text:
/* vector/gsl_vector_uint.h
 *
 * Copyright (C) 1996, 1997, 1998, 1999, 2000, 2007 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 3 of the License, or (at
 * your option) any later version.
 *
 * This program is distributed in the hope that it will be useful, but
 * WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 * General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this program; if not, write to the Free Software
 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
 */

#ifndef __GSL_VECTOR_UINT_H__
#define __GSL_VECTOR_UINT_H__

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

#include <stdlib.h>
#include <gsl/gsl_types.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_inline.h>
#include <gsl/gsl_check_range.h>
#include <gsl/gsl_block_uint.h>

#undef __BEGIN_DECLS
#undef __END_DECLS
#ifdef __cplusplus
# define __BEGIN_DECLS extern "C" {
# define __END_DECLS }
#else
# define __BEGIN_DECLS /* empty */
# define __END_DECLS /* empty */
#endif

__BEGIN_DECLS

typedef struct
{
  size_t size;
  size_t stride;
  unsigned int *data;
  gsl_block_uint *block;
  int owner;
}
gsl_vector_uint;

typedef struct
{
  gsl_vector_uint vector;
} _gsl_vector_uint_view;

typedef _gsl_vector_uint_view gsl_vector_uint_view;

typedef struct
{
  gsl_vector_uint vector;
} _gsl_vector_uint_const_view;

typedef const _gsl_vector_uint_const_view gsl_vector_uint_const_view;

/* Allocation */

GSL_FUN gsl_vector_uint *gsl_vector_uint_alloc (const size_t n);
GSL_FUN gsl_vector_uint *gsl_vector_uint_calloc (const size_t n);

GSL_FUN gsl_vector_uint *gsl_vector_uint_alloc_from_block (gsl_block_uint * b,
                                                     const size_t offset,
                                                     const size_t n,
                                                     const size_t stride);

GSL_FUN gsl_vector_uint *gsl_vector_uint_alloc_from_vector (gsl_vector_uint * v,
                                                      const size_t offset,
                                                      const size_t n,
                                                      const size_t stride);

GSL_FUN void gsl_vector_uint_free (gsl_vector_uint * v);

/* Views */

GSL_FUN _gsl_vector_uint_view
gsl_vector_uint_view_array (unsigned int *v, size_t n);

GSL_FUN _gsl_vector_uint_view
gsl_vector_uint_view_array_with_stride (unsigned int *base,
                                         size_t stride,
                                         size_t n);

GSL_FUN _gsl_vector_uint_const_view
gsl_vector_uint_const_view_array (const unsigned int *v, size_t n);

GSL_FUN _gsl_vector_uint_const_view
gsl_vector_uint_const_view_array_with_stride (const unsigned int *base,
                                               size_t stride,
                                               size_t n);

GSL_FUN _gsl_vector_uint_view
gsl_vector_uint_subvector (gsl_vector_uint *v,
                            size_t i,
                            size_t n);

GSL_FUN _gsl_vector_uint_view
gsl_vector_uint_subvector_with_stride (gsl_vector_uint *v,
                                        size_t i,
                                        size_t stride,
                                        size_t n);

GSL_FUN _gsl_vector_uint_const_view
gsl_vector_uint_const_subvector (const gsl_vector_uint *v,
                                  size_t i,
                                  size_t n);

GSL_FUN _gsl_vector_uint_const_view
gsl_vector_uint_const_subvector_with_stride (const gsl_vector_uint *v,
                                              size_t i,
                                              size_t stride,
                                              size_t n);

/* Operations */

GSL_FUN void gsl_vector_uint_set_zero (gsl_vector_uint * v);
GSL_FUN void gsl_vector_uint_set_all (gsl_vector_uint * v, unsigned int x);
GSL_FUN int gsl_vector_uint_set_basis (gsl_vector_uint * v, size_t i);

GSL_FUN int gsl_vector_uint_fread (FILE * stream, gsl_vector_uint * v);
GSL_FUN int gsl_vector_uint_fwrite (FILE * stream, const gsl_vector_uint * v);
GSL_FUN int gsl_vector_uint_fscanf (FILE * stream, gsl_vector_uint * v);
GSL_FUN int gsl_vector_uint_fprintf (FILE * stream, const gsl_vector_uint * v,
                              const char *format);

GSL_FUN int gsl_vector_uint_memcpy (gsl_vector_uint * dest, const gsl_vector_uint * src);

GSL_FUN int gsl_vector_uint_reverse (gsl_vector_uint * v);

GSL_FUN int gsl_vector_uint_swap (gsl_vector_uint * v, gsl_vector_uint * w);
GSL_FUN int gsl_vector_uint_swap_elements (gsl_vector_uint * v, const size_t i, const size_t j);

GSL_FUN unsigned int gsl_vector_uint_max (const gsl_vector_uint * v);
GSL_FUN unsigned int gsl_vector_uint_min (const gsl_vector_uint * v);
GSL_FUN void gsl_vector_uint_minmax (const gsl_vector_uint * v, unsigned int * min_out, unsigned int * max_out);

GSL_FUN size_t gsl_vector_uint_max_index (const gsl_vector_uint * v);
GSL_FUN size_t gsl_vector_uint_min_index (const gsl_vector_uint * v);
GSL_FUN void gsl_vector_uint_minmax_index (const gsl_vector_uint * v, size_t * imin, size_t * imax);

GSL_FUN int gsl_vector_uint_add (gsl_vector_uint * a, const gsl_vector_uint * b);
GSL_FUN int gsl_vector_uint_sub (gsl_vector_uint * a, const gsl_vector_uint * b);
GSL_FUN int gsl_vector_uint_mul (gsl_vector_uint * a, const gsl_vector_uint * b);
GSL_FUN int gsl_vector_uint_div (gsl_vector_uint * a, const gsl_vector_uint * b);
GSL_FUN int gsl_vector_uint_scale (gsl_vector_uint * a, const unsigned int x);
GSL_FUN int gsl_vector_uint_add_constant (gsl_vector_uint * a, const double x);
GSL_FUN int gsl_vector_uint_axpby (const unsigned int alpha, const gsl_vector_uint * x, const unsigned int beta, gsl_vector_uint * y);
GSL_FUN unsigned int gsl_vector_uint_sum (const gsl_vector_uint * a);

GSL_FUN int gsl_vector_uint_equal (const gsl_vector_uint * u,
                            const gsl_vector_uint * v);

GSL_FUN int gsl_vector_uint_isnull (const gsl_vector_uint * v);
GSL_FUN int gsl_vector_uint_ispos (const gsl_vector_uint * v);
GSL_FUN int gsl_vector_uint_isneg (const gsl_vector_uint * v);
GSL_FUN int gsl_vector_uint_isnonneg (const gsl_vector_uint * v);

GSL_FUN INLINE_DECL unsigned int gsl_vector_uint_get (const gsl_vector_uint * v, const size_t i);
GSL_FUN INLINE_DECL void gsl_vector_uint_set (gsl_vector_uint * v, const size_t i, unsigned int x);
GSL_FUN INLINE_DECL unsigned int * gsl_vector_uint_ptr (gsl_vector_uint * v, const size_t i);
GSL_FUN INLINE_DECL const unsigned int * gsl_vector_uint_const_ptr (const gsl_vector_uint * v, const size_t i);

#ifdef HAVE_INLINE

INLINE_FUN
unsigned int
gsl_vector_uint_get (const gsl_vector_uint * v, const size_t i)
{
#if GSL_RANGE_CHECK
  if (GSL_RANGE_COND(i >= v->size))
    {
      GSL_ERROR_VAL ("index out of range", GSL_EINVAL, 0);
    }
#endif
  return v->data[i * v->stride];
}

INLINE_FUN
void
gsl_vector_uint_set (gsl_vector_uint * v, const size_t i, unsigned int x)
{
#if GSL_RANGE_CHECK
  if (GSL_RANGE_COND(i >= v->size))
    {
      GSL_ERROR_VOID ("index out of range", GSL_EINVAL);
    }
#endif
  v->data[i * v->stride] = x;
}

INLINE_FUN
unsigned int *
gsl_vector_uint_ptr (gsl_vector_uint * v, const size_t i)
{
#if GSL_RANGE_CHECK
  if (GSL_RANGE_COND(i >= v->size))
    {
      GSL_ERROR_NULL ("index out of range", GSL_EINVAL);
    }
#endif
  return (unsigned int *) (v->data + i * v->stride);
}

INLINE_FUN
const unsigned int *
gsl_vector_uint_const_ptr (const gsl_vector_uint * v, const size_t i)
{
#if GSL_RANGE_CHECK
  if (GSL_RANGE_COND(i >= v->size))
    {
      GSL_ERROR_NULL ("index out of range", GSL_EINVAL);
    }
#endif
  return (const unsigned int *) (v->data + i * v->stride);
}
#endif /* HAVE_INLINE */

__END_DECLS

#endif /* __GSL_VECTOR_UINT_H__ */


text:
/* This file is necessary otherwise clang instantiate the bindgen tools
 * multiple time concurrently which produce non-deterministic results. */

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

ext: h

max_stars_repo_name: florv/scala-offheap

max_issues_repo_name: florv/scala-offheap

max_forks_repo_name: florv/scala-offheap


text:
/* ===============================================================================*/
/*   Version 1.0.             Cullan Howlett                                      */
/*   Copyright (c) 2017       International Centre for Radio Astronomy Research,  */
/*   The MIT License (MIT)    University of Western Australia                     */
/*                                                                                */
/* Permission is hereby granted, free of charge, to any person obtaining a copy   */
/* of this software and associated documentation files (the "Software"), to deal  */
/* in the Software without restriction, including without limitation the rights   */
/* to use, copy, modify, merge, publish, distribute, sublicense, and/or sell      */
/* copies of the Software, and to permit persons to whom the Software is          */
/* furnished to do so, subject to the following conditions:                       */
/*                                                                                */
/* The above copyright notice and this permission notice shall be included in     */
/* all copies or substantial portions of the Software.                            */
/*                                                                                */
/* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR     */
/* IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,       */
/* FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE    */
/* AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER         */
/* LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,  */
/* OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN      */
/* THE SOFTWARE.                                                                  */
/* ===============================================================================*/

#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_spline.h>
#include <gsl/gsl_linalg.h>
#include <gsl/gsl_integration.h>

// Fisher matrix calculation for surveys with velocity and density field measurements.
// ASSUMPTIONS:
//  -  Uncorrelated shot-noise between the density and velocity fields
//  -  I use the trapezium rule to integrate over r. There will be some error due to this, but this makes the most sense as
//     we are binning the number density anyway, and it makes hella difference in the speed of the code.
//  -  The redshift dependence of the non-linear matter and velocity divergence power spectra is captured using linear interpolation.
//  -  The PV error scales as a fixed pecentage of H0*r.
//  -  Flat LCDM cosmology (but not necessarily GR as gammaval can be changed).
//  -  The damping of the velocity and density fields due to non-linear RSD is redshift independent

// The parameters necessary for the calculation
static int nparams = 4;           // The number of free parameters (we can use any of beta, fsigma8, r_g, sigma_g, sigma_u)
static int Data[4] = {0,1,3,4};   // A vector of flags for the parameters we are interested in (0=beta, 1=fsigma8, 2=r_g, 3=sigma_g, 4=sigma_u). MAKE SURE THE LENGTH OF THIS VECTOR, NPARAMS AND THE ENTRIES AGREE/MAKE SENSE, OR YOU MIGHT GET NONSENSE RESULTS!!
static int nziter = 1;            // Now many bins in redshift between zmin and zmax we are considering
static double zmin = 0.0;         // The minimum redshift to consider (You must have power spectra that are within this range or GSL spline will error out)
static double zmax = 0.1;         // The maximum redshift to consider (You must have power spectra that are within this range or GSL spline will error out)
static double Om = 0.3121;        // The matter density at z=0
static double c = 299792.458;     // The speed of light in km/s
static double gammaval = 0.55;    // The value of gammaval to use in the forecasts (where f(z) = Om(z)^gammaval)
static double r_g = 1.0;          // The cross correlation coefficient between the velocity and density fields
static double beta0 = 0.393;      // The value of beta (at z=0, we'll modify this by the redshift dependent value of bias and f as required)
static double sigma80 = 0.8150;   // The value of sigma8 at z=0
static double sigma_u = 13.00;    // The value of the velocity damping parameter in Mpc/h. I use the values from Jun Koda's paper
static double sigma_g = 4.24;     // The value of the density damping parameter in Mpc/h. I use the values from Jun Koda's paper
static double kmax = 0.2;         // The maximum k to evaluate for dd, dv and vv correlations (Typical values are 0.1 - 0.2, on smaller scales the models are likely to break down).
static double survey_area[3] = {0.0, 0.0, 1.745};   // We need to know the survey area for each survey and the overlap area between the surveys (redshift survey only first, then PV survey only, then overlap.
                                                    // For fully overlapping we would have {0, 0, size_overlap}. For redshift larger than PV, we would have {size_red-size_overlap, 0, size_overlap}). Units are pi steradians, such that full sky is 4.0, half sky is 2.0 etc.
