text:
/* blas/blas.c
 *
 * Copyright (C) 1996, 1997, 1998, 1999, 2000, 2001 Gerard Jungman & 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 2 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.
 */

/* GSL implementation of BLAS operations for vectors and dense
 * matrices.  Note that GSL native storage is row-major.  */

#include <gsl/gsl_math.h>
#include <gsl/gsl_errno.h>
#include <gsl/blas/gsl_blas_types.h>
#include <gsl/blas/gsl_blas.h>
#include <gsl/blas/cblas.h>

/* ========================================================================
 * Level 1
 * ========================================================================
 */

/* CBLAS defines vector sizes in terms of int. GSL defines sizes in
   terms of size_t, so we need to convert these into integers.  There
   is the possibility of overflow here. FIXME: Maybe this could be
   caught */

#define INT(X) ((int)(X))

int gsl_blas_ddot (const gsl_vector * X, const gsl_vector * Y, double *result){
  if (X->size == Y->size){
	  *result = ddot (&INT (X->size), X->data, &INT (X->stride), Y->data, &INT (Y->stride));
      return GSL_SUCCESS;
    }
  else {
      GSL_ERROR ("invalid length", GSL_EBADLEN);
    }
}

/**
*  DNRM2 returns the euclidean norm of a vector via the function
*  name, so that
*
*     DNRM2 := sqrt( x'*x )
*/

double gsl_blas_dnrm2 (const gsl_vector * X){
  return dnrm2 (&INT (X->size), X->data, &INT (X->stride));
}

/**
*     interchanges two vectors.
*/
int gsl_blas_dswap (gsl_vector * X, gsl_vector * Y){
  if (X->size == Y->size){
      dswap (&INT (X->size), X->data, &INT (X->stride), Y->data, &INT (Y->stride));
      return GSL_SUCCESS;
    }
  else{
      GSL_ERROR ("invalid length", GSL_EBADLEN);
    };
}

/**
*     copies a vector, x, to a vector, y.
*/
int gsl_blas_dcopy (const gsl_vector * X, gsl_vector * Y){
  if (X->size == Y->size){
      dcopy (&INT (X->size), X->data, &INT (X->stride), Y->data, &INT (Y->stride));
      return GSL_SUCCESS;
    }
  else {
      GSL_ERROR ("invalid length", GSL_EBADLEN);
    }
}

/* ===========================================================================
 * Level 2
 * ===========================================================================
 */

int
gsl_blas_dgemv (CBLAS_TRANSPOSE_t TransA, double alpha, const gsl_matrix * A,
                const gsl_vector * X, double beta, gsl_vector * Y)
{
	//here the matrix is by default transposed, so to get the right behaviour we will transpose it when the user didn't want it tranposed, and leave it the way it is when it should be transposed

  //note that M and N are changed to show that the matrix is transposed by default.
  const size_t M = A->size1;
  const size_t N = A->size2;

  char aTransp = (TransA == CblasNoTrans) ? 'T' : 'N';

  if ((TransA == CblasNoTrans && N == X->size && M == Y->size)
      || (TransA == CblasTrans && M == X->size && N == Y->size))
    {
      dgemv (&aTransp, &INT (N), &INT (M), &alpha, A->data,
                   &INT (A->tda), X->data, &INT (X->stride), &beta, Y->data,
                   &INT (Y->stride));
      return GSL_SUCCESS;
    }
  else
    {
      GSL_ERROR ("invalid length", GSL_EBADLEN);
    }
}

/*
 * ===========================================================================
 * Prototypes for level 3 BLAS
 * ===========================================================================
 */

int gsl_blas_dgemm (CBLAS_TRANSPOSE_t TransA, CBLAS_TRANSPOSE_t TransB,
                double alpha, const gsl_matrix * A, const gsl_matrix * B,
                double beta, gsl_matrix * C)
{
/**
	Since BLAS uses column major matrices, when asked to perform C = A*B we will
	actually compute C = B*A, with the number of rows and cols switched up.
	(C' = B' * A' if C = A * B, and we're working with transposes by default)

	That's why we need to switch the order of the matrices as well as the order of M and N.

*/

  const size_t M = C->size1;
  const size_t N = C->size2;
  const size_t MA = (TransA == CblasNoTrans) ? A->size1 : A->size2;
  const size_t NA = (TransA == CblasNoTrans) ? A->size2 : A->size1;
  const size_t MB = (TransB == CblasNoTrans) ? B->size1 : B->size2;
  const size_t NB = (TransB == CblasNoTrans) ? B->size2 : B->size1;

  char aTransp = (TransA == CblasNoTrans) ? 'N' : 'T';
  char bTransp = (TransB == CblasNoTrans) ? 'N' : 'T';

  if (M == MA && N == NB && NA == MB)   /* [MxN] = [MAxNA][MBxNB] */
    {
	  dgemm (&bTransp, &aTransp, &INT (N), &INT (M), &INT (NA),
                   &alpha, B->data, &INT (B->tda), A->data, &INT (A->tda), &beta,
                   C->data, &INT (C->tda));
      return GSL_SUCCESS;
    }
  else
    {
      GSL_ERROR ("invalid length", GSL_EBADLEN);
    }
}

max_stars_repo_name: MontyThibault/centre-of-mass-awareness

max_issues_repo_name: MontyThibault/centre-of-mass-awareness

max_forks_repo_name: MontyThibault/centre-of-mass-awareness


text:
#include <algorithm>
#include <cstdio>
#include <random>
#include <gsl/gsl_multimin.h>

#include "dgp.h"
#include "vec.h"
#include "fobj_sph.h"

typedef struct
{
   double tau;
   double lambda;
   double rho;
   dgp_t *G;
} gsl_params_t;

void gsl_fdf(const gsl_vector *x, void *params, double *f, gsl_vector *g)
{
   gsl_params_t *p = (gsl_params_t *)params;
   *f = fobj_sph(p->lambda, p->tau, *(p->G), x->data, g->data);
}

void gsl_df(const gsl_vector *x, void *params, gsl_vector *g)
{
   gsl_params_t *p = (gsl_params_t *)params;
   fobj_sph(p->lambda, p->tau, *(p->G), x->data, g->data);
}

double gsl_f(const gsl_vector *x, void *params)
{
   gsl_params_t *p = (gsl_params_t *)params;
   return fobj_sph(p->lambda, p->tau, *(p->G), x->data, NULL);
}

double init_tau(dgp_t &G)
{
   std::vector<double> tau_all(G.m_nedges);

   for (int i = 0; i < G.m_nedges; i++)
   {
      tau_all[i] = (G.m_l[i] + G.m_u[i]) / 2.0;
   }

   std::sort(tau_all.begin(), tau_all.end());

   return tau_all[(int)(0.75 * G.m_nedges)];;
}

double gsl_vector_norm(gsl_vector *v)
{
   double v_nrm = 0.0;
   for (size_t i = 0; i < v->size; i++)
   {
      v_nrm += v->data[i] * v->data[i];
   }
   return sqrt(v_nrm);
}

void init_sol(dgp_t &G, double *x)
{
   double dij, lij, uij, *xi, *xj;
   double y[3], y_nrm;
   int *neighs, nneighs, j;
   std::vector<int> s(G.m_nnodes);
   std::vector<bool> b(G.m_nnodes, false);

   // s = [0, 1, 2, ..., nnodes-1]
   for (int k = 0; k < G.m_nnodes; k++)
   {
      s[k] = k;
   }
   std::random_shuffle(s.begin(), s.end());

   // x = zeros(3 * nnodes)
   for (int k = 0; k < (3 * G.m_nnodes); ++k)
   {
      x[k] = 0.0;
   }

   std::default_random_engine rndEngine;
   std::normal_distribution<double> gaussian(0.0, 1.0);

   for (int k = 0; k < G.m_nnodes; k++)
   {
      int i = s[k];
      xi = &(x[3 * i]);
      G.neighs(i, nneighs, &neighs);
      for (int jj = 0; jj < nneighs; jj++)
      {
         j = neighs[jj];
         if (b[j])
         {
            continue;
         }
         G.vals(i, j, NULL, &lij, &uij);
         dij = (lij + uij) / 2.0;
         xj = &(x[3 * j]);
         // y is randomly distributed inside the sphere((0,0,0),1)
         vec_set(y, gaussian(rndEngine), gaussian(rndEngine), gaussian(rndEngine));
         y_nrm = vec_norm(y);
         // xj = dij * (y / y_nrm) + xi  (randomly distributed over sphere(xi, dij))
         vec_axpby(dij / y_nrm, y, 1.0, xi, xj);
      }
   }
}

void sph_print_iter(int i, gsl_params_t &p, double f, double fs, gsl_vector *g, int k)
{
   if ((i % 10) == 0)
   {
      printf(" iter |    tau   |    f     |    fs    |  nrm(g)  | #f_calls \n");
   }
   printf("%5d | %3.2e | %3.2e | %3.2e | %3.2e | %5d\n", i, p.tau, f, fs, gsl_vector_norm(g), k);
}

bool sph(dgp_t &G, double ftol, double &f, double *x, bool verbose)
{
   gsl_params_t gsl_params;
   gsl_params.lambda = 0.5;
   gsl_params.rho = 0.99;
   gsl_params.tau = init_tau(G);
   gsl_params.G = &G;

   gsl_vector *y = gsl_vector_alloc(3 * G.m_nnodes);
   gsl_vector *g = gsl_vector_alloc(3 * G.m_nnodes);

   init_sol(G, y->data);

   const gsl_multimin_fdfminimizer_type *fdfmin_method = gsl_multimin_fdfminimizer_vector_bfgs2;
   gsl_multimin_fdfminimizer *fdfmin;

   gsl_multimin_function_fdf gsl_fobj;
   gsl_fobj.n = 3 * G.m_nnodes;
   gsl_fobj.f = gsl_f;
   gsl_fobj.df = gsl_df;
   gsl_fobj.fdf = gsl_fdf;
   gsl_fobj.params = (void *)(&gsl_params);

   double fs;
   bool solved = false;
   int maxit = 1000, k = 0, status;
   for (int i = 0; i < maxit; i++)
   {
      fs = fobj_sph(gsl_params.lambda, gsl_params.tau, G, y->data, g->data);
      f = fobj_sph(gsl_params.lambda, 1E-16, G, y->data, NULL);

      if (verbose)
      {
         sph_print_iter(i, gsl_params, f, fs, g, k);
      }

      if (f < ftol)
      {
         solved = true;
         break;
      }

      gsl_params.tau *= gsl_params.rho;
      fdfmin = gsl_multimin_fdfminimizer_alloc(fdfmin_method, 3 * G.m_nnodes);
      gsl_multimin_fdfminimizer_set(fdfmin, &gsl_fobj, y, 0.01, 1E-3);

      // call local optimization
      for (k = 0, status = GSL_CONTINUE; k < 1000 && status == GSL_CONTINUE; k++)
      {
         status = gsl_multimin_fdfminimizer_iterate(fdfmin);
         if (status)
         {
            break;
         }
         status = gsl_multimin_test_gradient(fdfmin->gradient, 1E-8);
      }

      // copy from fdfmin to y
      for (int j = 0; j < 3 * G.m_nnodes; j++)
      {
         y->data[j] = fdfmin->x->data[j];
      }

      gsl_multimin_fdfminimizer_free(fdfmin);
   }

   // copy solution from y to x
   for (int i = 0; i < 3 * G.m_nnodes; i++)
   {
      x[i] = y->data[i];
   }

   gsl_vector_free(y);

   return solved;
}

ext: h

max_stars_repo_name: michaelsouza/sph

max_issues_repo_name: michaelsouza/sph

max_forks_repo_name: michaelsouza/sph


text:
/* ---------------------------------------------------------------- *
 * hod.h						            *
 * Martin Kilbinger, Henry J. McCracken, Jean Coupon 2008-2013      *
 * ---------------------------------------------------------------- */

#ifndef __HOD_H
#define __HOD_H

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

#include <gsl/gsl_sf_erf.h>
#include <gsl/gsl_sf_gamma.h>
#include <gsl/gsl_integration.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_roots.h>

#include "io.h"
#include "errorlist.h"
#include "config.h"
#include "maths.h"

#include "cosmo.h"
#include "nofz.h"
#include "halomodel.h"

/* errors */
#define HOD_BASE        -50000
#define HOD_MSTAR_MAX  HOD_BASE  + 1
#define HOD_NGAL_FIT   HOD_BASE  + 2

#define GM_base      -400000
#define GM_nz_type   GM_base + 1
#define GM_pi_max    GM_base + 2
#define GM_SHMR      GM_base + 3
#define GM_SMF       GM_base + 4
#define GM_GGL       GM_base + 5
#define GM_MSTAR_MAX GM_base + 6
#define GM_ETA       GM_base + 7

/* Limits for xi(r), 1- and 2-halo, in Mpc.h */
#define RMIN1 0.001
#define RMAX1 5.0
#define RMIN2 0.1
#define RMAX2 400.0

#define RMIN  RMIN1
#define RMAX  RMAX2
#define NR    50

#define MAXCHAR 1024
#define NLINES 100

#define NFIELD 100
#define NCHAR 20

/* GG: galaxy-galaxy stuff. GM: galaxy-dar matter stuff (lensing) */
#define GG 0
#define GM 1

#define getDoubleValue(array,col)  atof(array+NCHAR*(col-1))
#define getIntValue(array,col)     atoi(array+NCHAR*(col-1))
#define getCharValue(array,col)    array+NCHAR*(col-1)
#define getLine(array,i)           array+NFIELD*NCHAR*i

#define Nhalodata_t 4
typedef enum {w_of_theta, wp_rp, deltaSigma, smf} halodata_t;
#define shalodata_t(i) ( \
  i==w_of_theta    ? "w_of_theta" : \
  i==wp_rp         ? "wp_rp"      : \
  i==deltaSigma    ? "deltaSigma" : \
  i==smf           ? "smf"        : \
  "")

#define Nhalomode_t 3
typedef enum {galcorr_var, galcorr_cov, galcorr_log} halomode_t;
#define shalomode_t(i) (   \
 i==galcorr_var ? "galcorr_var" : \
 i==galcorr_cov ? "galcorr_cov" : \
 i==galcorr_log ? "galcorr_log" : \
 "")

