text:
/* histogram/file2d.c
 *
 * Copyright (C) 1996, 1997, 1998, 1999, 2000, 2007 Brian Gough
 *
 * This program is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation; either version 3 of the License, or (at
 * your option) any later version.
 *
 * This program is distributed in the hope that it will be useful, but
 * WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 * General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this program; if not, write to the Free Software
 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
 */

#include <config.h>
#include <stdio.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_block.h>
#include <gsl/gsl_histogram2d.h>

int
gsl_histogram2d_fread (FILE * stream, gsl_histogram2d * h)
{
  int status = gsl_block_raw_fread (stream, h->xrange, h->nx + 1, 1);

  if (status)
    return status;

  status = gsl_block_raw_fread (stream, h->yrange, h->ny + 1, 1);

  if (status)
    return status;

  status = gsl_block_raw_fread (stream, h->bin, h->nx * h->ny, 1);

  return status;
}

int
gsl_histogram2d_fwrite (FILE * stream, const gsl_histogram2d * h)
{
  int status = gsl_block_raw_fwrite (stream, h->xrange, h->nx + 1, 1);

  if (status)
    return status;

  status = gsl_block_raw_fwrite (stream, h->yrange, h->ny + 1, 1);

  if (status)
    return status;

  status = gsl_block_raw_fwrite (stream, h->bin, h->nx * h->ny, 1);
  return status;
}

int
gsl_histogram2d_fprintf (FILE * stream, const gsl_histogram2d * h,
                         const char *range_format, const char *bin_format)
{
  size_t i, j;
  const size_t nx = h->nx;
  const size_t ny = h->ny;
  int status;

  for (i = 0; i < nx; i++)
    {
      for (j = 0; j < ny; j++)
        {
          status = fprintf (stream, range_format, h->xrange[i]);

          if (status < 0)
            {
              GSL_ERROR ("fprintf failed", GSL_EFAILED);
            }

          status = putc (' ', stream);

          if (status == EOF)
            {
              GSL_ERROR ("putc failed", GSL_EFAILED);
            }

          status = fprintf (stream, range_format, h->xrange[i + 1]);

          if (status < 0)
            {
              GSL_ERROR ("fprintf failed", GSL_EFAILED);
            }

          status = putc (' ', stream);

          if (status == EOF)
            {
              GSL_ERROR ("putc failed", GSL_EFAILED);
            }

          status = fprintf (stream, range_format, h->yrange[j]);

          if (status < 0)
            {
              GSL_ERROR ("fprintf failed", GSL_EFAILED);
            }

          status = putc (' ', stream);

          if (status == EOF)
            {
              GSL_ERROR ("putc failed", GSL_EFAILED);
            }

          status = fprintf (stream, range_format, h->yrange[j + 1]);

          if (status < 0)
            {
              GSL_ERROR ("fprintf failed", GSL_EFAILED);
            }

          status = putc (' ', stream);

          if (status == EOF)
            {
              GSL_ERROR ("putc failed", GSL_EFAILED);
            }

          status = fprintf (stream, bin_format, h->bin[i * ny + j]);

          if (status < 0)
            {
              GSL_ERROR ("fprintf failed", GSL_EFAILED);
            }

          status = putc ('\n', stream);

          if (status == EOF)
            {
              GSL_ERROR ("putc failed", GSL_EFAILED);
            }
        }
      status = putc ('\n', stream);

      if (status == EOF)
        {
          GSL_ERROR ("putc failed", GSL_EFAILED);
        }
    }

  return GSL_SUCCESS;
}

int
gsl_histogram2d_fscanf (FILE * stream, gsl_histogram2d * h)
{
  size_t i, j;
  const size_t nx = h->nx;
  const size_t ny = h->ny;
  double xupper, yupper;

  for (i = 0; i < nx; i++)
    {
      for (j = 0; j < ny; j++)
        {
          int status = fscanf (stream,
                               "%lg %lg %lg %lg %lg",
                               h->xrange + i, &xupper,
                               h->yrange + j, &yupper,
                               h->bin + i * ny + j);

          if (status != 5)
            {
              GSL_ERROR ("fscanf failed", GSL_EFAILED);
            }
        }
      h->yrange[ny] = yupper;
    }

  h->xrange[nx] = xupper;

  return GSL_SUCCESS;
}

max_stars_repo_name: ielomariala/Hex-Game

max_issues_repo_name: Brian-ning/HMNE

max_forks_repo_name: Brian-ning/HMNE


text:
/**
 *
 * @file testing_cgetri.c
 *
 *  PLASMA testing routines
 *  PLASMA is a software package provided by Univ. of Tennessee,
 *  Univ. of California Berkeley and Univ. of Colorado Denver
 *
 * @version 2.6.0
 * @author Hatem Ltaief
 * @date 2010-11-15
 * @generated c Tue Jan  7 11:45:17 2014
 *
 **/
#include <stdlib.h>
#include <stdio.h>
#include <string.h>
#include <math.h>

#include <plasma.h>
#include <cblas.h>
#include <lapacke.h>
#include <core_blas.h>
#include "testing_cmain.h"

static int check_factorization(int, PLASMA_Complex32_t*, PLASMA_Complex32_t*, int, int*, float);
static int check_inverse(int, PLASMA_Complex32_t *, PLASMA_Complex32_t *, int, int*, float);

int testing_cgetri(int argc, char **argv)
{

    /* Check for number of arguments*/
    if (argc != 2){
        USAGE("GETRI", "N LDA",
              "   - N    : the size of the matrix\n"
              "   - LDA  : leading dimension of the matrix A\n");
        return -1;
    }

    int N     = atoi(argv[0]);
    int LDA   = atoi(argv[1]);
    float eps;
    int info_inverse, info_factorization;
    int i, j;

    PLASMA_Complex32_t *A1   = (PLASMA_Complex32_t *)malloc(LDA*N*sizeof(PLASMA_Complex32_t));
    PLASMA_Complex32_t *A2   = (PLASMA_Complex32_t *)malloc(LDA*N*sizeof(PLASMA_Complex32_t));
    PLASMA_Complex32_t *WORK = (PLASMA_Complex32_t *)malloc(2*LDA*sizeof(PLASMA_Complex32_t));
    float *D                = (float *)malloc(LDA*sizeof(float));
    int *IPIV = (int *)malloc(N*sizeof(int));

    /* Check if unable to allocate memory */
    if ( (!A1) || (!A2) || (!IPIV) ){
        printf("Out of Memory \n ");
        return -2;
    }

    eps = LAPACKE_slamch_work('e');

    /*-------------------------------------------------------------
    *  TESTING CGETRI
    */

    /* Initialize A1 and A2 Matrix */
    PLASMA_cplrnt(N, N, A1, LDA, 3453);
    for ( i = 0; i < N; i++)
        for (  j = 0; j < N; j++)
            A2[LDA*j+i] = A1[LDA*j+i];

    printf("\n");
    printf("------ TESTS FOR PLASMA CGETRI ROUTINE -------  \n");
    printf("            Size of the Matrix %d by %d\n", N, N);
    printf("\n");
    printf(" The matrix A is randomly generated for each test.\n");
    printf("============\n");
    printf(" The relative machine precision (eps) is to be %e \n", eps);
    printf(" Computational tests pass if scaled residuals are less than 60.\n");

    /* PLASMA CGETRF */
    PLASMA_cgetrf(N, N, A2, LDA, IPIV);

    /* Check the factorization */
    info_factorization = check_factorization( N, A1, A2, LDA, IPIV, eps);

    /* PLASMA CGETRI */
    PLASMA_cgetri(N, A2, LDA, IPIV);

    /* Check the inverse */
    info_inverse = check_inverse(N, A1, A2, LDA, IPIV, eps);

    if ( (info_inverse == 0) && (info_factorization == 0) ) {
        printf("***************************************************\n");
        printf(" ---- TESTING CGETRI ..................... PASSED !\n");
        printf("***************************************************\n");
    }
    else {
        printf("***************************************************\n");
        printf(" - TESTING CGETRI ... FAILED !\n");
        printf("***************************************************\n");
    }

    free(A1); free(A2); free(IPIV); free(WORK); free(D);

    return 0;
}

/*------------------------------------------------------------------------
 *  Check the factorization of the matrix A2
 */
static int check_factorization(int N, PLASMA_Complex32_t *A1, PLASMA_Complex32_t *A2, int LDA, int *IPIV, float eps)
{
    int info_factorization;
    float Rnorm, Anorm, Xnorm, Bnorm, result;
    PLASMA_Complex32_t alpha, beta;

    alpha = 1.0;
    beta  = -1.0;

    PLASMA_Complex32_t *b  = (PLASMA_Complex32_t *)malloc(LDA*sizeof(PLASMA_Complex32_t));
    PLASMA_Complex32_t *x  = (PLASMA_Complex32_t *)malloc(LDA*sizeof(PLASMA_Complex32_t));

    LAPACKE_clarnv_work(1, ISEED, LDA, x);
    LAPACKE_clacpy_work(LAPACK_COL_MAJOR, 'A', N, 1, x, LDA, b, LDA);

    PLASMA_cgetrs( PlasmaNoTrans, N, 1, A2, LDA, IPIV, x, LDA );

    Xnorm = PLASMA_clange(PlasmaInfNorm, N, 1, x,  LDA);
    Anorm = PLASMA_clange(PlasmaInfNorm, N, N, A1, LDA);
    Bnorm = PLASMA_clange(PlasmaInfNorm, N, 1, b,  LDA);

    PLASMA_cgemm( PlasmaNoTrans, PlasmaNoTrans, N, 1, N,
                  alpha, A1, LDA, x, LDA, beta, b, LDA);

    Rnorm = PLASMA_clange(PlasmaInfNorm, N, 1, b, LDA);

    if (getenv("PLASMA_TESTING_VERBOSE"))
      printf( "||A||_oo=%f\n||X||_oo=%f\n||B||_oo=%f\n||A X - B||_oo=%e\n", Anorm, Xnorm, Bnorm, Rnorm );

    result = Rnorm / ( (Anorm*Xnorm+Bnorm)*N*eps ) ;
    printf("============\n");
    printf("Checking the Residual of the solution \n");
    printf("-- ||Ax-B||_oo/((||A||_oo||x||_oo+||B||_oo).N.eps) = %e \n", result);

    if (  isnan(Xnorm) || isinf(Xnorm) || isnan(result) || isinf(result) || (result > 60.0) ) {
        printf("-- The factorization is suspicious ! \n");
        info_factorization = 1;
     }
    else{
        printf("-- The factorization is CORRECT ! \n");
        info_factorization = 0;
    }
    free(x); free(b);

    return info_factorization;
}

/*------------------------------------------------------------------------
 *  Check the accuracy of the computed inverse
 */

