text:
#include <stdio.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_min.h>

double fn1 (double x, void * params)
{
  return cos(x) + 1.0;
}

int
main (void)
{
  int status;
  int iter = 0, max_iter = 100;
  const gsl_min_fminimizer_type *T;
  gsl_min_fminimizer *s;
  double m = 2.0, m_expected = M_PI;
  double a = 0.0, b = 6.0;
  gsl_function F;

  F.function = &fn1;
  F.params = 0;

  T = gsl_min_fminimizer_brent;
  s = gsl_min_fminimizer_alloc (T);
  gsl_min_fminimizer_set (s, &F, m, a, b);

  printf ("using %s method\n",
          gsl_min_fminimizer_name (s));

  printf ("%5s [%9s, %9s] %9s %10s %9s\n",
          "iter", "lower", "upper", "min",
          "err", "err(est)");

  printf ("%5d [%.7f, %.7f] %.7f %+.7f %.7f\n",
          iter, a, b,
          m, m - m_expected, b - a);

  do
    {
      iter++;
      status = gsl_min_fminimizer_iterate (s);

      m = gsl_min_fminimizer_x_minimum (s);
      a = gsl_min_fminimizer_x_lower (s);
      b = gsl_min_fminimizer_x_upper (s);

      status
        = gsl_min_test_interval (a, b, 0.001, 0.0);

      if (status == GSL_SUCCESS)
        printf ("Converged:\n");

      printf ("%5d [%.7f, %.7f] "
              "%.7f %+.7f %.7f\n",
              iter, a, b,
              m, m - m_expected, b - a);
    }
  while (status == GSL_CONTINUE && iter < max_iter);

  gsl_min_fminimizer_free (s);

  return status;
}

ext: c

max_stars_repo_name: ruslankuzmin/julia


text:
#include <math.h>
#include <stdlib.h>
#include <gsl/gsl_complex.h>
#include <gsl/gsl_complex_math.h>
#include <gsl/gsl_poly.h>
#include "userFunc.h"

#define G        6.67384e-8    /* 2010 CODATA value in CGS units */
#define ALPHAMIN 1.0e-3        /* Minimum alpha allowed for numerical reasons */

/**************************************************************************/
/* This defines userFunc routines for the Krumholz & Burkert (2010) GI    */
/* disk problem.                                                          */
/**************************************************************************/

void
userEOS(const double t, const double dt, const grid *grd,
	const double *col, const double *pres, const double *eInt,
	void *params,
	double *gamma, double *delta) {
  fprintf(stderr,
	  "Warning: userEOS function called but not implemented!\n");
  return;
}

void
userAlpha(const double t, const double dt, const grid *grd,
	  const double *col, const double *pres, const double *eInt,
	  const double *gamma, const double *delta,
	  void *params,
	  double *alpha) {

  gsl_poly_complex_workspace *wksp;
  double coef[6], sol[10];
  double omega, kappa, T1, sigma, kcrit, J, Q, d, tfast, tgrowth;
  int i, j;
  const int m=2;
  gsl_complex nu2, nu;
  double alphacoef = ((double *) params)[4];

  /* Allocate workspace */
  wksp = gsl_poly_complex_workspace_alloc(6);

  /* Store the coefficiencts of the dispersion relation polynomial
     that don't change */
  coef[2] = -8.0;
  coef[3] = coef[4] = 0.0;

  /* Loop over cells */
  for (i=0; i<grd->nr; i++) {

    /* Handle negative column densities, which may arise if the time
       step overshoots. We just need to make sure the code doesn't
       barf here, as the negative column densities will be fixed by
       the iterative solver. */
    if (col[i] <= 0.0) {
      alpha[i] = ALPHAMIN;
      continue;
    }

    /* The various quantities needed for the stability analysis */
    sigma = sqrt(pres[i]/col[i]);
    omega = grd->vphi_g[i+1]/grd->r_g[i+1];
    kappa = sqrt(2.0*(1.0+grd->beta_g[i+1]))*omega;
    T1 = -(2.0*m*omega/(kappa*grd->r_g[i+1])) *
      (2.0*m*omega/(kappa*grd->r_g[i+1])) * (grd->beta_g[i+1]-1);
    kcrit = kappa*kappa / (2.0*M_PI*G*col[i]);
    Q = kappa*sigma / (M_PI*G*col[i]);
    if (Q < 0)
      printf("Q < 0! kappa = %f, sigma = %f, col = %f\n", kappa, sigma, col[i]);
    J = sqrt(T1)/kcrit;

    /* Coefficients of the disperison relation polynomial */
    coef[0] = -Q*Q*Q*Q;
    coef[1] = 6.0*Q*Q;
    coef[5] = 16.0*J*J;

    /* Find roots, giving minima and maxima of D */
    gsl_poly_complex_solve(coef, 6, wksp, sol);

    /* For each root, compute the growth time of the instability */
    tfast = 1.0e30;
    for (j=0; j<5; j++) {

      /* Skip roots with negative real part, or imaginary part that is
	 greater than roundoff */
      if ((sol[2*j] < 0) || (fabs(sol[2*j+1]) > 1.0e-6)) continue;

      /* Compute D */
      d = (Q*Q/(4.0*sol[2*j]*sol[2*j]) - 1.0/sol[2*j]) *
	(Q*Q/(4.0*sol[2*j]*sol[2*j]) - 1.0/sol[2*j] -
	 4.0*J*J*sol[2*j]*sol[2*j]);

      /* Ignore stable roots */
      if (d > 0) continue;

      /* Get frequency and growth time */
      GSL_SET_COMPLEX(&nu2,
		      1.0 + 0.5*(Q*Q/(4.0*sol[2*j]*sol[2*j]) - 1.0/sol[2*j]),
		      0.5 * sqrt(-d));
      nu = gsl_complex_sqrt(nu2);
      tgrowth = 1.0 / (GSL_IMAG(nu)*kappa) / (2.0*M_PI/omega);

      /* Store minimum growth time */
      if (tgrowth < tfast) tfast = tgrowth;
    }

    /* Now check for gravitational instability */
    if (Q < 1) {
      GSL_SET_COMPLEX(&nu, 0.0, sqrt(1.0/(Q*Q)-1.0));
      tgrowth = 1.0 / (GSL_IMAG(nu)*kappa) / (2.0*M_PI/omega);
      if (tgrowth < tfast) tfast = tgrowth;
    }

    /* Compute alpha based on the fastest growing mode timescale */
    alpha[i] = alphacoef * exp(-tfast+1.0);
    if (alpha[i] > 1.0) alpha[i] = 1.0;
    if (alpha[i] < ALPHAMIN) alpha[i] = ALPHAMIN;
  }

  /* Free workspace */
  gsl_poly_complex_workspace_free(wksp);
}

void
userMassSrc(const double t, const double dt, const grid *grd,
	    const double *col, const double *pres, const double *eInt,
	    const double *gamma, const double *delta,
	    void *params,
	    double *massSrc) {
  /* Estimate Mdot by
     dSigma/dt = - eta_ML Sigma_SFR
               = - eta_ML eps_ff Sigma_g/t_ff,
  */

  int i;
  double shapefac = ((double *) params)[2];
  double zetad = ((double *) params)[3];
  double etaML = ((double *) params)[5];
  double epsffmax = ((double *) params)[6];
  double rhostar, a, b, c, h;
  double rhog, tff, alphavir, epsff;

  for (i=0; i<grd->nr; i++) {
    /* Stellar density */
    rhostar = shapefac * SQR(grd->vphi_g[i+1]) *
      (1.0+2.0*grd->beta_g[i+1]) /
      (4.0*M_PI*G*SQR(grd->r_g[i+1]));
    /* Get scale height from OML model */
    if (rhostar > 0.0) {
      a = 2.0*M_PI*zetad*G*rhostar*col[i];
      b = M_PI/2.0*G*SQR(col[i]);
      c = -pres[i];
      h = (-b + sqrt(b*b-4.0*a*c))/(2.0*a);
    } else {
      h = pres[i] / (M_PI/2.0*G*SQR(col[i]));
    }
    /* Gas density and free-fall time */
    rhog = col[i]/(2*h);
    tff = sqrt(3*M_PI/(32*G*rhog));
    /* Virial ratio */
    alphavir = pres[i]/(M_PI/2.0*G*col[i]*col[i]*h);
    /* epsff */
    if (alphavir <= 1.0) epsff = epsffmax;
    else epsff = epsffmax*exp(-alphavir);
    /* Mdot */
    massSrc[i] = -(1.0+etaML) * epsff * col[i] / tff;
  }
}

void
userIntEnSrc(const double t, const double dt, const grid *grd,
	     const double *col, const double *pres, const double *eInt,
	     const double *gamma, const double *delta,
	     void *params,
	     double *intEnSrc) {
  /* Cooling rate = eta Sigma sigma^2 Omega = eta P vphi/r */
  int i;
  double eta = ((double *) params)[0];
  double sigmath = ((double *) params)[1];
  double shapefac = ((double *) params)[2];
  double zetad = ((double *) params)[3];
  double sigma2, sigmaNT, rhostar, a, b, c, h;

  for (i=0; i<grd->nr; i++) {
    /* Gas velocity dispersion */
    sigma2 = pres[i]/col[i];
    if (sigma2 > SQR(sigmath)) {
      sigmaNT = sqrt(sigma2 - SQR(sigmath));
      /* Stellar density */
      rhostar = shapefac * SQR(grd->vphi_g[i+1]) *
	(1.0+2.0*grd->beta_g[i+1]) /
	(4.0*M_PI*G*SQR(grd->r_g[i+1]));
      /* Get scale height from OML model */
      if (rhostar > 0.0) {
	a = 2.0*M_PI*zetad*G*rhostar*col[i];
	b = M_PI/2.0*G*SQR(col[i]);
	c = -pres[i];
	h = (-b + sqrt(b*b-4.0*a*c))/(2.0*a);
      } else {
	h = pres[i] / (M_PI/2.0*G*SQR(col[i]));
      }
      /* Cooling rate = eta Sigma sigmaNT^2 / (h/sigmaNT) */
      intEnSrc[i] = -eta * col[i] * SQR(sigmaNT) / (h/sigmaNT);
    } else {
      intEnSrc[i] = 0.0;
    }
  }
}

void
userIBC(const double t, const double dt, const grid *grd,
	const double *col, const double *pres, const double *eInt,
	const double *gamma, const double *delta,
	const pres_bc_type ibc_pres, const enth_bc_type ibc_enth,
	void *params,
	double *ibc_pres_val, double *ibc_enth_val) {
  fprintf(stderr,
	  "Warning: userIBC function called but not implemented!\n");
  return;
}

void
userOBC(const double t, const double dt, const grid *grd,
	const double *col, const double *pres, const double *eInt,
	const double *gamma, const double *delta,
	const pres_bc_type obc_pres, const enth_bc_type obc_enth,
	void *params,
	double *obc_pres_val, double *obc_enth_val) {
  fprintf(stderr,
	  "Warning: userOBC function called but not implemented!\n");
  return;
}

void
userPreTimestep(const double t, const double dt,
		const grid *grd, double *col, double *pres,
		double *eInt, double *mBnd, double *eBnd,
		double *mSrc, double *eSrc,
		void *params, const unsigned long nUserOut,
		double *userOut) {
  fprintf(stderr,
	  "Warning: userPreTimestep function called but not implemented!\n");
  return;
}

void
userPostTimestep(const double t, const double dt,
		 const grid *grd, double *col, double *pres,
		 double *eInt, double *mBnd, double *eBnd,
		 double *mSrc, double *eSrc,
		 void *params, const unsigned long nUserOut,
		 double *userOut) {
  fprintf(stderr,
	  "Warning: userPostTimestep function called but not implemented!\n");
  return;
}

void
userCheckRead(
	      FILE *fp, grid *grd, const unsigned long nOut,
	      double *tOut, double *colOut,
	      double *presOut, double *eIntOut, double *mBndOut,
	      double *eBndOut, double *mSrcOut, double *eSrcOut,
	      const unsigned long nUserOut, double *userOut,
	      void *params
	      ) {
  fprintf(stderr,
	  "Warning: userCheckRead function called but not implemented!\n");
  return;
}

void
userCheckWrite(
	      FILE *fp,
	      const grid *grd, const unsigned long nOut,
	      const double *tOut, const double *colOut,
	      const double *presOut, const double *eIntOut,
	      const double *mBndOut, const double *eBndOut,
	      const double *mSrcOut, const double *eSrcOut,
	      const unsigned long nUserOut, const double *userOut,
	      const void *params
	      ) {
  fprintf(stderr,
	  "Warning: userCheckWrite function called but not implemented!\n");
  return;
}

max_stars_repo_name: franciscaconcha/amuse-vader

max_issues_repo_name: franciscaconcha/amuse-vader

max_forks_repo_name: franciscaconcha/amuse-vader


text:
//          Copyright Jean Pierre Cimalando 2019.
// Distributed under the Boost Software License, Version 1.0.
//    (See accompanying file LICENSE.md or copy at
//          http://www.boost.org/LICENSE_1_0.txt)

#pragma once
#include "keystate.h"
#include <gsl/gsl>
#include <string>
#include <cstdint>

struct Midi_Output {
    std::string id;
    std::string name;
};