#define Nngal_fit_t 5
typedef enum {ngal_log_fit, ngal_lin_fit, ngal_no_fit, ngal_lin_fit_only} ngal_fit_t;
#define sngal_fit_t(i) ( \
  i==ngal_log_fit  ? "ngal_log_fit"  : \
  i==ngal_lin_fit  ? "ngal_lin_fit"  : \
  i==ngal_no_fit   ? "ngal_no_fit"   : \
  i==ngal_lin_fit_only ? "ngal_lin_fit_only" : \
  "")

#define Nintconst_t 2
typedef enum {constant, random_file} intconst_t;
#define sintconst_t(i) ( \
  i==constant  ?    "constant"  : \
  i==random_file  ? "random_file"  : \
  "")

typedef struct {
  char WTHETA[500];
  char COVNAME[500];
  char DNDZ[500];
  char OUTFILE[500];
  int nbins;
  double ngal, ngalerr, area_deg, intconst, delta;
  double alpha_min,alpha_max,dalpha;
  double M1_min,M1_max,dM1;
  double Mmin_min,Mmin_max,dMmin;
  double Ngal_min,Ngal_max,dNgal;
  ngal_fit_t ngal_fit_type;
} config_hm_info;

/* ---------------------------------------------------------------- *
 * Global variables and functions                                   *
 * ---------------------------------------------------------------- */

#define ODD 0
#define EVEN 1

#define ABS(a) ((a) < 0 ? -(a) : (a))
#define PARITY(a) (a)%2 ? ODD : EVEN
#define FFTLog_SWAP(a,b) {FFTLog_TMP = (a); (a) = (b); (b) = FFTLog_TMP;}

#define CHANGE(fct) int change_##fct(cosmo_hm*, cosmo_hm*)
CHANGE(w_of_theta);
CHANGE(Pthg);
CHANGE(HOD);
CHANGE(ngd);
CHANGE(vc);
#undef CHANGE

/* ---------------------------------------------------------------- *
 * New types                                                        *
 * ---------------------------------------------------------------- */

/* This structure contains w_theta input data */
typedef struct {
  int nbins, nbins_RR;
  double *th;
  double *w;
  double *werr;   /* Error bars */
  double *wcov;   /* Covariance matrix */
  double *th_RR;  /* theta for random pairs (can be larger than for the data [even recommended]) */
  double *RR;     /* Random pairs (for the integral constraint) */
  double ngal;    /* Galaxy number density*/
  double ngal_err;
} wt_t;

typedef struct Nz_hist
{
  int nbins;
  double *z;
  double *n;
  double zmin;
  double zmax;
  double mean;
  double dz;
  double Norm;
  double max;
} Nz_hist;

/* ---------------------------------------------------------------- *
 * log-likelihood                                                   *
 * ---------------------------------------------------------------- */

double chi2_hm(cosmo_hm *model, halodata_t halodata, halomode_t halomode, const wt_t* data, ngal_fit_t ngal_fit_type,
	       double ngal, double ngalerr, intconst_t intconst_type, error **err);

/* ---------------------------------------------------------------- *
 * w(theta) and wp(rp)                                              *
 * ---------------------------------------------------------------- */

double *woftheta(cosmo_hm *model, pofk_t pofk, double *theta, int Ntheta, int i_bin, int j_bin, error **err);
double *woftheta_FFTLog_total(cosmo_hm *model, double ln_theta_min, double ln_delta_theta, int Ntheta,
			      double **theta, int i_bin, int j_bin, error **err);
double *wp(cosmo_hm *model, pofk_t pofk, const double *rp, int Nrp, double pi_max, int type, error **err);
double int_for_wp(double logr, void *params, error **err);

/* ---------------------------------------------------------------- *
 * Lensing functions (only for leauthaud11 model)                   *
 * ---------------------------------------------------------------- */

double *DeltaSigma(cosmo_hm *model, pofk_t pofk, const double *r, int N, error **err);
double int_for_Sigma(double logr, void *params, error **err);

/* ---------------------------------------------------------------- *
 * Stellar mass function (only for leauthaud11 model)               *
 * ---------------------------------------------------------------- */

double *dNdlogM10stellar(cosmo_hm *model, double *log10Mstellar, int N, error **err);
double *dNdlogM10stellar_c(cosmo_hm *model, double *log10Mstellar, int N, error **err);
double *dNdlogM10stellar_s(cosmo_hm *model, double *log10Mstellar, int N, error **err);

/* ---------------------------------------------------------------- *
 * HOD P(k) and xi(r)                                               *
 * ---------------------------------------------------------------- */

double *xiofr(cosmo_hm *model, pofk_t pofk, const double *r, int N, int type, error **err);
double *xi_1hcs(cosmo_hm *model,double a, const double *r, int N, int type, error **err);
double int_for_xi_1hcs(double logM, void *params, error **err);
double *xi_1hss(cosmo_hm *model,double a, const double *r, int N, int type, error **err);
double P1hss(double k, void *params);
double int_for_P1hss(double logM, void *params, error **err);
double concentration_sat(cosmo_hm *model, double Mh, double a, error **err);
double *xi_2h(cosmo_hm *model,double a, const double *r, int N, int type, error **err);
double *xi_P_NL(cosmo_hm *model, double a, const double *r, int N, error **err);
double FFTLog_P_NL(double k, void *params);
double P2h(double k, void *params);
double int_for_P2h(double logM, void *params, error **err);
double int_for_P2h_dm(double logM, void *params, error **err);
double bias_r(double M, void *params);

/* ---------------------------------------------------------------- *
 * HOD model functions                                              *
 * ---------------------------------------------------------------- */

double log10_fSHMR(double log10Mh, cosmo_hm *model);
double log10_fSHMR_inv_minus_x(double log10Mstar, void *p);
double log10_fSHMR_inv(double log10Mstar, cosmo_hm *model);
double Ngal_c(cosmo_hm *model, double Mh, double Mstellar_min, double Mstellar_max);
double int_for_Ngal_c(double Mstellar, void *params, error **err);
double phi_c_Mstellar(cosmo_hm *model, double log10Mstellar, double Mh, error **err);
double sigma_log_M(cosmo_hm *model, double log10Mstellar, error **err);
double eta_cen(cosmo_hm *model, double Mh, error **err);
double Ngal_s(cosmo_hm *model, double M, double Mstellar_min, double Mstellar_max);
double Ngal(cosmo_hm *model, double M, double Mstellar_min, double Mstellar_max);

double av_Mh_given_Mstar(cosmo_hm *model, double Mstellar, double a, error **err);
double int_for_phi_c_Mh(double log10Mh, void *params, error **err);
double int_for_phi_c_Mh_norm(double log10Mh, void *params, error **err);

/* ---------------------------------------------------------------- *
 * Deduced quantities                                               *
 * ---------------------------------------------------------------- */
double Mstar_tot_c(cosmo_hm *model, double Mh, double Mstellar_min, double Mstellar_max, error **err);
double int_for_Mstar_tot_c(double logMstar, void *params, error **err);
double Mstar_tot_s(cosmo_hm *model, double Mh, double Mstellar_min, double Mstellar_max, error **err);
double int_for_Mstar_tot_s(double logMstar, void *params, error **err);
double av_gal_bias(cosmo_hm *model, double a, error **err);
double int_for_av_gal_bias(double logM, void *intpar, error **err);
double av_halo_mass(cosmo_hm *model, double a, error **err);
double int_for_av_halo_mass(double logM, void *intpar, error **err);
double mass_weighted_av_stellar_mass(cosmo_hm *model, double a, error **err);
double int_for_mass_weighted_av_stellar_mass(double logM, void *params, error **err);
double int_for_mass_denum_weighted_av_stellar_mass(double logM, void *params, error **err);
double av_stellar_mass(cosmo_hm *model, double a, error **err);
double int_for_av_stellar_mass(double logM, void *params, error **err);
double av_frsat(cosmo_hm *model, double a, error **err);
double int_for_av_frsat(double logM, void *intpar, error **err);
double av_halo_bias(cosmo_hm *model, double a, error **err);
double int_for_av_halo_bias(double logM, void *intpar, error **err);
double dn_dlnM_integrand(double logM, void *intpar, error **err);

/* ---------------------------------------------------------------- *
 * Physical quantities                                              *
 * ---------------------------------------------------------------- */

double vc(cosmo_hm *model,double zmin, double zmax, error **err);
double int_for_vc(double z, void *params, error **err);
double logM_lim(cosmo_hm *model, double a, double r, error **err);
double ngal_triax(cosmo_hm *model, double a, double r,  error **err);
double ngal_den(cosmo_hm *model, double a, double logMlim, double Mstellar_min, double Mstellar_max, error **err);
double int_for_ngal_den(double logM, void *params,  error **err);
double ngal_den_c(cosmo_hm *model, double a, double logMlim, double Mstellar_min, double Mstellar_max, error **err);
double int_for_ngal_den_c(double logM, void *params,  error **err);
double ngal_den_s(cosmo_hm *model, double a, double logMlim, double Mstellar_min, double Mstellar_max, error **err);
double int_for_ngal_den_s(double logM, void *params,  error **err);
double ngal_den_vol(cosmo_hm *model, error **err);
double ngal_weighted(cosmo_hm *model, error **err);

/* ---------------------------------------------------------------- *
 * FFTLog                                                           *
 * ---------------------------------------------------------------- */

double xi_from_Pkr(gsl_function *ak, double r_prime, FFTLog_config *fc, error **err);
void FFTLog(FFTLog_config *fc, const gsl_function *ar_in,double *k, double *ak_out, int dir, error **err);
FFTLog_config *FFTLog_init(int N, double min, double max, double q, double mu);
void FFTLog_free(FFTLog_config *fc);
FFTLog_complex FFTLog_U_mu(double mu, FFTLog_complex z);

/* ---------------------------------------------------------------- *
 * Utils                                                            *
 * ---------------------------------------------------------------- */

void updateFrom_hm(cosmo_hm* avant, cosmo_hm* apres, error **err);
int change_vc(cosmo_hm *avant, cosmo_hm *apres);
int change_ngd(cosmo_hm *avant, cosmo_hm *apres);
int change_HOD(cosmo_hm *avant, cosmo_hm *apres);
int change_Pthg(cosmo_hm* avant, cosmo_hm* apres);
int change_w_of_theta(cosmo_hm *avant, cosmo_hm *apres);
wt_t *init_wt_t(double nlines, int wcov, error **err);
void free_wt_t(wt_t **model);
wt_t *read_wtheta(const char *data_file_name, intconst_t intconst_type, double delta, double intconst, const char *ran_file_name, error **err);
int getStrings(char *line, char *strings, char *delimit, size_t *N, error **err);
void read_config_hm_file(config_hm_info *config, char cname[], error **err);
Nz_hist *get_nz(redshift_t *nz_mk, int n_bin, error **err);
void  free_nz( Nz_hist *nz);
double normalize(double *data, size_t Nbins, double binsize,error **err);

/* ---------------------------------------------------------------- *
 * Obsolete stuff                                                   *
 * ---------------------------------------------------------------- */

double *woftheta_total_OBSOLETE(cosmo_hm *model, double log_theta_min, double log_delta_theta, int Ntheta, error **err);
nz_t *read_dndz_OBSOLETE(char *infile,error **err);
double Ngal_mean(cosmo_hm *model, double a, error **err);

double xir_mwolk(cosmo_hm *model, double a, double r, error **err);
wt_t *read_wtheta_mwolk(const char *name, double delta, double intconst, error **err);
double chi2_mwolk(cosmo_hm *model, const wt_t* data, ngal_fit_t ngal_fit_type,
		  double ngal, double ngalerr, double area, error **err);
double wp_mwolk(cosmo_hm *model, double rp, error **err);
double compute_chisq_wp(cosmo_hm *model, const wt_t *wth, double ngd_obs, double ngd_err,
			ngal_fit_t ngal_fit_type, double *ngd, int dologw, error **err);

#endif

max_stars_repo_name: danielgruen/ccv

max_issues_repo_name: danielgruen/ccv

max_forks_repo_name: danielgruen/ccv


text:
#include <iostream>
#include <vector>
#include <stdio.h>
#include <gsl\gsl_randist.h>
#include <gsl\gsl_rng.h>
#include <time.h>
#include <SFML\Graphics.hpp>
#ifndef _HEART_H
#define _HEART_H

class heart
{
	private:
		int L;  //sytem size length for now square

		std::vector<std::vector<int>> state;  //value indicates whether in excited state refactory or resting
		std::vector<std::vector<int>> new_state;  //using new_state and state to distinguish newly excited cells from previously excited
		std::vector<std::vector<int>> functional; //1 if functional 0 if disfunctional
		std::vector<std::vector<int>> bond; //1 if bond to site above current cell, 0 if no bond above site
		std::vector<std::vector<int>> unablated; // Defaults to zero initial value
		int heart_time; //internal clock
		//model parameters
		int tau; //refactory period
		int T; // pacemaker rate
		double nu, epsilom, delta; //probability vertical connection exists, probability of misfire given faulty and probability cell faulty

		//random number generator
		const gsl_rng_type*Type;
		gsl_rng*r;

		int total_bonds;
		int total_functional;
		int total_excited;

	public:

		heart()
		{

		}
		heart(const int _L,  int _T, int _tau ,double _nu, double _epsilom, double _delta)
		{
			L=_L;
			nu= _nu;
			tau=_tau;
			T=_T;
			epsilom=_epsilom;
			delta=_delta;
			heart_time = T;

			state.resize(L);
			new_state.resize(L);
			functional.resize(L);
			bond.resize(L);
			unablated.resize(L);
			//pacemaker action
			for(unsigned i=0; i<L; i++)
			{
				state[i].resize(L,0);  // all cells eligible for excitation
				new_state[i].resize(L,0);
				functional[i].resize(L,0);// none disfunctional at the moment
				bond[i].resize(L, 1); // all transverse bonds present at the moment
				unablated[i].resize(L,1);
			}
		//	total_functional = L*L;
		//	total_bonds = L*L;
			total_excited =0;//could be useful avoiding calling size of some vector

			 //HAVE NOT USED gsl_rng in class before need to make sure is working correctly

			//  environment variable GSL_RNG_TYPE
			gsl_rng_env_setup();

			Type = gsl_rng_mt19937;//random number generator using the Mersenne Twister algorithm
			r = gsl_rng_alloc (Type);

			gsl_rng_set(r,time(NULL));//random number generator seeded with time
			//gsl_rng_set(r,0);//random number generator seeded with time

			//assigning cells as having a bond above them or not and functional or not functional
			for(int i=0; i<L; i++)
			{
				for(int j=0; j<L; j++)
				{

					bond[i][j] = (gsl_rng_uniform(r)<nu);
					functional[i][j] = (gsl_rng_uniform(r)>delta);// inequality reversed
					//double x =  gsl_rng_uniform(r);
				//	std::cout <<x<< '\t' << (x<nu) << std::endl;

				}
			}

		}
		void operate()//one time step
		{
			total_excited =0;

			//Relax excited cells
				for(int i=0; i<L ;i++)
			{
				for(int j=0; j<L; j++)
				{
					if(state[i][j]==1)//searching 4 neighbours of excited cell trying to excite if resting and if don't misfire and if there is a bond
					{
						new_state[((i+1)<L)*(i+1)][j]+=(state[((i+1)<L)*(i+1)][j]<1)*bond[i][j]*(1-(1-functional[((i+1)<L)*(i+1)][j])*(gsl_rng_uniform(r)<epsilom))*unablated[((i+1)<L)*(i+1)][j];

						new_state[(i-1)*(i>0)+(L-1)*(i<1)][j]+=(state[(i-1)*(i>0)+(L-1)*(i<1)][j]<1)*bond[(i-1)*(i>0)+(L-1)*(i<1)][j]*(1-(1-functional[(i-1)*(i>0)+(L-1)*(i<1)][j])*(gsl_rng_uniform(r)<epsilom))*unablated[(i-1)*(i>0)+(L-1)*(i<1)][j];

						new_state[i][(j+1)*(j<L-1)] += (j<L-1)*(state[i][(j+1)*(j<L-1)]<1)*(1-(1-functional[i][(j+1)*(j<L-1)])*(gsl_rng_uniform(r)<epsilom))*unablated[i][(j+1)*(j<L-1)];

						new_state[i][(j-1)*(j>0)]+= (j>0)*(state[i][(j-1)*(j>0)]<1)*(1-(1-functional[i][(j-1)*(j>0)])*(gsl_rng_uniform(r)<epsilom))*unablated[i][(j-1)*(j>0)];
					}
				}
			}
				//end relaxing

			//evolving cells that were not excited previously or now (some subtlty to those possibly excited now)
			for(int i=0;i<L;i++)
			{
				for(int j=0; j<L; j++)
				{
					new_state[i][j] += (state[i][j]>0);
					//if state was excited or refactory before +1 states that were resting may have been newly excited or not but this ignores them
					new_state[i][j] *= (state[i][j]<tau);//promote to resting after end of refactory
					total_excited += (new_state[i][j]<2)*new_state[i][j];
				}
			}

				//SAN
			for(int i=0; i<L; i++)
			{
				new_state[i][0] += (state[i][0]<1)*(heart_time>T-1)*(1-(1-functional[i][0])*(gsl_rng_uniform(r)<epsilom));
					//will do something every T steps should also add in misfire
			}

			heart_time++;
			heart_time *= (heart_time<T+1);
			//end SAN

			state=new_state;

		}
		int get_state(int i, int j)//perhaps have alternative for 2d array
		{

			return state[i][j];

		}
		void set_state(int i, int j, int val)
		{

			state[i][j]=val;
		}
		int get_function(int i, int j)
		{
			return functional[i][j];
		}
		void set_function(int i, int j, int val)
		{
			functional[i][j]=val;
		}
		int get_bond(int i, int j)
		{
			return bond[i][j];
		}
		void set_bond(int i, int j, int val)
		{
			bond[i][j]=val;
		}
		int get_tau()
		{
			return tau;
		}
		int get_T()
		{
			return T;
		}
		void age_nu()
		{

		}
		void set_disfunction()
		{

		}
		int get_excited()
		{
			return total_excited;
		}
		double colour(int i, int j)
		{
			return (1-state[i][j]/(1.0*tau))*(state[i][j]>1);
			//return state[i][j]/(1.0*T);

		}
		void interact(sf::RenderWindow & renderWindow, int length)
		{
			sf::Event event;
			float x;
			float y;

			float Mx = sf::Mouse::getPosition(renderWindow).x;
			float My = sf::Mouse::getPosition(renderWindow).y;
			//std::cout << Mx << " " << My << std::endl;
			while (renderWindow.pollEvent(event))
			{
			switch(event.type)
			{
			case sf::Event::KeyPressed:

				state[My/length][Mx/length] = 0;
				new_state[My/length][Mx/length] = 0;
				unablated[My/length][Mx/length] = 0;

			}
			}
		}
};
#endif

max_stars_repo_path: ablation 1/heart.h

max_issues_repo_path: ablation 1/heart.h

max_forks_repo_path: ablation 1/heart.h


text:
#ifndef _H_INTERP2D
#define _H_INTERP2D

#include <gsl/gsl_spline.h>
#include <gsl/gsl_interp.h>

typedef struct Interpolator2D{
  double *x1;
  double *x2;
  double *y;
  int N1;
  int N2;
  double *y1;
  gsl_spline *x1dir_spline;
  gsl_interp_accel *x1dir_accel;
  gsl_spline **x2dir_splines;
  gsl_interp_accel **x2dir_accels;
} Interpolator2D;

typedef double (*interp2d_modifier_function)(double,double,double,void*);

Interpolator2D * init_interp_2d(double *x1, double *x2, double *y, int N1, int N2, const gsl_interp_type *T);
Interpolator2D * init_interp_2d_akima(double *x1, double *x2, double *y, int N1, int N2);
Interpolator2D * init_interp_2d_grid(double *x1, double *x2, double **y, int N1, int N2, const gsl_interp_type *T);
Interpolator2D * init_interp_2d_akima_grid(double *x1, double *x2, double **y, int N1, int N2);

Interpolator2D * init_interp_2d_function(double *x1, double *x2, double *y,
  int N1, int N2, const gsl_interp_type *T, interp2d_modifier_function function, void * args);
Interpolator2D * init_interp_2d_akima_function(double *x1, double *x2, double *y,
  int N1, int N2, interp2d_modifier_function function, void * args);
Interpolator2D * init_interp_2d_grid_function(double *x1, double *x2, double **y,
  int N1, int N2, const gsl_interp_type *T, interp2d_modifier_function function, void * args);
Interpolator2D * init_interp_2d_akima_grid_function(double *x1, double *x2,
  double **y, int N1, int N2, interp2d_modifier_function function, void * args);

void destroy_interp_2d(Interpolator2D * interp2d);
double interp_2d(double x1, double x2, Interpolator2D * interp2d);
Interpolator2D * load_interp_2d_file(const char * filename, int N1, int N2);

double interp_2d_xmin(Interpolator2D * interp2d);
double interp_2d_xmax(Interpolator2D * interp2d);
double interp_2d_ymin(Interpolator2D * interp2d);
double interp_2d_ymax(Interpolator2D * interp2d);

#endif


text:
/*
 */

#ifndef _aXe_GRISM_H
#define _aXe_GRISM_H

#include <gsl/gsl_matrix.h>
#include <gsl/gsl_vector.h>
#include "spc_trace_functions.h"
//#include "aXe_errors.h"

#define RELEASE "See github"

#define gsl_matrix              gsl_matrix_float
#define gsl_matrix_get          gsl_matrix_float_get
#define gsl_matrix_set          gsl_matrix_float_set
#define gsl_matrix_set_all      gsl_matrix_float_set_all
#define gsl_matrix_alloc        gsl_matrix_float_alloc
#define gsl_matrix_free         gsl_matrix_float_free
#define gsl_matrix_fprintf      gsl_matrix_float_fprintf
#define gsl_matrix_add_constant gsl_matrix_float_add_constant
#define gsl_matrix_scale        gsl_matrix_float_scale
#define gsl_matrix_min          gsl_matrix_float_min
#define gsl_matrix_max          gsl_matrix_float_max

/* WARNING: If you ever change PIXEL_T, make sure you also change
   the FITS loader */
#define PIXEL_T float           /* type of pixel data -- this should be in
                                   sync with the GSL definitions */
#define MAX_BEAMS 27            /* Max. number of spectrum orders we handle */
                                // now equals the alphabet...

#define BEAM(x) (char)((int)x+65)       /* Defines names for 3 different beams */

#define MAXCHAR 255             /* Number of characters in small strings and FITS keywords */

/**
        A set of two relative (pixel) offsets.
*/
typedef struct
{
  int dx0, dx1;
}
px_offset;

/**
        A set of two relative double offsets.
*/
typedef struct
{
  double dx0, dx1;
}
dx_offset;

/**
  A point with integer (pixel) coordinates.
*/
typedef struct
{
  int x, y;
}
px_point;

/**
  A point with double coordinates.
*/
typedef struct
{
  double x, y;
}
d_point;

/**
  An aggregation of the grism exposure, its errors, the background,
  and flatfields
*/
typedef struct
{
  gsl_matrix *grism;            /* The grism image */
  gsl_matrix *pixerrs;  /* (Absolute) errors associated with the grism image */
  gsl_matrix *dq;        /* Data quality associated with grism image */
}
observation;

/**
   A beam, i.e., a description of the image of a spectrum of a given order,
   including a reference point, the corners of a convex quadrangle
   bounding the image, and the parametrization of the spectrum trace
*/
typedef struct
{
  int ID;                       /* An integer ID for this beam
                                   (see enum beamtype */
  d_point  refpoint;            /* Pixel coordinates of the reference point */
  px_point corners[4];          /* Bounding box of spectrum as an arbitrary
                                   quadrangle */
  px_point bbox[2];             /* Bounding box of corners, will be filled
                                   in by spc_extract */
  trace_func *spec_trace;       /* Parametrization of the spectrum trace */
  double width;                 /* largest axis of object in pixels */
  double orient;                /* orientation of largest axis in radian,
                                   0==horizontal */
  int modspec;
  int modimage;

  double awidth;                /* the object width along the extraction direction */
  double bwidth;                /* the object width perpendicular to ewidth
                                   (which is closer to the dispersion direction */
  double aorient;               /* orientation of awidth (CCW) =0 along the x-axis */

  double slitgeom[4];          /* containst the values for the slit length, angle and
                                   so on */

  gsl_vector *flux;              /* the flux of the object at various wavelengths */

  int ignore;           /* This beam should be ignored if this is not set to 0 */
}
beam;

/**
  An object, describing the width and orientation of the object,
  the beams in which the spectrum images lie, and the images
  comprising an observation of the object.
*/
typedef struct
{
     int ID;                    /* An integer ID number for this object */
     beam beams[MAX_BEAMS];     /* table of beams (orders of spectrum) */
     int nbeams;                /* number of beams in beams member */
     observation *grism_obs;
}
object;

/**
  An aperture pixel, i.e., a pixel within a beam.  The structure
  contains the coordinates, the abscissa of the section point between
  the object axis and the spectrum trace, the distance between
  the pixel and the section point along the object axis, the path
  length of the section point along the trace, and the photon count for
  this pixel.  A table of these is produced by the
  make_spc_table function.
*/
typedef struct
{
  int p_x, p_y;         /* the absolute coordinates of the
                                   source pixel. A -1 in p_x signifies
                                   the end of a list of ap_pixels */
  double x, y;          /* the logical coordinates of this (logical)
                                   pixel relative to the beam's reference
                                   point */
  double dist;          /* distance between point and section */
  double xs;                    /* abscissa of section point relative to reference point */
  double ys;                    /* y coord of section point relative to reference point */

  // traditionally (at least aXe-1.3-1.6 "dxs" was called "Projected size of pixel along the trace".
  // however in the code (spc_extract.c, function "handle_one_pixel()" it becomes clear
  // that this quantity indeed is the local trace angle
  double dxs;

  double xi;                    /* path length of the section along the
                                   spectrum trace (computed from xs) */
  double lambda;                /* path length of the section along the
                                   spectrum trace (computed from xi) */
  double dlambda;               /* dpath length of the section along the
                                   spectrum trace */
  double count;         /* intensity for this pixel */
  double weight;     /* extraction weight for this pixel */
  double error;         /* absolute error of count */
  double contam;     /* flag whether this pixel is contaminated by another spectrum */
  double model;     /* the pixel flux computed for one of the emission models */
  long dq;           /* DQ bit information */
}
ap_pixel;

typedef struct
{
  ap_pixel *ap;                 /* A pointer to an ap_pixel array */
  object *obj;                  /* A pointer to the corresponding object
                                   structure */
  beam *b;                      /* A pointer to the beam in the object
                                   that this group of pixel correspond to */
}
PET_entry;

/* See comment to spectrum structure */
typedef enum
{
     SPC_W_CONTAM,
     SPC_W_BORDER
}
spec_warning;

/**
  A point within a spectrum, consisting of the minimum, maximum, and
  (weighted) mean lambda of the points that contributed to the
  bin, as well as an error estimate for the value (that's in there,
  too, of course)
*/
typedef struct
{
  double lambda_min, lambda_mean;
  double lambda_max, dlambda;/* wavelength information */
  double error;         /* error in counts in this bin */
  double count;         /* number of counts in this bin */
  double weight;                /* weight of this bin */
  double flux;          /* flux in physical units in this bin */
  double ferror;                /* error in physical units in this bin */
  double contam;     /* spectral contamination flag */
  //  int    contam;     /* spectral contamination flag */
  long   dq;        /* DQ information */
}
spc_entry;

/**
  A spectrum, consisting of a gsl_array containing the data, a start
  value (center wave length of lowest bin), a bin size, and warn flags,
  which may be

SPC_W_CONTAM -- Spectrum is contaminated

SPC_W_BORDER -- Spectrum lies partially beyond the image borders.

*/
typedef struct
{
  spc_entry *spec;      /* the actual spectrum */
  double lambdamin;
  double lambdamax;
  int spec_len;         /* Number of items in spec */
  spec_warning warning; /* problems with this spectrum */
}
spectrum;