static double error_rand = 300.0;    // The observational error due to random non-linear velocities (I normally use 300km/s as in Jun Koda's paper)
static double error_dist = 0.05;     // The percentage error on the distance indicator (Typically 0.05 - 0.10 for SNe IA, 0.2 or more for Tully-Fisher or Fundamental Plane)
static double verbosity = 0;         // How much output to give: 0 = only percentage errors on fsigma8, 1 = other useful info and nuisance parameters, 2 = full fisher and covariance matrices

// The number of redshifts and the redshifts themselves of the input matter and velocity divergence power spectra.
// These numbers are multiplied by 100, converted to ints and written in the form _z0p%02d which is then appended to the filename Pvel_file. See routine read_power.
static double nzin = 11;
static double zin[11] = {0.00, 0.05, 0.10, 0.15, 0.20, 0.25, 0.30, 0.35, 0.40, 0.45, 0.50};
char * Pvel_file = "./example_files/example_pk";                                                  // The file containing the velocity divergence power spectrum. Don't include .dat as we'll append the redshifts on read in

// The files containing the number density of the surveys. First is the PV survey, then the redshift survey. These files MUST have the same binning and redshift range,
// so that the sum over redshift bins works (would be fine if we used splines), i.e., if one survey is shallower then that file must contain rows with n(z)=0.
// I also typically save nbar x 10^6 in the input file to make sure I don't lose precision when outputting small nbar values to files. This is corrected when the nbar file
// is read in, so see the read_nz() routine!
char * nbar_file[300] = {"./example_files/example_nbar_vel.dat",
                         "./example_files/example_nbar_red.dat"};

// Other global parameters and arrays
int NK, * NRED;
double pkkmin;         // The minimum kmin to integrate over, based on the input power spectrum file
double pkkmax;         // The maximum k in the input power spectrum. The maximum k to integrate over is the smallest of this or kmax
double * zarray;
double * rarray;
double * deltararray;
double * growtharray;
double ** nbararray;
double * karray, * deltakarray;
double ** pmmarray, ** pmtarray, ** pttarray;
gsl_spline * growth_spline, * r_spline;
gsl_interp_accel * growth_acc, * r_acc;

// Prototypes
double zeff_integrand(double mu, void * pin);
double mu_integrand(double mu, void * pin);
double ezinv(double x, void *p);
double rz(double red);
double growthfunc(double x, void *p);
double growthz(double red);
void read_nz();
void read_power();

// Calculates the fished matrix for a velocity survey.
int main(int argc, char **argv) {

    FILE * fout;
    int i, j;

    // Read in the velocity divergence power spectrum output from the COPTER code (Carlson 2009)
    read_power();

    // Read in the number densities of the surveys
    read_nz();

    // Run some checks
    if (!((survey_area[0] > 0.0) || (survey_area[2] > 0.0))) {
        for (i=0; i<nparams; i++) {
            if (Data[i] == 2) {
                printf("ERROR: r_g is a free parameter, but there is no information in the density field (Fisher matrix will be singular)\n");
                exit(0);
            }
            if (Data[i] == 3) {
                printf("ERROR: sigma_g is a free parameter, but there is no information in the density field (Fisher matrix will be singular)\n");
                exit(0);
            }
        }
    }
    if (!((survey_area[1] > 0.0) || (survey_area[2] > 0.0))) {
        for (i=0; i<nparams; i++) {
            if (Data[i] == 4) {
                printf("ERROR: sigma_u is a free parameter, but there is no information in the velocity field (Fisher matrix will be singular)\n");
                exit(0);
            }
        }
    }
    if ((sizeof(Data)/sizeof(*Data)) != nparams) {
        printf("ERROR: Size of Data vector for parameters of interest must be equal to nparams\n");
        exit(0);
    }

    // Calculate the Fisher matrices for all bins.
    gsl_matrix * Fisher_Tot = gsl_matrix_alloc(nparams, nparams);
    for (i=0; i<nparams; i++) {
        for (j=0; j<nparams; j++) gsl_matrix_set(Fisher_Tot, i, j, 0.0);
    }

    printf("Evaluating the Fisher Matrix for %d bins between [z_min = %lf, z_max = %lf]\n", nziter, zmin, zmax);

    if (verbosity == 0) printf("#     zmin         zmax         zeff      fsigma8(z_eff)   percentage error(z_eff)\n");

    int ziter;
    for (ziter = 0; ziter<nziter; ziter++) {

        double zbinwidth = (zmax-zmin)/(nziter);
        double zmin_iter = ziter*zbinwidth + zmin;
        double zmax_iter = (ziter+1.0)*zbinwidth + zmin;

        double rzmax = gsl_spline_eval(r_spline, zmax_iter, r_acc);
        double kmin = M_PI/rzmax;

        if (verbosity > 0) printf("Evaluating the Fisher Matrix for [k_min = %lf, k_max = %lf] and [z_min = %lf, z_max = %lf]\n", kmin, kmax, zmin_iter, zmax_iter);

        // Calculate the effective redshift (which I base on the sum of the S/N for the density and velocity fields)
        int numk;
        double k_sum1 = 0.0, k_sum2 = 0.0;
        for (numk=0; numk<NK; numk++) {

            double k = karray[numk]+0.5*deltakarray[numk];
            double deltak = deltakarray[numk];
            if (k < kmin) continue;
            if (k > kmax) continue;

            double result, error;
            double params[4] = {numk, k, zmin_iter, zmax_iter};

            size_t nevals = 1000;
            gsl_function F;
            F.function = &zeff_integrand;
            F.params = &params;
            gsl_integration_workspace * w = gsl_integration_workspace_alloc(1000);

            gsl_integration_qags(&F, 0.0, 1.0, 0, 5e-3, nevals, w, &result, &error);
            gsl_integration_workspace_free(w);

            k_sum1 += k*k*deltak*result;
            k_sum2 += k*k*deltak;

        }
        double z_eff = k_sum1/k_sum2;
        if (verbosity > 0) printf("Effective redshift z_eff = %lf\n", z_eff);

        double growth_eff = gsl_spline_eval(growth_spline, z_eff, growth_acc);

        // Calculate the fisher matrix, integrating over k, then mu, then r (r is last as it means we are effectively integrating over effective volume).
        // As the input spectra are tabulated we'll just use the trapezium rule to integrate over k
        gsl_matrix * Fisher = gsl_matrix_alloc(nparams, nparams);
        for (i=0; i<nparams; i++) {
            for (j=i; j<nparams; j++) {

                double k_sum = 0.0;
                for (numk=0; numk<NK; numk++) {

                    double k = karray[numk]+0.5*deltakarray[numk];
                    double deltak = deltakarray[numk];
                    if (k < kmin) continue;
                    if (k > kmax) continue;

                    double result, error;
                    double params[6] = {numk, k, Data[i], Data[j], zmin_iter, zmax_iter};

                    size_t nevals = 1000;
                    gsl_function F;
                    F.function = &mu_integrand;
                    F.params = &params;
                    gsl_integration_workspace * w = gsl_integration_workspace_alloc(1000);

                    gsl_integration_qags(&F, 0.0, 1.0, 0, 5e-3, nevals, w, &result, &error);
                    gsl_integration_workspace_free(w);

                    k_sum += k*k*deltak*result;

                }
                //printf("%d, %d, %lf\n", i, j, k_sum/(4.0*M_PI));
                gsl_matrix_set(Fisher, i, j, k_sum/(4.0*M_PI));
                gsl_matrix_set(Fisher, j, i, k_sum/(4.0*M_PI));

            }
        }

        for (i=0; i<nparams; i++) {
            for (j=0; j<nparams; j++) {
                double val = gsl_matrix_get(Fisher_Tot, i, j) + gsl_matrix_get(Fisher, i, j);
                gsl_matrix_set(Fisher_Tot, i, j, val);
            }
        }

        if (verbosity == 2) {
            printf("Fisher Matrix\n======================\n");
            for (i=0; i<nparams; i++) {
                printf("[");
                for (j=0; j<nparams; j++) printf("%15.6lf,\t", gsl_matrix_get(Fisher, i, j));
                printf("],\n");
            }
        }

        // Now invert the Fisher matrix
        int s;
        gsl_permutation * p;
        gsl_matrix * Covariance = gsl_matrix_alloc(nparams, nparams);
        p = gsl_permutation_alloc(nparams);
        gsl_linalg_LU_decomp(Fisher, p, &s);
        gsl_linalg_LU_invert(Fisher, p, Covariance);
        gsl_permutation_free(p);

        double sigma8 = sigma80 * growth_eff;
        double Omz = Om*ezinv(z_eff,NULL)*ezinv(z_eff,NULL)*(1.0+z_eff)*(1.0+z_eff)*(1.0+z_eff);
        double f = pow(Omz, gammaval);
        double beta = f*beta0*growth_eff/pow(Om,0.55);

        if (verbosity == 0) {
            for (i=0; i<nparams; i++) {
                if (Data[i] == 1) printf("%12.6lf  %12.6lf  %12.6lf  %12.6lf  %12.6lf\n", zmin_iter, zmax_iter, z_eff, f*sigma8, 100.0*sqrt(gsl_matrix_get(Covariance, i, i))/(f*sigma8));
            }
        }

        if (verbosity > 0) {
            for (i=0; i<nparams; i++) {
                if (Data[i] == 0) {
                    printf("beta = %12.6lf +/- %12.6lf\n", beta, sqrt(gsl_matrix_get(Covariance, i, i)));
                    printf("%4.2lf percent error on beta\n", 100.0*sqrt(gsl_matrix_get(Covariance, i, i))/beta);
                }
                if (Data[i] == 1) {
                    printf("fsigma8 = %12.6lf +/- %12.6lf\n", f*sigma8, sqrt(gsl_matrix_get(Covariance, i, i)));
                    printf("%4.2lf percent error on fsigma8\n", 100.0*sqrt(gsl_matrix_get(Covariance, i, i))/(f*sigma8));
                }
                if (Data[i] == 2) {
                    printf("r_g = %12.6lf +/- %12.6lf\n", r_g, sqrt(gsl_matrix_get(Covariance, i, i)));
                    printf("%4.2lf percent error on r_g\n", 100.0*sqrt(gsl_matrix_get(Covariance, i, i))/r_g);
                }
                if (Data[i] == 3) {
                    printf("sigma_g = %12.6lf +/- %12.6lf\n", sigma_g, sqrt(gsl_matrix_get(Covariance, i, i)));
                    printf("%4.2lf percent error on sigma_g\n", 100.0*sqrt(gsl_matrix_get(Covariance, i, i))/sigma_g);
                }
                if (Data[i] == 4) {
                    printf("sigma_u = %12.6lf +/- %12.6lf\n", sigma_u, sqrt(gsl_matrix_get(Covariance, i, i)));
                    printf("%4.2lf percent error on sigma_u\n", 100.0*sqrt(gsl_matrix_get(Covariance, i, i))/sigma_u);
                }
            }
        }

        if (verbosity == 2) {
            printf("Covariance Matrix\n======================\n");
            for (i=0; i<nparams; i++) {
                printf("[");
                for (j=0; j<nparams; j++) printf("%15.6lf,\t", gsl_matrix_get(Covariance, i, j));
                printf("],\n");
            }
        }

        gsl_matrix_free(Fisher);
        gsl_matrix_free(Covariance);
    }

    // Now the full Fisher matrix over all redshifts if we had more than 1 redshift bin
    if (nziter > 1) {
        double rzmax = gsl_spline_eval(r_spline, zmax, r_acc);
        double kmin = M_PI/rzmax;

        if (verbosity > 0) printf("Finally, evaluating the Fisher Matrix for [k_min = %lf, k_max = %lf] and [z_min = %lf, z_max = %lf]\n", kmin, kmax, zmin, zmax);

        // Calculate the effective redshift
        int numk;
        double k_sum1 = 0.0, k_sum2 = 0.0;
        for (numk=0; numk<NK; numk++) {

            double k = karray[numk]+0.5*deltakarray[numk];
            double deltak = deltakarray[numk];
            if (k < kmin) continue;
            if (k > kmax) continue;

            double result, error;
            double params[4] = {numk, k, zmin, zmax};

            size_t nevals = 1000;
            gsl_function F;
            F.function = &zeff_integrand;
            F.params = &params;
            gsl_integration_workspace * w = gsl_integration_workspace_alloc(1000);

            gsl_integration_qags(&F, 0.0, 1.0, 0, 5e-3, nevals, w, &result, &error);
            gsl_integration_workspace_free(w);

            k_sum1 += k*k*deltak*result;
            k_sum2 += k*k*deltak;

        }
        double z_eff = k_sum1/k_sum2;
        if (verbosity > 0) printf("Effective redshift z_eff = %lf\n", z_eff);

        double growth_eff = gsl_spline_eval(growth_spline, z_eff, growth_acc);

        if (verbosity == 2) {
            printf("Fisher Matrix\n======================\n");
            for (i=0; i<nparams; i++) {
                printf("[");
                for (j=0; j<nparams; j++) printf("%15.6lf,\t", gsl_matrix_get(Fisher_Tot, i, j));
                printf("],\n");
            }
        }

        // Now invert the Fisher matrix
        int s;
        gsl_permutation * p;
        gsl_matrix * Covariance = gsl_matrix_alloc(nparams, nparams);
        p = gsl_permutation_alloc(nparams);
        gsl_linalg_LU_decomp(Fisher_Tot, p, &s);
        gsl_linalg_LU_invert(Fisher_Tot, p, Covariance);
        gsl_permutation_free(p);

        double sigma8 = sigma80 * growth_eff;
        double Omz = Om*ezinv(z_eff,NULL)*ezinv(z_eff,NULL)*(1.0+z_eff)*(1.0+z_eff)*(1.0+z_eff);
        double f = pow(Omz, gammaval);
        double beta = f*beta0*growth_eff/pow(Om,0.55);

        if (verbosity == 0) {
            printf("# Full redshift range:\n");
            for (i=0; i<nparams; i++) {
                if (Data[i] == 1) printf("%12.6lf  %12.6lf  %12.6lf  %12.6lf  %12.6lf\n", zmin, zmax, z_eff, f*sigma8, 100.0*sqrt(gsl_matrix_get(Covariance, i, i))/(f*sigma8));
            }
        }

        if (verbosity > 0) {
            for (i=0; i<nparams; i++) {
                if (Data[i] == 0) {
                    printf("beta = %12.6lf +/- %12.6lf\n", beta, sqrt(gsl_matrix_get(Covariance, i, i)));
                    printf("%4.2lf percent error on beta\n", 100.0*sqrt(gsl_matrix_get(Covariance, i, i))/beta);
                }
                if (Data[i] == 1) {
                    printf("fsigma8 = %12.6lf +/- %12.6lf\n", f*sigma8, sqrt(gsl_matrix_get(Covariance, i, i)));
                    printf("%4.2lf percent error on fsigma8\n", 100.0*sqrt(gsl_matrix_get(Covariance, i, i))/(f*sigma8));
                }
                if (Data[i] == 2) {
                    printf("r_g = %12.6lf +/- %12.6lf\n", r_g, sqrt(gsl_matrix_get(Covariance, i, i)));
                    printf("%4.2lf percent error on r_g\n", 100.0*sqrt(gsl_matrix_get(Covariance, i, i))/r_g);
                }
                if (Data[i] == 3) {
                    printf("sigma_g = %12.6lf +/- %12.6lf\n", sigma_g, sqrt(gsl_matrix_get(Covariance, i, i)));
                    printf("%4.2lf percent error on sigma_g\n", 100.0*sqrt(gsl_matrix_get(Covariance, i, i))/sigma_g);
                }
                if (Data[i] == 4) {
                    printf("sigma_u = %12.6lf +/- %12.6lf\n", sigma_u, sqrt(gsl_matrix_get(Covariance, i, i)));
                    printf("%4.2lf percent error on sigma_u\n", 100.0*sqrt(gsl_matrix_get(Covariance, i, i))/sigma_u);
                }
            }
        }

        if (verbosity == 2) {
            printf("Covariance Matrix\n======================\n");
            for (i=0; i<nparams; i++) {
                printf("[");
                for (j=0; j<nparams; j++) printf("%15.6lf,\t", gsl_matrix_get(Covariance, i, j));
                printf("],\n");
            }
        }
    }

    gsl_matrix_free(Fisher_Tot);
    gsl_spline_free(growth_spline);
    gsl_interp_accel_free(growth_acc);

    return 0;
}

