text:
/*=========================================================================

  Copyright (c) Yue Sheng, University of Manchester and David Worth, STFC
  All rights reserved.

     This software is distributed WITHOUT ANY WARRANTY; without even
     the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
     PURPOSE.  See the above copyright notice for more information.

=========================================================================*/

/** \file vvQuan3D.c
\brief This is a VolView plugin to calculate several characteristics from a
labelled image.

The following characteristics are calculated:
\li Volume by voxel counts
\li Equivalent sphere diameter by voxel counts
\li Bounding box diagonal
\li Principal Component Analysis
\li Ellipsoid fitting by PCA
\li Equivalent circle diameter by PCA
\li Isosurface by marching cube
\li Surface area
\li Surface volume
\li Equivalent sphere diameter from surface volume
\li Sphercity
\li Normalised surface area to volume ratio (Radius*Sa/Vol)
*/

#include "vvHelper.h"
#include "vtkVVPluginAPI.h"
#include <dlfcn.h>
#include <limits.h>
#include <math.h>
#include <string.h>
#include <stdlib.h>
#include <stdio.h>
#include <gsl/gsl_eigen.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_vector.h>

/**
Struct used when calculating the number of voxels with each voxel value
*/
typedef struct {
  /** The voxel value */
  int value;
  /** The number of voxels with this value */
  int count;
} vvQuan3DVoxelCount;

/**
Comparison function for sorting the image data.
\param a Integer value
\param b Integer value
\return 1 if \c a>b; -1 if \c a<b; 0 if \c a==b.
*/
int cmpfun(void *a, void *b)
{
  return ( *(int*)a - *(int*)b );
}

/**
\brief Calculate sphere diameter from the sphere volume.

Use the formula \f$ 2*(3v/(4\pi))^{\frac{1}{3}} \f$ where \f$ v \f$ is the
sphere volume.

\param volume The sphere volume
\return Diameter of the sphere
*/
double vol2dia(double volume)
{
  return 2.0 * pow((3.0*volume)/(4.0*3.14159265),(1.0/3.0));
}

/**
\brief Find the (i,j,k) coordinates for each voxel id and also calculate
maximum, minimum and mean coordinates.

The positions pointer must be allocated by the calling function as follows
\code
positions = malloc(3*sizeof(int)*voxelCount);
\endcode

\param voxelSubArray Array of voxel ids
\param voxelCount Number of voxels in voxelSubArray <b>NB This is not the
length of the array</b>
\param dim The dimsnesions of the image
\param positions Pointer to matrix which will hold the (i,j,k) coordinates
\param min  The minimum coordinate in each dimension
\param max  The maximum coordinate in each dimension
\param mean The mean of the coordinates in each dimension
*/
void getCoords(int *voxelSubArray, int voxelCount, int dim[],
               gsl_matrix *positions, int min[], int max[], int mean[])
{
  int i;
  int voxelId;
  int icoord, jcoord, kcoord, rem;
  min[0] = min[1] = min[2] = INT_MAX;
  max[0] = max[1] = max[2] = INT_MIN;
  mean[0] = mean[1] = mean[2] = 0;

  /* Loop over the voxels and calulate i,j,k coords for each one */
  for (i = 0; i < voxelCount; i++) {
    voxelId = voxelSubArray[i];
    kcoord = floor(voxelId / (dim[0]*dim[1]));
    rem = voxelId % (dim[0]*dim[1]);
    jcoord = floor(rem / dim[0]);
    icoord = rem % dim[0];
    //printf("voxel %d (%d,%d,%d)\n",voxelId,icoord,jcoord,kcoord);
    gsl_matrix_set(positions,i,0,icoord);
    gsl_matrix_set(positions,i,1,jcoord);
    gsl_matrix_set(positions,i,2,kcoord);
    if (icoord > max[0]) max[0] = icoord;
    if (jcoord > max[1]) max[1] = jcoord;
    if (kcoord > max[2]) max[2] = kcoord;
    if (icoord < min[0]) min[0] = icoord;
    if (jcoord < min[1]) min[1] = jcoord;
    if (kcoord < min[2]) min[2] = kcoord;
    mean[0] += icoord;
    mean[1] += jcoord;
    mean[2] += kcoord;
  }
  mean[0] /= voxelCount;
  mean[1] /= voxelCount;
  mean[2] /= voxelCount;
}

/**
\brief Calculate the data from the image

\author David Worth, STFC
\date May 2012
\todo This only works for unsigned short image type. Make it do other integer
datatypes.

\param inf Pointer to object that holds data for this plugin. The voxel size
can be obtained from it.
\param pds Pointer to object that carries information on data set to be
processed. Includes actual buffer of voxel data, number of voxels along each
dimension, voxel spacing and voxel type.
*/
static int ProcessData(void *inf, vtkVVProcessDataStruct *pds)
{
  /* Information about the plugin and the image */
  vtkVVPluginInfo *info = (vtkVVPluginInfo *)inf;
  /* Number of voxels in each direction of image */
  int *dim = info->InputVolumeDimensions;
  /* Size of datatype in the image */
  unsigned int rsize = info->InputVolumeScalarSize;
  /* Voxel size - can be set by user in VolView */
  double voxelSize = atof(info->GetGUIProperty(info, 0, VVP_GUI_VALUE));
  /* Minimum and maximum raw values from image. */
  double imgMin = info->InputVolumeScalarRange[0],
         imgMax = info->InputVolumeScalarRange[1];
  /* Array of voxel counts for each value. Array index is the value. */
  int maxVoxelVals = (int)ceil(imgMax)-(int)floor(imgMin)+1;
  int *voxelCounts = malloc(sizeof(int)*maxVoxelVals);
  /* Array of arrays - each sub-array holding the ids of voxels with that value */
  int* *voxelArray = malloc(sizeof(int*)*maxVoxelVals);
  /* Pointer to sub array used when creating voxelArray */
  int *voxelSubArray = NULL;
  /* Index array used to write into voxelSubArray */
  int *indexVoxelSubArray = malloc(sizeof(int)*maxVoxelVals);
  /* Number of the non-background voxels */
  int sumVoxels = 0;
  /* Number of different voxel values */
  int totalParticles = 0;
  /* Array of (i,j,k) coordinates for voxels with a particular value */
  gsl_matrix *positions = NULL;
  /* Stats for voxel data for each value */
  int max[3], min[3], mean[3], boundingBoxSize[3];
  double boundingBoxDiam;
  /* Multiplication factors based on physical voxel size */
  double areaFactor = pow(voxelSize,2), volumeFactor = pow(voxelSize,3);
  /* Output data */
  int voxelLabel, voxelCount;
  double voxelCountVolume, voxelCountDia;
  /* Data used when covariance matrix is calculated */
  double* *matrix = NULL;
  double *row = NULL;
  double* *covMatrix = NULL;
  /* Data for eigenanalysis of covariance matrix */
  gsl_matrix *gslCovMatrix = NULL;
  gsl_vector *eval = NULL;
  gsl_matrix *evecs = NULL;
  gsl_eigen_symmv_workspace *workspace = NULL;
  gsl_vector vector;
  int result;

  void *gslcblas_handle, *gsl_handle;

  /* Buffer for output to user */
  char buffer[1000];
  int i, j = 0;

  /* The raw image data */
  unsigned short *ptr = (unsigned short *)pds->inData;

  /* Check the data type */
  switch (info->InputVolumeScalarType) {
    case VTK_CHAR:
    case VTK_UNSIGNED_CHAR:
    case VTK_FLOAT:
    case VTK_DOUBLE: {
      sprintf(buffer, "Image datatype must be integer or boolean. This \
      algorithm cannot continue");
      info->SetProperty(info, VVP_ERROR, buffer);
      return -1;
      break;
    }
  }

  printf("Image has type %d size %d\n", info->InputVolumeScalarType,rsize);
  printf("The dimensions of the image are: %d, %d, %d \n",
          dim[0], dim[1], dim[2]);
  sprintf(buffer, "The dimensions of the image are: %d, %d, %d \n",
          dim[0], dim[1], dim[2]);
  printf("There are %d components\n",info->InputVolumeNumberOfComponents);
  printf("The minimum value in the image is %d\n", (int)imgMin);
  printf("and the maximum value is %d\n", (int)imgMax);

  /* Sort the image values */
  /*qsort(ptr, (size_t)(dim[0]*dim[1]*dim[2]), (size_t)info->InputVolumeScalarSize,
        (__compar_fn_t)cmpfun);*/

  /*
   * Create array of voxel counts for each voxel value except the min value
   * which we take to be the background
   */
  for (i = 0; i < maxVoxelVals; i++) voxelCounts[i]=0;
  for (i = 0; i < dim[0]*dim[1]*dim[2]; i++) {
    if (ptr[i] != imgMin) { /* Not the minimum value */
      if (ptr[i] > imgMax) printf("%d %d\n",i,ptr[i]); /* Something wrong */
      /* Not seen this value before so it's a new particle */
      if (voxelCounts[ptr[i]] == 0) totalParticles++;
      voxelCounts[ptr[i]]++;
      sumVoxels++;
    }
  }

  /* Create the array of arrays with voxel ids for each voxel value */
  for (i = 0; i < maxVoxelVals; i++) {
    if (voxelCounts[i] > 0) {
      /*printf("Non zero voxel count for i = %d - count is %d\n",i,voxelCounts[i]);*/
      voxelArray[i] = malloc(voxelCounts[i]*sizeof(int));
    }
    else voxelArray[i] = NULL;
    indexVoxelSubArray[i] = 0;
  }
  /* Now populate that array with values */
  for (i = 0; i < dim[0]*dim[1]*dim[2]; i++) {
    if (ptr[i] != imgMin) {
      //printf("add entry to voxel array %d - entry id %d\n",ptr[i],i);
      voxelSubArray = voxelArray[ptr[i]];
      voxelSubArray[ indexVoxelSubArray[ptr[i]] ] = i;
      indexVoxelSubArray[ptr[i]]++;
    }
  }

  printf("Porosity is %f\n",(double)sumVoxels/(double)(dim[0]*dim[1]*dim[2]));
  sprintf(buffer+strlen(buffer), "Porosity is %f\n",
          (double)sumVoxels/(double)(dim[0]*dim[1]*dim[2]));
  printf("Number of particles is %d\n",totalParticles);

  /* Load the GSL BLAS and GSL itself. RTLD_GLOBAL means we can use the
   * function names from GSL rather than having to get function pointers with
   * dlsym. */
  gslcblas_handle = dlopen("/usr/lib/libgslcblas.so", RTLD_LAZY|RTLD_GLOBAL);
  if (!gslcblas_handle) {
    printf("Error opening libgslcblas.so - %s\n",dlerror());
  }
  gsl_handle = dlopen("/usr/lib/libgsl.so", RTLD_LAZY|RTLD_GLOBAL);
  if (!gsl_handle) {
    printf("Error opening libgsl.so - %s\n",dlerror());
  }

  /* Allocate the GSL matrices and vectors for eigenanalysis */
  gslCovMatrix = gsl_matrix_alloc(3,3);
  eval = gsl_vector_alloc(3);
  evecs = gsl_matrix_alloc(3,3);
  workspace = gsl_eigen_symmv_alloc(3);

  /* Loop over the particles */
  for (i = 0; i < 1/*totalParticles*/; i++) {
    /* Find the next particle with a non-zero voxel count */
    while ((j < maxVoxelVals) && (voxelCounts[j] == 0)) j++;
    /* Do some simple stuff */
    voxelLabel = j;
    voxelCount = voxelCounts[j];
    voxelCountVolume = voxelCount * volumeFactor;
    voxelCountDia = voxelSize*vol2dia(voxelCount);
    printf("%d, %d, %f, %f",voxelLabel,voxelCount,voxelCountVolume,
      voxelCountDia);

    /* Find the (i,j,k) coordinates in the image of the voxels */
    positions = gsl_matrix_alloc(voxelCount,3);
    voxelSubArray = voxelArray[voxelLabel];
    getCoords(voxelSubArray, voxelCount, dim, positions, min, max, mean);
    /*printf(" (%d, %d, %d)",min[0],min[1],min[2]);
    printf(" (%d, %d, %d)",max[0],max[1],max[2]);
    printf(" (%d, %d, %d)",mean[0],mean[1],mean[2]);*/
    boundingBoxSize[0] = max[0]-min[0]+1;
    boundingBoxSize[1] = max[1]-min[1]+1;
    boundingBoxSize[2] = max[2]-min[2]+1;
    boundingBoxDiam = sqrt(boundingBoxSize[0]*boundingBoxSize[0] +
                           boundingBoxSize[1]*boundingBoxSize[1] +
                           boundingBoxSize[2]*boundingBoxSize[2])*voxelSize;
    printf(" (%d, %d, %d)",boundingBoxSize[0],boundingBoxSize[1],
                           boundingBoxSize[2]);
    printf(" %f", boundingBoxDiam);

    /* Calculate positions-mean */
    matrix = malloc(voxelCount*sizeof(double*));
    for (j = 0; j < voxelCount; j++) {
      row = malloc(3*sizeof(double));
      row[0] = gsl_matrix_get(positions,j,0)-mean[0];
      row[1] = gsl_matrix_get(positions,j,1)-mean[1];
      row[2] = gsl_matrix_get(positions,j,2)-mean[2];
      matrix[j] = row;
    }
    covMatrix = CovMat(matrix, voxelCount, 3);

    for (j = 0; j < voxelCount; j++) {
      free(matrix[j]);
    }
    free(matrix);

    /*printf("\nCovariance matrix\n");
    for (j = 0; j < 3; j++) {
      printf("%f %f %f\n", covMatrix[j][0], covMatrix[j][1], covMatrix[j][2]);
    }*/

    /* To do the eigenanalysis use GNU Scientific library so create matrix
     * in that format */
    for (j = 0; j < 3; j++) {
      gsl_matrix_set(gslCovMatrix,0,j,covMatrix[0][j]);
      gsl_matrix_set(gslCovMatrix,1,j,covMatrix[1][j]);
      gsl_matrix_set(gslCovMatrix,2,j,covMatrix[2][j]);
    }

    /* Calculate eigenvalues and vectors */
    result = gsl_eigen_symmv(gslCovMatrix, eval, evecs, workspace);

    printf("\nEigenvalues\n");
    gsl_vector_fprintf(stdout, eval, "%g");
    printf("\nEigenvectors\n");
    vector = gsl_matrix_const_column(evecs,0).vector;
    gsl_vector_fprintf(stdout, &vector, "%g");
    printf("\n");
    vector = gsl_matrix_const_column(evecs,1).vector;
    gsl_vector_fprintf(stdout, &vector, "%g");
    printf("\n");
    vector = gsl_matrix_const_column(evecs,2).vector;
    gsl_vector_fprintf(stdout, &vector, "%g");

    for (j = 0; j < 3; j++) {
      free(covMatrix[j]);
    }
    free(covMatrix);

    gsl_matrix_free(positions);

    printf("\n");
    /* Move on to next particle */
    j++;
  }

  dlclose(gsl_handle);
  dlclose(gslcblas_handle);

  info->SetProperty(info, VVP_REPORT_TEXT, buffer);

  free(voxelCounts);
  for (i = 0; i < maxVoxelVals; i++) {
    if (voxelArray[i] != NULL) free(voxelArray[i]);
  }

  gsl_matrix_free(gslCovMatrix);
  gsl_vector_free(eval);
  gsl_matrix_free(evecs);
  gsl_eigen_symmv_free(workspace);

  free(voxelArray);
  free(indexVoxelSubArray);
  return 0;
}