#define DEBUG_FILE 0x01
#define DEBUG 0x0000

extern void
print_ap_pixel_table (const ap_pixel * ap_p);

extern ap_pixel *
make_spc_table (object * const ob,
                const int beamorder, int *const flags);

extern ap_pixel *
make_gps_table (object * const ob,
                const int beamorder, int *const flags,
                int xval, int yval);

#endif /* !_aXe_GRISM_H */

max_stars_repo_name: sosey/pyaxe

max_issues_repo_name: sosey/pyaxe

max_forks_repo_name: sosey/pyaxe


text:
/*
 * ParticleFilter.h
 *
 *  Created on: Mar 21, 2012
 *      Author: nair
 */

#ifndef PARTICLEFILTER_H_
#define PARTICLEFILTER_H_

#include <gsl/gsl_rng.h>
#include <opencv2/opencv.hpp>
#include <Eigen/Dense>

namespace tracking1
{

class ParticleFilterParams
{
public:
    ParticleFilterParams();
    ~ParticleFilterParams();

    void parse(std::string filename);

public:

    int NUMBER_OF_PARTICLES;

    /* standard deviations for gaussian sampling in transition model */
    float TRANS_X_STD;
    float TRANS_Y_STD;
    float TRANS_SCALE_STD;

    /* autoregressive dynamics parameters for transition model */
    float DYNAMICS_A1;
    float DYNAMICS_A2;
    float DYNAMICS_B0;

    double LAMBDA;

    bool DYNAMIC_MODEL_CV;
};

struct Particle
{
    float x;          /**< current x coordinate */
    float y;          /**< current y coordinate */
    float s;          /**< current z coordinate */
    float xp;         /**< previous x coordinate */
    float yp;         /**< previous y coordinate */
    float sp;         /**< previous z coordinate */
    float x0;         /**< original x coordinate */
    float y0;         /**< original y coordinate */
    float s0;			/**< original z coordinate */
    float x_velocity; /**< current x_velocity coordinate */
    float y_velocity; /**< current y_velocity coordinate */
    float s_velocity; /**< current z_velocity coordinate */

    float w;          /**< weight */
};

class ParticleFilter
{
public:

    ParticleFilter();
    ~ParticleFilter();

    void initialize(const cv::Vec3d& initPose);

    void predict();

    void predict_with_simple_model();

    void predict_with_cv_model();

    void correct();

    void resizeParticleSet(int n);
    void setParticleWeight(int id, float w);

    double getMaxLikelihood()
    {
        double max=0.0;
        for(int i=0;i<mParams.NUMBER_OF_PARTICLES;++i)
        {
            if(max<mParticles[i].w)
            {
                max=mParticles[i].w;
            }
        }
        return max;
    }

    double getAvgLikelihood()
    {
        double avg=0.0;
        for(int i=0;i<mParams.NUMBER_OF_PARTICLES;++i)
        {

            avg+=mParticles[i].w;
        }
        return avg/mParams.NUMBER_OF_PARTICLES;
    }

    int nOfParticles()
    {
        return mParams.NUMBER_OF_PARTICLES;
    }

    static int particle_cmp(const void* p1,const void* p2 )
    {
        Particle* _p1 = (Particle*)p1;
        Particle* _p2 = (Particle*)p2;

        if( _p1->w > _p2->w )
            return -1;
        if( _p1->w < _p2->w )
            return 1;
        return 0;
    }

    void getParticlePose(int id, cv::Vec3d& pose);

    void getOutputPose(cv::Vec3d& pose);

    void getOutputAvgPose(cv::Vec3d& pose);

    void getOutputVelocity(cv::Vec3d& vel);

    double getLambda()
    {
        return mParams.LAMBDA;
    }

private:

    void normalize();

    void resample();

private:

    ParticleFilterParams mParams;

    Particle* mParticles;
    cv::Vec3d mOutputPose;

    Eigen::MatrixXd A_Matrix;
    Eigen::MatrixXd G_Matrix;
    Eigen::MatrixXd W_Matrix;

    Eigen::MatrixXd State_Matrix;
    Eigen::MatrixXd State_Prev_Matrix;

    gsl_rng* rng;

};

}
#endif /* PARTICLEFILTER_H_ */

max_stars_repo_name: cognitivesystems/smartcamera

max_issues_repo_name: cognitivesystems/smartcamera

max_forks_repo_name: cognitivesystems/smartcamera


text:
#include "stdio.h"
#include "gsl/gsl_integration.h"
#include "gsl/gsl_errno.h"
#include "non_limber.h"
#include <gsl/gsl_interp2d.h>
#include <gsl/gsl_spline2d.h>
#include <gsl/gsl_spline.h>
#include <gsl/gsl_monte_vegas.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_monte.h>
#include <gsl/gsl_sf_bessel.h>

// This is a workspace size for the gsl integrator
#define LIMBER_FIXED_TABLE_SIZE 4096

// data that is passed into the integrator
// This is everything we need to compute the
// integrand
struct ClIntegrandData {
	double ell;
	double chimin;
	double chimax;
	gsl_spline * WbX;
	gsl_spline * WbY;
    gsl_spline * D_chi;
	gsl_spline2d * P;
	gsl_interp_accel * accelerator_wx;
	gsl_interp_accel * accelerator_wy;
	gsl_interp_accel * accelerator_px;
	gsl_interp_accel * accelerator_py;
	gsl_interp_accel * accelerator_growth;
};

double cl_integrand(double *x, size_t dim, void *data_void)
{
	struct ClIntegrandData * data = (struct ClIntegrandData*) data_void;

	double chi1 = x[0];
	double log_delta = x[1];
	double log_nu = x[2];
	double nu = exp(log_nu);
	double delta = exp(log_delta);
	double chi2_p = chi1 * ( 1 + delta);
	double chi2_m = chi1 * (1 - delta);
	double logk = log_nu / chi1;

	if(chi2_m < data->chimin || chi2_p > data->chimax) return 0.0;

	double FX = gsl_spline_eval(data->WbX, chi1, data->accelerator_wx);
	double FY_p = gsl_spline_eval(data->WbY, chi2_p, data->accelerator_wy);
	double FY_m = gsl_spline_eval(data->WbY, chi2_m, data->accelerator_wy);

	double growth_1 = gsl_spline_eval(data->D_chi, chi1, data->accelerator_growth);
	double growth_2_p = gsl_spline_eval(data->D_chi, chi2_p, data->accelerator_growth);
	double growth_2_m = gsl_spline_eval(data->D_chi, chi2_m, data->accelerator_growth);

	// Get P(k,chi)
	// Deactivate error handling
	gsl_error_handler_t * old_handler = gsl_set_error_handler_off();

	double p;
	int status = gsl_spline2d_eval_e(data->P, chi1, logk, data->accelerator_px, data->accelerator_py, &p);

	// Restore the old error handler
	gsl_set_error_handler(old_handler);
	if (status)
	{
		//printf(stderr,'spline failed, for now, assume this is because k was out of range and return 0...should probably be more careful here though.');
		return 0.0;
	}
	double j1 = gsl_sf_bessel_jl(data->ell, nu);
	double j2_p = gsl_sf_bessel_jl(data->ell, nu*(1+delta));
	double j2_m = gsl_sf_bessel_jl(data->ell, nu*(1-delta));
	double prefac = 2 * delta * 2 * nu * nu * nu * p / growth_1 / M_PI / chi1 / chi1 ;
	return prefac * ( FY_p * growth_2_p * j2_p + FY_m * growth_2_m * j2_m );
}

static double inline gsl_spline_min_x(gsl_spline * s)
{
	return s->x[0];
}
static double inline gsl_spline_max_x(gsl_spline * s)
{
	return s->x[s->size-1];
}

int get_chi_mean_sigma(double chimin, double chimax, gsl_spline * WbX, gsl_spline * WbY,
				gsl_interp_accel * accelerator_wx,	gsl_interp_accel * accelerator_wy,
				double *chi_mean, double *chi_sigma)
{
    int nchi = 1000;
    double dchi = (chimax-chimin)/nchi;
    double chi = chimin;
    double wxwy_sum = 0.;
    double wxwy2_sum = 0.;
    double wx,wy;
    for (int i=0; i<nchi; i++){
    	wx = gsl_spline_eval(WbX, chi, accelerator_wx);
    	wy = gsl_spline_eval(WbY, chi, accelerator_wy);
    	wxwy_sum += wx*wy*chi;
    	wxwy2_sum += wx*wy*chi*chi;
    	chi += dchi;
    }
    *chi_mean = wxwy_sum*dchi;
    *chi_sigma = sqrt(wxwy2_sum*dchi - *chi_mean * *chi_mean);
    return 0;
}

int cl_integral(cl_config * config, gsl_spline * WbX,
			    gsl_spline * WbY, gsl_spline2d * P, gsl_spline * D_chi,
			    double * cl_out, double * cl_err_out)
{

    config->status = LIMBER_STATUS_ERROR;
    int any_parameter_error=0;
    if (WbX==NULL){
        fprintf(stderr, "NULL WX parameter in limber_integral\n");
        any_parameter_error = 1;
    }
    if (WbY==NULL){
        fprintf(stderr, "NULL WY parameter in limber_integral\n");
        any_parameter_error = 1;
    }
    if (P==NULL){
        fprintf(stderr, "NULL P parameter in limber_integral\n");
        any_parameter_error = 1;
    }
    if (D_chi==NULL){
        fprintf(stderr, "NULL D_chi parameter in limber_integral\n");
        any_parameter_error = 1;
    }
    if (config->n_ell<0){
        fprintf(stderr, "Negative n_ell parameter in limber_integral\n");
        any_parameter_error = 1;
    }
    if (any_parameter_error){
        return 1;
    }

	config->status = LIMBER_STATUS_OK;

    static int n_ell_zero_warning = 0;
    if (config->n_ell==0){
        if (n_ell_zero_warning==0){
            fprintf(stderr, "Warning: n_ell=0 in Limber. Will not be treated as an error. Warning once only per process.\n");
        }
        n_ell_zero_warning = 1;
        return 1;
    }

	// Get the appropriate ranges over which to integrate
	// It is assumed that (at least one of) the kernel
	// splines should go to zero in some kind of reasonable
	// place, so we just use the range they specify
	struct ClIntegrandData data;
	double chimin_x = gsl_spline_min_x(WbX);
	double chimin_y = gsl_spline_min_x(WbY);
	double chimax_x = gsl_spline_max_x(WbX);
	double chimax_y = gsl_spline_max_x(WbY);
	double c_ell, c_ell_error;

	// Workspaces
    const gsl_rng_type *T;
    gsl_rng *r;
    gsl_rng_env_setup ();
    T = gsl_rng_default;
    r = gsl_rng_alloc (T);
    size_t calls = 5000000;

	double reltol = config->relative_tolerance;
	double abstol = config->absolute_tolerance;

	// Take the smallest range since we want both the
	// splines to be valid there.
	// This range as well as all the data needed to compute
	// the integrand is put into a struct to be passed
	// through the integrator to the function above.
	double chimin = chimin_x>chimin_y ? chimin_x : chimin_y;
    double chimax = chimax_x<chimax_y ? chimax_x : chimax_y;
	data.chimin = chimin_y;
	data.chimax = chimax_y;
	data.WbX = WbX;
	data.WbY = WbY;
	data.P = P;
	data.accelerator_wx = gsl_interp_accel_alloc();
	data.accelerator_wy = gsl_interp_accel_alloc();
	data.accelerator_px = gsl_interp_accel_alloc();
	data.accelerator_py = gsl_interp_accel_alloc();
	data.accelerator_growth = gsl_interp_accel_alloc();
	data.D_chi = D_chi;

	double chi_mean, chi_sigma;
	get_chi_mean_sigma(chimin, chimax, WbX, WbY,
				data.accelerator_wx, data.accelerator_wy,
				&chi_mean, &chi_sigma);

	// save ell values for return
	double ell_vector[config->n_ell];
	double nu_0;

    gsl_monte_vegas_state *s = gsl_monte_vegas_alloc (3);

	// loop through ell values according to the input configuration
	for (int i_ell = 0; i_ell<config->n_ell; i_ell++){
		double ell = config->ell[i_ell];
		double log_nu_0 = log(ell+0.5);
		double log_nu_min = log_nu_0 - config->log_nu_range;
		double log_nu_max = log_nu_0 + config->log_nu_range;
		printf("delta_range_factor %f\n", config->delta_range_factor);
		double delta_min = 1. / ell / config->delta_range_factor;
		double delta_max = 1. * config->delta_range_factor / ell;

		double log_delta_min = log(delta_min);
		double log_delta_max = log(delta_max);
		printf("log_delta_min,max %f %f\n", delta_min, delta_max);
		double xl[3] = { chimin_x, log_delta_min, log_nu_min };
		double xu[3] = { chimax_x, log_delta_max, log_nu_max };
		data.ell=ell;
		gsl_monte_function G = {&cl_integrand, 3, &data};

        gsl_monte_vegas_integrate (&G, xl, xu, 3, 10000, r, s,
                               &c_ell, &c_ell_error);
        gsl_monte_vegas_integrate (&G, xl, xu, 3, calls/5, r, s,
                               &c_ell, &c_ell_error);
        gsl_monte_vegas_init (s);

		//Include the prefactor scaling
		c_ell *= config->prefactor;
		// Record the results into arrays
		cl_out[i_ell] = c_ell;
		cl_err_out[i_ell] = c_ell_error;
		ell_vector[i_ell] = ell;
	}
	gsl_monte_vegas_free (s);

	// Tidy up
	gsl_interp_accel_free(data.accelerator_wx);
	gsl_interp_accel_free(data.accelerator_wy);
	gsl_interp_accel_free(data.accelerator_px);
	gsl_interp_accel_free(data.accelerator_py);
	gsl_interp_accel_free(data.accelerator_growth);

	// These two are not deallocated because they are static and only initialized once.
	// gsl_integration_glfixed_table_free(table);
	// gsl_integration_workspace_free(W);

	// And that's it
	return 0;
}

ext: c


