text:
#include <config.h>
#include <stdio.h>
#include <math.h>

#include <gsl/gsl_errno.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_blas.h>
#include <gsl/gsl_test.h>

#include <gsl/gsl_wavelet.h>
#include <gsl/gsl_wavelet2d.h>

#define N_BS 11

double urand (void);

double
urand (void)
{
  static unsigned long int x = 1;
  x = (1103515245 * x + 12345) & 0x7fffffffUL;
  return x / 2147483648.0;
}

const size_t member[N_BS] =
  { 309, 307, 305, 303, 301, 208, 206, 204, 202, 105, 103 };

void
test_1d (size_t N, size_t stride, const gsl_wavelet_type * T, size_t member);

void
test_2d (size_t N, size_t tda, const gsl_wavelet_type * T, size_t member, int type);

int
main (int argc, char **argv)
{
  size_t i, N, stride, tda;
  const int S = 1, NS = 2;  /* Standard & Non-standard transforms */

  /* One-dimensional tests */

  for (N = 1; N <= 16384; N *= 2)
    {
      for (stride = 1; stride <= 5; stride++)
        {
          for (i = 0; i < N_BS; i++)
            {
              test_1d (N, stride, gsl_wavelet_bspline, member[i]);
              test_1d (N, stride, gsl_wavelet_bspline_centered, member[i]);
            }

          for (i = 4; i <= 20; i += 2)
            {
              test_1d (N, stride, gsl_wavelet_daubechies, i);
              test_1d (N, stride, gsl_wavelet_daubechies_centered, i);
            }

          test_1d (N, stride, gsl_wavelet_haar, 2);
          test_1d (N, stride, gsl_wavelet_haar_centered, 2);
        }
    }

  /* Two-dimensional tests */

  for (N = 1; N <= 64; N *= 2)
    {
      for (tda = N; tda <= N + 5; tda++)
        {
          for (i = 0; i < N_BS; i++)
            {
              test_2d (N, tda, gsl_wavelet_bspline, member[i], S);
              test_2d (N, tda, gsl_wavelet_bspline_centered, member[i], S);

              test_2d (N, tda, gsl_wavelet_bspline, member[i], NS);
              test_2d (N, tda, gsl_wavelet_bspline_centered, member[i], NS);
            }

          for (i = 4; i <= 20; i += 2)
            {
              test_2d (N, tda, gsl_wavelet_daubechies, i, S);
              test_2d (N, tda, gsl_wavelet_daubechies_centered, i, S);

              test_2d (N, tda, gsl_wavelet_daubechies, i, NS);
              test_2d (N, tda, gsl_wavelet_daubechies_centered, i, NS);
            }

          test_2d (N, tda, gsl_wavelet_haar, 2, S);
          test_2d (N, tda, gsl_wavelet_haar_centered, 2, S);

          test_2d (N, tda, gsl_wavelet_haar, 2, NS);
          test_2d (N, tda, gsl_wavelet_haar_centered, 2, NS);
        }
    }

  exit (gsl_test_summary ());
}

void
test_1d (size_t N, size_t stride, const gsl_wavelet_type * T, size_t member)
{
  gsl_wavelet_workspace *work;
  gsl_vector *v1, *v2, *vdelta;
  gsl_vector_view v;
  gsl_wavelet *w;

  size_t i;
  double *data = malloc (N * stride * sizeof (double));

  for (i = 0; i < N * stride; i++)
    data[i] = 12345.0 + i;

  v = gsl_vector_view_array_with_stride (data, stride, N);
  v1 = &(v.vector);

  for (i = 0; i < N; i++)
    {
      gsl_vector_set (v1, i, urand ());
    }

  v2 = gsl_vector_alloc (N);
  gsl_vector_memcpy (v2, v1);

  vdelta = gsl_vector_alloc (N);

  work = gsl_wavelet_workspace_alloc (N);

  w = gsl_wavelet_alloc (T, member);

  gsl_wavelet_transform_forward (w, v2->data, v2->stride, v2->size, work);
  gsl_wavelet_transform_inverse (w, v2->data, v2->stride, v2->size, work);

  for (i = 0; i < N; i++)
    {
      double x1 = gsl_vector_get (v1, i);
      double x2 = gsl_vector_get (v2, i);
      gsl_vector_set (vdelta, i, fabs (x1 - x2));
    }

  {
    double x1, x2;
    i = gsl_vector_max_index (vdelta);
    x1 = gsl_vector_get (v1, i);
    x2 = gsl_vector_get (v2, i);

    gsl_test (fabs (x2 - x1) > N * 1e-15,
              "%s(%d), n = %d, stride = %d, maxerr = %g",
              gsl_wavelet_name (w), member, N, stride, fabs (x2 - x1));
  }

  if (stride > 1)
    {
      int status = 0;

      for (i = 0; i < N * stride; i++)
        {
          if (i % stride == 0)
            continue;

          status |= (data[i] != (12345.0 + i));
        }

      gsl_test (status, "%s(%d) other data untouched, n = %d, stride = %d",
                gsl_wavelet_name (w), member, N, stride);
    }

  gsl_wavelet_workspace_free (work);
  gsl_wavelet_free (w);
  gsl_vector_free (vdelta);
  gsl_vector_free (v2);
  free (data);
}

void
test_2d (size_t N, size_t tda, const gsl_wavelet_type * T, size_t member, int type)
{
  gsl_wavelet_workspace *work;
  gsl_matrix *m2;
  gsl_wavelet *w;
  gsl_matrix *m1;
  gsl_matrix *mdelta;
  gsl_matrix_view m;
  size_t i;
  size_t j;

  double *data = malloc (N * tda * sizeof (double));

  const char * typename;

  typename = (type == 1) ? "standard" : "nonstd" ;

  for (i = 0; i < N * tda; i++)
    data[i] = 12345.0 + i;

  m = gsl_matrix_view_array_with_tda (data, N, N, tda);
  m1 = &(m.matrix);

  for (i = 0; i < N; i++)
    {
      for (j = 0; j < N; j++)
        {
            gsl_matrix_set (m1, i, j, urand());
        }
    }

  m2 = gsl_matrix_alloc (N, N);
  gsl_matrix_memcpy (m2, m1);

  mdelta = gsl_matrix_alloc (N, N);

  work = gsl_wavelet_workspace_alloc (N);

  w = gsl_wavelet_alloc (T, member);

  switch (type)
    {
    case 1:
      gsl_wavelet2d_transform_matrix_forward (w, m2, work);
      gsl_wavelet2d_transform_matrix_inverse (w, m2, work);
      break;
    case 2:
      gsl_wavelet2d_nstransform_matrix_forward (w, m2, work);
      gsl_wavelet2d_nstransform_matrix_inverse (w, m2, work);
      break;
    }

  for (i = 0; i < N; i++)
    {
      for (j = 0; j < N; j++)
        {
          double x1 = gsl_matrix_get (m1, i, j);
          double x2 = gsl_matrix_get (m2, i, j );
          gsl_matrix_set (mdelta, i, j, fabs (x1 - x2));
        }
    }

  {
    double x1, x2;
    gsl_matrix_max_index (mdelta, &i, &j);
    x1 = gsl_matrix_get (m1, i, j);
    x2 = gsl_matrix_get (m2, i, j);

    gsl_test (fabs (x2 - x1) > N * 1e-15,
              "%s(%d)-2d %s, n = %d, tda = %d, maxerr = %g",
              gsl_wavelet_name (w), member, typename, N, tda, fabs (x2 - x1));
  }

  if (tda > N)
    {
      int status = 0;

      for (i = 0; i < N ; i++)
        {
          for (j = N; j < tda; j++)
            {
              status |= (data[i*tda+j] != (12345.0 + (i*tda+j)));
            }
        }

      gsl_test (status, "%s(%d)-2d %s other data untouched, n = %d, tda = %d",
                gsl_wavelet_name (w), member, typename, N, tda);
    }

  free (data);
  gsl_wavelet_workspace_free (work);
  gsl_wavelet_free (w);
  gsl_matrix_free (m2);
  gsl_matrix_free (mdelta);
}

max_stars_repo_name: manggoguy/parsec-modified

max_issues_repo_name: manggoguy/parsec-modified

max_forks_repo_name: manggoguy/parsec-modified


text:
/*
* -----------------------------------------------------------------
*  Partially Stirred Reaction Library --- pasr_lib.h
*  Version: 1.6180
*  Date: Nov 16, 2010
* -----------------------------------------------------------------
*  Programmer: Americo Barbosa da Cunha Junior
*              americo.cunhajr@gmail.com
* -----------------------------------------------------------------
*  Copyright (c) 2010 by Americo Barbosa da Cunha Junior
*
*  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.
*
*  A copy of the GNU General Public License is available in
*  LICENSE.txt or http://www.gnu.org/licenses/.
* -----------------------------------------------------------------
*  This is the header file for the PaSR model.
* -----------------------------------------------------------------
*/

#ifndef __PaSR_LIB_H__
#define __PaSR_LIB_H__

#include <gsl/gsl_vector.h>
#include <gsl/gsl_rng.h>

/*
* ------------------------------------------------------------
*  Types: struct pasr, pasr_wrk
* ------------------------------------------------------------
* This structure contains fields of PaSR model.
*
* last update: Nov 5, 2009
* ------------------------------------------------------------
*/

typedef struct pasr
{
    int *pasr_idx;         /* particles index vector */
    gsl_vector **pasr_ps;  /* particle system        */
} pasr_wrk;
/*------------------------------------------------------------*/

/*
*------------------------------------------------------------
* function prototypes
*------------------------------------------------------------
*/

void pasr_title();

int pasr_input(unsigned long int *seed,unsigned int *Np,unsigned int *K,
			double *t0,double *delta_t,double *tau_res,
				double *tau_mix,double *Tmax,double *Tmin);

pasr_wrk *pasr_alloc();

int pasr_init(int Np,int Neq,pasr_wrk *pasr);

void pasr_free(int Np,void **pasr_bl);

int pasr_set_all(int Np,gsl_vector *phi,pasr_wrk *pasr);

int pasr_Nexc(int Np,double delta_t,double tau_res);

void pasr_meanvar(int Np,int k,gsl_vector **ps,double *mean,double *var);

void pasr_iem(int Np,int Neq,double delta_t,double tau_mix,
			    double phi_m,int *idx,gsl_vector **ps);

void pmsr_iops(int Np,int Nexc,int Npair,gsl_rng *r,
				    gsl_vector *phi,pmsr_wrk *pmsr);

#endif /* __PaSR_LIB_H__ */

ext: h

max_stars_repo_name: americocunhajr/CRFlowLib

max_issues_repo_name: americocunhajr/CRFlowLib

max_forks_repo_name: americocunhajr/CRFlowLib


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

#include <gsl/gsl_errno.h>
#include <gsl/gsl_odeiv2.h>
#include <gsl/gsl_spline.h>
#include <gsl/gsl_integration.h>
#include <gsl/gsl_roots.h>

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

/* --------- ROUTINE: h_over_h0 ---------
INPUT: scale factor, cosmology
TASK: Compute E(a)=H(a)/H0
*/
static double h_over_h0(double a, ccl_cosmology * cosmo, int *status)
{
  // Check if massive neutrinos are present - if not, we don't need to
  // compute their contribution
  double Om_mass_nu;
  if ((cosmo->params.N_nu_mass)>1e-12) {
    Om_mass_nu = ccl_Omeganuh2(
      a, cosmo->params.N_nu_mass, cosmo->params.mnu, cosmo->params.T_CMB,
      cosmo->data.accelerator, status) / (cosmo->params.h) / (cosmo->params.h);
    ccl_check_status(cosmo, status);
  }
  else {
    Om_mass_nu = 0;
  }

  /* Calculate h^2 using the formula (eqn 2 in the CCL paper):
    E(a)^2 = Omega_m a^-3 +
             Omega_l a^(-3*(1+w0+wa)) exp(3*wa*(a-1)) +
             Omega_k a^-2 +
             (Omega_g + Omega_n_rel) a^-4 +
             Om_mass_nu
  */
  return sqrt(
    (cosmo->params.Omega_m +
     cosmo->params.Omega_l *
       pow(a,-3*(cosmo->params.w0+cosmo->params.wa)) *
       exp(3*cosmo->params.wa*(a-1)) +
     cosmo->params.Omega_k * a +
     (cosmo->params.Omega_g + cosmo->params.Omega_n_rel) / a +
     Om_mass_nu * a*a*a) / (a*a*a));
}