enum Midi_Message_Flag {
    Midi_Message_Is_First = 1,
};

///
class Midi_Instrument {
public:
    Midi_Instrument();
    virtual ~Midi_Instrument() {}
    void send_message(const uint8_t *data, unsigned len, double ts, uint8_t flags);
    void initialize();
    void all_sound_off();

    const Keyboard_State &keyboard_state() const noexcept { return kbs_; }

    virtual void flush_events() {}

    virtual void open_midi_output(gsl::cstring_span id) = 0;
    virtual void close_midi_output() = 0;

    virtual bool is_synth() const { return false; }

protected:
    virtual void handle_send_message(const uint8_t *data, unsigned len, double ts, uint8_t flags) = 0;

private:
    Keyboard_State kbs_;
};

ext: h

max_stars_repo_name: jpcima/smf-dsp

max_issues_repo_name: jpcima/smf-dsp

max_forks_repo_name: jpcima/smf-dsp


text:
#include <stdio.h>
#include <gsl/gsl_poly.h>

int main() {
  double a = 1.0, b = -3.0, c = 2.0;
  double x0 = 0.0, x1 = 0.0;
  gsl_poly_solve_quadratic(a, b, c, &x0, &x1);
  printf("%4.1fx^2 + %4.1fx + %4.1f = 0 -> x = %4.1f, %4.1f\n",
      a, b, c, x0, x1);
  return 0;
}

max_stars_repo_name: julnamoo/practice-linux

max_issues_repo_name: julnamoo/practice-linux

max_forks_repo_name: julnamoo/practice-linux


text:
#ifndef OPENMC_TALLIES_FILTER_TIME_H
#define OPENMC_TALLIES_FILTER_TIME_H

#include <gsl/gsl-lite.hpp>

#include "openmc/tallies/filter.h"
#include "openmc/vector.h"

namespace openmc {

//==============================================================================
//! Bins the incident particle time.
//==============================================================================

class TimeFilter : public Filter {
public:
  //----------------------------------------------------------------------------
  // Constructors, destructors

  ~TimeFilter() = default;

  //----------------------------------------------------------------------------
  // Methods

  std::string type() const override { return "time"; }

  void from_xml(pugi::xml_node node) override;

  void get_all_bins(const Particle& p, TallyEstimator estimator,
    FilterMatch& match) const override;

  void to_statepoint(hid_t filter_group) const override;

  std::string text_label(int bin) const override;

  //----------------------------------------------------------------------------
  // Accessors

  const vector<double>& bins() const { return bins_; }
  void set_bins(gsl::span<const double> bins);

protected:
  //----------------------------------------------------------------------------
  // Data members

  vector<double> bins_;
};

} // namespace openmc
#endif // OPENMC_TALLIES_FILTER_ENERGY_H

max_stars_repo_name: RyotaroOKabe/openmc

max_issues_repo_name: RyotaroOKabe/openmc

max_forks_repo_name: RyotaroOKabe/openmc


text:
/*
	Simbicon 1.5 Controller Editor Framework,
	Copyright 2009 Stelian Coros, Philippe Beaudoin and Michiel van de Panne.
	All rights reserved. Web: www.cs.ubc.ca/~van/simbicon_cef

	This file is part of the Simbicon 1.5 Controller Editor Framework.

	Simbicon 1.5 Controller Editor Framework is free software: you can
	redistribute it and/or modify
	it under the terms of the GNU General Public License as published by
	the Free Software Foundation, either version 3 of the License, or
	(at your option) any later version.

	Simbicon 1.5 Controller Editor Framework is distributed in the hope
	that it will be useful, but WITHOUT ANY WARRANTY; without even the
	implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
	See the GNU General Public License for more details.

	You should have received a copy of the GNU General Public License
	along with Simbicon 1.5 Controller Editor Framework.
	If not, see <http://www.gnu.org/licenses/>.
*/

#pragma once

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

#include <MathLib/MathLib.h>
#include <MathLib/ThreeTuple.h>

#include <gsl/matrix/gsl_matrix.h>

#define MATRIX_AT(m, i, j) (*((m->data + ((i) * m->tda + (j)))))

class Vector3d;

/*====================================================================================================================================================================*
 | This class will be used to represent matrices of arbitrary sizes (m rows by n columns) that have elements of type double. The underlying data strucutre used by    |
 | this class is gsl's (Gnu Scientific Library) matrix class. This class also makes use of the ATLAS implementation of BLAS for some operations such as matrix-matrix |
 | multiplication. This class is meant to improve performance, not necessarily ease of use.                                                                           |
 *====================================================================================================================================================================*/

class Matrix  {
friend class Vector;
protected:

	//this data structure holds all the matrix data
	gsl_matrix *matrix;

public:
	/**
		constructor - creates an m rows by n columns matrix that is not initialized to any particular values
	*/
	Matrix(int m, int n);

	/**
		default constructor
	*/
	Matrix();

	/**
		copy constructor - performs a deep copy of the matrix passed in as a parameter.
	*/
	Matrix(const Matrix& other);

	/**
		destructor.
	*/
        virtual ~Matrix();

	/**
		loads the matrix with all zero values.
	*/
	void loadZero();

	/**
		loads the matrix with 1's on the diagonal, 0's everywhere else - note: the matrix doesn't have to be square.
	*/
	void loadIdentity();

	/**
		this method sets the current matrix to a 3x3 matrix that is equal to the outer product of the vectors a and b
	*/
	void setToOuterproduct(const Vector3d& a, const Vector3d& b);

	/**
		this method resizes the current matrix to have m rows and n cols. If not enough space is allocated for it, then a new matrix of
		correct dimensions is allocated. There is no guarantee with regards to the data that is contained in the matrix after a resize operation.
	*/
	void resizeTo(int m, int n);

	/**
		copy operator - performs a deep copy of the matrix passed in as a parameter.
	*/
	Matrix& operator=(const Matrix &other);

	/**
		this method performs a shallow copy of the matrix that is passed in as a parameter.
	*/
	void shallowCopy(const Matrix& other, int startRow = 0, int startCol = 0, int endRow = -1, int endCol = -1);

	/**
		this method performs a deep copy of the matrix that is passed in as a paramerer.
	*/
	void deepCopy(const Matrix& other);

	/**
		Returns the number of columns
	*/
	int getColumnCount() const;

	/**
		Returns the number of rows
	*/
	int getRowCount()const;

	/**
		Multiplies each element in the current matrix by a constant
	*/
	void multiplyBy(const double val);

	/**
		This method sets the current matrix to be equal to one of the products: a * b, a'*b, a*b' or a'*b'.
		The values of transA and transB indicate which of the matrices are tranposed and which ones are not.
	*/
	void setToProductOf(const Matrix& a, const Matrix& b, bool transA = false, bool transB = false);

	/**
		This method computes the inverse of the matrix a and writes it over the current matrix. The implementation
		for the inverse of a matrix was obtained from Graphite. The parameter t is used as a threshold value for
		determinants, etc so that we still get a result for the inverse of our matrix even if it is very poorly conditioned.
	*/
	void setToInverseOf(const Matrix &a, double t = 0);

        /**
                This method prints the contents of the matrix - testing purpose only.
        */
        void printMatrix() const;

        /**
                This method prints the contents of the matrix - testing purpose only.
        */
        std::string toString() const;

	/**
		This method sets the current matrix to be a sub-matrix (starting at (i,j) and ending at (i+rows, j+cols)
		of the one that is passed in as a parameter - shallow copy only.
	*/
	void setToSubmatrix(const Matrix &a, int i, int j, int rows, int cols);

	/**
		This method returns a copy of the value of the matrix at (i,j)
	*/
	double get(int i, int j) const;

	/**
		This method sets the value of the matrix at (i,j) to newVal.
	*/
	void set(int i, int j, double newVal);

	/**
		This method is used to set the values in the matrix to the ones that are passed in the array of doubles.
		It is assumed that the array contains the right number of elements and that there is no space between consecutive rows (tda == nrCols).
	*/
	void setValues(double* vals);

	/**
		this method returns a pointer to its internal matrix data structure
	*/
	gsl_matrix* getMatrixPointer() const;

	/**
		Implement this operator to have a quick way of multiplying 3x3 matrices by vectors - used for dynamics for instance
	*/
	Vector3d operator * (const Vector3d &other);

	/**
		Implement a potentially faster function for 3x3 matrix - vector3d multiplication. Result and v need to be different!!
	*/
	void postMultiplyVector(const Vector3d& v, Vector3d& result);

	/**
		Still need to do (if need be): get row/col vector. Initialize a matrix from a vector so that we can create a 1xn matrix easily.
		make sure dgemv works properly, implement vector class, etc, swap, add, scale matrices, etc
	*/

	/**
		this method adds the matrix that is passed in as a parameter to the current matrix. In addition, it scales, both matrices by the two
		numbers that are passed in as parameters:
		*this = a * *this + b * other.
	*/
	void add(const Matrix& other, double scaleA = 1.0, double scaleB = 1.0);

	/**
		this method subtracts the matrix that is passed in as a parameter from the current matrix. In addition, it scales, both matrices by the two
		numbers that are passed in as parameters:
		*this = a * *this - b * other.
	*/
	void sub(const Matrix& other, double scaleA = 1.0, double scaleB = 1.0);

};

void testMatrixClass();


text:
/*System includes*/
#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include <string.h>
#include <sys/stat.h>
#include <float.h>

/*GSL includes*/
#include <gsl/gsl_vector.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_sf.h>
#include <gsl/gsl_integration.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_roots.h>
#include <gsl/gsl_linalg.h>
#include <gsl/gsl_statistics_double.h>
#include <gsl/gsl_fft_complex.h>
#include <gsl/gsl_complex_math.h>
#include <gsl/gsl_multimin.h>
#include <gsl/gsl_deriv.h>
#include <gsl/gsl_randist.h>
#include <gsl/gsl_cblas.h>
#include <gsl/gsl_blas.h>
#include <pthread.h>
#include <gsl/gsl_sf_exp.h>

/*User includes*/
#include "c_vbgmm_fit.h"

void readInputData(const char *szFile, t_Data *ptData);

void c_vbgmm_fit (double* adX, int nN, int nD, int nK, int* anAssign, int debug, int bAssign);

void readAssigns(const char *szFile, int *anAssign, int nN);

int main()
{
    t_Data tData;
    int *anAssign = NULL;
    int i = 0, j = 0, nN = -1, nD = -1, nK = 0;
    double* adX = NULL;
    const gsl_rng_type * T;
    gsl_rng * r;

    gsl_rng_env_setup();

    T = gsl_rng_default;
    r = gsl_rng_alloc (T);

    readInputData("PCA_transformed_data_gt2000.csv", &tData);

    nN = tData.nN;
    nD = tData.nD;

    adX = (double *) malloc(nN*nD*sizeof(double));
    anAssign = (int *) malloc(nN*sizeof(int));

    readAssigns("clustering_gt2000.csv", anAssign, nN);

    for(i = 0; i < nN; i++){
        if(anAssign[i] > nK){
            nK = anAssign[i];
        }
        for(j = 0; j < nD; j++){
            adX[i*nD + j] = tData.aadX[i][j];
        }
    }

    nK = nK + 1;
    fprintf(stderr,"Run c_vbgmm_fit with %d clusters\n",nK);
    fflush(stderr);

    c_vbgmm_fit (adX, nN, nD, nK, anAssign, FALSE, TRUE);
    for(i = 0; i < nN; i++){
        printf("%d,%d\n",i,anAssign[i]);
    }
    free(adX);
    free(anAssign);
    return 0;
}

void readInputData(const char *szFile, t_Data *ptData)
{
  double  **aadX = NULL;
  int  i = 0, j = 0, nD = 0, nN = 0;
  char *szLine = (char *) malloc(sizeof(char)*MAX_LINE_LENGTH);
  FILE* ifp = NULL;

  if(!szLine)
    goto memoryError;

  ifp = fopen(szFile, "r");

  if(ifp){
    char* szTok   = NULL;
    char* pcError = NULL;

    if(fgets(szLine, MAX_LINE_LENGTH, ifp) == NULL)
      goto formatError;

    szTok = strtok(szLine, DELIM);
    /*count dimensions*/
    while(strtok(NULL, DELIM) != NULL){

      nD++;
    }
    /*count data points*/
    while(fgets(szLine, MAX_LINE_LENGTH, ifp) != NULL){
    	nN++;
    }
    fclose(ifp);

    /*reopen input file*/
    ifp = fopen(szFile, "r");

    if(fgets(szLine, MAX_LINE_LENGTH, ifp) == NULL)
      goto formatError;

    /*allocate memory for dimension names*/
    //ptData->aszDimNames = (char **) malloc(nD*sizeof(char*));
    //if(!ptData->aszDimNames)
      //goto memoryError;

    szTok = strtok(szLine, DELIM);
    /*read in dim names*/
    for(i = 0; i < nD; i++){
      szTok = strtok(NULL, DELIM);
     // ptData->aszDimNames[i] = strdup(szTok);
    }

    /*allocate memory for data matrix*/
    aadX = (double **) malloc(nN*sizeof(double*));
    if(!aadX)
      goto memoryError;
    for(i = 0; i < nN; i++){
      aadX[i] = (double *) malloc(nD*sizeof(double));
      if(!aadX[i])
	goto memoryError;
    }

    /*read in input data*/
    //ptData->aszSampleNames = (char **) malloc(nN*sizeof(char*));
    //if(!ptData->aszSampleNames)
      //goto memoryError;

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

      if(fgets(szLine, MAX_LINE_LENGTH, ifp) == NULL)
	goto formatError;

      szTok = strtok(szLine, DELIM);
     // ptData->aszSampleNames[i] = strdup(szTok);
      for(j = 0; j < nD; j++){
	szTok = strtok(NULL, DELIM);

	aadX[i][j] = strtod(szTok,&pcError);

	if(*pcError != '\0'){
	  goto formatError;
	}
      }
    }
  }
  else{
    fprintf(stderr, "Failed to open abundance data file %s aborting\n", szFile);
    fflush(stderr);
    exit(EXIT_FAILURE);
  }

  free(szLine);
  ptData->nD = nD;
  ptData->nN = nN;
  ptData->aadX = aadX;
  return;

 memoryError:
  fprintf(stderr, "Failed allocating memory in readInputData\n");
  fflush(stderr);
  exit(EXIT_FAILURE);

 formatError:
  fprintf(stderr, "Incorrectly formatted abundance data file\n");
  fflush(stderr);
  exit(EXIT_FAILURE);
}