// The integrand to calculate the effective redshift. I'm not actually sure how this is done in the case of
// density and velocity field measurements, but it seems logical to base it on the sum of the integral of the density spectra and the velocity power spectra
// weighted by their effective signal to noise. In this way the effective redshift is calculated in the same way as a redshift survey, but there is some dependence
// on the S/N in the velocity power spectrum too. In any case, the S/N of the density field measurement is always much higher (because there are
// no errors and the number density is higher) and so this dominates the effective redshift calculation.
double zeff_integrand(double mu, void * pin) {

    int i, j, m, q, u, surv;
    double * p = (double *)pin;

    int numk = (int)p[0];
    double k = p[1];
    double zminval = p[2];
    double zmaxval = p[3];
    gsl_interp_accel * Pmm_acc, * Pmt_acc, * Ptt_acc;
    gsl_spline * Pmm_spline, * Pmt_spline, * Ptt_spline;
    if (nzin > 1) {
        double * Pmm_array = (double *)malloc(nzin*sizeof(double));
        double * Pmt_array = (double *)malloc(nzin*sizeof(double));
        double * Ptt_array = (double *)malloc(nzin*sizeof(double));
        for (j=0; j<nzin; j++) {
            Pmm_array[j] = pmmarray[j][numk];
            Pmt_array[j] = pmtarray[j][numk];
            Ptt_array[j] = pttarray[j][numk];
        }
        Pmm_acc    = gsl_interp_accel_alloc();
        Pmm_spline = gsl_spline_alloc(gsl_interp_cspline, nzin);
        gsl_spline_init(Pmm_spline, zin, Pmm_array, nzin);
        free(Pmm_array);

        Pmt_acc    = gsl_interp_accel_alloc();
        Pmt_spline = gsl_spline_alloc(gsl_interp_cspline, nzin);
        gsl_spline_init(Pmt_spline, zin, Pmt_array, nzin);
        free(Pmt_array);

        Ptt_acc    = gsl_interp_accel_alloc();
        Ptt_spline = gsl_spline_alloc(gsl_interp_cspline, nzin);
        gsl_spline_init(Ptt_spline, zin, Ptt_array, nzin);
        free(Ptt_array);
    }

    double dendamp = sqrt(1.0/(1.0+0.5*(k*k*mu*mu*sigma_g*sigma_g)));     // This is unitless
    double veldamp = sin(k*sigma_u)/(k*sigma_u);                          // This is unitless

    double dVeff = 0.0, zdVeff = 0.0;
    for (i=0; i<NRED[0]; i++) {

        double zval = zarray[i];
        if (zval < zminval) continue;
        if (zval > zmaxval) break;

        double r_sum = 0.0;
        double r = rarray[i];
        double deltar = deltararray[i];

        double dd_prefac=0.0, vv_prefac=0.0;
        double P_gg=0.0, P_uu=0.0;

        double sigma8 = sigma80 * growtharray[i];

        // First lets calculate the relevant power spectra. Interpolate the power spectra linearly in redshift
        double Pmm, Pmt, Ptt;
        Pmm = gsl_spline_eval(Pmm_spline, zval, Pmm_acc);
        Pmt = gsl_spline_eval(Pmt_spline, zval, Pmm_acc);
        Ptt = gsl_spline_eval(Ptt_spline, zval, Pmm_acc);

        double Omz = Om*ezinv(zval,NULL)*ezinv(zval,NULL)*(1.0+zval)*(1.0+zval)*(1.0+zval);
        double f = pow(Omz, gammaval);
        double beta = f*beta0*growtharray[i]/pow(Om,0.55);

        vv_prefac  = 1.0e2*f*mu*veldamp/k;
        dd_prefac = (1.0/(beta*beta) + 2.0*r_g*mu*mu/beta + mu*mu*mu*mu)*f*f*dendamp*dendamp;
        P_gg = dd_prefac*Pmm;
        P_uu = vv_prefac*vv_prefac*Ptt;

        // We need to do the overlapping and non-overlapping parts of the redshifts and PV surveys separately
        for (surv=0; surv<3; surv++) {
            double surv_sum = 0.0;
            if (survey_area[surv] > 0.0) {
                double error_obs, error_noise, n_g = 0.0, n_u = 0.0;

                // Set the nbar for each section.
                if (surv == 0) {
                    n_g = nbararray[1][i];
                } else if (surv == 1) {
                    error_obs = 100.0*error_dist*r;                              // Percentage error * distance * H0 in km/s (factor of 100.0 comes from hubble parameter)
                    error_noise = error_rand*error_rand + error_obs*error_obs;   // Error_noise is in km^{2}s^{-2}
                    n_u = nbararray[0][i]/error_noise;
                } else {
                    error_obs = 100.0*error_dist*r;                              // Percentage error * distance * H0 in km/s (factor of 100.0 comes from hubble parameter)
                    error_noise = error_rand*error_rand + error_obs*error_obs;   // Error_noise is in km^{2}s^{-2}
                    n_u = nbararray[0][i]/error_noise;
                    n_g = nbararray[1][i];
                }

                double value1 = n_g/(1.0 + n_g*P_gg);
                double value2 = n_u/(1.0 + n_u*P_uu);
                surv_sum += value1*value1 + value2*value2;

                surv_sum *= survey_area[surv];
                r_sum += surv_sum;
            }
        }

        dVeff += r*r*deltar*r_sum;
        zdVeff += zval*r*r*deltar*r_sum;

    }

    gsl_spline_free(Pmm_spline);
    gsl_spline_free(Pmt_spline);
    gsl_spline_free(Ptt_spline);
    gsl_interp_accel_free(Pmm_acc);
    gsl_interp_accel_free(Pmt_acc);
    gsl_interp_accel_free(Ptt_acc);

    return zdVeff/dVeff;

}