/* --------- ROUTINE: ccl_omega_x ---------
INPUT: cosmology object, scale factor, species label
TASK: Compute the density relative to critical, Omega(a) for a given species.
Possible values for "label":
ccl_species_crit_label <- critical (physical)
ccl_species_m_label <- matter
ccl_species_l_label <- DE
ccl_species_g_label <- radiation
ccl_species_k_label <- curvature
ccl_species_ur_label <- massless neutrinos
ccl_species_nu_label <- massive neutrinos
*/
double ccl_omega_x(ccl_cosmology * cosmo, double a, ccl_species_x_label label, int *status)
{
  // If massive neutrinos are present, compute the phase-space integral and
  // get OmegaNuh2. If not, set OmegaNuh2 to zero.
  double OmNuh2;
  if ((cosmo->params.N_nu_mass) > 0.0001) {
    // Call the massive neutrino density function just once at this redshift.
    OmNuh2 = ccl_Omeganuh2(a, cosmo->params.N_nu_mass, cosmo->params.mnu,
		       cosmo->params.T_CMB, cosmo->data.accelerator, status);
    ccl_check_status(cosmo, status);
  }
  else {
    OmNuh2 = 0.;
  }

  double hnorm = h_over_h0(a, cosmo, status);

  switch(label) {
    case ccl_species_crit_label :
      return 1.;
    case ccl_species_m_label :
      return cosmo->params.Omega_m / (a*a*a) / hnorm / hnorm;
    case ccl_species_l_label :
      return
        cosmo->params.Omega_l *
        pow(a,-3 * (1 + cosmo->params.w0 + cosmo->params.wa)) *
        exp(3 * cosmo->params.wa * (a-1)) / hnorm / hnorm;
    case ccl_species_g_label :
      return cosmo->params.Omega_g / (a*a*a*a) / hnorm / hnorm;
    case ccl_species_k_label :
      return cosmo->params.Omega_k / (a*a) / hnorm / hnorm;
    case ccl_species_ur_label :
      return cosmo->params.Omega_n_rel / (a*a*a*a) / hnorm / hnorm;
    case ccl_species_nu_label :
      return OmNuh2 / (cosmo->params.h) / (cosmo->params.h) / hnorm / hnorm;
    default:
      *status = CCL_ERROR_PARAMETERS;
      ccl_cosmology_set_status_message(
        cosmo, "ccl_background.c: ccl_omega_x(): Species %d not supported\n", label);
      return NAN;
  }
}

/* --------- ROUTINE: ccl_rho_x ---------
INPUT: cosmology object, scale factor, species label
TASK: Compute rho_x(a), with x defined by species label.
Possible values for "label":
ccl_species_crit_label <- critical (physical)
ccl_species_m_label <- matter (physical)
ccl_species_l_label <- DE (physical)
ccl_species_g_label <- radiation (physical)
ccl_species_k_label <- curvature (physical)
ccl_species_ur_label <- massless neutrinos (physical)
ccl_species_nu_label <- massive neutrinos (physical)
*/
double ccl_rho_x(ccl_cosmology * cosmo, double a, ccl_species_x_label label, int is_comoving, int *status)
{
  double comfac;
  if (is_comoving) {
     comfac = a*a*a;
  } else {
      comfac = 1.0;
  }
  double hnorm = h_over_h0(a, cosmo, status);
  double rhocrit =
    RHO_CRITICAL *
    (cosmo->params.h) *
    (cosmo->params.h) * hnorm * hnorm * comfac;

  return rhocrit * ccl_omega_x(cosmo, a, label, status);
}

// Structure to hold parameters of chi_integrand
typedef struct {
  ccl_cosmology *cosmo;
  int * status;
} chipar;

/* --------- ROUTINE: chi_integrand ---------
INPUT: scale factor
TASK: compute the integrand of the comoving distance
*/
static double chi_integrand(double a, void * params_void)
{
  ccl_cosmology * cosmo = ((chipar *)params_void)->cosmo;
  int *status = ((chipar *)params_void)->status;

  return CLIGHT_HMPC/(a*a*h_over_h0(a, cosmo, status));
}

/* --------- ROUTINE: growth_ode_system ---------
INPUT: scale factor
TASK: Define the ODE system to be solved in order to compute the growth (of the density)
*/
static int growth_ode_system(double a,const double y[],double dydt[],void *params)
{
  int status = 0;
  ccl_cosmology * cosmo = params;

  double hnorm=h_over_h0(a,cosmo, &status);
  double om=ccl_omega_x(cosmo, a, ccl_species_m_label, &status);

  dydt[0]=y[1]/(a*a*a*hnorm);
  dydt[1]=1.5*hnorm*a*om*y[0];

  return status;
}

/* --------- ROUTINE: df_integrand ---------
INPUT: scale factor, spline object
TASK: Compute integrand from modified growth function
*/
static double df_integrand(double a,void * spline_void)
{
  if(a<=0)
    return 0;
  else {
    gsl_spline *df_a_spline=(gsl_spline *)spline_void;

    return gsl_spline_eval(df_a_spline,a,NULL)/a;
  }
}

/* --------- ROUTINE: growth_factor_and_growth_rate ---------
INPUT: scale factor, cosmology
TASK: compute the growth (D(z)) and the growth rate, logarithmic derivative (f?)
*/

static int  growth_factor_and_growth_rate(double a,double *gf,double *fg,ccl_cosmology *cosmo, int *stat)
{
  if(a<EPS_SCALEFAC_GROWTH) {
    *gf=a;
    *fg=1;
    return 0;
  }
  else {
    int gslstatus;
    double y[2];
    double ainit=EPS_SCALEFAC_GROWTH;
    gsl_odeiv2_system sys={growth_ode_system,NULL,2,cosmo};
    gsl_odeiv2_driver *d=
      gsl_odeiv2_driver_alloc_y_new(&sys,gsl_odeiv2_step_rkck,0.1*EPS_SCALEFAC_GROWTH,0,ccl_gsl->ODE_GROWTH_EPSREL);

    y[0]=EPS_SCALEFAC_GROWTH;
    y[1]=EPS_SCALEFAC_GROWTH*EPS_SCALEFAC_GROWTH*EPS_SCALEFAC_GROWTH*
      h_over_h0(EPS_SCALEFAC_GROWTH,cosmo, stat);

    gslstatus=gsl_odeiv2_driver_apply(d,&ainit,a,y);
    gsl_odeiv2_driver_free(d);

    if(gslstatus != GSL_SUCCESS) {
      ccl_raise_gsl_warning(gslstatus, "ccl_background.c: growth_factor_and_growth_rate():");
      return NAN;
    }

    *gf=y[0];
    *fg=y[1]/(a*a*h_over_h0(a,cosmo, stat)*y[0]);
    return 0;
  }
}

/* --------- ROUTINE: compute_chi ---------
INPUT: scale factor, cosmology
OUTPUT: chi -> radial comoving distance
TASK: compute radial comoving distance at a
*/
static void compute_chi(double a, ccl_cosmology *cosmo, double * chi, int * stat)
{
  int gslstatus;
  double result;
  chipar p;

  p.cosmo=cosmo;
  p.status=stat;

  gsl_integration_cquad_workspace * workspace = gsl_integration_cquad_workspace_alloc(ccl_gsl->N_ITERATION);
  gsl_function F;
  F.function = &chi_integrand;
  F.params = &p;
  //TODO: CQUAD is great, but slower than other methods. This could be sped up if it becomes an issue.
  gslstatus=gsl_integration_cquad(
    &F, a, 1.0, 0.0, ccl_gsl->INTEGRATION_DISTANCE_EPSREL, workspace, &result, NULL, NULL);
  *chi=result/cosmo->params.h;
  gsl_integration_cquad_workspace_free(workspace);

  if (gslstatus != GSL_SUCCESS) {
    ccl_raise_gsl_warning(gslstatus, "ccl_background.c: compute_chi():");
    *stat = CCL_ERROR_COMPUTECHI;
  }
}

//Root finding for a(chi)
typedef struct {
  double chi;
  ccl_cosmology *cosmo;
  int * status;
} Fpar;

static double fzero(double a,void *params)
{
  double chi,chia,a_use=a;

  chi=((Fpar *)params)->chi;
  compute_chi(a_use,((Fpar *)params)->cosmo,&chia, ((Fpar *)params)->status);

  return chi-chia;
}

static double dfzero(double a,void *params)
{
  ccl_cosmology *cosmo=((Fpar *)params)->cosmo;
  int *stat = ((Fpar *)params)->status;

  chipar p;
  p.cosmo=cosmo;
  p.status=stat;

  return chi_integrand(a,&p)/cosmo->params.h;
}

static void fdfzero(double a,void *params,double *f,double *df)
{
  *f=fzero(a,params);
  *df=dfzero(a,params);
}

/* --------- ROUTINE: a_of_chi ---------
INPUT: comoving distance chi, cosmology, stat, a_old, gsl_root_fdfsolver
OUTPUT: scale factor
TASK: compute the scale factor that corresponds to a given comoving distance chi
Note: This routine uses a root solver to find an a such that compute_chi(a) = chi.
The root solver uses the derivative of compute_chi (which is chi_integrand) and
the value itself.
*/
static void a_of_chi(double chi, ccl_cosmology *cosmo, int* stat, double *a_old, gsl_root_fdfsolver *s)
{
  if(chi==0) {
    *a_old=1;
  }
  else {
    Fpar p;
    gsl_function_fdf FDF;
    double a_previous,a_current=*a_old;

    p.cosmo=cosmo;
    p.chi=chi;
    p.status=stat;
    FDF.f=&fzero;
    FDF.df=&dfzero;
    FDF.fdf=&fdfzero;
    FDF.params=&p;
    gsl_root_fdfsolver_set(s,&FDF,a_current);

    int iter=0, gslstatus;
    do {
      iter++;
      gslstatus=gsl_root_fdfsolver_iterate(s);
      if(gslstatus!=GSL_SUCCESS) ccl_raise_gsl_warning(gslstatus, "ccl_background.c: a_of_chi():");
      a_previous=a_current;
      a_current=gsl_root_fdfsolver_root(s);
      gslstatus=gsl_root_test_delta(a_current, a_previous, 0, ccl_gsl->ROOT_EPSREL);
    } while(gslstatus==GSL_CONTINUE && iter <= ccl_gsl->ROOT_N_ITERATION);

    *a_old=a_current;

    // Allows us to pass a status to h_over_h0 for the neutrino integral calculation.
    if(gslstatus==GSL_SUCCESS) {
      *stat = *(p.status);
    }
    else {
      ccl_raise_gsl_warning(gslstatus, "ccl_background.c: a_of_chi():");
      *stat = CCL_ERROR_COMPUTECHI;
    }
  }
}

/* ----- ROUTINE: ccl_cosmology_compute_distances ------
INPUT: cosmology
TASK: if not already there, make a table of comoving distances and of E(a)
*/

