text:
#ifndef INCLUDED_TELEPATH_LAPACK_DRIVERS_LLS_H
#define INCLUDED_TELEPATH_LAPACK_DRIVERS_LLS_H

#include <lapacke.h>

namespace telepath{

    template< typename T >
    auto posv( int matrix_order, char uplo, std::size_t n, std::size_t nrhs,
        T* a, std::size_t lda, T* b, std::size_t ldb );
    template<>
    inline auto posv<float>( int matrix_order, char uplo, std::size_t n,
        std::size_t nrhs, float* a, std::size_t lda, float* b,
        std::size_t ldb ){
        return LAPACKE_sposv( matrix_order, uplo, n, nrhs, a, lda, b, ldb );
    }
    template<>
    inline auto posv<double>( int matrix_order, char uplo, std::size_t n,
        std::size_t nrhs, double* a, std::size_t lda, double* b,
        std::size_t ldb ){
        return LAPACKE_dposv( matrix_order, uplo, n, nrhs, a, lda, b, ldb );
    }
    template<>
    inline auto posv<std::complex<float> >( int matrix_order, char uplo,
        std::size_t n, std::size_t nrhs, std::complex<float>* a,
        std::size_t lda, std::complex<float>* b, std::size_t ldb ){
        return LAPACKE_cposv( matrix_order, uplo, n, nrhs, a, lda, b, ldb );
    }
    template<>
    inline auto posv<std::complex<double> >( int matrix_order, char uplo,
        std::size_t n, std::size_t nrhs, std::complex<double>* a,
        std::size_t lda, std::complex<double>* b, std::size_t ldb ){
        return LAPACKE_zposv( matrix_order, uplo, n, nrhs, a, lda, b, ldb );
    }
    inline auto posv( int matrix_order, char uplo, std::size_t n,
        std::size_t nrhs, double* a, std::size_t lda, double* b,
        std::size_t ldb, double* x, std::size_t ldx, lapack_int* iter ){
        return LAPACKE_dsposv( matrix_order, uplo, n, nrhs, a, lda, b, ldb,
            x, ldx, iter );
    }
    inline auto posv( int matrix_order, char uplo, std::size_t n,
        std::size_t nrhs, std::complex<double>* a, std::size_t lda,
        std::complex<double>* b, std::size_t ldb, std::complex<double>* x,
        std::size_t ldx, lapack_int* iter ){
        return LAPACKE_zcposv( matrix_order, uplo, n, nrhs, a, lda, b, ldb,
            x, ldx, iter );
    }

    template< typename T >
    auto gelsd( int matrix_order, std::size_t m, std::size_t n,
        std::size_t nrhs, T* a, std::size_t lda, T* b, std::size_t ldb,
        T* s, T rcond, lapack_int* rank );
    template<>
    inline auto gelsd<double>( int matrix_order, std::size_t m, std::size_t n,
        std::size_t nrhs, double* a, std::size_t lda, double* b,
        std::size_t ldb,double* s, double rcond, lapack_int* rank ){
        return LAPACKE_dgelsd( matrix_order, m, n, nrhs, a, lda, b, ldb,
            s, rcond, rank );
    }
    template<>
    inline auto gelsd<float>( int matrix_order, std::size_t m, std::size_t n,
        std::size_t nrhs, float* a, std::size_t lda, float* b,
        std::size_t ldb,float* s, float rcond, lapack_int* rank ){
        return LAPACKE_sgelsd( matrix_order, m, n, nrhs, a, lda, b, ldb,
            s, rcond, rank );
    }
    template<>
    inline auto gelsd<std::complex<double> >( int matrix_order, std::size_t m,
        std::size_t n, std::size_t nrhs, std::complex<double>* a,
        std::size_t lda, std::complex<double>* b, std::size_t ldb,
        std::complex<double>* s, std::complex<double> rcond, lapack_int* rank ){
        return LAPACKE_zgelsd( matrix_order, m, n, nrhs, a, lda, b, ldb,
            s, rcond, rank );
    }
    template<>
    inline auto gelsd<std::complex<float> >( int matrix_order, std::size_t m,
        std::size_t n, std::size_t nrhs, std::complex<float>* a,
        std::size_t lda, std::complex<float>* b, std::size_t ldb,
        std::complex<float>* s, std::complex<float> rcond, lapack_int* rank ){
        return LAPACKE_cgelsd( matrix_order, m, n, nrhs, a, lda, b, ldb,
            s, rcond, rank );
    }

} //namespace telepath

#endif

max_stars_repo_name: tbepler/telepath

max_issues_repo_name: tbepler/telepath

max_forks_repo_name: tbepler/telepath


text:
#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include <assert.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_integration.h>
#include <gsl/gsl_spline.h>
#include <gsl/gsl_sort.h>
#include <gsl/gsl_fit.h>

#include "cosmocalc.h"

static double linear_corrfunc_integ_funct(double k, void *p);

static double linear_corrfunc_integ_funct(double k, void *p)
{
  double r = ((double*)p)[0];
  double a = ((double*)p)[1];
  return linear_powspec(k,a)*k/r;
}

double linear_corrfunc_exact(double r, double a)
{
  double I0;
  double abserr,p[2];
  gsl_integration_workspace *workspace;
  gsl_function F;
  gsl_integration_qawo_table *wf;
  gsl_integration_workspace *cycle_workspace;

#define WORKSPACE_NUM 70
#define ABSERR 1e-3
  workspace = gsl_integration_workspace_alloc((size_t) WORKSPACE_NUM);
  cycle_workspace = gsl_integration_workspace_alloc((size_t) WORKSPACE_NUM);
  wf = gsl_integration_qawo_table_alloc(r,1.0,GSL_INTEG_SINE,(size_t) WORKSPACE_NUM);

  F.function = &linear_corrfunc_integ_funct;
  p[0] = r;
  p[1] = a;
  F.params = p;
  gsl_integration_qawf(&F,1e-8,ABSERR,(size_t) WORKSPACE_NUM,workspace,cycle_workspace,wf,&I0,&abserr);

  gsl_integration_qawo_table_free(wf);
  gsl_integration_workspace_free(cycle_workspace);
  gsl_integration_workspace_free(workspace);
#undef ABSERR
#undef WORKSPACE_NUM

  return I0/2.0/M_PI/M_PI;
}

double linear_corrfunc(double r, double a)
{
  static int initFlag = 1;
  static int currCosmoNum;
  static gsl_spline *cosmocalc_linear_corrfunc_spline = NULL;
  static gsl_interp_accel *cosmocalc_linear_corrfunc_acc = NULL;

  double linear_corrfunc_table[COSMOCALC_LINEAR_CORRFUNC_TABLE_LENGTH];
  double r_table[COSMOCALC_LINEAR_CORRFUNC_TABLE_LENGTH];
  long i;
  double gf;

  if(initFlag == 1 || currCosmoNum != cosmoData.cosmoNum)
    {
      initFlag = 0;
      currCosmoNum = cosmoData.cosmoNum;

      for(i=0;i<COSMOCALC_LINEAR_CORRFUNC_TABLE_LENGTH;++i)
	{
	  r_table[i] = log(CF_R_MAX/CF_R_MIN)/(COSMOCALC_LINEAR_CORRFUNC_TABLE_LENGTH-1.0)*((double) i) + log(CF_R_MIN);
	  linear_corrfunc_table[i] = linear_corrfunc_exact(exp(r_table[i]),1.0);
	}

      //init the spline and accelerators
      if(cosmocalc_linear_corrfunc_spline != NULL)
	gsl_spline_free(cosmocalc_linear_corrfunc_spline);
      cosmocalc_linear_corrfunc_spline = gsl_spline_alloc(GSL_SPLINE_TYPE,(size_t) (COSMOCALC_LINEAR_CORRFUNC_TABLE_LENGTH));
      gsl_spline_init(cosmocalc_linear_corrfunc_spline,r_table,linear_corrfunc_table,(size_t) (COSMOCALC_LINEAR_CORRFUNC_TABLE_LENGTH));
      if(cosmocalc_linear_corrfunc_acc != NULL)
	gsl_interp_accel_reset(cosmocalc_linear_corrfunc_acc);
      else
	cosmocalc_linear_corrfunc_acc = gsl_interp_accel_alloc();
    }

  gf = growth_function(a);

  if(r < CF_R_MIN)
    return 0.0;
  else if(r < CF_R_MAX)
    return gsl_spline_eval(cosmocalc_linear_corrfunc_spline,log(r),cosmocalc_linear_corrfunc_acc)*gf*gf;
  else
    return 0.0;
}

ext: c

max_stars_repo_name: beckermr/cosmocalc

max_issues_repo_name: beckermr/cosmocalc

max_forks_repo_name: beckermr/cosmocalc


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

#ifndef KATLA_POSIX_FILE_H
#define KATLA_POSIX_FILE_H

#include "katla/core/core.h"

#include <optional>
#include <gsl/span>

#include <fcntl.h>

namespace katla {

    class PosixFile {
    public:
        enum class OpenFlags : int {
            ReadOnly    = O_RDONLY,
            WriteOnly   = O_WRONLY,
            ReadWrite   = O_RDWR,
            Create      = O_CREAT,
            Exclusive   = O_EXCL,
            Truncate    = O_TRUNC,
            Append      = O_APPEND,
            NonBlocking = O_NONBLOCK,
            Sync        = O_SYNC,
            Async       = O_ASYNC,

#if defined(O_TMPFILE)
            TmpFile     = O_TMPFILE
#endif
        };

        PosixFile();
        ~PosixFile();

        outcome::result<void> create(std::string_view filePath, OpenFlags flags);
        outcome::result<void> open(std::string_view filePath, OpenFlags flags);
        outcome::result<void> close();

        outcome::result<ssize_t> read(gsl::span<std::byte> &buffer);
        outcome::result<ssize_t> write(gsl::span<std::byte> &buffer);

        static outcome::result<std::string> absolutePath(std::string path);

        outcome::result<size_t> size();

    private:
        int m_fd;
    };

    inline PosixFile::OpenFlags operator|(PosixFile::OpenFlags lhs, PosixFile::OpenFlags rhs) {
        using T = std::underlying_type_t<PosixFile::OpenFlags>;
        return static_cast<PosixFile::OpenFlags>(static_cast<T>(lhs) | static_cast<T>(rhs));
    }

    inline PosixFile::OpenFlags &operator|=(PosixFile::OpenFlags &lhs, PosixFile::OpenFlags rhs) {
        lhs = lhs | rhs;
        return lhs;
    }

}

#endif // KATLA_POSIX_FILE_H

ext: h

max_stars_repo_name: plok/katla

max_issues_repo_name: plok/katla

max_forks_repo_name: plok/katla


text:
/* ============================================================ *
 * halomodel.h							*
 * Martin Kilbinger 2006-2009					*
 * ============================================================ */

#ifndef __HALOMODEL_H
#define __HALOMODEL_H

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

#include <gsl/gsl_sf_erf.h>

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

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

#define hm_base     -1900
#define hm_hodtype   hm_base + 1
#define hm_Mmin      hm_base + 2
#define hm_pofk      hm_base + 3
#define hm_nfw       hm_base + 4
#define hm_par       hm_base + 5
#define hm_overflow  hm_base + 6
#define hm_io        hm_base + 7
#define hm_zbin      hm_base + 8
#define hm_alpha     hm_base + 9
#define hm_negative  hm_base + 10
#define hm_zmean_2h  hm_base + 11
#define hm_halo_bias hm_base + 12
#define hm_undef     hm_base + 13
#define hm_gsl_int   hm_base + 14

/* Ranges of interpolation tables */
#define k_max_HOD     3336.0

/* Present critical density [M_sol h^2 / Mpc^3] */
#define rho_c0  2.7754e11

/* Mass limits for integration over mass functions */
#define logMmin (3.0*log(10.0))
#define logMmax (16.0*log(10.0))

/* Number of steps for scale-factor-integration (redshift) */
#define Na_hm 20

/* Bit-coded power spectrum types */
typedef enum {pofk_undef=-1, pl=1, pnl=2, p1hdm=4, p2hdm=8, pthdm=16, p1hg=32, p2hg=64, pthg=128,
	      p1hgcs=256, p1hgss=512, pstellar=1024} pofk_t;

/* Halo mass function type */
typedef enum {ps, st, st2, j01} massfct_t;
#define smassfct_t(i) ( \
  i==ps  ? "ps" : \
  i==st  ? "st" : \
  i==st2 ? "st2" : \
  i==j01 ? "j01" : \
  "")
#define Nmassfct_t 4

/* Halo bias type */
typedef enum {halo_bias_sc, halo_bias_tinker05, halo_bias_tinker10} halo_bias_t;
#define shalo_bias_t(i) ( \
   i==halo_bias_sc       ? "halo_bias_sc" : \
   i==halo_bias_tinker05 ? "halo_bias_tinker05" : \
   i==halo_bias_tinker10 ? "halo_bias_tinker10" : \
   "")
#define Nhalo_bias_t 3

/* HOD (Halo occupation distribution) type */
#define Nhod_t 5
typedef enum {hod_none, hamana04, berwein02, berwein02_hexcl, leauthaud11} hod_t;
#define shod_t(i) (		\
 i==hod_none  ? "hod_none" :	\
 i==hamana04  ? "hamana04" :	\
 i==berwein02 ? "berwein02" :	\
 i==berwein02_hexcl ? "berwein02_hexcl" :\
 i==leauthaud11 ? "leauthaud11" :	\
 "")

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

double FFTLog_TMP;