// The integrand for the integral over mu in the Fisher matrix calculation.
// For each mu we need to create a 4x4 matrix of the relevant power spectra derivatives and the inverse of the power spectrum matrix.
// Because there are some regions where the number density goes to zero we have to work directly with the inverse as it is difficult to invert numerically
// but if we deal with the inverse only then we can just set the relevant parts to zero when the number density is zero.
double mu_integrand(double mu, void * pin) {

    int i, j, m, q, u, surv;
    double * p = (double *)pin;
    double result, error;

    int numk = (int)p[0];
    double k = p[1];
    double zminval = p[4];
    double zmaxval = p[5];
    gsl_interp_accel * Pmm_acc, * Pmt_acc, * Ptt_acc;
    gsl_spline * Pmm_spline, * Pmt_spline, * Ptt_spline;
    double * Pmm_array = (double *)malloc(nzin*sizeof(double));
    double * Pmt_array = (double *)malloc(nzin*sizeof(double));
    double * Ptt_array = (double *)malloc(nzin*sizeof(double));
    for (j=0; j<nzin; j++) {
        Pmm_array[j] = pmmarray[j][numk];
        Pmt_array[j] = pmtarray[j][numk];
        Ptt_array[j] = pttarray[j][numk];
    }
    Pmm_acc    = gsl_interp_accel_alloc();
    Pmm_spline = gsl_spline_alloc(gsl_interp_cspline, nzin);
    gsl_spline_init(Pmm_spline, zin, Pmm_array, nzin);
    free(Pmm_array);

    Pmt_acc    = gsl_interp_accel_alloc();
    Pmt_spline = gsl_spline_alloc(gsl_interp_cspline, nzin);
    gsl_spline_init(Pmt_spline, zin, Pmt_array, nzin);
    free(Pmt_array);

    Ptt_acc    = gsl_interp_accel_alloc();
    Ptt_spline = gsl_spline_alloc(gsl_interp_cspline, nzin);
    gsl_spline_init(Ptt_spline, zin, Ptt_array, nzin);
    free(Ptt_array);

    double dendamp = sqrt(1.0/(1.0+0.5*(k*k*mu*mu*sigma_g*sigma_g)));     // This is unitless
    double veldamp = sin(k*sigma_u)/(k*sigma_u);                          // This is unitless

    double result_sum = 0.0;
    for (i=0; i<NRED[0]; i++) {

        double zval = zarray[i];
        double r_sum = 0.0;
        double r = rarray[i];
        double deltar = deltararray[i];

        if (zval < zminval) continue;
        if (zval > zmaxval) break;

        double dd_prefac=0.0, dv_prefac=0.0, vv_prefac=0.0;
        double P_gg=0.0, P_ug=0.0, P_uu=0.0;

        double sigma8 = sigma80 * growtharray[i];

        // First lets calculate the relevant power spectra. Interpolate the power spectra linearly in redshift
        double Pmm, Pmt, Ptt;
        Pmm = gsl_spline_eval(Pmm_spline, zval, Pmm_acc);
        Pmt = gsl_spline_eval(Pmt_spline, zval, Pmm_acc);
        Ptt = gsl_spline_eval(Ptt_spline, zval, Pmm_acc);

        double Omz = Om*ezinv(zval,NULL)*ezinv(zval,NULL)*(1.0+zval)*(1.0+zval)*(1.0+zval);
        double f = pow(Omz, gammaval);
        double beta = f*beta0*growtharray[i]/pow(Om,0.55);

        vv_prefac  = 1.0e2*f*mu*veldamp/k;
        dd_prefac = (1.0/(beta*beta) + 2.0*r_g*mu*mu/beta + mu*mu*mu*mu)*f*f*dendamp*dendamp;
        dv_prefac = (r_g/beta + mu*mu)*f*dendamp;
        P_gg = dd_prefac*Pmm;
        P_ug = vv_prefac*dv_prefac*Pmt;
        P_uu = vv_prefac*vv_prefac*Ptt;

        // And now the derivatives. Need to create a matrix of derivatives for each of the two parameters of interest
        gsl_matrix * dPdt1 = gsl_matrix_calloc(2, 2);
        gsl_matrix * dPdt2 = gsl_matrix_calloc(2, 2);
        double value;
        switch((int)p[2]) {
            // Differential w.r.t betaA
            case 0:
                value = -2.0*(1.0/beta + r_g*mu*mu)*f*f*dendamp*dendamp*Pmm/(beta*beta);
                gsl_matrix_set(dPdt1, 0, 0, value);
                value = -(vv_prefac*f*r_g*dendamp*Pmt)/(beta*beta);
                gsl_matrix_set(dPdt1, 0, 1, value);
                gsl_matrix_set(dPdt1, 1, 0, value);
                break;
            // Differential w.r.t fsigma8
            case 1:
                value = 2.0*(f/(beta*beta) + 2.0*f*r_g*mu*mu/beta + f*mu*mu*mu*mu)*dendamp*dendamp*Pmm/sigma8;
                gsl_matrix_set(dPdt1, 0, 0, value);
                value = 2.0*vv_prefac*(r_g/beta + mu*mu)*dendamp*Pmt/sigma8;
                gsl_matrix_set(dPdt1, 0, 1, value);
                gsl_matrix_set(dPdt1, 1, 0, value);
                value = (2.0*P_uu)/(f*sigma8);
                gsl_matrix_set(dPdt1, 1, 1, value);
                break;
            // Differential w.r.t r_g
            case 2:
                value = 2.0*(1.0/beta)*mu*mu*f*f*dendamp*dendamp*Pmm;
                gsl_matrix_set(dPdt1, 0, 0, value);
                value = vv_prefac*(1.0/beta)*f*dendamp*Pmt;
                gsl_matrix_set(dPdt1, 0, 1, value);
                gsl_matrix_set(dPdt1, 1, 0, value);
                break;
            // Differential w.r.t sigma_g
            case 3:
                value = -k*k*mu*mu*dendamp*dendamp*sigma_g*P_gg;
                gsl_matrix_set(dPdt1, 0, 0, value);
                value = -0.5*k*k*mu*mu*dendamp*dendamp*sigma_g*P_ug;
                gsl_matrix_set(dPdt1, 0, 1, value);
                gsl_matrix_set(dPdt1, 1, 0, value);
                break;
            // Differential w.r.t sigma_u
            case 4:
                value = P_ug*(k*cos(k*sigma_u)/sin(k*sigma_u) - 1.0/sigma_u);
                gsl_matrix_set(dPdt1, 0, 1, value);
                gsl_matrix_set(dPdt1, 1, 0, value);
                value = 2.0*P_uu*(k*cos(k*sigma_u)/sin(k*sigma_u) - 1.0/sigma_u);
                gsl_matrix_set(dPdt1, 1, 1, value);
                break;
            default:
                break;
        }
        switch((int)p[3]) {
            // Differential w.r.t betaA
            case 0:
                value = -2.0*(1.0/beta + r_g*mu*mu)*f*f*dendamp*dendamp*Pmm/(beta*beta);
                gsl_matrix_set(dPdt2, 0, 0, value);
                value = -(vv_prefac*f*r_g*dendamp*Pmt)/(beta*beta);
                gsl_matrix_set(dPdt2, 0, 1, value);
                gsl_matrix_set(dPdt2, 1, 0, value);
                break;
            // Differential w.r.t fsigma8
            case 1:
                value = 2.0*(f/(beta*beta) + 2.0*f*r_g*mu*mu/beta + f*mu*mu*mu*mu)*dendamp*dendamp*Pmm/sigma8;
                gsl_matrix_set(dPdt2, 0, 0, value);
                value = 2.0*vv_prefac*(r_g/beta + mu*mu)*dendamp*Pmt/sigma8;
                gsl_matrix_set(dPdt2, 0, 1, value);
                gsl_matrix_set(dPdt2, 1, 0, value);
                value = (2.0*P_uu)/(f*sigma8);
                gsl_matrix_set(dPdt2, 1, 1, value);
                break;
            // Differential w.r.t r_g
            case 2:
                value = 2.0*(1.0/beta)*mu*mu*f*f*dendamp*dendamp*Pmm;
                gsl_matrix_set(dPdt2, 0, 0, value);
                value = vv_prefac*(1.0/beta)*f*dendamp*Pmt;
                gsl_matrix_set(dPdt2, 0, 1, value);
                gsl_matrix_set(dPdt2, 1, 0, value);
                break;
            // Differential w.r.t sigma_g
            case 3:
                value = -k*k*mu*mu*dendamp*dendamp*sigma_g*P_gg;
                gsl_matrix_set(dPdt2, 0, 0, value);
                value = -0.5*k*k*mu*mu*dendamp*dendamp*sigma_g*P_ug;
                gsl_matrix_set(dPdt2, 0, 1, value);
                gsl_matrix_set(dPdt2, 1, 0, value);
                break;
            // Differential w.r.t sigma_u
            case 4:
                value = P_ug*(k*cos(k*sigma_u)/sin(k*sigma_u) - 1.0/sigma_u);
                gsl_matrix_set(dPdt2, 0, 1, value);
                gsl_matrix_set(dPdt2, 1, 0, value);
                value = 2.0*P_uu*(k*cos(k*sigma_u)/sin(k*sigma_u) - 1.0/sigma_u);
                gsl_matrix_set(dPdt2, 1, 1, value);
                break;
            default:
                break;
        }

        // We need to do the overlapping and non-overlapping parts of the surveys separately
        for (surv=0; surv<3; surv++) {
            double surv_sum = 0.0;
            if (survey_area[surv] > 0.0) {
                double error_obs, error_noise, n_g = 0.0, n_u = 0.0;

                // Set the nbar for each section.
                if (surv == 0) {
                    n_g = nbararray[1][i];
                } else if (surv == 1) {
                    error_obs = 100.0*error_dist*r;                              // Percentage error * distance * H0 in km/s (factor of 100.0 comes from hubble parameter)
                    error_noise = error_rand*error_rand + error_obs*error_obs;   // Error_noise is in km^{2}s^{-2}
                    n_u = nbararray[0][i]/error_noise;
                } else {
                    error_obs = 100.0*error_dist*r;                              // Percentage error * distance * H0 in km/s (factor of 100.0 comes from hubble parameter)
                    error_noise = error_rand*error_rand + error_obs*error_obs;   // Error_noise is in km^{2}s^{-2}
                    n_u = nbararray[0][i]/error_noise;
                    n_g = nbararray[1][i];
                }

                //printf("%lf, %lf, %lf\n", r, n_g, 1.0e6*n_u);

                if (!((n_u > 0.0) || (n_g > 0.0))) continue;

                // First we need the determinant.
                double det = 1.0 + n_u*n_g*(P_gg*P_uu - P_ug*P_ug) + n_u*P_uu + n_g*P_gg;

                // Now the inverse matrix.
                gsl_matrix * iP = gsl_matrix_calloc(2, 2);
                value = n_u*n_g*P_uu + n_g;
                gsl_matrix_set(iP, 0, 0, value);
                value = n_g*n_u*P_gg + n_u;
                gsl_matrix_set(iP, 1, 1, value);
                value = - n_g*n_u*P_ug;
                gsl_matrix_set(iP, 0, 1, value);
                gsl_matrix_set(iP, 1, 0, value);

                // Finally we need to compute the Fisher integrand by summing over the inverse and differential matrices
                for (j=0; j<2; j++) {
                    for (m=0; m<2; m++) {
                        for (u=0; u<2; u++) {
                            for (q=0; q<2; q++) {
                                value = gsl_matrix_get(dPdt1, j, q)*gsl_matrix_get(iP, q, u)*gsl_matrix_get(dPdt2, u, m)*gsl_matrix_get(iP, m, j);
                                surv_sum += value;
                            }
                        }
                    }
                }
                surv_sum /= det*det;
                surv_sum *= survey_area[surv];
                r_sum += surv_sum;
                gsl_matrix_free(iP);
                //printf("%d, %lf, %lf, %lf, %lf\n", surv, k, mu, r, r_sum);

            }
        }
        //printf("%lf, %lf, %lf, %lf\n", k, mu, r, r_sum);

        result_sum += r*r*deltar*r_sum;

        gsl_matrix_free(dPdt1);
        gsl_matrix_free(dPdt2);
    }

    gsl_spline_free(Pmm_spline);
    gsl_spline_free(Pmt_spline);
    gsl_spline_free(Ptt_spline);
    gsl_interp_accel_free(Pmm_acc);
    gsl_interp_accel_free(Pmt_acc);
    gsl_interp_accel_free(Ptt_acc);

    return result_sum;
}