void ccl_cosmology_compute_distances(ccl_cosmology * cosmo, int *status)
{

  //Do nothing if everything is computed already
  if(cosmo->computed_distances)
    return;

  if(ccl_splines->A_SPLINE_MAX>1.) {
    *status = CCL_ERROR_COMPUTECHI;
    ccl_cosmology_set_status_message(cosmo, "ccl_background.c: scale factor cannot be larger than 1.\n");
    return;
  }

  // Create logarithmically and then linearly-spaced values of the scale factor
  int na = ccl_splines->A_SPLINE_NA+ccl_splines->A_SPLINE_NLOG-1;
  double * a = ccl_linlog_spacing(ccl_splines->A_SPLINE_MINLOG, ccl_splines->A_SPLINE_MIN,
				  ccl_splines->A_SPLINE_MAX, ccl_splines->A_SPLINE_NLOG,
				  ccl_splines->A_SPLINE_NA);
  // Allocate arrays for all three of E(a), chi(a), and a(chi)
  double *E_a = malloc(sizeof(double)*na);
  double *chi_a = malloc(sizeof(double)*na);
  // Allocate E(a) and chi(a) splines
  gsl_spline * E = gsl_spline_alloc(A_SPLINE_TYPE, na);
  gsl_spline * chi = gsl_spline_alloc(A_SPLINE_TYPE, na);
  // a(chi) spline allocated below

  //Check for too little memory
  if (a==NULL || E_a==NULL || chi_a==NULL){
    *status=CCL_ERROR_MEMORY;

    ccl_cosmology_set_status_message(cosmo, "ccl_background.c: ccl_cosmology_compute_distances(): ran out of memory\n");
  }

  //Check for messed up scale factor conditions
  if (!*status){
    if ((fabs(a[0]-ccl_splines->A_SPLINE_MINLOG)>1e-5) ||
	(fabs(a[na-1]-ccl_splines->A_SPLINE_MAX)>1e-5) ||
	(a[na-1]>1.0)) {
      *status = CCL_ERROR_LINSPACE;
      ccl_cosmology_set_status_message(cosmo, "ccl_background.c: ccl_cosmology_compute_distances(): Error creating first logarithmic and then linear spacing in a\n");
    }
  }

  // Fill in E(a) - note, this step cannot change the status variable
  if (!*status)
    for (int i=0; i<na; i++)
      E_a[i] = h_over_h0(a[i], cosmo, status);

  // Create a E(a) spline
  if (!*status){
    if (gsl_spline_init(E, a, E_a, na)){
      *status = CCL_ERROR_SPLINE;
      ccl_cosmology_set_status_message(cosmo, "ccl_background.c: ccl_cosmology_compute_distances(): Error creating  E(a) spline\n");
    }
  }

  // Compute chi(a)
  if (!*status){
    for (int i=0; i<na; i++)
      compute_chi(a[i], cosmo, &chi_a[i], status);
    if (*status){
      *status = CCL_ERROR_INTEG;
      ccl_cosmology_set_status_message(cosmo, "ccl_background.c: ccl_cosmology_compute_distances(): chi(a) integration error \n");
    }
  }

  // Initialize chi(a) spline
  if (!*status){
    if (gsl_spline_init(chi, a, chi_a, na)){//in Mpc
      *status = CCL_ERROR_SPLINE;
      ccl_cosmology_set_status_message(cosmo, "ccl_background.c: ccl_cosmology_compute_distances(): Error creating  chi(a) spline\n");
    }
  }

  if (*status){//If there was an error, free the GSL splines and return
    gsl_spline_free(E); //Note: you are allowed to call gsl_free() on NULL
    gsl_spline_free(chi);
  }

  // Set up the boundaries for the a(chi) spline
  double dchi, chi0, chif, a0, af;
  if(!*status){
    dchi=5.;
    chi0=chi_a[na-1];
    chif=chi_a[0];
    a0=a[na-1];
    af=a[0];
  }

  //TODO: The interval in chi (5. Mpc) should be made a macro
  free(a); //Free these, in preparation for making a(chi) splines
  free(E_a);
  free(chi_a);
  //Note: you are allowed to call free() on NULL

  ////////////////////////////////////////////////////////////////////////////////////////////////////
  //Below here na (length of some arrays) changes, so this function has to be split at this point.
  ////////////////////////////////////////////////////////////////////////////////////////////////////

  na = (int)((chif-chi0)/dchi);
  dchi  = (chif-chi0)/na; // <=5, since na is an integer
  //Allocate new arrays for a and chi(a)
  chi_a = ccl_linear_spacing(chi0, chif, na);
  a     = malloc(sizeof(double)*na);
  //Allocate space for GSL root finders
  const gsl_root_fdfsolver_type *T=gsl_root_fdfsolver_newton;
  gsl_root_fdfsolver *s=gsl_root_fdfsolver_alloc(T);
  gsl_spline *achi;

  //Check for too little memory
  if (!*status){
    achi=gsl_spline_alloc(A_SPLINE_TYPE, na);
    if (a==NULL || chi_a==NULL){
      *status=CCL_ERROR_MEMORY;
      ccl_cosmology_set_status_message(cosmo, "ccl_background.c: ccl_cosmology_compute_distances(): ran out of memory\n");
    }else if(fabs(chi_a[0]-chi0)>1e-5 || fabs(chi_a[na-1]-chif)>1e-5) { //Check for messed up chi conditions
      *status = CCL_ERROR_LINSPACE;
      ccl_cosmology_set_status_message(cosmo, "ccl_background.c: ccl_cosmology_compute_distances(): Error creating linear spacing in chi\n");
    }
  }

  // Calculate a(chi)
  if (!*status){
    a[0]=a0; a[na-1]=af;
    for(int i=1;i<na-1;i++) {
      // we are using the previous value as a guess here to help the root finder
      // as long as we use small steps in a this should be fine
      a_of_chi(chi_a[i],cosmo, status, &a0, s);
      a[i]=a0;
    }
    if(*status) {
      *status = CCL_ERROR_ROOT;
      ccl_cosmology_set_status_message(cosmo, "ccl_background.c: ccl_cosmology_compute_distances(): a(chi) root-finding error \n");
    }
  }

  // Initialize the a(chi) spline
  if (!*status){
    if(gsl_spline_init(achi, chi_a, a, na)){
      *status = CCL_ERROR_SPLINE;
      ccl_cosmology_set_status_message(cosmo, "ccl_background.c: ccl_cosmology_compute_distances(): Error creating  a(chi) spline\n");
    }
  }

  free(a);
  free(chi_a); //Note: you are allowed to call free() on NULL
  gsl_root_fdfsolver_free(s);
  if (*status){//If there was an error, free the GSL splines and return
    gsl_spline_free(E); //Note: you are allowed to call gsl_free() on NULL
    gsl_spline_free(chi);
    gsl_spline_free(achi);
    return;
  }

  //If there were no errors, attach the splines to the cosmo struct and end the function.
  if(cosmo->data.accelerator==NULL)
    cosmo->data.accelerator=gsl_interp_accel_alloc();
  cosmo->data.E             = E;
  cosmo->data.chi           = chi;
  cosmo->data.achi          = achi;
  cosmo->computed_distances = true;
}

/* ----- ROUTINE: ccl_cosmology_compute_growth ------
INPUT: cosmology
TASK: if not already there, make a table of growth function and growth rate
      normalize growth to input parameter growth0
*/

void ccl_cosmology_compute_growth(ccl_cosmology * cosmo, int * status)
{

  // This is not valid for massive neutrinos; if we have massive neutrinos, exit.
  if (cosmo->params.N_nu_mass>0){
	  *status = CCL_ERROR_NOT_IMPLEMENTED;
	  ccl_cosmology_set_status_message(cosmo, "ccl_background.c: ccl_cosmology_compute_growth(): Support for the growth rate in cosmologies with massive neutrinos is not yet implemented.\n");
	  return;
  }

  if(cosmo->computed_growth)
    return;

  // Create logarithmically and then linearly-spaced values of the scale factor
  int  chistatus = 0, na = ccl_splines->A_SPLINE_NA+ccl_splines->A_SPLINE_NLOG-1;
  double * a = ccl_linlog_spacing(ccl_splines->A_SPLINE_MINLOG, ccl_splines->A_SPLINE_MIN, ccl_splines->A_SPLINE_MAX, ccl_splines->A_SPLINE_NLOG, ccl_splines->A_SPLINE_NA);
  if (a==NULL ||
      (fabs(a[0]-ccl_splines->A_SPLINE_MINLOG)>1e-5) ||
      (fabs(a[na-1]-ccl_splines->A_SPLINE_MAX)>1e-5) ||
      (a[na-1]>1.0)
      ) {
    free(a);
    *status = CCL_ERROR_LINSPACE;
    ccl_cosmology_set_status_message(cosmo, "ccl_background.c: ccl_cosmology_compute_growth(): Error creating logarithmically and then linear spacing in a\n");
    return;
  }

  gsl_integration_cquad_workspace * workspace=NULL;
  gsl_function F;
  gsl_spline *df_a_spline=NULL;
  if(cosmo->params.has_mgrowth) {
    double *df_arr=malloc(na*sizeof(double));
    if(df_arr==NULL) {
      free(a);
      *status=CCL_ERROR_MEMORY;
      ccl_cosmology_set_status_message(cosmo, "ccl_background.c: ccl_cosmology_compute_distances(): ran out of memory\n");
      return;
    }
    //Generate spline for Delta f(z) that we will then interpolate into an array of a
    gsl_spline *df_z_spline=gsl_spline_alloc(A_SPLINE_TYPE,cosmo->params.nz_mgrowth);
    chistatus=gsl_spline_init(df_z_spline,cosmo->params.z_mgrowth,cosmo->params.df_mgrowth,
			      cosmo->params.nz_mgrowth);

    if(chistatus) {
      free(a);
      free(df_arr);
      gsl_spline_free(df_z_spline);
      *status = CCL_ERROR_SPLINE;
      ccl_cosmology_set_status_message(cosmo, "ccl_background.c: ccl_cosmology_compute_growth(): Error creating Delta f(z) spline\n");
      return;
    }
    for (int i=0; i<na; i++) {
      if(a[i]>0) {
        double z=1./a[i]-1.;

        if(z<=cosmo->params.z_mgrowth[0])
          df_arr[i]=cosmo->params.df_mgrowth[0];
        else if(z>cosmo->params.z_mgrowth[cosmo->params.nz_mgrowth-1])
          df_arr[i]=cosmo->params.df_mgrowth[cosmo->params.nz_mgrowth-1];
        else
	       chistatus |= gsl_spline_eval_e(df_z_spline,z,NULL,&df_arr[i]);
      } else {
	     df_arr[i]=0;
      }
    }
    if(chistatus) {
      free(a);
      free(df_arr);
      gsl_spline_free(df_z_spline);
      *status = CCL_ERROR_SPLINE;
      ccl_cosmology_set_status_message(cosmo, "ccl_background.c: ccl_cosmology_compute_growth(): Error evaluating Delta f(z) spline\n");
      return;
    }
    gsl_spline_free(df_z_spline);

    //Generate Delta(f) spline
    df_a_spline=gsl_spline_alloc(A_SPLINE_TYPE,na);
    chistatus=gsl_spline_init(df_a_spline,a,df_arr,na);
    free(df_arr);
    if (chistatus) {
      free(a);
      gsl_spline_free(df_a_spline);
      *status = CCL_ERROR_SPLINE;
      ccl_cosmology_set_status_message(cosmo, "ccl_background.c: ccl_cosmology_compute_growth(): Error creating Delta f(a) spline\n");
      return;
    }

    workspace=gsl_integration_cquad_workspace_alloc(ccl_gsl->N_ITERATION);
    F.function=&df_integrand;
    F.params=df_a_spline;
  }

  // allocate space for y, which will be all three
  // of E(a), chi(a), D(a) and f(a) in turn.
  int  status_mg=0, gslstatus;
  double growth0,fgrowth0;
  double *y = malloc(sizeof(double)*na);
  if(y==NULL) {
    free(a);
    *status=CCL_ERROR_MEMORY;
    ccl_cosmology_set_status_message(cosmo, "ccl_background.c: ccl_cosmology_compute_distances(): ran out of memory\n");
    return;
  }
  double *y2 = malloc(sizeof(double)*na);
  if(y2==NULL) {
    free(a);
    free(y);
    *status=CCL_ERROR_MEMORY;
    ccl_cosmology_set_status_message(cosmo, "ccl_background.c: ccl_cosmology_compute_distances(): ran out of memory\n");
    return;
  }

  chistatus|=growth_factor_and_growth_rate(1.,&growth0,&fgrowth0,cosmo, status);
  for(int i=0; i<na; i++) {
    chistatus|=growth_factor_and_growth_rate(a[i],&(y[i]),&(y2[i]),cosmo, status);
    if(cosmo->params.has_mgrowth) {
      if(a[i]>0) {
	double df,integ;
	//Add modification to f
	gslstatus = gsl_spline_eval_e(df_a_spline,a[i],NULL,&df);
  if(gslstatus != GSL_SUCCESS) {
    ccl_raise_gsl_warning(gslstatus, "ccl_background.c: ccl_cosmology_compute_growth():");
    status_mg |= gslstatus;
  }
	y2[i]+=df;
	//Multiply D by exp(-int(df))
	gslstatus = gsl_integration_cquad(&F,a[i],1.0,0.0,ccl_gsl->INTEGRATION_DISTANCE_EPSREL,workspace,&integ,NULL,NULL);
	if(gslstatus != GSL_SUCCESS) {
    ccl_raise_gsl_warning(gslstatus, "ccl_background.c: ccl_cosmology_compute_growth():");
    status_mg |= gslstatus;
  }
  y[i]*=exp(-integ);
      }
    }
    y[i]/=growth0;
  }
  if(chistatus || status_mg || *status) {
    free(a);
    free(y);
    free(y2);
    if(df_a_spline!=NULL)
      gsl_spline_free(df_a_spline);
    if(workspace!=NULL)
      gsl_integration_cquad_workspace_free(workspace);
    if (chistatus) {
      *status = CCL_ERROR_INTEG;
      ccl_cosmology_set_status_message(cosmo, "ccl_background.c: ccl_cosmology_compute_growth(): integral for linear growth factor didn't converge\n");
    }
    if(status_mg) {
      *status = CCL_ERROR_INTEG;
      ccl_cosmology_set_status_message(cosmo, "ccl_background.c: ccl_cosmology_compute_growth(): integral for MG growth factor didn't converge\n");
    }
    return;
  }

  if(cosmo->params.has_mgrowth) {
    gsl_spline_free(df_a_spline);
    gsl_integration_cquad_workspace_free(workspace);
  }

  gsl_spline * growth = gsl_spline_alloc(A_SPLINE_TYPE, na);
  chistatus = gsl_spline_init(growth, a, y, na);
  if(chistatus) {
    free(a);
    free(y);
    free(y2);
    gsl_spline_free(growth);
    *status = CCL_ERROR_SPLINE;
    ccl_cosmology_set_status_message(cosmo, "ccl_background.c: ccl_cosmology_compute_growth(): Error creating D(a) spline\n");
    return;
  }

  gsl_spline * fgrowth = gsl_spline_alloc(A_SPLINE_TYPE, na);
  chistatus = gsl_spline_init(fgrowth, a, y2, na);
  if(chistatus) {
    free(a);
    free(y);
    free(y2);
    gsl_spline_free(growth);
    gsl_spline_free(fgrowth);
    *status = CCL_ERROR_SPLINE;
    ccl_cosmology_set_status_message(cosmo, "ccl_background.c: ccl_cosmology_compute_growth(): Error creating f(a) spline\n");
    return;
  }

  // Initialize the accelerator which speeds the splines and
  // assign all the splines we've just made to the structure.
  if(cosmo->data.accelerator==NULL)
    cosmo->data.accelerator=gsl_interp_accel_alloc();
  cosmo->data.growth = growth;
  cosmo->data.fgrowth = fgrowth;
  cosmo->data.growth0 = growth0;
  cosmo->computed_growth = true;

  free(a);
  free(y);
  free(y2);

  return;
}