typedef struct FFTLog_complex
{
  double re;
  double im;
  double amp;
  double arg;
}  FFTLog_complex;

typedef struct {
  int N;
  fftw_plan p_forward;
  fftw_plan p_backward;
  fftw_complex *an;
  fftw_complex *ak;
  fftw_complex *cm;
  fftw_complex *um;
  fftw_complex *cmum;
  double min;
  double max;
  double q;
  double mu;
  double kr;
} FFTLog_config;

typedef struct {

  cosmo *cosmo;
  redshift_t *redshift;
  double zmin, zmax;

  /* Dark matter halo profile                                                      */
  double c0;			/* concentration parameter                         */
  double alpha_NFW;		/* density slope                                   */
  double beta_NFW;		/* concentration slope as fct of mass              */
  massfct_t massfct;            /* halo mass function				   */
  halo_bias_t halo_bias;        /* Halo bias                                       */

  /* Mass function parameters (Sheth&Torman). Do not set manually, they are set    *
   * in set_massfct() according to enum massfct.				   */
  double nmz_a;                 /* Called q in CS02                                */
  double nmz_p;                 /* a=1, p=1/2 is Press-Schechter mass fct.         */

  /* HOD (halo occupation distribution) parameters				   */
  hod_t hod;                    /* HOD type                                        */
  double M1, M0, sigma_log_M;
  double M_min;
  double alpha;
  double pi_max;
  double eta;                   /* central galaxy proportion                       */

  /* galaxy-galaxy lensing and wp(rp) */
  double log10Mhalo;
  double coord_phys;

  /* For Leauthaud11 model                                                         */
  double beta,delta,gamma,Mstar0;
  double beta_sat,B_sat,beta_cut,B_cut;
  double x;                     /* any parameter to propagate if needed            */
  double Mstellar_min, Mstellar_max;
  double fcen1, fcen2;

  /* Precomputed stuff */
  double A;			/* Mass function normalisation                     */
  double Mstar;                 /* M_*(a=1.0)					   */
  interTable2D *Pthdm;
  interTable *xir;
  interTable *xi_dm;
  interTable2D *rhohat;
  splineTable* sigRsqr;
  double a_xir;

  /* FFTLOG flag - OBSOLETE */
  int FFTLog;

} cosmo_hm;

typedef struct {

  cosmo    *cosmo;
  cosmo_hm *model;
  double   a, r, k, ng, ngp, eps, c;
  double  logMlim, bias_fac, Mh, Mstellar, Mstellar_min, Mstellar_max;
  double  M, r_vir, *kk;
  error    **err;

  double   logrmin, logrmax, rp, xi;
  gsl_interp_accel *acc;
  gsl_spline *spline;

  int i, j, type, asymptotic, logintegrate;

  double   (*bias_func)(double, void *);
} cosmo_hm_params;

typedef struct gsl_int_params
{
  void *params;
  funcwithpars func;
  error **err;

} gsl_int_params;

typedef struct {
  double *z;
  double *fac;
  double *ypn; /* for spline interpolation */
  double zm;   /* average weighted redshift*/
  int nbins;
} nz_t;

cosmo_hm* init_parameters_hm(double OMEGAM, double OMEGADE, double W0_DE, double W1_DE,
			     double *W_POLY_DE, int N_POLY_DE,
			     double H100, double OMEGAB, double OMEGANUMASS,
			     double NEFFNUMASS, double NORM, double NSPEC,
			     int Nzbin, const int *Nnz, const nofz_t *nofz, double *par_nz,
			     double zmin, double zmax,
			     nonlinear_t NONLINEAR, transfer_t TRANSFER, growth_t GROWTH,
			     de_param_t DEPARAM, norm_t normmode,
			     double C0, double ALPHANFW, double BETANFW, massfct_t MASSFCT, halo_bias_t HALO_BIAS,
			     double M_min, double M1, double M0, double sigma_log_M, double alpha,
			     double Mstar0, double beta, double delta, double gamma, double B_cut, double B_sat,
			     double beta_cut, double beta_sat, double Mstellar_min, double Mstellar_max, double eta,
			     double fcen1, double fcen2,
			     hod_t HOD, double pi_max, error **err);

cosmo_hm* copy_parameters_hm_only(cosmo_hm* source, error **err);
cosmo_hm *copy_parameters_hm(cosmo_hm *source, error **err);
void read_cosmological_parameters_hm(cosmo_hm **model, FILE *F, error **err);
cosmo_hm *set_cosmological_parameters_to_default_hm(error **err);
void free_parameters_hm(cosmo_hm** model);

void set_massfct(massfct_t massfct, double *nmz_a, double *nmz_p, error **err);
void dump_param_only_hm(cosmo_hm* model, FILE *F);
void dump_param_hm(cosmo_hm* model, FILE *F, error **err);

double sm2_rtbis(double (*func)(double, void *, error **), double x1, double x2,
		 double xacc, void *param, error **err);

/* From nrcomplex.h,c */
#ifndef _DCOMPLEX_DECLARE_T_
typedef struct DCOMPLEX {double r,i;} dcomplex;
#define _DCOMPLEX_DECLARE_T_
#endif /* _DCOMPLEX_DECLARE_T_ */
dcomplex Complex(double re, double im);
dcomplex Cadd(dcomplex a, dcomplex b);
dcomplex Cmul(dcomplex a, dcomplex b);
dcomplex Cdiv(dcomplex a, dcomplex b);
dcomplex RCmul(double x, dcomplex a);
void sm2_cisi(double x, double *ci, double *si, error **err);

double delta_c(cosmo *model, double a, error **err);
double bis_Mstar(double logM, void *param, error **err);
double bis_Mstar_a(double logM, void *param, error **err);
double Mstar(cosmo_hm *model, error **err);
double Mstar_a(cosmo_hm *model, double a, error **err);
double concentration(cosmo_hm *model, double Mh, double a, error **err);
double Delta_vir(cosmo_hm *model, double a);
double dsigma_R_sqr_dR(cosmo_hm *model, double R, error **err);
double nufnu(cosmo_hm *model, double nu, int asymptotic, error **err);
double nufnu_j01(double x);
double sigma_R_sqr(cosmo_hm *model, double R, error **err);
double sigmasqr_M(cosmo_hm *model, double M, error **err);

double dsigma_m1_dlnM(cosmo_hm *model, double M, error **err);
double dnu_dlnM(cosmo_hm *model, double M, double a, error **err);
double dn_dlnM_lnM(double logM, void *intpar, error **err);
double dn_dlnM_uf(double M, cosmo_hm *model, double a, error **err);
double dn_dlnM(double M, void *intpar, error **err);

double r_vir(cosmo_hm *model, double M, double a, error **err);
double M_vir(cosmo_hm *model, double r_vir, double a, error **err);
double Delta_h(cosmo_hm *model, double a, error **err);
double rho_crit(cosmo_hm *model, double a, error **err);
double rho_crit_halo(cosmo_hm *model, double a, error **err);
double Omega_m_halo(cosmo_hm *model, double a, error **err);
double rho_halo(cosmo_hm *model, double r, double a, double Mh, double c, error **err);

double DeltaSigma_WB2000(cosmo_hm *model, double r, const double a, const double M,  double c,  double Delta, error **err);
double g_inf(double x, error **err);
double g_sup(double x, error **err);

double int_for_rhohat(double, void *, error **err);
double rhohat_halo(cosmo_hm *model, double k, double M, double a, double c, error **err);

double halo_bias(cosmo_hm *model, double M, double a, int k, error **err);
double bias(cosmo_hm *model, double M, double a, int k, error **err);
double bias_tinker(cosmo_hm *model, double M, double a, error **err);
double bias_tinker10(cosmo_hm *model, double M, double a, error **err);
double int_for_bias_norm(double logM, void *intpar, error **err);
double bias_norm(cosmo_hm *model, double a, error **err);

double int_for_M_ij(double, void *, error **);
double M_ij(cosmo_hm *model, int i, int j, double a, const double *k, error **err);
double P1h_dm(cosmo_hm *model, double a, double k, error **err);
double P2h_dm(cosmo_hm *model, double a, double k, error **err);

double xi_dm_NL_OBSOLETE(cosmo_hm *model, double a, double r, error **err);  // non-linear DM xi
double int_for_xi_dm_NL_OBSOLETE(double k, void *intpar,  error **err);  // non-linear DM xi

#define CHANGE(fct) int change_##fct(cosmo_hm*, cosmo_hm*)

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

double int_gsl(funcwithpars func,void *params, double a, double b, double eps, error **err);
double integrand_gsl(double x,void *p);

CHANGE(massfct);
CHANGE(massfct_params);
CHANGE(halo_bias);
CHANGE(sigma_R_sqr);
CHANGE(Mstar);
CHANGE(rhohat_halo);
CHANGE(Pth);

#undef CHANGE

#endif

max_stars_repo_name: danielgruen/ccv

max_issues_repo_name: danielgruen/ccv

max_forks_repo_name: danielgruen/ccv


text:
/*
 * Copyright 2017 RIS.
 * Authored by Ron Dahlgren <ron@red83.net>
 */
#include "greatest.h"

#include "entities/entities.h"
#include <stdlib.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_vector.h>

// Alloc and free
TEST test_entity_alloc_and_free()
{
  double l, r;
  gsl_matrix* identity = gsl_matrix_alloc(4,4);
  gsl_matrix_set_identity(identity);

  entity_info* new_entity = NULL;
  // A newly alloc'd entity should have an identity transform
  new_entity = entities.entity_alloc();
  ASSERT_FALSE(new_entity == NULL);
  ASSERT_FALSE(new_entity->transform == NULL);
  for (int i = 0; i < 4; i++) {
    for (int j = 0; j < 4; j++) {
      l = gsl_matrix_get(identity, i, j);
      r = new_entity->transform[i * 4 + j];
      ASSERT_EQ(l, r);
    }
  }
  // function pointers should all be set to null
  new_entity->initialize = NULL;
  new_entity->terminate = NULL;
  new_entity->pre_update = NULL;
  new_entity->post_update = NULL;

  entities.entity_free(new_entity);
  PASS();
}

int pre_update_calls;
int post_update_calls;

void pre_update(double d, entity_info* info) {
  pre_update_calls++;
}

void post_update(double d, entity_info* info) {
  post_update_calls++;
}

TEST test_update_transforms()
{
  pre_update_calls = 0;
  post_update_calls = 0;
  // initialize the entity manager
  entities.initialize();
  // Allocate some entities
  entity_info* alfred = entities.entity_alloc();
  // alfred has no callbacks
  entity_info* bryce = entities.entity_alloc();
  // bryce has pre
  bryce->pre_update = pre_update;
  entity_info* charles = entities.entity_alloc();
  // charles has post
  charles->post_update = post_update;
  entity_info* danny = entities.entity_alloc();
  // danny has both
  danny->pre_update = pre_update;
  danny->post_update = post_update;
  // Associate them with the entity manager
  entities.add_entity(alfred);
  entities.add_entity(bryce);
  entities.add_entity(charles);
  entities.add_entity(danny);
  // do an update pass
  entities.update(1);
  // call backs should have fired
  ASSERT_EQ(pre_update_calls, 2);
  ASSERT_EQ(post_update_calls, 2);
  // remove them from the entity manager
  entities.remove_entity(alfred);
  entities.remove_entity(bryce);
  entities.remove_entity(charles);
  entities.remove_entity(danny);
  // do an update pass
  entities.update(1);
  // No callbacks should have fired
  ASSERT_EQ(pre_update_calls, 2);
  ASSERT_EQ(post_update_calls, 2);
  // free everything
  entities.entity_free(alfred);
  entities.entity_free(bryce);
  entities.entity_free(charles);
  entities.entity_free(danny);
  entities.terminate();
  PASS();
}

static void
init_entities(void *data)
{
  entities.initialize();
}

static void
destroy_entities(void *data)
{
  entities.terminate();
}

SUITE(entity_test_suite) {
  SET_SETUP(init_entities, NULL);
  SET_TEARDOWN(init_entities, NULL);
  RUN_TEST(test_entity_alloc_and_free);
  RUN_TEST(test_update_transforms);
}

GREATEST_MAIN_DEFS();

int main(int argc, char** argv) {
  GREATEST_MAIN_BEGIN();
  RUN_SUITE(entity_test_suite);
  GREATEST_MAIN_END();
}


text:
#ifndef __INC_MULTIVARSUPPORT__
#define __INC_MULTIVARSUPPORT__

#include "multi_modelstruct.h"
#include "emulator_struct.h"
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_vector.h>

struct multi_modelstruct;
struct emulator_struct;

typedef struct multi_emulator {
	int nt;
	int nr;
	int nparams;
	int nmodel_points;
	int nregression_fns;
	int nthetas;
	multi_modelstruct *model;
	emulator_struct **emu_struct_array;

} multi_emulator;

multi_emulator *alloc_multi_emulator(multi_modelstruct *model);

void free_multi_emulator(multi_emulator *e);

void estimate_multi(multi_modelstruct *m, FILE* outfp);

void emulate_point_multi(multi_emulator *emu, gsl_vector *the_point,
												 gsl_vector *the_mean,
												 gsl_vector *the_variance);

void emulate_point_multi_covar(multi_emulator *emu, gsl_vector *the_point,
															 gsl_vector *the_mean, gsl_vector *the_variance,
															 gsl_matrix *the_covar);

void emulate_point_multi_pca(multi_emulator *emu, gsl_vector *the_point,
												 gsl_vector *the_mean,
												 gsl_vector *the_variance);

#endif

max_stars_repo_name: jackdawjackdaw/emulator