/**
\brief Update the VolView GUI to display user parameters.

Sets one GUI parameter - the physical size of the voxels in the image.

\param inf Pointer to object that should be modified to set up GUI elements for
this plugin. It also contains details of the input and output images.
*/
static int UpdateGUI(void *inf)
{
  int i;
  vtkVVPluginInfo *info = (vtkVVPluginInfo *)inf;

  info->SetGUIProperty(info, 0, VVP_GUI_LABEL, "Voxel Size");
  info->SetGUIProperty(info, 0, VVP_GUI_TYPE, VVP_GUI_SCALE);
  info->SetGUIProperty(info, 0, VVP_GUI_DEFAULT , "1");
  info->SetGUIProperty(info, 0, VVP_GUI_HELP,
               "The physical size of a voxel in the image");

  /* what range should we show for possible output values */
  vvPluginSetGUIScaleRange(0);

  info->OutputVolumeScalarType = info->InputVolumeScalarType;
  info->OutputVolumeNumberOfComponents = info->InputVolumeNumberOfComponents;
  for (i = 0; i < 3; i++)
    {
    info->OutputVolumeDimensions[i] = info->InputVolumeDimensions[i];
    info->OutputVolumeSpacing[i] = info->InputVolumeSpacing[i];
    info->OutputVolumeOrigin[i] = info->InputVolumeOrigin[i];
    }

  return 1;
}

/**
\brief Initialise this plugin.

Sets the name and group for this plugin in the plugin list shown to the VolView
user and gives some documentation. Also defines properties so VolView can judge
the memory requirements and potential for undoing this plugin.

\param info Pointer to object that should be modified to give details about this
plugin.
*/
void VV_PLUGIN_EXPORT vvQuan3DCInit(vtkVVPluginInfo *info)
{
  /* always check the version */
  vvPluginVersionCheck();

  /* setup information that never changes */
  info->ProcessData = ProcessData;
  info->UpdateGUI = UpdateGUI;
  info->SetProperty(info, VVP_NAME, "Quantify 3D in C");
  info->SetProperty(info, VVP_GROUP, "Quantification");
  info->SetProperty(info, VVP_TERSE_DOCUMENTATION,
     "Quantify several characteristics from a labelled image");
  info->SetProperty(info, VVP_FULL_DOCUMENTATION,
    "Quantify the following characteristics from a labelled image: Volume by \
voxel counts, Equivalent sphere diameter by voxel counts, Bounding box \
diagonal, Principal Component Analysis, Ellipsoid fitting by PCA, Equivalent \
circle diameter by PCA, Isosurface by marching cube, Surface area, Surface \
volume, Equivalent sphere diameter from surface volume, Sphercity, Normalised \
surface area to volume ratio");

  info->SetProperty(info, VVP_SUPPORTS_IN_PLACE_PROCESSING, "1");
  info->SetProperty(info, VVP_SUPPORTS_PROCESSING_PIECES,   "1");
  info->SetProperty(info, VVP_REQUIRED_Z_OVERLAP,           "0");
  info->SetProperty(info, VVP_NUMBER_OF_GUI_ITEMS,          "1");
  info->SetProperty(info, VVP_REQUIRES_SERIES_INPUT,        "0");
  info->SetProperty(info, VVP_SUPPORTS_PROCESSING_SERIES_BY_VOLUMES, "0");
  info->SetProperty(info, VVP_PRODUCES_OUTPUT_SERIES, "0");
  info->SetProperty(info, VVP_PRODUCES_PLOTTING_OUTPUT, "0");
}

max_stars_repo_name: vais-ral/CCPi-Quantification

max_issues_repo_name: vais-ral/CCPi-Quantification

max_forks_repo_name: vais-ral/CCPi-Quantification


text:
#include <erl_nif.h>
#include <stdio.h>
#include <string.h>
#include <math.h>
#include <cblas.h>

/*
--------------------------------------------------------------------|
            ------------------------------                          |
            |          LAPACKE           |                          |
            ------------------------------                          |
--------------------------------------------------------------------|
*/

void dgetrf_(int* M, int *N, double* A, int* lda, int* IPIV, int* INFO);
void dgetri_(int* N, double* A, int* lda, int* IPIV, double* WORK, int* lwork, int* INFO);

/*
--------------------------------------------------------------------|
            ------------------------------                          |
            |           INIT FC          |                          |
            ------------------------------                          |
--------------------------------------------------------------------|
*/

ERL_NIF_TERM atom_nok;
ERL_NIF_TERM atom_true;
ERL_NIF_TERM atom_false;
ERL_NIF_TERM atom_matrix;

ErlNifResourceType *MULT_YIELDING_ARG = NULL;

int load(ErlNifEnv* env, void** priv_data, ERL_NIF_TERM load_info){
    atom_nok = enif_make_atom(env, "nok\0");
    atom_true = enif_make_atom(env, "true\0");
    atom_false = enif_make_atom(env, "false\0");
    atom_matrix = enif_make_atom(env, "matrix\0");
    return 0;
}

int upgrade(ErlNifEnv* caller_env, void** priv_data, void** old_priv_data, ERL_NIF_TERM load_info){
    return 0;
}

//Gives easier access to an ErlangBinary containing a matrix.
typedef struct{
    int n_rows;
    int n_cols;

    //Content of the matrix, in row major format.
    double* content;

    //Erlang binary containing the matrix.
    ErlNifBinary bin;
} Matrix;

//Access asked coordinates of matrix
double* matrix_at(int col, int row, Matrix m){
    return &m.content[row*m.n_cols + col];
}

//Allocates memory space of for matrix of dimensions n_rows, n_cols.
//The matrix_content can be modified, until a call to array_to_erl.
//Matrix content is stored in row major format.
Matrix matrix_alloc(int n_rows, int n_cols){
    ErlNifBinary bin;

    enif_alloc_binary(sizeof(double)*n_rows*n_cols, &bin);

    Matrix matrix;
    matrix.n_cols = n_cols;
    matrix.n_rows = n_rows;
    matrix.bin = bin;
    matrix.content = (double*) bin.data;

    return matrix;
}

//Creates a duplicate of a matrix.
//This duplicate can be modified until uploaded.
Matrix matrix_dup(Matrix m){
    Matrix d = matrix_alloc(m.n_rows, m.n_cols);
    memcpy(d.content, m.content, d.bin.size);
    return d;
}

//Free an allocated matrix that was not sent back to Erlang.
void matrix_free(Matrix m){
    enif_release_binary(&m.bin);
}

//Constructs a matrix erlang term.
//No modifications can be made afterwards to the matrix.
ERL_NIF_TERM matrix_to_erl(ErlNifEnv* env, Matrix m){
    ERL_NIF_TERM term = enif_make_binary(env, &m.bin);
    return enif_make_tuple4(env, atom_matrix, enif_make_int(env,m.n_rows), enif_make_int(env,m.n_cols), term);
}

int enif_is_matrix(ErlNifEnv* env, ERL_NIF_TERM term){
    int arity;
    const ERL_NIF_TERM* content;

    if(!enif_is_tuple(env, term))
        return 0;

    enif_get_tuple(env, term, &arity, &content);
    if(arity != 4)
        return 0;

    if(content[0] != atom_matrix
        || !enif_is_number(env, content[1])
        || !enif_is_number(env, content[2])
        || !enif_is_binary(env, content[3]))

            return 0;
    return 1;
}

//Reads an erlang term as a matrix.
//As such, no modifications can be made to the red matrix.
//Returns true if it was possible to read a matrix, false otherwise
int enif_get_matrix(ErlNifEnv* env, ERL_NIF_TERM term, Matrix *dest){

    int arity;
    const ERL_NIF_TERM* content;

    if(!enif_is_tuple(env, term))
        return 0;

    enif_get_tuple(env, term, &arity, &content);
    if(arity != 4)
        return 0;

    if(content[0] != atom_matrix
        || !enif_get_int(env, content[1], &dest->n_rows)
        || !enif_get_int(env, content[2], &dest->n_cols)
        || !enif_inspect_binary(env, content[3], &dest->bin))
    {
        return 0;
    }

    dest->content = (double*) (dest->bin.data);
    return 1;
}

//Used to translate at once a number of ERL_NIF_TERM.
//Data types are inferred via content of format string:
//  n: number (int or double) translated to double.
//  m: matrix
//  i: int
int enif_get(ErlNifEnv* env, const ERL_NIF_TERM* erl_terms, const char* format, ...){
    va_list valist;
    va_start(valist, format);
    int valid = 1;

    while(valid && *format != '\0'){
        switch(*format++){
            case 'n':
                //Read a number as a double.
                ;
                double *d = va_arg(valist, double*);
                int i;
                if(!enif_get_double(env, *erl_terms, d))
                    if(enif_get_int(env, *erl_terms, &i))
                        *d = (double) i;
                    else valid = 0;
                break;

            case 'm':
                //Reads a matrix.
                valid = enif_get_matrix(env, *erl_terms, va_arg(valist, Matrix*));
                break;

            case 'i':
                //Reads an int.
                valid = enif_get_int(env, *erl_terms, va_arg(valist, int*));
                break;

            default:
                //Unknown type... give an error.
                valid = 0;
                break;
        }
        erl_terms ++;
    }

    va_end(valist);
    return valid;
}

//----------------------------------------------------------------------------------------------------|
//                        ------------------------------------------                                  |
//                        |                   NIFS                 |                                  |
//                        ------------------------------------------                                  |
//----------------------------------------------------------------------------------------------------|

//@arg 0: List of Lists of numbers.
//@return: Matrix of dimension
ERL_NIF_TERM nif_matrix(ErlNifEnv *env, int argc, const ERL_NIF_TERM argv[]){
    unsigned n_rows, line_length, dest = 0, n_cols = -1;
    ERL_NIF_TERM list = argv[0], row, elem;
    Matrix m;

    //Reading incoming matrix.
    if(!enif_get_list_length(env, list, &n_rows) && n_rows > 0)
        return enif_make_badarg(env);

    for(int i = 0; enif_get_list_cell(env, list, &row, &list); i++){
        if(!enif_get_list_length(env, row, &line_length))
            return enif_make_badarg(env);

        if(n_cols == -1){
            //Allocate binary, make matrix accessor.
            n_cols = line_length;
            m = matrix_alloc(n_rows,n_cols);
        }

        if(n_cols != line_length)
            return enif_make_badarg(env);

        for(int j = 0; enif_get_list_cell(env, row, &elem, &row); j++){
            if(!enif_get_double(env, elem, &m.content[dest])){
                int i;
                if(enif_get_int(env, elem, &i)){
                    m.content[dest] = (double) i;
                }
                else{
                    return enif_make_badarg(env);
                }
            }
            dest++;
        }
    }

    return matrix_to_erl(env, m);
}

#define PRECISION 10

//Used for debug purpose.
void debug_write(char info[]){
    FILE* fp = fopen("debug.txt", "a");
    fprintf(fp, info);
    fprintf(fp, "\n");
    fclose(fp);
}

void debug_write_matrix(Matrix m){
    char *content = enif_alloc(sizeof(char)*((2*m.n_cols-1)*m.n_rows*PRECISION + m.n_rows*2 + 3));
    content[0] = '[';
    content[1] = '\0';
    char converted[PRECISION];

    for(int i=0; i<m.n_rows; i++){
        strcat(content, "[");
        for(int j = 0; j<m.n_cols; j++){
            snprintf(converted, PRECISION-1, "%.5lf", m.content[i*m.n_cols+j]);
            strcat(content, converted);
            if(j != m.n_cols-1)
                strcat(content, " ");
        }
        strcat(content, "]");
    }
    strcat(content, "]");
    debug_write(content);
    enif_free(content);

}

//@arg 0: matrix.
//@arg 1: int, coord m: row
//@arg 2: int, coord n: col
//@return: double, at cord Matrix(m,n).
ERL_NIF_TERM nif_get(ErlNifEnv * env, int argc, const ERL_NIF_TERM argv[]){
    int m,n;
    Matrix matrix;

    if(!enif_get(env, argv, "iim", &m, &n, &matrix))
        return enif_make_badarg(env);
    m--, n--;

    if(m < 0 || m >= matrix.n_rows || n < 0 || n >= matrix.n_cols)
        return enif_make_badarg(env);

    int index = m*matrix.n_cols+n;
    return enif_make_double(env, matrix.content[index]);
}