//Expansion rate normalized to 1 today

double ccl_h_over_h0(ccl_cosmology * cosmo, double a, int* status)
{

  if(!cosmo->computed_distances) {
    ccl_cosmology_compute_distances(cosmo,status);
    ccl_check_status(cosmo, status);
  }

  double h_over_h0;
  int gslstatus = gsl_spline_eval_e(cosmo->data.E, a, cosmo->data.accelerator,&h_over_h0);
  if(gslstatus != GSL_SUCCESS) {
    ccl_raise_gsl_warning(gslstatus, "ccl_background.c: ccl_h_over_h0():");
    *status = gslstatus;
    ccl_cosmology_set_status_message(cosmo, "ccl_background.c: ccl_h_over_h0(): Scale factor outside interpolation range.\n");
    return NAN;
  }

  return h_over_h0;
}

void ccl_h_over_h0s(ccl_cosmology * cosmo, int na, double a[], double output[], int * status)
{
  int _status;

  for (int i=0; i<na; i++) {
    _status = 0;
    output[i] = ccl_h_over_h0(cosmo, a[i], &_status);
    *status |= _status;
  }
}

// Distance-like function examples, all in Mpc
double ccl_comoving_radial_distance(ccl_cosmology * cosmo, double a, int * status)
{
  if((a > (1.0 - 1.e-8)) && (a<=1.0)) {
    return 0.;
  }
  else if(a>1.) {
    *status = CCL_ERROR_COMPUTECHI;
    ccl_cosmology_set_status_message(cosmo, "ccl_background.c: scale factor cannot be larger than 1.\n");
    ccl_check_status(cosmo,status);
    return NAN;
  }
  else {
    if(!cosmo->computed_distances) {
      ccl_cosmology_compute_distances(cosmo, status);
      ccl_check_status(cosmo,status);
    }

    double crd;
    int gslstatus = gsl_spline_eval_e(cosmo->data.chi, a, cosmo->data.accelerator, &crd);
    if(gslstatus != GSL_SUCCESS) {
      ccl_raise_gsl_warning(gslstatus, "ccl_background.c: ccl_comoving_radial_distance():");
      *status = gslstatus;
      ccl_cosmology_set_status_message(cosmo, "ccl_background.c: ccl_comoving_radial_distance(): Scale factor outside interpolation range.\n");
      return NAN;
    }
    return crd;
  }
}

void ccl_comoving_radial_distances(ccl_cosmology * cosmo, int na, double a[], double output[], int* status)
{
  int _status;

  for (int i=0; i<na; i++) {
    _status = 0;
    output[i] = ccl_comoving_radial_distance(cosmo, a[i], &_status);
    *status |= _status;
  }
}

double ccl_sinn(ccl_cosmology *cosmo, double chi, int *status)
{
  //////
  //         { sin(x)  , if k==1
  // sinn(x)={  x      , if k==0
  //         { sinh(x) , if k==-1
  switch(cosmo->params.k_sign) {
  case -1:
    return sinh(cosmo->params.sqrtk * chi) / cosmo->params.sqrtk;
  case 1:
    return sin(cosmo->params.sqrtk*chi) / cosmo->params.sqrtk;
  case 0:
    return chi;
  default:
    *status = CCL_ERROR_PARAMETERS;
    ccl_cosmology_set_status_message(cosmo, "ccl_background.c: ccl_sinn: ill-defined cosmo->params.k_sign = %d",
	    cosmo->params.k_sign);
    return NAN;
  }
}

double ccl_comoving_angular_distance(ccl_cosmology * cosmo, double a, int* status)
{
  if((a > (1.0 - 1.e-8)) && (a<=1.0)) {
    return 0.;
  }
  else if(a>1.) {
    *status = CCL_ERROR_COMPUTECHI;
    ccl_cosmology_set_status_message(cosmo, "ccl_background.c: scale factor cannot be larger than 1.\n");
    ccl_check_status(cosmo,status);
    return NAN;
  }
  else {
    if (!cosmo->computed_distances) {
      ccl_cosmology_compute_distances(cosmo, status);
      ccl_check_status(cosmo, status);
    }

    double chi;
    int gslstatus = gsl_spline_eval_e(cosmo->data.chi, a,
                                      cosmo->data.accelerator,&chi);
    if(gslstatus != GSL_SUCCESS) {
      ccl_raise_gsl_warning(gslstatus, "ccl_background.c: ccl_comoving_angular_distance():");
      *status |= gslstatus;
      ccl_cosmology_set_status_message(cosmo, "ccl_background.c: ccl_comoving_angular_distance(): Scale factor outside interpolation range.\n");
      return NAN;
    }
    return ccl_sinn(cosmo,chi,status);
  }
}

void ccl_comoving_angular_distances(ccl_cosmology * cosmo, int na, double a[],
                                    double output[], int* status)
{
  int _status;

  for (int i=0; i < na; i++) {
    _status = 0;
    output[i] = ccl_comoving_angular_distance(cosmo, a[i], &_status);
    *status |= _status;
  }
}

double ccl_luminosity_distance(ccl_cosmology * cosmo, double a, int* status)
{
  return ccl_comoving_angular_distance(cosmo, a, status) / a;
}

void ccl_luminosity_distances(ccl_cosmology * cosmo, int na, double a[], double output[], int * status)
{
  int _status;

  for (int i=0; i<na; i++) {
    _status = 0;
    output[i] = ccl_luminosity_distance(cosmo, a[i], &_status);
    *status |= _status;
  }
}

double ccl_distance_modulus(ccl_cosmology * cosmo, double a, int* status)
{
  if((a > (1.0 - 1.e-8)) && (a<=1.0)) {
    *status = CCL_ERROR_COMPUTECHI;
    ccl_cosmology_set_status_message(cosmo, "ccl_background.c: distance_modulus undefined for a=1.\n");
    ccl_check_status(cosmo,status);
    return NAN;
  } else if(a>1.) {
    *status = CCL_ERROR_COMPUTECHI;
    ccl_cosmology_set_status_message(cosmo, "ccl_background.c: scale factor cannot be larger than 1.\n");
    ccl_check_status(cosmo,status);
    return NAN;
  } else {
    if (!cosmo->computed_distances) {
      ccl_cosmology_compute_distances(cosmo, status);
      ccl_check_status(cosmo, status);
    }
    /* distance modulus = 5 * log10(d) - 5
       Since d in CCL is in Mpc, you get
         5*log10(10^6) - 5 = 30 - 5 = 25
      for the constant.
    */
    return 5 * log10(ccl_luminosity_distance(cosmo, a, status)) + 25;
  }
}

void ccl_distance_moduli(ccl_cosmology * cosmo, int na, double a[], double output[], int * status)
{
  int _status;

  for (int i=0; i<na; i++) {
    _status = 0;
    output[i] = ccl_distance_modulus(cosmo, a[i], &_status);
    *status |= _status;
  }
}

//Scale factor for a given distance
double ccl_scale_factor_of_chi(ccl_cosmology * cosmo, double chi, int * status)
{
  if((chi < 1.e-8) && (chi>=0.)) {
    return 1.;
  }
  else if(chi<0.) {
    *status = CCL_ERROR_COMPUTECHI;
    ccl_cosmology_set_status_message(cosmo, "ccl_background.c: distance cannot be smaller than 0.\n");
    ccl_check_status(cosmo,status);
    return NAN;
  }
  else {
    if (!cosmo->computed_distances) {
      ccl_cosmology_compute_distances(cosmo,status);
      ccl_check_status(cosmo,status);
    }
    double a;
    int gslstatus = gsl_spline_eval_e(cosmo->data.achi, chi,cosmo->data.accelerator_achi, &a);
    if(gslstatus != GSL_SUCCESS) {
      ccl_raise_gsl_warning(gslstatus, "ccl_background.c: ccl_scale_factor_of_chi():");
      *status |= gslstatus;
    }
    return a;
  }
}

//
void ccl_scale_factor_of_chis(ccl_cosmology * cosmo, int nchi, double chi[], double output[], int * status)
{
  int _status;

  for (int i=0; i<nchi; i++) {
    _status = 0;
    output[i] = ccl_scale_factor_of_chi(cosmo, chi[i], &_status);
    *status |= _status;
  }
}

double ccl_growth_factor(ccl_cosmology * cosmo, double a, int * status)
{
  if(a==1.){
    return 1.;
  }
  else if(a>1.) {
    *status = CCL_ERROR_COMPUTECHI;
    ccl_cosmology_set_status_message(cosmo, "ccl_background.c: scale factor cannot be larger than 1.\n");
    ccl_check_status(cosmo,status);
    return NAN;
  }
  else {
    if (!cosmo->computed_growth) {
      ccl_cosmology_compute_growth(cosmo, status);
      ccl_check_status(cosmo, status);
    }
    if (*status!= CCL_ERROR_NOT_IMPLEMENTED) {
      double D;
      int gslstatus = gsl_spline_eval_e(cosmo->data.growth, a, cosmo->data.accelerator,&D);
      if(gslstatus != GSL_SUCCESS) {
        ccl_raise_gsl_warning(gslstatus, "ccl_background.c: ccl_growth_factor():");
        *status |= gslstatus;
        ccl_cosmology_set_status_message(cosmo, "ccl_background.c: ccl_growth_factor(): Scale factor outside interpolation range.\n");
        return NAN;
      }
      return D;
    }
    else {
      return NAN;
    }
  }
}