max_issues_repo_name: jackdawjackdaw/emulator

max_forks_repo_name: jackdawjackdaw/emulator


text:
/* Siconos-Numerics, Copyright INRIA 2005-2018.
 * Siconos is a program dedicated to modeling, simulation and control
 * of non smooth dynamical systems.
 * Siconos is a free software; you can redistribute it and/or modify
 * it under the terms of the GNU Lesser General Public License as published by
 * the Free Software Foundation; either version 2 of the License, or
 * (at your option) any later version.
 * Siconos 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 Lesser General Public License for more details.
 *
 * You should have received a copy of the GNU Lesser General Public License
 * along with Siconos; if not, write to the Free Software
 * Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301  USA
 *
 * Contact: Vincent ACARY, siconos-team@lists.gforge.inria.fr
 */

#ifndef SiconosBlas_H
#define SiconosBlas_H
#include "SiconosConfig.h"

#if defined(__cplusplus)
extern "C"
{
#endif

// tells include-what-you-use to keep this file
// and not to suggest cblas.h or alike.
// IWYU pragma: begin_exports
#if defined(HAS_MKL_CBLAS)
#include <mkl_cblas.h>
#elif defined(HAS_MATLAB_BLAS)
#include <blas.h>
#define cblas_daxpy daxpy
#define cblas_dcopy dcopy
#define cblas_ddot ddot
#define cblas_dgemm dgemm
#define cblas_dgemv dgemv
#define cblas_dnrm2 dnrm2
#define cblas_dscal dscal
#else
#include <cblas.h>
#endif
// IWYU pragma: end_exports

#ifdef __cplusplus
}
#undef restrict
#define restrict __restrict
#endif

static inline double* NMD_row_rmajor(double* restrict mat, unsigned ncols, unsigned rindx)
{
  return &mat[rindx*ncols];
}

static inline void NMD_copycol_rmajor(int nrows, double* col, double* restrict mat, int ncols, unsigned cindx)
{
  cblas_dcopy(nrows, col, 1, &mat[cindx], ncols);
}

static inline void NMD_dense_gemv(int nrows, int ncols, double alpha, double* restrict mat, double* restrict y, double beta, double* restrict x)
{
  cblas_dgemv(CblasColMajor, CblasTrans, ncols, nrows, alpha, mat, ncols, y, 1, beta, x, 1);
}

#endif // SiconosBlas_H

max_stars_repo_name: fperignon/sandbox

max_issues_repo_name: fperignon/sandbox

max_forks_repo_name: fperignon/sandbox


text:
#ifndef __FUNCTIONLIB_H__
#define __FUNCTIONLIB_H__

#include <iostream>
#include <math.h>
#include <stdio.h>
#include <string.h>
#include <gsl/gsl_sf.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_linalg.h>
#include <gsl/gsl_permutation.h>
#include <gsl/gsl_complex.h>
#include <gsl/gsl_complex_math.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_eigen.h>
#include <gsl/gsl_poly.h>
#include <gsl/gsl_integration.h>
#include <gsl/gsl_sort.h>
#include <gsl/gsl_sort_vector.h>

#endif

///////////////////////// Definition of the structures, mainly used as parameters for functions ///////////////////////

///////////Meaning of Words:////////////////////
//
//  coef ---> coefficient
//  transl ---> tranlation
//  stddev ---> standard deviation
//  magnif ---> magnification of the variable,
//  e.g.: sin(2x) where 2 is the magnification
//
//  in all functions the combination of magnification
//  and translation is interpreted as
//  f( magnif( x - transl) ).
//
////////////////////////////////////////////////

//struct of chained list describing a compound function of several terms multiplying together, where ID refers to the type of funcion (in other words the 'true' struct the pointer 'target' is referring to) . 'flag' signifies this term is 'multiplying' (value 1) or 'dividing' (other value) the whole term.

typedef struct MulMap
{
    void * NEXT;
    void * target;
    int ID;
    int flag;
}MulMap;

//struct of chained list describing a series of multiplied compound function summed up together, where flag means the term is summing (value 1) or subtracting (other value) to the whole term (though coef in each struct describing the functions has the ability to describe this)

typedef struct SumMap
{
    void *NEXT;
    void *MulMaphead;
    int flag;
}SumMap;

typedef struct intresult
{
    double result;
    double error;
}intresult;

//ID: 1
typedef struct PolyFunc
{
    int length;
    double* coef;
    double transl;
    double magnif;
}PolyFunc;

//ID: 2
typedef struct GaussianFunc
{
    double mean;
    double stddev;
    double coef;
}GaussianFunc;

//ID: 3
typedef struct SinCosFunc
{
    double* sincoef;
    double* coscoef;
    double transl;
    double magnif;
    int length;
}SinCosFunc;

//ID: 4
typedef struct Exp
{
    double coef;
    double transl;
    double magnif;
}Exp;

//ID: 5
typedef struct Log
{
    double coef;
    double magnif;
    double base;
    double transl;
}Log;

//ID: 6
typedef struct LegendreP
{
    int l;
    double coef;
    double magnif;
    double transl;
}LegendreP;

//ID: 7
typedef struct Rect
{
    double a;
    double b;

    double height;
}Rect;

double IDMapping(double x, int ID, void * param);
double FunctionMultiply(double x, void * HeadMulMap);
double FunctionSum(double x,void * summap);

double Combination(int n,int r);

double PolyFunction(double x, void * param);
double GaussianFunction(double x,void * param);
double SinCosFunction(double x,void * param);
double Exponential(double x,void * param);
double Logarithm(double x,void * param);

//Series of functions to help create the structs
void FuncStructConstruct(FILE * inputfile,void * target, int ID);
void MulMapConstruct(FILE * inputfile, MulMap * headpointer);
void SumMapConstruct(FILE * inputfile, SumMap * headpointer);

//Series of functions to help free the memory taken by the structs
void FuncStructDelete(void* target, int ID);
void MulMapDelete(void * target);
void SumMapDelete(void * target);

void OrthoPolyCoefGen(gsl_matrix * target,int rows,int columns,double a,double b);
void PolynomialTranslationCoef(gsl_vector * target,PolyFunc * param,int vectorlength);
void LegendrePCoefficient(PolyFunc * target, LegendreP * param);

//GSL trifles
void gsl_vector_complex_convert(gsl_vector * source, gsl_vector_complex * target, int length);
void gsl_matrix_complex_convert(gsl_matrix * source, gsl_matrix_complex * target, int rows, int columns);
void gsl_vector_poly_convert(PolyFunc * target,gsl_vector * source, int vectorlength);
void gsl_vector_complex_extract(gsl_vector_complex * source, gsl_vector * real, gsl_vector * imag, int length);
void gsl_matrix_complex_extract(gsl_vector_complex * source, gsl_matrix * real,gsl_matrix * imag, int rows, int columns);
void gsl_vector_complex_combine(gsl_vector * real, gsl_vector * imag, gsl_vector_complex * target);
void gsl_matrix_complex_combine(gsl_matrix * real,gsl_matrix * imag, gsl_matrix_complex * target);
void gsl_matrix_diag(gsl_matrix * target, gsl_vector * diag, int length);
void gsl_matrix_complex_diag(gsl_matrix_complex * target, gsl_vector_complex * diag, int length);
void gsl_matrix_mul(gsl_matrix * A, gsl_matrix *B, gsl_matrix * Result,int Acolumn,int Arow,int Bcolumn);
void gsl_matrix_complex_mul(gsl_matrix_complex * A, gsl_matrix_complex *B, gsl_matrix_complex * Result,int Acolumn,int Arow,int Bcolumn);
double gsl_vector_inner_product(gsl_vector * A, gsl_vector * B,int length);
gsl_complex gsl_vector_complex_product(gsl_vector_complex * A, gsl_vector_complex * B, int length);
gsl_complex gsl_vector_complex_inner_product(gsl_vector_complex * A, gsl_vector_complex * B, int length);
void gsl_vector_transform(gsl_vector * vec,gsl_matrix * trf,int length);
void gsl_vector_complex_transform(gsl_vector_complex * vec, gsl_matrix_complex * trf,int length);
void gsl_matrix_unitmatrix(gsl_matrix * m,int length);
void gsl_matrix_complex_unitmatrix(gsl_matrix_complex * m,int length);
void gsl_vector_complex_conjugate(gsl_vector_complex * v, int length);
void gsl_matrix_complex_conjugate(gsl_matrix_complex * m, int rows, int columns);

void gsl_vector_poly_convert(PolyFunc * target,gsl_vector * source, int vectorlength);

max_stars_repo_name: Walter-Feng/HomebrewLib

max_issues_repo_name: Walter-Feng/HomebrewLib

max_forks_repo_name: Walter-Feng/HomebrewLib


text:
/*
 * Simultaneous ML fitting of N 3D Gaussian kernels against a volumetric image
 * MATLAB interface in df_mlfitN.c
 *
 * TODO:
 * - how to handle background, still as a constant over the
 * patch?
 * - bg not handled yet.
 *
 */

#include <assert.h>
#include <inttypes.h>
#include <stdint.h>
#include <stdio.h>
#include <string.h>

#include <gsl/gsl_multimin.h>
#include <gsl/gsl_statistics.h>

#include "blit3.h"
#include "mlfit.h"

#ifndef verbose
#define verbose 1
#endif

// When low on bugs, uncomment the following lines for speed
// Try with valgrind and cmocka first ...
//
// replace gsl_vector_get(v,i) by v->data[i*v->stride]
// #define GSL_RANGE_CHECK_OFF

// Globals
uint32_t maxIterations = 5000;
double convCriteria = 1e-6;

// Headers
int    localizeDot(double *, size_t, double *,  double *);
int    localize(double *, size_t, size_t, size_t, double *, size_t, size_t, double *);
double matrix_get(double * , size_t, size_t, size_t, double *);
double my_f (const gsl_vector *, void *);
int    unit_tests(void);

// Optimization constants
typedef struct {
  double * V; // Volumetric image
  double * W; // Temporary space for model during iterations
  double * D; // dots, x,y,z, nphot, sigmax, sigmay, sigmaz
  size_t Vm, Vn, Vp;
  size_t Dm, Dn;
  double * bg;
} optParams;

double my_f (const gsl_vector *v, void *params)
  // The function to optimize. mlsimplex controls/varies the variables
  // in v. params are to set up the context
{

  optParams *p = (optParams *) params;

  // Size of volumetric image
  const size_t Vm = p->Vm; // Size of image
  const size_t Vn = p->Vn;
  const size_t Vp = p->Vp;
  const size_t Dn = p->Dn; // number of dots
  const size_t Dm = p->Dm;  // and number of features
  double * D = p->D;
  const double * restrict V = p->V; // Volumetric data to fit against
  double * W = p->W; // Pre allocated memory for the model
  double * bg = p->bg;

#if verbose > 0
  printf("-> my_f, %d dots\n", (int) Dn);
#endif

  // Get the other parameters ...

  double photsum = 0;

  // Copy the cordinates and nphot at the current iteration to D
  for(size_t dd = 0; dd<Dn; dd++)
  {
#if verbose > 0
    printf("%lu: (%.1f %.1f %.1f) %.0f (%.1f %.1f %.1f)\n", dd, D[Dm*dd+0],D[Dm*dd+1], D[Dm*dd+2], D[Dm*dd+3], D[Dm*dd+4], D[Dm*dd+5], D[Dm*dd+6]);
#endif
    D[Dm*dd+0] = gsl_vector_get(v, 4*dd+0); // x
    D[Dm*dd+1] = gsl_vector_get(v, 4*dd+1); // y
    D[Dm*dd+2] = gsl_vector_get(v, 4*dd+2); // z
    D[Dm*dd+3] = gsl_vector_get(v, 4*dd+3); // nphot

#if verbose>0
    printf("%lu: xyz (%.1f %.1f %.1f) %.0f (%.1f %.1f %.1f)\n", dd, D[Dm*dd+0],D[Dm*dd+1], D[Dm*dd+2], D[Dm*dd+3], D[Dm*dd+4], D[Dm*dd+5], D[Dm*dd+6]);
#endif
    if(D[Dm*dd+3]<0)
      return INFINITY;

    photsum = photsum + D[7*dd+3];
  }

#if verbose > 0
  printf("phosum of dot parameters: %f\n", photsum);
#endif

  // Copy background map to model
  memcpy(W, bg, sizeof(double)*Vm*Vn*Vp);
  //  memset(W, 0, sizeof(double)*Vm*Vn*Vp); // W = 0;
  blit3g(W, Vm, Vn, Vp, D, Dn, mid_point,0);

  // set bg as well with blit3(g)
  for(size_t kk = 0; kk<Vn*Vm*Vp; kk++) {
    //W[kk]++;
  }

  double vmin = 10e99;
  double vmax = -10e99;
  double wmin = 10e99;
  double wmax = -10e99;
  double diffmax = 0;
  for(size_t kk = 0; kk<Vn*Vm*Vp; kk++)
  {
    vmin = GSL_MIN(vmin, V[kk]);
    vmax = GSL_MAX(vmax, V[kk]);
    wmin = GSL_MIN(wmin, W[kk]);
    wmax = GSL_MAX(wmax, W[kk]);
    diffmax = GSL_MAX(diffmax, fabs(V[kk]-W[kk]));
  }
#if verbose>0
  printf("V -- min: %f max %f\n", vmin, vmax);
  printf("W -- min: %f max %f\n", wmin, wmax);
  printf("diffmax: %f\n", diffmax);
#endif

  /*
   * from LL2PG.m
   *   model = x(3)+x(4)*gaussianInt2([x(1), x(2)], x(5), (size(patch, 1)-1)/2);
   *   mask = disk2d((size(patch,1)-1)/2);
   *   %L = -sum(sum(-(patch-model).^2./model - .5*log(model)));
   *   L = -sum(sum(mask.*(-(patch-model).^2./model - .5*log(model))));
   */

  double E = 0;
  for (size_t kk=0; kk<Vm*Vn*Vp; kk++) {
    //E+= (V[kk]-W[kk])*(V[kk]-W[kk])/V[kk] - .5*log(W[kk]);
    E+= (V[kk]-W[kk])*(V[kk]-W[kk]);
  }
  //  E = -E;

  /* Quadratic
   * for (size_t kk=0; kk<Rw*Rw; kk++)
   * E+= (GI[kk]-R[kk])*(GI[kk]-R[kk]);
   */
#if verbose >0
  printf("E: %f\n", E);
#endif
  return E;
}
double matrix_get(double * V, size_t Vm, size_t Vn, size_t Vp, double * D)
  // Safely get V(D), i.e. check bounds. Returns 0 if outside
  //
{

  double x = D[0];
  double y = D[1];
  double z = D[2];

  if(x<0)
    return 0;
  if(x>Vm-1)
    return 0;
  if(y<0)
    return 0;
  if(y>Vn-1)
    return 0;
  if(z<0)
    return 0;
  if(z>Vp-1)
    return 1;

  size_t pos = (size_t) nearbyint(x) + nearbyint(y)*Vm +  nearbyint(z)*Vm*Vn;
  assert(pos<Vm*Vn*Vp);
  return V[pos];
}