text:
/**
 *
 * @file testing_ssyev.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
 * @author Azzam Haidar
 * @date 2010-11-15
 * @generated s Tue Jan  7 11:45:18 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_smain.h"

#undef COMPLEX
#define REAL

static int check_orthogonality(int, int, float*, int, float);
static int check_reduction(int, int, int, float*, float*, int, float*, float);
static int check_solution(int, float*, float*, float);

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

    int    N     = atoi(argv[0]);
    int    LDA   = atoi(argv[1]);
    int    LDQ   = LDA;
    int    mode  = 4;
    float eps   = LAPACKE_slamch_work('e');
    float dmax  = 1.0;
    float rcond = 1.0e6;
    int INFO=-1;

    PLASMA_enum uplo = PlasmaUpper;
    PLASMA_enum vec  = PlasmaVec;
    int info_ortho     = 0;
    int info_solution  = 0;
    int info_reduction = 0;
    int LDAxN = LDA*N;
    int LDQxN = LDQ*N;

    float *A1   = NULL;
    float *A2   = (float *)malloc(LDAxN*sizeof(float));
    float *Q    = NULL;
    float             *W1   = (float *)malloc(N*sizeof(float));
    float             *W2   = (float *)malloc(N*sizeof(float));
    float *work = (float *)malloc(3*N* sizeof(float));
    PLASMA_desc *T;

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

    PLASMA_Disable(PLASMA_AUTOTUNING);
    PLASMA_Set(PLASMA_TILE_SIZE, 120);
    PLASMA_Set(PLASMA_INNER_BLOCK_SIZE, 20);

    PLASMA_Enable(PLASMA_WARNINGS);
    PLASMA_Enable(PLASMA_ERRORS);
    PLASMA_Alloc_Workspace_ssyev(N, N, &T);

    /*----------------------------------------------------------
    *  TESTING SSYEV
    */
    /* Initialize A1 */
    LAPACKE_slatms_work( LAPACK_COL_MAJOR, N, N,
                         lapack_const(PlasmaDistSymmetric), ISEED,
                         lapack_const(PlasmaHermGeev), W1, mode, rcond,
                         dmax, N, N,
                         lapack_const(PlasmaNoPacking), A2, LDA, work );

    /*
     * Sort the eigenvalue because when computing the tridiag
     * and then the eigenvalue of the DSTQR are sorted.
     * So to avoid testing fail when having good results W1 should be sorted
    */
    LAPACKE_slasrt_work( 'I', N, W1 );

    if ( vec == PlasmaVec ) {
        Q  = (float *)malloc(LDQxN*sizeof(float));
        A1 = (float *)malloc(LDAxN*sizeof(float));

        /* Copy A2 into A1 */
        LAPACKE_slacpy_work(LAPACK_COL_MAJOR, 'A', N, N, A2, LDA, A1, LDA);
    }
    /*
     * PLASMA SSYEV
     */
    INFO = PLASMA_ssyev(vec, uplo, N, A2, LDA, W2, T, Q, LDQ);
    if(INFO!=0){
            printf(" ERROR OCCURED INFO %d\n",INFO);
            goto fin;
    }

    printf("\n");
    printf("------ TESTS FOR PLASMA SSYEV 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");

    /* Check the orthogonality, reduction and the eigen solutions */
    if (vec == PlasmaVec) {
        info_ortho = check_orthogonality(N, N, Q, LDQ, eps);
        info_reduction = check_reduction(uplo, N, 1, A1, W2, LDA, Q, eps);
    }
    info_solution = check_solution(N, W1, W2, eps);

    if ( (info_solution == 0) & (info_ortho == 0) & (info_reduction == 0) ) {
        printf("***************************************************\n");
        printf(" ---- TESTING SSYEV ...................... PASSED !\n");
        printf("***************************************************\n");
    }
    else {
        printf("************************************************\n");
        printf(" - TESTING SSYEV ... FAILED !\n");
        printf("************************************************\n");
    }

fin:
    PLASMA_Dealloc_Handle_Tile(&T);
    free(A2);
    free(W1);
    free(W2);
    free(work);
    if (Q  != NULL) free(Q);
    if (A1 != NULL) free(A1);

    return 0;
}

/*-------------------------------------------------------------------
 * Check the orthogonality of Q
 */
static int check_orthogonality(int M, int N, float *Q, int LDQ, float eps)
{
    float  done  =  1.0;
    float  mdone = -1.0;
    float  normQ, result;
    int     info_ortho;
    int     minMN = min(M, N);
    float *work = (float *)malloc(minMN*sizeof(float));

    /* Build the idendity matrix */
    float *Id = (float *) malloc(minMN*minMN*sizeof(float));
    LAPACKE_slaset_work(LAPACK_COL_MAJOR, 'A', minMN, minMN, 0., 1., Id, minMN);

    /* Perform Id - Q'Q */
    if (M >= N)
        cblas_ssyrk(CblasColMajor, CblasUpper, CblasTrans, N, M, done, Q, LDQ, mdone, Id, N);
    else
        cblas_ssyrk(CblasColMajor, CblasUpper, CblasNoTrans,   M, N, done, Q, LDQ, mdone, Id, M);

    normQ = LAPACKE_slansy_work(LAPACK_COL_MAJOR, lapack_const(PlasmaInfNorm), 'U', minMN, Id, minMN, work);

    result = normQ / (minMN * eps);
    printf(" ======================================================\n");
    printf(" ||Id-Q'*Q||_oo / (N*eps)          : %15.3E \n",  result );
    printf(" ======================================================\n");

    if ( isnan(result) || isinf(result) || (result > 60.0) ) {
        printf("-- Orthogonality is suspicious ! \n");
        info_ortho=1;
    }
    else {
        printf("-- Orthogonality is CORRECT ! \n");
        info_ortho=0;
    }
    free(work); free(Id);
    return info_ortho;
}

/*------------------------------------------------------------
 *  Check the reduction
 */
static int check_reduction(int uplo, int N, int bw, float *A, float *D, int LDA, float *Q, float eps )
{
    float zone  =  1.0;
    float mzone = -1.0;
    float *TEMP     = (float *)malloc(N*N*sizeof(float));
    float *Residual = (float *)malloc(N*N*sizeof(float));
    float *work = (float *)malloc(N*sizeof(float));
    float Anorm, Rnorm, result;
    int info_reduction;
    int i;

    /* Compute TEMP =  Q * LAMBDA */
    LAPACKE_slacpy_work(LAPACK_COL_MAJOR, lapack_const(PlasmaUpperLower), N, N, Q, LDA, TEMP, N);
    for (i = 0; i < N; i++){
            cblas_sscal(N, D[i], &(TEMP[i*N]),1);
    }
    /* Compute Residual = A - Q * LAMBDA * Q^H */
    /* A is Hermetian but both upper and lower
     * are assumed valable here for checking
     * otherwise it need to be symetrized before
     * checking.
     */
    LAPACKE_slacpy_work(LAPACK_COL_MAJOR, lapack_const(PlasmaUpperLower), N, N, A, LDA, Residual, N);
    cblas_sgemm(CblasColMajor, CblasNoTrans, CblasTrans, N, N, N, (mzone), TEMP, N,  Q, LDA, (zone), Residual,     N);

    Rnorm = LAPACKE_slange_work(LAPACK_COL_MAJOR, lapack_const(PlasmaOneNorm), N, N, Residual, N,   work);
    Anorm = LAPACKE_slansy_work(LAPACK_COL_MAJOR, lapack_const(PlasmaOneNorm), lapack_const(uplo), N, A,       LDA, work);

    result = Rnorm / ( Anorm * N * eps);
    if ( uplo == PlasmaLower ){
        printf(" ======================================================\n");
        printf(" ||A-Q*LAMBDA*Q'||_oo/(||A||_oo.N.eps) : %15.3E \n",  result );
        printf(" ======================================================\n");
    }else{
        printf(" ======================================================\n");
        printf(" ||A-Q'*LAMBDA*Q||_oo/(||A||_oo.N.eps) : %15.3E \n",  result );
        printf(" ======================================================\n");
    }

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

    free(TEMP); free(Residual);
    free(work);

    return info_reduction;
}
/*------------------------------------------------------------
 *  Check the eigenvalues
 */
static int check_solution(int N, float *E1, float *E2, float eps)
{
    int info_solution, i;
    float resid;
    float maxtmp;
    float maxel = fabs( fabs(E1[0]) - fabs(E2[0]) );
    float maxeig = max( fabs(E1[0]), fabs(E2[0]) );
    for (i = 1; i < N; i++){
        resid   = fabs(fabs(E1[i])-fabs(E2[i]));
        maxtmp  = max(fabs(E1[i]), fabs(E2[i]));

        /* Update */
        maxeig = max(maxtmp, maxeig);
        maxel  = max(resid,  maxel );
    }

    maxel = maxel / (maxeig * N * eps);
    printf(" ======================================================\n");
    printf(" | D - eigcomputed | / (|D| * N * eps) : %15.3E \n",  maxel );
    printf(" ======================================================\n");

    if ( isnan(maxel) || isinf(maxel) || (maxel > 100) ) {
        printf("-- The eigenvalues are suspicious ! \n");
        info_solution = 1;
    }
    else{
        printf("-- The eigenvalues are CORRECT ! \n");
        info_solution = 0;
    }
    return info_solution;
}


text:
/**
 * @file bblas_zaccuracy.c
 *
 * @brief BBLAS accuracy tests for double_Complex precision.
 *
 *  BBLAS is a software package provided by Univ. of Manchester,
 *  Univ. of Tennessee.
 *
 * @version 1.0.0
 * @author  Samuel  D. Relton
 * @author  Pedro   V. Lara
 * @author  Mawussi Zounon
 * @date    2016-02-20
 *
 **/

#ifndef DOXYGEN_SHOULD_SKIP_THIS
/**
 * Code generation
 * @precisions normal z -> c d s
 **/
#endif

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

#define COMPLEX

/**
 * bblas_zgemm_tolerance computes the forward error
 * @f[\|C - C_{\mathrm{seq}}\|_{\infty}/((K|\alpha|\|A\|_{\infty}\|B\|_{\infty} +
 *           |\beta|\|C_{\mathrm{init}}\|_{\infty})\epsilon).@f]
 *
 **/

void bblas_zgemm_tolerance(bblas_ztest_t *test, BBLAS_Complex64_t **C_final)
{

    /*Local variables */
    int batch_count = test->batch_count;
    int batch_iter, first_index = 0;
    int Am, An, Bm, Bn, M, K, N, ldc;
    int max_work_size = max(test->maxK,max(test->maxM, test->maxN));
    double Anorm, Bnorm, Rnorm, result;
    double eps = LAPACKE_dlamch_work('e');
    BBLAS_Complex64_t alpha =-1;
    double *work = (double *)malloc(max_work_size*sizeof(double));

    /*Temporary buffer to save (test-arrayC -C_final) */
    BBLAS_Complex64_t **C_diff;
    C_diff = (BBLAS_Complex64_t **) malloc(batch_count*sizeof(BBLAS_Complex64_t *));

    /*Make a copy of test->arrayC in C_diff */
    bblas_zcopy_Cinit(test, C_diff);

    /*Variable size */
    if( test->batch_opts == BBLAS_VARIABLE )
    {
	for( batch_iter =0; batch_iter < batch_count ; batch_iter++)
	{
	    M = test->M[batch_iter];
	    K   = test->K[batch_iter];
	    N   = test->N[batch_iter];
	    ldc = test->ldc[batch_iter];

	    if (test->transA[batch_iter] == BblasNoTrans) {
		Am = M; An = K;
	    } else {
		Am = K; An = M;
	    }
	    if (test->transB[batch_iter] == BblasNoTrans) {
		Bm = K; Bn = N;
	    } else {
		Bm = N; Bn = K;
	    }

	    /*Compute the error C - C_finial */
	    cblas_zaxpy (ldc*N, CBLAS_SADDR(alpha), C_final[batch_iter], 1, C_diff[batch_iter], 1);

	    /*Compute  the infinity norm associated  with the error */
	    Rnorm = LAPACKE_zlange_work(LAPACK_COL_MAJOR, 'I', M, N, C_diff[batch_iter], ldc, work);

	    /*Compute the infinity norm of A */
	    Anorm = LAPACKE_zlange_work(LAPACK_COL_MAJOR, 'I', Am, An, test->arrayA[batch_iter], Am, work);

	    /*Compute the infinity norm of B */
	    Bnorm = LAPACKE_zlange_work(LAPACK_COL_MAJOR, 'I', Bm, Bn, test->arrayB[batch_iter], Bm, work);

	    result = Rnorm / ((K*cabs(test->alpha[batch_iter])*Anorm*Bnorm + cabs(test->beta[batch_iter])*test->Cinitnorm[batch_iter])*eps);

	    /*Compute the relative error */
	    switch(test->target)
	    {
		case BBLAS_MKL:
		    test->mkl_error[batch_iter] = result;
		    break;

		case BBLAS_CUBLAS:
		case BBLAS_MAGMA:
		    test->device_error[batch_iter] = result;
		    break;

		case BBLAS_OTHER:
		    test->other_error[batch_iter] = result;
		    break;
		default:
		    printf("In bblas_zcheck_Cfinal, Variable: Target no defined\n");
		    exit(EXIT_FAILURE);
	    }
	}

    }else  if( test->batch_opts == BBLAS_FIXED ) //fixed  size
    {
	M   = test->M[first_index];
	K   = test->K[first_index];
	ldc = test->ldc[first_index];
	N   = test->N[first_index];

	if (test->transA[first_index] == BblasNoTrans) {
	    Am = M; An = K;
	} else {
	    Am = K; An = M;
	}
	if (test->transB[first_index] == BblasNoTrans) {
	    Bm = K; Bn = N;
	} else {
	    Bm = N; Bn = K;
	}

	for( batch_iter =0; batch_iter < batch_count ; batch_iter++)
	{
		/*Compute the error C - C_finial */
	    cblas_zaxpy (ldc*N, CBLAS_SADDR(alpha), C_final[batch_iter], 1, C_diff[batch_iter], 1);
	    /*Compute  the infinity norm associated  with the error */
	    Rnorm = LAPACKE_zlange_work(LAPACK_COL_MAJOR, 'I', M, N, C_diff[batch_iter], ldc, work);

	    /*Compute the infinity norm of A */
	    Anorm = LAPACKE_zlange_work(LAPACK_COL_MAJOR, 'I', Am, An, test->arrayA[batch_iter], Am, work);

	    /*Compute the infinity norm of B */
	    Bnorm = LAPACKE_zlange_work(LAPACK_COL_MAJOR, 'I', Bm, Bn, test->arrayB[batch_iter], Bm, work);

	    /*Compute the relative error */
	    result = Rnorm / ((K*cabs(test->alpha[first_index])*Anorm*Bnorm + cabs(test->beta[first_index])*test->Cinitnorm[batch_iter])*eps);

	    switch(test->target)
	    {
		case BBLAS_MKL:
		    test->mkl_error[batch_iter] = result;
		    break;

		case BBLAS_CUBLAS:
		case BBLAS_MAGMA:
		    test->device_error[batch_iter] = result;
		    break;

		case BBLAS_OTHER:
		    test->other_error[batch_iter] = result;
		    break;
		default:
		    printf("In bblas_zcheck_Cfinal: Fixed, Target no defined\n");
		    exit(EXIT_FAILURE);
	    }
	}
    }else
    {
	bblas_error("bblas_ztesting.c", "wrong batch_opts value");
    }

/*Free allocated memory */
    for(int index =0; index < batch_count; index++)
    {
	free(C_diff[index]);
    }
    free(C_diff);
    free(work);
}