void readAssigns(const char *szFile, int *anAssign, int nN)
{
  int  i = 0;
  char *szLine = (char *) malloc(sizeof(char)*MAX_LINE_LENGTH);
  FILE* ifp = NULL;

  if(!szLine)
    goto memoryError;

  ifp = fopen(szFile, "r");

  if(ifp){
    char* szTok   = NULL;
    char* pcError = NULL;

    if(fgets(szLine, MAX_LINE_LENGTH, ifp) == NULL)
        goto formatError;

    for(i = 0; i < nN; i++){
        if(fgets(szLine, MAX_LINE_LENGTH, ifp) == NULL)
            goto formatError;

        szTok = strtok(szLine, DELIM);
        szTok = strtok(NULL, DELIM);
    	anAssign[i] = strtod(szTok,&pcError);

	    if(*pcError != '\0'){
	        goto formatError;
	    }
    }
  }
  else{
    fprintf(stderr, "Failed to open abundance data file %s aborting\n", szFile);
    fflush(stderr);
    exit(EXIT_FAILURE);
  }

  free(szLine);

  return;

 memoryError:
  fprintf(stderr, "Failed allocating memory in readInputData\n");
  fflush(stderr);
  exit(EXIT_FAILURE);

 formatError:
  fprintf(stderr, "Incorrectly formatted abundance data file\n");
  fflush(stderr);
  exit(EXIT_FAILURE);
}

max_stars_repo_name: merenlab/CONCOCT

max_issues_repo_name: merenlab/CONCOCT

max_forks_repo_name: merenlab/CONCOCT


text:
#include <stdlib.h>
#include <stdio.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_blas.h>
#include <gsl/gsl_multifit_nlinear.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>

struct data
{
  double *t;
  double *y;
  size_t n;
};

/* model function: a * exp( -1/2 * [ (t - b) / c ]^2 ) */
double
gaussian(const double a, const double b, const double c, const double t)
{
  const double z = (t - b) / c;
  return (a * exp(-0.5 * z * z));
}

int
func_f (const gsl_vector * x, void *params, gsl_vector * f)
{
  struct data *d = (struct data *) params;
  double a = gsl_vector_get(x, 0);
  double b = gsl_vector_get(x, 1);
  double c = gsl_vector_get(x, 2);
  size_t i;

  for (i = 0; i < d->n; ++i)
    {
      double ti = d->t[i];
      double yi = d->y[i];
      double y = gaussian(a, b, c, ti);

      gsl_vector_set(f, i, yi - y);
    }

  return GSL_SUCCESS;
}

int
func_df (const gsl_vector * x, void *params, gsl_matrix * J)
{
  struct data *d = (struct data *) params;
  double a = gsl_vector_get(x, 0);
  double b = gsl_vector_get(x, 1);
  double c = gsl_vector_get(x, 2);
  size_t i;

  for (i = 0; i < d->n; ++i)
    {
      double ti = d->t[i];
      double zi = (ti - b) / c;
      double ei = exp(-0.5 * zi * zi);

      gsl_matrix_set(J, i, 0, -ei);
      gsl_matrix_set(J, i, 1, -(a / c) * ei * zi);
      gsl_matrix_set(J, i, 2, -(a / c) * ei * zi * zi);
    }

  return GSL_SUCCESS;
}

int
func_fvv (const gsl_vector * x, const gsl_vector * v,
          void *params, gsl_vector * fvv)
{
  struct data *d = (struct data *) params;
  double a = gsl_vector_get(x, 0);
  double b = gsl_vector_get(x, 1);
  double c = gsl_vector_get(x, 2);
  double va = gsl_vector_get(v, 0);
  double vb = gsl_vector_get(v, 1);
  double vc = gsl_vector_get(v, 2);
  size_t i;

  for (i = 0; i < d->n; ++i)
    {
      double ti = d->t[i];
      double zi = (ti - b) / c;
      double ei = exp(-0.5 * zi * zi);
      double Dab = -zi * ei / c;
      double Dac = -zi * zi * ei / c;
      double Dbb = a * ei / (c * c) * (1.0 - zi*zi);
      double Dbc = a * zi * ei / (c * c) * (2.0 - zi*zi);
      double Dcc = a * zi * zi * ei / (c * c) * (3.0 - zi*zi);
      double sum;

      sum = 2.0 * va * vb * Dab +
            2.0 * va * vc * Dac +
                  vb * vb * Dbb +
            2.0 * vb * vc * Dbc +
                  vc * vc * Dcc;

      gsl_vector_set(fvv, i, sum);
    }

  return GSL_SUCCESS;
}

void
callback(const size_t iter, void *params,
         const gsl_multifit_nlinear_workspace *w)
{
  gsl_vector *f = gsl_multifit_nlinear_residual(w);
  gsl_vector *x = gsl_multifit_nlinear_position(w);
  double avratio = gsl_multifit_nlinear_avratio(w);
  double rcond;

  (void) params; /* not used */

  /* compute reciprocal condition number of J(x) */
  gsl_multifit_nlinear_rcond(&rcond, w);

  fprintf(stderr, "iter %2zu: a = %.4f, b = %.4f, c = %.4f, |a|/|v| = %.4f cond(J) = %8.4f, |f(x)| = %.4f\n",
          iter,
          gsl_vector_get(x, 0),
          gsl_vector_get(x, 1),
          gsl_vector_get(x, 2),
          avratio,
          1.0 / rcond,
          gsl_blas_dnrm2(f));
}

void
solve_system(gsl_vector *x, gsl_multifit_nlinear_fdf *fdf,
             gsl_multifit_nlinear_parameters *params)
{
  const gsl_multifit_nlinear_type *T = gsl_multifit_nlinear_trust;
  const size_t max_iter = 200;
  const double xtol = 1.0e-8;
  const double gtol = 1.0e-8;
  const double ftol = 1.0e-8;
  const size_t n = fdf->n;
  const size_t p = fdf->p;
  gsl_multifit_nlinear_workspace *work =
    gsl_multifit_nlinear_alloc(T, params, n, p);
  gsl_vector * f = gsl_multifit_nlinear_residual(work);
  gsl_vector * y = gsl_multifit_nlinear_position(work);
  int info;
  double chisq0, chisq, rcond;

  /* initialize solver */
  gsl_multifit_nlinear_init(x, fdf, work);

  /* store initial cost */
  gsl_blas_ddot(f, f, &chisq0);

  /* iterate until convergence */
  gsl_multifit_nlinear_driver(max_iter, xtol, gtol, ftol,
                              callback, NULL, &info, work);

  /* store final cost */
  gsl_blas_ddot(f, f, &chisq);

  /* store cond(J(x)) */
  gsl_multifit_nlinear_rcond(&rcond, work);

  gsl_vector_memcpy(x, y);

  /* print summary */

  fprintf(stderr, "NITER         = %zu\n", gsl_multifit_nlinear_niter(work));
  fprintf(stderr, "NFEV          = %zu\n", fdf->nevalf);
  fprintf(stderr, "NJEV          = %zu\n", fdf->nevaldf);
  fprintf(stderr, "NAEV          = %zu\n", fdf->nevalfvv);
  fprintf(stderr, "initial cost  = %.12e\n", chisq0);
  fprintf(stderr, "final cost    = %.12e\n", chisq);
  fprintf(stderr, "final x       = (%.12e, %.12e, %12e)\n",
          gsl_vector_get(x, 0), gsl_vector_get(x, 1), gsl_vector_get(x, 2));
  fprintf(stderr, "final cond(J) = %.12e\n", 1.0 / rcond);

  gsl_multifit_nlinear_free(work);
}

int
main (void)
{
  const size_t n = 300;  /* number of data points to fit */
  const size_t p = 3;    /* number of model parameters */
  const double a = 5.0;  /* amplitude */
  const double b = 0.4;  /* center */
  const double c = 0.15;  /* width */
  const gsl_rng_type * T = gsl_rng_default;
  gsl_vector *f = gsl_vector_alloc(n);
  gsl_vector *x = gsl_vector_alloc(p);
  gsl_multifit_nlinear_fdf fdf;
  gsl_multifit_nlinear_parameters fdf_params =
    gsl_multifit_nlinear_default_parameters();
  struct data fit_data;
  gsl_rng * r;
  size_t i;

  gsl_rng_env_setup ();
  r = gsl_rng_alloc (T);

  fit_data.t = malloc(n * sizeof(double));
  fit_data.y = malloc(n * sizeof(double));
  fit_data.n = n;

  /* generate synthetic data with noise */
  for (i = 0; i < n; ++i)
    {
      double t = (double)i / (double) n;
      double y0 = gaussian(a, b, c, t);
      double dy = gsl_ran_gaussian (r, 0.1 * y0);

      fit_data.t[i] = t;
      fit_data.y[i] = y0 + dy;
    }

  /* define function to be minimized */
  fdf.f = func_f;
  fdf.df = func_df;
  fdf.fvv = func_fvv;
  fdf.n = n;
  fdf.p = p;
  fdf.params = &fit_data;

  /* starting point */
  gsl_vector_set(x, 0, 1.0);
  gsl_vector_set(x, 1, 0.0);
  gsl_vector_set(x, 2, 1.0);

  fdf_params.trs = gsl_multifit_nlinear_trs_lmaccel;
  solve_system(x, &fdf, &fdf_params);

  /* print data and model */
  {
    double A = gsl_vector_get(x, 0);
    double B = gsl_vector_get(x, 1);
    double C = gsl_vector_get(x, 2);

    for (i = 0; i < n; ++i)
      {
        double ti = fit_data.t[i];
        double yi = fit_data.y[i];
        double fi = gaussian(A, B, C, ti);

        printf("%f %f %f\n", ti, yi, fi);
      }
  }

  gsl_vector_free(f);
  gsl_vector_free(x);
  gsl_rng_free(r);

  return 0;
}

max_stars_repo_name: zzpwahaha/Chimera-Control-Trim

max_issues_repo_name: zzpwahaha/Chimera-Control-Trim

max_forks_repo_name: zzpwahaha/Chimera-Control-Trim


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

/*
 * This generator is taken from
 *
 * Donald E. Knuth
 * The Art of Computer Programming
 * Volume 2
 * Third Edition
 * Addison-Wesley
 * Page 108
 *
 * It is called "Fishman"
 *
 * This implementation copyright (C) 2001 Carlo Perassi
 * and (C) 2003 Heiko Bauke.
 */

#include <config.h>
#include <stdlib.h>
#include <gsl/gsl_rng.h>

static inline unsigned long int ran_get (void *vstate);
static double ran_get_double (void *vstate);
static void ran_set (void *state, unsigned long int s);

static const long int m = 2147483647, a = 48271, q = 44488, r = 3399;

typedef struct
{
  unsigned long int x;
}
ran_state_t;

static inline unsigned long int
ran_get (void *vstate)
{
  ran_state_t *state = (ran_state_t *) vstate;

  const unsigned long int x = state->x;

  const long int h = x / q;
  const long int t = a * (x - h * q) - h * r;

  if (t < 0)
    {
      state->x = t + m;
    }
  else
    {
      state->x = t;
    }

  return state->x;
}

static double
ran_get_double (void *vstate)
{
  ran_state_t *state = (ran_state_t *) vstate;

  return ran_get (state) / 2147483647.0;
}

static void
ran_set (void *vstate, unsigned long int s)
{
  ran_state_t *state = (ran_state_t *) vstate;

  if ((s%m) == 0)
    s = 1;                      /* default seed is 1 */

  state->x = s & m;

  return;
}

static const gsl_rng_type ran_type = {
  "fishman20",                  /* name */
  2147483646,                   /* RAND_MAX */
  1,                            /* RAND_MIN */
  sizeof (ran_state_t),
  &ran_set,
  &ran_get,
  &ran_get_double
};

const gsl_rng_type *gsl_rng_fishman20 = &ran_type;

max_stars_repo_name: ielomariala/Hex-Game

max_issues_repo_name: Brian-ning/HMNE