ERL_NIF_TERM nif_at(ErlNifEnv * env, int argc, const ERL_NIF_TERM argv[]){
    int n;
    Matrix matrix;

    if(!enif_get(env, argv, "mi", &matrix, &n))
        return enif_make_badarg(env);
    n--;

    if( n < 0 || n >= matrix.n_cols * matrix.n_rows)
        return enif_make_badarg(env);

    return enif_make_double(env, matrix.content[n]);
}

//Matrix to flattened list of ints
ERL_NIF_TERM nif_mtfli(ErlNifEnv * env, int argc, const ERL_NIF_TERM argv[]){
    Matrix M;
    if(!enif_get(env, argv, "m", &M)){
        return enif_make_badarg(env);
    }

    int n_elems = M.n_cols * M.n_rows;
    ERL_NIF_TERM *arr = enif_alloc(sizeof(ERL_NIF_TERM)*n_elems);
    for(int i = 0; i<n_elems; i++){
        arr[i] = enif_make_int(env, (int)M.content[i]);
    }

    ERL_NIF_TERM result = enif_make_list_from_array(env, arr, n_elems);
    enif_free(arr);
    return result;
}

//Matrix to flattened list of ints
ERL_NIF_TERM nif_mtfl(ErlNifEnv * env, int argc, const ERL_NIF_TERM argv[]){
    Matrix M;
    if(!enif_get(env, argv, "m", &M)){
        return enif_make_badarg(env);
    }

    int n_elems = M.n_cols * M.n_rows;
    ERL_NIF_TERM *arr = enif_alloc(sizeof(ERL_NIF_TERM)*n_elems);
    for(int i = 0; i<n_elems; i++){
        arr[i] = enif_make_double(env, M.content[i]);
    }

    ERL_NIF_TERM result = enif_make_list_from_array(env, arr, n_elems);
    enif_free(arr);
    return result;
}

//Equal all doubles
//Compares whether all doubles are approximately the same.
int equal_ad(double* a, double* b, int size){
    for(int i = 0; i<size; i++){
        if(fabs(a[i] - b[i])> 1e-6)
            return 0;
    }
    return 1;
}

//@arg 0: Array.
//@arg 1: Array.
//@return: true if arrays share content, false if they have different content || size..
ERL_NIF_TERM nif_equals(ErlNifEnv *env, int argc, const ERL_NIF_TERM argv[]){
    Matrix a,b;

    if(!enif_get(env, argv, "mm", &a, &b))
        return atom_false;

    //Compare number of columns and rows
    if((a.n_cols != b.n_cols || a.n_rows != b.n_rows))
        return atom_false;

    //Compare content of arrays
    if(!equal_ad(a.content, b.content, a.n_cols*a.n_rows))
        return atom_false;

    return atom_true;
}

//@arg 0: int.
//@arg 1: Array.
//@return: returns an array, containing requested row.
ERL_NIF_TERM nif_row(ErlNifEnv * env, int argc, const ERL_NIF_TERM argv[]){
    int row_req;
    Matrix matrix;
    if(!enif_get(env, argv, "im", &row_req, &matrix))
        return enif_make_badarg(env);

    row_req--;
    if(row_req<0 || row_req >= matrix.n_rows)
        return enif_make_badarg(env);

    Matrix row = matrix_alloc(1, matrix.n_cols);
    memcpy(row.content, matrix.content + (row_req * matrix.n_cols), matrix.n_cols*sizeof(double));

    return matrix_to_erl(env, row);
}

//@arg 0: int.
//@arg 1: Array.
//@return: returns an array, containing requested col.
ERL_NIF_TERM nif_col(ErlNifEnv * env, int argc, const ERL_NIF_TERM argv[]){
    int col_req;
    Matrix matrix;
    if(!enif_get(env, argv, "im", &col_req, &matrix))
        return enif_make_badarg(env);

    col_req--;
    if(col_req<0 || col_req >= matrix.n_cols)
        return enif_make_badarg(env);

    Matrix col = matrix_alloc(matrix.n_rows, 1);

    for(int i = 0; i < matrix.n_rows; i++){
        col.content[i] = matrix.content[i * matrix.n_rows + col_req];
    }

    return matrix_to_erl(env, col);
}

//@arg 0: int.
//@arg 1: int.
//@return: empty Matrix of requested dimension..
ERL_NIF_TERM nif_zeros(ErlNifEnv *env, int argc, const ERL_NIF_TERM argv[]){
    int m,n;
    if(!enif_get(env, argv, "ii", &m, &n))
        return enif_make_badarg(env);

    Matrix a = matrix_alloc(m,n);
    memset(a.content, 0, sizeof(double)*m*n);
    return matrix_to_erl(env, a);
}

//@arg 0: int.
//@arg 1: int.
//@return: empty matrix of dimension [arg 0, arg 1]..
ERL_NIF_TERM nif_eye(ErlNifEnv *env, int argc, const ERL_NIF_TERM argv[]){
    int m;
    if(!enif_get_int(env, argv[0], &m))
        return enif_make_badarg(env);

    if(m <= 0)
        return enif_make_badarg(env);

    Matrix a = matrix_alloc(m,m);
    memset(a.content, 0, sizeof(double)*m*m);
    for(int i = 0; i<m; i++){
        a.content[i*m+i] = 1.0;
    }
    return matrix_to_erl(env, a);
}

//Either element-wise multiplication, or multiplication of a matrix by a number.
ERL_NIF_TERM nif_mult(ErlNifEnv *env, int argc, const ERL_NIF_TERM argv[]){

    Matrix a,b;
    double c;

    if(enif_get(env, argv, "mn", &a, &c)){
        Matrix d = matrix_alloc(a.n_rows, a.n_cols);

        for(int i = 0; i<a.n_cols*a.n_rows; i++){
            d.content[i] = a.content[i] * c;
        }
        return matrix_to_erl(env, d);
    }
    else
    if(enif_get(env, argv, "mm", &a, &b)
            && a.n_rows*a.n_cols == b.n_rows*b.n_cols){

        Matrix d = matrix_alloc(a.n_rows, a.n_cols);

        for(int i = 0; i < a.n_rows*a.n_cols; i++){
            d.content[i] = a.content[i] * b.content[i];
        }
        return matrix_to_erl(env, d);
    }

    return enif_make_badarg(env);
}

//Either element-wise addition, or addition of a matrix by a number.
ERL_NIF_TERM nif_add(ErlNifEnv *env, int argc, const ERL_NIF_TERM argv[]){

    Matrix a,b;
    double c;

    if(enif_get(env, argv, "mn", &a, &c)){
        Matrix d = matrix_alloc(a.n_rows, a.n_cols);

        for(int i = 0; i<a.n_cols*a.n_rows; i++){
            d.content[i] = a.content[i] + c;
        }
        return matrix_to_erl(env, d);
    }
    else
    if(enif_get(env, argv, "mm", &a, &b)
            && a.n_rows*a.n_cols == b.n_rows*b.n_cols){

        Matrix d = matrix_alloc(a.n_rows, a.n_cols);

        for(int i = 0; i < a.n_rows*a.n_cols; i++){
            d.content[i] = a.content[i] + b.content[i];
        }
        return matrix_to_erl(env, d);
    }

    return enif_make_badarg(env);
}

//Either element-wise subtraction, or subtraction of a matrix by a number.
ERL_NIF_TERM nif_sub(ErlNifEnv *env, int argc, const ERL_NIF_TERM argv[]){

    Matrix a,b;
    double c;

    if(enif_get(env, argv, "mn", &a, &c)){
        Matrix d = matrix_alloc(a.n_rows, a.n_cols);

        for(int i = 0; i<a.n_cols*a.n_rows; i++){
            d.content[i] = a.content[i] - c;
        }
        return matrix_to_erl(env, d);
    }
    else
    if(enif_get(env, argv, "mm", &a, &b)
            && a.n_rows*a.n_cols == b.n_rows*b.n_cols){

        Matrix d = matrix_alloc(a.n_rows, a.n_cols);

        for(int i = 0; i < a.n_rows*a.n_cols; i++){
            d.content[i] = a.content[i] - b.content[i];
        }
        return matrix_to_erl(env, d);
    }

    return enif_make_badarg(env);
}

//Either element-wise division, or division of a matrix by a number.
ERL_NIF_TERM nif_divide(ErlNifEnv *env, int argc, const ERL_NIF_TERM argv[]){

    Matrix a,b;
    double c;

    if(enif_get(env, argv, "mn", &a, &c)){

        if(c!=0.0){
            Matrix d = matrix_alloc(a.n_rows, a.n_cols);

            for(int i = 0; i<a.n_cols*a.n_rows; i++){
                d.content[i] = a.content[i] / c;
            }
            return matrix_to_erl(env, d);
        }
    }
    else if(enif_get(env, argv, "mm", &a, &b)
            && a.n_rows*a.n_cols == b.n_rows*b.n_cols){

        Matrix d = matrix_alloc(a.n_rows, a.n_cols);

        for(int i = 0; i < a.n_rows*a.n_cols; i++){
            if(b.content[i] == 0.0){
                return enif_make_badarg(env);
            }
            else{
                d.content[i] = a.content[i] / b.content[i];
            }
        }
        return matrix_to_erl(env, d);
    }

    return enif_make_badarg(env);
}

//Transpose of a matrix.
Matrix tr(Matrix a){
    Matrix result = matrix_alloc(a.n_cols ,a.n_rows);

    for(int j = 0; j < a.n_rows; j++){
        for(int i = 0; i < a.n_cols; i++){
            result.content[i*result.n_cols+j] = a.content[j*a.n_cols+i];
        }
    }
    return result;
}
ERL_NIF_TERM nif_transpose(ErlNifEnv *env, int argc, const ERL_NIF_TERM argv[]){

    Matrix a;
    if(!enif_get_matrix(env, argv[0], &a))
        return enif_make_badarg(env);

    return matrix_to_erl(env, tr(a));
}

//arg0: Matrix.
ERL_NIF_TERM nif_inv(ErlNifEnv *env, int argc, const ERL_NIF_TERM argv[]){

    Matrix a;
    if(!enif_get_matrix(env, argv[0], &a))
        return enif_make_badarg(env);

    Matrix inv = matrix_dup(a);

    int N = a.n_rows;
    int *IPIV = enif_alloc(sizeof(int)*N);
    int LWORK = N*N;
    double* WORK = enif_alloc(sizeof(int)*N*N);
    int INFO1, INFO2;

    dgetrf_(&N,&N,inv.content,&N,IPIV,&INFO1);
    dgetri_(&N,inv.content,&N,IPIV,WORK,&LWORK,&INFO2);

    enif_free(IPIV);
    enif_free(WORK);
    ERL_NIF_TERM result;

    if(INFO1 > 0 || INFO2 > 0){
        result = enif_raise_exception(env, enif_make_atom(env, "nif_inv: could not invert singular matrix."));
        matrix_free(inv);
    }
    else if(INFO1 < 0 || INFO2 < 0){
        result = enif_raise_exception(env, enif_make_atom(env, "nif_inv: LAPACK error."));
        matrix_free(inv);
    }
    else result = matrix_to_erl(env, inv);

    return result;
}

//----------------------------------------------------------------------------------------------------|
//                        ------------------------------------------                                  |
//                        |                   CBLAS                |                                  |
//                        ------------------------------------------                                  |
//----------------------------------------------------------------------------------------------------|
//Some CBLAS wrappers.

//Calculates the norm of input vector/matrix, aka the square root of the sum of its composants.
ERL_NIF_TERM nif_dnrm2(ErlNifEnv *env, int argc, const ERL_NIF_TERM argv[]){
    Matrix x;

    if(!enif_get(env, argv, "m", &x)){
        return enif_make_badarg(env);
    }

    double result = cblas_dnrm2(x.n_cols*x.n_rows, x.content, 1);

    return enif_make_double(env, result);
}

//Performs blas_ddot
//Input: two vectors / matrices
ERL_NIF_TERM nif_ddot(ErlNifEnv *env, int argc, const ERL_NIF_TERM argv[]){
    Matrix x,y;

    if(!enif_get(env, argv, "mm", &x, &y)){
        return enif_make_badarg(env);
    }

    int n = fmin(x.n_rows*x.n_cols, y.n_rows*y.n_cols);
    if(n <= 0){
        return enif_make_badarg(env);
    }

    double result = cblas_ddot(n, x.content, 1, y.content, 1);

    return enif_make_double(env, result);
}

//Performs blas_daxpy
//Input: a number, vectors X and Y
//Output: a vector of same dimension then Y, containing alpha X + Y
//------------------------------------------------------------------------
ERL_NIF_TERM nif_daxpy(ErlNifEnv *env, int argc, const ERL_NIF_TERM argv[]){
    Matrix x,y;
    int n;
    double alpha;

    if(!enif_get(env, argv, "inmm", &n, &alpha, &x, &y)){
        return enif_make_badarg(env);
    }

    if(fmin(x.n_rows, x.n_cols) * fmin(y.n_rows, y.n_cols) != 1){
        //We are not using vectors...
        return enif_make_badarg(env);
    }

    Matrix ny = matrix_dup(y);

    cblas_daxpy(n, alpha, x.content, 1, ny.content, 1);

    return matrix_to_erl(env, ny);
}