static int check_inverse(int N, PLASMA_Complex32_t *A1, PLASMA_Complex32_t *A2, int LDA, int *IPIV, float eps )
{
    int info_inverse;
    int i;
    float Rnorm, Anorm, Ainvnorm, result;
    PLASMA_Complex32_t alpha, beta, zone;
    PLASMA_Complex32_t *work = (PLASMA_Complex32_t *)malloc(N*N*sizeof(PLASMA_Complex32_t));

    alpha = -1.0;
    beta  = 0.0;
    zone = 1.0;

    PLASMA_cgemm( PlasmaNoTrans, PlasmaNoTrans, N, N, N, alpha, A2, LDA, A1, LDA, beta, work, N);

    /* Add the identity matrix to work */
    for(i=0; i<N; i++)
        *(work+i+i*N) = *(work+i+i*N) + zone;

    Rnorm    = PLASMA_clange(PlasmaInfNorm, N, N, work, N);
    Anorm    = PLASMA_clange(PlasmaInfNorm, N, N, A1,   LDA);
    Ainvnorm = PLASMA_clange(PlasmaInfNorm, N, N, A2,   LDA);

    if (getenv("PLASMA_TESTING_VERBOSE"))
      printf( "||A||_1=%f\n||Ainv||_1=%f\n||Id - A*Ainv||_1=%e\n", Anorm, Ainvnorm, Rnorm );

    result = Rnorm / ( (Anorm*Ainvnorm)*N*eps ) ;
    printf("============\n");
    printf("Checking the Residual of the inverse \n");
    printf("-- ||Id - A*Ainv||_1/((||A||_1||Ainv||_1).N.eps) = %e \n", result);

    if (  isnan(Ainvnorm) || isinf(Ainvnorm) || isnan(result) || isinf(result) || (result > 60.0) ) {
        printf("-- The inverse is suspicious ! \n");
        info_inverse = 1;
     }
    else{
        printf("-- The inverse is CORRECT ! \n");
        info_inverse = 0;
    }

    free(work);

    return info_inverse;
}

ext: c


text:
#ifndef __GSL_PERMUTE_VECTOR_H__
#define __GSL_PERMUTE_VECTOR_H__

#if !defined( GSL_FUN )
#  if !defined( GSL_DLL )
#    define GSL_FUN extern
#  elif defined( BUILD_GSL_DLL )
#    define GSL_FUN extern __declspec(dllexport)
#  else
#    define GSL_FUN extern __declspec(dllimport)
#  endif
#endif

#include <gsl/gsl_permute_vector_complex_long_double.h>
#include <gsl/gsl_permute_vector_complex_double.h>
#include <gsl/gsl_permute_vector_complex_float.h>

#include <gsl/gsl_permute_vector_long_double.h>
#include <gsl/gsl_permute_vector_double.h>
#include <gsl/gsl_permute_vector_float.h>

#include <gsl/gsl_permute_vector_ulong.h>
#include <gsl/gsl_permute_vector_long.h>

#include <gsl/gsl_permute_vector_uint.h>
#include <gsl/gsl_permute_vector_int.h>

#include <gsl/gsl_permute_vector_ushort.h>
#include <gsl/gsl_permute_vector_short.h>

#include <gsl/gsl_permute_vector_uchar.h>
#include <gsl/gsl_permute_vector_char.h>

#endif /* __GSL_PERMUTE_VECTOR_H__ */

ext: h

max_stars_repo_name: zzpwahaha/Chimera-Control-Trim

max_issues_repo_name: zzpwahaha/Chimera-Control-Trim

max_forks_repo_name: zzpwahaha/Chimera-Control-Trim


text:
#ifndef H_MATH_UTILS_GSL
#define H_MATH_UTILS_GSL

#include <gsl/gsl_math.h>
//#include <egsl/egsl.h>
#include </home/alecone/catkin_graph/src/csm/sm/lib/egsl/egsl.h>

#include "laser_data.h"

#define gvg gsl_vector_get
#define gvs gsl_vector_set

/* GSL stuff */
	const char* gsl_friendly_pose(gsl_vector*v);
	gsl_vector * vector_from_array(unsigned int n, double *x);
	void vector_to_array(const gsl_vector*v, double*);
	void copy_from_array(gsl_vector*v, double*);

	void oplus(const gsl_vector*x1,const gsl_vector*x2, gsl_vector*res);
	void ominus(const gsl_vector*x, gsl_vector*res);
	void pose_diff(const gsl_vector*pose2,const gsl_vector*pose1,gsl_vector*res);

	void transform(const gsl_vector* point2d, const gsl_vector* pose, gsl_vector*result2d);
	void gsl_vector_set_nan(gsl_vector*v);

	double distance(const gsl_vector* a,const gsl_vector* b);
	double distance_squared(const gsl_vector* a,const gsl_vector* b);

	/** Returns norm of 2D point p */
	double norm(const gsl_vector*p);
	const char* egsl_friendly_pose(val pose);
	const char* egsl_friendly_cov(val cov);

/** Returns Fisher's information matrix. You still have to multiply
    it by (1/sigma^2). */
val ld_fisher0(LDP ld);

#endif


text:
/* specfunc/bessel_Jnu.c
 *
 * Copyright (C) 1996, 1997, 1998, 1999, 2000 Gerard Jungman, 2017 Konrad Griessinger
 *
 * This program is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation; either version 3 of the License, or (at
 * your option) any later version.
 *
 * This program is distributed in the hope that it will be useful, but
 * WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 * General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this program; if not, write to the Free Software
 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
 */

/* Author:  G. Jungman */

#include <config.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_sf_bessel.h>
#include <gsl/gsl_sf_sincos_pi.h>

#include "error.h"

#include "bessel.h"
#include "bessel_olver.h"
#include "bessel_temme.h"

/* Evaluate at large enough nu to apply asymptotic
 * results and apply backward recurrence.
 */
#if 0
static
int
bessel_J_recur_asymp(const double nu, const double x,
                     gsl_sf_result * Jnu, gsl_sf_result * Jnup1)
{
  const double nu_cut = 25.0;
  int n;
  int steps = ceil(nu_cut - nu) + 1;

  gsl_sf_result r_Jnp1;
  gsl_sf_result r_Jn;
  int stat_O1 = gsl_sf_bessel_Jnu_asymp_Olver_e(nu + steps + 1.0, x, &r_Jnp1);
  int stat_O2 = gsl_sf_bessel_Jnu_asymp_Olver_e(nu + steps,       x, &r_Jn);
  double r_fe = fabs(r_Jnp1.err/r_Jnp1.val) + fabs(r_Jn.err/r_Jn.val);
  double Jnp1 = r_Jnp1.val;
  double Jn   = r_Jn.val;
  double Jnm1;
  double Jnp1_save;

  for(n=steps; n>0; n--) {
    Jnm1 = 2.0*(nu+n)/x * Jn - Jnp1;
    Jnp1 = Jn;
    Jnp1_save = Jn;
    Jn   = Jnm1;
  }

  Jnu->val = Jn;
  Jnu->err = (r_fe + GSL_DBL_EPSILON * (steps + 1.0)) * fabs(Jn);
  Jnup1->val = Jnp1_save;
  Jnup1->err = (r_fe + GSL_DBL_EPSILON * (steps + 1.0)) * fabs(Jnp1_save);

  return GSL_ERROR_SELECT_2(stat_O1, stat_O2);
}
#endif

/*-*-*-*-*-*-*-*-*-*-*-* Functions with Error Codes *-*-*-*-*-*-*-*-*-*-*-*/

int
gsl_sf_bessel_Jnupos_e(const double nu, const double x, gsl_sf_result * result)
{
  /* CHECK_POINTER(result) */

  if(x == 0.0) {
    if(nu == 0.0) {
      result->val = 1.0;
      result->err = 0.0;
    }
    else {
      result->val = 0.0;
      result->err = 0.0;
    }
    return GSL_SUCCESS;
  }
  else if(x*x < 10.0*(nu+1.0)) {
    return gsl_sf_bessel_IJ_taylor_e(nu, x, -1, 100, GSL_DBL_EPSILON, result);
  }
  else if(nu > 50.0) {
    return gsl_sf_bessel_Jnu_asymp_Olver_e(nu, x, result);
  }
  else if(x > 1000.0)
  {
    /* We need this to avoid feeding large x to CF1; note that
     * due to the above check, we know that n <= 50. See similar
     * block in bessel_Jn.c.
     */
    return gsl_sf_bessel_Jnu_asympx_e(nu, x, result);
  }
  else {
    /* -1/2 <= mu <= 1/2 */
    int N = (int)(nu + 0.5);
    double mu = nu - N;

    /* Determine the J ratio at nu.
     */
    double Jnup1_Jnu;
    double sgn_Jnu;
    const int stat_CF1 = gsl_sf_bessel_J_CF1(nu, x, &Jnup1_Jnu, &sgn_Jnu);

    if(x < 2.0) {
      /* Determine Y_mu, Y_mup1 directly and recurse forward to nu.
       * Then use the CF1 information to solve for J_nu and J_nup1.
       */
      gsl_sf_result Y_mu, Y_mup1;
      const int stat_mu = gsl_sf_bessel_Y_temme(mu, x, &Y_mu, &Y_mup1);

      double Ynm1 = Y_mu.val;
      double Yn   = Y_mup1.val;
      double Ynp1 = 0.0;
      int n;
      for(n=1; n<N; n++) {
        Ynp1 = 2.0*(mu+n)/x * Yn - Ynm1;
        Ynm1 = Yn;
        Yn   = Ynp1;
      }

      result->val = 2.0/(M_PI*x) / (Jnup1_Jnu*Yn - Ynp1);
      result->err = GSL_DBL_EPSILON * (N + 2.0) * fabs(result->val);
      return GSL_ERROR_SELECT_2(stat_mu, stat_CF1);
    }
    else {
      /* Recurse backward from nu to mu, determining the J ratio
       * at mu. Use this together with a Steed method CF2 to
       * determine the actual J_mu, and thus obtain the normalization.
       */
      double Jmu;
      double Jmup1_Jmu;
      double sgn_Jmu;
      double Jmuprime_Jmu;
      double P, Q;
      const int stat_CF2 = gsl_sf_bessel_JY_steed_CF2(mu, x, &P, &Q);
      double gamma;

      double Jnp1 = sgn_Jnu * GSL_SQRT_DBL_MIN * Jnup1_Jnu;
      double Jn   = sgn_Jnu * GSL_SQRT_DBL_MIN;
      double Jnm1;
      int n;
      for(n=N; n>0; n--) {
        Jnm1 = 2.0*(mu+n)/x * Jn - Jnp1;
        Jnp1 = Jn;
        Jn   = Jnm1;
      }
      Jmup1_Jmu = Jnp1/Jn;
      sgn_Jmu   = GSL_SIGN(Jn);
      Jmuprime_Jmu = mu/x - Jmup1_Jmu;

      gamma = (P - Jmuprime_Jmu)/Q;
      Jmu   = sgn_Jmu * sqrt(2.0/(M_PI*x) / (Q + gamma*(P-Jmuprime_Jmu)));

      result->val = Jmu * (sgn_Jnu * GSL_SQRT_DBL_MIN) / Jn;
      result->err = 2.0 * GSL_DBL_EPSILON * (N + 2.0) * fabs(result->val);

      return GSL_ERROR_SELECT_2(stat_CF2, stat_CF1);
    }
  }
}