max_forks_repo_name: Brian-ning/HMNE


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

#define N 1000000
#define REP 500

/* compilar: gcc fonte.c -o exec -lblas */

int main( void ) {
  float *x, *y, *y2;
  int i, j;
  double ti = 0, tf = 0; /* ti = tempo inicial // tf = tempo final */
  struct timeval tempo_inicio, tempo_fim;

  /* aloca vetores dianmicamente */

  x = ( float* ) malloc( N * sizeof( float ) );
  y = ( float* ) malloc( N * sizeof( float ) );
  y2 = ( float* ) malloc( N * sizeof( float ) );
  /* atribue valores aleatorios */
  for( i = 0; i < N; i++ ) {
    x[ i ] = random( );
    y2[ i ] = y[ i ] = random( );
  }
  /* marca tempo inicial */
  gettimeofday( &tempo_inicio, NULL );
  /* Calcula y=alfa*x+y */
  for( j = 0; j < REP; j++ ) {
    cblas_saxpy( N, 0.5, x, 1, y, 1 );
  }
  /* calcula o tempo gasto */
  gettimeofday( &tempo_fim, NULL );
  tf = ( double ) tempo_fim.tv_usec / 1000000.0 + tempo_fim.tv_sec;
  ti = ( double ) tempo_inicio.tv_usec / 1000000.0 + tempo_inicio.tv_sec;
  printf( "Tempo gasto em saxpy (BLAS) %.15f segundos \n", tf - ti );

  /* marca tempo inicial */
  gettimeofday( &tempo_inicio, NULL );
  /* Calcula y=alfa*x+y */
  for( j = 0; j < REP; j++ ) {
    for( i = 0; i < N; i++ ) {
      y2[ i ] = 0.5 * x[ i ] + y2[ i ];
    }
  }
  /* calcula o tempo gasto */
  gettimeofday( &tempo_fim, NULL );
  tf = ( double ) tempo_fim.tv_usec / 1000000.0 + tempo_fim.tv_sec;
  ti = ( double ) tempo_inicio.tv_usec / 1000000.0 + tempo_inicio.tv_sec;
  printf( "Tempo gasto direto %.15f segundos \n", tf - ti );

  return( 0 );
}

max_stars_repo_name: lellisls/PAD

max_issues_repo_name: lellisls/PAD

max_forks_repo_name: lellisls/PAD


text:
// --------------------------------------------------------------
// File: ex10.c
// Description: Example with the sections directive
// Software stack: Exercises suite for V. Keller lecture
// Version: 1.0
// License: BSD
// Author: Vincent Keller (Vincent.Keller@epfl.ch), CADMOS, 2012
//
// COMPILATION AND LINKING
//
// WITHOUT MKL :  gcc -O3 -ftree-vectorize dgemm.c -lgsl -lgslcblas -lm -o dgemm
// WITH MKL : export MKLROOT=/software/intel/mkl; icc -DMKLINUSE -DMKL_ILP64 -openmp -I${MKLROOT}/include -O3 -xHost dgemm.c  -L${MKLROOT}/lib/intel64 -lmkl_intel_ilp64 -lmkl_core -lmkl_intel_thread -lpthread -lm -o dgemm
// --------------------------------------------------------------

#include <stdio.h>
#include <stdlib.h>
#include <omp.h>

#if defined(MKLINUSE)
#include "mkl.h"
#else
#include <gsl/gsl_blas.h>
#endif

double verification(double ** array, int N){
        double ret;
        int i,j;
        for (i=0;i<N;i++){
                for (j=0;j<N;j++){
                        ret+=array[i][j];
                }
        }
        return ret;
}

double second()
{
        struct timeval tp;
        int i;
        i = gettimeofday(&tp,NULL);
        return ( (double) tp.tv_sec + (double) tp.tv_usec * 1.e-6 );
}

int main( )
{
        int N, myrank,i,j,k,chunk, mysize;
        double ** A;
        double ** B;
        double ** C;
        double t1, t2;

        double *A_dgemm, *B_dgemm, *C_dgemm;
        double alpha, beta;

        mysize = 1;
        N=2000;

#if defined(MKLINUSE)
        A_dgemm = (double *)mkl_malloc( N*N*sizeof( double ), 64 );
        B_dgemm = (double *)mkl_malloc( N*N*sizeof( double ), 64 );
        C_dgemm = (double *)mkl_malloc( N*N*sizeof( double ), 64 );
#else
        A_dgemm = (double *)malloc( N*N*sizeof( double ));
        B_dgemm = (double *)malloc( N*N*sizeof( double ));
        C_dgemm = (double *)malloc( N*N*sizeof( double ));
#endif

        for (i = 0; i < (N*N); i++) {
                A_dgemm[i] = 1.0;
        }

        for (i = 0; i < (N*N); i++) {
                B_dgemm[i] = 2.0;
        }

        for (i = 0; i < (N*N); i++) {
                C_dgemm[i] = 0.0;
        }

        alpha= 1.0;
        beta = 0.0;

        A = (double**)malloc(N*sizeof(double*));
        B = (double**)malloc(N*sizeof(double*));
        C = (double**)malloc(N*sizeof(double*));

        for (i=0;i<N;i++){
                A[i]=(double*)malloc(N*sizeof(double));
                B[i]=(double*)malloc(N*sizeof(double));
                C[i]=(double*)malloc(N*sizeof(double));
        }

        for (i=0;i<N;i++){
                for (j=0;j<N;j++){
                        A[i][j] = 1.0;
                        B[i][j] = 2.0;
                        C[i][j] = 0.0;
                }
        }

// ============================================================================================================
        t1 = second();
        for (i=0;i<N;i++){
                for (j=0;j<N;j++){
                        for (k=0;k<N;k++){
                                C[i][j]=C[i][j] + A[i][k]*B[k][j];
                        }
                }
        }
        t2 = second();

        printf("[IJK]   Compute time [s]   : %6.3f \n", (t2-t1));
        printf("[IJK]   Performance  [GF/s]: %6.3f \n", (((double)N*(double)N*(double)N*2.0)/((t2-t1)*1.e9)));
        printf("[IJK]   Verification       : %6.3f \n\n", verification(C,N));

// ============================================================================================================
        for (i=0;i<N;i++){
                for (j=0;j<N;j++){
                        A[i][j] = 1.0;
                        B[i][j] = 2.0;
                        C[i][j] = 0.0;
                }
        }
        t1 = second();
        for (i=0;i<N;i++){
                for (k=0;k<N;k++){
                        for (j=0;j<N;j++){
                                C[i][j]=C[i][j] + A[i][k]*B[k][j];
                        }
                }
        }
        t2 = second();

        printf("[IKJ]   Compute time [s]   : %6.3f \n", (t2-t1));
        printf("[IKJ]   Performance  [GF/s]: %6.3f \n", (((double)N*(double)N*(double)N*2.0)/((t2-t1)*1.e9)));
        printf("[IKJ]   Verification       : %6.3f \n\n", verification(C,N));

// ============================================================================================================
        for (i=0;i<N;i++){
                for (j=0;j<N;j++){
                        A[i][j] = 1.0;
                        B[i][j] = 2.0;
                        C[i][j] = 0.0;
                }
        }
        t1 = second();
        for (j=0;j<N;j++){
                for (i=0;i<N;i++){
                        for (k=0;k<N;k++){
                                C[i][j]=C[i][j] + A[i][k]*B[k][j];
                        }
                }
        }
        t2 = second();

        printf("[JIK]   Compute time [s]   : %6.3f \n", (t2-t1));
        printf("[JIK]   Performance  [GF/s]: %6.3f \n", (((double)N*(double)N*(double)N*2.0)/((t2-t1)*1.e9)));
        printf("[JIK]   Verification       : %6.3f \n\n", verification(C,N));

// ============================================================================================================
        for (i=0;i<N;i++){
                for (j=0;j<N;j++){
                        A[i][j] = 1.0;
                        B[i][j] = 2.0;
                        C[i][j] = 0.0;
                }
        }
        t1 = second();
        for (j=0;j<N;j++){
                for (k=0;k<N;k++){
                        for (i=0;i<N;i++){
                                C[i][j]=C[i][j] + A[i][k]*B[k][j];
                        }
                }
        }
        t2 = second();

        printf("[JKI]   Compute time [s]   : %6.3f \n", (t2-t1));
        printf("[JKI]   Performance  [GF/s]: %6.3f \n", (((double)N*(double)N*(double)N*2.0)/((t2-t1)*1.e9)));
        printf("[JKI]   Verification       : %6.3f \n\n", verification(C,N));

// ============================================================================================================
        for (i=0;i<N;i++){
                for (j=0;j<N;j++){
                        A[i][j] = 1.0;
                        B[i][j] = 2.0;
                        C[i][j] = 0.0;
                }
        }
        t1 = second();
        for (k=0;k<N;k++){
                for (i=0;i<N;i++){
                        for (j=0;j<N;j++){
                                C[i][j]=C[i][j] + A[i][k]*B[k][j];
                        }
                }
        }
        t2 = second();

        printf("[KIJ]   Compute time [s]   : %6.3f \n", (t2-t1));
        printf("[KIJ]   Performance  [GF/s]: %6.3f \n", (((double)N*(double)N*(double)N*2.0)/((t2-t1)*1.e9)));
        printf("[KIJ]   Verification       : %6.3f \n\n", verification(C,N));

// ============================================================================================================
        for (i=0;i<N;i++){
                for (j=0;j<N;j++){
                        A[i][j] = 1.0;
                        B[i][j] = 2.0;
                        C[i][j] = 0.0;
                }
        }
        t1 = second();
        for (k=0;k<N;k++){
                for (j=0;j<N;j++){
                        for (i=0;i<N;i++){
                                C[i][j]=C[i][j] + A[i][k]*B[k][j];
                        }
                }
        }
        t2 = second();

        printf("[KJI]   Compute time [s]   : %6.3f \n", (t2-t1));
        printf("[KJI]   Performance  [GF/s]: %6.3f \n", (((double)N*(double)N*(double)N*2.0)/((t2-t1)*1.e9)));
        printf("[KJI]   Verification       : %6.3f \n\n", verification(C,N));

// ============================================================================================================
        t1 = second();
        cblas_dgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans, N,N,N, alpha, A_dgemm, N, B_dgemm, N, beta, C_dgemm, N);
        t2 = second();

        printf("[DGEMM OPT]   Compute time [s]   : %6.3f \n", (t2-t1));
        printf("[DGEMM OPT]   Performance  [GF/s]: %6.3f \n", (((double)N*(double)N*(double)N*2.0)/((t2-t1)*1.e9)));
        printf("[DGEMM OPT]   Verification       : %6.3f \n\n", verification(C,N));

        for (i=0;i<N;i++){
                free(A[i]);
                free(B[i]);
                free(C[i]);
        }

        free(A);
        free(B);
        free(C);

        return 0;
}


text:
/******************************************************************************
CosmoLike Configuration Space Covariances for Projected Galaxy 2-Point Statistics
https://github.com/CosmoLike/CosmoCov
by CosmoLike developers
******************************************************************************/

#include <math.h>
#include <stdlib.h>
#include <stdio.h>
#include <assert.h>
// #include "../class/include/class.h"

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

//void omega_a(double aa,double *om_m,double *om_v);
double omv_vareos(double a);
static inline double hoverh0(double a);
double growfac(double a);
//int func_for_growfac(double a,const double y[],double f[],void *params);
double Tsqr_EH_wiggle(double khoverMPC);
//double int_for_sigma_r_sqr(double k, void * args);
double sigma_r_sqr();
//double Delta_L_wiggle(double k);
//double Delta_lin_wiggle(double k,double a);
double p_lin(double k,double a);
//double int_sig_R_knl(double logk, void *args);
//double int_neff(double lnk, void *args);
//double int_cur(double lnk, void *args);
//void nonlin_scale(double amp, double *R_NL, double *neff, double *Curv);
//double Halofit(double k, double amp, double omm, double omv,double w_z, double R_NL, double neff,double Curv, double P_delta_Lin);
//void Delta_halofit(double **table_P_NL,double logkmin, double logkmax, double dk, double da);
//double Delta_NL_Halofit(double k_NL, double a); //k in h/Mpc
double Pdelta(double k_NL,double a); //k in coverH0 units

//double int_for_chi(double a,void * args);
double f_K(double chi);
double chi(double a);
double a_chi(double chi1);
// extern void emu(double *xstar, double *ystar, double *kstar);

//c%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
//variable Omega_v
//c%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

double omv_vareos(double a)
{
  return(cosmology.Omega_v*exp(-3.*((cosmology.w0+cosmology.wa+1.)*log(a)+cosmology.wa*(1.-a))));
}

//c%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
//c evolution of omega matter and omega lamda with expansion factor

void omega_a(double aa,double *om_m,double *om_v)
{
  double a2,omega_curv;
  a2=aa*aa;
  omega_curv=1.0-cosmology.Omega_m- cosmology.Omega_v;
  *om_m=cosmology.Omega_m /(cosmology.Omega_m +aa*(omv_vareos(aa) *a2 +omega_curv));
  *om_v=omv_vareos(aa)*a2*aa/(cosmology.Omega_m+aa*(a2*omv_vareos(aa) +omega_curv));
}
//c%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
//growth factor including Dark energy parameters w0, wa
//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