// Arguments: alpha, A, x, beta, y
// Performs alpha*A*x + beta*y.
// alpha, beta are numbers
// A, x, y are matrices (x and y being vectors)
// x and y are expected to have a length of A.n_cols
ERL_NIF_TERM nif_dgemv(ErlNifEnv *env, int argc, const ERL_NIF_TERM argv[]){
    Matrix A,x,y;
    double alpha, beta;

    if(!enif_get(env, argv, "nmmnm", &alpha, &A, &x, &beta, &y)){
        enif_make_badarg(env);
    }

    //Check dimensions compatibility
    int vec_length = fmin(fmax(x.n_cols, x.n_rows), fmax(y.n_cols, y.n_rows));
    if(vec_length < A.n_cols || fmin(x.n_cols, x.n_rows) != 1 || fmin(y.n_cols, y.n_rows) != 1){
        enif_make_badarg(env);
    }

    Matrix ny = matrix_dup(y);

    cblas_dgemv(CblasRowMajor, CblasNoTrans, A.n_rows, A.n_cols, alpha, A.content, A.n_rows, x.content, 1,  beta, ny.content, 1);

    return matrix_to_erl(env, ny);
}

//Arguments: double alpha, matrix A, matrix B, double beta, matrix C
ERL_NIF_TERM nif_dgemm(ErlNifEnv *env, int argc, const ERL_NIF_TERM argv[]){
    Matrix A, B;
    double alpha = 1.0;
    double beta = 0.0;

    if(!enif_get(env, argv, "mm", &A, &B)
            || A.n_cols != B.n_rows){

        return enif_make_badarg(env);
    }

    Matrix C = matrix_alloc(A.n_rows, B.n_cols);

    cblas_dgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans,
    A.n_rows, B.n_cols, A.n_cols,
    alpha, A.content, A.n_cols, B.content, B.n_rows,
    1.0, C.content, C.n_cols);

    return matrix_to_erl(env, C);
}

ErlNifFunc nif_funcs[] = {
    {"matrix", 1, nif_matrix},
    {"get", 3, nif_get},
    {"at", 2, nif_at},
    {"mtfli", 1, nif_mtfli},
    {"mtfl", 1, nif_mtfl},
    {"equals", 2, nif_equals},
    {"row", 2, nif_row},
    {"col", 2, nif_col},
    {"zeros", 2, nif_zeros},
    {"eye", 1, nif_eye},
    {"mult", 2, nif_mult},
    {"add",  2, nif_add},
    {"sub",  2, nif_sub},
    {"divide", 2, nif_divide},
    {"transpose", 1, nif_transpose},
    {"inv", 1, nif_inv},

    //--- BLAS----------
    {"nrm2", 1, nif_dnrm2},
    {"vec_dot", 2, nif_ddot},
    {"dot", 2, nif_dgemm}
};

ERL_NIF_INIT(numerl, nif_funcs, load, NULL, upgrade, NULL)

ext: c

max_stars_repo_name: tanguyl/numerl

max_issues_repo_name: tanguyl/numerl

max_forks_repo_name: tanguyl/numerl


text:
#include <stdio.h>
#include <stdlib.h>
#include <stdarg.h>
#include <math.h>
#include <string.h>

#include <gsl/gsl_errno.h>
#include <gsl/gsl_odeiv.h>
#include <gsl/gsl_spline.h>
#include <gsl/gsl_integration.h>

#include "ccl.h"
#include "ccl_params.h"

//
// Macros for replacing relative paths
#define EXPAND_STR(s) STRING(s)
#define STRING(s) #s

const ccl_configuration default_config = {ccl_boltzmann_class, ccl_halofit, ccl_nobaryons, ccl_tinker10, ccl_duffy2008, ccl_emu_strict};

const ccl_gsl_params default_gsl_params = {GSL_EPSREL,                          // EPSREL
                                           GSL_N_ITERATION,                     // N_ITERATION
                                           GSL_INTEGRATION_GAUSS_KRONROD_POINTS,// INTEGRATION_GAUSS_KRONROD_POINTS
                                           GSL_EPSREL,                          // INTEGRATION_EPSREL
                                           GSL_INTEGRATION_GAUSS_KRONROD_POINTS,// INTEGRATION_LIMBER_GAUSS_KRONROD_POINTS
                                           GSL_EPSREL,                          // INTEGRATION_LIMBER_EPSREL
                                           GSL_EPSREL_DIST,                     // INTEGRATION_DISTANCE_EPSREL
                                           GSL_EPSREL_DNDZ,                     // INTEGRATION_DNDZ_EPSREL
                                           GSL_EPSREL_SIGMAR,                   // INTEGRATION_SIGMAR_EPSREL
                                           GSL_EPSREL_NU,                       // INTEGRATION_NU_EPSREL
                                           GSL_EPSABS_NU,                       // INTEGRATION_NU_EPSABS
                                           GSL_EPSREL,                          // ROOT_EPSREL
                                           GSL_N_ITERATION,                     // ROOT_N_ITERATION
                                           GSL_EPSREL_GROWTH                    // ODE_GROWTH_EPSREL
                                          };

/* ------- ROUTINE: ccl_cosmology_read_config ------
   INPUTS: none, but will look for ini file in include/ dir
   TASK: fill out global variables of splines with user defined input.
   The variables are defined in ccl_params.h.

   The following are the relevant global variables:
*/

ccl_spline_params * ccl_splines=NULL; // Global variable
ccl_gsl_params * ccl_gsl=NULL; // Global variable

void ccl_cosmology_read_config(void)
{

  int CONFIG_LINE_BUFFER_SIZE=100;
  int MAX_CONFIG_VAR_LEN=100;
  FILE *fconfig;
  char buf[CONFIG_LINE_BUFFER_SIZE];
  char var_name[MAX_CONFIG_VAR_LEN];
  char* rtn;
  double var_dbl;

  // Get parameter .ini filename from environment variable or default location
  const char* param_file;
  const char* param_file_env = getenv("CCL_PARAM_FILE");
  if (param_file_env != NULL) {
    param_file = param_file_env;
  }
  else {
    // Use default ini file
    param_file = EXPAND_STR(__CCL_DATA_DIR__) "/ccl_params.ini";
  }
  if ((fconfig=fopen(param_file, "r")) == NULL) {
    ccl_raise_exception(CCL_ERROR_MISSING_CONFIG_FILE, "ccl_core.c: Failed to open config file: %s", param_file);
    return;
  }

  if(ccl_splines == NULL) {
    ccl_splines = malloc(sizeof(ccl_spline_params));
  }
  if(ccl_gsl == NULL) {
    ccl_gsl = malloc(sizeof(ccl_gsl_params));
    memcpy(ccl_gsl, &default_gsl_params, sizeof(ccl_gsl_params));
  }

  /* Exit gracefully if we couldn't allocate memory */
  if(ccl_splines==NULL || ccl_gsl==NULL) {
    ccl_raise_exception(CCL_ERROR_MEMORY, "ccl_core.c: Failed to allocate memory for config file data.");
    return;
  }

#define MATCH(s, action) if (0 == strcmp(var_name, s)) { action ; continue;} do{} while(0)

  int lineno = 0;
  while(! feof(fconfig)) {
    rtn = fgets(buf, CONFIG_LINE_BUFFER_SIZE, fconfig);
    lineno ++;

    if (buf[0]==';' || buf[0]=='[' || buf[0]=='\n') {
      continue;
    }
    else {
      sscanf(buf, "%99[^=]=%le\n",var_name, &var_dbl);

      // Spline parameters
      MATCH("A_SPLINE_NA", ccl_splines->A_SPLINE_NA=(int) var_dbl);
      MATCH("A_SPLINE_NLOG", ccl_splines->A_SPLINE_NLOG=(int) var_dbl);
      MATCH("A_SPLINE_MINLOG", ccl_splines->A_SPLINE_MINLOG=var_dbl);
      MATCH("A_SPLINE_MIN", ccl_splines->A_SPLINE_MIN=var_dbl);
      MATCH("A_SPLINE_MINLOG_PK", ccl_splines->A_SPLINE_MINLOG_PK=var_dbl);
      MATCH("A_SPLINE_MIN_PK", ccl_splines->A_SPLINE_MIN_PK=var_dbl);
      MATCH("A_SPLINE_MAX", ccl_splines->A_SPLINE_MAX=var_dbl);
      MATCH("LOGM_SPLINE_DELTA", ccl_splines->LOGM_SPLINE_DELTA=var_dbl);
      MATCH("LOGM_SPLINE_NM", ccl_splines->LOGM_SPLINE_NM=(int) var_dbl);
      MATCH("LOGM_SPLINE_MIN", ccl_splines->LOGM_SPLINE_MIN=var_dbl);
      MATCH("LOGM_SPLINE_MAX", ccl_splines->LOGM_SPLINE_MAX=var_dbl);
      MATCH("A_SPLINE_NA_PK", ccl_splines->A_SPLINE_NA_PK=(int) var_dbl);
      MATCH("A_SPLINE_NLOG_PK", ccl_splines->A_SPLINE_NLOG_PK=(int) var_dbl);
      MATCH("K_MAX_SPLINE", ccl_splines->K_MAX_SPLINE=var_dbl);
      MATCH("K_MAX", ccl_splines->K_MAX=var_dbl);
      MATCH("K_MIN", ccl_splines->K_MIN=var_dbl);
      MATCH("N_K", ccl_splines->N_K=(int) var_dbl);

      // 3dcorr parameters
      MATCH("N_K_3DCOR", ccl_splines->N_K_3DCOR=(int) var_dbl);

      // GSL parameters
      MATCH("GSL_EPSREL", ccl_gsl->EPSREL=var_dbl);
      MATCH("GSL_N_ITERATION", ccl_gsl->N_ITERATION=(size_t) var_dbl);
      MATCH("GSL_INTEGRATION_GAUSS_KRONROD_POINTS", ccl_gsl->INTEGRATION_GAUSS_KRONROD_POINTS=(int) var_dbl);
      MATCH("GSL_INTEGRATION_EPSREL", ccl_gsl->INTEGRATION_EPSREL=var_dbl);
      MATCH("GSL_INTEGRATION_DISTANCE_EPSREL", ccl_gsl->INTEGRATION_DISTANCE_EPSREL=var_dbl);
      MATCH("GSL_INTEGRATION_DNDZ_EPSREL", ccl_gsl->INTEGRATION_DNDZ_EPSREL=var_dbl);
      MATCH("GSL_INTEGRATION_SIGMAR_EPSREL", ccl_gsl->INTEGRATION_SIGMAR_EPSREL=var_dbl);
      MATCH("GSL_INTEGRATION_NU_EPSREL", ccl_gsl->INTEGRATION_NU_EPSREL=var_dbl);
      MATCH("GSL_INTEGRATION_NU_EPSABS", ccl_gsl->INTEGRATION_NU_EPSABS=var_dbl);
      MATCH("GSL_INTEGRATION_LIMBER_GAUSS_KRONROD_POINTS", ccl_gsl->INTEGRATION_LIMBER_GAUSS_KRONROD_POINTS=(int) var_dbl);
      MATCH("GSL_INTEGRATION_LIMBER_EPSREL", ccl_gsl->INTEGRATION_LIMBER_EPSREL=var_dbl);
      MATCH("GSL_ROOT_EPSREL", ccl_gsl->ROOT_EPSREL=var_dbl);
      MATCH("GSL_ROOT_N_ITERATION", ccl_gsl->ROOT_N_ITERATION=(int) var_dbl);
      MATCH("GSL_ODE_GROWTH_EPSREL", ccl_gsl->ODE_GROWTH_EPSREL=var_dbl);

      ccl_raise_exception(CCL_ERROR_MISSING_CONFIG_FILE, "ccl_core.c: Failed to parse config file at line %d: %s", lineno, buf);
    }
  }
#undef MATCH

  fclose(fconfig);
}

/* ------- ROUTINE: ccl_cosmology_create ------
INPUTS: ccl_parameters params
        ccl_configuration config
TASK: creates the ccl_cosmology struct and passes some values to it
DEFINITIONS:
chi: comoving distance [Mpc]
growth: growth function (density)
fgrowth: logarithmic derivative of the growth (density) (dlnD/da?)
E: E(a)=H(a)/H0
accelerator: interpolation accelerator for functions of a
accelerator_achi: interpolation accelerator for functions of chi
growth0: growth at z=0, defined to be 1
sigma: ?
p_lin: linear matter power spectrum at z=0?
p_lnl: nonlinear matter power spectrum at z=0?
computed_distances, computed_growth,
computed_power, computed_sigma: store status of the computations
*/
ccl_cosmology * ccl_cosmology_create(ccl_parameters params, ccl_configuration config)
{
  ccl_cosmology * cosmo = malloc(sizeof(ccl_cosmology));
  cosmo->params = params;
  cosmo->config = config;

  cosmo->data.chi = NULL;
  cosmo->data.growth = NULL;
  cosmo->data.fgrowth = NULL;
  cosmo->data.E = NULL;
  cosmo->data.accelerator=NULL;
  cosmo->data.accelerator_achi=NULL;
  cosmo->data.accelerator_m=NULL;
  cosmo->data.accelerator_d=NULL;
  cosmo->data.accelerator_k=NULL;
  cosmo->data.growth0 = 1.;
  cosmo->data.achi=NULL;

  cosmo->data.logsigma = NULL;
  cosmo->data.dlnsigma_dlogm = NULL;

  // hmf parameter for interpolation
  cosmo->data.alphahmf = NULL;
  cosmo->data.betahmf = NULL;
  cosmo->data.gammahmf = NULL;
  cosmo->data.phihmf = NULL;
  cosmo->data.etahmf = NULL;

  cosmo->data.p_lin = NULL;
  cosmo->data.p_nl = NULL;
  //cosmo->data.nu_pspace_int = NULL;
  cosmo->computed_distances = false;
  cosmo->computed_growth = false;
  cosmo->computed_power = false;
  cosmo->computed_sigma = false;
  cosmo->computed_hmfparams = false;
  cosmo->status = 0;
  ccl_cosmology_set_status_message(cosmo, "");

  return cosmo;
}