// Routine to read in the number density as a function of redshift. We need a file containing the left-most edge of each redshift bin and teh number density in that bin.
// From this we create arrays to store the bin centre, the bin width, the comoving distance and growth factor at the bin centre and the number density.
// The last bin width and bin centre is constructed from the last row of the input and the value of zmax at the top of the code.
// ITS VERY IMPORTANT THAT THE NUMBER OF ROWS AND THE REDSHIFTS OF BOTH THE DENSITY AND PV NUMBER DENSITIES MATCH AS THE INTEGRATION OVER Z IS DONE USING THE TRAPEZIUM RULE.
// ALSO MAKE NOTE OF THE FACTOR OF 1.0e-6 ON LINE 827. THIS IS BECAUSE I TYPICALLY SAVE THE VALUE OF NBAR x 10^6 IN THE INPUT FILES< SO THAT I DON'T LOSE PRECISION
// WHEN SMALL VALUES OF THE NUMBER DENSITY ARE WRITTEN TO A FILE!
void read_nz() {

    FILE * fp;
    char buf[500];
    int i, nsamp;

    NRED = (int *)calloc(2, sizeof(int));
    nbararray = (double **)calloc(2, sizeof(double*));
    double * zinarray;

    for (nsamp = 0; nsamp < 2; nsamp++) {

        if(!(fp = fopen(nbar_file[nsamp], "r"))) {
            printf("\nERROR: Can't open nbar file '%s'.\n\n", nbar_file[nsamp]);
            exit(0);
        }

        NRED[nsamp] = 0;
        while(fgets(buf,500,fp)) {
            if(strncmp(buf,"#",1)!=0) {
                double tz, tnbar;
                if(sscanf(buf, "%lf %lf\n", &tz, &tnbar) != 2) {printf("nbar read error\n"); exit(0);};
                if (tz > zmax) break;
                NRED[nsamp]++;
            }
        }
        fclose(fp);

        if (nsamp == 0) zinarray = (double *)calloc(NRED[nsamp], sizeof(double));
        nbararray[nsamp] = (double *)calloc(NRED[nsamp], sizeof(double));

        NRED[nsamp] = 0;
        fp = fopen(nbar_file[nsamp], "r");
        while(fgets(buf,500,fp)) {
            if(strncmp(buf,"#",1)!=0) {
                double tz, tnbar;
                if(sscanf(buf, "%lf %lf\n", &tz, &tnbar) != 2) {printf("nbar read error\n"); exit(0);};
                if (tz > zmax) break;
                if (nsamp == 0) zinarray[NRED[nsamp]] = tz;
                nbararray[nsamp][NRED[nsamp]] = 1.0e-6*tnbar;
                NRED[nsamp]++;
            }
        }
        fclose(fp);
    }

    if (NRED[1] != NRED[0]) {
        printf("ERROR: The number of redshift bins for each sample must match\n");
        exit(0);
    }

    zarray = (double *)calloc(NRED[0], sizeof(double));
    rarray = (double *)calloc(NRED[0], sizeof(double));
    deltararray = (double *)calloc(NRED[0], sizeof(double));
    growtharray = (double *)calloc(NRED[0], sizeof(double));

    for (i=0; i<NRED[0]-1; i++) {
        zarray[i] = (zinarray[i+1]+zinarray[i])/2.0;
        rarray[i] = rz(zarray[i]);
        deltararray[i] = rz(zinarray[i+1]) - rz(zinarray[i]);
        growtharray[i] = growthz(zarray[i])/growthz(0.0);
        //printf("%12.6lf %12.6lf %12.6lf %12.6lf %12.6lf %12.6lf\n", zarray[i], rarray[i], deltararray[i], growtharray[i], nbararray[0][i], nbararray[1][i]);
    }
    zarray[NRED[0]-1] = (zmax+zinarray[NRED[0]-1])/2.0;
    rarray[NRED[0]-1] = rz(zarray[NRED[0]-1]);
    deltararray[NRED[0]-1] = rz(zmax) - rz(zinarray[NRED[0]-1]);
    growtharray[NRED[0]-1] = growthz(zarray[NRED[0]-1])/growthz(0.0);
    //printf("%12.6lf %12.6lf %12.6lf %12.6lf %12.6lf %12.6lf\n", zarray[NRED[0]-1], rarray[NRED[0]-1], deltararray[NRED[0]-1], growtharray[NRED[0]-1], nbararray[0][NRED[0]-1], nbararray[1][NRED[0]-1]);

    growth_acc    = gsl_interp_accel_alloc();
    growth_spline = gsl_spline_alloc(gsl_interp_cspline, NRED[0]);
    gsl_spline_init(growth_spline, zarray, growtharray, NRED[0]);

    free(zinarray);

    // Also create a simple redshift-distance spline
    int nbins = 400;
    double REDMIN = 0.0;
    double REDMAX = 2.0;
    double redbinwidth = (REDMAX-REDMIN)/(double)(nbins-1);
    double RMIN = rz(REDMIN);
    double RMAX = rz(REDMAX);
    double * ztemp = (double *)malloc(nbins*sizeof(double));
    double * rtemp = (double *)malloc(nbins*sizeof(double));
    for (i=0;i<nbins;i++) {
        ztemp[i] = i*redbinwidth+REDMIN;
        rtemp[i] = rz(ztemp[i]);
    }
    r_acc = gsl_interp_accel_alloc();
    r_spline = gsl_spline_alloc(gsl_interp_cspline, nbins);
    gsl_spline_init(r_spline, ztemp, rtemp, nbins);

    free(ztemp);
    free(rtemp);

    return;
}

// Routine to read in the velocity power spectrum.
void read_power() {

    FILE * fp;
    char buf[500];
    int i, j;

    pmmarray = (double**)malloc(nzin*sizeof(double*));
    pmtarray = (double**)malloc(nzin*sizeof(double*));
    pttarray = (double**)malloc(nzin*sizeof(double*));

    for (i = 0; i<nzin; i++) {

        char Pvel_file_in[500];
        sprintf(Pvel_file_in, "%s_z0p%02d.dat", Pvel_file, (int)(100.0*zin[i]));

        if(!(fp = fopen(Pvel_file_in, "r"))) {
            printf("\nERROR: Can't open power file '%s'.\n\n", Pvel_file_in);
            exit(0);
        }

        NK = 0;
        while(fgets(buf,500,fp)) {
            if(strncmp(buf,"#",1)!=0) {
                double tk, pkdelta, pkdeltavel, pkvel;
                if(sscanf(buf, "%lf %lf %lf %lf\n", &tk, &pkdelta, &pkdeltavel, &pkvel) != 4) {printf("Pvel read error\n"); exit(0);};
                NK++;
            }
        }
        fclose(fp);

        if (i == 0) {
            karray = (double *)calloc(NK, sizeof(double));
            deltakarray = (double *)calloc(NK-1, sizeof(double));
        }
        pmmarray[i] = (double *)calloc(NK, sizeof(double));
        pmtarray[i] = (double *)calloc(NK, sizeof(double));
        pttarray[i] = (double *)calloc(NK, sizeof(double));

        NK = 0;
        fp = fopen(Pvel_file_in, "r");
        while(fgets(buf,500,fp)) {
            if(strncmp(buf,"#",1)!=0) {
                double tk, pkdelta, pkdeltavel, pkvel;
                if(sscanf(buf, "%lf %lf %lf %lf\n", &tk, &pkdelta, &pkdeltavel, &pkvel) != 4) {printf("Pvel read error\n"); exit(0);};
                if (i == 0) karray[NK] = tk;
                pttarray[i][NK] = pkvel;
                pmmarray[i][NK] = pkdelta;
                pmtarray[i][NK] = pkdeltavel;
                NK++;
            }
        }
        fclose(fp);
    }

    for (i=0; i<NK-1; i++) deltakarray[i] = karray[i+1]-karray[i];

    pkkmin = karray[0];
    pkkmax = karray[NK-1];

    if (pkkmax < kmax) {
        printf("ERROR: The maximum k in the input power spectra id less than k_max\n");
        exit(0);
    }

    return;
}

// Integrand for the comoving distance
double ezinv(double x, void *p) {
  return 1.0/sqrt(Om*(1.0+x)*(1.0+x)*(1.0+x)+(1.0-Om));
}

// Calculates the comoving distance from the redshift
double rz(double red) {
  double result, error;
  gsl_function F;
  gsl_integration_workspace * w = gsl_integration_workspace_alloc(1000);
  F.function = &ezinv;
  gsl_integration_qags(&F, 0.0, red, 0, 1e-7, 1000, w, &result, &error);
  gsl_integration_workspace_free(w);
  return c*result/100.0;
}

// The integrand for the normalised growth factor
double growthfunc(double x, void *p) {
    double red = 1.0/x - 1.0;
    double Omz = Om*ezinv(red,NULL)*ezinv(red,NULL)/(x*x*x);
    double f = pow(Omz, gammaval);
    return f/x;
}

// Calculates the normalised growth factor as a function of redshift given a value of gammaval
double growthz(double red) {
  double result, error;
  gsl_function F;
  gsl_integration_workspace * w = gsl_integration_workspace_alloc(1000);
  F.function = &growthfunc;
  double a = 1.0/(1.0+red);
  gsl_integration_qags(&F, a, 1.0, 0, 1e-7, 1000, w, &result, &error);
  gsl_integration_workspace_free(w);
  return exp(-result);
}

ext: c


text:
/**
@file

Copyright John Reid 2006

*/

#include "biopsy/defs.h"

#include <gsl/gsl_rng.h>

namespace biopsy {

gsl_rng *
get_gsl_rng();

} //namespace biopsy

max_stars_repo_name: JohnReid/biopsy

max_issues_repo_name: JohnReid/biopsy

max_forks_repo_name: JohnReid/biopsy


text:
/*
MIT License

Copyright (c) 2021 Shunji Uno <shunji_uno@iscpc.jp>

Permission is hereby granted, free of charge, to any person obtaining a copy
of this software and associated documentation files (the "Software"), to deal
in the Software without restriction, including without limitation the rights
to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
copies of the Software, and to permit persons to whom the Software is
furnished to do so, subject to the following conditions:

The above copyright notice and this permission notice shall be included in all
copies or substantial portions of the Software.

THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
SOFTWARE.
*/

#pragma once
#include <stdint.h>
#include <math.h>
#include <float.h>

#ifdef MKL
#include <mkl.h>
#ifdef MKL_SBLAS
#include <mkl_spblas.h>
#endif /* MKL_SBLAS */
#else
#include <cblas.h>
#endif /* MKL */

#ifdef SXAT
#include <sblas.h>
#endif

#define DNRM2(N,X) cblas_dnrm2(N,X,1)
#define DDOT(N,X,Y) cblas_ddot(N,X,1,Y,1)

#define MATRIX_TYPE_INDEX0     (0)
#define MATRIX_TYPE_INDEX1     (1)
#define MATRIX_TYPE_ASYMMETRIC (0<<4)
#define MATRIX_TYPE_SYMMETRIC  (1<<4)
#define MATRIX_TYPE_LOWER      (1<<5)
#define MATRIX_TYPE_UPPER      (0)
#define MATRIX_TYPE_UNIT       (1<<6)
#define MATRIX_TYPE_NON_UNIT   (0)
#define MATRIX_TYPE_CSC        (0<<8)
#define MATRIX_TYPE_CSR        (1<<8)
#define MATRIX_TYPE_ELLPACK    (2<<8)
#define MATRIX_TYPE_JAD        (3<<8)
#define MATRIX_TYPE_DCSC       (8<<8)
#define MATRIX_TYPE_DCSR       (9<<8)
#define MATRIX_TYPE_MASK       (0xf<<8)

#define MATRIX_INDEX_TYPE(A)    (((A)->flags)&0xf)
#define MATRIX_IS_SYMMETRIC(A)  (((A)->flags)&MATRIX_TYPE_SYMMETRIC)
#define MATRIX_IS_LOWER(A)      (((A)->flags)&MATRIX_TYPE_LOWER)
#define MATRIX_IS_UNIT(A)       (((A)->flags)&MATRIX_TYPE_UNIT)
#define MATRIX_IS_CSC(A)        ((((A)->flags)&MATRIX_TYPE_MASK) == MATRIX_TYPE_CSC)
#define MATRIX_IS_CSR(A)        ((((A)->flags)&MATRIX_TYPE_MASK) == MATRIX_TYPE_CSR)