int localize(double * V, size_t Vm, size_t Vn, size_t Vp,
    double * D, size_t Dm, size_t Dn, double * F)
  // Localization for a dot roughly centered in V of size Vm x Vm
  // D[0], D[1], D[3] are the global coordinates of the dot
  // F are the fitted coordinates
{

  // avoid crashing matlab
  gsl_set_error_handler_off();

#if verbose > 0
  printf("-> Localize\n");
  printf("V: %lu %lu %lu\n", Vm, Vn, Vp);
  printf("D: %lu %lu\n", Dm, Dn);
#endif
  assert(Dm == 7); // x y z nphot sigmax sigmay sigmaz

  // Non-optimized parameters
  optParams par;
  par.V = V;
  par.Vm = Vm;
  par.Vn = Vn;
  par.Vp = Vp;
  par.Dn = Dn;
  par.Dm = Dm;
  par.W = malloc(Vm*Vn*Vp*sizeof(double));
  assert(par.W != NULL);
  par.D = D;
  par.bg = malloc(Vm*Vn*Vp*sizeof(double));

#if verbose > 0
  printf("Estimating local background\n");
#endif
  // Estimating a constant from the surroundings of the first dot
  // printf("Vm: %lu %Vn: %lu Vp: %lu, D: %f %f %f\n", Vm, Vn, Vp, D[0], D[1], D[2]);
  double bg = estimateBGV(V, Vm, Vn, Vp, D);
#if verbose>0
  printf("Background estimated to %f\n", bg);
#endif
  printf("Background estimated to %f\n", bg);

  for(size_t kk=0; kk<Vm*Vn*Vp; kk++)
    par.bg[kk] = bg;

  const gsl_multimin_fminimizer_type *T =
    gsl_multimin_fminimizer_nmsimplex2;
  gsl_multimin_fminimizer *s = NULL;
  gsl_vector *ss, *x;
  gsl_multimin_function minex_func;

  size_t iter = 0;
  int status;
  double size;

  /* Starting point */
#if verbose > 0
  printf("Setting up starting vector\n");
#endif
  const int nf = 4; // # of features to be optimized: x, y, z, nphot
  int nParameters = nf*Dn;
  x = gsl_vector_alloc (nParameters);
  for(size_t dd = 0; dd<Dn; dd++)
  {
    gsl_vector_set(x, dd*nf+0, D[dd*Dm+0]); // x position
    gsl_vector_set(x, dd*nf+1, D[dd*Dm+1]); // y position
    gsl_vector_set(x, dd*nf+2, D[dd*Dm+2]); // z position

    // number of photons
    // Since the blitting is scaled by the central pixel, the number
    // of photons can be estimated as the central value of each peak -
    // background from the surrounding
    double nphot = gsl_max(0, matrix_get(V, Vm, Vn, Vp, D+ dd*Dm)-par.bg[0]);
    printf("nphot: %f\n", nphot);
    gsl_vector_set(x, dd*nf+3, nphot);
  }
#if verbose > 0
  printf("Setting up initial step sizes\n");
#endif
  /* Set initial step sizes */
  ss = gsl_vector_alloc(nParameters);
  gsl_vector_set_all(ss, 0.01);
  for(size_t dd = 0; dd<Dn; dd++)
  {
    gsl_vector_set(ss, dd*nf+3, 5); // Number of photons
  }

#if verbose > 0
  printf("Testing the fitting with the initial data\n");
#endif
  my_f(x,(void*) &par);

  /* Initialize method and iterate */
  minex_func.n = nParameters;
  minex_func.f = my_f;
  minex_func.params = ∥
#if verbose > 0
  printf("Allocating for the solver\n");
#endif
  s = gsl_multimin_fminimizer_alloc (T, nParameters);
#if verbose > 0
  printf("Initializing solver\n");
#endif
  gsl_multimin_fminimizer_set (s, &minex_func, x, ss);

  do
  {
    iter++;
#if verbose>0
    printf("----------------------------- Iteration: %lu\n", iter);
#endif
    status = gsl_multimin_fminimizer_iterate(s);

    if (status)
      break;

    size = gsl_multimin_fminimizer_size(s);
    status = gsl_multimin_test_size(size, convCriteria);

    if (status == GSL_SUCCESS)
    {
#if verbose > 0
      printf ("converged to minimum at\n");
      printf ("%5lu x:%10.3e y:%10.3e z:%10.3e NP:%6.1f f() = %7.3f size = %10.3e\n",
          iter,
          gsl_vector_get (s->x, 0),
          gsl_vector_get (s->x, 1),
          gsl_vector_get (s->x, 2),
          gsl_vector_get (s->x, 3),
          s->fval, size);
#endif
    }
  }
  while (status == GSL_CONTINUE && iter < maxIterations);

  for(size_t dd = 0; dd<Dn; dd++)
  {
    F[5*dd+0] = gsl_vector_get (s->x, nf*dd+0); // x
    F[5*dd+1] = gsl_vector_get (s->x, nf*dd+1); // y
    F[5*dd+2] = gsl_vector_get (s->x, nf*dd+2); // z
    F[5*dd+3] = gsl_vector_get (s->x, nf*dd+3); // nphot
    F[5*dd+4] = status;                            // status
  }

  gsl_vector_free(x);
  gsl_vector_free(ss);
  gsl_multimin_fminimizer_free (s);

  free(par.W);
  free(par.bg);
  return status;
}

#ifndef _MAIN
#define _MAIN
int unit_tests(){
  double * V; // image
  int Vm = 102; int Vn = 107; int Vp = 60;
  double * D; // list of dots
  size_t Dm = 7; // number of features of the dots
  size_t Dn = 10; // number of dots
  double * F; // fitted dots

  printf("Image size: %dx%dx%d\n", Vm, Vn, Vp);
  printf("Localizing %lu dots\n", Dn);

  V = malloc(Vm*Vn*Vp*sizeof(double));
  D = malloc(Dm*Dn*sizeof(double));
  F = malloc(5*Dn*sizeof(double));

  // Initialize the data
  for(int kk=0; kk<Vm*Vn*Vp; kk++)
    V[kk] = rand_range(2,9);

  for(uint32_t kk=0; kk<Dn; kk++)
  {
    size_t pos = kk*Dm;
    D[pos] =   rand_range(6, Vm-7);
    D[pos+1] = rand_range(6, Vn-7);
    D[pos+2] = rand_range(6, Vp-7);
    D[pos+3] = 1000;
    D[pos+4] = 1.2;
    D[pos+5] = 1.2;
    D[pos+6] = 1.4;

    if(D[pos] < 6)
      D[pos] = 6;
    if(D[pos+1] < 6)
      D[pos+1] = 6;
    if(D[pos+2] < 6)
      D[pos+2] = 6;
#if verbose > 0
    printf("D %03d %f %f %f\n", kk, D[pos], D[pos+1], D[pos+2]);
#endif
  }

  for(uint32_t kk=0; kk<Dn; kk++)
  {
    size_t pos = nearbyint(D[kk*7]) + nearbyint(D[kk*7+1])*Vm + nearbyint(D[kk*7+2])*Vm*Vn;
    assert(pos<(size_t) Vm*Vn*Vp);
    assert(pos>0);
    V[pos] = 5;
  }

  // Run the optimization
  localize(V, Vm, Vn, Vp, D, Dm, Dn, F);

  // In next version, also supply clustering information

#if verbose >0
  for(size_t kk = 0; kk<Dm; kk++)
  {
    size_t Dpos = kk*7;
    size_t Fpos = kk*5;
    printf("%6lu (%f, %f, %f) -> (%f, %f, %f)\n", kk,
        D[Dpos], D[Dpos+1], D[Dpos+2],
        F[Fpos], F[Fpos+1], F[Fpos+2]);
  }
#endif
  free(F);
  free(D);
  free(V);
  return 0;
}

int main(int argc, char ** argv)
  // For testing, not used any more, see the MATLAB interface in
  // df_mlfit.c
{

  printf("%s\n", argv[0]);

  if(argc == 1)
    return unit_tests();

  return 0;
}
#endif

max_stars_repo_name: elgw/dotter

max_issues_repo_name: elgw/dotter

max_forks_repo_name: elgw/dotter


text:
#include <gsl/gsl_math.h>
#include <gsl/gsl_cblas.h>
#include "cblas.h"

void
cblas_zgeru (const enum CBLAS_ORDER order, const int M, const int N,
             const void *alpha, const void *X, const int incX, const void *Y,
             const int incY, void *A, const int lda)
{
#define BASE double
#include "source_geru.h"
#undef BASE
}

max_stars_repo_name: manggoguy/parsec-modified

max_issues_repo_name: manggoguy/parsec-modified

max_forks_repo_name: manggoguy/parsec-modified


text:
/* -*- Mode: c++ -*- */
// copyright (c) 2006-2010 by Christos Dimitrakakis
// <christos.dimitrakakis@gmail.com>
/***************************************************************************
 *                                                                         *
 *   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.                                   *
 *                                                                         *
 ***************************************************************************/

#ifndef MATRIX_H
#define MATRIX_H

#include <gsl/gsl_cblas.h>
#include <gsl/gsl_linalg.h>
#include <gsl/gsl_matrix.h>
#include <exception>
#include <iostream>
#include <stdexcept>
#include <vector>
#include "Vector.h"
#include "real.h"
//#include <atlas/cblas.h>

#undef REFERENCE_ACCESS

#define ACCURACY_LIMIT 1e-12

/** \brief An n-by-m dimensional matrix.

    TODO: Use the BLAS interface for some / most routines.
 */
class Matrix {
 public:
  /// Check bounds?
  enum BoundsCheckingStatus { NO_CHECK_BOUNDS = 0, CHECK_BOUNDS = 1 };
#ifdef NDEBUG
  static Matrix Null(int rows, int columns,
                     enum BoundsCheckingStatus check = NO_CHECK_BOUNDS);
  static Matrix Unity(int rows, int columns,
                      enum BoundsCheckingStatus check = NO_CHECK_BOUNDS);

  Matrix(int rows = 1, int columns = 1,
         enum BoundsCheckingStatus check = NO_CHECK_BOUNDS);
  Matrix(int rows, int columns, real* y,
         enum BoundsCheckingStatus check = NO_CHECK_BOUNDS);
  explicit Matrix(const Vector& v,
                  enum BoundsCheckingStatus check_ = NO_CHECK_BOUNDS);
#else
  static Matrix Null(int rows, int columns,
                     enum BoundsCheckingStatus check = CHECK_BOUNDS);
  static Matrix Unity(int rows, int columns,
                      enum BoundsCheckingStatus check = CHECK_BOUNDS);

  Matrix(int rows = 1, int columns = 1,
         enum BoundsCheckingStatus check = CHECK_BOUNDS);
  Matrix(int rows, int columns, real* y,
         enum BoundsCheckingStatus check = CHECK_BOUNDS);
  explicit Matrix(const Vector& v,
                  enum BoundsCheckingStatus check_ = CHECK_BOUNDS);
#endif
  Matrix(const Matrix& rhs, bool clone = true);
  ~Matrix();
  void Resize(int rows_, int columns_);
  Matrix AddRow(const Vector& rhs);
  Matrix AddColumn(const Vector& rhs);
  Matrix& operator=(const Matrix& rhs);
  bool operator==(const Matrix& rhs) const;
  bool operator!=(const Matrix& rhs) const;
  Matrix operator+(const Matrix& rhs);
  Matrix& operator+=(const Matrix& rhs);
  Matrix operator-(const Matrix& rhs);
  Matrix& operator-=(const Matrix& rhs);
  Matrix& operator*=(const real& rhs);
  Matrix operator*(const Matrix& rhs);
  Matrix operator*(const real& rhs);
  Matrix operator+(const real& rhs);
  Matrix operator-(const real& rhs);
  Matrix operator/(const real& rhs);
  /// Matrix inversion (defaults to GSL with LU)
  Matrix Inverse(real epsilon = ACCURACY_LIMIT) const {
    return GSL_Inverse();
    // return Inverse_LU();
  }
  Vector SVD_Solve(const Vector& b) const;