/* ------ ROUTINE: ccl_parameters_fill_initial -------
INPUT: ccl_parameters: params
TASK: fill parameters not set by ccl_parameters_create with some initial values
DEFINITIONS:
Omega_g = (Omega_g*h^2)/h^2 is the radiation parameter; "g" is for photons, as in CLASS
T_CMB: CMB temperature in Kelvin
Omega_l: Lambda
A_s: amplitude of the primordial PS, enforced here to initially set to NaN
sigma8: variance in 8 Mpc/h spheres for normalization of matter PS, enforced here to initially set to NaN
z_star: recombination redshift
 */
void ccl_parameters_fill_initial(ccl_parameters * params, int *status)
{
  // Fixed radiation parameters
  // Omega_g * h**2 is known from T_CMB
  params->T_CMB =  TCMB;
  // kg / m^3
  double rho_g = 4. * STBOLTZ / pow(CLIGHT, 3) * pow(params->T_CMB, 4);
  // kg / m^3
  double rho_crit = RHO_CRITICAL * SOLAR_MASS/pow(MPC_TO_METER, 3) * pow(params->h, 2);
  params->Omega_g = rho_g/rho_crit;

  // Get the N_nu_rel from Neff and N_nu_mass
  params->N_nu_rel = params->Neff - params->N_nu_mass * pow(TNCDM, 4) / pow(4./11.,4./3.);

  // Temperature of the relativistic neutrinos in K
  double T_nu= (params->T_CMB) * pow(4./11.,1./3.);
  // in kg / m^3
  double rho_nu_rel = params->N_nu_rel* 7.0/8.0 * 4. * STBOLTZ / pow(CLIGHT, 3) * pow(T_nu, 4);
  params-> Omega_n_rel = rho_nu_rel/rho_crit;

  // If non-relativistic neutrinos are present, calculate the phase_space integral.
  if((params->N_nu_mass)>0) {
    // Pass NULL for the accelerator here because we don't have our cosmology object defined yet.
    params->Omega_n_mass = ccl_Omeganuh2(1.0, params->N_nu_mass, params->mnu, params->T_CMB, NULL, status) / ((params->h)*(params->h));
    ccl_check_status_nocosmo(status);
  }
  else{
    params->Omega_n_mass = 0.;
  }

  params->Omega_m = params->Omega_b + params-> Omega_c;
  params->Omega_l = 1.0 - params->Omega_m - params->Omega_g - params->Omega_n_rel -params->Omega_n_mass- params->Omega_k;
  // Initially undetermined parameters - set to nan to trigger
  // problems if they are mistakenly used.
  if (isfinite(params->A_s)) {params->sigma8 = NAN;}
  if (isfinite(params->sigma8)) {params->A_s = NAN;}
  params->z_star = NAN;

  if(fabs(params->Omega_k)<1E-6)
    params->k_sign=0;
  else if(params->Omega_k>0)
    params->k_sign=-1;
  else
    params->k_sign=1;
  params->sqrtk=sqrt(fabs(params->Omega_k))*params->h/CLIGHT_HMPC;
}

/* ------ ROUTINE: ccl_parameters_create -------
INPUT: numbers for the basic cosmological parameters needed by CCL
TASK: fill params with some initial values provided by the user
DEFINITIONS:
Omega_c: cold dark matter
Omega_b: baryons
Omega_m: matter
Omega_k: curvature
little omega_x means Omega_x*h^2
Neff : Effective number of neutrino speces
mnu : Pointer to either sum of neutrino masses or list of three masses.
mnu_type : how the neutrino mass(es) should be treated
w0: Dark energy eq of state parameter
wa: Dark energy eq of state parameter, time variation
H0: Hubble's constant in km/s/Mpc.
h: Hubble's constant divided by (100 km/s/Mpc).
A_s: amplitude of the primordial PS
n_s: index of the primordial PS

 */
ccl_parameters ccl_parameters_create(
                     double Omega_c, double Omega_b, double Omega_k,
				     double Neff, double* mnu, ccl_mnu_convention mnu_type,
				     double w0, double wa, double h, double norm_pk,
				     double n_s, double bcm_log10Mc, double bcm_etab,
				     double bcm_ks, int nz_mgrowth, double *zarr_mgrowth,
				     double *dfarr_mgrowth, int *status)
{
  #ifndef USE_GSL_ERROR
    gsl_set_error_handler_off ();
  #endif

  ccl_parameters params;
  // Initialize params
  params.mnu = NULL;
  params.z_mgrowth=NULL;
  params.df_mgrowth=NULL;
  params.sigma8 = NAN;
  params.A_s = NAN;
  params.Omega_c = Omega_c;
  params.Omega_b = Omega_b;
  params.Omega_k = Omega_k;
  params.Neff = Neff;

  // Set the sum of neutrino masses
  params.sum_nu_masses = *mnu;
  double mnusum = *mnu;
  double *mnu_in = NULL;

  /* Check whether ccl_splines and ccl_gsl exist. If either is not set yet, load
     parameters from the config file. */
  if(ccl_splines==NULL || ccl_gsl==NULL) {
    ccl_cosmology_read_config();
  }

  // Decide how to split sum of neutrino masses between 3 neutrinos. We use
  // a Newton's rule numerical solution (thanks M. Jarvis).

  if (mnu_type==ccl_mnu_sum){
	  // Normal hierarchy

	  mnu_in = malloc(3*sizeof(double));

	  // Check if the sum is zero
	  if (*mnu<1e-15){
		  mnu_in[0] = 0.;
		  mnu_in[1] = 0.;
		  mnu_in[2] = 0.;
	  } else{

	      mnu_in[0] = 0.; // This is a starting guess.

	      double sum_check;
	      // Check that sum is consistent
	      mnu_in[1] = sqrt(DELTAM12_sq);
	      mnu_in[2] = sqrt(DELTAM13_sq_pos);
	      sum_check = mnu_in[0] + mnu_in[1] + mnu_in[2];
	      if (ccl_mnu_sum < sum_check){
		      *status = CCL_ERROR_MNU_UNPHYSICAL;
          }

          double dsdm1;
          // This is the Newton's method
          while (fabs(*mnu - sum_check) > 1e-15){

              dsdm1 = 1. + mnu_in[0] / mnu_in[1] + mnu_in[0] / mnu_in[2];
              mnu_in[0] = mnu_in[0] - (sum_check - *mnu) / dsdm1;
              mnu_in[1] = sqrt(mnu_in[0]*mnu_in[0] + DELTAM12_sq);
              mnu_in[2] = sqrt(mnu_in[0]*mnu_in[0] + DELTAM13_sq_pos);
              sum_check = mnu_in[0] + mnu_in[1] + mnu_in[2];
          }
	  }

  } else if (mnu_type==ccl_mnu_sum_inverted){
	  // Inverted hierarchy

	  mnu_in = malloc(3*sizeof(double));

	  	  // Check if the sum is zero
	  if (*mnu<1e-15){
		  mnu_in[0] = 0.;
		  mnu_in[1] = 0.;
		  mnu_in[2] = 0.;
	  } else{

	      mnu_in[0] = 0.; // This is a starting guess.

	      double sum_check;
	      // Check that sum is consistent
	      mnu_in[1] = sqrt(-1.* DELTAM13_sq_neg - DELTAM12_sq);
	      mnu_in[2] = sqrt(-1.* DELTAM13_sq_neg);
	      sum_check = mnu_in[0] + mnu_in[1] + mnu_in[2];
	      if (ccl_mnu_sum < sum_check){
		      *status = CCL_ERROR_MNU_UNPHYSICAL;
          }

          double dsdm1;
          // This is the Newton's method
          while (fabs(*mnu- sum_check) > 1e-15){
              dsdm1 = 1. + (mnu_in[0] / mnu_in[1]) + (mnu_in[0] / mnu_in[2]);
              mnu_in[0] = mnu_in[0] - (sum_check - *mnu) / dsdm1;
              mnu_in[1] = sqrt(mnu_in[0]*mnu_in[0] + DELTAM12_sq);
              mnu_in[2] = sqrt(mnu_in[0]*mnu_in[0] + DELTAM13_sq_neg);
              sum_check = mnu_in[0] + mnu_in[1] + mnu_in[2];
          }

      }

  } else if (mnu_type==ccl_mnu_sum_equal){
	    // Split the sum of masses equally
	    mnu_in = malloc(3*sizeof(double));
	    mnu_in[0] = params.sum_nu_masses / 3.;
	    mnu_in[1] = params.sum_nu_masses / 3.;
	    mnu_in[2] = params.sum_nu_masses / 3.;
  } else if (mnu_type == ccl_mnu_list){
      // A list of neutrino masses was already passed in
	  params.sum_nu_masses = mnu[0] + mnu[1] + mnu[2];
	  mnu_in = malloc(3*sizeof(double));
	  for(int i=0; i<3; i++) mnu_in[i] = mnu[i];
  } else {
	  *status = CCL_ERROR_NOT_IMPLEMENTED;
  }
  // Check for errors in the neutrino set up (e.g. unphysical mnu)
  ccl_check_status_nocosmo(status);

  // Check which of the neutrino species are non-relativistic today
  int N_nu_mass = 0;
  for(int i = 0; i<3; i=i+1){
  	if (mnu_in[i] > 0.00017){ // Limit taken from Lesgourges et al. 2012
  		N_nu_mass = N_nu_mass + 1;
  	}
  }
  params.N_nu_mass = N_nu_mass;

  // Fill the array of massive neutrinos
  if (N_nu_mass>0){
  	params.mnu = malloc(params.N_nu_mass*sizeof(double));
  	int relativistic[3] = {0, 0, 0};
	for (int i = 0; i < N_nu_mass; i = i + 1){
		for (int j = 0; j<3; j = j +1){
			if ((mnu_in[j]>0.00017) && (relativistic[j]==0)){
				relativistic[j]=1;
				params.mnu[i] = mnu_in[j];
				break;
			}
		} // end loop over neutrinos
	} // end loop over massive neutrinos
  } else{
	  params.mnu = malloc(sizeof(double));
	  params.mnu[0] = 0.;
  }
  // Free mnu_in
  if (mnu_in != NULL) free(mnu_in);

  // Dark Energy
  params.w0 = w0;
  params.wa = wa;

  // Hubble parameters
  params.h = h;
  params.H0 = h*100;

  // Primordial power spectra
  if(norm_pk<1E-5)
    params.A_s=norm_pk;
  else
    params.sigma8=norm_pk;
  params.n_s = n_s;

  //Baryonic params
  if(bcm_log10Mc<0)
    params.bcm_log10Mc=log10(1.2e14);
  else
    params.bcm_log10Mc=bcm_log10Mc;
  if(bcm_etab<0)
    params.bcm_etab=0.5;
  else
    params.bcm_etab=bcm_etab;
  if(bcm_ks<0)
    params.bcm_ks=55.0;
  else
    params.bcm_ks=bcm_ks;

  // Set remaining standard and easily derived parameters
  ccl_parameters_fill_initial(&params, status);

  //Trigger modified growth function if nz>0
  if(nz_mgrowth>0) {
    params.has_mgrowth=true;
    params.nz_mgrowth=nz_mgrowth;
    params.z_mgrowth=malloc(params.nz_mgrowth*sizeof(double));
    params.df_mgrowth=malloc(params.nz_mgrowth*sizeof(double));
    memcpy(params.z_mgrowth,zarr_mgrowth,params.nz_mgrowth*sizeof(double));
    memcpy(params.df_mgrowth,dfarr_mgrowth,params.nz_mgrowth*sizeof(double));
  }
  else {
    params.has_mgrowth=false;
    params.nz_mgrowth=0;
    params.z_mgrowth=NULL;
    params.df_mgrowth=NULL;
  }

  return params;
}

/* ------- ROUTINE: ccl_parameters_create_flat_lcdm --------
INPUT: some cosmological parameters needed to create a flat LCDM model
TASK: call ccl_parameters_create to produce an LCDM model
*/
ccl_parameters ccl_parameters_create_flat_lcdm(double Omega_c, double Omega_b, double h,
                                               double norm_pk, double n_s, int *status)
{
  double Omega_k = 0.0;
  double Neff = 3.046;
  double w0 = -1.0;
  double wa = 0.0;
  double *mnu;
  double mnuval = 0.;  // a pointer to the variable is not kept past the lifetime of this function
  mnu = &mnuval;
  ccl_mnu_convention mnu_type = ccl_mnu_sum;

  ccl_parameters params = ccl_parameters_create(Omega_c, Omega_b, Omega_k, Neff,
        mnu, mnu_type, w0, wa, h, norm_pk, n_s, -1, -1, -1, -1, NULL, NULL, status);

  return params;

}

/**
 * Write a cosmology parameters object to a file in yaml format.
 * @param cosmo Cosmological parameters
 * @param f FILE* pointer opened for reading
 * @return void
 */