#ifdef __cplusplus
extern "C" {
#endif

typedef struct ellpack_info {
    int32_t nw;
    double* COEF;
    int32_t* ICOL;
} ellpack_info_t;

typedef struct jad_info {
    int32_t* perm;
    int32_t* ptr;
    int32_t* index;
    double* value;
    int32_t maxnzr;
} jad_info_t;

typedef struct Matrix {
    int NROWS;
    int NNZ;
    uint32_t flags;
    int* pointers;
    int* indice;
    double* values;
    ellpack_info_t _ellpack;
    jad_info_t _jad;
    void* info;
#ifdef MKL
    sparse_matrix_t hdl;
#endif
#ifdef SXAT
    sblas_handle_t hdl;
#endif
#ifdef SSL2
    double *w;
    int *iw;
#endif
    int optimized;
} Matrix_t;

/*
 * Generic Matrix operations
 */
void Matrix_init_generic(Matrix_t *A);
void Matrix_free_generic(Matrix_t *A);
int Matrix_optimize_generic(Matrix_t *A);
int Matrix_MV_generic(const Matrix_t *A, const double alpha, const double* x, const double beta, const double* y, double* z);

/*
 * Architecture dependent Matrix operations (prototype for override)
 */
void Matrix_init(Matrix_t *A);
void Matrix_free(Matrix_t *A);
int Matrix_optimize(Matrix_t *A);
int Matrix_MV(const Matrix_t *A, const double alpha, const double* x, const double beta, const double* y, double* z);

/*
 * Architecture independent functions
 */
void Matrix_setMatrixCSR(Matrix_t *A, const int nrows, const int nnz, const int *aptr,
    const int *aind, const double *aval, const uint32_t flags);
Matrix_t* Matrix_duplicate(const Matrix_t* A);
int Matrix_convert_index(Matrix_t* A, int base);
int Matrix_transpose(Matrix_t* A);
int Matrix_extract_symmetric(Matrix_t* A);
int Matrix_create_ellpack(Matrix_t* A);
int Matrix_create_jad(Matrix_t* A);

#ifdef __cplusplus
}
#endif

max_stars_repo_name: ISCPC/vesolver

max_issues_repo_name: ISCPC/vesolver

max_forks_repo_name: ISCPC/vesolver


text:
#ifndef CONTROLLERS_UPDATERECTIFICATIONCONTROLLER_H
#define CONTROLLERS_UPDATERECTIFICATIONCONTROLLER_H

#include "imageoperationactioncontroller.h"

#include "s3d/cv/image_operation/image_operation.h"

#include <QAction>

#include <gsl/gsl>

class ImageOperationTriggerController : public ImageOperationActionController {
Q_OBJECT

public:
  ImageOperationTriggerController(gsl::not_null<QAction*> action,
                                  gsl::not_null<s3d::image_operation::ImageOperation*> imageOperation);

  void onActionTriggered() override;
};

#endif //CONTROLLERS_UPDATERECTIFICATIONCONTROLLER_H


text:
#if !defined(compressibleFlow_h)
#define compressibleFlow_h
#include <petsc.h>
#include "fluxDifferencer.h"

typedef enum {RHO, RHOE, RHOU, RHOV, RHOW, TOTAL_COMPRESSIBLE_FLOW_COMPONENTS} CompressibleFlowComponents;
typedef enum {T, VEL, TOTAL_COMPRESSIBLE_AUX_COMPONENTS} CompressibleAuxComponents;

struct _FlowData_CompressibleFlow{
    /*Courant Friedrichs Lewy*/
    PetscReal cfl;
    /* thermal conductivity*/
    PetscReal k;
    /* dynamic viscosity*/
    PetscReal mu;

    /* store method used for flux differencer */
    FluxDifferencerFunction fluxDifferencer;

    // EOS function calls
    PetscErrorCode (*decodeStateFunction)(const PetscReal* yi, PetscInt dim, PetscReal density, PetscReal totalEnergy, const PetscReal* velocity, PetscReal* internalEnergy, PetscReal* a, PetscReal* p, void* ctx);
    void* decodeStateFunctionContext;
    PetscErrorCode (*computeTemperatureFunction)(const PetscReal* yi, PetscInt dim, PetscReal density, PetscReal totalEnergy, const PetscReal* massFlux, PetscReal* T, void* ctx);
    void* computeTemperatureContext;

    PetscBool automaticTimeStepCalculator;
} ;

typedef struct _FlowData_CompressibleFlow* FlowData_CompressibleFlow;

typedef PetscErrorCode (*FVAuxFieldUpdateFunction)(FlowData_CompressibleFlow flowData, PetscReal time, PetscInt dim, const PetscFVCellGeom *cellGeom, const PetscScalar* conservedValues, PetscScalar* auxField);

PETSC_EXTERN PetscErrorCode FVFlowUpdateAuxFieldsFV(DM dm, DM auxDM, PetscReal time, Vec locXVec, Vec locAuxField, PetscInt numberUpdateFunctions, FVAuxFieldUpdateFunction* updateFunctions, FlowData_CompressibleFlow data);
PETSC_EXTERN PetscErrorCode CompressibleFlowDiffusionSourceRHSFunctionLocal(DM dm, DM auxDM, PetscReal time, Vec locXVec, Vec locAuxVec, Vec globFVec, FlowData_CompressibleFlow flowParameters);
PETSC_EXTERN void CompressibleFlowComputeEulerFlux(PetscInt dim, PetscInt Nf, const PetscReal *qp, const PetscReal *area, const PetscReal *xL, const PetscReal *xR, PetscInt numConstants, const PetscScalar constants[], PetscReal *flux, void* ctx);
PETSC_EXTERN PetscErrorCode CompressibleFlowComputeStressTensor(PetscInt dim, PetscReal mu, const PetscReal* gradVelL, const PetscReal * gradVelR, PetscReal* tau);

#endif

max_stars_repo_name: pakserep/ablate

max_issues_repo_name: pakserep/ablate

max_forks_repo_name: pakserep/ablate


text:
#include <stdlib.h>
#include <utility>
#include <gsl/gsl_multiroots.h>
#include <gsl/gsl_sf.h>
#include <string>
#include <random>
#include <array>

using namespace std;
struct mass_eff_scalar_density_params
{
    vector<double> scalar_coeff;
    vector<double> scalar_exp;
    double number_density;
    double nucleon_mass;
    double degeneracy;
    double temperature;
    vector<double> vec_coeff;
    vector<double> vec_exp;
};
//struct for the parameters for root solving for the coupling constants
struct coeff_params
{
    double saturation_densitiy;
    double nucleon_mass;
    double binding_energy;
    double degeneracy;
    vector<double> scalar_exp;
    vector<double> vec_exp;
    unsigned int* terms;
};
//struct for the parameters for finding the critical point
struct crit_params
{
    vector<double> scalar_coeff;
    vector<double> scalar_exp;
    vector<double> vec_coeff;
    vector<double> vec_exp;
    double nucleon_mass;
    double degeneracy;
};
//struct for the parameters for finding the interaction terms for 1 CP
struct interaction_params
{
    double nucleon_mass;
    double degeneracy;
    double saturation_density;
    double binding_energy;
    double critical_temperature;
    double critical_density;
    unsigned int* terms;
};
//struct for the parameters for finding the interaction terms for 2 CPs
struct interaction_params_2crit
{
    double nucleon_mass;
    double degeneracy;
    double saturation_density;
    double binding_energy;
    double critical_temperature_lg;
    double critical_density_lg;
    double critical_temperature_qgp;
    double critical_density_qgp;
    double spinodial_l_density;
    double spinodial_r_density;
    unsigned int* terms;
};
int mass_eff_scalar_density_root(const gsl_vector * x, void *params, gsl_vector * f);
int print_state_mass_eff_scalar_density (size_t iter, gsl_multiroot_fsolver * s);
pair<double, double> get_mass_eff_scalar_density(double degeneracy, double nucleon_mass,vector<double> scalar_coeff ,vector<double>  scalar_exp,double number_density,bool print);
int T_mass_eff_mu_eff_scalar_density_root(const gsl_vector * x, void *params, gsl_vector * f);
int print_state_T_mass_eff_mu_eff_scalar_density (size_t iter, gsl_multiroot_fsolver * s);
tuple<double, double,double,bool> get_T_mass_eff_mu_eff_scalar_density(vector<double> scalar_coeff, vector<double> scalar_exp, vector<double> vec_coeff, vector<double> vec_exp, double nucleon_mass,double number_density, double temperature, double degeneracy, bool print);
int coeff_root(const gsl_vector * x, void *params, gsl_vector * f);
int print_state_coeff (size_t iter, gsl_multiroot_fsolver * s);
pair<double, double> get_coeff(double nucleon_mass, double saturation_density, double binding_energy, double degeneracy, vector<double> scalar_exp, vector<double> vec_exp, bool print, unsigned int terms[2], bool * success=NULL);
int crit_root(const gsl_vector * x, void *params, gsl_vector * f);
pair<double, double> get_crit(double nucleon_mass,double saturation_density, double binding_energy, double degeneracy, vector<double> scalar_exp, vector<double> vec_exp, bool print, double crit_T=19.09, double crit_density=540000, vector<double> scalar_coeff={0.0}, vector<double> vec_coeff={0.0});
int interaction_root(const gsl_vector * x, void *params, gsl_vector * f);
int print_state_interaction (size_t iter, gsl_multiroot_fsolver * s);
tuple<double, double,double,double, double, bool> get_interaction_4D(void * p, bool print,vector<double> init_exp={2.05,2.05},vector<double> init_coeff={10*4*M_PI/pow(550,2.0),10*4*M_PI/pow(783,2.0)});
int interaction_root_2crit(const gsl_vector * x, void *params, gsl_vector * f);
int print_state_interaction_2crit (size_t iter, gsl_multiroot_fsolver * s);
tuple<double, vector<double>, vector<double>, bool> get_interaction_2crit(void * p, bool print,vector<double> init_exp={2.05,2.05,3,3},vector<double> init_coeff={10*4*M_PI/pow(550,2.0),10*4*M_PI/pow(783,2.0),0,0});
double energy_pp_minus_mass_solv(double degeneracy, double nucleon_mass,double number_density, vector<double> scalar_coeff,vector<double> scalar_exp, vector<double> vec_coeff, vector<double> vec_exp);
double energy_pp_minus_mass_dn_solv(double degeneracy, double nucleon_mass,double number_density,vector<double> scalar_coeff,vector<double> scalar_exp, vector<double> vec_coeff, vector<double> vec_exp);
double press_solv(double degeneracy, double nucleon_mass,double number_density,vector<double> scalar_coeff,vector<double> scalar_exp, vector<double> vec_coeff, vector<double> vec_exp);
double press_dn_solv(double degeneracy, double nucleon_mass,double number_density,vector<double> scalar_coeff,vector<double> scalar_exp, vector<double> vec_coeff, vector<double> vec_exp);
double incsolv(double degeneracy, double nucleon_mass,double number_density,vector<double> scalar_coeff,vector<double> scalar_exp, vector<double> vec_coeff, vector<double> vec_exp,double delta=1);
double T_press_solv(double degeneracy, double nucleon_mass,double number_density,vector<double> scalar_coeff,vector<double> scalar_exp, vector<double> vec_coeff, vector<double> vec_exp, double temperature);
double  T_press_dmu_solv(double degeneracy, double nucleon_mass,double number_density,vector<double> scalar_coeff,vector<double> scalar_exp, vector<double> vec_coeff, vector<double> vec_exp, double temperature, double delta=1.);
double  T_press_dT_solv(double degeneracy, double nucleon_mass,double number_density,vector<double> scalar_coeff,vector<double> scalar_exp, vector<double> vec_coeff, vector<double> vec_exp, double temperature, double delta=1e-2);
double  T_press_dn_solv(double degeneracy, double nucleon_mass,double number_density,vector<double> scalar_coeff,vector<double> scalar_exp, vector<double> vec_coeff, vector<double> vec_exp, double temperature, double delta=1e3);
double  T_press_dn2_solv(double degeneracy, double nucleon_mass,double number_density,vector<double> scalar_coeff,vector<double> scalar_exp, vector<double> vec_coeff, vector<double> vec_exp,double temperature, double delta=1e3);
pair<double, double> T_press_dn12_solv(double degeneracy, double nucleon_mass,double number_density,vector<double> scalar_coeff,vector<double> scalar_exp, vector<double> vec_coeff, vector<double> vec_exp, double temperature, double delta=1e3);
double T_eps_solv(double degeneracy, double nucleon_mass,double number_density,vector<double> scalar_coeff,vector<double> scalar_exp, vector<double> vec_coeff, vector<double> vec_exp, double temperature);
void gsl_handler (const char * reason, const char * file, int line, int gsl_errno);
bool validate_interaction(double error, vector<double> coeff_guess, vector<double> exp_guess,vector<double> coeff, vector<double> exp, string filename, int timestamp, void* params);
void interaction_4D_grid(double nucleon_mass,double critical_temperature, double critical_density, double binding_energy, double saturation_density, double degeneracy, double boundaries  [4][3], unsigned int terms [2], string filename, bool print, int num_sol=0);
void interaction_4D_crit_grid(double nucleon_mass, double binding_energy, double saturation_density, double degeneracy, double boundaries_model  [4][3], double boundaries_crit [2][3], unsigned int terms [2], string filename, bool print, int num_sol);
double round_to_n_digits(double x, int n);