//function for growfac (DGL)
int func_for_growfac(double a,const double y[],double f[],void *params)
{
  //double *p=(double *)params;
  if (a == 0) {
    printf("a=0 in function 'func_for_growfac'!\n");
    exit(1);
  }
  double aa=a*a;
  double omegam=cosmology.Omega_m/(aa*a);
  double omegav=omv_vareos(a);
  double hub = hoverh0(a);
  double one_plus_mg_mu = 1.;
  hub = hub*hub;
  f[0]=y[1];
  if(cosmology.MGmu != 0){
    one_plus_mg_mu += cosmology.MGmu*omegav/hub/cosmology.Omega_v;
  }
  f[1]=y[0]*3.*cosmology.Omega_m/(2.*hub*aa*aa*a)*one_plus_mg_mu-y[1]/a*(2.-(omegam+(3.*(cosmology.w0+cosmology.wa*(1.-a))+1)*omegav)/(2.*hub));
  return GSL_SUCCESS;
}

static inline double hoverh0(double a){
  return sqrt(cosmology.Omega_m /(a*a*a) + (1.-cosmology.Omega_m -cosmology.Omega_v )/(a*a) + omv_vareos(a) );
}

double growfac(double a)
{
  const double MINA=1.e-8;
  static cosmopara C;
  static double *ai;
  static double *table;
  double res;

  gsl_interp *intf=gsl_interp_alloc(gsl_interp_linear,Ntable.N_a);
  gsl_interp_accel *acc=gsl_interp_accel_alloc();

  if (recompute_expansion(C))
  {

    if(table!=0) free_double_vector(table,0, Ntable.N_a-1);
    if(ai!=0) free_double_vector(ai,0, Ntable.N_a-1);
    ai=create_double_vector(0, Ntable.N_a-1);
    table=create_double_vector(0, Ntable.N_a-1);

    int i;

    const gsl_odeiv_step_type *T=gsl_odeiv_step_rkf45;
    gsl_odeiv_step *s=gsl_odeiv_step_alloc(T,2);
    gsl_odeiv_control *c=gsl_odeiv_control_y_new(1.e-6,0.0);
    gsl_odeiv_evolve *e=gsl_odeiv_evolve_alloc(2);

    double t=MINA;            //start a
    double t1=1.1;                //final a
    double h=1.e-6;              //initial step size
    double y[2]={MINA,MINA};   //initial conditions
    double norm;
    double par[0]={};
    gsl_odeiv_system sys={func_for_growfac,NULL,2,&par};

    for (i=1;i<=Ntable.N_a;i++) {
      ai[i-1]=i*t1/(1.*Ntable.N_a);
      while(t<ai[i-1])
        gsl_odeiv_evolve_apply(e,c,s,&sys,&t,ai[i-1],&h,y);
      if (i==1) norm=y[0]/ai[i-1];
      table[i-1]=y[0]/norm;
    }

    gsl_odeiv_evolve_free(e);
    gsl_odeiv_control_free(c);
    gsl_odeiv_step_free(s);
    update_cosmopara(&C);

  }
  gsl_interp_init(intf,ai,table,Ntable.N_a);
  res=gsl_interp_eval(intf,ai,table,a,acc);
  gsl_interp_accel_free(acc);
  gsl_interp_free(intf);
  return(res);
}

// ---------------------------- Transfer Function from EH98 ----------------------
//Input: k -- Wavenumber at which to calculate transfer function, in Mpc^-1. Output: Returns the value of the full transfer function fitting formula. This is the form given in Section 3 of Eisenstein & Hu (1997). Notes: Units are Mpc, not h^-1 Mpc.

double Tsqr_EH_wiggle(double khoverMPC)
{
  static double omhh=-123.;
  static double obhh=-123.;
  static double OMEGA_V = -123.;
  static double f_baryon;

  static double k_equality;
  static double sound_horizon;
  static double beta_c;
  static double alpha_c;
  static double beta_node;
  static double alpha_b;
  static double beta_b;
  static double k_silk;

  //if (omhh != cosmology.Omega_m*cosmology.h0*cosmology.h0 || obhh != cosmology.omb*cosmology.h0*cosmology.h0|| OMEGA_V != cosmology.Omega_v){
  double theta_cmb,z_equality,z_drag,R_drag,R_equality;
  double z_drag_b1, z_drag_b2;
  double alpha_c_a1, alpha_c_a2, beta_c_b1, beta_c_b2, alpha_b_G, y;

  omhh = cosmology.Omega_m*cosmology.h0*cosmology.h0;
  obhh = cosmology.omb*cosmology.h0*cosmology.h0;
  OMEGA_V = cosmology.Omega_v;
  f_baryon = obhh/omhh;
    //printf("%le\n",f_baryon);

    theta_cmb=2.728/2.7;// Tcmb in units of 2.7 K
    z_equality= 2.50e4*omhh/POW4(theta_cmb);//Redshift of matter-radiation equality, really 1+z

    k_equality=0.0746*omhh/SQR(theta_cmb);//Scale of equality, in Mpc^-1
    z_drag_b1 = 0.313*pow(omhh,-0.419)*(1+0.607*pow(omhh,0.674));
    z_drag_b2 = 0.238*pow(omhh,0.223);
    z_drag=1291*pow(omhh,0.251)/(1+0.659*pow(omhh,0.828))*(1+z_drag_b1*pow(obhh,z_drag_b2));//Redshift of drag epoch
    R_drag=31.5*obhh/POW4(theta_cmb)*(1000/(1+z_drag));//Photon-baryon ratio at drag epoch

    R_equality= 31.5*obhh/POW4(theta_cmb)*(1000/z_equality);//Photon-baryon ratio at equality epoch
    sound_horizon=2./3./k_equality*sqrt(6./R_equality)*log((sqrt(1+R_drag)+sqrt(R_drag+R_equality))/(1+sqrt(R_equality)));//Sound horizon at drag epoch, in Mpc
    k_silk= 1.6*pow(obhh,0.52)*pow(omhh,0.73)*(1+pow(10.4*omhh,-0.95));//Silk damping scale, in Mpc^-1
    alpha_c_a1 = pow(46.9*omhh,0.670)*(1+pow(32.1*omhh,-0.532));
    alpha_c_a2 = pow(12.0*omhh,0.424)*(1+pow(45.0*omhh,-0.582));
    alpha_c=pow(alpha_c_a1,-f_baryon)*pow(alpha_c_a2,-CUBE(f_baryon)); //CDM suppression

    beta_c_b1 = 0.944/(1+pow(458*omhh,-0.708));
    beta_c_b2 = pow(0.395*omhh, -0.0266);
    beta_c=1.0/(1+beta_c_b1*(pow(1-f_baryon, beta_c_b2)-1));//CDM log shift

    y = z_equality/(1+z_drag);
    alpha_b_G = y*(-6.*sqrt(1+y)+(2.+3.*y)*log((sqrt(1+y)+1)/(sqrt(1+y)-1)));
    alpha_b=2.07*k_equality*sound_horizon*pow(1+R_drag,-0.75)*alpha_b_G;//Baryon suppression
    beta_node = 8.41*pow(omhh, 0.435);//Sound horizon shift
    beta_b= 0.5+f_baryon+(3.-2.*f_baryon)*sqrt(pow(17.2*omhh,2.0)+1);//Baryon envelope shift
  //}
  // Start of TFfit_onek from the original tf_fit.c at http://background.uchicago.edu/~whu/transfer/transferpage.html
    double T_c_ln_beta,T_c_ln_nobeta , T_c_C_alpha, T_c_C_noalpha;
    double q,qsqr, xx, xx_tilde;
    double T_c_f, T_c, s_tilde, T_b_T0, T_b;

  double k=khoverMPC*cosmology.h0; //internally this routine uses Mpc^-1 not h/Mpc
  q = k/13.41/k_equality;
  qsqr =SQR(q);
  xx = k*sound_horizon;

  T_c_ln_beta = log(2.718282+1.8*beta_c*q);
  T_c_ln_nobeta = log(2.718282+1.8*q);

  T_c_ln_beta = log(2.718282+1.8*beta_c*q);
  T_c_C_alpha = 14.2/alpha_c + 386.0/(1+69.9*pow(q,1.08));
  T_c_C_noalpha = 14.2 + 386.0/(1+69.9*pow(q,1.08));

  T_c_f = 1.0/(1.0+POW4(xx/5.4));
  T_c = T_c_f*T_c_ln_beta/(T_c_ln_beta+T_c_C_noalpha*qsqr) +(1-T_c_f)*T_c_ln_beta/(T_c_ln_beta+T_c_C_alpha*qsqr);

  s_tilde = sound_horizon*pow(1+CUBE(beta_node/xx),-1./3.);
  xx_tilde = k*s_tilde;

  T_b_T0 = T_c_ln_nobeta/(T_c_ln_nobeta+T_c_C_noalpha*qsqr);

  T_b = sin(xx_tilde)/(xx_tilde)*(T_b_T0/(1+SQR(xx/5.2))+alpha_b/(1+CUBE(beta_b/xx))*exp(-pow(k/k_silk,1.4)));

  return SQR(f_baryon*T_b + (1-f_baryon)*T_c);
}

//Calculate Normalization see Cosmology Notes 8.105
double int_for_sigma_r_sqr(double k, void * args)
{
  double kR, res, x;
  kR = k*8.; // r=8 Mpc/h
  x = (sin(kR) - kR*cos(kR))/(kR*kR*kR);
  res = pow(k,2.+cosmology.n_spec+ 0.5*cosmology.alpha_s*log(k/0.05))*Tsqr_EH_wiggle(k)*x*x;
  return res;
}

double sigma_r_sqr()
{
  static double res = -123.;
  static cosmopara C;
  double integral,array[1];

  if (recompute_Delta(C)) //strictly speaking, this is recomputed unnecessarily if only sigma_8 changes
  {
    integral = int_gsl_integrate_medium_precision(int_for_sigma_r_sqr,(void*)array,1e-4,1e6,NULL,512);
    res = 9.0*integral;   //see Peackock97, eq. 29
    update_cosmopara(&C);

  }
  if (!(res>0.0)){
    fprintf(stderr,"failed with sigma_r_sqr = %le\n", res);
  }
  assert(res>0.0);
  return res;
}

double Delta_L_wiggle(double k)
{
  static cosmopara C;

  static double *table_P;
  static double dk = .0, logkmin = .0, logkmax = .0;

  double klog,f1,norm;
  int i;

  if (k < limits.k_min_mpc || k > limits.k_max_mpc){
    norm=cosmology.sigma_8*cosmology.sigma_8/sigma_r_sqr();

    return norm*pow(k,cosmology.n_spec+ 0.5*cosmology.alpha_s*log(k/0.05)+3.0)*Tsqr_EH_wiggle(k);
    //printf("outside Delta_L_tab\n");
  }
  else{
    if (recompute_Delta(C))
    {
      if (cosmology.M_nu > 0){
        printf("Implementation of EH transfer function does not support massive neutrinos\n EXIT\n");
      }
      update_cosmopara(&C);
      norm=cosmology.sigma_8*cosmology.sigma_8/sigma_r_sqr();

      if(table_P!=0) free_double_vector(table_P,0, Ntable.N_k_lin-1);
      table_P=create_double_vector(0, Ntable.N_k_lin-1);

      logkmin = log(limits.k_min_mpc);
      logkmax = log(limits.k_max_mpc);
      dk = (logkmax - logkmin)/(Ntable.N_k_lin-1.);
      klog = logkmin;

      for (i=0; i<Ntable.N_k_lin; i++, klog += dk) {
       table_P[i]=log(norm*pow(exp(klog),cosmology.n_spec+ 0.5*cosmology.alpha_s*log(k/0.05)+3.0)*Tsqr_EH_wiggle(exp(klog)));
     }
      //printf("finished Delta_L_wiggle\n");
   }
 }
 klog=log(k);
 f1=interpol(table_P, Ntable.N_k_lin, logkmin, logkmax, dk,klog, 1.0,1.0 );
 return exp(f1);
}

// linear power spectrum routine with k in units H_0/c; used in covariances.c for beat coupling and in halo.c
double p_lin(double k,double a)
{
  static cosmopara C;
  static double **table_P_Lz = 0;
  static double logkmin = 0., logkmax = 0., dk = 0., da = 0.;
  int status;
  // if (strcmp(pdeltaparams.runmode,"CLASS")==0 || strcmp(pdeltaparams.runmode,"class")==0) return p_class(k,a,0, &status);

  double amp,ampsqr,grow0,aa,klog,val;

  int i,j;
  if (a >= 0.99999){a =0.99999;}
  if (recompute_cosmo3D(C)){
    update_cosmopara(&C);
    if (table_P_Lz!=0) free_double_matrix(table_P_Lz,0, Ntable.N_a-1, 0, Ntable.N_k_lin-1);
    table_P_Lz = create_double_matrix(0, Ntable.N_a-1, 0, Ntable.N_k_lin-1);
    grow0=growfac(1.);
    da = (1. - limits.a_min)/(Ntable.N_a-1.);
    aa = limits.a_min;
    for (i=0; i<Ntable.N_a; i++, aa +=da) {
      if(aa>1.0) aa=1.0;
      amp=growfac(aa)/grow0;
      ampsqr=amp*amp;

      logkmin = log(limits.k_min_mpc);
      logkmax = log(limits.k_max_mpc);
      dk = (logkmax - logkmin)/(Ntable.N_k_lin-1.);
      klog = logkmin;
      for (j=0; j<Ntable.N_k_lin; j++, klog += dk) {
        table_P_Lz[i][j] = log(ampsqr*Delta_L_wiggle(exp(klog)));
        //printf("%le %le",exp(klog),Delta_L_wiggle(exp(klog));
      }
    }
  }
  if (k/cosmology.coverH0 > exp(logkmax)) return 0.0;
  klog = log(k/cosmology.coverH0);
  val = interpol2d(table_P_Lz, Ntable.N_a, limits.a_min, 1., da, a, Ntable.N_k_lin, logkmin, logkmax, dk, klog, 3.0+cosmology.n_spec, 0.0);
  if(isnan(val) || (k==0)) return 0.0;
  return 2.0*constants.pi_sqr*exp(val)/k/k/k;
}