int
gsl_sf_bessel_Jnu_e(const double nu, const double x, gsl_sf_result * result)
{
  /* CHECK_POINTER(result) */

  if(x <= 0.0) {
    DOMAIN_ERROR(result);
  }
  else if (nu < 0.0) {
    int Jstatus = gsl_sf_bessel_Jnupos_e(-nu, x, result);
    double Jval = result->val;
    double Jerr = result->err;
    int Ystatus = gsl_sf_bessel_Ynupos_e(-nu, x, result);
    double Yval = result->val;
    double Yerr = result->err;
    /* double s = sin(M_PI*nu), c = cos(M_PI*nu); */
    int sinstatus = gsl_sf_sin_pi_e(nu, result);
    double s = result->val;
    double serr = result->err;
    int cosstatus = gsl_sf_cos_pi_e(nu, result);
    double c = result->val;
    double cerr = result->err;
    result->val = s*Yval + c*Jval;
    result->err = fabs(c*Yerr) + fabs(s*Jerr) + fabs(cerr*Yval) + fabs(serr*Jval);
    return GSL_ERROR_SELECT_4(Jstatus, Ystatus, sinstatus, cosstatus);
  }
  else return gsl_sf_bessel_Jnupos_e(nu, x, result);
}

/*-*-*-*-*-*-*-*-*-* Functions w/ Natural Prototypes *-*-*-*-*-*-*-*-*-*-*/

#include "eval.h"

double gsl_sf_bessel_Jnu(const double nu, const double x)
{
  EVAL_RESULT(gsl_sf_bessel_Jnu_e(nu, x, &result));
}

max_issues_repo_name: ielomariala/Hex-Game

max_forks_repo_name: karanbirsandhu/nu-sense


text:
#pragma once

#include "options.h"
#include <cstdio>
#include <fmt/color.h>
#include <fmt/printf.h>
#include <gsl/gsl-lite.hpp>

namespace angonoka::cli {
/**
    Print text to stdout if not in quiet mode.

    @param options CLI options
*/
template <typename... T>
void print(
    const Options& options,
    fmt::format_string<T...> fmt,
    T&&... args)
{
    if (options.quiet) return;
    fmt::print(fmt, std::forward<T>(args)...);
}

/**
    Prints red text.

    Conditionally disables the color depending on
    CLI options.

    @param options CLI options
*/
template <typename... T>
void print_error(
    const Options& options,
    std::FILE* output,
    fmt::format_string<T...> fmt,
    T&&... args)
{
    if (options.color) {
        fmt::print(
            output,
            fg(fmt::terminal_color::red),
            fmt::string_view{fmt},
            std::forward<T>(args)...);
    } else {
        fmt::print(output, fmt, std::forward<T>(args)...);
    }
}

/**
    Critical error message.

    Used for progress messages with ellipsis like

    Progress message... <die()>Error
    An error has occured.

    @param options CLI options
*/
template <typename... T>
void die(
    const Options& options,
    fmt::format_string<T...> fmt,
    T&&... args)
{

    if (!options.quiet) print_error(options, stdout, "Error\n");
    print_error(options, stderr, fmt, std::forward<T>(args)...);
}
} // namespace angonoka::cli

max_stars_repo_name: coffee-lord/angonoka

max_issues_repo_name: coffee-lord/angonoka

max_forks_repo_name: coffee-lord/angonoka


text:
/*
How to use blas and lapack with the standard
interfaces provided by their respective projects, repectively through
`cblas.h` and `lapacke.h`
*/

#include <assert.h>
#include <math.h>
#include <stdbool.h>
#include <stdio.h>
#include <stdlib.h>

#include <cblas.h>
#include <lapacke.h>

/**
 * assert two integers are equal
 * if not, print them to stderr and assert false
 */
void assert_eqi(int i1, int i2) {
    if (i1 != i2) {
        fprintf(stderr, "%d\n%d\n\n\n", i1, i2);
        assert(false);
    }
}

/**
 * assert two doubles are equal within err precision
 * if not, print them to stderr
 */
void assert_eqd(double d1, double d2, double err) {
    if (fabs(d1 - d2) > err) {
        fprintf(stderr, "%f\n%f\n\n\n", d1, d2);
        assert(false);
    }
}

/** print an array of doubles to stderr */
void print_vecd(int n, double * v) {
    int i;
    for (i=0; i<n; i++) {
        fprintf(stderr, "%.2f ", v[i]);
    }
    printf("\n");
}

/*
bool eq_vecd2(int n, double * v1, double * v2, double err)
{
    if(fabs(d1 - d2) > err)
        return false
    return true;
}

void assert_eqvd(int n, double * v1, double * v2, double err){
    if(eq_vecd(n, v1, v2, err) != true){
        print_vecd(n, v1);
        print_vecd(n, v2);
    }
}
*/

int main(void) {
    int info, ipiv2[2];
    float err = 1e-6;
    float x2[2], b2[2], c2[2];
    float a2x2[2][2];

	/* cblas */
	{
        x2[0] =  1.0;
        x2[1] = -2.0;
        assert_eqd(cblas_snrm2(2, x2, 1), sqrt(5.0), err);
    }

	/* lapacke */
	{
		/*
		sgesv

		Matrix vector multiply.
		*/
		{
			a2x2[0][0] = 1.0;
			a2x2[1][0] = 2.0;
			a2x2[0][1] = 3.0;
			a2x2[1][1] = 4.0;
			b2[0] = 5.;
			b2[1] = 11.;

			info = LAPACKE_sgesv(
				LAPACK_COL_MAJOR,
				2,
				1,
				&a2x2[0][0],
				2,
				&ipiv2[0],
				&b2[0],
				2
			);
			c2[0] = 1.0;
			c2[1] = 2.0;
			assert_eqi(info, 0);
			assert_eqd(b2[0], c2[0], err);
			assert_eqd(b2[1], c2[1], err);
		}
    }

    return EXIT_SUCCESS;
}

max_stars_repo_name: liujiamingustc/phd

max_issues_repo_name: liujiamingustc/phd

max_forks_repo_name: liujiamingustc/phd


text:
/*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*
**                                                                                **
** This file forms part of the Underworld geophysics modelling application.         **
**                                                                                **
** For full license and copyright information, please refer to the LICENSE.md file  **
** located at the project root, or contact the authors.                             **
**                                                                                **
**~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*/
#include <petsc.h>
#include <petscmat.h>

#include "common-driver-utils.h"

#include <StGermain/libStGermain/src/StGermain.h>
#include <StgDomain/libStgDomain/src/StgDomain.h>
#include <StgFEM/libStgFEM/src/StgFEM.h>
#include <PICellerator/libPICellerator/src/PICellerator.h>
#include <Underworld/libUnderworld/src/Underworld.h>
#include "Solvers/SLE/src/SLE.h" /* to give the AugLagStokes_SLE type */
#include "Solvers/KSPSolvers/src/KSPSolvers.h" /* for __KSP_COMMON */
#include "BSSCR.h"
#include "writeMatVec.h"

#if( (PETSC_VERSION_MAJOR==2) && (PETSC_VERSION_MINOR==3) && (PETSC_VERSION_SUBMINOR==0) )
#define FILE_OPTION PETSC_FILE_RDONLY
#endif
#if( (PETSC_VERSION_MAJOR==2) && (PETSC_VERSION_MINOR==3) && (PETSC_VERSION_SUBMINOR>=2) )
#define FILE_OPTION FILE_MODE_READ
#endif
#if( PETSC_VERSION_MAJOR==3 )
#define FILE_OPTION FILE_MODE_READ
#endif

PetscErrorCode BSSCR_FormSchurApproximation1( Mat A11, Mat A12, Mat A21, Mat A22, Mat *_Shat, PetscTruth sym );
PetscErrorCode BSSCR_FormSchurApproximationDiag( Mat A11, Mat A12, Mat A21, Mat A22, Mat *_Shat, PetscTruth sym );

PetscErrorCode BSSCR_BSSCR_StokesReadPCSchurMat_binary( MPI_Comm comm, Mat *S )
{
	PetscViewer mat_view_file;
	char	op_name[PETSC_MAX_PATH_LEN];
	PetscTruth flg;

	PetscOptionsGetString( PETSC_NULL,"-stokes_Smat",op_name,PETSC_MAX_PATH_LEN-1,&flg );
	*S = PETSC_NULL;
	if (flg) {
		if (!S)	Stg_SETERRQ(1,"Memory space for Smat is NULL");

		PetscViewerBinaryOpen( comm, op_name,  FILE_OPTION, &mat_view_file );
		Stg_MatLoad( mat_view_file, MATAIJ, S );
		Stg_PetscViewerDestroy(&mat_view_file );
	}

	PetscFunctionReturn(0);
}

PetscErrorCode BSSCR_BSSCR_StokesReadPCSchurMat_ascii( MPI_Comm comm, Mat *S )
{
	char	op_name[PETSC_MAX_PATH_LEN];
	PetscTruth flg;

	PetscOptionsGetString( PETSC_NULL,"-stokes_Smat",op_name,PETSC_MAX_PATH_LEN-1,&flg );
	*S = PETSC_NULL;
	if (flg) {
		if (!S)	Stg_SETERRQ(1,"Memory space for Smat is NULL");
		Stg_SETERRQ(1,"Currently Disabled");
		//MatAIJLoad_MatrixMarket( comm, op_name, S );
	}

	PetscFunctionReturn(0);
}

PetscErrorCode BSSCR_StokesReadPCSchurMat( MPI_Comm comm, Mat *S )
{
	char	op_name[PETSC_MAX_PATH_LEN];
	PetscTruth flg;

	PetscOptionsGetString( PETSC_NULL,"-stokes_ascii",op_name,PETSC_MAX_PATH_LEN-1,&flg );
	if (flg==PETSC_TRUE) {
		BSSCR_BSSCR_StokesReadPCSchurMat_ascii( comm, S );
	}
	else {
		BSSCR_BSSCR_StokesReadPCSchurMat_binary( comm, S );
	}

	PetscFunctionReturn(0);
}

PetscErrorCode BSSCR_StokesCreatePCSchur( Mat K, Mat G, PC pc_S )
{
	char pc_type[PETSC_MAX_PATH_LEN];
	PetscTruth flg;

	PetscOptionsGetString( PETSC_NULL, "-Q22_pc_type", pc_type, PETSC_MAX_PATH_LEN-1, &flg );
	if( !flg ) {
		Stg_SETERRQ( PETSC_ERR_SUP, "OPTION: -Q22_pc_type must be set" );
	}

	/* 1. define S pc to be "none" */

	if( strcmp(pc_type,"none")==0 ) { /* none */
		PetscPrintf(PETSC_COMM_WORLD,"  Setting schur_pc to \"none\" \n" );
		PCSetType( pc_S, "none" );
	}
	else if( strcmp(pc_type,"uw")==0 ) { /* diag */
		Mat S, Amat;
		MatStructure mstruct;
		MPI_Comm comm;

		PetscPrintf(PETSC_COMM_WORLD,"  Setting schur_pc to \"uw\" \n" );

		PetscObjectGetComm( (PetscObject)pc_S, &comm );
		S = PETSC_NULL;
		BSSCR_StokesReadPCSchurMat( comm, &S );
		if (!S) {	Stg_SETERRQ(1,"Must indicate location of file for SchurPC matrix with the stokes_Smat option");	}

		Stg_PCGetOperators( pc_S, &Amat, &Amat, &mstruct );
		Stg_PCSetOperators( pc_S, Amat, S, SAME_NONZERO_PATTERN );

		MatView( S, PETSC_VIEWER_STDOUT_WORLD );
	}
	else { /* not valid option */
		Stg_SETERRQ( PETSC_ERR_SUP, "OPTION: -seg_schur_pc_type is not valid" );
	}

	PetscFunctionReturn(0);
}