max_stars_repo_name: NGoetz/CWalecka

max_issues_repo_name: NGoetz/CWalecka

max_forks_repo_name: NGoetz/CWalecka


text:
/***************************************************************************
    File                 : MyParser.h
    Project              : QtiPlot
    --------------------------------------------------------------------
    Copyright            : (C) 2006 by Ion Vasilief
    Email (use @ for *)  : ion_vasilief*yahoo.fr
    Description          : Parser class based on muParser

 ***************************************************************************/

/***************************************************************************
 *                                                                         *
 *  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                                           *
 *                                                                         *
 ***************************************************************************/
#ifndef MYPARSER_H
#define MYPARSER_H

#include <muParser.h>
#include <muParserDLL.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_sf.h>

#include <qstringlist.h>

using namespace mu;

/*!\brief Mathematical parser class based on muParser.
 *
 * \section future Future Plans
 * Eliminate in favour of Script/ScriptingEnv.
 * This will allow you to use e.g. Python's global variables and functions everywhere.
 * Before this happens, a cleaner and more generic solution for accessing the current ScriptingEnv
 * should be implemented (maybe by making it a property of Project; see ApplicationWindow).
 */
class MyParser : public Parser
{
public:
	MyParser();

	static QStringList functionsList();
	static QString explainFunction(int index);

	static double bessel_J0(double x)
		{
		return gsl_sf_bessel_J0 (x);
		}

	static double bessel_J1(double x)
		{
		return gsl_sf_bessel_J1 (x);
		}

	static double bessel_Jn(double x, double n)
		{
		return gsl_sf_bessel_Jn ((int)n, x);
		}

	static double bessel_Y0(double x)
		{
		return gsl_sf_bessel_Y0 (x);
		}

	static double bessel_Y1(double x)
		{
		return gsl_sf_bessel_Y1 (x);
		}
	static double bessel_Yn(double x, double n)
		{
		return gsl_sf_bessel_Yn ((int)n, x);
		}
	static double beta(double a, double b)
		{
		return gsl_sf_beta (a, b);
		}
	static double erf(double x)
		{
		return gsl_sf_erf (x);
		}
	static double erfc(double x)
		{
		return gsl_sf_erfc (x);
		}
	static double erfz(double x)
		{
		return gsl_sf_erf_Z (x);
		}
	static double erfq(double x)
		{
		return gsl_sf_erf_Q (x);
		}
	static double gamma(double x)
		{
		return gsl_sf_gamma (x);
		}
	static double gammaln(double x)
		{
		return gsl_sf_lngamma (x);
		}
	static double hazard(double x)
		{
		return gsl_sf_hazard (x);
		}
};

#endif

max_stars_repo_name: hoehnp/SpaceDesignTool

max_issues_repo_name: hoehnp/SpaceDesignTool

max_forks_repo_name: hoehnp/SpaceDesignTool


text:
/*
 * Copyright (C) 2015  University of Oregon
 *
 * You may distribute under the terms of either the GNU General Public
 * License or the Apache License, as specified in the LICENSE file.
 *
 * For more information, see the LICENSE file.
 */
/* aslmirtime.h */
/*---------------------------------------------------------------------------*/
/*                                                                           */
/* aslmirtime.h: Global includes and defines                                 */
/*                                                                           */
/* Copyright (C) 2011 Paul Kinchesh                                          */
/*                                                                           */
/* This file is part of aslmirtime.                                          */
/*                                                                           */
/* aslmirtime is free software: you can redistribute it and/or modify        */
/* it under the terms of the GNU General Public License as published by      */
/* the Free Software Foundation, either version 3 of the License, or         */
/* (at your option) any later version.                                       */
/*                                                                           */
/* aslmirtime 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 aslmirtime. If not, see <http://www.gnu.org/licenses/>.        */
/*                                                                           */
/*---------------------------------------------------------------------------*/
/**/

/*--------------------------------------------------------*/
/*---- Wrap the header to prevent multiple inclusions ----*/
/*--------------------------------------------------------*/
#ifndef _H_aslmirtime_H
#define _H_aslmirtime_H

/*----------------------------------*/
/*---- Include standard headers ----*/
/*----------------------------------*/
#include <stdio.h>
#include <math.h>
#include <string.h>
#include <stdlib.h>
#include <sys/param.h>
#include <sys/stat.h>
#include <stddef.h>
#include <sys/time.h>

/*-----------------------------*/
/*---- Include GSL Headers ----*/
/*-----------------------------*/
#include <gsl/gsl_multimin.h>

/*-------------------------------------------*/
/*---- Make sure __FUNCTION__ is defined ----*/
/*-------------------------------------------*/
#ifndef __FUNCTION__
#define __FUNCTION__ __func__
#endif

/*------------------------------------------------------*/
/*---- Floating point comparison macros (as in SGL) ----*/
/*------------------------------------------------------*/
/* EPSILON is the largest allowable deviation due to floating point storage */
#define EPSILON 1e-9
#define FP_LT(A,B)  (((A)<(B)) && (fabs((A)-(B))>EPSILON))  /* A less than B */
#define FP_GT(A,B)  (((A)>(B)) && (fabs((A)-(B))>EPSILON))  /* A greater than B */
#define FP_EQ(A,B)  (fabs((A)-(B))<=EPSILON)                /* A equal to B */
#define FP_NEQ(A,B) (!FP_EQ(A,B))                           /* A not equal to B */
#define FP_GTE(A,B) (FP_GT(A,B) || FP_EQ(A,B))              /* A greater than or equal to B */
#define FP_LTE(A,B) (FP_LT(A,B) || FP_EQ(A,B))              /* A less than or equal to B */

/*---------------------------------------------------------------*/
/*---- So we can simply include aslmirtime.h in aslmirtime.c ----*/
/*---------------------------------------------------------------*/
#ifdef LOCAL
#define EXTERN
#else
#define EXTERN extern
#endif

/*-----------------------------------*/
/*---- Functions in aslmirtime.c ----*/
/*-----------------------------------*/
double my_f (const gsl_vector *v, void *params);
void my_df (const gsl_vector *v, void *params, gsl_vector *df);
void my_fdf (const gsl_vector *x, void *params, double *f, gsl_vector *df);
void getinput(double **pars,int argc,char *argv[]);
int nomem(char *file,const char *function,int line);
int usage_terminate();

/*----------------------------*/
/*---- End of header wrap ----*/
/*----------------------------*/
#endif

max_stars_repo_name: DanIverson/OpenVnmrJ

max_issues_repo_name: timburrow/openvnmrj-source

max_forks_repo_name: timburrow/openvnmrj-source


text:
/*
   Copyright [2021] [IBM Corporation]
   Licensed under the Apache License, Version 2.0 (the "License");
   you may not use this file except in compliance with the License.
   You may obtain a copy of the License at
       http://www.apache.org/licenses/LICENSE-2.0
   Unless required by applicable law or agreed to in writing, software
   distributed under the License is distributed on an "AS IS" BASIS,
   WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
   See the License for the specific language governing permissions and
   limitations under the License.
*/

#ifndef MCAS_HSTORE_HEAP_MC_SHIM_H
#define MCAS_HSTORE_HEAP_MC_SHIM_H

#include <mm_plugin_itf.h>

#include <ccpm/interfaces.h> /* ownership_callback, (IHeap_expandable, region_vector_t) */
#include <common/byte_span.h>
#include <common/string_view.h>
#include <gsl/span>

#include <cstddef> /* size_t */
#include <functional> /* function */

struct heap_mc_shim
	: public ccpm::IHeap_expandable
{
private:
	MM_plugin_wrapper _mm;
public:
	heap_mc_shim(common::string_view path, ccpm::persister *pe, gsl::span<common::byte_span> range, std::function<bool(const void *)> callee_owns);
	heap_mc_shim(common::string_view path, ccpm::persister *pe);
	heap_mc_shim(MM_plugin_wrapper &&pw);

	bool reconstitute(
		ccpm::region_span regions
		, ccpm::ownership_callback_t resolver
		, bool force_init
	) override;

	status_t allocate(void * & ptr
		, std::size_t bytes
		, std::size_t alignment
	) override;

	status_t free(void * & ptr
		, std::size_t bytes
	) override;

	status_t remaining(std::size_t& out_size) const override;

	ccpm::region_vector_t get_regions() const override;

	void add_regions(ccpm::region_span regions) override;

	bool includes(const void *ptr) const override;

	bool is_crash_consistent() const;
};

#endif

max_stars_repo_name: omriarad/mcas

max_issues_repo_name: omriarad/mcas

max_forks_repo_name: omriarad/mcas


text:
#pragma once

#include "Cesium3DTilesSelection/BoundingVolume.h"
#include "Cesium3DTilesSelection/Library.h"
#include "Cesium3DTilesSelection/Tile.h"
#include "Cesium3DTilesSelection/TileContext.h"
#include "Cesium3DTilesSelection/TileID.h"
#include "Cesium3DTilesSelection/TileRefine.h"
#include "Cesium3DTilesSelection/TilesetOptions.h"

#include <gsl/span>
#include <spdlog/fwd.h>

#include <cstddef>
#include <memory>

namespace Cesium3DTilesSelection {
/**
 * @brief The information that is passed to a {@link TileContentLoader} to
 * create a {@link TileContentLoadResult}.
 *
 * For many types of tile content, only the `data` field is required. The other
 * members are used for content that can generate child tiles, like external
 * tilesets or composite tiles. These members are usually initialized from
 * the corresponding members of the {@link Tile} that the content belongs to.
 */
struct CESIUM3DTILESSELECTION_API TileContentLoadInput {

  /**
   * @brief Creates a new, uninitialized instance for the given tile.
   *
   * The `data`, `contentType` and `url` will have default values,
   * and have to be initialized before this instance is passed to
   * one of the loader functions.
   *
   * @param pLogger The logger that will be used
   * @param tile The {@link Tile} that the content belongs to
   */
  TileContentLoadInput(
      const std::shared_ptr<spdlog::logger> pLogger,
      const Tile& tile);

  /**
   * @brief Creates a new instance for the given tile.
   *
   * @param pLogger The logger that will be used
   * @param data The actual data that the tile content will be created from
   * @param contentType The content type, if the data was received via a
   * network response
   * @param url The URL that the data was loaded from
   * @param tile The {@link Tile} that the content belongs to
   */
  TileContentLoadInput(
      const std::shared_ptr<spdlog::logger> pLogger,
      const gsl::span<const std::byte>& data,
      const std::string& contentType,
      const std::string& url,
      const Tile& tile);

  /**
   * @brief Creates a new instance.
   *
   * For many types of tile content, only the `data` field is required. The
   * other parameters are used for content that can generate child tiles, like
   * external tilesets or composite tiles.
   *
   * @param pLogger The logger that will be used
   * @param data The actual data that the tile content will be created from
   * @param contentType The content type, if the data was received via a
   * network response
   * @param url The URL that the data was loaded from
   * @param tileID The {@link TileID}
   * @param tileBoundingVolume The tile {@link BoundingVolume}
   * @param tileContentBoundingVolume The tile content {@link BoundingVolume}
   * @param tileRefine The {@link TileRefine} strategy
   * @param tileGeometricError The geometric error of the tile
   * @param tileTransform The tile transform
   */
  TileContentLoadInput(
      const std::shared_ptr<spdlog::logger> pLogger,
      const gsl::span<const std::byte>& data,
      const std::string& contentType,
      const std::string& url,
      const TileID& tileID,
      const BoundingVolume& tileBoundingVolume,
      const std::optional<BoundingVolume>& tileContentBoundingVolume,
      TileRefine tileRefine,
      double tileGeometricError,
      const glm::dmat4& tileTransform,
      const TilesetContentOptions& contentOptions);