void ccl_growth_factors(ccl_cosmology * cosmo, int na, double a[], double output[], int * status)
{
  int _status;

  for (int i=0; i<na; i++) {
    _status = 0;
    output[i] = ccl_growth_factor(cosmo, a[i], &_status);
    *status |= _status;
  }
}

double ccl_growth_factor_unnorm(ccl_cosmology * cosmo, double a, int * status)
{

  if (!cosmo->computed_growth) {
    ccl_cosmology_compute_growth(cosmo, status);
    ccl_check_status(cosmo, status);
  }

  if (*status != CCL_ERROR_NOT_IMPLEMENTED) {
      return cosmo->data.growth0 * ccl_growth_factor(cosmo, a, status);
  } else {
    return NAN;
  }
}

void ccl_growth_factors_unnorm(ccl_cosmology * cosmo, int na, double a[], double output[], int * status)
{
  int _status;

  for (int i=0; i<na; i++) {
    _status = 0;
    output[i] = ccl_growth_factor_unnorm(cosmo, a[i], &_status);
    *status |= _status;
  }
}

double ccl_growth_rate(ccl_cosmology * cosmo, double a, int * status)
{
  if(a>1.) {
    *status = CCL_ERROR_COMPUTECHI;
    ccl_cosmology_set_status_message(cosmo, "ccl_background.c: scale factor cannot be larger than 1.\n");
    ccl_check_status(cosmo,status);
    return NAN;
  } else {
    if (!cosmo->computed_growth) {
      ccl_cosmology_compute_growth(cosmo, status);
      ccl_check_status(cosmo, status);
    }
    if(*status != CCL_ERROR_NOT_IMPLEMENTED) {
      double g;
      int gslstatus = gsl_spline_eval_e(cosmo->data.fgrowth, a, cosmo->data.accelerator,&g);
      if(gslstatus != GSL_SUCCESS) {
        ccl_raise_gsl_warning(gslstatus, "ccl_background.c: ccl_growth_rate():");
        *status |= gslstatus;
        ccl_cosmology_set_status_message(cosmo, "ccl_background.c: ccl_growth_rate(): Scale factor outside interpolation range.\n");
        return NAN;
      }
      return g;
    } else {
	    return NAN;
	  }
  }
}

void ccl_growth_rates(ccl_cosmology * cosmo, int na, double a[], double output[], int * status)
{
  int _status;

  for (int i=0; i<na; i++) {
    _status = 0;
    output[i] = ccl_growth_rate(cosmo, a[i], &_status);
    *status |= _status;
  }
}

ext: c

max_stars_repo_name: vrastil/CCL

max_issues_repo_name: vrastil/CCL

max_forks_repo_name: vrastil/CCL


text:
#ifndef EXECUTE_INCLUDE
#define EXECUTE_INCLUDE

#include <filesystem>
#include <string>
#include <vector>

#include <gsl/string_span>

#include "path.h"

namespace execHelper {
namespace test {
namespace baseUtils {
namespace execution {
using CommandLine = std::vector<std::string>;

int execute(const CommandLine& commandLine,
            const Path& workingDir = std::filesystem::current_path()) noexcept;
} // namespace execution
} // namespace baseUtils
} // namespace test
} // namespace execHelper

#endif /* EXECUTE_INCLUDE */

max_stars_repo_name: exec-helper/source

max_issues_repo_name: exec-helper/source

max_forks_repo_name: exec-helper/source


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

  Library   : Image Registration Toolkit (IRTK)
  Module    : $Id$
  Copyright : Imperial College, Department of Computing
              Visual Information Processing (VIP), 2008 onwards
  Date      : $Date$
  Version   : $Revision$
  Changes   : $Author$

Copyright (c) 1999-2014 and onwards, Imperial College London
All rights reserved.
See LICENSE for details

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

#ifndef _IRTKMATRIX_H

#define _IRTKMATRIX_H

#include <float.h>

#ifdef USE_VXL
#include <vnl/vnl_matrix.h>
#else
#include <gsl/gsl_matrix.h>
#endif

#include <irtkObject.h>

#include <iostream>
#include <cmath>

/**

  Matrix class.

*/

class irtkMatrix : public irtkObject
{

protected:

  /// Number of rows
  int _rows;

  /// Number of colums
  int _cols;

  /// Data
  double **_matrix;

public:

  /// Default constructor
  irtkMatrix();

  /// Constructor for given number of rows and columns
  irtkMatrix(int, int);

  /// Copy constructor
  irtkMatrix(const irtkMatrix &);

  /// Destructor
  ~irtkMatrix();

  /// Initialize matrix with number of rows and columns
  void Initialize(int, int);

  //
  // Matrix access functions
  //

  /// Returns number of rows
  int Rows() const;

  /// Returns number of columns
  int Cols() const;

  /// Puts matrix value
  void   Put(int, int, double);

  /// Gets matrix value
  double Get(int, int) const;

  //
  // Operators for matrix access
  //

  /// Puts matrix value
  double& operator()(int, int);

  /// Gets matrix value
  double  operator()(int, int) const;

  /// Access matrix get operator
  irtkMatrix  operator()(int, int, int, int);

  /// Access matrix put operator
  void    operator()(irtkMatrix &, int, int);

  //
  // Matrix operators for doubles
  //

  /// Subtraction of a double
  irtkMatrix& operator-=(const double&);

  /// Addition of a double
  irtkMatrix& operator+=(const double&);

  /// Multiplication with a double
  irtkMatrix& operator*=(const double&);

  /// Division by a double
  irtkMatrix& operator/=(const double&);

  /// Return result of subtraction of a double
  irtkMatrix  operator- (const double&);

  /// Return result of addition of a double
  irtkMatrix  operator+ (const double&);

  /// Return result of multiplication with a double
  irtkMatrix  operator* (const double&);

  /// Return result of division by a double
  irtkMatrix  operator/ (const double&);

  //
  // Matrix operators for matrices
  //

  /// Matrix copy operator
  irtkMatrix& operator =(const irtkMatrix&);

  /// Matrix subtraction operator
  irtkMatrix& operator-=(const irtkMatrix&);

  /// Matrix addition operator
  irtkMatrix& operator+=(const irtkMatrix&);

  /// Matrix multiplication operator
  irtkMatrix& operator*=(const irtkMatrix&);

  /// Return result of matrix subtraction
  irtkMatrix  operator- (const irtkMatrix&);

  /// Return result of matrix addition
  irtkMatrix  operator+ (const irtkMatrix&);

  /// Return result of matrix multiplication
  irtkMatrix  operator* (const irtkMatrix&);

  /// Matrix inversion operator
  irtkMatrix  operator! (void);

  /// Matrix transpose operator
  irtkMatrix  operator~ (void);

  /// Matrix comparison operator =
  int operator==(const irtkMatrix &);

  // Matrix exponential via Pade approximation.
  // See Golub and Van Loan, Matrix Computations, Algorithm 11.3-1.
  friend irtkMatrix expm(irtkMatrix);

  // Matrix logarithm.
  friend irtkMatrix logm(irtkMatrix);

  // Matrix square root.
  friend irtkMatrix sqrtm(irtkMatrix);

  friend irtkMatrix FrechetMean(irtkMatrix *, int, int = 10);

  friend irtkMatrix FrechetMean(irtkMatrix *, double *, int, int = 10);

#ifndef USE_STL
  /// Comparison operator != (if USE_STL is defined, negate == operator)
  int operator!=(const irtkMatrix &);
#endif

  //
  // Matrix functions
  //

  /// Calculate norm of matrix
  double Norm(void) const;

  /// Calculate trace of matrix
  double Trace(void) const;

  // The infinity norm is the maximum of the absolute value row sums.
  double InfinityNorm(void) const;

  /// Calculate determinant of matrix
  double Det() const;

  /* Calculate SVD of matrix
   * now does the same as Eigenvalues (but eigenvalues are in square root)
   * algorithm chosen according to matrix properties (symmetry, diagonizability)
   */
  void   SVD(irtkMatrix &, irtkVector &, irtkMatrix &);

  /// Identity matrix
  void   Ident();

  /// Returns true if the matrix is an identity matrix.
  bool IsIdentity() const;

  /// Invert of matrix
  void   Invert();

  /// Adjugate of matrix and return determine;
  void   Adjugate(double &);

  /// Transpose matrix
  void   Transpose();

  /* Calculate eigenvalues and eigenvectors of matrix
   * now does the same as SVD (but eigenvalues are NOT in square root)
   * algorithm chosen according to matrix properties (symmetry, diagonizability)
   */
  void   Eigenvalues(irtkMatrix &, irtkVector &, irtkMatrix &);

  /// Calculate least square fit via SVD
  void   LeastSquaresFit(const irtkVector &, irtkVector &);

  //
  // Matrix in- and output
  //

  /// Interface to output stream
  friend std::ostream& operator<< (std::ostream&, const irtkMatrix&);

  /// Interface to input stream
  friend std::istream& operator>> (std::istream&, irtkMatrix&);

  /// Print matrix
  void Print();

  /// Read matrix from file
  void Read (char *);

  /// Write matrix to file
  void Write(char *);

  /// Import matrix from text file (requires no. of expected rows and cols)
  void Import (char *, int, int);

#ifdef USE_VXL

  /// Conversion to VNL matrix
  template <class T> void Matrix2Vnl(vnl_matrix<T> *m) const;

  /// Conversion from VNL matrix
  template <class T> void Vnl2Matrix(vnl_matrix<T> *m);

#else

  /// Conversion to GSL matrix (memory must be allocated)
  void Matrix2GSL(gsl_matrix *) const;

  /// Conversion from GSL matrix
  void GSL2Matrix(gsl_matrix *);

  //
  // Matrix operators for vectors
  //

  /// Return result of multiplication of matrix and vector
  irtkVector  operator* (const irtkVector&);

#endif
};

//
// Access operators
//

inline int irtkMatrix::Rows() const
{
  return _rows;
}

inline int irtkMatrix::Cols() const
{
  return _cols;
}

inline void irtkMatrix::Put(int rows, int cols, double matrix)
{
#ifdef NO_BOUNDS
  _matrix[cols][rows] = matrix;
#else
  if ((rows >= 0) && (rows < _rows) && (cols >= 0) && (cols < _cols)) {
    _matrix[cols][rows] = matrix;
  } else {
    cout << "irtkMatrix::Put: parameter out of range\n";
  }
#endif
}

inline double irtkMatrix::Get(int rows, int cols) const
{
#ifdef NO_BOUNDS
  return _matrix[cols][rows];
#else
  if ((rows >= 0) && (rows < _rows) && (cols >= 0) && (cols < _cols)) {
    return _matrix[cols][rows];
  } else {
    cout << "irtkMatrix::Get: parameter out of range\n";
    return 0;
  }
#endif
}

inline double &irtkMatrix::operator()(int rows, int cols)
{
#ifdef NO_BOUNDS
  return _matrix[cols][rows];
#else
  if ((rows >= 0) && (rows < _rows) && (cols >= 0) && (cols < _cols)) {
    return _matrix[cols][rows];
  } else {
    cout << "irtkMatrix::operator(): parameter out of range\n";
    return _matrix[0][0];
  }
#endif
}

inline double irtkMatrix::operator()(int rows, int cols) const
{
#ifdef NO_BOUNDS
  return _matrix[cols][rows];
#else
  if ((rows >= 0) && (rows < _rows) && (cols >= 0) && (cols < _cols)) {
    return _matrix[cols][rows];
  } else {
    cout << "irtkMatrix::operator(): parameter out of range\n";
    return 0;
  }
#endif
}