PetscErrorCode BSSCR_BSSCR_StokesCreatePCSchur2(
   Mat K, Mat G, Mat D, Mat C, Mat Smat, PC pc_S,
   PetscTruth sym, KSP_BSSCR * bsscrp )
{
	char pc_type[PETSC_MAX_PATH_LEN];
	PetscTruth flg;

	PetscOptionsGetString( PETSC_NULL, "-Q22_pc_type", pc_type, PETSC_MAX_PATH_LEN-1, &flg );
	if( !flg ) {
	    strcpy(pc_type, "uw");
	    //Stg_SETERRQ( PETSC_ERR_SUP, "OPTION: -Q22_pc_type must be set" );
	}

	/* 1. define S pc to be "none" */

	if( strcmp(pc_type,"none")==0 ) { /* none */
		/* Mat Amat,Pmat; */
		/* MatStructure mstruct; */

		PetscPrintf(PETSC_COMM_WORLD,"  Setting schur_pc to \"none\" \n" );
		/* PCGetOperators( pc_S, &Amat, &Pmat, &mstruct ); */
		/* PCSetOperators( pc_S, Amat, Amat, SAME_NONZERO_PATTERN ); */
		PCSetType( pc_S, "none" );
	}
	else if( strcmp(pc_type,"uw")==0 ) { /* diag */
	    Mat Amat,Pmat;
		MatStructure mstruct;

		PetscPrintf(PETSC_COMM_WORLD,"  Setting schur_pc to \"uw\" \n" );

		if (!Smat) {	Stg_SETERRQ(1,"Smat cannot be NULL if -Q22_pc_type = uw");	}

		Stg_PCGetOperators( pc_S, &Amat, &Pmat, &mstruct );
		Stg_PCSetOperators( pc_S, Amat, Smat, SAME_NONZERO_PATTERN );
	}
	else if( strcmp(pc_type,"uwscale")==0 ) { /* diag */
	    Mat Amat, Shat, Pmat;
		MatStructure mstruct;

		PetscPrintf(PETSC_COMM_WORLD,"  Setting schur_pc to \"uwscale\" \n" );

		if (!Smat) {	Stg_SETERRQ(1,"Smat cannot be NULL if -Q22_pc_type = uwscale");	}

		BSSCR_FormSchurApproximation1( K, G, D, C, &Shat, sym );

		Stg_PCGetOperators( pc_S, &Amat, &Pmat, &mstruct );
		Stg_PCSetOperators( pc_S, Amat, Shat, SAME_NONZERO_PATTERN );

		Stg_MatDestroy(&Shat);
	}
	else if( strcmp(pc_type,"gkgdiag")==0 ) { /* diag */
	    Mat Amat, Shat, Pmat;
		MatStructure mstruct;

		PetscPrintf(PETSC_COMM_WORLD,"  Setting schur_pc to \"gkgdiag\" \n" );

		if (!Smat) {	Stg_SETERRQ(1,"Smat cannot be NULL if -Q22_pc_type = uwscale");	}

		BSSCR_FormSchurApproximationDiag( K, G, D, C, &Shat, sym );

		Stg_PCGetOperators( pc_S, &Amat, &Pmat, &mstruct );
		Stg_PCSetOperators( pc_S, Amat, Shat, SAME_NONZERO_PATTERN );

		Stg_MatDestroy(&Shat);
	}
	else if( strcmp(pc_type,"gtkg")==0 ) { /* GtKG */
	    PetscPrintf(PETSC_COMM_WORLD,"  Setting schur_pc to \"gtkg\" \n" );

	    /* Build the schur pc GtKG */
	    PCSetType( pc_S, "gtkg" );
	    //AugLagStokes_SLE * stokesSLE = (AugLagStokes_SLE*)bsscrp->st_sle;
        StokesBlockKSPInterface* Solver = bsscrp->solver;
	    Mat M=0;
	    if (Solver->vmStiffMat){
		M = Solver->vmStiffMat->matrix;
	    }
	    BSSCR_PCGtKGSet_Operators( pc_S, K, G, M );
	    //BSSCR_PCGtKGAttachNullSpace( pc_S );
	}
	else { /* not valid option */
		Stg_SETERRQ( PETSC_ERR_SUP, "OPTION: -Q22_pc_type is not valid" );
	}

	PetscFunctionReturn(0);
}

/*

Recomended usage:

Stg_KSPSetOperators( ksp_S, S, Shat, SAME_NONZERO_PATTERN );

KSPGetPC( ksp_S, &pc_S );
PCSetType( pc_S, "jacobi" );

NOTE: In pracise "jacobi" works better than keep the entire matrix and
factoring it using something like "cholesky"

*/
PetscErrorCode BSSCR_FormSchurApproximation1( Mat A11, Mat A12, Mat A21, Mat A22, Mat *_Shat, PetscTruth sym )
{
   Mat Shat, A21_cpy;
   Vec diag;

   MatGetVecs( A11, &diag, PETSC_NULL );
   MatGetDiagonal( A11, diag );
   VecReciprocal( diag );

/*    if( sym ) { */
/* #if( PETSC_VERSION_MAJOR <= 2 ) */
/*       MatTranspose( A12, &A21_cpy ); */
/* #else */
/*       MatTranspose( A12, MAT_INITIAL_MATRIX, &A21_cpy ); */
/* #endif */
/*       MatDiagonalScale(A21_cpy, PETSC_NULL, diag ); */
/*    } */
/*    else { */
      MatDuplicate( A21, MAT_COPY_VALUES, &A21_cpy );
      MatDiagonalScale(A21_cpy, PETSC_NULL, diag );
   /* } */

   MatMatMult( A21_cpy, A12, MAT_INITIAL_MATRIX, 1.2, &Shat );  /* A21 diag(K)^{-1} A12 */

   if( A22 != PETSC_NULL )
      MatAXPY( Shat, -1.0, A22, DIFFERENT_NONZERO_PATTERN ); /* S <- -C + A21 diag(K)^{-1} A12 */

   *(void**)_Shat = (void*)Shat;

   Stg_MatDestroy(&A21_cpy);
   Stg_VecDestroy(&diag);

   PetscFunctionReturn(0);
}

PetscErrorCode BSSCR_FormSchurApproximationDiag( Mat A11, Mat A12, Mat A21, Mat A22, Mat *_Shat, PetscTruth sym )
{
   Mat Shat, A21_cpy;
   Vec diag;

   MatGetVecs( A11, &diag, PETSC_NULL );
   MatGetDiagonal( A11, diag );
   VecReciprocal( diag );

/*    if( sym ) { */
/* #if( PETSC_VERSION_MAJOR <= 2 ) */
/*       MatTranspose( A12, &A21_cpy ); */
/* #else */
/*       MatTranspose( A12, MAT_INITIAL_MATRIX, &A21_cpy ); */
/* #endif */
/*       MatDiagonalScale(A21_cpy, PETSC_NULL, diag ); */
/*    } */
/*    else { */
      MatDuplicate( A21, MAT_COPY_VALUES, &A21_cpy );
      MatDiagonalScale(A21_cpy, PETSC_NULL, diag );
   /* } */

   MatMatMult( A21_cpy, A12, MAT_INITIAL_MATRIX, 1.2, &Shat );  /* A21 diag(K)^{-1} A12 */

   if( A22 != PETSC_NULL )
      MatAXPY( Shat, -1.0, A22, DIFFERENT_NONZERO_PATTERN ); /* S <- -C + A21 diag(K)^{-1} A12 */

   Stg_MatDestroy(&A21_cpy);
   Stg_VecDestroy(&diag);

   MatGetVecs( Shat, &diag, PETSC_NULL );
   MatGetDiagonal( Shat, diag );

   MatZeroEntries( Shat );
   MatDiagonalSet( Shat, diag, INSERT_VALUES );

   *(void**)_Shat = (void*)Shat;

   Stg_VecDestroy(&diag);

   PetscFunctionReturn(0);
}


text:
#include <stdio.h>
#include <stdlib.h>
#include <time.h>
#include <math.h>
#include <string.h>
#include <gsl/gsl_cblas.h>
#include <gsl/gsl_randist.h>

#include "allvars.h"
#include "proto.h"

void mcmc_con_run()
{
  char fname_mcmc[100];
  strcpy(fname_mcmc, "data/mcmc_con.txt");

  printf("*******con mcmc\n");
  mcmc_con_init();
  if(flag_mcmc==1)
  {
    mcmc_sampling(fname_mcmc, &probability_con);
    mcmc_stats(fname_mcmc);
  }
  else
  {
    nc = nc_lim_low;
    printf("reading par.txt\n");
    read_input();
  }

  reconstruct_con();

  memcpy(theta_best_con, theta_best, ntheta*sizeof(double));
  memcpy(theta_best_var_con, theta_best_var, ntheta*2*sizeof(double));
}

void mcmc_con_init()
{
  int i;

  ntheta = 2;

  i=0;
  theta_range[i][0] = log(1.0e-6);
  theta_range[i++][1] = log(10.0);
  theta_range[i][0] = log(1.0e-2);
  theta_range[i++][1] = log(1.0e5);

  i = 0;
  theta_input[i++] = log(0.1);
  theta_input[i++] = log(100.0);

  i = 0;
  sigma_input[i++] = 0.1;
  sigma_input[i++] = 0.5;

  i = 0;
  theta_fixed[i++] = 0;
  theta_fixed[i++] = 0;

}
/*
 * probability for continuum mcmc.
 */