double int_sig_R_knl(double lnk, void *args) // tak12 A4
{
  double krsqr;

  double *params= (double *) args;
  double Rscale=params[0];
  //printf("Rscale %le k %le\n",Rscale,exp(lnk));
  krsqr= SQR(exp(lnk)*Rscale);
  return Delta_L_wiggle(exp(lnk))*exp(-krsqr);
}

double int_neff(double lnk, void *args) //tak12 A5
{
  double krsqr;
  double *params= (double *) args;
  double Rscale=params[0];
  krsqr= SQR(exp(lnk)*Rscale);
  return Delta_L_wiggle(exp(lnk))*2.0*krsqr*exp(-krsqr); //see S03 eq. 59
}

double int_cur(double lnk, void *args) //tak12 A5
{
  double krsqr;
  double *params= (double *) args;
  double Rscale=params[0];
  krsqr= SQR(exp(lnk)*Rscale);
  return Delta_L_wiggle(exp(lnk))*4.0*krsqr*(1.0-krsqr)*exp(-krsqr); // S03 eq.60
}

//iterative calculation of the nonlinear scale as defined in tak12 A4
void nonlin_scale(double amp, double *R_NL, double *neff, double *Curv)
{
  double sig_R,kmax,logkmax,sig_R_noamp,neffplus3;
  int iterstep;
  const int itermax  = 40;
  int converged=0;
  double array[1];
  double logRmin = -3.0;
  double logRmax =  4.0;

  iterstep=0;
  while(converged==0)
  {
    array[0]=pow(10.,(logRmin+logRmax)/2.0);

    //flexible upper limit of integration depending on drop-off of filter function
    kmax  = sqrt(5.*log(10.))/array[0];
    if (kmax<8000.0) logkmax = log(8000.0);
    sig_R_noamp=sqrt(int_gsl_integrate_medium_precision(int_sig_R_knl,(void*)array,-4.5,logkmax,NULL,512)); //integral goes over ln k exponent correspond to k_min~0.011

    sig_R=amp*sig_R_noamp;
    if (sig_R>1.0)  logRmin=log10(array[0]);
    if (sig_R<1.0)  logRmax=log10(array[0]);
    iterstep=iterstep+1;
    if(fabs(sig_R-1.0) < 0.0001 || iterstep>itermax) converged=1;
  }
  *R_NL=array[0]; //R where sig_R==1
  neffplus3=int_gsl_integrate_medium_precision(int_neff,(void*)array,-4.5,logkmax,NULL,512)/sig_R_noamp/sig_R_noamp;
  *neff= neffplus3 - 3.0;
  *Curv= int_gsl_integrate_medium_precision(int_cur,(void*)array,-4.5,logkmax,NULL,512)/sig_R_noamp/sig_R_noamp + SQR(neffplus3);

  //printf("%d %le\n",iterstep,amp);
}

double Halofit(double k, double amp, double omm, double omv,double w_z, double R_NL, double neff,double Curv, double P_delta_Lin)
{
  double y_scale,n2eff,n3eff,n4eff;
  double a_n,b_n,c_n,gamma_n,alpha_n,beta_n,nu_n,f1,f2,f3;

  double Delta_H,Delta_H_Prime,Delta_Q;
  //determine nonlinear scale, neff and curvature, see tak12 A4, A5
  y_scale=k*R_NL;

  n2eff=neff*neff;
  n3eff=n2eff*neff;
  n4eff=n2eff*n2eff;

  //calculate coefficients
  a_n = pow(10.,1.5222+2.8553*neff + 2.3706*n2eff+0.9903*n3eff+0.2250*n4eff-0.6038*Curv+0.1749*omv*(1.0+w_z));
  b_n = pow(10., -0.5642+0.5864*neff + 0.5716*n2eff-1.5474*Curv +0.2279*omv*(1.0+w_z));
  c_n = pow(10., 0.3698+ 2.0404*neff + 0.8161*n2eff+0.5869*Curv);
  gamma_n = 0.1971-0.0843*neff + 0.8460*Curv;
  alpha_n = fabs(6.0835 + 1.3373*neff - 0.1959*n2eff - 5.5274*Curv);
  beta_n = 2.0379 - 0.7354*neff + 0.3157*n2eff + 1.2490*n3eff + 0.3980*n4eff - 0.1682*Curv;
  nu_n = pow(10,5.2105+3.6902*neff);

  f1 = pow(omm,(-0.0307));
  f2 = pow(omm,(-0.0585));
  f3 = pow(omm,(0.0743));

  //TwoHaloTerm
  Delta_Q=P_delta_Lin*(pow((1.0+P_delta_Lin),beta_n)/(1.0+alpha_n*P_delta_Lin))*exp(-(y_scale/4.0+y_scale*y_scale/8.0));
  //OneHaloterm
  Delta_H_Prime=(a_n*pow(y_scale,3.0*f1))/(1.0+b_n*pow(y_scale,f2)+pow(c_n*f3*y_scale,3.0-gamma_n));
  Delta_H=Delta_H_Prime/(1.0+nu_n*pow(y_scale,-2.0)); // using mu=0.0 Tak A12
  //printf("Delta_Q %le Delta_H %le\n",Delta_Q,Delta_H);
  return Delta_H+Delta_Q;
}

void Delta_halofit(double **table_P_NL,double logkmin, double logkmax, double dk, double da)
{
  double rk,omm,omv,w_z,amp,grow0,aa,klog;
  double R_NL,Curv,neff,P_delta,P_delta_Lin;
  int i,j;

  grow0=growfac(1.);
  aa = limits.a_min;
  //binning in k and a must be the same as in emu
  for (i=0; i<Ntable.N_a; i++, aa +=da) {
    if(aa>1.0) aa=1.0;
    omega_a(aa,&omm,&omv);
    w_z=cosmology.w0+cosmology.wa*(1.-aa);
    amp=growfac(aa)/grow0;
    nonlin_scale(amp, &R_NL, &neff, &Curv);
    //printf("%le %le %le %le\n",aa,R_NL,neff,Curv);
    klog = logkmin;
    for (j=0; j<Ntable.N_k_nlin; j++, klog += dk) {
      rk=exp(klog);
      P_delta_Lin=amp*amp*Delta_L_wiggle(rk);
      P_delta=Halofit(rk, amp, omm, omv, w_z, R_NL, neff, Curv, P_delta_Lin);
      table_P_NL[i][j]=log(P_delta);
    }
  }
}

double Delta_NL_Halofit(double k_NL, double a)
{
  static cosmopara C;
  static double logkmin = 0., logkmax = 0., dk = 0., da = 0.;

  static double **table_P_NL=0;
  double klog,val;

  if (recompute_cosmo3D(C)){
    update_cosmopara(&C);
    if (table_P_NL!=0) free_double_matrix(table_P_NL,0, Ntable.N_a-1, 0, Ntable.N_k_nlin-1);
    table_P_NL = create_double_matrix(0, Ntable.N_a-1, 0, Ntable.N_k_nlin-1);

    da = (1. - limits.a_min)/(Ntable.N_a-1.);
    logkmin = log(limits.k_min_mpc);
    logkmax = log(limits.k_max_mpc);
    dk = (logkmax - logkmin)/(Ntable.N_k_nlin-1.);

    Delta_halofit(table_P_NL,logkmin, logkmax, dk, da);
  }
  klog = log(k_NL);
  val = interpol2d(table_P_NL, Ntable.N_a, limits.a_min, 1., da, a, Ntable.N_k_nlin, logkmin, logkmax, dk, klog, cosmology.n_spec, 0.0);
  return exp(val);
  // returns the dimensionless power spectrum as a function of scale factor a and k
}

double nonlinear_scale_computation(double a)
{
  static cosmopara C;
  static double da = 0.;
  static double *table=0;

  double omm,omv,amp,grow0,aa,res;
  double R_NL,Curv,neff;
  int i;

  if (recompute_cosmo3D(C)){
    update_cosmopara(&C);
    grow0=growfac(1.);
    da = (1.-limits.a_min)/(Ntable.N_a-1.);
    aa = limits.a_min;
    if (table!=0) free_double_vector(table, 0, Ntable.N_a-1);
    table   = create_double_vector(0, Ntable.N_a-1);
    for (i=0; i<Ntable.N_a; i++, aa+=da) {
      if(aa>1.0) aa=1.0;
      omega_a(aa,&omm,&omv);
      amp=growfac(aa)/grow0;
      nonlin_scale(amp, &R_NL, &neff, &Curv);
      table[i] = 1./R_NL;
    }
  }
  res = interpol(table, Ntable.N_a, limits.a_min, 1., da, a, 0.0, 0.0);
  return res;
}

/*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*/
//Pdelta is called with k in units H0/c since the comoving distance chi is in units c/H0. Upstream Pdelta all routines are in h/mpc
double Pdelta(double k_NL,double a)
{
  static int P_type = -1;
  if (P_type == -1){
    if (strcmp(pdeltaparams.runmode,"Halofit")==0) P_type = 0;
    if (strcmp(pdeltaparams.runmode,"halofit")==0) P_type = 0;
    // if (strcmp(pdeltaparams.runmode,"emu")==0) P_type = 1;
    // if (strcmp(pdeltaparams.runmode,"emu_only")==0) P_type = 2;
    if (strcmp(pdeltaparams.runmode,"linear")==0) P_type = 3;
    // if (strcmp(pdeltaparams.runmode,"CLASS")==0) P_type = 4;
    // if (strcmp(pdeltaparams.runmode,"class")==0) P_type = 4;
    if (strcmp(pdeltaparams.runmode,"cosmo_sim_test") ==0) P_type = 5;

  }

  //printf("%s set\n",pdeltaparams.runmode);
  double pdelta = 0.,kintern=k_NL/cosmology.coverH0,error,k_nonlin,res;
  int status;
  switch (P_type){
    case 0: pdelta=2.0*constants.pi_sqr*Delta_NL_Halofit(kintern,a)/k_NL/k_NL/k_NL; break;
    // case 1: pdelta=2.0*constants.pi_sqr*Delta_NL_emu(kintern,a)/k_NL/k_NL/k_NL; break;
    // case 2: pdelta=2.0*constants.pi_sqr*Delta_NL_emu_only(kintern,a)/k_NL/k_NL/k_NL; break;
    case 3: pdelta=p_lin(k_NL,a); break;
    // case 4: pdelta=p_class(k_NL,a,1, &status); break;
    case 5: k_nonlin=nonlinear_scale_computation(a);
    if (kintern<0.01) pdelta=2.0*constants.pi_sqr*Delta_NL_Halofit(kintern,a)/k_NL/k_NL/k_NL;
    else{
      error=0.01*pow((pdeltaparams.DIFF_A*kintern/k_nonlin),pdeltaparams.DIFF_n);
      pdelta=2.0*constants.pi_sqr*Delta_NL_Halofit(kintern,a)*(1.0+error)/k_NL/k_NL/k_NL;
    }
    break;
    default:
    printf("cosmo3D:Pdelta: %s Pdelta runmode not defined\n",pdeltaparams.runmode);
    printf("using Halofit (standard)\n");
    pdelta=2.0*constants.pi_sqr*Delta_NL_Halofit(kintern,a)/k_NL/k_NL/k_NL;
    break;
  }

//      double z = 1./a -1 ;
//	  if (z > 4.) {
//		 printf("	z:%lf    a:%lf   	k_Mpc:%lf    PkRatio:%lf \n",z,a,kintern,PkRatio_baryons(kintern, a));
//	  }

  // if (bary.isPkbary==1) pdelta = pdelta*PkRatio_baryons(kintern, a);
  return pdelta;
}

/*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*/

/*============================================================
 *see BS 2.41 bzw Logbook for detailed calculation of chi from a.*/

double int_for_chi(double a, void * args){
  //double res,asqr;
  //asqr=a*a;
  //res= 1./sqrt(a*cosmology.Omega_m + asqr*(1.-cosmology.Omega_m -cosmology.Omega_v ) + asqr*asqr*omv_vareos(a));
  //return res;
  return 1./(a*a*hoverh0(a)); //changed to call of hoverh0 to be ready for other parametrizations
}