/**
 * bblas_zstarmm_tolerance computes the forward error
 * @f[\|C - C_{\mathrm{seq}}\|_{\infty}/((K|\alpha|\|A\|_{\infty}\|B\|_{\infty} +
 * |\beta|\|C_{\mathrm{init}}\|_{\infty})\epsilon).@f]
 *
 **/

void bblas_zstarmm_tolerance(bblas_ztest_t *test, BBLAS_Complex64_t **C_final)
{
    /*Local variables */
    int batch_count = test->batch_count;
    int batch_iter, first_index = 0;
    int Am, Bm, Bn, N, ldc;
    int max_work_size = max(test->maxM, test->maxN);
    double Anorm, Bnorm, Rnorm, result;
    double eps = LAPACKE_dlamch_work('e');
    BBLAS_Complex64_t alpha =-1;
    double *work = (double *)malloc(max_work_size*sizeof(double));

    /*Temporary buffer to save (test-arrayC -C_final) */
    BBLAS_Complex64_t **C_diff;
    C_diff = (BBLAS_Complex64_t **) malloc(batch_count*sizeof(BBLAS_Complex64_t *));

    /*Make a copy of test->arrayC in C_diff */
    bblas_zcopy_Cinit(test, C_diff);

    /*Variable size */
    if( test->batch_opts == BBLAS_VARIABLE )
    {
	for( batch_iter =0; batch_iter < batch_count ; batch_iter++)
	{
	    N = test->N[batch_iter];
	    ldc = test->ldc[batch_iter];

	    if(test->side[batch_iter] == BblasLeft )
	    {
		Am = test->M[batch_iter];
	    } else {
		Am = test->N[batch_iter];
	    }
	    Bm = test->ldb[batch_iter];
	    Bn = test->N[batch_iter];
	    /*Compute the error C - C_finial */
	    cblas_zaxpy (ldc*N, CBLAS_SADDR(alpha), C_final[batch_iter], 1, C_diff[batch_iter], 1);

		/*Compute  the infinity norm associated  with the error */
	    Rnorm = LAPACKE_zlange_work(LAPACK_COL_MAJOR, 'I', ldc, N, C_diff[batch_iter], ldc, work);

	    /*Compute the infinity norm of A */
	    Anorm = LAPACKE_zlange_work(LAPACK_COL_MAJOR, 'I', Am, Am, test->arrayA[batch_iter], Am, work);

	    /*Compute the infinity norm of B */
	    Bnorm = LAPACKE_zlange_work(LAPACK_COL_MAJOR, 'I', Bm, Bn, test->arrayB[batch_iter], Bm, work);

	    result = Rnorm / ((N*cabs(test->alpha[batch_iter])*Anorm*Bnorm + cabs(test->beta[batch_iter])*test->Cinitnorm[batch_iter])*eps);

	    /*Compute the relative error */
	    switch(test->target)
	    {
		case BBLAS_MKL:
		    test->mkl_error[batch_iter] = result;
		    break;

		case BBLAS_CUBLAS:
		case BBLAS_MAGMA:
		    test->device_error[batch_iter] = result;
		    break;

		case BBLAS_OTHER:
		    test->other_error[batch_iter] = result;
		    break;
		default:
		    printf("In bblas_zcheck_Cfinal, Variable: Target no defined\n");
		    exit(EXIT_FAILURE);
	    }
	}

    }else  if( test->batch_opts == BBLAS_FIXED ) //fixed  size
    {
	N = test->N[first_index];
	ldc = test->ldc[first_index];

	if(test->side[first_index] == BblasLeft )
	{
	    Am = test->M[first_index];
	} else {
	    Am = test->N[first_index];
	}
	Bm = test->ldb[first_index];
	Bn = test->N[first_index];

	for( batch_iter =0; batch_iter < batch_count ; batch_iter++)
	{
	    /*Compute the error C - C_finial */
	    cblas_zaxpy (ldc*N, CBLAS_SADDR(alpha), C_final[batch_iter], 1, C_diff[batch_iter], 1);

	    /*Compute  the infinity norm associated  with the error */
	    Rnorm = LAPACKE_zlange_work(LAPACK_COL_MAJOR, 'I', ldc, N, C_diff[batch_iter], ldc, work);

	    /*Compute the infinity norm of A */
	    Anorm = LAPACKE_zlange_work(LAPACK_COL_MAJOR, 'I', Am, Am, test->arrayA[batch_iter], Am, work);

	    /*Compute the infinity norm of B */
	    Bnorm = LAPACKE_zlange_work(LAPACK_COL_MAJOR, 'I', Bm, Bn, test->arrayB[batch_iter], Bm, work);

	    /*Compute the relative error */
	    result = Rnorm / ((N*cabs(test->alpha[first_index])*Anorm*Bnorm + cabs(test->beta[first_index])*test->Cinitnorm[batch_iter])*eps);

	    switch(test->target)
	    {
		case BBLAS_MKL:
		    test->mkl_error[batch_iter] = result;
		    break;

		case BBLAS_CUBLAS:
		case BBLAS_MAGMA:
		    test->device_error[batch_iter] = result;
		    break;

		case BBLAS_OTHER:
		    test->other_error[batch_iter] = result;
		    break;
		default:
		    printf("In bblas_zcheck_Cfinal: Fixed, Target no defined\n");
		    exit(EXIT_FAILURE);
	    }
	}
    }else
    {
	bblas_error("bblas_ztesting.c", "wrong batch_opts value");
    }

/*Free allocated memory */
    for(int index =0; index < batch_count; index++)
    {
	free(C_diff[index]);
    }
    free(C_diff);
    free(work);
}

/**
 * bblas_zstark_tolerance computes the forward error
 * @f[\|C - C_{\mathrm{seq}}\|_{\infty}/((K|\alpha|\|A\|_{\infty}\|A\|_1 +
 * |\beta|\|C_{\mathrm{init}}\|_{\infty})\epsilon).@f]
 *
 **/

void bblas_zstark_tolerance(bblas_ztest_t *test, BBLAS_Complex64_t **C_final)
{
    /*Local variables */
    int batch_count = test->batch_count;
    int batch_iter, first_index = 0;
    int Am, An, N, K, ldc;
    char u;
    int max_work_size = max(test->maxK, test->maxN);
    double Anorm, ATnorm, Rnorm, result;
    double eps = LAPACKE_dlamch_work('e');
    double alpha_norm, beta_norm;
    BBLAS_Complex64_t alpha =-1;
    double *work = (double *)malloc(max_work_size*sizeof(double));

    /*Temporary buffer to save (test-arrayC -C_final) */
    BBLAS_Complex64_t **C_diff;
    C_diff = (BBLAS_Complex64_t **) malloc(batch_count*sizeof(BBLAS_Complex64_t *));

    /*Make a copy of test->arrayC in C_diff */
    bblas_zcopy_Cinit(test, C_diff);

    /*Variable size */
    if( test->batch_opts == BBLAS_VARIABLE )
    {
	for( batch_iter =0; batch_iter < batch_count ; batch_iter++)
	{
	    N = test->N[batch_iter];
	    K = test->K[batch_iter];
	    ldc = test->ldc[batch_iter];

	    if (test->trans[batch_iter] == BblasNoTrans) {
		Am = N; An = K;
	    } else {
		Am = K; An = N;
	    }

	    if(test->uplo[batch_iter] == BblasLower )
	      {
		u = 'L';
	      }else
	      {
		u = 'U';
	      }
	    /*Compute the error C - C_finial */
	    cblas_zaxpy (ldc*N, CBLAS_SADDR(alpha), C_final[batch_iter], 1, C_diff[batch_iter], 1);

	    /*Compute  the infinity norm associated  with the error */
	    Rnorm = LAPACKE_zlanhe_work(LAPACK_COL_MAJOR, 'I', u, N, C_diff[batch_iter], ldc, work);

	    /*Compute the infinity norm of A */
	    Anorm = LAPACKE_zlange_work(LAPACK_COL_MAJOR, 'I', Am, An, test->arrayA[batch_iter], Am, work);

	    /*Compute the one norm of A */
	    ATnorm = LAPACKE_zlange_work(LAPACK_COL_MAJOR, 'O', Am, An, test->arrayA[batch_iter], Am, work);

	    if (test->routine == BBLAS_HERK)
	    {
		beta_norm   = cabs(test->beta_herk[batch_iter]);
		alpha_norm  = cabs(test->alpha_herk[batch_iter]);
	    }else{
		beta_norm  = cabs(test->beta[batch_iter]);
		alpha_norm  = cabs(test->alpha[batch_iter]);
	    }

	    result = Rnorm / ((K*alpha_norm*Anorm*ATnorm + beta_norm*test->Cinitnorm[batch_iter])*eps);

	    /*Compute the relative error */
	    switch(test->target)
	    {
		case BBLAS_MKL:
		    test->mkl_error[batch_iter] = result;
		    break;

		case BBLAS_CUBLAS:
		case BBLAS_MAGMA:
		    test->device_error[batch_iter] = result;
		    break;

		case BBLAS_OTHER:
		    test->other_error[batch_iter] = result;
		    break;
		default:
		    printf("In bblas_zcheck_Cfinal, Variable: Target no defined\n");
		    exit(EXIT_FAILURE);
	    }
	}

    }else  if( test->batch_opts == BBLAS_FIXED ) //fixed  size
    {
	    N = test->N[first_index];
	    K = test->K[first_index];
	    ldc = test->ldc[first_index];

	    if (test->trans[first_index] == BblasNoTrans) {
		Am = N; An = K;
	    } else {
		Am = K; An = N;
	    }
	    if(test->uplo[first_index] == BblasLower )
	      {
		u = 'L';
	      }else
	      {
		u = 'U';
	      }

	    for( batch_iter =0; batch_iter < batch_count ; batch_iter++)
	      {
		/*Compute the error C - C_finial */
		cblas_zaxpy (ldc*N, CBLAS_SADDR(alpha), C_final[batch_iter], 1, C_diff[batch_iter], 1);

		Rnorm = LAPACKE_zlanhe_work(LAPACK_COL_MAJOR, 'I', u, N, C_diff[batch_iter], ldc, work);

		/*Compute the infinity norm of A */
		Anorm = LAPACKE_zlange_work(LAPACK_COL_MAJOR, 'I', Am, An, test->arrayA[batch_iter], Am, work);

		/*Compute the one norm of A */
		ATnorm = LAPACKE_zlange_work(LAPACK_COL_MAJOR, 'O', Am, An, test->arrayA[batch_iter], Am, work);

		if (test->routine == BBLAS_HERK)
		{
		    beta_norm   = cabs(test->beta_herk[first_index]);
		    alpha_norm  = cabs(test->alpha_herk[first_index]);
		}else{
		    beta_norm   = cabs(test->beta[first_index]);
		    alpha_norm  = cabs(test->alpha[first_index]);
		}

		/*Compute the relative error */
		result = Rnorm / ((K*alpha_norm*Anorm*ATnorm +
				   beta_norm*test->Cinitnorm[batch_iter])*eps);

		switch(test->target)
		{
		    case BBLAS_MKL:
			test->mkl_error[batch_iter] = result;
			break;

		case BBLAS_CUBLAS:
		case BBLAS_MAGMA:
		  test->device_error[batch_iter] = result;
		  break;

		    case BBLAS_OTHER:
			test->other_error[batch_iter] = result;
			break;
		    default:
			printf("In bblas_zcheck_Cfinal: Fixed, Target no defined\n");
			exit(EXIT_FAILURE);
		}
	    }
    }else
    {
	bblas_error("bblas_ztesting.c", "wrong batch_opts value");
    }

/*Free allocated memory */
    for(int index =0; index < batch_count; index++)
    {
	free(C_diff[index]);
    }
    free(C_diff);
    free(work);
}

/**
 * bblas_zstar2k_tolerance computes the forward error
 * @f[\|C - C_{\mathrm{seq}}\|_{\infty}/((K|\alpha|\|A\|_{\infty}\|B\|_1 +
 * K|\alpha|\|B\|_{\infty}\|A\|_1  +
 * |\beta|\|C_{\mathrm{init}}\|_{\infty})\epsilon).@f]
 *
 **/