  /**
   * @brief The logger that receives details of loading errors and warnings.
   */
  std::shared_ptr<spdlog::logger> pLogger;

  /**
   * @brief The raw input data.
   *
   * The {@link TileContentFactory} will try to determine the type of the
   * data using the first four bytes (i.e. the "magic header"). If this
   * does not succeed, it will try to determine the type based on the
   * `contentType` field.
   */
  gsl::span<const std::byte> data;

  /**
   * @brief The content type.
   *
   * If the data was obtained via a HTTP response, then this will be
   * the `Content-Type` of that response. The {@link TileContentFactory}
   * will try to interpret the data based on this content type.
   *
   * If the data was not directly obtained from an HTTP response, then
   * this may be the empty string.
   */
  std::string contentType;

  /**
   * @brief The source URL.
   */
  std::string url;

  /**
   * @brief The {@link TileID}.
   */
  TileID tileID;

  /**
   * @brief The tile {@link BoundingVolume}.
   */
  BoundingVolume tileBoundingVolume;

  /**
   * @brief Tile content {@link BoundingVolume}.
   */
  std::optional<BoundingVolume> tileContentBoundingVolume;

  /**
   * @brief The {@link TileRefine}.
   */
  TileRefine tileRefine;

  /**
   * @brief The geometric error.
   */
  double tileGeometricError;

  /**
   * @brief The tile transform
   */
  glm::dmat4 tileTransform;

  /**
   * @brief Options for parsing content and creating Gltf models.
   */
  TilesetContentOptions contentOptions;
};
} // namespace Cesium3DTilesSelection

max_stars_repo_name: zrkcode/cesium-native

max_issues_repo_name: zrkcode/cesium-native

max_forks_repo_name: zrkcode/cesium-native


text:
/*
 * Copyright 2020 Makani Technologies LLC
 *
 * Licensed under the Apache License, Version 2.0 (the "License");
 * you may not use this file except in compliance with the License.
 * You may obtain a copy of the License at
 *
 *      http://www.apache.org/licenses/LICENSE-2.0
 *
 * Unless required by applicable law or agreed to in writing, software
 * distributed under the License is distributed on an "AS IS" BASIS,
 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 * See the License for the specific language governing permissions and
 * limitations under the License.
 */

// The GSL extra library provides linear algebra functions using a
// similar interface to the rest of the GSL library.

#ifndef COMMON_C_MATH_GSL_LINALG_EXTRA_H_
#define COMMON_C_MATH_GSL_LINALG_EXTRA_H_

#include <stdint.h>

#include <gsl/gsl_blas_types.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_vector.h>

#ifdef __cplusplus
extern "C" {
#endif

int32_t GslTriangularSolve(CBLAS_SIDE_t side, CBLAS_UPLO_t uplo,
                           CBLAS_TRANSPOSE_t transpose, const gsl_matrix *T,
                           const gsl_matrix *B, gsl_matrix *X);
void GslTrapezoidalToTriangular(const gsl_matrix *R, gsl_matrix *T,
                                gsl_vector *tau);
void GslTrapezoidalToTriangularZTMat(const gsl_matrix *T, const gsl_vector *tau,
                                     const gsl_matrix *X, gsl_matrix *Zt_X);
int32_t GslMatrixDivide(CBLAS_SIDE_t side, const gsl_matrix *A,
                        const gsl_matrix *B, gsl_matrix *X);

#ifdef __cplusplus
}  // extern "C"
#endif

#endif  // COMMON_C_MATH_GSL_LINALG_EXTRA_H_


text:
#ifndef __GEMM_DRIVER_H
#define __GEMM_DRIVER_H

#include "util.h"

#include <stddef.h>
#include <stdlib.h>
#include <string.h>
#include <tuple>
#include <algorithm>
#include <assert.h>
#include <iostream>
#include <sstream>
#include <vector>
// openblas
#include <cblas.h>

#ifndef CEIL
#define CEIL(value, divider)  (   ((value)-1)/(divider)+1  )
#endif

#ifndef CEIL_WRAP
#define CEIL_WRAP(value, divider)  ( CEIL((value), (divider)) * (divider)  )
#endif

typedef enum {
    LAYOUT_ROW_MAJOR = 0,
    LAYOUT_COL_MAJOR
}layout_t;

typedef enum {
    TRANS_NO_TRANS = 0,
    TRANS_TRANS,
    TRANS_CONJ_TRANS,
    TRANS_CONJ_NO_TRANS,
}trans_t;

typedef enum {
    IDENT_A_MATRIX = 0,
    IDENT_B_MATRIX
}identifier_t;

// cblas helper function
static inline CBLAS_ORDER to_blas_layout(layout_t layout){
    if(layout == LAYOUT_ROW_MAJOR)
        return CblasRowMajor;
    //if(layout == LAYOUT_COL_MAJOR)
    //TODO: validation
    return CblasColMajor;
}
// cblas helper function
static inline CBLAS_TRANSPOSE to_blas_transpose(trans_t trans){
    switch(trans){
        case  TRANS_NO_TRANS:
            return CblasNoTrans;
        case  TRANS_TRANS:
            return CblasTrans;
        case TRANS_CONJ_TRANS:
            return CblasConjTrans;
        case TRANS_CONJ_NO_TRANS:
            return CblasConjNoTrans;
        default:
            return CblasConjNoTrans;
    }
}
static inline const char * to_layout_str(layout_t layout){
    if(layout == LAYOUT_ROW_MAJOR)
        return "CblasRowMajor";
    if(layout == LAYOUT_COL_MAJOR)
        return "CblasColMajor";
    return "n/a major";
}

static inline const char * to_trans_str(trans_t trans){
    if(trans == TRANS_NO_TRANS)
        return "CblasNoTrans";
    if(trans == TRANS_TRANS)
        return "CblasTrans";
    if(trans == TRANS_CONJ_TRANS)
        return "CblasConjTrans";
    if(trans == TRANS_CONJ_NO_TRANS)
        return "CblasConjNoTrans";
    return "n/a trans";
}

class gemm_context_t {
public:
// matrix descriptors
    layout_t    layout;
    trans_t     trans_a;
    trans_t     trans_b;
    size_t      m;
    size_t      n;
    size_t      k;

    size_t      lda;
    size_t      ldb;
    size_t      ldc;
    double      alpha;
    double      beta;

//
    size_t      alignment;  // used for alloc A/B/C, indeed not so useful

// blocking parameters
    size_t      mc;
    size_t      nc;
    size_t      kc;
    size_t      mr;
    size_t      nr;

// hw parameters
    //size_t      cpu_id;
    size_t      l1_size;
    size_t      l2_size;
    size_t      l3_size;
    size_t      tlb_entry_l1d;
    size_t      cacheline_size;
    size_t      page_size;
    std::vector<int>    cpu_list;

    double      frequency;  // MHz

    bool        cur_use_tuned {false};  // use tuned param from db or default, only a flag

    void serialize_layout_trans(std::ostream & os){
        std::string al, bl, cl;
        if(layout == LAYOUT_ROW_MAJOR){
            cl = "cr";
            if(trans_a == TRANS_NO_TRANS || trans_a == TRANS_CONJ_NO_TRANS){
                al = "ar";
            }else{
                al = "ac";
            }
            if(trans_b == TRANS_NO_TRANS || trans_b == TRANS_CONJ_NO_TRANS){
                bl = "br";
            }else{
                bl = "bc";
            }
        }else{
            cl = "cc";
            if(trans_a == TRANS_NO_TRANS || trans_a == TRANS_CONJ_NO_TRANS){
                al = "ac";
            }else{
                al = "ar";
            }
            if(trans_b == TRANS_NO_TRANS || trans_b == TRANS_CONJ_NO_TRANS){
                bl = "bc";
            }else{
                bl = "br";
            }
        }
        os<<al<<"-"<<bl<<"-"<<cl;
    }

    void deserialize_layout_trans(const std::istream & is){
        std::string al, bl, cl;
        std::stringstream tmp;
        tmp << is.rdbuf();
        std::string s( tmp.str() );

        al = s.substr(0, 2);
        bl = s.substr(3, 2);
        cl = s.substr(6, 2);

        if(cl == "cr"){
            layout = LAYOUT_ROW_MAJOR;
            if(al == "ar"){
                trans_a = TRANS_NO_TRANS;
            }else{
                trans_a = TRANS_TRANS;
            }
            if(bl == "br"){
                trans_b = TRANS_NO_TRANS;
            }else{
                trans_b = TRANS_TRANS;
            }
        }else if(cl == "cc"){
            layout = LAYOUT_COL_MAJOR;
            if(al == "ar"){
                trans_a = TRANS_TRANS;
            }else{
                trans_a = TRANS_NO_TRANS;
            }
            if(bl == "br"){
                trans_b = TRANS_TRANS;
            }else{
                trans_b = TRANS_NO_TRANS;
            }
        }
    }

    // TODO: fix more param
    void serialize(std::ostream & os){
        serialize_layout_trans(os);
        os<<"-"<<m<<"-"<<n<<"-"<<k;
    }
    void serialize(std::string & str){
        std::ostringstream oss;
        serialize(oss);
        str = oss.str();
    }
    void deserialize(std::istream & is){
        char _d;

        deserialize_layout_trans(is);
        is>>_d>>m>>_d>>n>>_d>>k;
    }
    void deserialize(std::string & str){
        std::istringstream iss;
        iss.str(str);
        deserialize(iss);
    }
};

// https://software.intel.com/en-us/mkl-developer-reference-c-cblas-gemm

class matrix_elem_t {
public:
    size_t operator() (size_t row, size_t col, size_t ldim, layout_t layout, trans_t trans){
        if(layout == LAYOUT_ROW_MAJOR){
            if(trans == TRANS_NO_TRANS || trans == TRANS_CONJ_NO_TRANS){
                assert(ldim>=col);
                return row * ldim;
            }
            else{
                assert(ldim>=row);
                return col * ldim;
            }
        }else{
            if(trans == TRANS_NO_TRANS || trans == TRANS_CONJ_NO_TRANS){
                assert(ldim>=row);
                return col * ldim;
            }
            else{
                assert(ldim>=col);
                return row * ldim;
            }
        }
    }
};

template<typename T>
class matrix_t{
public:
    matrix_t(size_t row_, size_t col_, size_t ldim_,
        layout_t layout_, trans_t trans_, size_t alignment_)
    {
        this->row = row_;
        this->col = col_;
        this->ldim = ldim_;
        // TODO ldim should be multiple of alignment!
        this->layout = layout_;
        this->trans = trans_;
        this->alignment = alignment_;

        size_t elements = matrix_elem_t()(row_, col_, ldim_, layout_, trans_);
        this->data = (T *) __aligned_malloc(sizeof(T)*elements, alignment_);
        rand_vector(this->data, elements);
    }
    ~matrix_t(){
        if(data)
            __aligned_free(this->data);
    }

    size_t dtype_size() const {
        return sizeof(T);
    }
    void __copy(const matrix_t<T> & rhs){
        this->row = rhs.row;
        this->col = rhs.col;
        this->ldim = rhs.ldim;

        this->layout = rhs.layout;
        this->trans = rhs.trans;
        this->alignment = rhs.alignment;
        size_t elements = matrix_elem_t()(rhs.row, rhs.col, rhs.ldim, rhs.layout, rhs.trans);
        this->data = (T *) __aligned_malloc(sizeof(T)*elements, alignment);
        memcpy(this->data, rhs.data, sizeof(T)*elements);
    }

    matrix_t(const matrix_t<T> & rhs){
        __copy(rhs);
    }
    matrix_t& operator =(const matrix_t<T> & rhs){
        __copy(rhs);
        return *this;
    }

    T           *data {nullptr};
    size_t      row;
    size_t      col;
    size_t      ldim;   // leading dimension

    layout_t    layout;
    trans_t     trans;
    size_t      alignment;
};

//typedef matrix_t<float> matrix_fp32_t;
//typedef matrix_t<double> matrix_fp64_t;

#endif