  Matrix GSL_Inverse() const;

  /** Matrix inversion using the Cholesky decomposition.

      First call Cholesky with accuracy epsilon to obtain $\fL,
      L'\f$ such that \f$A = LL'\f$.  Then solve for \f$X = A^{-1}\f$:
      \f[
      LL'X = I,
      \f]
      by dynamic programming.
  */
  Matrix Inverse_Cholesky(real epsilon = ACCURACY_LIMIT) const {
    Matrix U = Cholesky(epsilon);
    Matrix L(U, false);
    L.Transpose();
    return Inverse(L, U);
  }
  /** Matrix inversion using the LU decomposition.

      First call LUDecomposition with accuracy epsilon to obtain $\fL,
      U\f$ such that \f$A = LU\f$.  Then solve for \f$X = A^{-1}\f$:
      \f[
      LUX = I,
      \f]
      by dynamic programming.
  */
  Matrix Inverse_LU(real epsilon = ACCURACY_LIMIT) const {
    real det;
    Matrix tmp(*this);
    std::vector<Matrix> A = tmp.LUDecomposition(det, epsilon);
    return Inverse(A[0], A[1]);
  }

  Matrix Inverse(const Matrix& L, const Matrix& U) const;
  bool isSymmetric() const;
  bool isTriangular() const;
  bool isUpperTriangular() const;
  bool isLowerTriangular() const;
  real det() const;
  real tr() const;
  real L1Norm() const;
  real L2Norm() const;
  Matrix Multiple(const Matrix& rhs) const;
  Matrix Kron(const Matrix& rhs) const;
  real ColumnSum(int c) const;
  real RowSum(int r) const;
  real compute_det_triangular() const;
  real gaussian_elimination_forward(real epsilon = ACCURACY_LIMIT);
  std::vector<Matrix> QRDecomposition() const;
  std::vector<Matrix> LUDecomposition(real& determinant,
                                      real epsilon = ACCURACY_LIMIT);
  Matrix Cholesky(real epsilon = ACCURACY_LIMIT) const;
  void Cholesky(Matrix& chol, real epsilon = ACCURACY_LIMIT) const;
  void Clear();
  void Transpose();
  Vector getColumn(int c) const;
  Vector getRow(int r) const;
  void setColumn(int c, const Vector& x);
  void setRow(int r, const Vector& x);
  void SortRow(int r);
  void SortColumn(int r);
  int Rows() const;
  int Columns() const;
  Vector Vec() const;
  void Vec(const Vector& x);
  Vector RowMax() const;
  Vector ColumnMax() const;
  Vector RowMin() const;
  Vector ColumnMin() const;
  real& operator()(int i, int j);
  const real& operator()(int i, int j) const;
  void print(FILE* f) const;
  friend Matrix operator*(const real& lhs, const Matrix& rhs);
  friend Matrix operator*(const Vector& lhs, const Matrix& rhs);
  friend Vector operator*(const Matrix& lhs, const Vector& rhs);

 protected:
  int rows;     ///< number of rows in the matrix
  int columns;  ///< number of columns in the matrix
  real* x;      ///< data
#ifdef REFERENCE_ACCESS
  real** x_list;  ///< data pointers
  void MakeReferences();
#endif
  const enum BoundsCheckingStatus checking_bounds;
  bool transposed;
  bool clear_data;
  const real& qGet(int i, int j);
  void qSet(int i, int j, real v);
};

Matrix operator*(const real& lhs, const Matrix& rhs);
Matrix operator*(const Vector& lhs, const Matrix& rhs);
Vector operator*(const Matrix& lhs, const Vector& rhs);
real Mahalanobis2(const Vector& x, const Matrix& A, const Vector& y);

/// Kronecker product
inline Matrix Kron(const Matrix& lhs, const Matrix& rhs) {
  return lhs.Kron(rhs);
}

inline Matrix OuterProduct(const Vector& lhs, const Vector& rhs) {
  Matrix R(lhs.n, rhs.n);
  MatrixProduct(lhs, rhs, R);
  return R;
}

Matrix Transpose(const Matrix& rhs);

/// In-place transpose matrix
inline void Matrix::Transpose() { transposed = !transposed; }

inline int Matrix::Rows() const {
  if (transposed) {
    return columns;
  }
  return rows;
}

inline int Matrix::Columns() const {
  if (transposed) {
    return rows;
  }
  return columns;
}

inline real& Matrix::operator()(int i, int j) {
  if (transposed) {
    int tmp = i;
    i = j;
    j = tmp;
  }

  if (checking_bounds) {
    if (i < 0 || j < 0 || i >= rows || j >= columns) {
      fprintf(stderr, "bad access: (%d, %d) on a %d x %d matrix\n", i, j, rows,
              columns);
      fflush(stderr);
      throw std::out_of_range("matrix index out of range");
    }
  }
#ifdef REFERENCE_ACCESS
  return x_list[i][j];
#else
  return x[i * columns + j];
#endif
}

inline const real& Matrix::operator()(int i, int j) const {
  if (transposed) {
    int tmp = i;
    i = j;
    j = tmp;
  }

  if (checking_bounds) {
    if (i < 0 || j < 0 || i >= rows || j >= columns) {
      fprintf(stderr, "bad access: (%d, %d) on a %d x %d matrix\n", i, j, rows,
              columns);
      fflush(stderr);
      throw std::out_of_range("matrix index out of range");
    }
  }
#ifdef REFERENCE_ACCESS
  return x_list[i][j];
#else
  return x[i * columns + j];
#endif
}

#if 0
inline const real& Matrix::qGet(int i, int j) {
#ifdef REFERENCE_ACCESS
  return x_list[i][j];
#else
  return x[i * columns + j];
#endif
}

inline void Matrix::qSet(int i, int j, real v) {
#ifdef REFERENCE_ACCESS
  x_list[i][j] = v;
#else
  x[i * columns + j] = v;
#endif
}
#endif

#endif

max_stars_repo_name: litlpoet/beliefbox

max_issues_repo_name: litlpoet/beliefbox

max_forks_repo_name: litlpoet/beliefbox


text:
#pragma once

#include <type_traits>
#include <vector>

#include <gsl-lite/gsl-lite.hpp>

namespace thrustshift {

/*! \brief Make multiple spans from a range of objects which contain ranges.
 *
 *  \param r
 *  \param get_range This functor must **not** return a copy of any object
 *
 *  ```cpp
 *  struct S {
 *
 *    std::vector<int> v;
 *  };
 *  ...
 *
 *  std::vector<S> ss;
 *  // fill ss
 *
 *  //// WRONG
 *  auto spans_with_invalid_ptrs = make_spans_from_ranges(ss, [](const auto& e) { return e.v; });
 *
 *  //// CORRECT
 *  auto spans = make_spans_from_ranges(ss, [](const auto& e) { return gsl_lite::make_span(e.v); });
 *  ```
 *
 *  The wrong lambda is returning a copy of member vector `v`, which is then deallocated.
 *
 */
template <class T, class RangeOfRanges, class GetRangeCallback>
std::vector<gsl_lite::span<T>> make_spans_from_ranges(
    RangeOfRanges&& r,
    GetRangeCallback&& get_range) {
	std::vector<gsl_lite::span<T>> spans(r.size());
	for (size_t i = 0; i < spans.size(); ++i) {
		spans[i] =
		    gsl_lite::make_span(get_range(r[i]).data(), get_range(r[i]).size());
	}
	return spans;
}

template <class T, class RangeOfRanges>
std::vector<gsl_lite::span<T>> make_spans_from_ranges(RangeOfRanges&& r) {
	auto get_range = [](auto& c) { return gsl_lite::make_span(c); };
	return make_spans_from_ranges<T>(std::forward<RangeOfRanges>(r), get_range);
}

template <class T, class RangeOfPtrs>
std::vector<gsl_lite::span<T>> make_spans_from_ptrs(RangeOfPtrs&& r, size_t N) {

	std::vector<gsl_lite::span<T>> spans(r.size());
	for (size_t i = 0; i < spans.size(); ++i) {
		spans[i] = gsl_lite::make_span(r[i].get(), N);
	}
	return spans;
}

template <typename T>
gsl_lite::span<T> subtract(gsl_lite::span<T>& pool, size_t size_of_new_span) {
	gsl_Expects(pool.size() >= size_of_new_span);
	auto piece = pool.first(size_of_new_span);
	pool = pool.subspan(size_of_new_span);
	return piece;
}

} // namespace thrustshift

max_stars_repo_name: pauleonix/thrustshift

max_issues_repo_name: pauleonix/thrustshift

max_forks_repo_name: pauleonix/thrustshift


text:
/******************************************\
* OpenMP/SSE/BLAS Optimized Matrix Library *
*                                          *
* by                                       *
*   Elliott Forney                         *
*   3.9.2010                               *
\******************************************/

/*
 *  Libraries
 */

#include <stdio.h>
#include <stdlib.h>
#include <time.h>
#include <math.h>

// cblas
#include <cblas.h>

// OpenMP
#include <omp.h>

// SSE 4.1 intrinsics
#include <smmintrin.h>

// SSE 2 intrinsics
// #include <pmmintrin.h>

//
#include "errcheck.h"
#include "matrix.h"

/*
 *  Macros
 */

//** Debugging level
#define DEBUG 0

//** Initialization

// max number of chars per line in input table
#define line_buff_size 1024

//** Access

// row major indexing
// expects stride and data to be set
#define data_at(r,c) data[r*stride+c]

//** Addition/Subtraction

// when number of values in matrix is
// larger than add_lim use big kernel
#define add_lim 10000

//** Transpose

#define trans_lim 22500

//#define trans_tile 64

/*
 *  Global variables
 */

unsigned rand_state = 8;

/*
 *  Function prototypes
 */

//
void build_drc(matrix *m);

// sigmoid function
float phi(const float v);

// derivative of sigmoid function given phi(v)
float phi_prime(const float z);

// addition kernel for small matrices
void add_small(float *a, float *b, float *c, const unsigned n);

// addition kernel for big matrices
void add_big(float *m, float *b, float *c, const unsigned n);

// subtraction kernel for small matrices
void sub_small(float *a, float *b, float *c, const unsigned n);

// subtraction kernel for big matrices
void sub_big(float *a, float *b, float *c, const unsigned n);

// pointwise square kernel for small matrices
void sqr_small(float *a, float *b, const unsigned n);

// pointwise square kernel for big matrices
void sqr_big(float *a, float *b, const unsigned n);

// scalar multiply kernel for small matrices
void scl_small(float *a, float b, float *c, const unsigned n);

// scalar multiply kernel for big matrices
void scl_big(float *a, float b, float *c, const unsigned n);

// scalar multiplication & pointwise addition kernel for small matrices
void scl_add_small(float *a, float b, float *c, float *d, const unsigned n);

// scalar multiplication & pointwise addition kernel for big matrices
void scl_add_big(float *a, float b, float *c, float *d, const unsigned n);

// pointwise multiplication for small matrices
void pmult_small(float *a, float *b, float *c, const unsigned n);

// pointwise multiplication for big matrices
void pmult_big(float *a, float *b, float *c, const unsigned n);

/*
 *  Function bodies
 */

// build row-col indices for double index
void build_drc(matrix *m)
{
  // not done/needed yet
}

// sigmoid function
float phi(const float v)
{
  // recommended by Lecun, find citation!!
  const float scl = 2.0f/3.0f;
  return 1.7159f * tanh( scl * v );
}

// derivative of sigmoid function given phi(v)
float phi_prime(const float z)
{
  const float scl = 2.0f/3.0f;
  return scl * (1.7159f - z*z);
}

// addition kernel for small matrices
void add_small(float *a, float *b, float *c, const unsigned n)
{
  // pointer to last destination
  const float *a_end = a+n;

  // loop through each value in destination
  while (a < a_end)
  {
    // four values from b and c into sse registers
    __m128 mm_v1 = _mm_load_ps(b);
    __m128 mm_v2 = _mm_load_ps(c);

    // add our sse vectors
    mm_v1 = _mm_add_ps(mm_v1, mm_v2);

    // store result into a
    _mm_store_ps(a, mm_v1);

    // increment pointers by 4
    a += 4; b += 4; c += 4;
  }
}

// addition kernel for big matrices
void add_big(float *a, float *b, float *c, const unsigned n)
{
  unsigned i;

  // loop through a, b and c by 4 in parallel
  #pragma omp parallel for
  for (i = 0; i < n; i += 4)
  {
    // load four values into sse registers
    __m128 mm_v1 = _mm_load_ps(b+i); // what if they have different strides!!! :(
    __m128 mm_v2 = _mm_load_ps(c+i);

    // add our sse vectors
    mm_v1 = _mm_add_ps(mm_v1, mm_v2);

    // store result into a
    _mm_store_ps(a+i, mm_v1);
  }
}

// subtraction kernel for small matrices
void sub_small(float *a, float *b, float *c, const unsigned n)
{
  const float *a_end = a+n;

  while (a < a_end)
  {
    __m128 mm_v1 = _mm_load_ps(b);
    __m128 mm_v2 = _mm_load_ps(c);

    mm_v1 = _mm_sub_ps(mm_v1, mm_v2);

    _mm_store_ps(a, mm_v1);

    a += 4; b += 4; c += 4;
  }
}