void bblas_zstar2k_tolerance(bblas_ztest_t *test, BBLAS_Complex64_t **C_final)
{
    /*Local variables */
    int batch_count = test->batch_count;
    int batch_iter, first_index = 0;
    int Am, An, Bn, Bm, N, K, ldc;
    int max_work_size = max(test->maxK, test->maxN);
    double Anorm, Bnorm, ATnorm, BTnorm,  Rnorm, result;
    double beta_norm;
    double eps = LAPACKE_dlamch_work('e');
    BBLAS_Complex64_t alpha =-1;
    double *work = (double *)malloc(max_work_size*sizeof(double));

    /*Temporary buffer to save (test-arrayC -C_final) */
    BBLAS_Complex64_t **C_diff;
    C_diff = (BBLAS_Complex64_t **) malloc(batch_count*sizeof(BBLAS_Complex64_t *));

    /*Make a copy of test->arrayC in C_diff */
    bblas_zcopy_Cinit(test, C_diff);

    /*Variable size */
    if( test->batch_opts == BBLAS_VARIABLE )
    {
	for( batch_iter =0; batch_iter < batch_count ; batch_iter++)
	{
	    N = test->N[batch_iter];
	    K = test->K[batch_iter];
	    ldc = test->ldc[batch_iter];

	    if (test->trans[batch_iter] == BblasNoTrans) {
		Am = N; An = K;
		Bn = K;
	    } else {
		Am = K; An = N;
		Bn = N;
	    }
	    Bm = test->ldb[batch_iter];

	    /*Compute the error C - C_finial */
	    cblas_zaxpy (ldc*N, CBLAS_SADDR(alpha), C_final[batch_iter], 1, C_diff[batch_iter], 1);

		/*Compute  the infinity norm associated  with the error */
	    Rnorm = LAPACKE_zlange_work(LAPACK_COL_MAJOR, 'I', ldc, N, C_diff[batch_iter], ldc, work);

	    /*Compute the infinity norm of A */
	    Anorm = LAPACKE_zlange_work(LAPACK_COL_MAJOR, 'I', Am, An, test->arrayA[batch_iter], Am, work);

	    /*Compute the one norm of A */
	    ATnorm = LAPACKE_zlange_work(LAPACK_COL_MAJOR, 'O', Am, An, test->arrayA[batch_iter], Am, work);

	    /*Compute the infinity norm of B */
	    Bnorm = LAPACKE_zlange_work(LAPACK_COL_MAJOR, 'I', Bm, Bn, test->arrayB[batch_iter], Bm, work);

	    /*Compute the one norm of B */
	    BTnorm = LAPACKE_zlange_work(LAPACK_COL_MAJOR, 'O', Bm, Bn, test->arrayB[batch_iter], Bm, work);

	    if(test->routine == BBLAS_HER2K)
	    {
		beta_norm = cabs(test->beta_herk[batch_iter]);
	    }else
	    {
		beta_norm = cabs(test->beta[batch_iter]);
	    }

	    result = Rnorm / ((K*cabs(test->alpha[batch_iter])*Anorm*BTnorm +
			       K*cabs(test->alpha[batch_iter])*Bnorm*ATnorm +
			       beta_norm*test->Cinitnorm[batch_iter])*eps);

	    /*Compute the relative error */
	    switch(test->target)
	    {
		case BBLAS_MKL:
		    test->mkl_error[batch_iter] = result;
		    break;

		case BBLAS_CUBLAS:
		case BBLAS_MAGMA:
		    test->device_error[batch_iter] = result;
		    break;

		case BBLAS_OTHER:
		    test->other_error[batch_iter] = result;
		    break;
		default:
		    printf("In bblas_zcheck_Cfinal, Variable: Target no defined\n");
		    exit(EXIT_FAILURE);
	    }
	}

    }else  if( test->batch_opts == BBLAS_FIXED ) //fixed  size
    {
	N = test->N[first_index];
	K = test->K[first_index];
	ldc = test->ldc[first_index];

	if (test->trans[first_index] == BblasNoTrans) {
		Am = N; An = K;
		Bn = K;
	} else {
	    Am = K; An = N;
	    Bn = N;
	}
	Bm = test->ldb[first_index];

	    for( batch_iter =0; batch_iter < batch_count ; batch_iter++)
	    {
		/*Compute the error C - C_finial */
		cblas_zaxpy (ldc*N, CBLAS_SADDR(alpha), C_final[batch_iter], 1, C_diff[batch_iter], 1);

		/*Compute  the infinity norm associated  with the error */
		Rnorm = LAPACKE_zlange_work(LAPACK_COL_MAJOR, 'I', ldc, N, C_diff[batch_iter], ldc, work);

		/*Compute the infinity norm of A */
		Anorm = LAPACKE_zlange_work(LAPACK_COL_MAJOR, 'I', Am, An, test->arrayA[batch_iter], Am, work);

		/*Compute the one norm of A */
		ATnorm = LAPACKE_zlange_work(LAPACK_COL_MAJOR, 'O', Am, An, test->arrayA[batch_iter], Am, work);

		/*Compute the infinity norm of B */
		Bnorm = LAPACKE_zlange_work(LAPACK_COL_MAJOR, 'I', Bm, Bn, test->arrayB[batch_iter], Bm, work);

		/*Compute the one norm of B */
		BTnorm = LAPACKE_zlange_work(LAPACK_COL_MAJOR, 'O', Bm, Bn, test->arrayB[batch_iter], Bm, work);

		if(test->routine == BBLAS_HER2K)
		{
		    beta_norm = cabs(test->beta_herk[first_index]);
		}else
		{
		    beta_norm = cabs(test->beta[first_index]);
		}

		/*Compute the relative error */
		result = Rnorm / ((K*cabs(test->alpha[first_index])*Anorm*BTnorm +
				   K*cabs(test->alpha[first_index])*Bnorm*ATnorm +
				   beta_norm*test->Cinitnorm[batch_iter])*eps);

		switch(test->target)
		{
		case BBLAS_MKL:
		  test->mkl_error[batch_iter] = result;
		  break;

		case BBLAS_CUBLAS:
		case BBLAS_MAGMA:
		  test->device_error[batch_iter] = result;
		  break;

		case BBLAS_OTHER:
		  test->other_error[batch_iter] = result;
		  break;
		default:
		  printf("In bblas_zcheck_Cfinal: Fixed, Target no defined\n");
		  exit(EXIT_FAILURE);
		}
	    }
    }else
    {
	bblas_error("bblas_ztesting.c", "wrong batch_opts value");
    }

/*Free allocated memory */
    for(int index =0; index < batch_count; index++)
    {
	free(C_diff[index]);
    }
    free(C_diff);
    free(work);
}

/**
 * bblas_ztrmm_tolerance computes the forward error
 * @f[\|C - C_{\mathrm{seq}}\|_{\infty}/((|\alpha|\|A\|_{\infty}\|B\|_{\infty}N\epsilon).@f]
 *
 **/

void bblas_ztrmm_tolerance(bblas_ztest_t *test, BBLAS_Complex64_t **B_final)
{
    /*Local variables */
    int batch_count = test->batch_count;
    int batch_iter, first_index = 0;
    int Am, An, ldb, N, M;
    int max_work_size = max(test->maxK, test->maxN);
    char u, diag;
    double Anorm, Rnorm, result;
    double eps = LAPACKE_dlamch_work('e');
    BBLAS_Complex64_t alpha =-1;
    double *work = (double *)malloc(max_work_size*sizeof(double));

    /*Temporary buffer to save (test-arrayB -B_final) */
    BBLAS_Complex64_t **B_diff;
    B_diff = (BBLAS_Complex64_t **) malloc(batch_count*sizeof(BBLAS_Complex64_t *));

    /*Make a copy of test->arrayB in B_diff */
    bblas_zcopy_Binit(test, B_diff);

    /*Variable size */
    if( test->batch_opts == BBLAS_VARIABLE )
    {
	for( batch_iter =0; batch_iter < batch_count ; batch_iter++)
	{
	    N = test->N[batch_iter];
	    M = test->M[batch_iter];
	    ldb = test->ldb[batch_iter];

	    Am = test->lda[batch_iter];

	    if(test->side[batch_iter] == BblasLeft )
	    {
		An = M;
	    }else
	    {
		An = N;
	    }

	    if(test->uplo[batch_iter] == BblasLower )
	    {
	       u = 'L';
	    }else
	    {
	       u = 'U';
	    }

	    if(test->diag[batch_iter] == BblasUnit )
	    {
		diag = 'U';
	    }else
	    {
		diag = 'N';
	    }

	    /*Compute the error B - B_finial */
	    cblas_zaxpy (ldb*N, CBLAS_SADDR(alpha), B_final[batch_iter], 1, B_diff[batch_iter], 1);

		/*Compute  the infinity norm associated  with the error */
	    Rnorm = LAPACKE_zlange_work(LAPACK_COL_MAJOR, 'I', ldb, N, B_diff[batch_iter], ldb, work);

	    /*Infinity norm of the triangular matrix A */
	    Anorm = LAPACKE_zlantr_work(LAPACK_COL_MAJOR, 'I', u, diag, Am,
					An, test->arrayA[batch_iter], Am, work);

		result = Rnorm/(N*cabs(test->alpha[batch_iter])*Anorm*test->Binitnorm[batch_iter]*eps);
	    /*Compute the relative error */
	    switch(test->target)
	    {
		case BBLAS_MKL:
		    test->mkl_error[batch_iter] = result;
		    break;

		case BBLAS_CUBLAS:
		case BBLAS_MAGMA:
		    test->device_error[batch_iter] = result;
		    break;

		case BBLAS_OTHER:
		    test->other_error[batch_iter] = result;
		    break;
		default:
		    printf("In bblas_zcheck_Cfinal, Variable: Target no defined\n");
		    exit(EXIT_FAILURE);
	    }

	}
    }else  if( test->batch_opts == BBLAS_FIXED )
    {

	N = test->N[first_index];
	M = test->M[first_index];
	ldb = test->ldb[first_index];

	Am = test->lda[first_index];

	if(test->side[first_index] == BblasLeft )
	{
	    An = M;
	}else
	{
	    An = N;
	}

	if(test->uplo[first_index] == BblasLower )
	{
	    u = 'L';
	}else
	{
	    u = 'U';
	}

	if(test->diag[first_index] == BblasUnit )
	{
	    diag = 'U';
	}else
	{
	    diag = 'N';
	}

	for( batch_iter =0; batch_iter < batch_count ; batch_iter++)
	{
	    /*Compute the error B - B_finial */
	    cblas_zaxpy (ldb*N, CBLAS_SADDR(alpha), B_final[batch_iter], 1, B_diff[batch_iter], 1);

		/*Compute  the infinity norm associated  with the error */
	    Rnorm = LAPACKE_zlange_work(LAPACK_COL_MAJOR, 'I', ldb, N, B_diff[batch_iter], ldb, work);

	    /*Infinity norm of the triangular matrix A */
	    Anorm = LAPACKE_zlantr_work(LAPACK_COL_MAJOR, 'I', u, diag, Am,
					An, test->arrayA[batch_iter], Am, work);

	    result = Rnorm/(N*cabs(test->alpha[first_index])*Anorm*test->Binitnorm[batch_iter]*eps);

	    /*Compute the relative error */
	    switch(test->target)
	    {
		case BBLAS_MKL:
		    test->mkl_error[batch_iter] = result;
		    break;

		case BBLAS_CUBLAS:
		case BBLAS_MAGMA:
		    test->device_error[batch_iter] = result;
		    break;

		case BBLAS_OTHER:
		    test->other_error[batch_iter] = result;
		    break;
		default:
		    printf("In bblas_zcheck_Cfinal, Variable: Target no defined\n");
		    exit(EXIT_FAILURE);
	    }
	}
    }
    else
    {
	bblas_error("bblas_ztesting.c", "wrong batch_opts value");
    }

    /*Free allocated memory */
    for(int index =0; index < batch_count; index++)
    {
	free(B_diff[index]);
    }
    free(B_diff);
    free(work);
}

/**
 * bblas_ztrsm_tolerance computes the backward error
 * @f[\|\alpha B - AX \|_{\infty}/((\|A\|_{\infty}\|X\| + |\alpha|\|B\|_{\infty})N\epsilon).@f]
 *
 **/