double probability_con(double *theta)//double sigmahat, double taud, double alpha)
{
  double sigmahat, taud, alpha;
  double prob, prior, lndet, lndet_ICq, sigma;
  double *ybuf, *Larr, *Cq, *ICq, *ave, *yave, *ysub;
  int i, nq, info;

  nq = (flag_detrend + 1);
  ybuf = workspace;
  Larr = workspace+ncon_data;
  Cq = Larr + nq*ncon_data;
  ICq = Cq + nq*nq;
  ave = ICq + nq*nq;
  yave = ave + nq;
  ysub = yave + ncon_data;

  alpha = 1.0;
  sigmahat = exp(theta[0]);
  taud = exp(theta[1]);

  sigma = sigmahat * sqrt(taud / 2.0);

  set_covar_mat_con(sigma, taud, alpha);

  for(i=0;i<ncon_data*ncon_data; i++)
  {
    Cmat[i] = Smat[i] + Nmat[i];
  }
  memcpy(ICmat, Cmat, ncon_data*ncon_data*sizeof(double));
  inverse_mat(ICmat, ncon_data, &info);

  for(i=0;i<ncon_data;i++)
  {
    if(flag_detrend==0)
    {
      Larr[i]=1.0;
    }
    else
    {
      Larr[i*nq + 0] = 1.0;
      Larr[i*nq + 1] = Tcon_data[i];
    }
  }
  //multiply_matvec(ICmat, Larr, ncon_data, ybuf);
  multiply_mat_MN(ICmat, Larr, Tmat1, ncon_data, nq, ncon_data);
  multiply_mat_MN_transposeA(Larr, Tmat1, ICq, nq, nq, ncon_data);
  memcpy(Cq, ICq, nq*nq*sizeof(double));
  inverse_mat(Cq, nq, &info);

  multiply_mat_MN_transposeA(Larr, ICmat, Tmat1, nq, ncon_data, ncon_data);
  multiply_mat_MN(Cq, Tmat1, Tmat2, nq, ncon_data, nq);
  multiply_mat_MN(Tmat2, Fcon_data, ave, nq, 1, ncon_data);
  multiply_mat_MN(Larr, ave, yave, ncon_data, 1, nq);

//  lambda = cblas_ddot(ncon_data, Larr, 1, ybuf, 1);
//  multiply_matvec(ICmat, Fcon_data, ncon_data, ybuf);
//  ave_con = cblas_ddot(ncon_data, Larr, 1, ybuf, 1);
//  ave_con /=lambda;

  for(i=0;i<ncon_data;i++)ysub[i] = Fcon_data[i] - yave[i];
  multiply_matvec(ICmat, ysub, ncon_data, ybuf);
  prob = -0.5 * cblas_ddot(ncon_data, ysub, 1, ybuf, 1);

  lndet = lndet_mat(Cmat, ncon_data, &info);
  lndet_ICq = lndet_mat(ICq, nq, &info);
  prob = prob - 0.5*lndet - 0.5*lndet_ICq;

/* penalize on larger tau than the length of continuum */
  prior = 0.0;
  prior += - theta[0];
  if(theta[1] > log(cad_con) )
  {
    prior = ( log(cad_con) - theta[1]);
  }
  else
  {
    prior +=  - (log(cad_con) - theta[1]);
  }
  return prob + prior;
}

/* calculate covariance matrix */
void set_covar_mat_con(double sigma, double tau, double alpha)
{
  double t1, t2, nerr;
  int i, j;

  for(i=0; i<ncon_data; i++)
  {
    t1 = Tcon_data[i];
    for(j=0; j<=i; j++)
    {

      t2 = Tcon_data[j];
      Smat[i*ncon_data+j] = sigma*sigma * exp (- pow (fabs(t1-t2) / tau, alpha));
      Smat[j*ncon_data+i] = Smat[i*ncon_data+j];

      Nmat[i*ncon_data+j] = Nmat[j*ncon_data+i] = 0.0;
    }
    nerr = Fcerrs_data[i];
    Nmat[i*ncon_data+i] = nerr*nerr;
  }
  return;
}

void set_covar_Umat_con(double sigma, double tau, double alpha)
{
  double t1, t2;
  int i, j;

  for(i=0; i<ncon; i++)
  {
    t1 = Tcon[i];
    for(j=0; j<ncon_data; j++)
    {
      t2 = Tcon_data[j];
      USmat[i*ncon_data+j] = sigma*sigma * exp (- pow (fabs(t1-t2) / tau, alpha) );
    }
  }
  return;
}

/* reconstruction */
void reconstruct_con()
{
  double sigma, tau, alpha, sigmahat;
  double *ybuf, *Larr, *Larr_rec, *yave, *yave_rec, *Cq, *ICq, *ysub, *ave;
  double *Tmp1, *Tmp2, *Tmp3, *Tmp4;

  int i, nq, info, nmax;

  nq = 1+flag_detrend;

  yave = workspace;
  ysub = yave + ncon_data;
  Larr_rec = ysub + ncon_data;    //Larr_rec is a nq*ncon matrix
  yave_rec = Larr_rec + nq*ncon;  //yave_rec is ncon matrix
  ybuf = yave_rec + ncon;
  Larr = ybuf + ncon_data;   // Larr is a nq*n matrix
  Cq = Larr + nq*ncon_data;
  ICq = Cq + nq*nq;
  ave = ICq + nq*nq;

  nmax = max(ncon, ncon_data);
  Tmp1 = array_malloc(nmax*nmax);     // temporary matrixes
  Tmp2 = array_malloc(nmax*nmax);
  Tmp3 = array_malloc(nmax*nmax);
  Tmp4 = array_malloc(nmax*nmax);

  sigmahat = exp(theta_best[0]);
  tau = exp(theta_best[1]);
  sigma = sigmahat * sqrt(tau/2.0);
  alpha = 1.0;

  set_covar_mat_con(sigma, tau, alpha);
  set_covar_Umat_con(sigma, tau, alpha);

  for(i=0;i<ncon_data*ncon_data; i++)
  {
    Cmat[i] = Smat[i] + Nmat[i];
  }
  memcpy(ICmat, Cmat, ncon_data*ncon_data*sizeof(double));
  inverse_mat(ICmat, ncon_data, &info);

  /* the best estimate for average */
  for(i=0;i<ncon_data;i++)
  {
    if(flag_detrend==0)
    {
      Larr[i]=1.0;
    }
    else
    {
      Larr[i*nq + 0] = 1.0;
      Larr[i*nq + 1] = Tcon_data[i];
    }
  }

  for(i=0;i<ncon;i++)
  {
    if(flag_detrend==0)
    {
      Larr_rec[i]=1.0;
    }
    else
    {
      Larr_rec[i*nq + 0] = 1.0;
      Larr_rec[i*nq + 1] = Tcon[i];
    }
  }

// calculate the covariance matrix Cq for q.
// Cq = (L^T x C^-1 x L)^-1
  multiply_mat_MN(ICmat, Larr, Tmat1, ncon_data, nq, ncon_data);
  multiply_mat_MN_transposeA(Larr, Tmat1, ICq, nq, nq, ncon_data);
  memcpy(Cq, ICq, nq*nq*sizeof(double));
  inverse_mat(Cq, nq, &info);
// calculate the best estimate of q.
// q = Cq x L^T x C^-1 x y
  multiply_mat_MN_transposeA(Larr, ICmat, Tmat1, nq, ncon_data, ncon_data);
  multiply_mat_MN(Cq, Tmat1, Tmat2, nq, ncon_data, nq);
  multiply_mat_MN(Tmat2, Fcon_data, ave, nq, 1, ncon_data);
// subtract the linear trend
  multiply_mat_MN(Larr, ave, yave, ncon_data, 1, nq);
  for(i=0;i<ncon_data;i++)ysub[i] = Fcon_data[i] - yave[i];
// calculate the best estimate for s.
// s = S x C^-1 x (y - Lxq)
  multiply_matvec(ICmat, ysub, ncon_data, ybuf);
  multiply_matvec_MN(USmat, ncon, ncon_data, ybuf, Fcon);
  multiply_mat_MN(Larr_rec, ave, yave_rec, ncon, 1, nq);

  for(i=0; i<ncon; i++)Fcon[i] = Fcon[i] + yave_rec[i];

  multiply_mat_MN(USmat, ICmat, Tmp1, ncon, ncon_data, ncon_data);
  multiply_mat_MN_transposeB(Tmp1, USmat, Tmp2, ncon, ncon, ncon_data);
  multiply_mat_MN(Tmp1, Larr, Tmp3, ncon, nq, ncon_data);
  for(i=0; i<ncon*nq; i++)Tmp3[i] -= Larr_rec[i];
  multiply_mat_MN(Tmp3, Cq, Tmp1, ncon, nq, nq);
  multiply_mat_MN_transposeB(Tmp1, Tmp3, Tmp4, ncon, ncon, nq);

  for(i=0; i<ncon; i++)
  {
    Fcerrs[i] = sqrt(sigma*sigma - Tmp2[i*ncon+i] + Tmp4[i*ncon+i]);
  }

  FILE *fp;
  fp = fopen("data/scon.txt", "w");
  for(i=0; i<ncon; i++)
  {
    fprintf(fp, "%f %f %f\n", Tcon[i], Fcon[i]*scale_con, Fcerrs[i]*scale_con);
  }
  fclose(fp);
  free(Tmp1);
  free(Tmp2);
  free(Tmp3);
  free(Tmp4);
}

max_stars_repo_name: LiyrAstroph/MICA

max_issues_repo_name: LiyrAstroph/MICA

max_forks_repo_name: LiyrAstroph/MICA


text:
#ifndef main_H_INCLUDED
#define main_H_INCLUDED

#include "wf.h"
#include "glasma.h"
#include "jet.h"
#include "bfkl.h"
#include "time.h"
#include "single.h"
#include <gsl/gsl_errno.h>

char tag[256];

//*
//*
// This section of code initializes and reads in the
// following global variables from file
//*
//*
#define X_FIELDS \
   X(int, wfTAG, "%d") \
   X(int, A1, "%d") \
   X(int, A2, "%d") \
   X(double, rts, "%lf") \

#include "ReadParameters.cxx"
//

#endif

max_stars_repo_name: kdusling/mpc

max_issues_repo_name: kdusling/mpc

max_forks_repo_name: kdusling/mpc


text:
// Geocode a DEM stored in ASF .meta format.

// Standard libraries.
#include <limits.h>
#include <math.h>
#include <string.h>

// Libraries from packages outside ASF.
#include <gsl/gsl_math.h>

// Libraries developed at ASF.
#include <asf.h>
#include <asf_meta.h>
#include <asf_raster.h>
#include <float_image.h>
#include <libasf_proj.h>

// Headers defined by this library.
#include "asf_geocode.h"

// Return true iff (x, y) falls in the region
// ([0, image->size_x - 1), [0, image->size_y - 1))
static int
in_image (FloatImage *image, ssize_t x, ssize_t y)
{
  g_assert (image->size_x <= SSIZE_MAX && image->size_y <= SSIZE_MAX);
  ssize_t sz_x = image->size_x, sz_y = image->size_y;

  return 0 <= x && x < sz_x && 0 <= y && y < sz_y;
}

// Given a projected image with metadata md, an unproject_input
// function, and pixel indicies (x, y), return the latitude and
// longitude of x, y (relative to whatever libasf_proj input spheroid
// is in effect at the time this function is called.  Probably easiest
// to see the purpose of this routine from the context in which it is
// called.
static int
get_pixel_lat_long (const meta_parameters *md, projector_t unproject_input,
		    size_t x, size_t y, double *lat, double *lon)
{
  // Convenience aliases.
  meta_projection *pb = md->projection;
  project_parameters_t *pp = &(md->projection->param);

  double xpc = pb->startX + pb->perX * x;
  double ypc = pb->startY + pb->perY * y;

  // Pseudoprojected images are in degrees and the projection and
  // unprojection function deal in radians.
  if ( md->projection->type == LAT_LONG_PSEUDO_PROJECTION ) {
    xpc *= D2R;
    ypc *= D2R;
  }

  int return_code = unproject_input (pp, xpc, ypc, ASF_PROJ_NO_HEIGHT,
				     lat, lon, NULL, md->projection->datum);
  g_assert (return_code);

  return return_code;
}