// subtraction kernel for big matrices
void sub_big(float *a, float *b, float *c, const unsigned n)
{
  unsigned i;

  #pragma omp parallel for
  for (i = 0; i < n; i += 4)
  {
    __m128 mm_v1 = _mm_load_ps(b+i);
    __m128 mm_v2 = _mm_load_ps(c+i);

    mm_v1 = _mm_sub_ps(mm_v1, mm_v2);

    _mm_store_ps(a+i, mm_v1);
  }
}

// pointwise square kernel for small matrices
void sqr_small(float *a, float *b, const unsigned n)
{
  const float *a_end = a+n;

  while (a < a_end)
  {
    __m128 mm_v = _mm_load_ps(b);

    mm_v = _mm_mul_ps(mm_v, mm_v);

    _mm_store_ps(a, mm_v);

    a += 4; b += 4;
  }
}

// pointwise square kernel for big matrices
void sqr_big(float *a, float *b, const unsigned n)
{
  unsigned i;

  #pragma omp parallel for
  for (i = 0; i < n; i += 4)
  {
    __m128 mm_v = _mm_load_ps(b+i);

    mm_v = _mm_mul_ps(mm_v, mm_v);

    _mm_store_ps(a+i, mm_v);
  }
}

// scalar multiplication kernel for small matrices
void scl_small(float *a, float b, float *c, const unsigned n)
{
  const float *a_end = a+n;

  __m128 mm_s = _mm_set1_ps(b);

  while (a < a_end)
  {
    __m128 mm_v = _mm_load_ps(c);

    mm_v = _mm_mul_ps(mm_s, mm_v);

    _mm_store_ps(a, mm_v);

    a += 4; c += 4;
  }
}

// scalar multiplication kernel for big matrices
void scl_big(float *a, float b, float *c, const unsigned n)
{
  unsigned i;

  const __m128 mm_s = _mm_set1_ps(b);

  #pragma omp parallel for
  for (i = 0; i < n; i += 4)
  {
    __m128 mm_v = _mm_load_ps(c+i);

    mm_v = _mm_mul_ps(mm_s, mm_v);

    _mm_store_ps(a+i, mm_v);
  }
}

// scalar multiplication & pointwise addition kernel for small matrices
void scl_add_small(float *a, float b, float *c, float *d, unsigned n)
{
  const float *a_end = a+n;

  __m128 mm_s = _mm_set1_ps(b);

  while (a < a_end)
  {
    __m128 mm_v1 = _mm_load_ps(c);
    __m128 mm_v2 = _mm_load_ps(d);

    mm_v1 = _mm_mul_ps(mm_s,  mm_v1);
    mm_v1 = _mm_add_ps(mm_v1, mm_v2);

    _mm_store_ps(a, mm_v1);

    a += 4; c += 4; d += 4;
  }
}

// scalar multiplication & pointwise addition kernel for big matrices
void scl_add_big(float *a, float b, float *c, float *d, unsigned n)
{
  unsigned i;

  const __m128 mm_s = _mm_set1_ps(b);

  #pragma omp parallel for
  for (i = 0; i < n; i += 4)
  {
    __m128 mm_v1 = _mm_load_ps(c+i);
    __m128 mm_v2 = _mm_load_ps(d+i);

    mm_v1 = _mm_mul_ps(mm_s,  mm_v1);
    mm_v1 = _mm_add_ps(mm_v1, mm_v2);

    _mm_store_ps(a+i, mm_v1);
  }
}

// pointwise multiplication kernel for small matrices
void pmult_small(float *a, float *b, float *c, const unsigned n)
{
  // pointer to last destination
  const float *a_end = a+n;

  // loop through each value in destination
  while (a < a_end)
  {
    // four values from b and c into sse registers
    __m128 mm_v1 = _mm_load_ps(b);
    __m128 mm_v2 = _mm_load_ps(c);

    // mult our sse vectors
    mm_v1 = _mm_mul_ps(mm_v1, mm_v2);

    // store result into a
    _mm_store_ps(a, mm_v1);

    // increment pointers by 4
    a += 4; b += 4; c += 4;
  }
}

// pointwise multiplication kernel for big matrices
void pmult_big(float *a, float *b, float *c, const unsigned n)
{
  unsigned i;

  // loop through a, b and c by 4 in parallel
  #pragma omp parallel for
  for (i = 0; i < n; i += 4)
  {
    // load four values into sse registers
    __m128 mm_v1 = _mm_load_ps(b+i);
    __m128 mm_v2 = _mm_load_ps(c+i);

    // mult our sse vectors
    mm_v1 = _mm_mul_ps(mm_v1, mm_v2);

    // store result into a
    _mm_store_ps(a+i, mm_v1);
  }
}

/*
 *  External function bodies
 */

//** Initialization & Destruction

// initialize matrix m of size r by c
void matrix_init(matrix *m, const unsigned r, const unsigned c)
{
  // set matrix dimensions
  m->r = r;
  m->c = c;

  // zero pad to nearest mult of 64, lcm 4 for sse, 64 for trans
  // leave one extra row/col to append bias terms
//  m->rstride = r + (64 - (r % 64)); // lcm(trans_tile, 4)
//  m->cstride = c + (64 - (c % 64));

  // 64 tiles speeds up transpose but hurts everything else
  m->rstride = r + (4 - (r % 4));
  m->cstride = c + (4 - (c % 4));

  // allocate space for matrix, 16 byte aligned
  m->cpu_data = (float*)_mm_malloc((m->rstride)*(m->cstride)*
                                   sizeof(float), 16);
  if (m->cpu_data == NULL)
    errcheck_cpu();

  // allocate space for last column holder
  m->cv = (float*)malloc((m->r)*sizeof(float));
  if (m->cv == NULL)
    errcheck_cpu();

  // zero out data for padding
  memset(m->cpu_data, '\0', (m->rstride)*(m->cstride)*sizeof(float));
  memset(m->cv, '\0', (m->r)*sizeof(float));

  // build row-col indices
  build_drc(m);
}

// initialize matrix m from ascii table
void matrix_init_file(matrix *m, char *file_name)
{
  // will need to do this dynamically if
  // we ever want to handle big tables!!!
  float table_data[line_buff_size][line_buff_size];
  unsigned r = 0, c = 0;

  // open file & setup input buffer
  FILE *table = fopen(file_name, "r");
  char  line_buffer[line_buff_size];

  // check if we were even able to open table file
  if (file_name == NULL)
    errcheck_cpu();

  // buckle-up, don't drink and code
  memset(table_data, '\0', line_buff_size*line_buff_size*sizeof(float));

  // for each line in table (row)
  while (fgets(line_buffer, line_buff_size, table))
  {
    if (DEBUG > 5)
      printf("line buffer: %s\n", line_buffer);

    // set up string tokenizer on input buffer
    char *cur_val = strtok(line_buffer, " "); // Note reentrant or thread safe!!

    // don't increment num rows on blank line
    if (cur_val != NULL)
    {
      c = 0; // new row, reset col counter

      // for each token (col)
      while (cur_val != NULL)
      {
        // convert from char to float
        table_data[r][c] = atof(cur_val);

        if (DEBUG > 5)
          printf("converting %d %d %f\n", r, c, table_data[r][c]);

        // get next token
        cur_val = strtok(NULL, " ");
        ++c; // increment num cols
      }
      ++r; // increment num rows
    }
  }

  // close file descriptor
  fclose(table);

  // initialize m
  matrix_init(m, r, c);

  // setup for data_at
  unsigned mr, mc;
  float *data = m->cpu_data;
  const unsigned stride = m->cstride;

  // loop through collected data and put into m
  for (mr = 0; mr < r; ++mr)
    for (mc = 0; mc < c; ++mc)
      data_at(mr,mc) = table_data[mr][mc];
}

// load zeros into matrix m
void matrix_load_zero(matrix m)
{
  float *data = m.cpu_data;
  const float *data_end = data + (m.rstride*m.cstride);

  // loop through all values and set to zero
  // memset?
  while (data < data_end)
    *(data++) = 0.0f;
}

// load values from an array
void matrix_load_array(matrix m, float *v)
{
  unsigned r, c, i = 0;

  // set data and stride for data_at
  float *data = m.cpu_data;
  const unsigned stride = m.cstride;

  for (r = 0; r < m.r; ++r)
    for (c = 0; c < m.c; ++c)
      data_at(r,c) = v[i++];
}

// load values from the random uniform distribution
void matrix_load_runif(matrix m, float min, float max)
{
  unsigned r, c;

  // set data and stride for data_at
  float *data = m.cpu_data;
  const unsigned stride = m.cstride;

  const float range = max - min;

  // needs work to be multithreaded!

  // set each non-padding value from random uniform
  for (r = 0; r < m.r; ++r)
    for (c = 0; c < m.c; ++c)
      data_at(r,c) = rand_r(&rand_state) * range / RAND_MAX + min;
}

void matrix_load_testa(matrix m, unsigned n)
{
  unsigned r, c;

  // set data and stride for data_at
  float *data = m.cpu_data;
  const unsigned stride = m.cstride;

  #pragma omp parallel for private(c)
  for (r = 0; r < m.r; ++r)
    for (c = 0; c < m.c; ++c)
      //data_at(r,c) = (float)n * (float)r*c;
      data_at(r,c) = (float)n;
}

void matrix_load_testb(matrix m)
{
  unsigned r, c;

  // set data and stride for data_at
  float *data = m.cpu_data;
  const unsigned stride = m.cstride;

  #pragma omp parallel for private(c)
  for (r = 0; r < m.r; ++r)
    for (c = 0; c < m.c; ++c)
      //data_at(r,c) = (float)r;
      data_at(r,c) = 1.0f;
}

void matrix_load_testc(matrix m)
{
  unsigned r, c;

  // set data and stride for data_at
  float *data = m.cpu_data;
  const unsigned stride = m.cstride;

  #pragma omp parallel for private(c)
  for (r = 0; r < m.r; ++r)
    for (c = 0; c < m.c; ++c)
      //data_at(r,c) = (float)c;
      data_at(r,c) = 1.0f;
}

// copy values to an array
void matrix_unload_array(matrix m, float *v)
{
  unsigned r, c, i = 0;

  float *data = m.cpu_data;
  const unsigned stride = m.cstride;

  for (r = 0; r < m.r; ++r)
    for (c = 0; c < m.c; ++c)
      v[i++] = data_at(r,c);
}

// write values to a file
void matrix_unload_file(matrix m, char *file_name)
{
  unsigned r, c;

  float *data = m.cpu_data;
  const unsigned stride = m.cstride;

  // open file for writing
  FILE *table = fopen(file_name, "w");

  // check if we were even able to open table file
  if (file_name == NULL)
    errcheck_cpu();

  for (r = 0; r < m.r; ++r)
  {
    for (c = 0; c < m.c-1; ++c)
      fprintf(table, "%.16f ", data_at(r,c));
    fprintf(table, "%.16f\n", data_at(r,c));
  }

  fclose(table);
}

// destroy matrix m
void matrix_dest(matrix *m)
{
  // free matrix data
  free(m->cpu_data);
  free(m->cv);

  // free row/col indices
  //free(m->drc);

  // set everything to zero for safety
  m->r = 0;
  m->c = 0;
  m->rstride = 0;
  m->cstride = 0;
  m->cpu_data = NULL;
}

//** Synchronization

// wait for asynchronous calls to finish
void matrix_wait()
{
  // nada
}

//** Addition/Subtraction

// a = b + c
void matrix_add(matrix a, matrix b, matrix c)
{
  #ifndef NO_ERRCHECK
    if ( (a.r != b.r) || (b.r != c.r) ) {
      fprintf(stderr, __FILE__ " %d: "
              "Rows don't match for addition!\n", __LINE__);
      exit(CPU_ERROR);
    }

    if ( (a.c != b.c) || (b.c != c.c) ) {
      fprintf(stderr, __FILE__ " %d: "
              "Cols don't match for addition!\n", __LINE__);
      exit(CPU_ERROR);
    }
  #endif

  const unsigned n = a.r * a.cstride;

/* naive
  unsigned i;
  float *aval = a.cpu_data;
  float *bval = b.cpu_data;
  float *cval = c.cpu_data;
  for (i = 0; i < n; ++i)
    aval[i] = bval[i] + cval[i]; */

  if (n < add_lim)
    add_small(a.cpu_data, b.cpu_data, c.cpu_data, n);
  else
    add_big(a.cpu_data, b.cpu_data, c.cpu_data, n);
}

// a = b - c
void matrix_sub(matrix a, matrix b, matrix c)
{
  #ifndef NO_ERRCHECK
    if ( (a.r != b.r) || (b.r != c.r) ) {
      fprintf(stderr, __FILE__ " %d: "
              "Rows don't match for subtraction!\n", __LINE__);
      exit(CPU_ERROR);
    }

    if ( (a.c != b.c) || (b.c != c.c) ) {
      fprintf(stderr, __FILE__ " %d: "
              "Cols don't match for subtraction!\n", __LINE__);
      exit(CPU_ERROR);
    }
  #endif

  const unsigned n = a.r * a.cstride;

  if (n < add_lim)
    sub_small(a.cpu_data, b.cpu_data, c.cpu_data, n);
  else
    sub_big(a.cpu_data, b.cpu_data, c.cpu_data, n);
}

//** Multiplication

// a = b * c
void matrix_mult(matrix a, matrix b, matrix c)
{
  #ifndef NO_ERRCHECK
    if ( (b.c != c.r) || (a.r != b.r) || (a.c != c.c) ) {
      fprintf(stderr, __FILE__ " %d: "
              "Dimensions don't match for multiplication!\n", __LINE__);
      exit(CPU_ERROR);
    }
  #endif

  // use cblas for multiplication, not going to beat this.. for now ;)
  cblas_sgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans,
              b.r, c.c, b.c, 1.0f,
              b.cpu_data, b.cstride, c.cpu_data, c.cstride,
              0.0f, a.cpu_data, a.cstride);
}