void bblas_ztrsm_tolerance(bblas_ztest_t *test, BBLAS_Complex64_t **B_final)
{
    /*Local variables */
    int batch_count = test->batch_count;
    int batch_iter, first_index = 0;
    int Am, An, ldb, N, M;
    int max_work_size = max(test->maxK, test->maxN);
    char u, diag;
    double Anorm, Rnorm, result, Xnorm;
    double eps = LAPACKE_dlamch_work('e');
    BBLAS_Complex64_t alpha;
    double *work = (double *)malloc(max_work_size*sizeof(double));

    /*Variable to save the residual alphaB -AX */
    BBLAS_Complex64_t **residual;
    residual = (BBLAS_Complex64_t **) malloc(batch_count*sizeof(BBLAS_Complex64_t *));

    /*Make a copy of test->arrayB (X) in residual */
    bblas_zcopy_Binit(test, residual);

    /*Variable size */
    if( test->batch_opts == BBLAS_VARIABLE )
    {
	for( batch_iter =0; batch_iter < batch_count ; batch_iter++)
	{
	    N = test->N[batch_iter];
	    M = test->M[batch_iter];
	    ldb = test->ldb[batch_iter];

	    Am = test->lda[batch_iter];

	    if(test->side[batch_iter] == BblasLeft )
	    {
		An = M;
	    }else
	    {
		An = N;
	    }

	    if(test->uplo[batch_iter] == BblasLower )
	    {
	       u = 'L';
	    }else
	    {
	       u = 'U';
	    }

	    if(test->diag[batch_iter] == BblasUnit )
	    {
		diag = 'U';
	    }else
	    {
		diag = 'N';
	    }

	    /*Compute residual <- AX  (ztrmm) */
	    alpha = 1;
	    cblas_ztrmm(BblasColMajor, test->side[batch_iter], test->uplo[batch_iter], test->transA[batch_iter],
			test->diag[batch_iter], M, N, CBLAS_SADDR(alpha), test->arrayA[batch_iter],
			test->lda[batch_iter], residual[batch_iter], ldb);

	    /*Compute Binit = alpha Binit */
	    for (int index=0;  index < ldb*N; index++)
	    {
		test->Binit[batch_iter][index] = test->alpha[batch_iter]*test->Binit[batch_iter][index];
	    }

	    /*Compute the  residual <- alpha B - Ax  */
	    alpha= -1;
	    cblas_zaxpy (ldb*N, CBLAS_SADDR(alpha), test->Binit[batch_iter], 1, residual[batch_iter], 1);
		/*Compute  the infinity norm associated  with the residual */
	    Rnorm = LAPACKE_zlange_work(LAPACK_COL_MAJOR, 'I', ldb, N, residual[batch_iter], ldb, work);

	    /*Infinity norm of the triangular matrix A */
	    Anorm = LAPACKE_zlantr_work(LAPACK_COL_MAJOR, 'I', u, diag, Am,
					An, test->arrayA[batch_iter], Am, work);

	    /*Compute  the infinity norm associated  with X (B) */
	    Xnorm = LAPACKE_zlange_work(LAPACK_COL_MAJOR, 'I', ldb, N, test->arrayB[batch_iter], ldb, work);

	    result = Rnorm/(Anorm*Xnorm*N*eps);

	    /*Compute the relative error */
	    switch(test->target)
	    {
		case BBLAS_MKL:
		    test->mkl_error[batch_iter] = result;
		    break;

		case BBLAS_CUBLAS:
		case BBLAS_MAGMA:
		    test->device_error[batch_iter] = result;
		    break;

		case BBLAS_OTHER:
		    test->other_error[batch_iter] = result;
		    break;
		default:
		    printf("In bblas_zcheck_Cfinal, Variable: Target no defined\n");
		    exit(EXIT_FAILURE);
	    }

	}
    }else  if( test->batch_opts == BBLAS_FIXED )
    {

	N = test->N[first_index];
	M = test->M[first_index];
	ldb = test->ldb[first_index];

	Am = test->lda[first_index];

	if(test->side[first_index] == BblasLeft )
	{
	    An = M;
	}else
	{
	    An = N;
	}

	if(test->uplo[first_index] == BblasLower )
	{
	    u = 'L';
	}else
	{
	    u = 'U';
	}

	if(test->diag[first_index] == BblasUnit )
	{
	    diag = 'U';
	}else
	{
	    diag = 'N';
	}

	for( batch_iter =0; batch_iter < batch_count ; batch_iter++)
	{

	    /*Compute residual <- AX  (ztrmm) */
	    alpha = 1;
	    cblas_ztrmm(BblasColMajor, test->side[first_index], test->uplo[first_index], test->transA[first_index],
			test->diag[first_index], M, N, CBLAS_SADDR(alpha), (const BBLAS_Complex64_t*) test->arrayA[batch_iter],
			test->lda[first_index], residual[batch_iter], ldb);

	    /*Compute Binit = alpha Binit */
	    for (int index=0;  index < ldb*N; index++)
	    {
		test->Binit[batch_iter][index] = test->alpha[first_index]*test->Binit[batch_iter][index];
	    }

	    /*Compute the  residual <- alpha B - Ax  */
	    alpha = -1;
	    cblas_zaxpy (ldb*N, CBLAS_SADDR(alpha), test->Binit[batch_iter], 1, residual[batch_iter], 1);

		/*Compute  the infinity norm associated  with the residual */
	    Rnorm = LAPACKE_zlange_work(LAPACK_COL_MAJOR, 'I', ldb, N, residual[batch_iter], ldb, work);

	    /*Infinity norm of the triangular matrix A */
	    Anorm = LAPACKE_zlantr_work(LAPACK_COL_MAJOR, 'I', u, diag, Am,
					An, test->arrayA[batch_iter], Am, work);

	    /*Compute  the infinity norm associated  with X (B) */
	    Xnorm = LAPACKE_zlange_work(LAPACK_COL_MAJOR, 'I', ldb, N, test->arrayB[batch_iter], ldb, work);

	    result = Rnorm/(Anorm*Xnorm*N*eps);

	    /*Compute the relative error */
	    switch(test->target)
	    {
		case BBLAS_MKL:
		    test->mkl_error[batch_iter] = result;
		    break;

		case BBLAS_CUBLAS:
		case BBLAS_MAGMA:
		    test->device_error[batch_iter] = result;
		    break;

		case BBLAS_OTHER:
		    test->other_error[batch_iter] = result;
		    break;
		default:
		    printf("In bblas_zcheck_Cfinal, Variable: Target no defined\n");
		    exit(EXIT_FAILURE);
	    }
	}
    }
    else
    {
	bblas_error("bblas_ztesting.c", "wrong batch_opts value");
    }

    /*Free allocated memory */
    for(int index =0; index < batch_count; index++)
    {
	free(residual[index]);
    }
    free(residual);
    free(work);
}

/**
 * Calls bblas_zstar2k_tolerance.
 **/
#ifdef COMPLEX
void bblas_zher2k_tolerance(bblas_ztest_t *test, BBLAS_Complex64_t **C_final)
{
    bblas_zstar2k_tolerance(test, C_final);

}
#endif
/**
 * Calls bblas_zstar2k_tolerance.
 **/

void bblas_zsyr2k_tolerance(bblas_ztest_t *test, BBLAS_Complex64_t **C_final)
{
    bblas_zstar2k_tolerance(test, C_final);
}

/**
 * Calls bblas_zstark_tolerance.
 **/

#ifdef COMPLEX
void bblas_zherk_tolerance(bblas_ztest_t *test, BBLAS_Complex64_t **C_final)
{
    bblas_zstark_tolerance(test, C_final);
}
#endif
/**
 * Calls bblas_zstark_tolerance.
 **/

void bblas_zsyrk_tolerance(bblas_ztest_t *test, BBLAS_Complex64_t **C_final)
{
    bblas_zstark_tolerance(test, C_final);
}

/**
 * Calls bblas_zstarmm_tolerance.
 **/

#ifdef COMPLEX
void bblas_zhemm_tolerance(bblas_ztest_t *test, BBLAS_Complex64_t **C_final)
{
    bblas_zstarmm_tolerance(test, C_final);
}
#endif
/**
 * Calls bblas_zstamm_tolerance.
 **/

void bblas_zsymm_tolerance(bblas_ztest_t *test, BBLAS_Complex64_t **C_final)
{
    bblas_zstarmm_tolerance(test, C_final);
}

max_stars_repo_name: NLAFET/BBLAS-ref

max_issues_repo_name: NLAFET/BBLAS-ref

max_forks_repo_name: NLAFET/BBLAS-ref


text:
#include <stdlib.h>
#include <cblas.h>
#include <lapacke.h>

#include "../../runtime/pcp.h"
#include "tasks.h"
#include "utils.h"

#define BLOCK_SIZE 32
#define MAX_MATRIX_SIZE 1024
#define MAX_NUM_ITER (MAX_MATRIX_SIZE / BLOCK_SIZE)
#define MAX_NUM_THREADS 16

// The matrix partitioning for TRSM and SYRK in each iteration.
//
// A partitioning x of a range 0, ..., n - 1 over p threads is
// represented as an array x[] of length p + 1 such that x[k] is the
// start of the block assigned to thread k and x[k + 1] - x[k] is the
// size of the block assigned to thread k, for k = 0, ..., p - 1.
struct partition
{
    // The matrix size.
    int matrix_size;

    // The block size.
    int block_size;

    // The number of threads.
    int num_threads;

    // The number of iterations = ceil(matrix_size, block_size).
    int num_iter;

    // The partitionings for TRSM.
    int trsm[MAX_NUM_ITER][MAX_NUM_THREADS + 1];

    // The partitionings for SYRK.
    int syrk[MAX_NUM_ITER][MAX_NUM_THREADS + 1];
};

// Barrier used to synchronize the workers between iterations.
static pcp_barrier_t *barrier = NULL;

static struct partition *partition_db[MAX_MATRIX_SIZE + 1][MAX_NUM_THREADS + 1] = {NULL};
static struct partition *partition = NULL;

static double *A;
static int ldA;
static int n;
#define A(i,j) A[(i) + (j) * ldA]

// Initializes the partitioning for TRSM for a given iteration.
//
// Uses uniform row blocks.
static void part_init_trsm(int iter)
{
    const int last = partition->matrix_size;
    const int first = min(last, (iter + 1) * partition->block_size);
    const int size = last - first;
    const int p = partition->num_threads;
    const int chunk = iceil(size, p);
    partition->trsm[iter][0] = first;
    partition->trsm[iter][p] = last;
    for (int th = 1; th < p; ++th) {
        partition->trsm[iter][th] = min(last, partition->trsm[iter][th - 1] + chunk);
    }
}

// Initializes the partitioning for SYRK for a given iteration.
static void part_init_syrk(int iter)
{
    // Balance the load by flops.
    const int last = partition->matrix_size;
    const int first = min(last, (iter + 1) * partition->block_size);
    const int size = last - first;
    const int p = partition->num_threads;
    const int total_work = size * (size + 1) / 2;
    const int ideal_part_work = total_work / p;
    partition->syrk[iter][0] = first;
    partition->syrk[iter][p] = last;
    for (int k = 1; k < p; ++k) {
        partition->syrk[iter][k] = partition->syrk[iter][k - 1];
        int work = 0;
        while (work < ideal_part_work && partition->syrk[iter][k] < last) {
            partition->syrk[iter][k] += 1;
            work += partition->syrk[iter][k] - first;
        }
        if (k == 1) {
            partition->syrk[iter][k] = max(partition->syrk[iter][k], partition->syrk[iter][k - 1] + partition->block_size);
            partition->syrk[iter][k] = min(partition->syrk[iter][k], last);
        }
    }
}

// Returns the start of a TRSM block.
static int part_block_start_trsm(int iter, int me)
{
    return partition->trsm[iter][me];
}

// Returns the size of a TRSM block.
static int part_block_size_trsm(int iter, int me)
{
    return part_block_start_trsm(iter, me + 1) - part_block_start_trsm(iter, me);
}

// Returns the start of a SYRK block.
static int part_block_start_syrk(int iter, int me)
{
    return partition->syrk[iter][me];
}

// Returns the size of a SYRK block.
static int part_block_size_syrk(int iter, int me)
{
    return part_block_start_syrk(iter, me + 1) - part_block_start_syrk(iter, me);
}

// Initializes the partitioning for given parameters.
//
// Follows this protocol:
// 1) Re-use existing if it matches.
// 2) Load from file if it matches.
// 3) Create from scratch.
static void part_init(int matrix_size, int block_size, int num_threads)
{
    // Check if the partitioning is available in the databse.
    if (partition_db[matrix_size][num_threads] == NULL) {
        // No, create from scratch and save in the database.
        partition = malloc(sizeof(*partition));
        partition->matrix_size = matrix_size;
        partition->block_size = block_size;
        partition->num_threads = num_threads;
        partition->num_iter = iceil(matrix_size, block_size);
        for (int iter = 0; iter < partition->num_iter; ++iter) {
            part_init_trsm(iter);
            part_init_syrk(iter);
        }
        partition_db[matrix_size][num_threads] = partition;
    } else {
        // Yes, reuse from database.
        partition = partition_db[matrix_size][num_threads];
    }
}

void chol_task_par_reconfigure(int nth)
{
    // Choose current size of the worker pool for barrier.
    if (barrier != NULL) {
        pcp_barrier_destroy(barrier);
    }
    barrier = pcp_barrier_create(nth);
}

void chol_task_par_finalize(void)
{
    if (barrier != NULL) {
        pcp_barrier_destroy(barrier);
        barrier = NULL;
    }
}

static void krnl_chol(int iter)
{
    int jp      = iter * partition->block_size;
    int n       = min(partition->matrix_size - jp, partition->block_size);
    double *A11 = &A(jp,jp);

    LAPACKE_dpotrf(LAPACK_COL_MAJOR, 'L', n, A11, ldA);
}

static void krnl_syrk(int iter, int me)
{
    //   j1        j2
    // i g g g g g s
    //   g g g g g s s
    //   g g g g g s s s
    //   g g g g g s s s s
    //
    // g = gemm
    // s = syrk

    int jp = iter * partition->block_size;
    int i  = part_block_start_syrk(iter, me);
    int m  = part_block_size_syrk(iter, me);
    int j1 = jp + partition->block_size;
    int j2 = i;
    int n1 = j2 - j1;
    int k  = partition->block_size;

    double *A21;
    double *A21a, *A21b;
    double *A22;

    /* SYRK */
    if (m > 0) {
        A22 = &A(i,j2);
        A21 = &A(i,jp);
        cblas_dsyrk(CblasColMajor, CblasLower, CblasNoTrans,
                    m, k,
                    -1.0, A21, ldA,
                     1.0, A22, ldA);
    }

    /* GEMM */
    if (m > 0 && n1 > 0) {
        A22  = &A(i,j1);
        A21a = &A(i,jp);
        A21b = &A(j1,jp);
        cblas_dgemm(CblasColMajor, CblasNoTrans, CblasTrans,
                    m, n1, k,
                    -1.0, A21a, ldA,
                          A21b, ldA,
                     1.0, A22,  ldA);
    }
}

static void krnl_trsm(int iter, int me)
{
    int jp      = iter * partition->block_size;
    int i       = part_block_start_trsm(iter, me);
    int m       = part_block_size_trsm(iter, me);
    int n       = min(partition->matrix_size - jp, partition->block_size);
    double *A11 = &A(jp,jp);
    double *A21 = &A(i,jp);

    if (m > 0) {
        cblas_dtrsm(CblasColMajor, CblasRight, CblasLower, CblasTrans, CblasNonUnit,
                    m, n,
                    1.0, A11, ldA,
                         A21, ldA);
    }
}

void chol_task_par(void *ptr, int nth, int me)
{
    struct chol_task_arg *arg = (struct chol_task_arg*) ptr;

    n   = arg->n;
    A   = arg->A;
    ldA = arg->ldA;

    if (me == 0) {
        part_init(n, BLOCK_SIZE, nth);
    }

    if (me == 0) {
        // chol(0)
        krnl_chol(0);
    }

    // Synchronize.
    pcp_barrier_wait(barrier);

    // Loop over iterations.
    for (int iter = 0; iter < partition->num_iter; ++iter) {
        const int final_iteration = (iter == partition->num_iter - 1);

        // Synchronize
        pcp_barrier_wait(barrier);

        //
        // Phase I: trsm(iter) in //
        //

        krnl_trsm(iter, me);

        // Synchronize
        pcp_barrier_wait(barrier);

        //
        // Phase II: syrk(iter) in // plus chol(iter + 1)
        //

        krnl_syrk(iter, me);
        if (me == 0 && !final_iteration) {
            krnl_chol(iter + 1);
        }
    }

    // Synchronize
    pcp_barrier_wait(barrier);
}

max_stars_repo_name: NLAFET/pcp-runtime

max_issues_repo_name: NLAFET/pcp-runtime

max_forks_repo_name: NLAFET/pcp-runtime