//
// Matrix operators for doubles
//

inline irtkMatrix& irtkMatrix::operator-=(const double &x)
{
  int i, j;

  for (j = 0; j < _cols; j++) {
    for (i = 0; i < _rows; i++) {
      _matrix[j][i] -= x;
    }
  }
  return *this;
}

inline irtkMatrix& irtkMatrix::operator+=(const double &x)
{
  int i, j;

  for (j = 0; j < _cols; j++) {
    for (i = 0; i < _rows; i++) {
      _matrix[j][i] += x;
    }
  }
  return *this;
}

inline irtkMatrix& irtkMatrix::operator*=(const double &x)
{
  int i, j;

  for (j = 0; j < _cols; j++) {
    for (i = 0; i < _rows; i++) {
      _matrix[j][i] *= x;
    }
  }
  return *this;
}

inline irtkMatrix& irtkMatrix::operator/=(const double &x)
{
  int i, j;

  for (j = 0; j < _cols; j++) {
    for (i = 0; i < _rows; i++) {
      _matrix[j][i] /= x;
    }
  }
  return *this;
}

inline irtkMatrix irtkMatrix::operator-(const double &x)
{
  int i, j;
  irtkMatrix m;

  m = *this;
  for (j = 0; j < _cols; j++) {
    for (i = 0; i < m._rows; i++) {
      m._matrix[j][i] = _matrix[j][i] - x;
    }
  }
  return m;
}

inline irtkMatrix irtkMatrix::operator+(const double &x)
{
  int i, j;
  irtkMatrix m;

  m = *this;
  for (j = 0; j < _cols; j++) {
    for (i = 0; i < m._rows; i++) {
      m._matrix[j][i] = _matrix[j][i] + x;
    }
  }
  return m;
}

inline irtkMatrix irtkMatrix::operator*(const double &x)
{
  int i, j;
  irtkMatrix m;

  m = *this;
  for (j = 0; j < _cols; j++) {
    for (i = 0; i < m._rows; i++) {
      m._matrix[j][i] = _matrix[j][i] * x;
    }
  }
  return m;
}

inline irtkMatrix irtkMatrix::operator/(const double &x)
{
  int i, j;
  irtkMatrix m;

  m = *this;
  for (j = 0; j < _cols; j++) {
    for (i = 0; i < m._rows; i++) {
      m._matrix[j][i] = _matrix[j][i] / x;
    }
  }
  return m;
}

//
// Matrix operators for matrices
//

inline int irtkMatrix::operator==(const irtkMatrix& m)
{
  int i, j;

  if ((m._rows != _rows) || (m._cols != _cols)) return 0;
  for (j = 0; j < _cols; j++) {
    for (i = 0; i < m._rows; i++) {
      if (m._matrix[j][i] != _matrix[j][i]) return 0;
    }
  }
  return 1;
}

#ifndef USE_STL
inline int irtkMatrix::operator!=(const irtkMatrix& m)
{
  return !(*this == m);
}
#endif

inline double irtkMatrix::Trace(void) const
{
    int i, j;
    double trace = 0;

    if(_rows == _cols){

        // The trace of a matrix
        for (j = 0; j < _cols; j++) {
            for (i = 0; i < _rows; i++) {
                trace += _matrix[j][i];
            }
        }
        return trace;
    }else{
        std::cerr << "irtkMatrix::Trace() matrix number of col != row" << std::endl;
        return 0;
    }
}

inline double irtkMatrix::Norm(void) const
{
  int i, j;
  double norm = 0;

  // The norm of a matrix M is defined as trace(M * M~)
  for (j = 0; j < _cols; j++) {
    for (i = 0; i < _rows; i++) {
      norm += _matrix[j][i]*_matrix[j][i];
    }
  }
  return std::sqrt(norm);
}

inline double irtkMatrix::InfinityNorm(void) const
{
  int i, j;
  double normInf = -1.0 * DBL_MAX;
  double sum;

  for (i = 0; i < _rows; ++i) {
    sum = 0;
    for (j = 0; j < _cols; ++j) {
      sum += abs(_matrix[j][i]);
    }
    if (sum > normInf)
      normInf = sum;
  }
  return normInf;
}

#endif

max_stars_repo_name: gordon-n-stevenson/fetalReconstruction

max_stars_repo_licenses: Zlib Unlicense Intel MIT

max_issues_repo_name: gordon-n-stevenson/fetalReconstruction

max_issues_repo_licenses: Zlib Unlicense Intel MIT

max_forks_repo_name: gordon-n-stevenson/fetalReconstruction

max_forks_repo_licenses: Zlib Unlicense Intel MIT


text:
/*
 * Copyright 2014 Google Inc. All rights reserved.
 *
 * Licensed under the Apache License, Version 2.0 (the "License");
 * you may not use this file except in compliance with the License.
 * You may obtain a copy of the License at
 *
 *    http://www.apache.org/licenses/LICENSE-2.0
 *
 * Unless required by applicable law or agreed to in writing, software
 * distributed under the License is distributed on an "AS IS" BASIS,
 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 * See the License for the specific language governing permissions and
 * limitations under the License.
 *
 * ---------------------------------------------------------------------------
 * Copyright 2014 Nervana Systems Inc.  All rights reserved.
 *
 * * Added sumabs, maxabs, save function headers
 * * Fixed numpy deprecation warnings
 * ---------------------------------------------------------------------------
 */

#ifndef MATRIX_H_
#define MATRIX_H_

#include "matrix_funcs.h"
#include <sstream>
#ifdef NUMPY_INTERFACE
#define NPY_NO_DEPRECATED_API NPY_1_7_API_VERSION
#include <Python.h>
#include <arrayobject.h>
#endif
#include <limits>
#include <assert.h>
#include <stdio.h>
#include <string.h>
#include <math.h>
#include <vector>
#include <cuda_runtime.h>

extern "C" {
#include <cblas.h>
}

#ifdef DOUBLE_PRECISION
#define CBLAS_GEMM cblas_dgemm
#define CBLAS_SCAL cblas_dscal
#define CBLAS_AXPY cblas_daxpy
#else
#define CBLAS_GEMM cblas_sgemm
#define CBLAS_SCAL cblas_sscal
#define CBLAS_AXPY cblas_saxpy
#endif /* DOUBLE_PRECISION */

#define MTYPE_MAX numeric_limits<MTYPE>::max()

typedef long long int int64;

class Matrix {
private:
    MTYPE* _data;
    bool _ownsData;
    int64 _numRows, _numCols;
    int64 _numElements;
    CBLAS_TRANSPOSE _trans;

    void _init(MTYPE* data, int64 numRows, int64 numCols, bool transpose, bool ownsData);
    void _tileTo2(Matrix& target) const;
    void _copyAllTo(Matrix& target) const;
    MTYPE _sum_column(int64 col) const;
    MTYPE _sum_row(int64 row) const;
    MTYPE _aggregate(MTYPE(*agg_func)(MTYPE, MTYPE), MTYPE initialValue) const;
    void _aggregate(int64 axis, Matrix& target, MTYPE(*agg_func)(MTYPE, MTYPE), MTYPE initialValue) const;
    MTYPE _aggregateRow(int64 row, MTYPE(*agg_func)(MTYPE, MTYPE), MTYPE initialValue) const;
    MTYPE _aggregateCol(int64 row, MTYPE(*agg_func)(MTYPE, MTYPE), MTYPE initialValue) const;
    void _updateDims(int64 numRows, int64 numCols);
    void _applyLoop(MTYPE(*func)(MTYPE));
    void _applyLoop(MTYPE (*func)(MTYPE), Matrix& target);
    void _applyLoop2(const Matrix& a, MTYPE(*func)(MTYPE, MTYPE), Matrix& target) const;
    void _applyLoop2(const Matrix& a, MTYPE (*func)(MTYPE,MTYPE, MTYPE), MTYPE scalar, Matrix& target) const;
    void _applyLoop2(const Matrix& a, MTYPE (*func)(MTYPE,MTYPE, MTYPE, MTYPE), MTYPE scalar1, MTYPE scalar2, Matrix& target) const;
    void _applyLoopScalar(const MTYPE scalar, MTYPE(*func)(MTYPE, MTYPE), Matrix& target) const;
    void _checkBounds(int64 startRow, int64 endRow, int64 startCol, int64 endCol) const;
    void _divideByVector(const Matrix& vec, Matrix& target);
    inline int64 _getNumColsBackEnd() const {
        return _trans == CblasNoTrans ? _numCols : _numRows;
    }
public:
    enum FUNCTION {
        TANH, RECIPROCAL, SQUARE, ABS, EXP, LOG, ZERO, ONE, LOGISTIC1, LOGISTIC2, SIGN
    };
    Matrix();
    Matrix(int64 numRows, int64 numCols);
    Matrix(int64 numRows, int64 numCols, bool transpose);
#ifdef NUMPY_INTERFACE
    Matrix(PyArrayObject *src);
    void _initpy(PyArrayObject *src);
#endif
    Matrix(const Matrix &like);
    Matrix(MTYPE* data, int64 numRows, int64 numCols);
    Matrix(float *buffer, int64 numRows, int64 numCols, int nil);
    Matrix(MTYPE* data, int64 numRows, int64 numCols, bool transpose);
    ~Matrix();

    inline MTYPE& getCell(int64 i, int64 j) const {
        assert(i >= 0 && i < _numRows);
        assert(j >= 0 && j < _numCols);
        if (_trans == CblasTrans) {
            return _data[j * _numRows + i];
        }
        return _data[i * _numCols + j];
    }

    MTYPE& operator()(int64 i, int64 j) const {
        return getCell(i, j);
    }

    inline MTYPE* getData() const {
        return _data;
    }

    inline bool isView() const {
        return !_ownsData;
    }

    inline int64 getNumRows() const {
        return _numRows;
    }

    inline int64 getNumCols() const {
        return _numCols;
    }

    inline int64 getNumDataBytes() const {
        return _numElements * sizeof(MTYPE);
    }

    inline int64 getNumElements() const {
        return _numElements;
    }

    inline int64 getLeadingDim() const {
        return _trans == CblasTrans ? _numRows : _numCols;
    }

    inline int64 getFollowingDim() const {
        return _trans == CblasTrans ? _numCols : _numRows;
    }

    inline CBLAS_TRANSPOSE getBLASTrans() const {
        return _trans;
    }

    inline bool isSameDims(const Matrix& a) const {
        return a.getNumRows() == getNumRows() && a.getNumCols() == getNumCols();
    }

    inline bool isTrans() const {
        return _trans == CblasTrans;
    }

    /*
     * Only use if you know what you're doing!
     * Does not update any dimensions. Just flips the _trans flag.
     *
     * Use transpose() if you want to get the transpose of this matrix.
     */
    inline void setTrans(bool trans) {
        assert(isTrans() == trans || !isView());
        _trans = trans ? CblasTrans : CblasNoTrans;
    }