// a = phi(b * c)
void matrix_mult_phi(matrix a, matrix b, matrix c)
{
  #ifndef NO_ERRCHECK
    if ( (b.c != c.r) || (a.r != b.r) || (a.c != c.c) ) {
      fprintf(stderr, __FILE__ " %d: "
              "Dimensions don't match for mult_phi!\n", __LINE__);
      exit(CPU_ERROR);
    }
  #endif

  // use cblas for multiplication, not going to beat this.. for now ;)
  cblas_sgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans,
              b.r, c.c, b.c, 1.0f,
              b.cpu_data, b.cstride, c.cpu_data, c.cstride,
              0.0f, a.cpu_data, a.cstride);

  // apply phi

  unsigned row, col;

  float *adata = a.cpu_data;
  const unsigned astride = a.cstride;

  // need to optimize small kern and put into functions
  if ((a.r*a.cstride) < add_lim)
  {
    const __m128 mm_s1 = _mm_set1_ps(1.7159f);
    const __m128 mm_s2 = _mm_set1_ps(2.0f/3.0f);
    for (row = 0; row < a.r; ++row)
      for (col = 0; col < a.c; col += 4)
      {
        unsigned i;
        const unsigned base = row*astride+col;

        // inside scalar mult
        __m128 mm_v = _mm_load_ps(adata+base);
        mm_v = _mm_mul_ps(mm_v, mm_s2);
        _mm_store_ps(adata+base, mm_v);

        // apply tanh
        for (i = base; i < base+4; ++i)
          adata[i] = tanh(adata[i]);

        // outside scalar mult
        mm_v = _mm_load_ps(adata+base);
        mm_v = _mm_mul_ps(mm_v, mm_s1);
        _mm_store_ps(adata+base, mm_v);
      }
  }
  else
  {
    const __m128 mm_s1 = _mm_set1_ps(1.7159f);
    const __m128 mm_s2 = _mm_set1_ps(2.0f/3.0f);
    #pragma omp parallel for private(col)
    for (row = 0; row < a.r; ++row)
      for (col = 0; col < a.c; col += 4)
      {
        unsigned i;
        const unsigned base = row*astride+col;

        // inside scalar mult
        __m128 mm_v = _mm_load_ps(adata+base);
        mm_v = _mm_mul_ps(mm_v, mm_s2);
        _mm_store_ps(adata+base, mm_v);

        // apply tanh
        for (i = base; i < base+4; ++i)
          adata[i] = tanh(adata[i]);

        // outside scalar mult
        mm_v = _mm_load_ps(adata+base);
        mm_v = _mm_mul_ps(mm_v, mm_s1);
        _mm_store_ps(adata+base, mm_v);
      }
  }

/*
  unsigned i;
  #pragma omp parallel for
  for (i = 0; i < a.r*a.cstride; ++i)
    adata[i] = phi(adata[i]);
*/
}

//** Transpose

// a = b^T
void matrix_trans(matrix a, matrix b)
{
  #ifndef NO_ERRCHECK
    if ( (a.r != b.c) || (a.c != b.r) ) {
      fprintf(stderr, __FILE__ " %d: "
              "Source and destination dimensions don't match for transpose!\n", __LINE__);
      exit(CPU_ERROR);
    }
  #endif

  unsigned r, c;

  float *adata = a.cpu_data;
  float *bdata = b.cpu_data;

  const unsigned astride = a.cstride;
  const unsigned bstride = b.cstride;

  // need to seperate into smaller functions
  if ((a.r*a.c) < trans_lim)
  {
    for (r = 0; r < a.r; r += 4)
      for (c = 0; c < a.c; c += 4)
      {
        // load 4x4 tile
        float *base = bdata + c*bstride+r;
        __m128 mm_v1 = _mm_load_ps(base            );
        __m128 mm_v2 = _mm_load_ps(base +   bstride);
        __m128 mm_v3 = _mm_load_ps(base + 2*bstride);
        __m128 mm_v4 = _mm_load_ps(base + 3*bstride);

        // transpose 4x4 tile
        _MM_TRANSPOSE4_PS(mm_v1, mm_v2, mm_v3, mm_v4);

        // store 4x4 tile back into a
        base = adata + r*astride+c;
        _mm_store_ps(base,             mm_v1);
        _mm_store_ps(base +   astride, mm_v2);
        _mm_store_ps(base + 2*astride, mm_v3);
        _mm_store_ps(base + 3*astride, mm_v4);
      }
  }
  else
  {
    #pragma omp parallel for private(c)
    for (r = 0; r < a.r; r += 4)
      for (c = 0; c < a.c; c += 4)
      {
        // load 4x4 tile
        float *base = bdata + c*bstride+r;
        __m128 mm_v1 = _mm_load_ps(base            );
        __m128 mm_v2 = _mm_load_ps(base +   bstride);
        __m128 mm_v3 = _mm_load_ps(base + 2*bstride);
        __m128 mm_v4 = _mm_load_ps(base + 3*bstride);

        // transpose 4x4 tile
        _MM_TRANSPOSE4_PS(mm_v1, mm_v2, mm_v3, mm_v4);

        // store 4x4 tile back into a
        base = adata + r*astride+c;
        _mm_store_ps(base,             mm_v1);
        _mm_store_ps(base +   astride, mm_v2);
        _mm_store_ps(base + 2*astride, mm_v3);
        _mm_store_ps(base + 3*astride, mm_v4);
      }
  }

  /* tiled version, fast but big tiles hurt everything else
    #pragma omp parallel for private(tc,r,c)
    for (tr = 0; tr < a.r; tr += trans_tile)
      for (tc = 0; tc < a.c; tc += trans_tile)
        for (r = tr; r < tr+trans_tile; r += 4)
          for (c = tc; c < tc+trans_tile; c += 4)
          {
            // load 4x4 tile
            float *base = bdata + c*bstride+r;
            __m128 mm_v1 = _mm_load_ps(base            );
            __m128 mm_v2 = _mm_load_ps(base +   bstride);
            __m128 mm_v3 = _mm_load_ps(base + 2*bstride);
            __m128 mm_v4 = _mm_load_ps(base + 3*bstride);

            // transpose 4x4 tile
            _MM_TRANSPOSE4_PS(mm_v1, mm_v2, mm_v3, mm_v4);

            // store 4x4 tile back into a
            base = adata + r*astride+c;
            _mm_store_ps(base,             mm_v1);
            _mm_store_ps(base +   astride, mm_v2);
            _mm_store_ps(base + 2*astride, mm_v3);
            _mm_store_ps(base + 3*astride, mm_v4);
          } */
}

//** Pointwise miscellaneous

// a = b^2
void matrix_sqr(matrix a, matrix b)
{
  #ifndef NO_ERRCHECK
    if ( (a.r != b.r) || (a.c != b.c) ) {
      fprintf(stderr, __FILE__ " %d: "
              "Source and destination dimensions don't match for square!\n", __LINE__);
      exit(CPU_ERROR);
    }
  #endif

  const unsigned n = a.r * a.cstride;

  if (n < add_lim)
    sqr_small(a.cpu_data, b.cpu_data, n);
  else
    sqr_big(a.cpu_data, b.cpu_data, n);
}

// scalar multiplication a = b*c
void matrix_scl(matrix a, float b, matrix c)
{
  #ifndef NO_ERRCHECK
    if ( (a.r != c.r) || (a.c != c.c) ) {
      fprintf(stderr, __FILE__ " %d: "
              "Source and destination dimensions don't match for scale!\n", __LINE__);
      exit(CPU_ERROR);
    }
  #endif

  const unsigned n = a.r * a.cstride;

  if (n < add_lim)
    scl_small(a.cpu_data, b, c.cpu_data, n);
  else
    scl_big(a.cpu_data, b, c.cpu_data, n);
}

// scalar multiplication & pointwise addition a = b.*c+d
void matrix_scl_add(matrix a, float b, matrix c, matrix d)
{
  #ifndef NO_ERRCHECK
    if ( (a.r != c.r) || (a.c != c.c) ) {
      fprintf(stderr, __FILE__ " %d: "
              "Source and destination dimensions don't match for scl_add!\n", __LINE__);
      exit(CPU_ERROR);
    }
    if ( (a.r != d.r) || (a.c != d.c) ) {
      fprintf(stderr, __FILE__ " %d: "
              "Source and destination dimensions don't match for scl_add!\n", __LINE__);
      exit(CPU_ERROR);
    }
  #endif

  const unsigned n = a.r * a.cstride;

  if (n < add_lim)
    scl_add_small(a.cpu_data, b, c.cpu_data, d.cpu_data, n);
  else
    scl_add_big(a.cpu_data, b, c.cpu_data, d.cpu_data, n);
}

// pointwise multiplication a = b.*c
void matrix_pmult(matrix a, matrix b, matrix c)
{
  #ifndef NO_ERRCHECK
    if ( (a.r != b.r) || (b.r != c.r) ) {
      fprintf(stderr, __FILE__ " %d: "
              "Rows don't match for pmult!\n", __LINE__);
      exit(CPU_ERROR);
    }

    if ( (a.c != b.c) || (b.c != c.c) ) {
      fprintf(stderr, __FILE__ " %d: "
              "Cols don't match for pmult!\n", __LINE__);
      exit(CPU_ERROR);
    }
  #endif

  const unsigned n = a.r * a.cstride;

  if (n < add_lim)
    pmult_small(a.cpu_data, b.cpu_data, c.cpu_data, n);
  else
    pmult_big(a.cpu_data, b.cpu_data, c.cpu_data, n);
}

//
void matrix_phi_prime(matrix a, matrix b)
{
  unsigned r, c;

  const unsigned astride = a.cstride;
  const unsigned bstride = b.cstride;

  float *adata = a.cpu_data;
  float *bdata = b.cpu_data;

  const __m128 mm_s1 = _mm_set1_ps(2.0f/3.0f);
  const __m128 mm_s2 = _mm_set1_ps(1.7159f);

  // need to optimize small kern and put into functions
  if ((a.r*a.cstride) < add_lim)
  {
    for (r = 0; r < a.r; ++r)
      for (c = 0; c < a.c; c += 4)
      {
        __m128 mm_v1 = _mm_load_ps(bdata + r*bstride+c);

        mm_v1 = _mm_mul_ps(mm_v1, mm_v1);

        mm_v1 = _mm_sub_ps(mm_s2, mm_v1);

        mm_v1 = _mm_mul_ps(mm_s1, mm_v1);

        _mm_store_ps(adata + r*astride+c, mm_v1);
      }
  }
  else
  {
    #pragma omp parallel for private(c)
    for (r = 0; r < a.r; ++r)
      for (c = 0; c < a.c; c += 4)
      {
        __m128 mm_v1 = _mm_load_ps(bdata + r*bstride+c);

        mm_v1 = _mm_mul_ps(mm_v1, mm_v1);

        mm_v1 = _mm_sub_ps(mm_s2, mm_v1);

        mm_v1 = _mm_mul_ps(mm_s1, mm_v1);

        _mm_store_ps(adata + r*astride+c, mm_v1);
      }
  }

//  const float scl = 2.0f/3.0f;
//  return scl * (1.7159f - z*z);

/*
  #pragma omp parallel for private(c)
  for (r = 0; r < a.r; ++r)
    for (c = 0; c < a.c; ++c)
      adata[r*astride + c] =
        phi_prime(bdata[r*bstride + c]);
*/
}

// a = b .* phi_prime(c)
void matrix_pmult_phi_prime(matrix a, matrix b, matrix c)
{
  unsigned i;

  float *adata = a.cpu_data;
  float *bdata = b.cpu_data;
  float *cdata = c.cpu_data;

  const __m128 mm_s1 = _mm_set1_ps(2.0f/3.0f);
  const __m128 mm_s2 = _mm_set1_ps(1.7159f);

  #pragma omp parallel for
  for (i = 0; i < a.r*a.cstride; i += 4)
    {
      __m128 mm_v1 = _mm_load_ps(bdata+i);
      __m128 mm_v2 = _mm_load_ps(cdata+i);

      mm_v2 = _mm_mul_ps(mm_v2, mm_v2);

      mm_v2 = _mm_sub_ps(mm_s2, mm_v2);

      mm_v2 = _mm_mul_ps(mm_s1, mm_v2);

      mm_v2 = _mm_mul_ps(mm_v1, mm_v2);

      _mm_store_ps(adata + i, mm_v2);
    }
}