void ccl_parameters_write_yaml(ccl_parameters * params, const char * filename, int *status)
{

  FILE * f = fopen(filename, "w");

  if (!f){
    *status = CCL_ERROR_FILE_WRITE;
    return;
  }

#define WRITE_DOUBLE(name) fprintf(f, #name ": %le\n",params->name)
#define WRITE_INT(name) fprintf(f, #name ": %d\n",params->name)

  // Densities: CDM, baryons, total matter, curvature
  WRITE_DOUBLE(Omega_c);
  WRITE_DOUBLE(Omega_b);
  WRITE_DOUBLE(Omega_m);
  WRITE_DOUBLE(Omega_k);
  WRITE_INT(k_sign);

  // Dark Energy
  WRITE_DOUBLE(w0);
  WRITE_DOUBLE(wa);

  // Hubble parameters
  WRITE_DOUBLE(H0);
  WRITE_DOUBLE(h);

  // Neutrino properties
  WRITE_DOUBLE(Neff);
  WRITE_INT(N_nu_mass);
  WRITE_DOUBLE(N_nu_rel);

  if (params->N_nu_mass>0){
    fprintf(f, "mnu: [");
    for (int i=0; i<params->N_nu_mass; i++){
      fprintf(f, "%le, ", params->mnu[i]);
    }
    fprintf(f, "]\n");
  }

  WRITE_DOUBLE(sum_nu_masses);
  WRITE_DOUBLE(Omega_n_mass);
  WRITE_DOUBLE(Omega_n_rel);

  // Primordial power spectra
  WRITE_DOUBLE(A_s);
  WRITE_DOUBLE(n_s);

  // Radiation parameters
  WRITE_DOUBLE(Omega_g);
  WRITE_DOUBLE(T_CMB);

  // BCM baryonic model parameters
  WRITE_DOUBLE(bcm_log10Mc);
  WRITE_DOUBLE(bcm_etab);
  WRITE_DOUBLE(bcm_ks);

  // Derived parameters
  WRITE_DOUBLE(sigma8);
  WRITE_DOUBLE(Omega_l);
  WRITE_DOUBLE(z_star);

  WRITE_INT(has_mgrowth);
  WRITE_INT(nz_mgrowth);

  if (params->has_mgrowth){
    fprintf(f, "z_mgrowth: [");
    for (int i=0; i<params->nz_mgrowth; i++){
      fprintf(f, "%le, ", params->z_mgrowth[i]);
    }
    fprintf(f, "]\n");

    fprintf(f, "df_mgrowth: [");
    for (int i=0; i<params->nz_mgrowth; i++){
      fprintf(f, "%le, ", params->df_mgrowth[i]);
    }
    fprintf(f, "]\n");
  }

#undef WRITE_DOUBLE
#undef WRITE_INT

  fclose(f);

}

/**
 * Write a cosmology parameters object to a file in yaml format.
 * @param cosmo Cosmological parameters
 * @param f FILE* pointer opened for reading
 * @return void
 */
ccl_parameters ccl_parameters_read_yaml(const char * filename, int *status)
{

  FILE * f = fopen(filename, "r");

  if (!f){
    *status = CCL_ERROR_FILE_READ;
    ccl_parameters bad_params;

    ccl_raise_exception(CCL_ERROR_FILE_READ, "ccl_core.c: Failed to read parameters from file.");

    return bad_params;
  }

#define READ_DOUBLE(name) double name; *status |= (0==fscanf(f, #name ": %le\n",&name));
#define READ_INT(name) int name; *status |= (0==fscanf(f, #name ": %d\n",&name))

  // Densities: CDM, baryons, total matter, curvature
  READ_DOUBLE(Omega_c);
  READ_DOUBLE(Omega_b);
  READ_DOUBLE(Omega_m);
  READ_DOUBLE(Omega_k);
  READ_INT(k_sign);

  // Dark Energy
  READ_DOUBLE(w0);
  READ_DOUBLE(wa);

  // Hubble parameters
  READ_DOUBLE(H0);
  READ_DOUBLE(h);

  // Neutrino properties
  READ_DOUBLE(Neff);
  READ_INT(N_nu_mass);
  READ_DOUBLE(N_nu_rel);

  double mnu[3] = {0.0, 0.0, 0.0};
  if (N_nu_mass>0){
    *status |= (0==fscanf(f, "mnu: ["));
    for (int i=0; i<N_nu_mass; i++){
      *status |= (0==fscanf(f, "%le, ", mnu+i));
    }
    *status |= (0==fscanf(f, "]\n"));
  }

  READ_DOUBLE(sum_nu_masses);
  READ_DOUBLE(Omega_n_mass);
  READ_DOUBLE(Omega_n_rel);

  // Primordial power spectra
  READ_DOUBLE(A_s);
  READ_DOUBLE(n_s);

  // Radiation parameters
  READ_DOUBLE(Omega_g);
  READ_DOUBLE(T_CMB);

  // BCM baryonic model parameters
  READ_DOUBLE(bcm_log10Mc);
  READ_DOUBLE(bcm_etab);
  READ_DOUBLE(bcm_ks);

  // Derived parameters
  READ_DOUBLE(sigma8);
  READ_DOUBLE(Omega_l);
  READ_DOUBLE(z_star);

  READ_INT(has_mgrowth);
  READ_INT(nz_mgrowth);

  double *z_mgrowth;
  double *df_mgrowth;

  if (has_mgrowth){
    z_mgrowth = malloc(nz_mgrowth*sizeof(double));
    df_mgrowth = malloc(nz_mgrowth*sizeof(double));
    *status |= (0==fscanf(f, "z_mgrowth: ["));
    for (int i=0; i<nz_mgrowth; i++){
      *status |= (0==fscanf(f, "%le, ", z_mgrowth+i));
    }
    *status |= (0==fscanf(f, "]\n"));

    *status |= (0==fscanf(f, "df_mgrowth: ["));
    for (int i=0; i<nz_mgrowth; i++){
      *status |= (0==fscanf(f, "%le, ", df_mgrowth+i));
    }
    *status |= (0==fscanf(f, "]\n"));
  }
  else{
    z_mgrowth = NULL;
    df_mgrowth = NULL;
  }

#undef READ_DOUBLE
#undef READ_INT

  fclose(f);

  if (status){
    char msg[256];
    snprintf(msg, 256, "ccl_core.c: Structure of YAML file incorrect: %s", filename);
    ccl_raise_exception(*status, msg);
  }

  double norm_pk;

  if (isnan(A_s)){
    norm_pk = sigma8;
  }
  else{
   norm_pk = A_s;
  }

  ccl_parameters params = ccl_parameters_create(
    Omega_c, Omega_b, Omega_k,
    Neff, mnu, ccl_mnu_list,
    w0, wa, h, norm_pk,
    n_s, bcm_log10Mc, bcm_etab,
    bcm_ks, nz_mgrowth, z_mgrowth,
    df_mgrowth, status);

  if(z_mgrowth) free(z_mgrowth);
  if (df_mgrowth) free(df_mgrowth);

  return params;

}

/* ------- ROUTINE: ccl_data_free --------
INPUT: ccl_data
TASK: free the input data
*/
void ccl_data_free(ccl_data * data)
{
  //We cannot assume that all of these have been allocated
  //TODO: it would actually make more sense to do this within ccl_cosmology_free,
  //where we could make use of the flags "computed_distances" etc. to figure out
  //what to free up
  gsl_spline_free(data->chi);
  gsl_spline_free(data->growth);
  gsl_spline_free(data->fgrowth);
  gsl_interp_accel_free(data->accelerator);
  gsl_interp_accel_free(data->accelerator_achi);
  gsl_spline_free(data->E);
  gsl_spline_free(data->achi);
  gsl_spline_free(data->logsigma);
  gsl_spline_free(data->dlnsigma_dlogm);
  gsl_spline2d_free(data->p_lin);
  gsl_spline2d_free(data->p_nl);
  gsl_spline_free(data->alphahmf);
  gsl_spline_free(data->betahmf);
  gsl_spline_free(data->gammahmf);
  gsl_spline_free(data->phihmf);
  gsl_spline_free(data->etahmf);
  gsl_interp_accel_free(data->accelerator_d);
  gsl_interp_accel_free(data->accelerator_m);
  gsl_interp_accel_free(data->accelerator_k);
}

/* ------- ROUTINE: ccl_cosmology_set_status_message --------
INPUT: ccl_cosmology struct, status_string
TASK: set the status message safely.
*/
void ccl_cosmology_set_status_message(ccl_cosmology * cosmo, const char * message, ...)
{
  const int trunc = 480; /* must be < 500 - 4 */
  va_list va;
  va_start(va, message);
  vsnprintf(cosmo->status_message, trunc, message, va);
  va_end(va);

  /* if truncation happens, message[trunc - 1] is not NULL, ... will show up. */
  strcpy(&cosmo->status_message[trunc], "...");
}

/* ------- ROUTINE: ccl_parameters_free --------
INPUT: ccl_parameters struct
TASK: free allocated quantities in the parameters struct
*/
void ccl_parameters_free(ccl_parameters * params)
{
  if (params->mnu != NULL){
    free(params->mnu);
    params->mnu = NULL;
  }
  if (params->z_mgrowth != NULL){
    free(params->z_mgrowth);
    params->z_mgrowth = NULL;
  }
  if (params->df_mgrowth != NULL){
    free(params->df_mgrowth);
    params->df_mgrowth = NULL;
  }
}

/* ------- ROUTINE: ccl_cosmology_free --------
INPUT: ccl_cosmology struct
TASK: free the input data and the cosmology struct
*/
void ccl_cosmology_free(ccl_cosmology * cosmo)
{
  ccl_data_free(&cosmo->data);
  free(cosmo);
}

max_stars_repo_name: Russell-Jones-OxPhys/CCL

max_issues_repo_name: Russell-Jones-OxPhys/CCL

max_forks_repo_name: Russell-Jones-OxPhys/CCL


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

#include <config.h>
#include <stdlib.h>
#include <math.h>

#include <gsl/gsl_math.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_spmatrix.h>

#include "avl.c"

static int compare_triplet(const void *pa, const void *pb, void *param);
static void *avl_spmalloc (size_t size, void *param);
static void avl_spfree (void *block, void *param);

static struct libavl_allocator avl_allocator_spmatrix =
{
  avl_spmalloc,
  avl_spfree
};

/*
gsl_spmatrix_alloc()
  Allocate a sparse matrix in triplet representation

Inputs: n1 - number of rows
        n2 - number of columns

Notes: if (n1,n2) are not known at allocation time, they can each be
set to 1, and they will be expanded as elements are added to the matrix
*/

gsl_spmatrix *
gsl_spmatrix_alloc(const size_t n1, const size_t n2)
{
  const double density = 0.1; /* estimate */
  size_t nzmax = (size_t) floor(n1 * n2 * density);

  if (nzmax == 0)
    nzmax = 10;

  return gsl_spmatrix_alloc_nzmax(n1, n2, nzmax, GSL_SPMATRIX_TRIPLET);
} /* gsl_spmatrix_alloc() */

/*
gsl_spmatrix_alloc_nzmax()
  Allocate a sparse matrix with given nzmax

Inputs: n1     - number of rows
        n2     - number of columns
        nzmax  - maximum number of matrix elements
        sptype - type of matrix (triplet, CCS, CRS)

Notes: if (n1,n2) are not known at allocation time, they can each be
set to 1, and they will be expanded as elements are added to the matrix
*/

gsl_spmatrix *
gsl_spmatrix_alloc_nzmax(const size_t n1, const size_t n2,
                         const size_t nzmax, const size_t sptype)
{
  gsl_spmatrix *m;

  if (n1 == 0)
    {
      GSL_ERROR_NULL ("matrix dimension n1 must be positive integer",
                      GSL_EINVAL);
    }
  else if (n2 == 0)
    {
      GSL_ERROR_NULL ("matrix dimension n2 must be positive integer",
                      GSL_EINVAL);
    }

  m = calloc(1, sizeof(gsl_spmatrix));
  if (!m)
    {
      GSL_ERROR_NULL("failed to allocate space for spmatrix struct",
                     GSL_ENOMEM);
    }

  m->size1 = n1;
  m->size2 = n2;
  m->nz = 0;
  m->nzmax = GSL_MAX(nzmax, 1);
  m->sptype = sptype;

  m->i = malloc(m->nzmax * sizeof(size_t));
  if (!m->i)
    {
      gsl_spmatrix_free(m);
      GSL_ERROR_NULL("failed to allocate space for row indices",
                     GSL_ENOMEM);
    }

  if (sptype == GSL_SPMATRIX_TRIPLET)
    {
      m->tree_data = malloc(sizeof(gsl_spmatrix_tree));
      if (!m->tree_data)
        {
          gsl_spmatrix_free(m);
          GSL_ERROR_NULL("failed to allocate space for AVL tree struct",
                         GSL_ENOMEM);
        }

      m->tree_data->n = 0;

      /* allocate tree data structure */
      m->tree_data->tree = avl_create(compare_triplet, (void *) m,
                                      &avl_allocator_spmatrix);
      if (!m->tree_data->tree)
        {
          gsl_spmatrix_free(m);
          GSL_ERROR_NULL("failed to allocate space for AVL tree",
                         GSL_ENOMEM);
        }

      /* preallocate nzmax tree nodes */
      m->tree_data->node_array = malloc(m->nzmax * sizeof(struct avl_node));
      if (!m->tree_data->node_array)
        {
          gsl_spmatrix_free(m);
          GSL_ERROR_NULL("failed to allocate space for AVL tree nodes",
                         GSL_ENOMEM);
        }

      m->p = malloc(m->nzmax * sizeof(size_t));
      if (!m->p)
        {
          gsl_spmatrix_free(m);
          GSL_ERROR_NULL("failed to allocate space for column indices",
                         GSL_ENOMEM);
        }
    }
  else if (sptype == GSL_SPMATRIX_CCS)
    {
      m->p = malloc((n2 + 1) * sizeof(size_t));
      m->work = malloc(GSL_MAX(n1, n2) *
                       GSL_MAX(sizeof(size_t), sizeof(double)));
      if (!m->p || !m->work)
        {
          gsl_spmatrix_free(m);
          GSL_ERROR_NULL("failed to allocate space for column pointers",
                         GSL_ENOMEM);
        }
    }
  else if (sptype == GSL_SPMATRIX_CRS)
    {
      m->p = malloc((n1 + 1) * sizeof(size_t));
      m->work = malloc(GSL_MAX(n1, n2) *
                       GSL_MAX(sizeof(size_t), sizeof(double)));
      if (!m->p || !m->work)
        {
          gsl_spmatrix_free(m);
          GSL_ERROR_NULL("failed to allocate space for row pointers",
                         GSL_ENOMEM);
        }
    }

  m->data = malloc(m->nzmax * sizeof(double));
  if (!m->data)
    {
      gsl_spmatrix_free(m);
      GSL_ERROR_NULL("failed to allocate space for data",
                     GSL_ENOMEM);
    }

  return m;
} /* gsl_spmatrix_alloc_nzmax() */