int
geocode_dem (projection_type_t projection_type,	// What we are projection to.
	     project_parameters_t *pp,    // Parameters we project to.
	     datum_type_t datum,                // Datum we project to.
	     // Pixel size of output image, in output projection units
	     // (meters or possibly degrees, if we decide to support
	     // projecting to pseudoprojected form).
	     double pixel_size,
	     resample_method_t resample_method,	// How to resample pixels.
	     const char *input_image, // Base name of input image.
	     const meta_parameters *imd, // Input DEM image metadata.
	     const char *output_image  // Base name of output image.
	     )
{
  int return_code;		// Holds return codes from functions.

  // Function to use to project or unproject between latlon and input
  // or output coordinates.
  projector_t project_input; 	// latlon => input image map projection
  projector_t unproject_input;	// input image_map_projection => latlon
  projector_t project_output;	// latlon => output image map projection
  projector_t unproject_output;	// output image map projection => latlon
  // Like the above, but act on arrays.
  array_projector_t array_project_input, array_unproject_input;
  array_projector_t array_project_output, array_unproject_output;

  // We only deal with reprojection map projected DEMs.
  g_assert (imd->projection != NULL);

  // FIXME: what to do with background value is something that still
  // needs to be determined (probably in consultation with the guys
  // working on terrain correction).
  const float background_value = 0.0;

  // Geocoding to pseudoprojected form presents issues, for example
  // with the meaning of the pixel_size argument, which is taken as a
  // distance in map projection coordinates for all other projections
  // (deciding how to interpret it when projecting to pseudoprojected
  // form is tough), and since there probably isn't much need, we
  // don't allow it.
  g_assert (projection_type != LAT_LONG_PSEUDO_PROJECTION);

  // Get the functions we want to use for projecting and unprojecting.
  set_projection_functions (imd->projection->type, &project_input,
			    &unproject_input, &array_project_input,
			    &array_unproject_input);
  set_projection_functions (projection_type, &project_output,
			    &unproject_output, &array_project_output,
			    &array_unproject_output);

  // Input image dimensions in pixels in x and y directions.
  size_t ii_size_x = imd->general->sample_count;
  size_t ii_size_y = imd->general->line_count;

  // Convenience aliases.
  meta_projection *ipb = imd->projection;
  project_parameters_t *ipp = &imd->projection->param;

  // First we march around the entire outside of the image and compute
  // projection coordinates for every pixel, keeping track of the
  // minimum and maximum projection coordinates in each dimension.
  // This lets us determine the exact extent of the DEM in
  // output projection coordinates.
  asfPrintStatus ("Determining input image extent in projection coordinate "
		  "space... ");

  double min_x = DBL_MAX;
  double max_x = -DBL_MAX;
  double min_y = DBL_MAX;
  double max_y = -DBL_MAX;

  // In going around the edge, we are just trying to determine the
  // extent of the image in the horizontal, so we don't care about
  // height yet.
  { // Scoping block.
    // Number of pixels in the edge of the image.
    size_t edge_point_count = 2 * ii_size_x + 2 * ii_size_y - 4;
    double *lats = g_new0 (double, edge_point_count);
    double *lons = g_new0 (double, edge_point_count);
    size_t current_edge_point = 0;
    size_t ii = 0, jj = 0;
    for ( ; ii < ii_size_x - 1 ; ii++ ) {
      return_code = get_pixel_lat_long (imd, unproject_input, ii, jj,
					&(lats[current_edge_point]),
					&(lons[current_edge_point]));
      g_assert (return_code);
      current_edge_point++;
    }
    for ( ; jj < ii_size_y - 1 ; jj++ ) {
      return_code = get_pixel_lat_long (imd, unproject_input, ii, jj,
					&(lats[current_edge_point]),
					&(lons[current_edge_point]));
      g_assert (return_code);
      current_edge_point++;
    }
    for ( ; ii > 0 ; ii-- ) {
      return_code = get_pixel_lat_long (imd, unproject_input, ii, jj,
					&(lats[current_edge_point]),
					&(lons[current_edge_point]));
      g_assert (return_code);
      current_edge_point++;
    }
    for ( ; jj > 0 ; jj-- ) {
      return_code = get_pixel_lat_long (imd, unproject_input, ii, jj,
					&(lats[current_edge_point]),
					&(lons[current_edge_point]));
      g_assert (return_code);
      current_edge_point++;
    }
    g_assert (current_edge_point == edge_point_count);
    // Pointers to arrays of projected coordinates to be filled in.
    // The projection function will allocate this memory itself.
    double *x = NULL, *y = NULL;
    // Project all the edge pixels.
    return_code = array_project_output (pp, lats, lons, NULL, &x, &y, NULL,
					edge_point_count, datum);
    g_assert (return_code == TRUE);
    // Find the extents of the image in projection coordinates.
    for ( ii = 0 ; ii < edge_point_count ; ii++ ) {
      if ( x[ii] < min_x ) { min_x = x[ii]; }
      if ( x[ii] > max_x ) { max_x = x[ii]; }
      if ( y[ii] < min_y ) { min_y = y[ii]; }
      if ( y[ii] > max_y ) { max_y = y[ii]; }
    }

    free (y);
    free (x);
    g_free (lons);
    g_free (lats);
  }

  asfPrintStatus ("done.\n\n");

  // Issue a warning when the chosen pixel size is smaller than the
  // input pixel size.  FIXME: this condition will really never fire
  // for pseudoprojected image, since the pixels size of the input is
  // tiny (degrees per pixel) and the pixel_size has already been
  // computed in asf_geocode function itself as an arc length on the
  // ground.
  if ( GSL_MIN(imd->general->x_pixel_size,
	       imd->general->y_pixel_size) > pixel_size ) {
    asfPrintWarning
      ("Requested pixel size %lf is smaller then the input image resolution "
       "(%le meters).\n", pixel_size,
       GSL_MIN (imd->general->x_pixel_size, imd->general->y_pixel_size));
  }

  // The pixel size requested by the user better not oversample by the
  // factor of 2.  Specifying --force will skip this check.  FIXME:
  // same essential problem as the above condition, but in this case
  // it always goes off.
  //  if (!force_flag && GSL_MIN(imd->general->x_pixel_size,
  //	       imd->general->y_pixel_size) > (2*pixel_size) ) {
  //    report_func
  //      ("Requested pixel size %lf is smaller then the minimum implied by half \n"
  //       "the input image resolution (%le meters), this is not supported.\n",
  //       pixel_size, GSL_MIN (imd->general->x_pixel_size,
  //			    imd->general->y_pixel_size));
  //  }

  asfPrintStatus ("Opening input DEM image... ");
  char *input_data_file = (char *) MALLOC(sizeof(char)*(strlen(input_image)+5));
  sprintf(input_data_file, "%s.img", input_image);
  FloatImage *iim
    = float_image_new_from_file (ii_size_x, ii_size_y, input_data_file, 0,
				 FLOAT_IMAGE_BYTE_ORDER_BIG_ENDIAN);
  FREE(input_data_file);
  asfPrintStatus ("done.\n\n");

  // Maximum pixel indicies in output image.
  size_t oix_max = ceil ((max_x - min_x) / pixel_size);
  size_t oiy_max = ceil ((max_y - min_y) / pixel_size);

  // Output image dimensions.
  size_t oi_size_x = oix_max + 1;
  size_t oi_size_y = oiy_max + 1;

  // Output image.
  FloatImage *oim = float_image_new (oi_size_x, oi_size_y);

  // Translate the command line notion of the resampling method into
  // the lingo known by the float_image class.  The compiler is
  // reassured with a default.
  float_image_sample_method_t float_image_sample_method
    = FLOAT_IMAGE_SAMPLE_METHOD_BILINEAR;
  switch ( resample_method ) {
  case RESAMPLE_NEAREST_NEIGHBOR:
    float_image_sample_method = FLOAT_IMAGE_SAMPLE_METHOD_NEAREST_NEIGHBOR;
    break;
  case RESAMPLE_BILINEAR:
    float_image_sample_method = FLOAT_IMAGE_SAMPLE_METHOD_BILINEAR;
    break;
  case RESAMPLE_BICUBIC:
    float_image_sample_method = FLOAT_IMAGE_SAMPLE_METHOD_BICUBIC;
    break;
  default:
    g_assert_not_reached ();
  }

  // We need to find the z coordinates in the output projection of all
  // the pixels in the input DEM.  We store these values in their own
  // FloatImage instance.

  //FloatImage *x_coords = float_image_new (ii_size_x, ii_size_y);
  //FloatImage *y_coords = float_image_new (ii_size_x, ii_size_y);
  FloatImage *z_coords = float_image_new (ii_size_x, ii_size_y);

  // We transform the points using the array transformation function
  // for efficiency, but we don't want to do them all at once, since
  // that would require huge gobs of memory.
  const size_t max_transform_chunk_pixels = 5000000;
  size_t rows_per_chunk = max_transform_chunk_pixels / ii_size_x;
  size_t chunk_pixels = rows_per_chunk * ii_size_x;
  double *chunk_x = g_new0 (double, chunk_pixels);
  double *chunk_y = g_new0 (double, chunk_pixels);
  double *chunk_z = g_new0 (double, chunk_pixels);
  double *lat = g_new0 (double, chunk_pixels);
  double *lon = g_new0 (double, chunk_pixels);
  double *height = g_new0 (double, chunk_pixels);

  asfPrintStatus ("Determining Z coordinates of input pixels in output "
		  "projection space... ");

  // Transform all the chunks, storing results in the z coordinate image.
  size_t ii, jj, kk;		// Index variables.
  for ( ii = 0 ; ii < ii_size_y ; ) {
    size_t rows_remaining = ii_size_y - ii;
    size_t rows_to_load
      = rows_per_chunk < rows_remaining ? rows_per_chunk : rows_remaining;
    for ( jj = 0 ; jj < rows_to_load ; jj++ ) {
      size_t current_image_row = ii + jj;
      for ( kk = 0 ; kk < ii_size_x ; kk++ ) {
	size_t current_chunk_pixel = jj * ii_size_x + kk;
	chunk_x[current_chunk_pixel] = ipb->startX + kk * ipb->perX;
	chunk_y[current_chunk_pixel]
	  = ipb->startY + current_image_row * ipb->perY;
	if ( imd->projection->type == LAT_LONG_PSEUDO_PROJECTION ) {
	  chunk_x[current_chunk_pixel] *= D2R;
	  chunk_y[current_chunk_pixel] *= D2R;
	}
	chunk_z[current_chunk_pixel]
	  = float_image_get_pixel (iim, kk, current_image_row);
      }
    }
    long current_chunk_pixels = rows_to_load * ii_size_x;
    array_unproject_input (ipp, chunk_x, chunk_y, chunk_z, &lat, &lon,
			   &height, current_chunk_pixels, ipb->datum);
    array_project_output (pp, lat, lon, height, &chunk_x, &chunk_y, &chunk_z,
			  current_chunk_pixels, datum);
    for ( jj = 0 ; jj < rows_to_load ; jj++ ) {
      size_t current_image_row = ii + jj;
      for ( kk = 0 ; kk < ii_size_x ; kk++ ) {
	size_t current_chunk_pixel = jj * ii_size_x + kk;
	// Current pixel x, y, z coordinates.
	//float cp_x = (float) chunk_x[current_chunk_pixel];
	//float cp_y = (float) chunk_y[current_chunk_pixel];
	float cp_z = (float) chunk_z[current_chunk_pixel];
	//float_image_set_pixel (x_coords, kk, current_image_row, cp_x);
	//float_image_set_pixel (y_coords, kk, current_image_row, cp_y);
	float_image_set_pixel (z_coords, kk, current_image_row, cp_z);
      }
    }

    ii += rows_to_load;
  }

  asfPrintStatus ("done.\n\n");

#ifdef DEBUG_GEOCODE_DEM_Z_COORDS_IMAGE_AS_JPEG
  // Take a look at the z_coordinate image (for debugging).
  float_image_export_as_jpeg_with_mask_interval (z_coords, "z_coords.jpg",
						 GSL_MAX (z_coords->size_x,
							  z_coords->size_y),
						 -FLT_MAX, -100);
#endif

  g_free (chunk_x);
  g_free (chunk_y);
  g_free (chunk_z);
  g_free (lat);
  g_free (lon);
  g_free (height);

  // Now we want to determine the pixel coordinates in the input which
  // correspond to each of the output pixels.  We can then sample the
  // new height value already computed for that input pixel to
  // determine the pixel value to use as output.

  // We want to proceed in chunks as we did when going in the other
  // direction.
  rows_per_chunk = max_transform_chunk_pixels / oi_size_x;
  chunk_pixels = rows_per_chunk * oi_size_x;
  chunk_x = g_new0 (double, chunk_pixels);
  chunk_y = g_new0 (double, chunk_pixels);
  // We don't have height information in this direction, nor do we care.
  chunk_z = NULL;
  lat = g_new0 (double, chunk_pixels);
  lon = g_new0 (double, chunk_pixels);
  // We don't have height information in this direction, nor do we care.
  height = NULL;

  asfPrintStatus ("Sampling Z coordinates to form pixels in output projection "
		  "space... ");

  // Transform all the chunks, using the results to form the output image.
  for ( ii = 0 ; ii < oi_size_y ; ) {
    size_t rows_remaining = oi_size_y - ii;
    size_t rows_to_load
      = rows_per_chunk < rows_remaining ? rows_per_chunk : rows_remaining;
    for ( jj = 0 ; jj < rows_to_load ; jj++ ) {
      size_t current_image_row = ii + jj;
      for ( kk = 0 ; kk < oi_size_x ; kk++ ) {
	size_t current_chunk_pixel = jj * oi_size_x + kk;
	chunk_x[current_chunk_pixel] = min_x + kk * pixel_size;
	chunk_y[current_chunk_pixel] = max_y - current_image_row * pixel_size;
      }
    }
    long current_chunk_pixels = rows_to_load * oi_size_x;
    array_unproject_output (pp, chunk_x, chunk_y, NULL, &lat, &lon, NULL,
			    current_chunk_pixels, datum);
    array_project_input (ipp, lat, lon, NULL, &chunk_x, &chunk_y, NULL,
			 current_chunk_pixels, ipb->datum);
    if ( imd->projection->type == LAT_LONG_PSEUDO_PROJECTION ) {
      ssize_t ll;     // For (semi)clarity we don't reuse index variable :)
      for ( ll = 0 ; ll < current_chunk_pixels ; ll++ ) {
	chunk_x[ll] *= R2D;
	chunk_y[ll] *= R2D;
      }
    }

    for ( jj = 0 ; jj < rows_to_load ; jj++ ) {
      size_t current_image_row = ii + jj;
      for ( kk = 0 ; kk < oi_size_x ; kk++ ) {
	size_t current_chunk_pixel = jj * oi_size_x + kk;

	// Compute pixel coordinates in input image.
	ssize_t in_x
	  = (chunk_x[current_chunk_pixel] - ipb->startX) / ipb->perX;
	ssize_t in_y
	  = (chunk_y[current_chunk_pixel] - ipb->startY) / ipb->perY;

	if ( in_image (z_coords, in_x, in_y) ) {
	  // FIXME: something needs to be done somewhere about
	  // propogating no data values.
	  float_image_set_pixel (oim, kk, current_image_row,
				 float_image_sample (z_coords, in_x, in_y,
						     resample_method));
	}
	else {
	  float_image_set_pixel (oim, kk, current_image_row, background_value);
	}
      }
    }

    ii += rows_to_load;
  }

  asfPrintStatus ("done.\n\n");

  g_free (chunk_x);
  g_free (chunk_y);
  g_free (lat);
  g_free (lon);

#ifdef DEBUG_GEOCODE_DEM_OUTPUT_IMAGE_AS_JPEG
  // Take a look at the output image (for debugging).
  float_image_export_as_jpeg_with_mask_interval (oim, "oim.jpg",
						 GSL_MAX (oim->size_x,
							  oim->size_y),
						 -FLT_MAX, -100);
#endif

  // Store the output image.
  asfPrintStatus ("Storing output image... ");
  char *output_data_file =
    (char *) MALLOC(sizeof(char)*(strlen(output_image)+5));
  sprintf(output_data_file, "%s.img", output_image);
  return_code = float_image_store (oim, output_data_file,
				   FLOAT_IMAGE_BYTE_ORDER_BIG_ENDIAN);
  g_assert (return_code == 0);
  asfPrintStatus ("done.\n\n");

  // Now we need some metadata for the output image.  We will just
  // start with the metadata from the input image and add the
  // geocoding parameters.

  char *input_meta_file = (char *) MALLOC(sizeof(char)*(strlen(input_image)+6));
  sprintf(input_meta_file, "%s.meta", input_image);

  char *output_meta_file =
    (char *) MALLOC(sizeof(char)*(strlen(output_image)+6));
  sprintf(output_meta_file, "%s.meta", output_image);

  meta_parameters *omd = meta_read (input_meta_file);

  // Adjust the metadata to correspond to the output image instead of
  // the input image.

  omd->general->x_pixel_size = pixel_size;
  omd->general->y_pixel_size = pixel_size;
  omd->general->line_count = oi_size_y;
  omd->general->sample_count = oi_size_x;

  // SAR block is not really appropriate for map projected images, but
  // since it ended up with this value that can signify map
  // projectedness in it somehow, we fill it in for safety.
  omd->sar->image_type = 'P';

  // Note that we have already verified that the input image is
  // projected, and since we initialize the output metadata from there
  // we know we will have a projection block.
  omd->projection->type = projection_type;
  omd->projection->startX = min_x;
  omd->projection->startY = max_y;
  omd->projection->perX = pixel_size;
  omd->projection->perY = -pixel_size;
  strcpy (omd->projection->units, "meters");

  // Set the spheroid axes lengths as appropriate for the output datum.
  spheroid_axes_lengths (datum_spheroid (datum), &(omd->projection->re_major),
			 &(omd->projection->re_minor));

  // What the heck, might as well set the ones in the general block as
  // well.
  spheroid_axes_lengths (datum_spheroid (datum), &(omd->general->re_major),
			 &(omd->general->re_minor));

  // Latitude and longitude at center of the output image.  We will
  // set these relative to the spheroid underlying the datum in use
  // for the projected image.  Yeah, that seems appropriate.
  double lat_0, lon_0;
  double center_x = omd->projection->startX + (omd->projection->perX
					       * omd->general->line_count / 2);
  double center_y = (omd->projection->startY
		     + (omd->projection->perY
			* omd->general->sample_count / 2));
  unproject_output (pp, center_x, center_y, ASF_PROJ_NO_HEIGHT, &lat_0, &lon_0,
		    NULL, datum);
  omd->general->center_latitude = R2D * lat_0;
  omd->general->center_longitude = R2D * lon_0;

  // FIXME: We are ignoring the meta_location fields for now since I'm
  // not sure whether they are supposed to refer to the corner pixels
  // or the corners of the data itself.

  if ( lat_0 > 0.0 ) {
    omd->projection->hem = 'N';
  }
  else {
    omd->projection->hem = 'S';
  }

  // Convert the projection parameter values back into degrees.
  to_degrees (projection_type, pp);
  omd->projection->param = *pp;
  meta_write (omd, output_meta_file);

  float_image_free (oim);
  FREE(output_data_file);
  meta_free (omd);
  FREE(input_meta_file);
  FREE(output_meta_file);

  return 0;
}