    void apply(FUNCTION f);
    void apply(Matrix::FUNCTION f, Matrix& target);
    void subtractFromScalar(MTYPE scalar);
    void subtractFromScalar(MTYPE scalar, Matrix &target) const;
    void biggerThanScalar(MTYPE scalar);
    void smallerThanScalar(MTYPE scalar);
    void equalsScalar(MTYPE scalar);
    void biggerThanScalar(MTYPE scalar, Matrix& target) const;
    void smallerThanScalar(MTYPE scalar, Matrix& target) const;
    void equalsScalar(MTYPE scalar, Matrix& target) const;
    void biggerThan(Matrix& a);
    void biggerThan(Matrix& a, Matrix& target) const;
    void smallerThan(Matrix& a);
    void smallerThan(Matrix& a, Matrix& target) const;
    void minWith(Matrix &a);
    void minWith(Matrix &a, Matrix &target) const;
    void maxWith(Matrix &a);
    void maxWith(Matrix &a, Matrix &target) const;
    void equals(Matrix& a);
    void equals(Matrix& a, Matrix& target) const;
    void notEquals(Matrix& a) ;
    void notEquals(Matrix& a, Matrix& target) const;
    void add(const Matrix &m);
    void add(const Matrix &m, MTYPE scale);
    void add(const Matrix &m, MTYPE scaleThis, MTYPE scaleM);
    void add(const Matrix &m, Matrix& target);
    void add(const Matrix &m, MTYPE scaleM, Matrix &target);
    void add(const Matrix &m, MTYPE scaleThis, MTYPE scaleM, Matrix &target);
    void subtract(const Matrix &m);
    void subtract(const Matrix &m, Matrix& target);
    void subtract(const Matrix &m, MTYPE scale);
    void subtract(const Matrix &m, MTYPE scale, Matrix& target);
    void addVector(const Matrix& vec, MTYPE scale);
    void addVector(const Matrix& vec, MTYPE scale, Matrix& target);
    void addVector(const Matrix& vec);
    void addVector(const Matrix& vec, Matrix& target);
    void addScalar(MTYPE scalar);
    void addScalar(MTYPE scalar, Matrix& target) const;
    void maxWithScalar(MTYPE scalar);
    void maxWithScalar(MTYPE scalar, Matrix &target) const;
    void minWithScalar(MTYPE scalar);
    void minWithScalar(MTYPE scalar, Matrix &target) const;
    void eltWiseMultByVector(const Matrix& vec);
    void eltWiseMultByVector(const Matrix& vec, Matrix& target);
    void eltWiseDivideByVector(const Matrix& vec);
    void eltWiseDivideByVector(const Matrix& vec, Matrix& target);
    void resize(int64 newNumRows, int64 newNumCols);
    void resize(const Matrix& like);
    Matrix& slice(int64 startRow, int64 endRow, int64 startCol, int64 endCol) const;
    void slice(int64 startRow, int64 endRow, int64 startCol, int64 endCol, Matrix &target) const;
    Matrix& sliceRows(int64 startRow, int64 endRow) const;
    void sliceRows(int64 startRow, int64 endRow, Matrix& target) const;
    Matrix& sliceCols(int64 startCol, int64 endCol) const;
    void sliceCols(int64 startCol, int64 endCol, Matrix& target) const;
    void rightMult(const Matrix &b, MTYPE scale);
    void rightMult(const Matrix &b, Matrix &target) const;
    void rightMult(const Matrix &b);
    void rightMult(const Matrix &b, MTYPE scaleAB, Matrix &target) const;
    void addProduct(const Matrix &a, const Matrix &b, MTYPE scaleAB, MTYPE scaleThis);
    void addProduct(const Matrix& a, const Matrix& b);
    void eltWiseMult(const Matrix& a);
    void eltWiseMult(const Matrix& a, Matrix& target) const;
    void eltWiseDivide(const Matrix& a);
    void eltWiseDivide(const Matrix& a, Matrix &target) const;
    Matrix& transpose() const;
    Matrix& transpose(bool hard) const;
    Matrix& tile(int64 timesY, int64 timesX) const;
    void tile(int64 timesY, int64 timesX, Matrix& target) const;
    void copy(Matrix &dest, int64 srcStartRow, int64 srcEndRow, int64 srcStartCol, int64 srcEndCol, int64 destStartRow, int64 destStartCol) const;
    Matrix& copy() const;
    void copy(Matrix& target) const;
    Matrix& sum(int64 axis) const;
    void sum(int64 axis, Matrix &target) const;
    MTYPE sum() const;
    MTYPE sumabs() const;
    MTYPE max() const;
    MTYPE maxabs() const;
    Matrix& max(int64 axis) const;
    void max(int64 axis, Matrix& target) const;
    MTYPE min() const;
    Matrix& min(int64 axis) const;
    void min(int64 axis, Matrix& target) const;
    MTYPE norm() const;
    MTYPE norm2() const;
    void scale(MTYPE scale);
    void scale(MTYPE alpha, Matrix& target);
    void reshape(int64 numRows, int64 numCols);
    Matrix& reshaped(int64 numRows, int64 numCols);
    void printShape(const char* name) const;
    bool hasNan() const;
    bool hasInf() const;

    void randomizeNormal(MTYPE mean, MTYPE stdev);
    void randomizeUniform();
    void randomizeNormal();
    void print() const;
    void print(int64 startRow,int64 rows, int64 startCol,int64 cols) const;
    void print(int64 rows, int64 cols) const;
};

typedef std::vector<Matrix*> MatrixV;

#endif /* MATRIX_H_ */

max_stars_repo_name: TimDettmers/dlearndb

max_issues_repo_name: TimDettmers/dlearndb

max_forks_repo_name: TimDettmers/dlearndb


text:
#include <gsl/span>

namespace engine
{
	template <class ElementType, std::ptrdiff_t Extent = gsl::dynamic_extent>
	using span = gsl::span<ElementType, Extent>;
}


text:
/*
* -----------------------------------------------------------------
*  isat_lib.h
*  In Situ Adaptive Tabulation Library
*  Version: 2.0
*  Last Update: Nov 3, 2019
*
*  Programmer: Americo Barbosa da Cunha Junior
*              americo.cunhajr@gmail.com
* -----------------------------------------------------------------
*  Copyright (c) 2010-2019, Americo Barbosa da Cunha Junior
*  All rights reserved.
* -----------------------------------------------------------------
*  This is the header file for ISAT_LIB module, a computational
*  library with In Situ Adaptive Tabulation (ISAT) algorithm
*  routines.
* -----------------------------------------------------------------
*/

#ifndef __ISAT_LIB_H__
#define __ISAT_LIB_H__

#include <gsl/gsl_blas.h>
#include <gsl/gsl_math.h>

#include "../include/bst_lib.h"

/*
*------------------------------------------------------------
*   structure for ISAT workspace
*
*   last update: Oct 9, 2019
*------------------------------------------------------------
*/

typedef struct isat_struct
{
    bst_node *root;         /* binary search tree root       */
    unsigned int lf;        /* # of leaves                   */
    unsigned int nd;        /* # of nodes                    */
    unsigned int add;       /* # of adds                     */
    unsigned int grw;       /* # of grows                    */
    unsigned int rtv;       /* # of retrieves                */
    unsigned int dev;       /* # of direct evaluations       */
    unsigned int hgt;       /* tree height                   */
    unsigned int maxleaves; /* max. # of leaves in the tree  */
    clock_t time_add;       /* CPU clock spent by additions  */
    clock_t time_grw;       /* CPU clock spent by growths    */
    clock_t time_rtv;       /* CPU clock spent by retrieves  */
    clock_t time_dev;       /* CPU clock spent by dir. eval. */

} isat_wrk;
/*------------------------------------------------------------*/

/*
*------------------------------------------------------------
* function prototypes
*------------------------------------------------------------
*/

isat_wrk *isat_alloc();

void isat_free(void **isat_bl);

int isat_bst_init(isat_wrk *isat);

void isat_statistics(isat_wrk *isat);

int isat_input(unsigned int *max_lf,
               double *etol);

void isat_eoa_mtrx(gsl_matrix *A,
                   double etol,
                   gsl_matrix *L);

double isat_lerror(gsl_vector *Rphi,
                   gsl_vector *Rphi0,
		           gsl_matrix *A,
                   gsl_vector *phi,
                   gsl_vector *phi0);

int isat4(isat_wrk *isat_mem,
          void *cvode_mem,
	      double etol,
          double t0,
          double delta_t,
          gsl_vector *phi,
          gsl_matrix *A,
          gsl_matrix *L,
          gsl_vector *Rphi);

#endif /* __ISAT_LIB_H__ */


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

/* Author:  G. Jungman */

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

#include "error.h"

#include "bessel.h"

/* i_{l+1}/i_l
 */
static
int
bessel_il_CF1(const int l, const double x, const double threshold, double * ratio)
{
  const int kmax = 2000;
  double tk   = 1.0;
  double sum  = 1.0;
  double rhok = 0.0;
  int k;

  for(k=1; k<=kmax; k++) {
    double ak = (x/(2.0*l+1.0+2.0*k)) * (x/(2.0*l+3.0+2.0*k));
    rhok = -ak*(1.0 + rhok)/(1.0 + ak*(1.0 + rhok));
    tk  *= rhok;
    sum += tk;
    if(fabs(tk/sum) < threshold) break;
  }

  *ratio = x/(2.0*l+3.0) * sum;

  if(k == kmax)
    GSL_ERROR ("error", GSL_EMAXITER);
  else
    return GSL_SUCCESS;
}

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

int gsl_sf_bessel_i0_scaled_e(const double x, gsl_sf_result * result)
{
  double ax = fabs(x);

  /* CHECK_POINTER(result) */

  if(x == 0.0) {
    result->val = 1.0;
    result->err = 0.0;
    return GSL_SUCCESS;
  }
  else if(ax < 0.2) {
    const double eax = exp(-ax);
    const double y = ax*ax;
    const double c1 = 1.0/6.0;
    const double c2 = 1.0/120.0;
    const double c3 = 1.0/5040.0;
    const double c4 = 1.0/362880.0;
    const double c5 = 1.0/39916800.0;
    const double sum = 1.0 + y*(c1 + y*(c2 + y*(c3 + y*(c4 + y*c5))));
    result->val = eax * sum;
    result->err = 2.0 * GSL_DBL_EPSILON * result->val;
  }
  else if(ax < -0.5*GSL_LOG_DBL_EPSILON) {
    result->val = (1.0 - exp(-2.0*ax))/(2.0*ax);
    result->err = 2.0 * GSL_DBL_EPSILON * result->val;
  }
  else {
    result->val = 1.0/(2.0*ax);
    result->err = 2.0 * GSL_DBL_EPSILON * result->val;
  }
  return GSL_SUCCESS;
}

int gsl_sf_bessel_i1_scaled_e(const double x, gsl_sf_result * result)
{
  double ax = fabs(x);

  /* CHECK_POINTER(result) */

  if(x == 0.0) {
    result->val = 0.0;
    result->err = 0.0;
    return GSL_SUCCESS;
  }
  else if(ax < 3.0*GSL_DBL_MIN) {
    UNDERFLOW_ERROR(result);
  }
  else if(ax < 0.25) {
    const double eax = exp(-ax);
    const double y  = x*x;
    const double c1 = 1.0/10.0;
    const double c2 = 1.0/280.0;
    const double c3 = 1.0/15120.0;
    const double c4 = 1.0/1330560.0;
    const double c5 = 1.0/172972800.0;
    const double sum = 1.0 + y*(c1 + y*(c2 + y*(c3 + y*(c4 + y*c5))));
    result->val = eax * x/3.0 * sum;
    result->err = 2.0 * GSL_DBL_EPSILON * fabs(result->val);
    return GSL_SUCCESS;
  }
  else {
    double ex = exp(-2.0*ax);
    result->val = 0.5 * (ax*(1.0+ex) - (1.0-ex)) / (ax*ax);
    result->err = 2.0 * GSL_DBL_EPSILON * fabs(result->val);
    if(x < 0.0) result->val = -result->val;
    return GSL_SUCCESS;
  }
}

int gsl_sf_bessel_i2_scaled_e(const double x, gsl_sf_result * result)
{
  double ax = fabs(x);

  /* CHECK_POINTER(result) */

  if(x == 0.0) {
    result->val = 0.0;
    result->err = 0.0;
    return GSL_SUCCESS;
  }
  else if(ax < 4.0*GSL_SQRT_DBL_MIN) {
    UNDERFLOW_ERROR(result);
  }
  else if(ax < 0.25) {
    const double y = x*x;
    const double c1 = 1.0/14.0;
    const double c2 = 1.0/504.0;
    const double c3 = 1.0/33264.0;
    const double c4 = 1.0/3459456.0;
    const double c5 = 1.0/518918400.0;
    const double sum = 1.0 + y*(c1 + y*(c2 + y*(c3 + y*(c4 + y*c5))));
    const double pre = exp(-ax) * x*x/15.0;
    result->val = pre * sum;
    result->err = 2.0 * GSL_DBL_EPSILON * fabs(result->val);
    return GSL_SUCCESS;
  }
  else {
    double ex = exp(-2.0*ax);
    double x2 = x*x;
    result->val = 0.5 * ((3.0+x2)*(1.0-ex) - 3.0*ax*(1.0+ex))/(ax*ax*ax);
    result->err = 2.0 * GSL_DBL_EPSILON * fabs(result->val);
    return GSL_SUCCESS;
  }
}