/*
gsl_spmatrix_free()
  Free sparse matrix object
*/

void
gsl_spmatrix_free(gsl_spmatrix *m)
{
  if (m->i)
    free(m->i);

  if (m->p)
    free(m->p);

  if (m->data)
    free(m->data);

  if (m->work)
    free(m->work);

  if (m->tree_data)
    {
      if (m->tree_data->tree)
        avl_destroy(m->tree_data->tree, NULL);

      if (m->tree_data->node_array)
        free(m->tree_data->node_array);

      free(m->tree_data);
    }

  free(m);
} /* gsl_spmatrix_free() */

/*
gsl_spmatrix_realloc()
  As elements are added to the sparse matrix, its possible that they
will exceed the previously specified nzmax - reallocate the matrix
with a new nzmax
*/

int
gsl_spmatrix_realloc(const size_t nzmax, gsl_spmatrix *m)
{
  int s = GSL_SUCCESS;
  void *ptr;

  if (nzmax < m->nz)
    {
      GSL_ERROR("new nzmax is less than current nz", GSL_EINVAL);
    }

  ptr = realloc(m->i, nzmax * sizeof(size_t));
  if (!ptr)
    {
      GSL_ERROR("failed to allocate space for row indices", GSL_ENOMEM);
    }

  m->i = (size_t *) ptr;

  if (GSL_SPMATRIX_ISTRIPLET(m))
    {
      ptr = realloc(m->p, nzmax * sizeof(size_t));
      if (!ptr)
        {
          GSL_ERROR("failed to allocate space for column indices", GSL_ENOMEM);
        }

      m->p = (size_t *) ptr;
    }

  ptr = realloc(m->data, nzmax * sizeof(double));
  if (!ptr)
    {
      GSL_ERROR("failed to allocate space for data", GSL_ENOMEM);
    }

  m->data = (double *) ptr;

  /* rebuild binary tree */
  if (GSL_SPMATRIX_ISTRIPLET(m))
    {
      size_t n;

      /* reset tree to empty state, but don't free root tree ptr */
      avl_empty(m->tree_data->tree, NULL);
      m->tree_data->n = 0;

      ptr = realloc(m->tree_data->node_array, nzmax * sizeof(struct avl_node));
      if (!ptr)
        {
          GSL_ERROR("failed to allocate space for AVL tree nodes", GSL_ENOMEM);
        }

      m->tree_data->node_array = ptr;

      /*
       * need to reinsert all tree elements since the m->data addresses
       * have changed
       */
      for (n = 0; n < m->nz; ++n)
        {
          ptr = avl_insert(m->tree_data->tree, &m->data[n]);
          if (ptr != NULL)
            {
              GSL_ERROR("detected duplicate entry", GSL_EINVAL);
            }
        }
    }

  /* update to new nzmax */
  m->nzmax = nzmax;

  return s;
} /* gsl_spmatrix_realloc() */

int
gsl_spmatrix_set_zero(gsl_spmatrix *m)
{
  m->nz = 0;

  if (GSL_SPMATRIX_ISTRIPLET(m))
    {
      /* reset tree to empty state and node index pointer to 0 */
      avl_empty(m->tree_data->tree, NULL);
      m->tree_data->n = 0;
    }

  return GSL_SUCCESS;
} /* gsl_spmatrix_set_zero() */

size_t
gsl_spmatrix_nnz(const gsl_spmatrix *m)
{
  return m->nz;
} /* gsl_spmatrix_nnz() */

/*
gsl_spmatrix_compare_idx()
  Comparison function for searching binary tree in triplet
representation.

To detect duplicate elements in the tree, we want to determine
if there already exists an entry for (i,j) in the tree. Since
the actual tree node stores only the data elements data[n],
we will do pointer arithmetic to get from the given data[n]
to the row/column indices i[n] and j[n].

This compare function will sort the tree first by row i,
and for equal rows, it will then sort by column j

Inputs: ia - row index of element a
        ja - column index of element a
        ib - row index of element b
        jb - column index of element b

Return:
  -1 if pa < pb: (ia,ja) < (ib,jb)
  +1 if pa > pb: (ia,ja) > (ib,jb)
   0 if pa = pb: (ia,ja) == (ib,jb)
*/

int
gsl_spmatrix_compare_idx(const size_t ia, const size_t ja,
                         const size_t ib, const size_t jb)
{
  if (ia < ib)
    return -1;
  else if (ia > ib)
    return 1;
  else
    {
      /* row indices are equal, sort by column index */
      if (ja < jb)
        return -1;
      else if (ja > jb)
        return 1;
      else
        return 0; /* row and column indices are equal */
    }
}

/*
gsl_spmatrix_tree_rebuild()
  When reading a triplet matrix from disk, or when
copying a triplet matrix, it is necessary to rebuild the
binary tree for element searches.

Inputs: m - triplet matrix
*/

int
gsl_spmatrix_tree_rebuild(gsl_spmatrix * m)
{
  if (!GSL_SPMATRIX_ISTRIPLET(m))
    {
      GSL_ERROR("m must be in triplet format", GSL_EINVAL);
    }
  else
    {
      size_t n;

      /* reset tree to empty state, but don't free root tree ptr */
      avl_empty(m->tree_data->tree, NULL);
      m->tree_data->n = 0;

      /* insert all tree elements */
      for (n = 0; n < m->nz; ++n)
        {
          void *ptr = avl_insert(m->tree_data->tree, &m->data[n]);
          if (ptr != NULL)
            {
              GSL_ERROR("detected duplicate entry", GSL_EINVAL);
            }
        }

      return GSL_SUCCESS;
    }
}

/*
compare_triplet()
  Comparison function for searching binary tree in triplet
representation.

To detect duplicate elements in the tree, we want to determine
if there already exists an entry for (i,j) in the tree. Since
the actual tree node stores only the data elements data[n],
we will do pointer arithmetic to get from the given data[n]
to the row/column indices i[n] and j[n].

This compare function will sort the tree first by row i,
and for equal rows, it will then sort by column j

Inputs: pa    - element 1 for comparison (double *)
        pb    - element 2 for comparison (double *)
        param - parameter (gsl_spmatrix)

Return:
  -1 if pa < pb: (ia,ja) < (ib,jb)
  +1 if pa > pb: (ia,ja) > (ib,jb)
   0 if pa = pb: (ia,ja) == (ib,jb)
*/

static int
compare_triplet(const void *pa, const void *pb, void *param)
{
  gsl_spmatrix *m = (gsl_spmatrix *) param;

  /* pointer arithmetic to find indices in data array */
  const size_t idxa = (const double *) pa - m->data;
  const size_t idxb = (const double *) pb - m->data;

  return gsl_spmatrix_compare_idx(m->i[idxa], m->p[idxa],
                                  m->i[idxb], m->p[idxb]);
} /* compare_triplet() */

static void *
avl_spmalloc (size_t size, void *param)
{
  gsl_spmatrix *m = (gsl_spmatrix *) param;

  if (size != sizeof(struct avl_node))
    {
      GSL_ERROR_NULL("attemping to allocate incorrect node size", GSL_EBADLEN);
    }

  /*
   * return the next available avl_node slot; index
   * m->tree_data->n keeps track of next open slot
   */
  if (m->tree_data->n < m->nzmax)
    {
      /* cast to char* for pointer arithmetic */
      unsigned char *node_ptr = (unsigned char *) m->tree_data->node_array;

      /* offset in bytes for next node slot */
      size_t offset = (m->tree_data->n)++ * sizeof(struct avl_node);

      return node_ptr + offset;
    }
  else
    {
      /*
       * we should never get here - gsl_spmatrix_realloc() should
       * be called before exceeding nzmax nodes
       */
      GSL_ERROR_NULL("attemping to allocate tree node past nzmax", GSL_EINVAL);
    }
}

static void
avl_spfree (void *block, void *param)
{
  (void)block;
  (void)param;
  /*
   * do nothing - instead of allocating/freeing individual nodes,
   * we malloc and free nzmax nodes at a time
   */
}

max_issues_repo_name: Brian-ning/HMNE

max_forks_repo_name: Brian-ning/HMNE


text:
#include <math.h>
#include <stdio.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_blas.h>

int
main ()
{
  size_t i,j;

  gsl_matrix *m = gsl_matrix_alloc (10, 10);

  for (i = 0; i < 10; i++)
    for (j = 0; j < 10; j++)
      gsl_matrix_set (m, i, j, sin (i) + cos (j));

  for (j = 0; j < 10; j++)
    {
      gsl_vector_view column = gsl_matrix_column (m, j);
      double d;

      d = gsl_blas_dnrm2 (&column.vector);

      printf ("matrix column %d, norm = %g\n", j, d);
    }

  gsl_matrix_free (m);
}

max_stars_repo_name: zzpwahaha/Chimera-Control-Trim

max_issues_repo_name: zzpwahaha/Chimera-Control-Trim

max_forks_repo_name: zzpwahaha/Chimera-Control-Trim


text:
/* linalg/trimult_complex.c
 *
 * Copyright (C) 2019 Patrick Alken
 *
 * This 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, or (at your option) any
 * later version.
 *
 * This source 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.
 *
 * This module contains code to compute L^T L where L is a lower triangular matrix
 */

#include <config.h>

#include <gsl/gsl_math.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_blas.h>
#include <gsl/gsl_linalg.h>
#include <gsl/gsl_complex.h>
#include <gsl/gsl_complex_math.h>

#include "recurse.h"

static int triangular_multherm_L2(CBLAS_UPLO_t Uplo, gsl_matrix_complex * T);
static int triangular_multherm_L3(CBLAS_UPLO_t Uplo, gsl_matrix_complex * T);
static int triangular_mult_L2(CBLAS_UPLO_t Uplo, gsl_matrix_complex * LU);
static int triangular_mult_L3(CBLAS_UPLO_t Uplo, gsl_matrix_complex * A);
static void complex_conj_vector(gsl_vector_complex * v);

int
gsl_linalg_complex_tri_LHL(gsl_matrix_complex * L)
{
  return triangular_multherm_L3(CblasLower, L);
}

int
gsl_linalg_complex_tri_UL(gsl_matrix_complex * LU)
{
  return triangular_mult_L3(CblasUpper, LU);
}

/*
triangular_multherm_L2()
  Compute L^H L or U U^H

Inputs: Uplo - CblasUpper or CblasLower
        T    - on output the upper (or lower) part of T
               is replaced by L^H L or U U^H

Return: success/error

Notes:
1) Based on LAPACK routine ZLAUU2 using Level 2 BLAS
*/

static int
triangular_multherm_L2(CBLAS_UPLO_t Uplo, gsl_matrix_complex * T)
{
  const size_t N = T->size1;

  if (N != T->size2)
    {
      GSL_ERROR ("matrix must be square", GSL_ENOTSQR);
    }
  else
    {
      size_t i;

      if (Uplo == CblasUpper)
        {
        }
      else
        {
          for (i = 0; i < N; ++i)
            {
              gsl_complex * Tii = gsl_matrix_complex_ptr(T, i, i);
              gsl_complex z0 = *Tii;

              if (i < N - 1)
                {
                  gsl_vector_complex_view v = gsl_matrix_complex_subcolumn(T, i, i + 1, N - i - 1);
                  double norm = gsl_blas_dznrm2(&v.vector);

                  GSL_REAL(*Tii) = gsl_complex_abs2(*Tii) + norm * norm;

                  if (i > 0)
                    {
                      gsl_vector_complex_view w = gsl_matrix_complex_subrow(T, i, 0, i);
                      gsl_matrix_complex_view m = gsl_matrix_complex_submatrix(T, i + 1, 0, N - i - 1, i);

                      complex_conj_vector(&w.vector);
                      gsl_blas_zgemv(CblasConjTrans, GSL_COMPLEX_ONE, &m.matrix, &v.vector, z0, &w.vector);
                      complex_conj_vector(&w.vector);
                    }
                }
              else
                {
                  gsl_vector_complex_view w = gsl_matrix_complex_row(T, i);
                  gsl_blas_zdscal(GSL_REAL(z0), &w.vector);
                }

              GSL_IMAG(*Tii) = 0.0;
            }
        }

      return GSL_SUCCESS;
    }
}

/*
triangular_multherm_L3()
  Compute L^H L or U U^H

Inputs: Uplo - CblasUpper or CblasLower
        T    - on output the upper (or lower) part of T
               is replaced by L^H L or U U^H

Return: success/error

Notes:
1) Based on ReLAPACK using Level 3 BLAS
*/

static int
triangular_multherm_L3(CBLAS_UPLO_t Uplo, gsl_matrix_complex * T)
{
  const size_t N = T->size1;

  if (N != T->size2)
    {
      GSL_ERROR ("matrix must be square", GSL_ENOTSQR);
    }
  else if (N <= CROSSOVER_TRIMULT)
    {
      /* use Level 2 algorithm */
      return triangular_multherm_L2(Uplo, T);
    }
  else
    {
      /*
       * partition matrix:
       *
       * T11 T12
       * T21 T22
       *
       * where T11 is N1-by-N1
       */
      int status;
      const size_t N1 = GSL_LINALG_SPLIT_COMPLEX(N);
      const size_t N2 = N - N1;
      gsl_matrix_complex_view T11 = gsl_matrix_complex_submatrix(T, 0, 0, N1, N1);
      gsl_matrix_complex_view T12 = gsl_matrix_complex_submatrix(T, 0, N1, N1, N2);
      gsl_matrix_complex_view T21 = gsl_matrix_complex_submatrix(T, N1, 0, N2, N1);
      gsl_matrix_complex_view T22 = gsl_matrix_complex_submatrix(T, N1, N1, N2, N2);

      /* recursion on T11 */
      status = triangular_multherm_L3(Uplo, &T11.matrix);
      if (status)
        return status;

      if (Uplo == CblasLower)
        {
          /* T11 += T21^T T21 */
          gsl_blas_zherk(Uplo, CblasConjTrans, 1.0, &T21.matrix, 1.0, &T11.matrix);

          /* T21 = T22^T * T21 */
          gsl_blas_ztrmm(CblasLeft, Uplo, CblasConjTrans, CblasNonUnit, GSL_COMPLEX_ONE, &T22.matrix, &T21.matrix);
        }
      else
        {
          /* T11 += T12 T12^T */
          gsl_blas_zherk(Uplo, CblasNoTrans, 1.0, &T12.matrix, 1.0, &T11.matrix);

          /* T12 = T12 * T22^T */
          gsl_blas_ztrmm(CblasRight, Uplo, CblasConjTrans, CblasNonUnit, GSL_COMPLEX_ONE, &T22.matrix, &T12.matrix);
        }

      /* recursion on T22 */
      status = triangular_multherm_L3(Uplo, &T22.matrix);
      if (status)
        return status;

      return GSL_SUCCESS;
    }
}