/*for the calculation of chi we have to integrate from a(z2)=a up to a(z1)=1, which means todays expansion factor*/
double chi(double a)
{
  static cosmopara C;
  static double *table;
  static double da = 0.;
  double aa,res;
  int i;
  double array[1];

  if (recompute_expansion(C)){
    update_cosmopara(&C);
    da = (1.-limits.a_min)/(Ntable.N_a-1.);
    aa = limits.a_min;
    if (table!=0) free_double_vector(table, 0, Ntable.N_a-1);
    table   = create_double_vector(0, Ntable.N_a-1);
    for (i=0; i<Ntable.N_a-1; i++, aa+=da) {
      table[i] = int_gsl_integrate_medium_precision(int_for_chi,(void*)array, aa, 1.,NULL,1000);
    }
    table[Ntable.N_a-1] =.0;
  }
  res = interpol(table, Ntable.N_a, limits.a_min, 1., da, a, 0.0, 0.0); // comoving distance in c/H_0
  if (res < 0){printf ("interpolation error in chi(%e)\n",a); res=0.01;}
  return res;
}
//auxilary function to look up a(chi)
double a_chi(double chi1){
  static gsl_spline * a_spline = NULL;
  static gsl_interp_accel * a_accel = NULL;
  static cosmopara C;
  static double chi_max =-1.;
  if (!a_spline){
    a_spline = gsl_spline_alloc(gsl_interp_cspline, Ntable.N_a);
    a_accel = gsl_interp_accel_alloc();
  }
  if (recompute_cosmo3D(C)){
    update_cosmopara(&C);
    double *table_a,*table_chi;
    table_a  = create_double_vector(0, Ntable.N_a-1);
    table_chi  = create_double_vector(0, Ntable.N_a-1);
    for (int i = 0; i < Ntable.N_a; i++){
      table_a[i] = 1.0 - (1.0 - 0.99*limits.a_min)/(Ntable.N_a-1.)*(double)i;
      table_chi[i] = int_gsl_integrate_medium_precision(int_for_chi,NULL, table_a[i], 1.,NULL,1000);
   //   printf("%d %e %e\n",i,table_a[i],table_chi[i]);
    }
    chi_max = int_gsl_integrate_medium_precision(int_for_chi,NULL, limits.a_min, 1.,NULL,1000);
    gsl_spline_init(a_spline, table_chi, table_a, Ntable.N_a);
    free_double_vector(table_a,0, Ntable.N_a-1);
    free_double_vector(table_chi,0, Ntable.N_a-1);
  }
  if (chi1 <=0.0){return 1.0;}
  if (chi1 > chi_max){printf("called a_chi(chi) with chi > chi(limits.a_min\nEXIT\n");exit(1);}
  return gsl_spline_eval(a_spline,chi1,a_accel);
}

/*===============================calculating the angular diameter distance f_K BS01 2.4, 2.30: f_K is a radial function that, depending on the curvature of the Universe, is a trigonometric, linear, or hyperbolic function of chi  */
double f_K(double chi)
{
  double K, K_h, f;
  K = (cosmology.Omega_m   + cosmology.Omega_v  - 1.);
  if (K > precision.medium) {           /* open */
    K_h = sqrt(K); // K in units H0/c see BS eq. 2.30
    f = 1./K_h*sin(K_h*chi);
    //printf("open\n");
  } else if (K < -precision.medium) {   /* closed */
    K_h = sqrt(-K);
    f = 1./K_h*sinh(K_h*chi);
    //printf("closed K=%le %le %le\n",K,cosmology.Omega_m,cosmology.Omega_v);
  } else {                     /* flat */
    f = chi;
    //printf("flatK=%le %le %le\n",K,cosmology.Omega_m,cosmology.Omega_v);
  }
  return f;
}

max_stars_repo_name: joezuntz/CosmoCov

max_issues_repo_name: joezuntz/CosmoCov

max_forks_repo_name: joezuntz/CosmoCov


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

#include <stdlib.h>

#include <gsl/gsl_errno.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_multifit.h>
#include <gsl/gsl_blas.h>

#include "gsl_multifit_ndlinear.h"

static int ndlinear_construct_row(const gsl_vector *d, gsl_vector *x,
                                  gsl_multifit_ndlinear_workspace *w);

/*
gsl_multifit_ndlinear_alloc()
  Allocate a ndlinear workspace

Inputs: n_dim  - dimension of fit function
        N      - number of terms in each sum; N[i] = N_i, 0 <= i < n_dim
        u      - basis functions to call
                 u[j] = u^{(j)}, 0 <= j < n_dim
        params - parameters to pass to basis functions

Return: pointer to new workspace

Notes: the supplied basis functions 'u[j]' must accept three
       arguments:

int uj(double x, double y[], void *params)

and fill the y[] vector so that y[i] = u_{i}^{(j)}(x) (the ith
basis function for the jth parameter evaluated at x)
*/

gsl_multifit_ndlinear_workspace *
gsl_multifit_ndlinear_alloc(size_t n_dim, size_t N[],
                            int (**u)(double x, double y[], void *p),
                            void *params)
{
  gsl_multifit_ndlinear_workspace *w;
  size_t n_coeffs; /* total number of fit coefficients */
  size_t i, idx;
  size_t sum_N;

  if (n_dim == 0)
    {
      GSL_ERROR_NULL("n_dim must be at least 1", GSL_EINVAL);
    }

  w = calloc(1, sizeof(gsl_multifit_ndlinear_workspace));
  if (!w)
    {
      GSL_ERROR_NULL("failed to allocate space for workspace", GSL_ENOMEM);
    }

  w->N = calloc(n_dim, sizeof(size_t));
  if (!w->N)
    {
      gsl_multifit_ndlinear_free(w);
      GSL_ERROR_NULL("failed to allocate space for N vector", GSL_ENOMEM);
    }

  n_coeffs = 1;
  sum_N = 0;
  for (i = 0; i < n_dim; ++i)
    {
      if (N[i] == 0)
        {
          gsl_multifit_ndlinear_free(w);
          GSL_ERROR_NULL("one of the sums is empty", GSL_EINVAL);
        }

      /* The total number of coefficients is: N_1 * N_2 * ... * N_n */
      n_coeffs *= N[i];
      w->N[i] = N[i];
      sum_N += N[i];
    }

  w->n_dim = n_dim;
  w->n_coeffs = n_coeffs;

  w->work = gsl_vector_alloc(n_coeffs);
  w->work2 = gsl_vector_alloc(sum_N);
  if (!w->work || !w->work2)
    {
      gsl_multifit_ndlinear_free(w);
      GSL_ERROR_NULL("failed to allocate space for basis vector",
                     GSL_ENOMEM);
    }

  w->v = calloc(n_dim, sizeof(gsl_vector_view));
  if (!w->v)
    {
      gsl_multifit_ndlinear_free(w);
      GSL_ERROR_NULL("failed to allocate space for basis vector",
                     GSL_ENOMEM);
    }

  w->u = calloc(n_dim, sizeof(int *));
  if (!w->u)
    {
      gsl_multifit_ndlinear_free(w);
      GSL_ERROR_NULL("failed to allocate space for basis functions",
                     GSL_ENOMEM);
    }

  idx = 0;
  for (i = 0; i < n_dim; ++i)
    {
      w->v[i] = gsl_vector_subvector(w->work2, idx, N[i]);
      idx += N[i];

      w->u[i] = u[i];
    }

  w->params = params;

  return (w);
} /* gsl_multifit_ndlinear_alloc() */

/*
gsl_multifit_ndlinear_free()
  Free workspace w
*/

void
gsl_multifit_ndlinear_free(gsl_multifit_ndlinear_workspace *w)
{
  if (w->N)
    free(w->N);

  if (w->work)
    gsl_vector_free(w->work);

  if (w->work2)
    gsl_vector_free(w->work2);

  if (w->v)
    free(w->v);

  if (w->u)
    free(w->u);

  free(w);
} /* gsl_multifit_ndlinear_free() */

/*
gsl_multifit_ndlinear_design()
  This function constructs the coefficient design matrix 'X'

Inputs: vars  - independent variable vectors for matrix X
                vars is a ndata-by-n_dim matrix where the ith row
                specifies the n_dim independent variables for the
                ith observation, so that
                vars_{ij} = (x_i)_j, the jth element of the
                ith input variable vector
        X     - (output) design matrix (must be ndata-by-w->n_coeffs)
        w     - workspace

Return: success or error
*/

int
gsl_multifit_ndlinear_design(const gsl_matrix *vars, gsl_matrix *X,
                             gsl_multifit_ndlinear_workspace *w)
{
  const size_t ndata = vars->size1;

  if ((X->size1 != ndata) || (X->size2 != w->n_coeffs))
    {
      GSL_ERROR("X matrix has wrong dimensions", GSL_EBADLEN);
    }
  else
    {
      size_t i; /* looping */
      int s;

      for (i = 0; i < ndata; ++i)
        {
          gsl_vector_const_view d = gsl_matrix_const_row(vars, i);
          gsl_vector_view xv = gsl_matrix_row(X, i);

          s = ndlinear_construct_row(&d.vector, &xv.vector, w);
          if (s != GSL_SUCCESS)
            return s;
        }

      return GSL_SUCCESS;
    }
} /* gsl_multifit_ndlinear_design() */

/*
gsl_multifit_ndlinear_est()
  Compute the model function at a given data point with errors

Inputs: x     - data point (w->n_dim elements)
        c     - coefficient vector
        cov   - covariance matrix
        y     - where to store fit function result
        y_err - standard deviation of fit
        w     - workspace

Return: success or error
*/

int
gsl_multifit_ndlinear_est(const gsl_vector *x, const gsl_vector *c,
                          const gsl_matrix *cov, double *y, double *y_err,
                          gsl_multifit_ndlinear_workspace *w)
{
  if (c->size != w->n_coeffs)
    {
      GSL_ERROR("c vector has wrong size", GSL_EBADLEN);
    }
  else
    {
      int s;

      s = ndlinear_construct_row(x, w->work, w);
      if (s != GSL_SUCCESS)
        return s;

      /*
       * Now w->work contains the appropriate basis functions
       * evaluated at the given point - compute the function value
       */
      s = gsl_multifit_linear_est(w->work, c, cov, y, y_err);

      return s;
    }
} /* gsl_multifit_ndlinear_est() */

/*
gsl_multifit_ndlinear_calc()
  Compute the model function at a given data point

Inputs: x - data point (w->n_dim elements)
        c - coefficient vector
        w - workspace

Return: model value
*/

double
gsl_multifit_ndlinear_calc(const gsl_vector *x, const gsl_vector *c,
                           gsl_multifit_ndlinear_workspace *w)
{
  if (c->size != w->n_coeffs)
    {
      GSL_ERROR_VAL("c vector has wrong size", GSL_EBADLEN, 0.0);
    }
  else
    {
      double y;
      int s;

      s = ndlinear_construct_row(x, w->work, w);
      if (s != GSL_SUCCESS)
        {
          GSL_ERROR_VAL("constructing matrix row failed", s, 0.0);
        }

      gsl_blas_ddot(w->work, c, &y);

      return y;
    }
} /* gsl_multifit_ndlinear_calc() */

/*
gsl_multifit_ndlinear_ncoeffs()
  Return the total number of fit coefficients
*/

size_t
gsl_multifit_ndlinear_ncoeffs(gsl_multifit_ndlinear_workspace *w)
{
  return w->n_coeffs;
} /* gsl_multifit_ndlinear_ncoeffs() */

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

/*
ndlinear_construct_row()

  Compute a row of the design matrix X:

X(:,j) = u_{r_0}^{(0)}(d_0) * u_{r_1)^{(1)}(d_1) * ... *
         u_{r_{n-1}}^{(n-1)}(d_{n-1})

where 'd' is the corresponding data vector for that row

Inputs: d - data vector of length w->n_dim
        x - (output) where to store row of design matrix X
        w - workspace

Return: success or error
*/

static int
ndlinear_construct_row(const gsl_vector *d, gsl_vector *x,
                       gsl_multifit_ndlinear_workspace *w)
{
  size_t j;
  int k, s;
  size_t denom, rk;
  double melement;

  /* compute basis functions for this data point */
  for (j = 0; j < w->n_dim; ++j)
    {
      s = w->u[j](gsl_vector_get(d, j), w->v[j].vector.data, w->params);

      if (s != GSL_SUCCESS)
        return s;
    }

  for (j = 0; j < w->n_coeffs; ++j)
    {
      /*
       * The (:,j) element of the matrix X will be:
       *
       * X_{:,j} = u_{r_0}^{(0)}(d_0) *
       *           u_{r_1)^{(1)}(d_1) *
       *           ... *
       *           u_{r_{n-1}}^{(n-1)}(d_{n-1})
       *
       * with the basis function indices r_k given by
       *
       * r_k = floor(j / Prod_{i=(k+1)..(n-1)} [ N_i ]) (mod N_k)
       *
       * In the case where N_i = N for all i,
       *
       * r_k = floor(j / N^{n - k - 1}) (mod N)
       *
       * n: dimension of fit function (w->n_dim)
       * N_i: number of terms in sum i of fit function
       */

      /* calculate the r_k and the matrix element X_{:,j} */

      denom = 1;
      melement = 1.0;
      for (k = (int)(w->n_dim - 1); k >= 0; --k)
        {
          rk = (j / denom) % w->N[k];
          denom *= w->N[k];

          melement *= gsl_vector_get(&(w->v[k]).vector, rk);
        }

      /* set the matrix element */
      gsl_vector_set(x, j, melement);
    }

  return GSL_SUCCESS;
} /* ndlinear_construct_row() */