int gsl_sf_bessel_il_scaled_e(const int l, double x, gsl_sf_result * result)
{
  double sgn = 1.0;
  double ax  = fabs(x);

  if(x < 0.0) {
    /* i_l(-x) = (-1)^l i_l(x) */
    sgn = ( GSL_IS_ODD(l) ? -1.0 : 1.0 );
    x = -x;
  }

  if(l < 0) {
    DOMAIN_ERROR(result);
  }
  else if(x == 0.0) {
    result->val = ( l == 0 ? 1.0 : 0.0 );
    result->err = 0.0;
    return GSL_SUCCESS;
  }
  else if(l == 0) {
    gsl_sf_result il;
    int stat_il = gsl_sf_bessel_i0_scaled_e(x, &il);
    result->val = sgn * il.val;
    result->err = il.err;
    return stat_il;
  }
  else if(l == 1) {
    gsl_sf_result il;
    int stat_il = gsl_sf_bessel_i1_scaled_e(x, &il);
    result->val = sgn * il.val;
    result->err = il.err;
    return stat_il;
  }
  else if(l == 2) {
    gsl_sf_result il;
    int stat_il = gsl_sf_bessel_i2_scaled_e(x, &il);
    result->val = sgn * il.val;
    result->err = il.err;
    return stat_il;
  }
  else if(x*x < 10.0*(l+1.5)/M_E) {
    gsl_sf_result b;
    int stat = gsl_sf_bessel_IJ_taylor_e(l+0.5, x, 1, 50, GSL_DBL_EPSILON, &b);
    double pre   = exp(-ax) * sqrt((0.5*M_PI)/x);
    result->val  = sgn * pre * b.val;
    result->err  = pre * b.err;
    result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
    return stat;
  }
  else if(l < 150) {
    gsl_sf_result i0_scaled;
    int stat_i0  = gsl_sf_bessel_i0_scaled_e(ax, &i0_scaled);
    double rat;
    int stat_CF1 = bessel_il_CF1(l, ax, GSL_DBL_EPSILON, &rat);
    double iellp1 = rat * GSL_SQRT_DBL_MIN;
    double iell   = GSL_SQRT_DBL_MIN;
    double iellm1;
    int ell;
    for(ell = l; ell >= 1; ell--) {
      iellm1 = iellp1 + (2*ell + 1)/x * iell;
      iellp1 = iell;
      iell   = iellm1;
    }
    result->val  = sgn * i0_scaled.val * (GSL_SQRT_DBL_MIN / iell);
    result->err  = i0_scaled.err * (GSL_SQRT_DBL_MIN / iell);
    result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
    return GSL_ERROR_SELECT_2(stat_i0, stat_CF1);
  }
  else if(GSL_MIN(0.29/(l*l+1.0), 0.5/(l*l+1.0+x*x)) < 0.5*GSL_ROOT3_DBL_EPSILON) {
    int status = gsl_sf_bessel_Inu_scaled_asymp_unif_e(l + 0.5, x, result);
    double pre = sqrt((0.5*M_PI)/x);
    result->val *= sgn * pre;
    result->err *= pre;
    return status;
  }
  else {
    /* recurse down from safe values */
    double rt_term = sqrt((0.5*M_PI)/x);
    const int LMAX = 2 + (int) (1.2 / GSL_ROOT6_DBL_EPSILON);
    gsl_sf_result r_iellp1;
    gsl_sf_result r_iell;
    int stat_a1 = gsl_sf_bessel_Inu_scaled_asymp_unif_e(LMAX + 1 + 0.5, x, &r_iellp1);
    int stat_a2 = gsl_sf_bessel_Inu_scaled_asymp_unif_e(LMAX     + 0.5, x, &r_iell);
    double iellp1 = r_iellp1.val;
    double iell   = r_iell.val;
    double iellm1 = 0.0;
    int ell;
    iellp1 *= rt_term;
    iell   *= rt_term;
    for(ell = LMAX; ell >= l+1; ell--) {
      iellm1 = iellp1 + (2*ell + 1)/x * iell;
      iellp1 = iell;
      iell   = iellm1;
    }
    result->val  = sgn * iellm1;
    result->err  = fabs(result->val)*(GSL_DBL_EPSILON + fabs(r_iellp1.err/r_iellp1.val) + fabs(r_iell.err/r_iell.val));
    result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);

    return GSL_ERROR_SELECT_2(stat_a1, stat_a2);
  }
}

int gsl_sf_bessel_il_scaled_array(const int lmax, const double x, double * result_array)
{
  if(x == 0.0) {
    int ell;
    result_array[0] = 1.0;
    for (ell = lmax; ell >= 1; ell--) {
      result_array[ell] = 0.0;
    };
    return GSL_SUCCESS;
  } else {
    int ell;
    gsl_sf_result r_iellp1;
    gsl_sf_result r_iell;
    int stat_0 = gsl_sf_bessel_il_scaled_e(lmax+1, x, &r_iellp1);
    int stat_1 = gsl_sf_bessel_il_scaled_e(lmax,   x, &r_iell);
    double iellp1 = r_iellp1.val;
    double iell   = r_iell.val;
    double iellm1;
    result_array[lmax] = iell;
    for(ell = lmax; ell >= 1; ell--) {
      iellm1 = iellp1 + (2*ell + 1)/x * iell;
      iellp1 = iell;
      iell   = iellm1;
      result_array[ell-1] = iellm1;
    }
    return GSL_ERROR_SELECT_2(stat_0, stat_1);
  }
}

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

#include "eval.h"

double gsl_sf_bessel_i0_scaled(const double x)
{
  EVAL_RESULT(gsl_sf_bessel_i0_scaled_e(x, &result));
}

double gsl_sf_bessel_i1_scaled(const double x)
{
  EVAL_RESULT(gsl_sf_bessel_i1_scaled_e(x, &result));
}

double gsl_sf_bessel_i2_scaled(const double x)
{
  EVAL_RESULT(gsl_sf_bessel_i2_scaled_e(x, &result));
}

double gsl_sf_bessel_il_scaled(const int l, const double x)
{
  EVAL_RESULT(gsl_sf_bessel_il_scaled_e(l, x, &result));
}


text:
/*
 * BSD 3-Clause License
 *
 * Copyright (c) 2018, mtezych
 * All rights reserved.
 *
 * Redistribution and use in source and binary forms, with or without
 * modification, are permitted provided that the following conditions are met:
 *
 * * Redistributions of source code must retain the above copyright notice,
 *   this list of conditions and the following disclaimer.
 *
 * * Redistributions in binary form must reproduce the above copyright notice,
 *   this list of conditions and the following disclaimer in the documentation
 *   and/or other materials provided with the distribution.
 *
 * * Neither the name of the copyright holder nor the names of its
 *   contributors may be used to endorse or promote products derived from
 *   this software without specific prior written permission.
 *
 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
 * AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO,
 * THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
 * PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR
 * CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
 * EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
 * PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS;
 * OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY,
 * WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR
 * OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE,
 * EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
 */

#ifndef CL_COMMAND_QUEUE
#define CL_COMMAND_QUEUE

#ifdef __APPLE__
	#include <OpenCL/cl.h>
#else
	#include <CL/cl.h>
#endif

#include <cl/Kernel.h>
#include <cl/Event.h>
#include <cl/Memory.h>

#include <util/util.h>

#include <gsl/span>

#include <array>
#include <vector>
#include <cstddef>

namespace cl
{
	struct Context;
	struct Device;

	enum NDRange : cl_uint
	{
		Range1D = 1,
		Range2D = 2,
		Range3D = 3,
	};

	struct Range
	{
		std::size_t offset;
		std::size_t size;
	};

	enum class ExecMode : cl_command_queue_properties
	{
		   InOrder = 0,
		OutOfOrder = CL_QUEUE_OUT_OF_ORDER_EXEC_MODE_ENABLE,
	};

	enum class Profiling : cl_command_queue_properties
	{
		 Enable = CL_QUEUE_PROFILING_ENABLE,
		Disable = 0,
	};

	struct CommandQueue
	{
		cl_command_queue clCommandQueue;

		CommandQueue
		(
			const Context& context, const Device& device,
			ExecMode  execMode  = ExecMode::InOrder,
			Profiling profiling = Profiling::Disable
		);

		~CommandQueue ();

		CommandQueue (CommandQueue&& commandQueue);
		CommandQueue (const CommandQueue& commandQueue) = delete;

		CommandQueue& operator = (CommandQueue&& commandQueue);
		CommandQueue& operator = (const CommandQueue& commandQueue) = delete;

		struct Properties
		{
			ExecMode  execMode;
			Profiling profiling;
		};

		enum class Info : cl_command_queue_info
		{
			ReferenceCount = CL_QUEUE_REFERENCE_COUNT,
			Properties     = CL_QUEUE_PROPERTIES,
		};

		template <Info info>
		auto GetInfo() const
		{
			auto infoSize = std::size_t { 0 };
			auto result = clGetCommandQueueInfo
			(
				clCommandQueue, util::enum_cast(info),
				0, nullptr, &infoSize
			);
			assert(result == CL_SUCCESS);

			auto infoBytes = std::vector<std::byte>
			{
				infoSize, std::byte { 0x00 }
			};
			result = clGetCommandQueueInfo
			(
				clCommandQueue, util::enum_cast(info),
				infoBytes.size(), infoBytes.data(), nullptr
			);
			assert(result == CL_SUCCESS);

			return InfoResult<info>::FromBytes(infoBytes);
		}

		template <NDRange WorkDimension>
		Event EnqueueKernel
		(
			const Kernel& kernel,
			const std::array<size_t, WorkDimension>& globalWorkSize,
			const std::array<size_t, WorkDimension>&  localWorkSize,
			const std::vector<Event>& waitEvents
		)
		{
			static_assert(sizeof(Event) == sizeof(cl_event));

			const auto globalWorkOffset = std::array<size_t, WorkDimension> { };

			auto signalEvent = cl_event { };

			auto result = clEnqueueNDRangeKernel
			(
				clCommandQueue, kernel.clKernel,
				WorkDimension,
				globalWorkOffset.data(),
				globalWorkSize.data(), localWorkSize.data(),
				static_cast<cl_uint>(waitEvents.size()),
				reinterpret_cast<const cl_event*>(util::data_or_null(waitEvents)),
				&signalEvent
			);
			assert(result == CL_SUCCESS);

			return Event { signalEvent };
		}

		template <typename ValueType>
		Event EnqueueReadBuffer
		(
			const Memory<ValueType>& memory,
			const gsl::span<ValueType> buffer,
			const Range& range,
			const std::vector<Event>& waitEvents
		)
		{
			assert(buffer.size() >= static_cast<std::ptrdiff_t>(range.size));

			auto signalEvent = cl_event { };

			auto result = clEnqueueReadBuffer
			(
				clCommandQueue, memory.clMemory,
				cl_bool { CL_BLOCKING },
				range.offset, sizeof(ValueType) * range.size,
				buffer.data(),
				static_cast<cl_uint>(waitEvents.size()),
				reinterpret_cast<const cl_event*>(util::data_or_null(waitEvents)),
				&signalEvent
			);
			assert(result == CL_SUCCESS);

			return Event { signalEvent };
		}

	private:

		template <Info info>
		struct InfoResult;
	};

	template <>
	struct CommandQueue::InfoResult<CommandQueue::Info::ReferenceCount>
	{
		static cl_uint
		FromBytes (const std::vector<std::byte>& infoBytes);
	};

	template <>
	struct CommandQueue::InfoResult<CommandQueue::Info::Properties>
	{
		static CommandQueue::Properties
		FromBytes (const std::vector<std::byte>& infoBytes);
	};
}

#endif

max_stars_repo_name: mtezych/cpp

max_issues_repo_name: mtezych/cpp

max_forks_repo_name: mtezych/cpp