/********************************************
 *           INTERNAL ROUTINES              *
 ********************************************/

/*
triangular_mult_L2()
  Compute U L or L U

Inputs: Uplo - CblasUpper or CblasLower (first triangular factor)
        LU   - on input, matrix in LU form;
               on output U*L or L*U

Return: success/error
*/

static int
triangular_mult_L2(CBLAS_UPLO_t Uplo, gsl_matrix_complex * LU)
{
  const size_t N = LU->size1;

  if (N != LU->size2)
    {
      GSL_ERROR ("matrix must be square", GSL_ENOTSQR);
    }
  else
    {
      size_t i;

      /* quick return */
      if (N == 1)
        return GSL_SUCCESS;

      if (Uplo == CblasUpper)
        {
          /* compute U*L and store in LU */

          for (i = 0; i < N; ++i)
            {
              gsl_complex * Aii = gsl_matrix_complex_ptr(LU, i, i);
              gsl_complex Uii = *Aii;

              if (i < N - 1)
                {
                  gsl_vector_complex_view lb = gsl_matrix_complex_subcolumn(LU, i, i + 1, N - i - 1);
                  gsl_vector_complex_view ur = gsl_matrix_complex_subrow(LU, i, i + 1, N - i - 1);
                  gsl_complex dot;

                  gsl_blas_zdotu(&lb.vector, &ur.vector, &dot);
                  *Aii = gsl_complex_add(*Aii, dot);

                  if (i > 0)
                    {
                      gsl_matrix_complex_view U_TR = gsl_matrix_complex_submatrix(LU, 0, i + 1, i, N - i - 1);
                      gsl_matrix_complex_view L_BL = gsl_matrix_complex_submatrix(LU, i + 1, 0, N - i - 1, i);
                      gsl_vector_complex_view ut = gsl_matrix_complex_subcolumn(LU, i, 0, i);
                      gsl_vector_complex_view ll = gsl_matrix_complex_subrow(LU, i, 0, i);

                      gsl_blas_zgemv(CblasTrans, GSL_COMPLEX_ONE, &L_BL.matrix, &ur.vector, Uii, &ll.vector);
                      gsl_blas_zgemv(CblasNoTrans, GSL_COMPLEX_ONE, &U_TR.matrix, &lb.vector, GSL_COMPLEX_ONE, &ut.vector);
                    }
                }
              else
                {
                  gsl_vector_complex_view v = gsl_matrix_complex_subrow(LU, i, 0, i);
                  gsl_blas_zscal(Uii, &v.vector);
                }
            }
        }
      else
        {
        }

      return GSL_SUCCESS;
    }
}

/*
triangular_mult_L3()
  Compute U L or L U

Inputs: Uplo - CblasUpper or CblasLower (for the first triangular factor)
        A    - on input, matrix in LU format;
               on output, U L or L U

Return: success/error
*/

static int
triangular_mult_L3(CBLAS_UPLO_t Uplo, gsl_matrix_complex * A)
{
  const size_t N = A->size1;

  if (N != A->size2)
    {
      GSL_ERROR ("matrix must be square", GSL_ENOTSQR);
    }
  else if (N <= CROSSOVER_TRIMULT)
    {
      return triangular_mult_L2(Uplo, A);
    }
  else
    {
      /* partition matrix:
       *
       * A11 A12
       * A21 A22
       *
       * where A11 is N1-by-N1
       */
      int status;
      const size_t N1 = GSL_LINALG_SPLIT_COMPLEX(N);
      const size_t N2 = N - N1;
      gsl_matrix_complex_view A11 = gsl_matrix_complex_submatrix(A, 0, 0, N1, N1);
      gsl_matrix_complex_view A12 = gsl_matrix_complex_submatrix(A, 0, N1, N1, N2);
      gsl_matrix_complex_view A21 = gsl_matrix_complex_submatrix(A, N1, 0, N2, N1);
      gsl_matrix_complex_view A22 = gsl_matrix_complex_submatrix(A, N1, N1, N2, N2);

      /* recursion on A11 */
      status = triangular_mult_L3(Uplo, &A11.matrix);
      if (status)
        return status;

      if (Uplo == CblasLower)
        {
        }
      else
        {
          /* form U * L */

          /* A11 += A12 A21 */
          gsl_blas_zgemm(CblasNoTrans, CblasNoTrans, GSL_COMPLEX_ONE, &A12.matrix, &A21.matrix, GSL_COMPLEX_ONE, &A11.matrix);

          /* A12 = A12 * L22 */
          gsl_blas_ztrmm(CblasRight, CblasLower, CblasNoTrans, CblasUnit, GSL_COMPLEX_ONE, &A22.matrix, &A12.matrix);

          /* A21 = U22 * A21 */
          gsl_blas_ztrmm(CblasLeft, CblasUpper, CblasNoTrans, CblasNonUnit, GSL_COMPLEX_ONE, &A22.matrix, &A21.matrix);
        }

      /* recursion on A22 */
      status = triangular_mult_L3(Uplo, &A22.matrix);
      if (status)
        return status;

      return GSL_SUCCESS;
    }
}

static void
complex_conj_vector(gsl_vector_complex * v)
{
  size_t i;

  for (i = 0; i < v->size; ++i)
    {
      gsl_complex * vi = gsl_vector_complex_ptr(v, i);
      GSL_IMAG(*vi) = -GSL_IMAG(*vi);
    }
}

max_stars_repo_name: ielomariala/Hex-Game

max_issues_repo_name: ielomariala/Hex-Game

max_forks_repo_name: karanbirsandhu/nu-sense


text:
#ifndef TEST_MATRIX_UTIL
#define TEST_MATRIX_UTIL
#include "matrix_util.h"
//#include <cblas.h>
//int test_blas(int result);

int test_my_zgemm(int result);

int test_inv2by2(int result);

int test_inv2by2inv2by2_overwrite(int result);
int test_outer(int result);

int test_my_zaxpy(int result);

int test_my_zdot(int result);

#endif

ext: h

max_stars_repo_name: ThayaFluss/cnl

max_issues_repo_name: ThayaFluss/cnl

max_forks_repo_name: ThayaFluss/cnl


text:
#ifndef TPEM_H
#define TPEM_H

#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <errno.h>
#ifndef NO_GSL
#include <gsl/gsl_integration.h>
#include <gsl/gsl_errno.h>
#endif

#ifdef	__cplusplus
extern	"C" {
#endif

#define	COMPLEX	1

#ifdef	COMPLEX
typedef	struct	{ double real, imag; } tpem_t;
#else
typedef	double	tpem_t;
#endif	/* COMPLEX */

void	*tpem_malloc	(size_t);
void	tpem_print	(tpem_t a, FILE *b);
tpem_t	tpem_cast	(double, double);
double	tpem_real	(tpem_t);
double tpem_imag (tpem_t z);
double tpem_chebyshev(int m,double x);

/* tpem_sparse.c
 *
 * struct tpem_sparse implements a sparse matrix in compressed row storage
 * format.
 * The CRS format puts the subsequent nonzero elements of the matrix rows
 * in contiguous memory locations. We create 3 vectors: one for tpem_t numbers
 * (*values) and the other two for integers (*colind, *rowptr).
 *
 * The *values vector stores the values of the nonzero elements of the matrix,
 * as they are traversed in a row-wise fashion.
 *
 * The *colind vector stores the column indices of the elements of the *values
 * vector. That is, if
 *	values[k] = a[i][j] then colind[k] = j
 *
 * The *rowptr vector stores the locations in the *values vector that start
 * a row, that is
 *	values[k] = a[i][j] if rowptr[i] <= i < rowptr[i + 1]
 *
 * By convention, we define rowptr[n] = #(nz), the number of nonzeros in the
 * matrix. The storage savings of this approach is significant. Instead of
 * storing n^2 elements, we need only 2#(nz) + n + 1 storage locations.
 *
 * As an example, consider the nonsymmetric matrix defined by
 *	[10  0  0  0 -2  0 ]
 *	[ 3  9  0  0  0  3 ]
 * A =	[ 0  7  8  7  0  0 ]
 *	[ 3  0  8  7  5  0 ]
 *	[ 0  8  0  9  9 13 ]
 *	[ 0  4  0  0  2 -1 ]
 * The CRS format for this matrix is then specified by the arrays given below
 * values = [10 -2  3  9  3  7  8  7  3 ... 9 13  4  2 -1 ]
 * colind = [ 0  4  0  1  5  1  2  3  0 ... 4  5  1  4  5 ]
 * rowptr = [ 0  2  5  8 12 16 19 ]
 *
 * Here's an example illustrating matrix vector product:
 *
 * void tpem_sparse_mult (sparse *matrix, tpem_t *dest, tpem_t *src)
 * {
 *	size_t row, col, k;
 *	tpem_t tmp;
 *
 *	for (row = 0; row < matrix->rank; row++) {
 *		sum = 0.0;
 *		for (k = matrix->rowptr[row]; k < matrix->rowptr[row+1]; k++) {
 *			col  = matrix->colind[k];
 *			sum += matrix->values[k] * src[col];
 *		}
 *		dest[row] = sum;
 *	}
 * }
 *
 * new_tpem_sparse (rank, size) is a constructor
 *	rank:	number of rows (we don't consider degenerate matrices)
 *	size:	#(nz), the total number of nonzero elements
 *
 * tpem_sparse_element (this, row, col)
 *	this:	pointer to the matrix (think pathetic OO)
 *	row:	row index
 *	col:	column index
 *	returns a pointer to a[row][col] or freaks
 *
 */
typedef	struct {
	size_t	rank;		/* number of rows */
	size_t	*rowptr;
	size_t	*colind;
	tpem_t	*values;
} tpem_sparse;

tpem_sparse	*new_tpem_sparse	(size_t, size_t);
void		tpem_sparse_free	(tpem_sparse *);
void		tpem_sparse_copy	(tpem_sparse *, tpem_sparse *);
tpem_t		*tpem_sparse_element	(tpem_sparse *, size_t, size_t);
/* void		tpem_sparse_print	(tpem_sparse *, FILE *); */

void		tpem_matrix_free	(tpem_t **);
tpem_t		**new_tpem_matrix	(size_t, size_t);

/* subspace.c
 *
 * struct tpem_subspace is opaque to the user, don't worry about it. But if
 * you must know, the boolean *flags vector tells whether the i^th state has
 * been used so far, think of the macro
 *	#define	isused(i)	(flags[i] != 0)
 * *stack and *top implement a stack of all states that have been used since
 * the last call to tpem_subspace_reset. The rest should be self-explanatory.
 *
 */
typedef	struct	{
	size_t	size;
	int	*flags;
	size_t	*stack, *top;
} tpem_subspace;

tpem_subspace	*new_tpem_subspace	(size_t);
void		tpem_subspace_free	(tpem_subspace *);
void		tpem_subspace_reset	(tpem_subspace *);
void		tpem_subspace_push	(tpem_subspace *, size_t);
/* void		tpem_subspace_print	(tpem_subspace *, FILE *);
*/

/* tpem.c */
typedef	double	(*tpem_func_ptr)	(double);
void	tpem_off_diagonal_element_tpem	(tpem_sparse *, size_t, tpem_t **,
					size_t, double);
void	tpem_calculate_moment_tpem		(tpem_sparse *, size_t, double *,
					double);
void	tpem_calculate_moment_pem		(tpem_sparse *, size_t, double *);
void	tpem_calculate_moment_diff_tpem	(tpem_sparse *, tpem_sparse *,
					size_t, double *,
					size_t, size_t *, double, double);
void	tpem_calculate_moment_diff_pem	(tpem_sparse *, tpem_sparse *,
					size_t, double *);
void	tpem_calculate_coeffs	(size_t, double *,tpem_func_ptr f);
void	tpem_calculate_coeffs_alt	(size_t, double *,tpem_func_ptr f);
void 	tpem_calculate_coeffs_dos 	(size_t , double *, size_t , size_t , size_t, size_t );
double	tpem_expansion			(size_t, double *, double *);
void	tpem_finalize			(void);
void 	tpem_init			(void);

/*
 * 	CALL TREE:
 * 	----------
 *
 *	tpem_calculate_moment:
 *		\tpem_diagonal_element:
 *			\tpem_sparse_product
 *
 *	tpem_calculate_moment_diff:
 *		|tpem_subspace_for_trace:
 *		|	\tpem_diagonal_element:
 *		|		\tpem_sparse_product
 *		\tpem_diagonal_element:
 *			\tpem_sparse_product
 *
 *	tpem_finalize:
 *		|tpem_diagonal_element
 *		|tpem_calculate_moment
 *		\tpem_calculate_moment_diff
 */

#ifdef	__cplusplus
}
#endif

#endif

max_stars_repo_name: Greatliao/skyrmions

max_issues_repo_name: Greatliao/skyrmions

max_forks_repo_name: Greatliao/skyrmions