max_stars_repo_name: glshort/MapReady

max_issues_repo_name: glshort/MapReady

max_forks_repo_name: glshort/MapReady


text:
#include <stdio.h>
#include <gsl/gsl_matrix.h>
#include "bstrlib/bstrlib.h"
#include "paths.h"
#include "ParseSCF.h"
#include "HTightBinding.h"
#include "BandEnergy.h"
#include "SpinOrbit.h"

int main(int argc, char *argv[]) {
    if (argc < 8) {
        printf("SOC-induced anisotropy calculation -- invoke with:\n");
        printf("anisotropy.out (system_name) (na) (nb) (nc) (soc_strength) (theta1) (phi1) (theta2) (phi2)\n");
        return 1;
    }
    char *system_name = argv[1];
    int na = atoi(argv[2]);
    int nb = atoi(argv[3]);
    int nc = atoi(argv[4]);
    double soc_strength = atof(argv[5]);
    double theta1 = atof(argv[6]);
    double phi1 = atof(argv[7]);
    double theta2 = atof(argv[8]);
    double phi2 = atof(argv[9]);

    char *scf_path = cwannier_data_path(system_name, "wannier", "\0", "scf.out");
    double num_electrons, alat;
    gsl_matrix *R = gsl_matrix_alloc(3, 3);

    int err = ParseSCF(scf_path, &num_electrons, &alat, R);
    if (err != CWANNIER_PARSESCF_OK) {
        printf("Error code = %d returned from ParseSCF\n", err);
        return err;
    }
    bcstrfree(scf_path);

    char *hr_up_path = cwannier_data_path(system_name, "wannier", system_name, "_up_hr.dat");
    HTightBinding *Hrs_up = ExtractHTightBinding(hr_up_path);
    bcstrfree(hr_up_path);

    char *hr_dn_path = cwannier_data_path(system_name, "wannier", system_name, "_dn_hr.dat");
    HTightBinding *Hrs_dn = ExtractHTightBinding(hr_dn_path);
    bcstrfree(hr_dn_path);

    printf("Hamiltonian loaded.\n");

    HTightBinding *Hrs_soc_1 = HamiltonianWithSOC(soc_strength, theta1, phi1, Hrs_up, Hrs_dn);
    HTightBinding *Hrs_soc_2 = HamiltonianWithSOC(soc_strength, theta2, phi2, Hrs_up, Hrs_dn);

    bool use_cache = true;
    double E_Fermi_1 = 0.0;
    double E_Fermi_2 = 0.0;

    double energy1 = BandEnergy(&E_Fermi_1, Hrs_soc_1, R, num_electrons, na, nb, nc, use_cache);
    printf("Got energy1 = %f\n", energy1);

    double energy2 = BandEnergy(&E_Fermi_2, Hrs_soc_2, R, num_electrons, na, nb, nc, use_cache);
    printf("Got energy2 = %f\n", energy2);

    printf("energy1 - energy2 = %e\n", energy1 - energy2);

    FreeHTightBinding(Hrs_up);
    FreeHTightBinding(Hrs_dn);
    FreeHTightBinding(Hrs_soc_1);
    FreeHTightBinding(Hrs_soc_2);
    return 0;
}