max_stars_repo_name: rennhak/Keyposes

max_issues_repo_name: rennhak/Keyposes

max_forks_repo_name: rennhak/Keyposes


text:
#ifndef structs_h
#define structs_h

#include <ceed.h>
#include <petsc.h>
#include "../problems/cl-problems.h"

// -----------------------------------------------------------------------------
// Command Line Options
// -----------------------------------------------------------------------------
// Forcing function options
typedef enum {
  FORCE_NONE = 0, FORCE_CONST = 1, FORCE_MMS = 2
} forcingType;
static const char *const forcing_types[] = {"none",
                                            "constant",
                                            "mms",
                                            "forcingType","FORCE_",0
                                           };
static const char *const forcing_types_for_disp[] = {"None",
                                                     "Constant",
                                                     "Manufactured solution"
                                                    };

// Multigrid options
typedef enum {
  MULTIGRID_LOGARITHMIC = 0, MULTIGRID_UNIFORM = 1, MULTIGRID_NONE = 2
} multigridType;
static const char *const multigrid_types [] = {"logarithmic",
                                               "uniform",
                                               "none",
                                               "multigridType","MULTIGRID",0
                                              };
static const char *const multigrid_types_for_disp[] = {"P-multigrid, logarithmic coarsening",
                                                       "P-multigrind, uniform coarsening",
                                                       "No multigrid"
                                                      };

// -----------------------------------------------------------------------------
// Application data structs
// -----------------------------------------------------------------------------
// Units
typedef struct Units_ *Units;
struct Units_ {
  // Fundamental units
  PetscScalar meter;
  PetscScalar kilogram;
  PetscScalar second;
  // Derived unit
  PetscScalar Pascal;
};

// Application context from user command line options
typedef struct AppCtx_ *AppCtx;
struct AppCtx_ {
  const char    *name, *name_for_disp;                 // problem name
  char          ceed_resource[PETSC_MAX_PATH_LEN];     // libCEED backend
  char          mesh_file[PETSC_MAX_PATH_LEN];         // exodusII mesh file
  char          output_dir[PETSC_MAX_PATH_LEN];
  PetscBool     test_mode;
  PetscBool     view_soln;
  PetscBool     view_final_soln;
  PetscViewer   energy_viewer;
  problemType   problem_choice;
  forcingType   forcing_choice;
  multigridType multigrid_choice;
  PetscInt      degree;
  PetscInt      q_extra;
  PetscInt      num_levels;
  PetscInt      *level_degrees;
  PetscInt      num_increments;                        // Number of steps
  PetscInt      bc_clamp_count;
  PetscInt      bc_clamp_faces[16];
  // [translation; 3] [rotation axis; 3] [rotation magnitude c_0, c_1]
  // The rotations are (c_0 + c_1 s) \pi, where s = x · axis
  PetscScalar   bc_clamp_max[16][8];
  PetscInt      bc_traction_count;
  PetscInt      bc_traction_faces[16];
  PetscScalar   bc_traction_vector[16][3];
  PetscScalar   forcing_vector[3];
  PetscReal     test_tol;
  PetscReal     expect_final_strain;
};

// Forcing function data
typedef struct {
  CeedQFunctionUser setup_forcing;
  const char        *setup_forcing_loc;
} forcingData;

extern forcingData forcing_options[3];

// Data for PETSc Matshell
typedef struct UserMult_ *UserMult;
struct UserMult_ {
  MPI_Comm        comm;
  DM              dm;
  Vec             X_loc, Y_loc, neumann_bcs;
  CeedVector      x_ceed, y_ceed;
  CeedOperator    op;
  CeedQFunction   qf;
  Ceed            ceed;
  PetscScalar     load_increment;
  CeedQFunctionContext ctx_phys, ctx_phys_smoother;
};

// Data for Jacobian setup routine
typedef struct FormJacobCtx_ *FormJacobCtx;
struct FormJacobCtx_ {
  UserMult     *jacob_ctx;
  PetscInt     num_levels;
  SNES         snes_coarse;
  Mat          *jacob_mat, jacob_mat_coarse;
  Vec          u_coarse;
};

// Data for PETSc Prolongation/Restriction Matshell
typedef struct UserMultProlongRestr_ *UserMultProlongRestr;
struct UserMultProlongRestr_ {
  MPI_Comm     comm;
  DM           dm_c, dm_f;
  Vec          loc_vec_c, loc_vec_f;
  CeedVector   ceed_vec_c, ceed_vec_f;
  CeedOperator op_prolong, op_restrict;
  Ceed         ceed;
};

#define SOLIDS_MAX_NUMBER_FIELDS 16

// libCEED data struct for level
typedef struct CeedData_ *CeedData;
struct CeedData_ {
  Ceed                ceed;
  CeedBasis           basis_x, basis_u, basis_c_to_f, basis_energy,
                      basis_diagnostic;
  CeedElemRestriction elem_restr_x, elem_restr_u, elem_restr_geo_data_i,
                      elem_restr_energy, elem_restr_diagnostic,
                      elem_restr_geo_data_diagnostic_i,
                      elem_restr_stored_fields_i[SOLIDS_MAX_NUMBER_FIELDS];
  CeedQFunction       qf_residual, qf_jacobian, qf_energy, qf_diagnostic;
  CeedOperator        op_residual, op_jacobian, op_restrict, op_prolong,
                      op_energy,
                      op_diagnostic;
  CeedVector          geo_data, geo_data_diagnostic, x_ceed, y_ceed,
                      true_soln, stored_fields[SOLIDS_MAX_NUMBER_FIELDS];
};

typedef struct {
  CeedQFunctionUser setup_geo, residual, jacobian, energy,
                    diagnostic, true_soln;
  const char *setup_geo_loc, *residual_loc, *jacobian_loc, *energy_loc,
        *diagnostic_loc, *true_soln_loc;
  CeedQuadMode quadrature_mode;
  CeedInt geo_data_size, number_fields_stored;
  CeedInt *field_sizes;
  const char *const *field_names;
} ProblemData;

#endif // structs_h

max_stars_repo_name: LeilaGhaffari/libCEED

max_issues_repo_name: LeilaGhaffari/libCEED

max_forks_repo_name: LeilaGhaffari/libCEED


text:
#pragma once

#include <memory>
#include <string>
#include <vector>

#include <gsl/span>

struct VfsSearchResult {
	std::string filename;
	bool dir;
	size_t sizeInBytes;
	uint32_t lastModified;
};

enum class SeekDir {
	Start = 0,
	Current = 1,
	End = 2
};

/*
	Abstractions for accessing files in the virtual file system.
	The virtual file system is backed by TIO in a normal game,
	which means that the files are scattered across the .dat files
	and data directories.
*/
class Vfs {
public:
	virtual ~Vfs() {
	};

	static Vfs* CreateStdIoVfs();

	/*
		Reads a file fully into a string.
	*/
	std::string ReadAsString(std::string_view filename);

	/*
	Reads a binary file fully into a vector of uint8_t.
	*/
	std::vector<uint8_t> ReadAsBinary(std::string_view filename);

	/**
	 * Does the file exist?
	 */
	virtual bool FileExists(std::string_view path) = 0;

	/**
	 * Does the directory exist?
	 */
	virtual bool DirExists(std::string_view path) = 0;

	/**
	 * Creates a directory.
	 */
	virtual bool MkDir(std::string_view path) = 0;

	/**
	 * Will return all files and directories that match the given glob pattern, which
	 * is not recursive.
	 * Example: mes\*.mes will return the names of all files in mes (not full paths) that
	            end with .mes.
	 */
	virtual std::vector<VfsSearchResult> Search(std::string_view globPattern) = 0;

	/**
	 * Removes an empty directory.
	 */
	virtual bool RemoveDir(std::string_view path) = 0;

	/**
	 * Removes a file (no directories).
	 */
	virtual bool RemoveFile(std::string_view path) = 0;

	/**
	* Deletes all files and dirctories within the given directory.
	*/
	bool CleanDir(std::string_view path);

	/**
	 * Returns true if the given directory does not contain anything.
	 */
	bool IsDirEmpty(std::string_view path);

	/**
	 * Writes binary data to a file.
	 */
	void WriteBinaryFile(std::string_view path, gsl::span<uint8_t> data);

	using FileHandle = void*;
	virtual FileHandle Open(std::string_view name, std::string_view mode) = 0;
	virtual size_t Read(void* buffer, size_t size, FileHandle handle) = 0;
	virtual size_t Write(const void* buffer, size_t size, FileHandle handle) = 0;
	virtual size_t Length(FileHandle handle) = 0;
	virtual size_t Tell(FileHandle handle) = 0;
	virtual void Seek(FileHandle handle, int position, SeekDir dir = SeekDir::Start) = 0;
	virtual void Close(FileHandle handle) = 0;

};

class VfsPath {
public:
	VfsPath() = delete;

	static bool IsFileSystem(std::string_view path);
	static std::string Concat(std::string_view a, std::string_view b);
};

extern std::unique_ptr<Vfs> vfs;

max_stars_repo_name: edoipi/TemplePlus

max_issues_repo_name: edoipi/TemplePlus

max_forks_repo_name: edoipi/TemplePlus


text:
/*
 * Copyright 2016 Maikel Nadolski
 *
 * 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.
 */

#include <iostream>
#include <gsl_assert.h>

#ifndef HMM_ALGORITHM_BACKWARD_H_
#define HMM_ALGORITHM_BACKWARD_H_

namespace maikel { namespace hmm {

template <class I, class J, class T>
    class backward_range_fn {
      public:
        using model       = hidden_markov_model<T>;
        using row_vector  = typename model::row_vector;
        using matrix      = typename model::matrix;
        using symbol_type = typename std::iterator_traits<I>::value_type;

        backward_range_fn() = delete;

        backward_range_fn(I seq_it, I seq_end, J scaling_it, model const& hmm)
        :  hmm_{&hmm}, seq_it_{seq_it}, seq_end_{seq_end}, scaling_it_{scaling_it},
          beta_(hmm.states()), next_beta_(hmm.states())
        {
          if (seq_it != seq_end)
            initial_coefficients(*scaling_it);
        }

        class iterator
        : public boost::iterator_facade<
              iterator, row_vector, std::input_iterator_tag, row_vector const&
         > {
          public:
            iterator() = default;
          private:
            friend class backward_range_fn;
            friend class boost::iterator_core_access;

            iterator(backward_range_fn& parent)
            : parent_{parent ? &parent : nullptr} {}

            backward_range_fn* parent_ = nullptr;

            row_vector const& dereference() const
            {
              Expects(parent_ && *parent_);
              return parent_->beta_;
            }

            void increment()
            {
              Expects(parent_ && *parent_);
              if (!parent_->next())
                parent_ = nullptr;
            }

            bool equal(iterator other) const noexcept
            {
              return parent_ == other.parent_;
            }
        };

        operator bool() const noexcept
        {
          return seq_it_ != seq_end_;
        }

        iterator begin() noexcept
        {
          return {*this};
        }

        iterator end() noexcept
        {
          return {};
        }

      private:
        model const* hmm_; // not owning
        I seq_it_, seq_end_;
        J scaling_it_;
        row_vector beta_;
        row_vector next_beta_;

        void initial_coefficients(T scaling) noexcept
        {
          Expects(beta_.size() == hmm_->states());
          beta_.fill(scaling);
        }

        void recursion_advance(symbol_type s, T scaling) noexcept
        {
          using size_type = typename model::size_type;
          matrix const& A = hmm_->transition_matrix();
          matrix const& B = hmm_->symbol_probabilities();
          next_beta_.swap(beta_);

          // check pre conditions
          size_type states = A.rows();
          Expects(A.cols() == states);
          Expects(B.rows() == states);
          Expects(next_beta_.size() == states);
          Expects(beta_.size() == states);
          size_type ob = gsl::narrow<size_type>(s);
          Expects(0 <= ob && ob < B.cols());

          // recursion formula
          for (size_type i = 0; i < states; ++i) {
            beta_(i) = 0.0;
            for (size_type j = 0; j < states; ++j)
              beta_(i) += A(i,j)*B(j,ob)*next_beta_(j);
            beta_(i) *= scaling;
          }
        }

        bool next()
        {
          Expects(*this);
          symbol_type s = *seq_it_++;
          ++scaling_it_;
          if (seq_it_ == seq_end_)
            return false;
          recursion_advance(s, *scaling_it_);
          return true;
        }
    };

  template <class I, class J, class T>
    backward_range_fn<I, J, T>
    backward(I begin, I end, J scaling, hidden_markov_model<T> const& hmm)
    {
      return {begin, end, scaling, hmm};
    }

} // namespace hmm
} // namespace maikel

#endif /* HMM_ALGORITHM_BACKWARD_H_ */

max_stars_repo_name: maikel/hidden-markov-model

max_issues_repo_name: maikel/Hidden-Markov-Model

max_forks_repo_name: maikel/Hidden-Markov-Model