//
void matrix_delta(matrix delta, matrix y, matrix g)
{
  unsigned r, c;

  #ifndef NO_ERRCHECK
    if ( (delta.r != y.r) || (y.r != g.r) ) {
      fprintf(stderr, __FILE__ " %d: "
              "Rows don't match for delta!\n", __LINE__);
      exit(CPU_ERROR);
    }

    if ( (delta.c != y.c) || (y.c != g.c) ) {
      fprintf(stderr, __FILE__ " %d: "
              "Cols don't match for delta!\n", __LINE__);
      exit(CPU_ERROR);
    }
  #endif

  const unsigned dstride = delta.cstride;
  const unsigned gstride = g.cstride;
  const unsigned ystride = y.cstride;

  float *ddata = delta.cpu_data;
  float *gdata = g.cpu_data;
  float *ydata = y.cpu_data;

  // need to optimize small kern and put into functions
  if ((delta.r*delta.cstride) < add_lim)
  {
    const __m128 mm_s = _mm_set1_ps(2.0f/(g.r*g.c));
    for (r = 0; r < delta.r; ++r)
      for (c = 0; c < delta.c; c += 4)
      {
        __m128 mm_v1 = _mm_load_ps(gdata + r*gstride+c);
        __m128 mm_v2 = _mm_load_ps(ydata + r*ystride+c);

        mm_v1 = _mm_sub_ps(mm_v2, mm_v1);
        mm_v1 = _mm_mul_ps(mm_v1, mm_s);

        _mm_store_ps(ddata + r*dstride+c, mm_v1);
      }
  }
  else
  {
    const __m128 mm_s = _mm_set1_ps(2.0f/(g.r*g.c));
    #pragma omp parallel for private (c)
    for (r = 0; r < delta.r; ++r)
      for (c = 0; c < delta.c; c += 4)
      {
        __m128 mm_v1 = _mm_load_ps(gdata + r*gstride+c);
        __m128 mm_v2 = _mm_load_ps(ydata + r*ystride+c);

        mm_v1 = _mm_sub_ps(mm_v2, mm_v1);
        mm_v1 = _mm_mul_ps(mm_v1, mm_s);

        _mm_store_ps(ddata + r*dstride+c, mm_v1);
      }
  }

/*
  const float denom = 2.0f / (g.r*g.c);

  #pragma omp parallel for private(c)
  for (r = 0; r < delta.r; ++r)
    for (c = 0; c < delta.c; ++c)
    {
      const float gval = gdata[r*gstride + c];
      const float yval = ydata[r*ystride + c];

      // wow, sse could be used all over here!
      ddata[r*dstride + c] =
        (yval - gval) * denom;
    }
*/
}

//** Combination/Separation

// remove last row of m
void matrix_r0(matrix *m)
{
  --(m->r);

  float *data = (m->cpu_data)+((m->r)*(m->cstride));
  const float *data_end = data + m->c;

  while (data < data_end)
    *(data++) = 0.0f;
}

// append a row of 1's to m
void matrix_r1(matrix *m)
{
  float *data = (m->cpu_data)+((m->r)*(m->cstride));
  const float *data_end = data + m->c;

  while (data < data_end)
    *(data++) = 1.0f;

  ++(m->r);
}

//** Error measurement

// rmse between values of actual and approx
float matrix_rmse(matrix actual, matrix approx)
{
  // if dims don't match throw error!

  unsigned r, c;
  const unsigned stride = actual.cstride;
  const unsigned len = actual.r*actual.c;
  float *d1 = actual.cpu_data;
  float *d2 = approx.cpu_data;
  float err = 0.0f;

  // need to use omp, maybe reduction?
  for (r = 0; r < actual.r; ++r)
    for (c = 0; c < actual.c; ++c)
    {
      const unsigned i = r*stride+c;
      err += (d1[i] - d2[i])*(d1[i] - d2[i]);
    }

  return sqrt(err/len);
}

// remove and save last col of m
void matrix_c0v(matrix *m)
{
  unsigned r;
  const unsigned stride = m->cstride;
  const unsigned c = --(m->c);
  float *data = m->cpu_data;

  for (r = 0; r < m->r; ++r)
  {
    float *cur = data + r*stride+c;
    m->cv[r] = *cur;
    *cur = 0.0f;
  }
}

// restore last col of m
void matrix_cv(matrix *m)
{
  unsigned r;

  const unsigned stride = m->cstride;
  const unsigned c = (m->c)++;
  float *data = m->cpu_data;

  for (r = 0; r < m->r; ++r)
    data[r*stride+c] = m->cv[r];
}

//** Output

// print m to standard out
void matrix_print(matrix m)
{
  unsigned r, c;
  const unsigned stride = m.cstride;
  float *data = m.cpu_data;

  // print each value to stdout
  for (r = 0; r < m.r; ++r)
  {
    for (c = 0; c < m.c; ++c)
      printf("%f ", data_at(r,c));
    printf("\n");
  }
}

// print m including padding to standard out
void matrix_print_padded(matrix m)
{
  unsigned r, c;

  // set data and stride for data_at
  const unsigned stride = m.cstride;
  float *data = m.cpu_data;

  // print each value to stdout
  for (r = 0; r < m.rstride; ++r)
  {
    for (c = 0; c < m.cstride; ++c)
      printf("%f ", data_at(r,c));
    printf("\n");
  }
}

max_stars_repo_name: idfah/badger

max_issues_repo_name: idfah/badger

max_forks_repo_name: idfah/badger


text:
/*

Copyright (C) 2019-2021  Jan BOON (Kaetemi) <jan.boon@kaetemi.be>
All rights reserved.

Redistribution and use in source and binary forms, with or without
modification, are permitted provided that the following conditions are met:
1. Redistributions of source code must retain the above copyright notice, this
   list of conditions and the following disclaimer.
2. 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.
3. 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.

*/

/*

Polyverse OÜ base include file for C++.

*/

#pragma once
#ifndef SEV_PLATFORM_H
#define SEV_PLATFORM_H

// Use C math defines for M_PI
#define _USE_MATH_DEFINES

#ifdef _WIN32
// Include Win32.
// Ensure malloc is included before anything else,
// There are some odd macros that may conflict otherwise.
// Include STL algorithm to ensure std::min and std::max
// are used everywhere, instead of min and max macros.
#define WIN32_LEAN_AND_MEAN
#define _WIN32_WINNT 0x0600
#ifdef __cplusplus
#define NOMINMAX
#endif /* __cplusplus */
#include <malloc.h>
#ifdef __cplusplus
#include <algorithm>
using std::max;
using std::min;
#endif /* __cplusplus */
#include <Windows.h>
#ifdef _MSC_VER
// #include <codeanalysis\sourceannotations.h>
#endif
#endif /* _WIN32 */

// C++
#ifdef __cplusplus

// Require C++17
#if defined(_MSC_VER) && (!defined(_HAS_CXX17) || !_HAS_CXX17)
static_assert(false, "C++17 is required");
#endif

// Define null, with color highlight
#ifndef null
constexpr decltype(nullptr) null = nullptr;
#define null null
#endif

// Include STL string and allow string literals.
// Always use sv suffix when declaring string literals.
// Ideally, assign them as `constexpr std::string_view`.
#include <string>
#include <string_view>
using namespace std::string_literals;
using namespace std::string_view_literals;

// Include GSL
// auto _ = gsl::finally([&] { delete xyz; });
#include <gsl/util>

// The usual
#include <functional>

// Force inline
#ifdef _MSC_VER
#	define SEV_FORCE_INLINE __forceinline
#else
#	define SEV_FORCE_INLINE inline __attribute__((always_inline))
#endif

#if defined(_DEBUG) && !defined(NDEBUG)
#define SEV_DEBUG
#else
#define SEV_RELEASE
#endif

// Include debug_break
#include <debugbreak.h>
#define SEV_RELEASE_BREAK() debug_break()

#ifdef SEV_DEBUG
#define SEV_DEBUG_BREAK() debug_break()
#define SEV_ASSERT(cond) do { if (!(cond)) SEV_DEBUG_BREAK(); } while (false)
#define SEV_VERIFY(cond) do { if (!(cond)) SEV_DEBUG_BREAK(); } while (false)
#else
#define SEV_DEBUG_BREAK() do { } while (false)
#define SEV_ASSERT(cond) do { } while (false)
#define SEV_VERIFY(cond) do { cond; } while (false)
#endif

// Library export decl
#ifdef _MSC_VER
#	define SEV_DECL_EXPORT __declspec(dllexport)
#	define SEV_DECL_IMPORT __declspec(dllimport)
#else
#	define SEV_DECL_EXPORT
#	define SEV_DECL_IMPORT
#endif

#if defined(SEV_LIB_EXPORT)
#ifdef _MSC_VER
#endif
#  define SEV_LIB SEV_DECL_EXPORT
#elif defined(SEV_LIB_STATIC)
#  define SEV_LIB
#else
#  define SEV_LIB SEV_DECL_IMPORT
#endif

#endif /* __cplusplus */

#endif /* SEV_PLATFORM_H */

/* end of file */

max_stars_repo_name: kaetemi/libsev

max_issues_repo_name: kaetemi/libsev

max_forks_repo_name: kaetemi/libsev


text:
//! \file openmc_driver.h
//! Driver to initialize and run OpenMC in stages
#ifndef ENRICO_OPENMC_DRIVER_H
#define ENRICO_OPENMC_DRIVER_H

#include "cell_instance.h"
#include "geom.h"
#include "neutronics_driver.h"

#include "openmc/tallies/tally.h"
#include <gsl/gsl>
#include <mpi.h>

#include <vector>

namespace enrico {

//! Driver to initialize and run OpenMC in stages
class OpenmcDriver : public NeutronicsDriver {
public:
  //! One-time initalization of OpenMC and member variables
  //! \param comm An existing MPI communicator used to inialize OpenMC
  explicit OpenmcDriver(MPI_Comm comm);

  //! One-time finalization of OpenMC
  ~OpenmcDriver();

  //! Create energy production tallies for a list of materials
  //! \param[in] materials  Indices into OpenMC's materials array
  void create_tallies(gsl::span<int32_t> materials);

  xt::xtensor<double, 1> heat_source(double power) const final;

  //! Initialization required in each Picard iteration
  void init_step() final;

  //! Runs OpenMC for one Picard iteration
  void solve_step() final;

  //! Writes OpenMC output for given timestep and iteration
  //! \param timestep timestep index
  //! \param iteration iteration index
  void write_step(int timestep, int iteration) final;

  //! Finalization required in each Picard iteration
  void finalize_step() final;

  // Data
  openmc::Tally* tally_;            //!< Fission energy deposition tally
  int32_t index_filter_;            //!< Index in filters arrays for material filter
  std::vector<CellInstance> cells_; //!< Array of cell instances
};

} // namespace enrico

#endif // ENRICO_OPENMC_DRIVER_H

max_stars_repo_name: sphamil/enrico

max_issues_repo_name: sphamil/enrico

max_forks_repo_name: sphamil/enrico


text:
#pragma once

#include "datastructures/MortonIndex.h"

#include <fstream>
#include <gsl/gsl>
#include <iostream>
#include <string>

struct OctreeIndexFileHeader
{
  const char magic[4] = { 'i', 'n', 'd', 'x' };
  uint32_t levels_per_index;
  size_t num_indices;
};

template<unsigned int MaxLevels>
void
write_octree_indices_to_file(const std::string& file_path,
                             gsl::span<MortonIndex<MaxLevels>> indices)
{
  OctreeIndexFileHeader header;
  header.levels_per_index = MaxLevels;
  header.num_indices = indices.size();

  std::ofstream fs{ file_path, std::ios::out | std::ios::binary };
  if (!fs.is_open()) {
    std::cerr << "Could not write index file " << file_path << std::endl;
    return;
  }

  fs.write(reinterpret_cast<char const*>(&header), sizeof(OctreeIndexFileHeader));
  for (auto idx : indices) {
    const auto idx_val = idx.get();
    fs.write(reinterpret_cast<char const*>(&idx_val), sizeof(idx_val));
  }

  fs.flush();
  fs.close();
}

template<unsigned int MaxLevels>
std::vector<MortonIndex<MaxLevels>>
read_octree_indices_from_file(const std::string& file_path)
{
  std::ifstream fs{ file_path, std::ios::in | std::ios::binary };
  if (!fs.is_open()) {
    std::cerr << "Could not read index file " << file_path << std::endl;
    return {};
  }

  fs.unsetf(std::ios::skipws);

  fs.seekg(0, std::ios::end);
  auto fileSize = fs.tellg();
  fs.seekg(0, std::ios::beg);

  std::vector<char> rawData;
  rawData.reserve(fileSize);

  rawData.insert(rawData.begin(), std::istream_iterator<char>(fs), std::istream_iterator<char>());

  const auto& header = *reinterpret_cast<OctreeIndexFileHeader const*>(rawData.data());

  if (header.levels_per_index != MaxLevels) {
    std::cerr << "Reading octree index file with indices that contain " << header.levels_per_index
              << " levels but requested to read into MortonIndexs with " << MaxLevels
              << " levels instead!" << std::endl;
    return {};
  }

  auto indices_begin = rawData.data() + sizeof(OctreeIndexFileHeader);
  const auto indices_end = rawData.data() + rawData.size();

  using IndexType_t = typename MortonIndex<MaxLevels>::Store_t;
  const auto index_binary_size = sizeof(IndexType_t);
  std::vector<MortonIndex<MaxLevels>> indices;
  indices.reserve(header.num_indices);
  for (; indices_begin < indices_end; indices_begin += index_binary_size) {
    indices.emplace_back(*reinterpret_cast<IndexType_t const*>(indices_begin));
  }

  return indices;
}

max_stars_repo_name: igd-geo/schwarzwald

max_issues_repo_name: igd-geo/schwarzwald

max_forks_repo_name: igd-geo/schwarzwald


text:
#ifdef HAVE_CONFIG_H
#include "config.h"
#endif

#include <php.h>
#include <cblas.h>
#include "kernel/operators.h"

/**
 * Sets the number of threads to use when parallel processesing.
 *
 * @param return_value
 * @param threads
 */
void tensor_set_num_threads(zval * return_value, zval * threads)
{
    int n = zephir_get_intval(threads);

    openblas_set_num_threads(n);

    RETURN_TRUE;
}

/**
 * Return the number of threads to use when parallel processesing.
 *
 * @param return_value
 */
void tensor_get_num_threads(zval * return_value)
{
    long threads = openblas_get_num_threads();

    RETURN_LONG(threads);
}