max_stars_repo_name: tflovorn/cwannier

max_issues_repo_name: tflovorn/cwannier

max_forks_repo_name: tflovorn/cwannier


text:
/**
 *
 * @file core_zpemv.c
 *
 *  PLASMA core_blas kernel
 *  PLASMA is a software package provided by Univ. of Tennessee,
 *  Univ. of California Berkeley and Univ. of Colorado Denver
 *
 * @version 2.6.0
 * @author Dulceneia Becker
 * @date 2011-06-29
 * @precisions normal z -> c d s
 *
 **/
#include <cblas.h>
#include <lapacke.h>
#include "common.h"

/***************************************************************************//**
 *
 * @ingroup CORE_PLASMA_Complex64_t
 *
 *  CORE_zpemv  performs one of the matrix-vector operations
 *
 *     y = alpha*op( A )*x + beta*y
 *
 *  where  op( A ) is one of
 *
 *     op( A ) = A   or   op( A ) = A**T   or   op( A ) = A**H,
 *
 *  alpha and beta are scalars, x and y are vectors and A is a
 *  pentagonal matrix (see further details).
 *
 *
 *  Arguments
 *  ==========
 *
 * @param[in] storev
 *
 *         @arg PlasmaColumnwise :  array A stored columwise
 *         @arg PlasmaRowwise    :  array A stored rowwise
 *
 * @param[in] trans
 *
 *         @arg PlasmaNoTrans   :  y := alpha*A*x    + beta*y.
 *         @arg PlasmaTrans     :  y := alpha*A**T*x + beta*y.
 *         @arg PlasmaConjTrans :  y := alpha*A**H*x + beta*y.
 *
 * @param[in] M
 *         Number of rows of the matrix A.
 *         M must be at least zero.
 *
 * @param[in] N
 *         Number of columns of the matrix A.
 *         N must be at least zero.
 *
 * @param[in] L
 *         Order of triangle within the matrix A (L specifies the shape
 *         of the matrix A; see further details).
 *
 * @param[in] ALPHA
 *         Scalar alpha.
 *
 * @param[in] A
 *         Array of size LDA-by-N.  On entry, the leading M by N part
 *         of the array A must contain the matrix of coefficients.
 *
 * @param[in] LDA
 *         Leading dimension of array A.
 *
 * @param[in] X
 *         On entry, the incremented array X must contain the vector x.
 *
 * @param[in] INCX
 *         Increment for the elements of X. INCX must not be zero.
 *
 * @param[in] BETA
 *         Scalar beta.
 *
 * @param[in,out] Y
 *         On entry, the incremented array Y must contain the vector y.
 *
 * @param[out] INCY
 *         Increment for the elements of Y. INCY must not be zero.
 *
 * @param[out] WORK
 *         Workspace array of size at least L.
 *
 *  Further Details
 *  ===============
 *
 *               |     N    |
 *            _   ___________   _
 *               |          |
 *     A:        |          |
 *          M-L  |          |
 *               |          |  M
 *            _  |.....     |
 *               \    :     |
 *            L    \  :     |
 *            _      \:_____|  _
 *
 *               |  L | N-L |
 *
 *
 *******************************************************************************
 *
 * @return
 *          \retval PLASMA_SUCCESS successful exit
 *          \retval <0 if -i, the i-th argument had an illegal value
 *
 ******************************************************************************/

#if defined(PLASMA_HAVE_WEAK)
#pragma weak CORE_zpemv = PCORE_zpemv
#define CORE_zpemv PCORE_zpemv
#endif
int CORE_zpemv(PLASMA_enum trans, int storev,
               int M, int N, int L,
               PLASMA_Complex64_t ALPHA,
               const PLASMA_Complex64_t *A, int LDA,
               const PLASMA_Complex64_t *X, int INCX,
               PLASMA_Complex64_t BETA,
               PLASMA_Complex64_t *Y, int INCY,
               PLASMA_Complex64_t *WORK)
{

   /*
    *  y = alpha * op(A) * x + beta * y
    */

    int K;
    static PLASMA_Complex64_t zzero = 0.0;

    /* Check input arguments */
    if ((trans != PlasmaNoTrans) && (trans != PlasmaTrans) && (trans != PlasmaConjTrans)) {
        coreblas_error(1, "Illegal value of trans");
        return -1;
    }
    if ((storev != PlasmaColumnwise) && (storev != PlasmaRowwise)) {
        coreblas_error(2, "Illegal value of storev");
        return -2;
    }
    if (!( ((storev == PlasmaColumnwise) && (trans != PlasmaNoTrans)) ||
           ((storev == PlasmaRowwise) && (trans == PlasmaNoTrans)) )) {
        coreblas_error(2, "Illegal values of trans/storev");
        return -2;
    }
    if (M < 0) {
        coreblas_error(3, "Illegal value of M");
        return -3;
    }
    if (N < 0) {
        coreblas_error(4, "Illegal value of N");
        return -4;
    }
    if (L > min(M ,N)) {
        coreblas_error(5, "Illegal value of L");
        return -5;
    }
    if (LDA < max(1,M)) {
        coreblas_error(8, "Illegal value of LDA");
        return -8;
    }
    if (INCX < 1) {
        coreblas_error(10, "Illegal value of INCX");
        return -10;
    }
    if (INCY < 1) {
        coreblas_error(13, "Illegal value of INCY");
        return -13;
    }

    /* Quick return */
    if ((M == 0) || (N == 0))
        return PLASMA_SUCCESS;
    if ((ALPHA == zzero) && (BETA == zzero))
        return PLASMA_SUCCESS;

    /* If L < 2, there is no triangular part */
    if (L == 1) L = 0;

    /* Columnwise */
    if (storev == PlasmaColumnwise) {
        /*
         *        ______________
         *        |      |     |    A1: A[ 0 ]
         *        |      |     |    A2: A[ M-L ]
         *        |  A1  |     |    A3: A[ (N-L) * LDA ]
         *        |      |     |
         *        |______| A3  |
         *        \      |     |
         *          \ A2 |     |
         *            \  |     |
         *              \|_____|
         *
         */

        /* Columnwise / NoTrans */
        if (trans == PlasmaNoTrans) {
            coreblas_error(1, "The case PlasmaNoTrans / PlasmaColumnwise is not yet implemented");
            return -1;
        }
        /* Columnwise / [Conj]Trans */
        else {

            /* L top rows of y */
            if (L > 0) {

                /* w = A_2' * x_2 */
                cblas_zcopy(
                            L, &X[INCX*(M-L)], INCX, WORK, 1);
                cblas_ztrmv(
                            CblasColMajor, (CBLAS_UPLO)PlasmaUpper,
                            (CBLAS_TRANSPOSE)trans,
                            (CBLAS_DIAG)PlasmaNonUnit,
                            L, &A[M-L], LDA, WORK, 1);

                if (M > L) {

                    /* y_1 = beta * y_1 [ + alpha * A_1 * x_1 ] */
                    cblas_zgemv(
                                CblasColMajor, (CBLAS_TRANSPOSE)trans,
                                M-L, L, CBLAS_SADDR(ALPHA), A, LDA,
                                X, INCX, CBLAS_SADDR(BETA), Y, INCY);

                    /* y_1 = y_1 + alpha * w */
                    cblas_zaxpy(L, CBLAS_SADDR(ALPHA), WORK, 1, Y, INCY);

                } else {

                    /* y_1 = y_1 + alpha * w */
                    if (BETA == zzero) {
                        cblas_zscal(L, CBLAS_SADDR(ALPHA), WORK, 1);
                        cblas_zcopy(L, WORK, 1, Y, INCY);
                    } else {
                        cblas_zscal(L, CBLAS_SADDR(BETA), Y, INCY);
                        cblas_zaxpy(L, CBLAS_SADDR(ALPHA), WORK, 1, Y, INCY);
                    }
                }
            }

            /* N-L bottom rows of Y */
            if (N > L) {
                K = N - L;
                cblas_zgemv(
                            CblasColMajor, (CBLAS_TRANSPOSE)trans,
                            M, K, CBLAS_SADDR(ALPHA), &A[LDA*L], LDA,
                            X, INCX, CBLAS_SADDR(BETA), &Y[INCY*L], INCY);
            }
        }
    }
    /* Rowwise */
    else {
        /*
         * --------------
         * |            | \           A1:  A[ 0 ]
         * |    A1      |   \         A2:  A[ (N-L) * LDA ]
         * |            | A2  \       A3:  A[ L ]
         * |--------------------\
         * |        A3          |
         * ----------------------
         *
         */

        /* Rowwise / NoTrans */
        if (trans == PlasmaNoTrans) {
            /* L top rows of A and y */
            if (L > 0) {

                /* w = A_2 * x_2 */
                cblas_zcopy(
                            L, &X[INCX*(N-L)], INCX, WORK, 1);
                cblas_ztrmv(
                            CblasColMajor, (CBLAS_UPLO)PlasmaLower,
                            (CBLAS_TRANSPOSE)PlasmaNoTrans,
                            (CBLAS_DIAG)PlasmaNonUnit,
                            L, &A[LDA*(N-L)], LDA, WORK, 1);

                if (N > L) {

                    /* y_1 = beta * y_1 [ + alpha * A_1 * x_1 ] */
                    cblas_zgemv(
                                CblasColMajor, (CBLAS_TRANSPOSE)PlasmaNoTrans,
                                L, N-L, CBLAS_SADDR(ALPHA), A, LDA,
                                X, INCX, CBLAS_SADDR(BETA), Y, INCY);

                    /* y_1 = y_1 + alpha * w */
                    cblas_zaxpy(L, CBLAS_SADDR(ALPHA), WORK, 1, Y, INCY);

                } else {

                    /* y_1 = y_1 + alpha * w */
                    if (BETA == zzero) {
                        cblas_zscal(L, CBLAS_SADDR(ALPHA), WORK, 1);
                        cblas_zcopy(L, WORK, 1, Y, INCY);
                    } else {
                        cblas_zscal(L, CBLAS_SADDR(BETA), Y, INCY);
                        cblas_zaxpy(L, CBLAS_SADDR(ALPHA), WORK, 1, Y, INCY);
                    }
                }
            }

            /* M-L bottom rows of Y */
            if (M > L) {
                cblas_zgemv(
                        CblasColMajor, (CBLAS_TRANSPOSE)PlasmaNoTrans,
                        M-L, N, CBLAS_SADDR(ALPHA), &A[L], LDA,
                        X, INCX, CBLAS_SADDR(BETA), &Y[INCY*L], INCY);
            }
        }
        /* Rowwise / [Conj]Trans */
        else {
            coreblas_error(1, "The case Plasma[Conj]Trans / PlasmaRowwise is not yet implemented");
            return -1;
        }
    }

    return PLASMA_SUCCESS;
}