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

#ifndef __GSL_POLY_H__
#define __GSL_POLY_H__

#include <stdlib.h>
#include <gsl/gsl_inline.h>
#include <gsl/gsl_complex.h>

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

__BEGIN_DECLS

/* Evaluate polynomial
 *
 * c[0] + c[1] x + c[2] x^2 + ... + c[len-1] x^(len-1)
 *
 * exceptions: none
 */

/* real polynomial, real x */
INLINE_DECL double gsl_poly_eval(const double c[], const int len, const double x);

/* real polynomial, complex x */
INLINE_DECL gsl_complex gsl_poly_complex_eval (const double c [], const int len, const gsl_complex z);

/* complex polynomial, complex x */
INLINE_DECL gsl_complex gsl_complex_poly_complex_eval (const gsl_complex c [], const int len, const gsl_complex z);

int gsl_poly_eval_derivs(const double c[], const size_t lenc, const double x, double res[], const size_t lenres);

#ifdef HAVE_INLINE
INLINE_FUN
double
gsl_poly_eval(const double c[], const int len, const double x)
{
  int i;
  double ans = c[len-1];
  for(i=len-1; i>0; i--) ans = c[i-1] + x * ans;
  return ans;
}

INLINE_FUN
gsl_complex
gsl_poly_complex_eval(const double c[], const int len, const gsl_complex z)
{
  int i;
  gsl_complex ans;
  GSL_SET_COMPLEX (&ans, c[len-1], 0.0);
  for(i=len-1; i>0; i--) {
    /* The following three lines are equivalent to
       ans = gsl_complex_add_real (gsl_complex_mul (z, ans), c[i-1]);
       but faster */
    double tmp = c[i-1] + GSL_REAL (z) * GSL_REAL (ans) - GSL_IMAG (z) * GSL_IMAG (ans);
    GSL_SET_IMAG (&ans, GSL_IMAG (z) * GSL_REAL (ans) + GSL_REAL (z) * GSL_IMAG (ans));
    GSL_SET_REAL (&ans, tmp);
  }
  return ans;
}

INLINE_FUN
gsl_complex
gsl_complex_poly_complex_eval(const gsl_complex c[], const int len, const gsl_complex z)
{
  int i;
  gsl_complex ans = c[len-1];
  for(i=len-1; i>0; i--) {
    /* The following three lines are equivalent to
       ans = gsl_complex_add (c[i-1], gsl_complex_mul (x, ans));
       but faster */
    double tmp = GSL_REAL (c[i-1]) + GSL_REAL (z) * GSL_REAL (ans) - GSL_IMAG (z) * GSL_IMAG (ans);
    GSL_SET_IMAG (&ans, GSL_IMAG (c[i-1]) + GSL_IMAG (z) * GSL_REAL (ans) + GSL_REAL (z) * GSL_IMAG (ans));
    GSL_SET_REAL (&ans, tmp);
  }
  return ans;
}
#endif /* HAVE_INLINE */

/* Work with divided-difference polynomials, Abramowitz & Stegun 25.2.26 */

int
gsl_poly_dd_init (double dd[], const double x[], const double y[],
                  size_t size);

INLINE_DECL double
gsl_poly_dd_eval (const double dd[], const double xa[], const size_t size, const double x);

#ifdef HAVE_INLINE
INLINE_FUN
double
gsl_poly_dd_eval(const double dd[], const double xa[], const size_t size, const double x)
{
  size_t i;
  double y = dd[size - 1];
  for (i = size - 1; i--;) y = dd[i] + (x - xa[i]) * y;
  return y;
}
#endif /* HAVE_INLINE */

int
gsl_poly_dd_taylor (double c[], double xp,
                    const double dd[], const double x[], size_t size,
                    double w[]);

/* Solve for real or complex roots of the standard quadratic equation,
 * returning the number of real roots.
 *
 * Roots are returned ordered.
 */
int gsl_poly_solve_quadratic (double a, double b, double c,
                              double * x0, double * x1);

int
gsl_poly_complex_solve_quadratic (double a, double b, double c,
                                  gsl_complex * z0, gsl_complex * z1);

/* Solve for real roots of the cubic equation
 * x^3 + a x^2 + b x + c = 0, returning the
 * number of real roots.
 *
 * Roots are returned ordered.
 */
int gsl_poly_solve_cubic (double a, double b, double c,
                          double * x0, double * x1, double * x2);

int
gsl_poly_complex_solve_cubic (double a, double b, double c,
                              gsl_complex * z0, gsl_complex * z1,
                              gsl_complex * z2);

/* Solve for the complex roots of a general real polynomial */

typedef struct
{
  size_t nc ;
  double * matrix ;
}
gsl_poly_complex_workspace ;

gsl_poly_complex_workspace * gsl_poly_complex_workspace_alloc (size_t n);
void gsl_poly_complex_workspace_free (gsl_poly_complex_workspace * w);

int
gsl_poly_complex_solve (const double * a, size_t n,
                        gsl_poly_complex_workspace * w,
                        gsl_complex_packed_ptr z);

__END_DECLS

#endif /* __GSL_POLY_H__ */

max_stars_repo_name: ielomariala/Hex-Game

max_issues_repo_name: ielomariala/Hex-Game

max_forks_repo_name: ielomariala/Hex-Game


text:
//setzt constanten auf
#define EXTERN_VARIABLES
#include <stdio.h>
#include <math.h>
#include <gsl/gsl_sf_bessel.h>
#include "global.h"
#include "utility.h"
#include "steppmethods.h"
#include <gsl/gsl_cblas.h>
#include <time.h>
#include <sys/types.h>
#include <sys/stat.h>
#include <unistd.h>
#include <omp.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>

void Set_Target_Length(double * TARGET_LENGTH, double volfrac,  void (*Stepper_Method)(int N,  double  x[N/2] , double v[N/2], double a[N/2], double *t,
					void (*deriv) (double *y, double *ans, double t,int N)),
					void (*Poti_Handle) (double *y, double *ans, double t,int N))
{	int FLAG = 1;
	*TARGET_LENGTH = 0.0;
	if (Stepper_Method == VVerlet_Step_deriv)
	{
		 *TARGET_LENGTH = 0.0;
		printf(" No  Targets\n");
		FLAG = 0;
	}
	//if (Stepper_Method == VVerlet_Step_Target_Circle)
	//{
	//	*TARGET_LENGTH = LATTICE_SPACING * sqrt(volfrac/ M_PI);
	//	printf(" Targets set to hard circles\n");
	//}

	if (Stepper_Method == VVerlet_Step_hardsphere_reflect)
	{
		 *TARGET_LENGTH = LATTICE_SPACING * sqrt(volfrac/ M_PI);
		printf(" Targets set to hard circles\n");
		FLAG = 0;
	}

	if (Stepper_Method == VVerlet_Step_hardsphere_reflect_Kupf)
	{
		 *TARGET_LENGTH = LATTICE_SPACING * sqrt(volfrac/ M_PI);
		printf(" Targets set to hard circles Kupferman \n");
		FLAG = 0;
	}

	if (Stepper_Method == VVerlet_Step_hardsphere_reflect_all)
	{
		 *TARGET_LENGTH = LATTICE_SPACING * sqrt(volfrac/ M_PI);
		printf(" Targets set to hard circles, reflect ALL particles\n");
		FLAG = 0;
	}

	if (Stepper_Method == VVerlet_Step_deriv_Kupf)
	{
		*TARGET_LENGTH = 0.0;
		printf(" No  Targets\n, Kupfermann bath");
		FLAG = 0;
	}

	if (Stepper_Method == VVerlet_Step_deriv_Box)
	{
		 //*TARGET_LENGTH = LATTICE_SPACING * sqrt(volfrac/ M_PI/5.0);
		//printf(" Targets set to soft Yukawa circles in 3 per pore corner\n");
		*TARGET_LENGTH = 0;
		FLAG = 0;
		if (DIM > 1)
		{
			printf("wrong Dimension, please switch to 1d");
			exit(1);
		}
		if (volfrac > 0.0)
		{
			printf("no need for volume fractions >0 \n");
			exit(1);
		}
		printf("box set up with hard wall \n");
	}
	if(FLAG)
	{
		printf("Kein bekannter Stepper in global, cannot set tarhet length!\n");
		exit(1);
	}
}

void Constants(){
// call to setup the bathh parameters and give rudimentary output for the run
	int i;
	OUTPUT_FLAG = 1;
	printf ("OSSZI = %d  DIM =%d  ORDER =%d  \n", OSSZI, DIM, ORDER);
	printf ("KBOLTZ = %e \nTIME_STEPS =%f  \nTIME_END =%f  \n", KBOLTZ, TIME_STEPS, TIME_END);
	for(i=0;i<OSSZI;i++){
		coupling[i] = 0.0;
	}
	for(i=0;i<OSSZI;i++){
		ommega[i] = 0.0;
	}
	for(i=0;i<OSSZI;i++){
		massq[i] = 1.0;
	}
	for(i=0;i<ORDER;i++){
		y[i] = 1.0;
	}
	printf("ALPHA = %1.2E  GAMMA = %1.2E  mass  = %1.2E \n", ALPHA, GAMMA, mass);
	printf("NUMBER_TO_INF = %d  \nTIME_ARRIVAL =%1.2E  \n ", NUMBER_TO_INF, TIME_ARRIVAL);
	printf("Startbedingungen = ");
	printf("%s",LABEL);
	printf("\n");
	if (!(OMMEGA_READ_BINARY))
	{
		// random Number Preparation for Taus GSL genrator
		gsl_rng * r;
		const gsl_rng_type * T;
		gsl_rng_env_setup();
		T = gsl_rng_taus2;
		r = gsl_rng_alloc (T);
		gsl_rng_set(r, time(NULL)); // Seed with time
		// Random number end //
	    Ommega_Setup(r);
	    gsl_rng_free (r);
	}else	// read in from ommega.dat
	{
		FILE *fp;
		fp = fopen("ommega.dat", "r");
		for (i=0;i<OSSZI;i++)
		{
			if (fscanf(fp, "%lf", &ommega[i]) != 1)
			{
	            printf("ERROT READING ommega.dat");
	            break;
        	}
		}
		fclose(fp);
		printf("read in ommega.dat\n");
	}
	if(VIRTUAL_FLAG) printf("Virtuelles Teilchen benutzt, ohne Badrückkopplung");
  	Gamma_Setup(coupling);

  	double ommega_max = ommega[cblas_idamax(OSSZI, ommega, 1)];
  	double coupling_max = coupling[cblas_idamax(OSSZI, coupling, 1)];
  	printf("MAX coupling =  %1.2e\n",coupling_max);
  	if(!(MAXSTEPSIZE_BINARY))
	{
  		MAXSTEPSIZE = (2*M_PI/ommega[cblas_idamax(OSSZI, ommega, 1)] * MAXSTEPSIZEMULT);
  		printf("MAXSTEPSIZE =  %1.2e set according to smallest period times %3.3e\n",MAXSTEPSIZE, MAXSTEPSIZEMULT);
  	}else
  	{
  		MAXSTEPSIZE = (MAXSTEPSIZEMULT / (double) OSSZI);
  		printf("MAXSTEPSIZE =  %1.2e set according to OSSZInr times Stepsize dt = %3.3e\n",MAXSTEPSIZE, MAXSTEPSIZEMULT);
  	}
  	Sum_Method =  SUMMATION_METHOD;
  	Stepper_Method = STEPPER_METHOD;
  	Poti_Handle = POTI_HANDLE;

  	int m = ipow(NUMBER_TO_INF,DIM); int n = DIM;
  	// Setze Gitter auf auf globale Matrix POSITIONS
	POSITIONS = (double **) malloc(m * sizeof(int *));
	POSITIONS[0] = (double *) malloc(m* n * sizeof(int));
	for (i=1; i<m; i++) POSITIONS[i] = POSITIONS[0] + n *i;
	vec_zero_i(lattice_position,DIM);
	printf("lattice allocated, current lattice cell set to zero\n");
	TARGET_CONTACT = 0;
	OUTPUT_FLAG = 0; // subsequent setup calls should not display text
	printf("Running on %d threads, optimally DIM = %d\n", THREADS, DIM);

	// set masses for chosen type of coupling
	if( !(KUPF_BINARY) )
	{
		for (int osi = 0; osi < OSSZI; osi ++)
		{
			massq[osi] 		= 1.0;
		}
		printf("Set all masses équal\n");
	}else
	{
		for (int osi = 0; osi < OSSZI; osi ++)
		{
			massq[osi] 		= coupling[osi] / ommega[osi] / ommega[osi];

		}
		printf("Set all masses according to k = m w^2l\n");
	}
	// Setup masses of bath for simpsons rule, split ends, odd and even parts

	if (SIMPSON_BINARY)
	{
		for (int osi = 1; osi < OSSZI/2-1; osi ++)
		{
			massq[osi * 2] 		*= 2.0/3.0;
			massq[osi * 2 + 1] 	*= 4.0/3.0;
		}
		massq[0] 		*= 1.0/3.0;;
		massq[OSSZI-1] 	*= 1.0/3.0;

		if (KUPF_BINARY)
		{
			for (int osi = 0; osi < OSSZI; osi ++)
			{
				coupling[osi] = massq[osi] * ommega[osi] * ommega[osi];
			}
		}
	}
}

ext: c


text:
#include <math.h>
#include <stdlib.h>
#include <stdio.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_integration.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_odeiv.h>
#include "growthfactor.h"

//Code to calculate the linear growth factor D, and linear growth rate, f. Where a linear perturbation delta grows as
//delta(a') = delta(a)*(D(a')/D(a))^2 and f = dlnD/dlna
//Note anyone using Komatsu's CRL library should note: growth_factor_crl = D *(1+z) and growth_rate_crl = f/(1+z)

double D;
double linear_f;
double w0,wa;
double omega_m,omega_lambda;

#define LIMIT_SIZE 1000

int get_growthfactor(double a,double om, double w, double w2, double *gf)
{
	w0 = w;
	wa = w2;
	omega_m = om;
	omega_lambda = 1.0 - omega_m;
        growth_de(a,gf);
	return 0;
}

// 2 parameter dark energy equation of state

double w (double a)
{
return w0 + (1.0-a)*wa;
}

double w_int(double z, void *param)
{
	return (1. + w(1./(z+1)) )/( 1. + z);
}

double Xde_int (double a,void *params )
{
    if (a == 0.) a = 1e-3;
    double Xde_int = w(a)/a;
    return Xde_int;
}

double Xde (double a)
{
    gsl_integration_workspace * worksp  = gsl_integration_workspace_alloc (LIMIT_SIZE);
    gsl_function F;
    F.function = &Xde_int;
    F.params =0;
    double integral,error;
    if (a == 0.)
		a = 1e-3;
    gsl_integration_qags (&F, a, 1,  1.0e-20, 1.0e-10, LIMIT_SIZE,worksp, &integral, &error);
    gsl_integration_workspace_free (worksp);

    return omega_m/(1.-omega_m)*exp(-3.*integral);
}

int func (double t, const double y[], double f[], void *params)
{
	double mu = *(double *)params;
	f[0] = y[1];
	f[1] = -( 3.5 - 1.5 * w(t)/( 1 + Xde(t) ) )*y[1]/t - 1.5 * ( 1 - w(t) )/( 1 + Xde(t))*(y[0]/t/t);
	return GSL_SUCCESS;
}

int jac (double t, const double y[], double *dfdy, double dfdt[], void *params)
{
	double mu = *(double *)params;
	gsl_matrix_view dfdy_mat = gsl_matrix_view_array (dfdy, 2, 2);
	gsl_matrix * m = &dfdy_mat.matrix;
	gsl_matrix_set (m, 0, 0, 0.0); //dy1/dx1
	gsl_matrix_set (m, 0, 1, 1.0);
	gsl_matrix_set (m, 1, 0, - 1.5 * ( 1 - w(t) )/( 1 + Xde(t))*(1./t/t));
	gsl_matrix_set (m, 1, 1, -( 3.5 - 1.5 * w(t)/( 1 + Xde(t) ) )/t);
	dfdt[0] = 0.0;
	dfdt[1] = 0.0;
	return GSL_SUCCESS;
}

double growth_de (double a, double *gf)
{
	const gsl_odeiv_step_type * T
		= gsl_odeiv_step_rk4;

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

	double mu = 10;
	gsl_odeiv_system sys = {func, jac, 2, &mu};

	double t =1.e-3, t1 = a;
	double h = 1e-6;
	double y[2] = {1.,0.};

	while (t < t1)
	{
		int status = gsl_odeiv_evolve_apply (e, c, s, &sys, &t, t1, &h, y);
		if (status != GSL_SUCCESS)
			break;
	}
	gsl_odeiv_evolve_free (e);
	gsl_odeiv_control_free (c);
	gsl_odeiv_step_free (s);

//	return y[0]*a;  //D growth factor
//	return y[1]*a*a/(y[0]*a) +1.; // f = d ln D/ d ln a
	gf[0] = y[0]*a;
	gf[1] = y[1]*a*a/(y[0]*a) +1.;
	return y[0]*a;

}


text:
#pragma once
//#include <gsl/string_span>
namespace Util {
	//static gsl::string_span<> toLower(gsl::string_span<> const&);
}; // namespace util

ext: h

max_issues_repo_name: TimScriptov/MCPELauncher

max_forks_repo_name: TimScriptov/MCPELauncher


text:
#ifndef LIC_METHOD_DEFINITION
#define LIC_METHOD_DEFINITION

#include <gsl/gsl>

#include "Metadata.h"

namespace lic
{

class MethodDefinition : Metadata
{
public:
    MethodDefinition(Assembly& assembly, MetadataTable table, size_t rid, gsl::byte* row);
    ~MethodDefinition();

    const char* Name() const;
    const gsl::byte* Code() const;
};

}

#endif // !LIC_METHOD_DEFINITION

max_stars_repo_name: roberthusak/lic

max_issues_repo_name: roberthusak/lic

max_forks_repo_name: roberthusak/lic


text:
#ifndef __GSL_MATRIX_H__
#define __GSL_MATRIX_H__

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

#include <gsl/gsl_matrix_complex_long_double.h>
#include <gsl/gsl_matrix_complex_double.h>
#include <gsl/gsl_matrix_complex_float.h>

#include <gsl/gsl_matrix_long_double.h>
#include <gsl/gsl_matrix_double.h>
#include <gsl/gsl_matrix_float.h>

#include <gsl/gsl_matrix_ulong.h>
#include <gsl/gsl_matrix_long.h>

#include <gsl/gsl_matrix_uint.h>
#include <gsl/gsl_matrix_int.h>

#include <gsl/gsl_matrix_ushort.h>
#include <gsl/gsl_matrix_short.h>

#include <gsl/gsl_matrix_uchar.h>
#include <gsl/gsl_matrix_char.h>

#endif /* __GSL_MATRIX_H__ */

max_stars_repo_name: berkus/music-cs

max_issues_repo_name: berkus/music-cs

max_forks_repo_name: berkus/music-cs


text:
/* examples/C/ssmfe/shift_invert.c */
/* Laplacian on a rectangular grid by shift-invert via LDLT factorization */
#include "spral.h"
#include <math.h>
#include <stdlib.h>
#include <stdio.h>
#include <cblas.h>

/* Headers that implements Laplacian and preconditioners and LDLT support */
#include "laplace2d.h"
#include "ldltf.h"

int main(void) {
   const int nx = 8;          /* grid points along x */
   const int ny = 8;          /* grid points along y */
   const int n = nx*ny;       /* problem size */
   const double sigma = 1.0;  /* shift */

   int ipiv[n];               /* LDLT pivot index */
   double lambda[n];          /* eigenvalues */
   double X[n][n];            /* eigenvectors */
   double A[n][n];            /* matrix */
   double LDLT[n][n];         /* factors */
   double work[n][n];         /* work array for dsytrf */
   struct spral_ssmfe_options options;    /* eigensolver options */
   struct spral_ssmfe_inform inform;      /* information */
   struct spral_ssmfe_rcid rci;           /* reverse communication data */
   void *keep;                            /* private data */

   /* Initialize options to default values */
   spral_ssmfe_default_options(&options);

   /* Set up then perform LDLT factorization of the shifted matrix */
   set_laplacian_matrix(nx, ny, n, A);
   for(int j=0; j<n; j++)
   for(int i=0; i<n; i++)
      LDLT[j][i] = (i==j) ? A[j][i] - sigma
                          : A[j][i];
   cwrap_dsytrf('L', n, &LDLT[0][0], n, ipiv, &work[0][0], n*n);

   /* Main loop */
   int left = num_neg_D(n, n, LDLT, ipiv);   /* all evalues to left of sigma */
   int right = 5;                            /* 5 evalues to right of sigma */
   rci.job = 0; keep = NULL;
   while(true) {
      spral_ssmfe_standard_shift_double(&rci, sigma, left, right, n, lambda,
            n, &X[0][0], n, &keep, &options, &inform);
      switch( rci.job ) {
      case 1:
         cblas_dgemm(CblasColMajor, CblasNoTrans, CblasNoTrans, n, rci.nx, n,
            1.0, &A[0][0], n, rci.x, n, 0.0, rci.y, n);
         break;
      case 2:
         // No preconditioning
         break;
      case 9:
         cblas_dcopy(n*rci.nx, rci.x, 1, rci.y, 1);
         cwrap_dsytrs('L', n, rci.nx, &LDLT[0][0], n, ipiv, rci.y, n);
         break;
      default:
         goto finished;
      }
   }
finished:
   printf("Eigenvalues near %e (took %d iterations)\n", sigma, inform.iteration);
   for(int i=0; i<inform.left+inform.right; i++)
      printf(" lambda[%1d] = %13.7e\n", i, lambda[i]);
   spral_ssmfe_free_double(&keep, &inform);

   /* Success */
   return 0;
}

max_stars_repo_name: mjacobse/spral

max_issues_repo_name: mjacobse/spral

max_forks_repo_name: mjacobse/spral


text:
///////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
//
//
// 			     //////////////////////////////////////////////////////////
//  			     //							     //
// 			     //   	        hybridMANTIS v1.0		     //
// 			     //   	       fastDETECT2 - C code  		     //
//			     //		   (optical photons transport)		     //
//			     //							     //
//			     //////////////////////////////////////////////////////////
//
//
//
//
// ****Disclaimer****
//  This software and documentation (the "Software") were developed at the Food and Drug Administration (FDA) by employees of the Federal Government in
//  the course of their official duties. Pursuant to Title 17, Section 105 of the United States Code, this work is not subject to copyright protection
//  and is in the public domain. Permission is hereby granted, free of charge, to any person obtaining a copy of the Software, to deal in the Software
//  without restriction, including without limitation the rights to use, copy, modify, merge, publish, distribute, sublicense, or sell copies of the
//  Software or derivatives, and to permit persons to whom the Software is furnished to do so. FDA assumes no responsibility whatsoever for use by other
//  parties of the Software, its source code, documentation or compiled executables, and makes no guarantees, expressed or implied, about its quality,
//  reliability, or any other characteristic. Further, use of this code in no way implies endorsement by the FDA or confers any advantage in regulatory
//  decisions. Although this software can be redistributed and/or modified freely, we ask that any derivative works bear some notice that they are
//  derived from it, and any modified versions bear some notice that they have been modified.
//
//
//	Associated publication: Sharma Diksha, Badal Andreu and Badano Aldo, "hybridMANTIS: a CPU-GPU Monte Carlo method for modeling indirect x-ray detectors with
//				columnar scintillators". Physics in Medicine and Biology, 57(8), pp. 2357–2372 (2012)
//
//
//	File:   	hybridMANTIS_c_ver1_0.c
//	Author: 	Diksha Sharma (US Food and Drug Administration)
//	Email: 		diksha.sharma@fda.hhs.gov
//	Last updated:  	Apr 13, 2012
//
///////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////

/////////////////////////////////////////
//
//      Header libraries
//
/////////////////////////////////////////

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

/////////////////////////////////////////
//
//      Global variables
//
/////////////////////////////////////////

	#define max_photon_per_EDE 900000	// maximum number of optical photons that can be generated per energy deposition event (EDE)

	#ifndef USING_CUDA
		#define mybufsize 2304000	// CPU buffer size: # of events sent to the CPU
	#endif

/////////////////////////////////////////
//
//      Include kernel program
//
/////////////////////////////////////////
	#include "kernel_cuda_c_ver1_0.cu"

////////////////////////////////////////////////////////////////////////////
//				MAIN PROGRAM			          //
////////////////////////////////////////////////////////////////////////////

#ifndef USING_CUDA

///////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
//
// cpuoptical(): Performs optical transport using CPU
//	  	 Input arguments: gpusize
//
// 		 gpusize - buffer size already processed in the GPU. CPU runs optical transport for rest of the buffer.
//
///////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////

	void cpuoptical_(int *gpusize)
	{

		// command line arguments
		float xdetector, ydetector, radius, height, n_C, n_IC, top_absfrac, bulk_abscoeff, beta, d_min, lbound_x, lbound_y, ubound_x, ubound_y, d_max, yield, sensorRefl;
		int pixelsize, num_primary, min_optphotons, max_optphotons, num_bins;

		float dcos[3]={0}; 		// directional cosines
		float normal[3]={0}; 		// normal to surface in case of TIR
		float pos[3] = {0}; 		// intersecting point coordinates
		float old_pos[3] = {0};  	// source coordinates

		int nbytes = (optical_.myctropt - (*gpusize))*sizeof(struct start_info);
		struct start_info *structa;
		structa = (struct start_info*) malloc(nbytes);
		if( structa == NULL )
			printf("\n Struct start_info array CANNOT BE ALLOCATED - %d !!", (optical_.myctropt-(*gpusize)));

		// get cpu time
		clock_t start, end;
		float num_sec;

		// get current time stamp to initialize seed input for RNG
		time_t seconds;
		seconds = time (NULL);
		struct timeval tv;
		gettimeofday(&tv,NULL);

		float rr=0.0f, theta=0.0f;
		float r=0.0f;			// random number
		float norm=0.0f;
		int jj=0;
		int my_index=0;
		int penindex = 0;		// equal to *gpusize
		int result_algo = 0;
		unsigned long long int *num_rebound;

		// initialize random number generator (RANECU)
		int seed_input = 271828182 ; 		// ranecu seed input
		int seed[2];

		// GNU scientific library (gsl) variables
		const gsl_rng_type * Tgsl;
		gsl_rng * rgsl;
		double mu_gsl;

		// output image variables
		int xdim = 0;
		int ydim = 0;
		int indexi=0, indexj=0;

		// copy to local variables from PENELOPE buffers
		xdetector = inputargs_.detx;		// x dimension of detector (in um). x in (0,xdetector)
		ydetector = inputargs_.dety;		// y dimension of detector (in um). y in (0,ydetector)
		height = inputargs_.detheight;		// height of column and thickness of detector (in um). z in range (-H/2, H/2)
		radius = inputargs_.detradius;		// radius of column (in um).
		n_C = inputargs_.detnC;			// refractive index of columns
		n_IC = inputargs_.detnIC;		// refractive index of intercolumnar material
		top_absfrac = inputargs_.dettop;	// column's top surface absorption fraction (0.0, 0.5, 0.98)
		bulk_abscoeff = inputargs_.detbulk;	// column's bulk absorption coefficient (in um^-1) (0.001, 0.1 cm^-1)
		beta = inputargs_.detbeta;		// roughness coefficient of column walls
		d_min = inputargs_.detdmin;		// minimum distance a photon can travel when transmitted from a column
		d_max = inputargs_.detdmax;
		lbound_x = inputargs_.detlboundx;	// x lower bound of region of interest of output image (in um)
		lbound_y = inputargs_.detlboundy;	// y lower bound (in um)
		ubound_x = inputargs_.detuboundx;	// x upper bound (in um)
		ubound_y = inputargs_.detuboundy;	// y upper bound (in um)
		yield = inputargs_.detyield;		// yield (/eV)
		pixelsize = inputargs_.detpixel;	// 1 pixel = pixelsize microns (in um)
		sensorRefl = inputargs_.detsensorRefl;	// Non-Ideal sensor reflectivity (%)
		num_primary = inputargs_.mynumhist;	// total number of primaries to be simulated
		min_optphotons = inputargs_.minphotons;	// minimum number of optical photons detected to be included in PHS
		max_optphotons = inputargs_.maxphotons;	// maximum number of optical photons detected to be included in PHS
		num_bins = inputargs_.mynumbins;	// number of bins for genrating PHS

	      	// create a generator chosen by the environment variable GSL_RNG_TYPE
	       	gsl_rng_env_setup();
	       	Tgsl = gsl_rng_default;
	       	rgsl = gsl_rng_alloc (Tgsl);

	       	// dimensions of PRF image
		xdim = ceil((ubound_x - lbound_x)/pixelsize);
		ydim = ceil((ubound_y - lbound_y)/pixelsize);
		unsigned long long int myimage[xdim][ydim];

		//initialize the output image 2D array
		for(indexi = 0; indexi < xdim; indexi++)
		{
		 for(indexj = 0; indexj < ydim; indexj++)
		 {
		  myimage[indexi][indexj] = 0;
		 }
		}

		// memory for storing histogram of # photons detected/primary
		int *h_num_detected_prim = 0;
		h_num_detected_prim = (int*)malloc(sizeof(int)*num_primary);

		for(indexj=0; indexj < num_primary; indexj++)
		  h_num_detected_prim[indexj] = 0;

		// memory for storing histogram of # photons detected/primary
		int *h_histogram = 0;
		h_histogram = (int*)malloc(sizeof(int)*num_bins);

		for(indexj=0; indexj < num_bins; indexj++)
		  h_histogram[indexj] = 0;

	penindex = *gpusize;

	for(my_index = 0; my_index < (optical_.myctropt-(*gpusize)); my_index++)		// iterate over x-rays
	{

		// reset the global counters
		num_generated = 0;
		num_detect=0;
		num_abs_top=0;
		num_abs_bulk=0;
		num_lost=0;
		num_outofcol=0;
		num_theta1=0;
		photon_distance=0.0f;

		//re-initialize the output image 2D array
		for(indexi = 0; indexi < xdim; indexi++)
		 for(indexj = 0; indexj < ydim; indexj++)
		   myimage[indexi][indexj] = 0;

		// copying PENELOPE buffer into *structa

		// units in the penelope output file are in cm. Convert them to microns.
		structa[my_index].str_x = optical_.xbufopt[penindex+my_index] * 10000.0f;	// x-coordinate of interaction event
		structa[my_index].str_y = optical_.ybufopt[penindex+my_index] * 10000.0f;	// y-coordinate
		structa[my_index].str_z = optical_.zbufopt[penindex+my_index] * 10000.0f;	// z-coordinate
		structa[my_index].str_E = optical_.debufopt[penindex+my_index];			// energy deposited
		structa[my_index].str_histnum = optical_.nbufopt[penindex+my_index];		// x-ray history

		// sample # optical photons based on light yield and energy deposited for this interaction event (using Poisson distribution)
		mu_gsl = (double)structa[my_index].str_E * yield;
		structa[my_index].str_N = gsl_ran_poisson(rgsl,mu_gsl);

		if(structa[my_index].str_N > max_photon_per_EDE)
		{
			printf("\n\n str_n exceeds max photons. program is exiting - %d !! \n\n",structa[my_index].str_N);
			exit(0);
		}

		num_rebound = (unsigned long long int*) malloc(structa[my_index].str_N*sizeof(unsigned long long int));
		if(num_rebound == NULL)
		printf("\n Error allocating num_rebound memory !\n");

		// start the clock
		start = clock();

		// initialize the RANECU generator in a position far away from the previous history:
		seed_input = (int)(seconds/3600+tv.tv_usec);			// seed input=seconds passed since 1970+current time in micro secs
		init_PRNG(my_index, 50000, seed_input, seed);      		// intialize RNG

		for(jj=0; jj<structa[my_index].str_N; jj++)
			num_rebound[jj] = 0;

		// reset the directional cosine and normal vectors
		dcos[0]=0.0f; dcos[1]=0.0f; dcos[2]=0.0f;
		normal[0]=0.0f; normal[1]=0.0f; normal[2]=0.0f;

		// re-initialize myimage
		for(indexi = 0; indexi < xdim; indexi++)
		 for(indexj = 0; indexj < ydim; indexj++)
			  myimage[indexi][indexj] = 0;

		// set starting location of photon
		pos[0] = structa[my_index].str_x; pos[1] = structa[my_index].str_y; pos[2] = structa[my_index].str_z;
		old_pos[0] = structa[my_index].str_x; old_pos[1] = structa[my_index].str_y; old_pos[2] = structa[my_index].str_z;

		// initializing the direction cosines for the first particle in each core
		r = (ranecu(seed) * 2.0f) - 1.0f; 	// random number between (-1,1)

		while(fabs(r) <= 0.01f)
		 {
		   	r = (ranecu(seed) * 2.0f) - 1.0f;
		 }

		dcos[2] = r;		// random number between (-1,1)
		rr = sqrt(1.0f-r*r);
		theta=ranecu(seed)*twopipen;
		dcos[0]=rr*cos(theta);
		dcos[1]=rr*sin(theta);

		norm = sqrt(dcos[0]*dcos[0] + dcos[1]*dcos[1] + dcos[2]*dcos[2]);

		if ((norm < (1.0f - epsilon)) || (norm > (1.0f + epsilon)))	// normalize
		 {
			dcos[0] = dcos[0]/norm;
			dcos[1] = dcos[1]/norm;
			dcos[2] = dcos[2]/norm;
		 }

		local_counter=0;	// total number of photons terminated (either detected at sensor, absorbed at the top or in the bulk) [global variable]
		while(local_counter < structa[my_index].str_N)		// until all the optical photons are not transported
		 {

			absorbed = 0;
			detect = 0;
			bulk_abs = 0;

			// set starting location of photon
			pos[0] = structa[my_index].str_x; pos[1] = structa[my_index].str_y; pos[2] = structa[my_index].str_z;
			old_pos[0] = structa[my_index].str_x; old_pos[1] = structa[my_index].str_y; old_pos[2] = structa[my_index].str_z;
			num_generated++;
			result_algo = 0;

			while(result_algo == 0)
			 {
			  	result_algo = algo(normal, old_pos, pos, dcos, num_rebound, seed, structa[my_index], &myimage[0][0], xdetector, ydetector, radius, height, n_C, n_IC, top_absfrac, bulk_abscoeff, beta, d_min, pixelsize, lbound_x, lbound_y, ubound_x, ubound_y, sensorRefl, d_max, ydim, h_num_detected_prim);
			 }

		 }

		// end the clock
		end = clock();

		num_sec = (((float)end - start)/CLOCKS_PER_SEC);

		 // type cast unsigned long long int to double
		 double cast_num_generated;
		 double cast_num_detect;
		 double cast_num_abs_top;
		 double cast_num_abs_bulk;
		 double cast_num_lost;
		 double cast_num_outofcol;
		 double cast_num_theta1;
		 double cast_gputime;

		 cast_num_generated = (double)num_generated;
		 cast_num_detect    = (double)num_detect;
		 cast_num_abs_top   = (double)num_abs_top;
		 cast_num_abs_bulk  = (double)num_abs_bulk;
		 cast_num_lost      = (double)num_lost;
		 cast_num_outofcol  = (double)num_outofcol;
		 cast_num_theta1    = (double)num_theta1;
		 cast_gputime	    = (double)(num_sec*1000.0);		// convert in millisecond

		 // save to global counters
		 optstats_.glgen      = optstats_.glgen      + cast_num_generated;
		 optstats_.gldetect   = optstats_.gldetect   + cast_num_detect;
		 optstats_.glabstop   = optstats_.glabstop   + cast_num_abs_top;
		 optstats_.glabsbulk  = optstats_.glabsbulk  + cast_num_abs_bulk;
		 optstats_.gllost     = optstats_.gllost     + cast_num_lost;
		 optstats_.gloutofcol = optstats_.gloutofcol + cast_num_outofcol;
		 optstats_.gltheta1   = optstats_.gltheta1   + cast_num_theta1;
		 optstats_.glgputime  = optstats_.glgputime  + cast_gputime;

		// release resources
		free(num_rebound);

	}	// my_index loop ends

	// make histogram of number of detected photons/primary for num_bins
	int binsize=0, newbin=0;
	int bincorr=0;

	binsize = floor((max_optphotons-min_optphotons)/num_bins);	// calculate size of each bin. Assuming equally spaced bins.
	bincorr = floor(min_optphotons/binsize);			// correction in bin number if min_optphotons > 0.

	for(indexi = 0; indexi < num_primary; indexi++)
	 {
	 	newbin = floor(h_num_detected_prim[indexi]/binsize) - bincorr;	// find bin #

	 	if(h_num_detected_prim[indexi] > 0)	// store only non-zero bins
	 	{
		 	if(h_num_detected_prim[indexi] <= min_optphotons)	// # detected < minimum photons given by user, add to the first bin
		 		h_histogram[0]++;
		 	else if(h_num_detected_prim[indexi] >= max_optphotons)	// # detected > maximum photons given by user, then add to the last bin
		 		h_histogram[num_bins-1]++;
		 	else
		 		h_histogram[newbin]++;
		}
	 }

	// add num_detected_primary to gldetprimary array in PENELOPE
	for(indexi = 0; indexi < num_bins; indexi++)
		outputdetprim_.gldetprimary[indexi] = outputdetprim_.gldetprimary[indexi] + h_histogram[indexi];

	// release resources
	free(structa);
	free(h_num_detected_prim);
	free(h_histogram);

		return;
	}	// C main() ends

#endif

max_stars_repo_name: diamfda/hybridmantis

max_issues_repo_name: diamfda/hybridmantis

max_forks_repo_name: diamfda/hybridmantis


text:
#pragma once

#include <imageview/ContinuousImageView.h>
#include <imageview/internal/ImageViewStorage.h>
#include <imageview/internal/PixelRef.h>

#include <gsl/span>

#include <cstddef>

namespace imageview {

template <class PixelFormat, bool Mutable = false>
class ImageView {
 public:
  static_assert(IsPixelFormat<PixelFormat>::value, "Not a PixelFormat.");

  using byte_type = std::conditional_t<Mutable, std::byte, const std::byte>;
  using value_type = typename PixelFormat::color_type;
  // TODO: consider using 'const value_type' instead of 'value_type' for immutable views.
  using reference = std::conditional_t<Mutable, detail::PixelRef<PixelFormat>, value_type>;

  // Construct an empty view.
  template <class Enable = std::enable_if_t<std::is_default_constructible_v<PixelFormat>>>
  constexpr ImageView() noexcept(noexcept(std::is_nothrow_default_constructible_v<PixelFormat>)) {}

  // Construct a view into an image.
  // \param height - height of the image.
  // \param width - width of the image.
  // \param stride -
  // \param data - bitmap data.
  constexpr ImageView(unsigned int height, unsigned int width, unsigned int stride, gsl::span<byte_type> data);

  // Construct a view into an image.
  // \param height - height of the image.
  // \param width - width of the image.
  // \param stride -
  // \param data - bitmap data.
  // \param pixel_format - PixelFormat instance to use.
  constexpr ImageView(unsigned int height, unsigned int width, unsigned int stride, gsl::span<byte_type> data,
                      const PixelFormat& pixel_format);
  constexpr ImageView(unsigned int height, unsigned int width, unsigned int stride, gsl::span<byte_type> data,
                      PixelFormat&& pixel_format);

  // Construct a read-only view from a mutable view.
  template <class Enable = std::enable_if_t<!Mutable>>
  constexpr ImageView(ImageView<PixelFormat, !Mutable> other);

  constexpr ImageView(ContinuousImageView<PixelFormat, Mutable> image);
  // Construct a read-only view from a mutable contiguous view.
  template <class Enable = std::enable_if_t<!Mutable>>
  constexpr ImageView(ContinuousImageView<PixelFormat, !Mutable> image);

  // Return the height of the image
  constexpr unsigned int height() const noexcept;

  constexpr unsigned int width() const noexcept;

  constexpr unsigned int stride() const noexcept;

  constexpr unsigned int area() const noexcept;

  // Returns true if the image has zero area, false otherwise.
  constexpr bool empty() const noexcept;

  // Returns the pixel format used by this image.
  constexpr const PixelFormat& pixelFormat() const noexcept;

  constexpr gsl::span<byte_type> data() const noexcept;

  constexpr reference operator()(unsigned int y, unsigned int x) const;

  constexpr ImageRowView<PixelFormat, Mutable> row(unsigned int y) const;

 private:
  detail::ImageViewStorage<PixelFormat, Mutable> storage_;
  unsigned int height_ = 0;
  unsigned int width_ = 0;
  unsigned int stride_ = 0;
};

template <class PixelFormat, bool Mutable>
constexpr ImageView<PixelFormat, Mutable>::ImageView(unsigned int height, unsigned int width, unsigned int stride,
                                                     gsl::span<byte_type> data)
    : storage_(data.data()), height_(height), width_(width), stride_(stride) {
  Expects(width <= stride);
  const std::size_t expected_size = (height == 0) ? 0 : ((height - 1) * stride + width) * PixelFormat::kBytesPerPixel;
  Expects(data.size() == expected_size);
}

template <class PixelFormat, bool Mutable>
constexpr ImageView<PixelFormat, Mutable>::ImageView(unsigned int height, unsigned int width, unsigned int stride,
                                                     gsl::span<byte_type> data, const PixelFormat& pixel_format)
    : storage_(data.data(), pixel_format), height_(height), width_(width), stride_(stride) {
  Expects(width <= stride);
  const std::size_t expected_size = (height == 0) ? 0 : ((height - 1) * stride + width) * PixelFormat::kBytesPerPixel;
  Expects(data.size() == expected_size);
}

template <class PixelFormat, bool Mutable>
constexpr ImageView<PixelFormat, Mutable>::ImageView(unsigned int height, unsigned int width, unsigned int stride,
                                                     gsl::span<byte_type> data, PixelFormat&& pixel_format)
    : storage_(data.data(), std::move(pixel_format)), height_(height), width_(width), stride_(stride) {
  Expects(width <= stride);
  const std::size_t expected_size = (height == 0) ? 0 : ((height - 1) * stride + width) * PixelFormat::kBytesPerPixel;
  Expects(data.size() == expected_size);
}

template <class PixelFormat, bool Mutable>
template <class Enable>
constexpr ImageView<PixelFormat, Mutable>::ImageView(ImageView<PixelFormat, !Mutable> other)
    : ImageView(other.height(), other.width(), other.stride(), other.data(), other.pixelFormat()) {}

template <class PixelFormat, bool Mutable>
constexpr ImageView<PixelFormat, Mutable>::ImageView(ContinuousImageView<PixelFormat, Mutable> image)
    : ImageView(image.height(), image.width(), image.width(), image.data()) {}

template <class PixelFormat, bool Mutable>
template <class Enable>
constexpr ImageView<PixelFormat, Mutable>::ImageView(ContinuousImageView<PixelFormat, !Mutable> image)
    : ImageView(ContinuousImageView<PixelFormat, Mutable>(image)) {}

template <class PixelFormat, bool Mutable>
constexpr unsigned int ImageView<PixelFormat, Mutable>::height() const noexcept {
  return height_;
}

template <class PixelFormat, bool Mutable>
constexpr unsigned int ImageView<PixelFormat, Mutable>::width() const noexcept {
  return width_;
}

template <class PixelFormat, bool Mutable>
constexpr unsigned int ImageView<PixelFormat, Mutable>::stride() const noexcept {
  return stride_;
}

template <class PixelFormat, bool Mutable>
constexpr unsigned int ImageView<PixelFormat, Mutable>::area() const noexcept {
  return height_ * width_;
}

template <class PixelFormat, bool Mutable>
constexpr bool ImageView<PixelFormat, Mutable>::empty() const noexcept {
  return height_ == 0 || width_ == 0;
}

template <class PixelFormat, bool Mutable>
constexpr const PixelFormat& ImageView<PixelFormat, Mutable>::pixelFormat() const noexcept {
  return storage_.pixelFormat();
}

template <class PixelFormat, bool Mutable>
constexpr auto ImageView<PixelFormat, Mutable>::data() const noexcept -> gsl::span<byte_type> {
  const std::size_t data_size = (height_ == 0) ? 0 : ((height_ - 1) * stride_ + width_) * PixelFormat::kBytesPerPixel;
  return gsl::span<byte_type>(storage_.data(), data_size);
}

template <class PixelFormat, bool Mutable>
constexpr auto ImageView<PixelFormat, Mutable>::operator()(unsigned int y, unsigned int x) const -> reference {
  Expects(y < height_);
  Expects(x < width_);
  const gsl::span<byte_type, PixelFormat::kBytesPerPixel> pixel_data(
      storage_.data() + (y * stride_ + x) * PixelFormat::kBytesPerPixel, PixelFormat::kBytesPerPixel);
  if constexpr (Mutable) {
    return detail::PixelRef<PixelFormat>(pixel_data, pixelFormat());
  } else {
    return pixelFormat().read(pixel_data);
  }
}

template <class PixelFormat, bool Mutable>
constexpr ImageRowView<PixelFormat, Mutable> ImageView<PixelFormat, Mutable>::row(unsigned int y) const {
  Expects(y < height_);
  const gsl::span<byte_type> row_data(storage_.data() + y * stride_ * PixelFormat::kBytesPerPixel,
                                      width_ * PixelFormat::kBytesPerPixel);
  return ImageRowView<PixelFormat>(row_data, width_, pixelFormat());
}

}  // namespace imageview

max_stars_repo_name: alexanderbelous/imageview

max_issues_repo_name: alexanderbelous/imageview

max_forks_repo_name: alexanderbelous/imageview


text:
/**
 *
 * @file core_ztstrf.c
 *
 *  PLASMA core_blas kernel
 *  PLASMA is a software package provided by Univ. of Tennessee,
 *  Univ. of California Berkeley and Univ. of Colorado Denver
 *
 * @version 2.6.0
 * @author Hatem Ltaief
 * @author Mathieu Faverge
 * @author Jakub Kurzak
 * @date 2010-11-15
 * @precisions normal z -> c d s
 *
 **/
#include "common.h"
#include <cblas.h>
#include <math.h>

/***************************************************************************//**
 *
 * @ingroup CORE_PLASMA_Complex64_t
 *
 *  CORE_ztstrf computes an LU factorization of a complex matrix formed
 *  by an upper triangular NB-by-N tile U on top of a M-by-N tile A
 *  using partial pivoting with row interchanges.
 *
 *  This is the right-looking Level 2.5 BLAS version of the algorithm.
 *
 *******************************************************************************
 *
 * @param[in] M
 *         The number of rows of the tile A.  M >= 0.
 *
 * @param[in] N
 *         The number of columns of the tile A.  N >= 0.
 *
 * @param[in] IB
 *         The inner-blocking size.  IB >= 0.
 *
 * @param[in] NB
 *
 * @param[in,out] U
 *         On entry, the NB-by-N upper triangular tile.
 *         On exit, the new factor U from the factorization
 *
 * @param[in] LDU
 *         The leading dimension of the array U.  LDU >= max(1,NB).
 *
 * @param[in,out] A
 *         On entry, the M-by-N tile to be factored.
 *         On exit, the factor L from the factorization
 *
 * @param[in] LDA
 *         The leading dimension of the array A.  LDA >= max(1,M).
 *
 * @param[in,out] L
 *         On entry, the IB-by-N lower triangular tile.
 *         On exit, the interchanged rows form the tile A in case of pivoting.
 *
 * @param[in] LDL
 *         The leading dimension of the array L.  LDL >= max(1,IB).
 *
 * @param[out] IPIV
 *         The pivot indices; for 1 <= i <= min(M,N), row i of the
 *         tile U was interchanged with row IPIV(i) of the tile A.
 *
 * @param[in,out] WORK
 *
 * @param[in] LDWORK
 *         The leading dimension of the array WORK.
 *
 * @param[out] INFO
 *
 *******************************************************************************
 *
 * @return
 *         \retval PLASMA_SUCCESS successful exit
 *         \retval <0 if INFO = -k, the k-th argument had an illegal value
 *         \retval >0 if INFO = k, U(k,k) is exactly zero. The factorization
 *              has been completed, but the factor U is exactly
 *              singular, and division by zero will occur if it is used
 *              to solve a system of equations.
 *
 ******************************************************************************/

#if defined(PLASMA_HAVE_WEAK)
#pragma weak CORE_ztstrf = PCORE_ztstrf
#define CORE_ztstrf PCORE_ztstrf
#define CORE_zssssm PCORE_zssssm
int  CORE_zssssm(int M1, int N1, int M2, int N2, int K, int IB,
                 PLASMA_Complex64_t *A1, int LDA1,
                 PLASMA_Complex64_t *A2, int LDA2,
                 const PLASMA_Complex64_t *L1, int LDL1,
                 const PLASMA_Complex64_t *L2, int LDL2,
                 const int *IPIV);
#endif
int CORE_ztstrf(int M, int N, int IB, int NB,
                PLASMA_Complex64_t *U, int LDU,
                PLASMA_Complex64_t *A, int LDA,
                PLASMA_Complex64_t *L, int LDL,
                int *IPIV,
                PLASMA_Complex64_t *WORK, int LDWORK,
                int *INFO)
{
    static PLASMA_Complex64_t zzero = 0.0;
    static PLASMA_Complex64_t mzone =-1.0;

    PLASMA_Complex64_t alpha;
    int i, j, ii, sb;
    int im, ip;

    /* Check input arguments */
    *INFO = 0;
    if (M < 0) {
        coreblas_error(1, "Illegal value of M");
        return -1;
    }
    if (N < 0) {
        coreblas_error(2, "Illegal value of N");
        return -2;
    }
    if (IB < 0) {
        coreblas_error(3, "Illegal value of IB");
        return -3;
    }
    if ((LDU < max(1,NB)) && (NB > 0)) {
        coreblas_error(6, "Illegal value of LDU");
        return -6;
    }
    if ((LDA < max(1,M)) && (M > 0)) {
        coreblas_error(8, "Illegal value of LDA");
        return -8;
    }
    if ((LDL < max(1,IB)) && (IB > 0)) {
        coreblas_error(10, "Illegal value of LDL");
        return -10;
    }

    /* Quick return */
    if ((M == 0) || (N == 0) || (IB == 0))
        return PLASMA_SUCCESS;

    /* Set L to 0 */
    memset(L, 0, LDL*N*sizeof(PLASMA_Complex64_t));

    ip = 0;
    for (ii = 0; ii < N; ii += IB) {
        sb = min(N-ii, IB);

        for (i = 0; i < sb; i++) {
            im = cblas_izamax(M, &A[LDA*(ii+i)], 1);
            IPIV[ip] = ii+i+1;

            if (cabs(A[LDA*(ii+i)+im]) > cabs(U[LDU*(ii+i)+ii+i])) {
                /*
                 * Swap behind.
                 */
                cblas_zswap(i, &L[LDL*ii+i], LDL, &WORK[im], LDWORK );
                /*
                 * Swap ahead.
                 */
                cblas_zswap(sb-i, &U[LDU*(ii+i)+ii+i], LDU, &A[LDA*(ii+i)+im], LDA );
                /*
                 * Set IPIV.
                 */
                IPIV[ip] = NB + im + 1;

                for (j = 0; j < i; j++) {
                    A[LDA*(ii+j)+im] = zzero;
                }
            }

            if ((*INFO == 0) && (cabs(A[LDA*(ii+i)+im]) == zzero)
                && (cabs(U[LDU*(ii+i)+ii+i]) == zzero)) {
                *INFO = ii+i+1;
            }

            alpha = ((PLASMA_Complex64_t)1. / U[LDU*(ii+i)+ii+i]);
            cblas_zscal(M, CBLAS_SADDR(alpha), &A[LDA*(ii+i)], 1);
            cblas_zcopy(M, &A[LDA*(ii+i)], 1, &WORK[LDWORK*i], 1);
            cblas_zgeru(
                CblasColMajor, M, sb-i-1,
                CBLAS_SADDR(mzone), &A[LDA*(ii+i)], 1,
                &U[LDU*(ii+i+1)+ii+i], LDU,
                &A[LDA*(ii+i+1)], LDA );
            ip = ip+1;
        }
        /*
         * Apply the subpanel to the rest of the panel.
         */
        if(ii+i < N) {
            for(j = ii; j < ii+sb; j++) {
                if (IPIV[j] <= NB) {
                    IPIV[j] = IPIV[j] - ii;
                }
            }

            CORE_zssssm(
                NB, N-(ii+sb), M, N-(ii+sb), sb, sb,
                &U[LDU*(ii+sb)+ii], LDU,
                &A[LDA*(ii+sb)], LDA,
                &L[LDL*ii], LDL,
                WORK, LDWORK, &IPIV[ii]);

            for(j = ii; j < ii+sb; j++) {
                if (IPIV[j] <= NB) {
                    IPIV[j] = IPIV[j] + ii;
                }
            }
        }
    }
    return PLASMA_SUCCESS;
}


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

#ifndef COMANCHE_HSTORE_POOL_MANAGER_H
#define COMANCHE_HSTORE_POOL_MANAGER_H

#include <api/kvstore_itf.h> /* status_t */

#include "alloc_key.h"

#include <common/logging.h> /* log_source */
#include <nupm/region_descriptor.h>
#include <gsl/pointers>
#include <sys/uio.h>
#include <cstddef>
#include <functional>
#include <string>
#include <system_error>

struct pool_path;

#pragma GCC diagnostic push
// #pragma GCC diagnostic ignored "-Wunused-parameter"

enum class pool_ec
{
  pool_fail,
  pool_unsupported_mode,
  region_fail,
  region_fail_general_exception,
  region_fail_api_exception,
};

struct pool_category
  : public std::error_category
{
  const char* name() const noexcept override { return "pool_category"; }
  std::string message( int condition ) const noexcept override
  {
    switch ( condition )
    {
    case int(pool_ec::pool_fail):
      return "default pool failure";
    case int(pool_ec::pool_unsupported_mode):
      return "pool unsupported flags";
    case int(pool_ec::region_fail):
      return "region-backed pool failure";
    case int(pool_ec::region_fail_general_exception):
      return "region-backed pool failure (General_exception)";
    case int(pool_ec::region_fail_api_exception):
      return "region-backed pool failure (API_exception)";
    default:
      return "unknown pool failure";
    }
  }
};

namespace
{
  pool_category pool_error_category;
}

struct pool_error
  : public std::error_condition
{
  std::string _msg;
public:
  pool_error(const std::string &msg_, pool_ec val_)
    : std::error_condition(int(val_), pool_error_category)
    , _msg(msg_)
  {}

};

struct dax_manager;

template <typename Pool>
  struct pool_manager
    : protected common::log_source
  {
    pool_manager(unsigned debug_level_) : common::log_source(debug_level_) {}
    virtual ~pool_manager() {}

    virtual void pool_create_check(const std::size_t size_) = 0;

    virtual void pool_close_check(const std::string &) = 0;

    virtual nupm::region_descriptor pool_get_regions(const Pool &) const = 0;

    /*
     * throws pool_error if create_region fails
     */
    virtual auto pool_create_1(
      const pool_path &path_
      , std::size_t size_
    ) -> nupm::region_descriptor = 0;

    virtual auto pool_create_2(
      AK_FORMAL
      const nupm::region_descriptor & rac
      , component::IKVStore::flags_t flags
      , std::size_t expected_obj_count
    ) -> std::unique_ptr<Pool> = 0;

    virtual auto pool_open_1(
      const pool_path &path_
    ) -> nupm::region_descriptor = 0;

    virtual auto pool_open_2(
      AK_FORMAL
      const nupm::region_descriptor & access_
      , component::IKVStore::flags_t flags_
    ) -> std::unique_ptr<Pool> = 0;

    virtual void pool_delete(const pool_path &path) = 0;
    virtual const std::unique_ptr<dax_manager> & get_dax_manager() const = 0;
  };
#pragma GCC diagnostic pop

#endif


text:
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>
#include <gsl/gsl_qrng.h>
#include <gsl/gsl_histogram2d.h>
#include <gsl/gsl_sf_laguerre.h>

#include <courier.h>
#include "data.h"

static const char* localhost = "127.0.0.1";
extern void splineval_init(int* nx, double *x, double *fcn_1d);
extern void spline_dump(int* nx);
extern void throw_darts(const int* ndarts, const int *nx, double *x, double *y, int* hits);
extern void classify(const int *nx, double *x, double *y, int* hits);

void report(size_t ndarts, double* x, double* y, int* hits, Point* rowvec) {
    int i;
}
void explore(const int ndarts, size_t nreps, size_t nx, size_t ny,
        double xmin, double xmax, double ymin, double ymax) {
    double u, v;
    int i, j;
    const char* deps[] = {"Point"};
    const gsl_rng_type * T = gsl_rng_default;
    gsl_rng * r = gsl_rng_alloc (T);
    gsl_histogram2d *h = gsl_histogram2d_alloc(nx,ny);
    gsl_histogram2d_pdf *p = gsl_histogram2d_pdf_alloc(nx,ny);

    gsl_histogram2d_set_ranges_uniform(h, xmin, xmax, ymin, ymax);
    gsl_histogram2d_shift(h, 100.0);

    double *x = malloc(sizeof(double)*ndarts),
           *y = malloc(sizeof(double)*ndarts);
    int *c = malloc(sizeof(int)*ndarts);
    Point *rowvec = malloc(sizeof(Point)*ndarts);
    int breakout = 0;
    for (j = 0; j < nreps; j++) {
        gsl_histogram2d_pdf_init(p,h);

        for (i = 0; i < ndarts; i++) {
            u = gsl_rng_uniform(r);
            v = gsl_rng_uniform(r);
            gsl_histogram2d_pdf_sample (p, u, v, &x[i], &y[i]);
        }
        classify(&ndarts,x,y,c);
        for (i = 0; i < ndarts; i++) {
            rowvec[i].x = x[i];
            rowvec[i].y = y[i];
            rowvec[i].category = c[i];
        }
        courier_write_table("Point", rowvec, sizeof(Point), ndarts, 0);
        char* res = courier_request_analysis(localhost, "Update", deps, 1);
        {
            char* token = strtok(res, " ");
            while (token != NULL)
            {
                size_t idx = strtod(token, NULL) - 1;
                token = strtok(NULL, " ");
                double w = strtod(token, NULL);
                h->bin[idx] = 1;//fmax(0,10 / w);
                token = strtok(NULL, " ");
            }
            if (gsl_histogram2d_sum(h) <= nx*ny) {
                printf("Stopping due to exhausted histogram\n");
                breakout = 1;
            }
        }
        courier_free_data(res);
        if (breakout) {
            printf("Iteration %d\n",j);
            break;
        }
    }
    FILE *hf = fopen("hist.tab", "w");
    gsl_histogram2d_fprintf(hf,h,"%g","%g");
    fclose(hf);

    free(rowvec);
    free(c);
    free(x);
    free(y);
    gsl_rng_free (r);
    gsl_histogram2d_pdf_free(p);
    gsl_histogram2d_free(h);
}
int main(int argc, char** argv) {
    int nx = 100;
    int i, j;
    double x[nx], fcn_1d[nx];
    const gsl_rng_type * T;
    gsl_rng * r;

    gsl_rng_env_setup();

    T = gsl_rng_default;
    r = gsl_rng_alloc (T);
    fcn_1d[0] = 0;
    x[0] = 0;
    for (i = 1; i < nx; i++)
    {
        x[i] = i + 1.0;
        double h = gsl_ran_levy_skew(r, 1, 1.0, 1.0);
        if (h < 1 && h > -1) { h = -5; }
        fcn_1d[i] = fmax(25,fmin(75,fcn_1d[i-1]+h));
    }
    splineval_init(&nx,x,fcn_1d);
    explore(1000, 1000, 64, 64, 0, 100, 0, 100);
    spline_dump(&nx);
    gsl_rng_free (r);
    return 0;
}

max_stars_repo_name: sashalemon/spark-coupling

max_issues_repo_name: sashalemon/spark-coupling

max_forks_repo_name: sashalemon/spark-coupling


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

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

#include "ccl.h"

// Global variable to hold the neutrino phase-space spline
gsl_spline* nu_spline = NULL;

// these are NOT adjustable
// this phase space integral is only done once and the following is accurate
// enough according to tests done by the devs
/**
 * Absolute precision in neutrino root finding
 */
#define GSL_EPSABS_NU 1E-7

/**
 * Relative precision in neutrino root finding
 */
#define GSL_EPSREL_NU 1E-7

/**
 * Number of iterations for neutrino root finding
 */
#define GSL_N_ITERATION_NU 1000

/* ------- ROUTINE: nu_integrand ------
INPUTS: x: dimensionless momentum, *r: pointer to a dimensionless mass / temperature
TASK: Integrand of phase-space massive neutrino integral
*/
static double nu_integrand(double x, void *r) {
  double rat = *((double*)(r));
  double x2 = x*x;
  return sqrt(x2 + rat*rat) / (exp(x)+1.0) * x2;
}

/* ------- ROUTINE: ccl_calculate_nu_phasespace_spline ------
TASK: Get the spline of the result of the phase-space integral required for massive neutrinos.
*/

gsl_spline* calculate_nu_phasespace_spline(int *status) {
  int N = CCL_NU_MNUT_N;
  double *mnut = NULL;
  double *y = NULL;
  gsl_spline* spl = NULL;
  gsl_integration_cquad_workspace * workspace = NULL;
  int stat = 0, gslstatus;
  gsl_function F;

  mnut = ccl_linear_spacing(log(CCL_NU_MNUT_MIN), log(CCL_NU_MNUT_MAX), N);
  y = malloc(sizeof(double)*CCL_NU_MNUT_N);
  if ((y == NULL) || (mnut == NULL)) {
    // Not setting a status_message here because we can't easily pass a
    // cosmology to this function - message printed in ccl_error.c.
    *status = CCL_ERROR_NU_INT;
  }

  if (*status == 0) {
    workspace = gsl_integration_cquad_workspace_alloc(GSL_N_ITERATION_NU);
    if (workspace == NULL)
      *status = CCL_ERROR_NU_INT;
  }

  if (*status == 0) {
    F.function = &nu_integrand;
    for (int i=0; i < CCL_NU_MNUT_N; i++) {
      double mnut_ = exp(mnut[i]);
      F.params = &(mnut_);
      gslstatus = gsl_integration_cquad(&F, 0, 1000.0,
                                        GSL_EPSABS_NU,
                                        GSL_EPSREL_NU,
                                        workspace, &y[i], NULL, NULL);
      if (gslstatus != GSL_SUCCESS) {
        ccl_raise_gsl_warning(gslstatus, "ccl_neutrinos.c: calculate_nu_phasespace_spline():");
        stat |= gslstatus;
      }
    }

    double renorm = 1./y[0];
    for (int i=0; i < CCL_NU_MNUT_N; i++)
      y[i] *= renorm;

    if (stat) {
      *status = CCL_ERROR_NU_INT;
    }
  }

  if (*status == 0) {
    spl = gsl_spline_alloc(gsl_interp_akima, CCL_NU_MNUT_N);
    if (spl == NULL)
      *status = CCL_ERROR_NU_INT;
  }

  if (*status == 0) {
    stat |= gsl_spline_init(spl, mnut, y, CCL_NU_MNUT_N);
    if (stat) {
      ccl_raise_gsl_warning(gslstatus, "ccl_neutrinos.c: calculate_nu_phasespace_spline():");
      *status = CCL_ERROR_NU_INT;
    }
  }

  // Check for errors in creating the spline
  if (stat || (*status)) {
    // Not setting a status_message here because we can't easily pass a
    // cosmology to this function - message printed in ccl_error.c.
    *status = CCL_ERROR_NU_INT;
    gsl_spline_free(spl);
  }

  gsl_integration_cquad_workspace_free(workspace);
  free(mnut);
  free(y);

  return spl;
}

/* ------- ROUTINE: ccl_nu_phasespace_intg ------
INPUTS: mnuOT: the dimensionless mass / temperature of a single massive neutrino
TASK: Get the value of the phase space integral at mnuOT
*/
double nu_phasespace_intg(double mnuOT, int* status)
{
  // Check if the global variable for the phasespace spline has been defined yet:
  if (nu_spline == NULL)
    nu_spline = calculate_nu_phasespace_spline(status);

  if (*status) {
    return NAN;
  }

  double integral_value = 0.;

  // First check the cases where we are in the limits.
  if (mnuOT < CCL_NU_MNUT_MIN)
    return 7./8.;
  else if (mnuOT > CCL_NU_MNUT_MAX)
    return 0.2776566337 * mnuOT;

  int gslstatus = gsl_spline_eval_e(nu_spline, log(mnuOT), NULL, &integral_value);
  if (gslstatus != GSL_SUCCESS) {
    ccl_raise_gsl_warning(gslstatus, "ccl_neutrinos.c: nu_phasespace_intg():");
    *status |= gslstatus;
  }
  return integral_value * 7./8.;
}

/* -------- ROUTINE: Omeganuh2 ---------
INPUTS: a: scale factor, Nnumass: number of massive neutrino species,
        mnu: total mass in eV of neutrinos, T_CMB: CMB temperature, status: pointer to status integer.
TASK: Compute Omeganu * h^2 as a function of time.
!! To all practical purposes, Neff is simply N_nu_mass !!
*/
double ccl_Omeganuh2 (double a, int N_nu_mass, double* mnu, double T_CMB, int* status) {
  double Tnu, a4, prefix_massless, mnuone, OmNuh2;
  double Tnu_eff, mnuOT, intval, prefix_massive;
  double total_mass; // To check if this is the massless or massive case.

  // First check if N_nu_mass is 0
  if (N_nu_mass == 0) return 0.0;

  Tnu = T_CMB*pow(4./11.,1./3.);
  a4 = a*a*a*a;
  // Check if mnu=0. We assume that in the massless case mnu is a
  // pointer to a single element and that element is 0. This should in principle never be called.
  if (mnu[0] < 0.00017) {  // Limit taken from Lesgourges et al. 2012
    prefix_massless = NU_CONST  * Tnu * Tnu * Tnu * Tnu;
    return N_nu_mass*prefix_massless*7./8./a4;
  }

  // And the remaining massive case.
  // Tnu_eff is used in the massive case because CLASS uses an effective
  // temperature of nonLCDM components to match to mnu / Omeganu =93.14eV. Tnu_eff = T_ncdm * T_CMB = 0.71611 * T_CMB
  Tnu_eff = Tnu * ccl_constants.TNCDM / (pow(4./11.,1./3.));

  // Define the prefix using the effective temperature (to get mnu / Omega = 93.14 eV) for the massive case:
  prefix_massive = NU_CONST * Tnu_eff * Tnu_eff * Tnu_eff * Tnu_eff;

  OmNuh2 = 0.; // Initialize to 0 - we add to this for each massive neutrino species.
  for(int i=0; i < N_nu_mass; i++) {
    // Get mass over T (mass (eV) / ((kb eV/s/K) Tnu_eff (K))
    // This returns the density normalized so that we get nuh2 at a=0
    mnuOT = mnu[i] / (Tnu_eff/a) * (ccl_constants.EV_IN_J / (ccl_constants.KBOLTZ));

    // Get the value of the phase-space integral
    intval = nu_phasespace_intg(mnuOT, status);
    OmNuh2 = intval*prefix_massive/a4 + OmNuh2;
  }

  return OmNuh2;
}

/* -------- ROUTINE: Omeganuh2_to_Mnu ---------
INPUTS: OmNuh2: neutrino mass density today Omeganu * h^2,
        label: how you want to split up the masses, see ccl_neutrinos.h for options,
        T_CMB: CMB temperature, status: pointer to status integer.
TASK: Given Omeganuh2 today, the method of splitting into masses, and the temperature of the
      CMB, output a pointer to the array of neutrino masses (may be length 1 if label asks for sum)
*/
double* ccl_nu_masses(double OmNuh2, ccl_neutrino_mass_splits mass_split,
                      double T_CMB,  int* status) {
  double sumnu;
  double *mnu = NULL;

  sumnu = 93.14 * OmNuh2;

  // Now split the sum up into three masses depending on the label given:
  if (mass_split == ccl_nu_normal) {
    mnu = malloc(3*sizeof(double));
    if (mnu == NULL) {
      *status = CCL_ERROR_MEMORY;
    }

    if (*status == 0) {
      // See CCL note for how we get these expressions for the neutrino masses in
      // normal and inverted hierarchy.
      mnu[0] = (
        2./3.* sumnu - 1./6. * pow(-6. * ccl_constants.DELTAM12_sq + 12. * ccl_constants.DELTAM13_sq_pos + 4. * sumnu*sumnu, 0.5)
        - 0.25 * ccl_constants.DELTAM12_sq / (2./3.* sumnu - 1./6. * pow(-6. * ccl_constants.DELTAM12_sq + 12. *
        ccl_constants.DELTAM13_sq_pos + 4. * sumnu*sumnu, 0.5)));
      mnu[1] = (
        2./3.* sumnu - 1./6. * pow(-6. * ccl_constants.DELTAM12_sq + 12. * ccl_constants.DELTAM13_sq_pos + 4. * sumnu*sumnu, 0.5)
        + 0.25 * ccl_constants.DELTAM12_sq / (2./3.* sumnu - 1./6. * pow(-6. * ccl_constants.DELTAM12_sq + 12. *
        ccl_constants.DELTAM13_sq_pos + 4. * sumnu*sumnu, 0.5)));
      mnu[2] = (
        -1./3. * sumnu + 1./3 * pow(-6. * ccl_constants.DELTAM12_sq +
          12. * ccl_constants.DELTAM13_sq_pos + 4. * sumnu*sumnu, 0.5));

      if (mnu[0] < 0 || mnu[1] < 0 || mnu[2] < 0) {
        // The user has provided a sum that is below the physical limit.
        if (sumnu < 1e-14) {
          mnu[0] = 0.;
          mnu[1] = 0.;
          mnu[2] = 0.;
        }
        else {
          *status = CCL_ERROR_MNU_UNPHYSICAL;
          ccl_raise_warning(
              *status,
              "CCL_ERROR_MNU_UNPHYSICAL: Sum of neutrinos masses for this Omeganu "
              "value is incompatible with the requested mass hierarchy.");
        }
      }
    }
  }
  else if (mass_split == ccl_nu_inverted) {
    mnu = malloc(3*sizeof(double));
    if (mnu == NULL) {
      *status = CCL_ERROR_MEMORY;
    }

    if (*status == 0) {
      mnu[0] = (
        2./3.* sumnu - 1./6. * pow(-6. * ccl_constants.DELTAM12_sq + 12. * ccl_constants.DELTAM13_sq_neg + 4. * sumnu * sumnu, 0.5)
        - 0.25 * ccl_constants.DELTAM12_sq / (2./3.* sumnu - 1./6. * pow(-6. * ccl_constants.DELTAM12_sq + 12. *
        ccl_constants.DELTAM13_sq_neg + 4. * sumnu * sumnu, 0.5)));
      mnu[1] = (
        2./3.* sumnu - 1./6. * pow(-6. * ccl_constants.DELTAM12_sq + 12. * ccl_constants.DELTAM13_sq_neg + 4. * sumnu*sumnu, 0.5)
        + 0.25 * ccl_constants.DELTAM12_sq / (2./3.* sumnu - 1./6. * pow(-6. *
        ccl_constants.DELTAM12_sq + 12. * ccl_constants.DELTAM13_sq_neg + 4. * sumnu*sumnu, 0.5)));
      mnu[2] = (
        -1./3. * sumnu + 1./3 * pow(-6. * ccl_constants.DELTAM12_sq +
          12. * ccl_constants.DELTAM13_sq_neg + 4. * sumnu*sumnu, 0.5));

      if (mnu[0] < 0 || mnu[1] < 0 || mnu[2] < 0) {
        // The user has provided a sum that is below the physical limit.
        if (sumnu < 1e-14) {
          mnu[0] = 0.;
          mnu[1] = 0.;
          mnu[2] = 0.;;
        }
        else {
          *status = CCL_ERROR_MNU_UNPHYSICAL;
          ccl_raise_warning(
              *status,
              "CCL_ERROR_MNU_UNPHYSICAL: Sum of neutrinos masses for this Omeganu "
              "value is incompatible with the requested mass hierarchy.");
        }
      }
    }
  }
  else if (mass_split == ccl_nu_equal) {
    mnu = malloc(3*sizeof(double));
    if (mnu == NULL) {
      *status = CCL_ERROR_MEMORY;
    }

    if (*status == 0) {
      mnu[0] = sumnu/3.;
      mnu[1] = sumnu/3.;
      mnu[2] = sumnu/3.;
    }
  }
  else if (mass_split == ccl_nu_sum) {
    mnu = malloc(sizeof(double));
    if (mnu == NULL) {
      *status = CCL_ERROR_MEMORY;
    }

    if (*status == 0) {
      mnu[0] = sumnu;
    }
  }
  else {
    *status = CCL_ERROR_MNU_UNPHYSICAL;
    ccl_raise_warning(
        *status,
        "mass option = %d not yet supported!", mass_split);
  }

  return mnu;
}

#undef GSL_EPSABS_NU
#undef GSL_EPSREL_NU
#undef GSL_N_ITERATION_NU

max_stars_repo_name: benediktdiemer/CCL

max_issues_repo_name: benediktdiemer/CCL

max_forks_repo_name: benediktdiemer/CCL


text:
#include <stdio.h>
#include <gsl/gsl_integration.h>
#include "function.h"

extern double f (double, void*);

int
main (void)
{
  /* Claim memory needed for adaptive integration */
  gsl_integration_workspace * w
    = gsl_integration_workspace_alloc (1000);

  /* Construct GSL function object to describe solution */
  gsl_function F;
  F.function = &f;
  double alpha = 1.0;
  F.params = NULL;

  /* Parameters to control quadrature routine */
  double a, b, tol_abs, tol_rel;
  int limit;
  a = 0.0;
  b = 3.0;
  tol_abs = 0;
  tol_rel = 1.0e-7;
  limit = 1000;

  // Test that the function is implemented properly
  //double test = GSL_FN_EVAL(&F,2);
  //printf("TEST RESULT: %.18f\n", test);

  /* Actual call to quadrature routine */
  double result, error;
  gsl_integration_qag (&F, a, b, tol_abs, tol_rel, limit,
                       GSL_INTEG_GAUSS15, w, &result, &error);
  double mathematica_result = 1.29646778572437307899;

  printf ("result                = % .18f\n", result);
  printf ("mathematica result    = % .18f\n", mathematica_result);
  printf ("estimated error       = % .18f\n", error);
  printf ("error vs. mathematica = % .18f\n",
           fabs(result - mathematica_result));
  printf ("# intervals     = %d\n", (int) w->size);

  /* Return memory allocated for solver */
  gsl_integration_workspace_free (w);

  return 0;
}


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

#include <config.h>
#include <stdlib.h>
#include <math.h>
#include <string.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_rstat.h>
#include <gsl/gsl_sort.h>
#include <gsl/gsl_statistics.h>
#include <gsl/gsl_test.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>
#include <gsl/gsl_ieee_utils.h>

double *
random_data(const size_t n, gsl_rng *r)
{
  size_t i;
  double *data = malloc(n * sizeof(double));

  for (i = 0; i < n; ++i)
    data[i] = gsl_rng_uniform(r);

  return data;
}

void
test_basic(const size_t n, const double data[], const double tol)
{
  gsl_rstat_workspace *rstat_workspace_p = gsl_rstat_alloc();
  const double expected_mean = gsl_stats_mean(data, 1, n);
  const double expected_var = gsl_stats_variance(data, 1, n);
  const double expected_sd = gsl_stats_sd(data, 1, n);
  const double expected_sd_mean = expected_sd / sqrt((double) n);
  const double expected_skew = gsl_stats_skew(data, 1, n);
  const double expected_kurtosis = gsl_stats_kurtosis(data, 1, n);
  double expected_rms = 0.0;
  double mean, var, sd, sd_mean, rms, skew, kurtosis;
  size_t i, num;
  int status;

  /* compute expected rms */
  for (i = 0; i < n; ++i)
    expected_rms += data[i] * data[i];

  expected_rms = sqrt(expected_rms / n);

  /* add data to rstat workspace */
  for (i = 0; i < n; ++i)
    gsl_rstat_add(data[i], rstat_workspace_p);

  mean     = gsl_rstat_mean(rstat_workspace_p);
  var      = gsl_rstat_variance(rstat_workspace_p);
  sd       = gsl_rstat_sd(rstat_workspace_p);
  sd_mean  = gsl_rstat_sd_mean(rstat_workspace_p);
  rms      = gsl_rstat_rms(rstat_workspace_p);
  skew     = gsl_rstat_skew(rstat_workspace_p);
  kurtosis = gsl_rstat_kurtosis(rstat_workspace_p);
  num      = gsl_rstat_n(rstat_workspace_p);

  gsl_test_int(num, n, "n n=%zu" , n);
  gsl_test_rel(mean, expected_mean, tol, "mean n=%zu", n);
  gsl_test_rel(var, expected_var, tol, "variance n=%zu", n);
  gsl_test_rel(sd, expected_sd, tol, "stddev n=%zu", n);
  gsl_test_rel(sd_mean, expected_sd_mean, tol, "stddev_mean n=%zu", n);
  gsl_test_rel(rms, expected_rms, tol, "rms n=%zu", n);
  gsl_test_rel(skew, expected_skew, tol, "skew n=%zu", n);
  gsl_test_rel(kurtosis, expected_kurtosis, tol, "kurtosis n=%zu", n);

  status = gsl_rstat_reset(rstat_workspace_p);
  gsl_test_int(status, GSL_SUCCESS, "rstat returned success");
  num = gsl_rstat_n(rstat_workspace_p);

  gsl_test_int(num, 0, "n n=%zu" , n);

  gsl_rstat_free(rstat_workspace_p);
}

void
test_quantile(const double p, const double data[], const size_t n,
              const double expected, const double tol, const char *desc)
{
  gsl_rstat_quantile_workspace *w = gsl_rstat_quantile_alloc(p);
  double result;
  size_t i;

  for (i = 0; i < n; ++i)
    gsl_rstat_quantile_add(data[i], w);

  result = gsl_rstat_quantile_get(w);

  if (fabs(expected) < 1.0e-4)
    gsl_test_abs(result, expected, tol, "%s p=%g", desc, p);
  else
    gsl_test_rel(result, expected, tol, "%s p=%g", desc, p);

  gsl_rstat_quantile_free(w);
}

int
main()
{
  gsl_rng *r = gsl_rng_alloc(gsl_rng_default);
  const double tol1 = 1.0e-8;
  const double tol2 = 1.0e-3;

  gsl_ieee_env_setup();

  {
    const size_t N = 2000000;
    double *data = random_data(N, r);
    double data2[] = { 4.0, 7.0, 13.0, 16.0, -5.0 };
    size_t i;

    test_basic(2, data, tol1);
    test_basic(100, data, tol1);
    test_basic(1000, data, tol1);
    test_basic(10000, data, tol1);
    test_basic(50000, data, tol1);
    test_basic(80000, data, tol1);
    test_basic(1500000, data, tol1);
    test_basic(2000000, data, tol1);

    for (i = 0; i < 5; ++i)
      data2[i] += 1.0e9;

    test_basic(5, data2, tol1);

    free(data);
  }

  {
    /* dataset from Jain and Chlamtac paper */
    const size_t n_jain = 20;
    const double data_jain[] = {  0.02,  0.15,  0.74,  3.39,  0.83,
                                  22.37, 10.15, 15.43, 38.62, 15.92,
                                  34.60, 10.28,  1.47,  0.40,  0.05,
                                  11.39,  0.27,  0.42,  0.09, 11.37 };
    double expected_jain = 4.44063435326;

    test_quantile(0.5, data_jain, n_jain, expected_jain, tol1, "jain");
  }

  {
    size_t n = 1000000;
    double *data = malloc(n * sizeof(double));
    double *sorted_data = malloc(n * sizeof(double));
    gsl_rstat_workspace *rstat_workspace_p = gsl_rstat_alloc();
    double p;
    size_t i;

    for (i = 0; i < n; ++i)
      {
        data[i] = gsl_ran_gaussian_tail(r, 1.3, 1.0);
        gsl_rstat_add(data[i], rstat_workspace_p);
      }

    memcpy(sorted_data, data, n * sizeof(double));
    gsl_sort(sorted_data, 1, n);

    /* test quantile calculation */
    for (p = 0.1; p <= 0.9; p += 0.1)
      {
        double expected = gsl_stats_quantile_from_sorted_data(sorted_data, 1, n, p);
        test_quantile(p, data, n, expected, tol2, "gauss");
      }

    /* test mean, variance */
    {
      const double expected_mean = gsl_stats_mean(data, 1, n);
      const double expected_var = gsl_stats_variance(data, 1, n);
      const double expected_sd = gsl_stats_sd(data, 1, n);
      const double expected_skew = gsl_stats_skew(data, 1, n);
      const double expected_kurtosis = gsl_stats_kurtosis(data, 1, n);
      const double expected_median = gsl_stats_quantile_from_sorted_data(sorted_data, 1, n, 0.5);

      const double mean = gsl_rstat_mean(rstat_workspace_p);
      const double var = gsl_rstat_variance(rstat_workspace_p);
      const double sd = gsl_rstat_sd(rstat_workspace_p);
      const double skew = gsl_rstat_skew(rstat_workspace_p);
      const double kurtosis = gsl_rstat_kurtosis(rstat_workspace_p);
      const double median = gsl_rstat_median(rstat_workspace_p);

      gsl_test_rel(mean, expected_mean, tol1, "mean");
      gsl_test_rel(var, expected_var, tol1, "variance");
      gsl_test_rel(sd, expected_sd, tol1, "stddev");
      gsl_test_rel(skew, expected_skew, tol1, "skew");
      gsl_test_rel(kurtosis, expected_kurtosis, tol1, "kurtosis");
      gsl_test_abs(median, expected_median, tol2, "median");
    }

    free(data);
    free(sorted_data);
    gsl_rstat_free(rstat_workspace_p);
  }

  gsl_rng_free(r);

  exit (gsl_test_summary());
}

max_forks_repo_name: karanbirsandhu/nu-sense


text:
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <assert.h>
#include <time.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>
#include <gsl/gsl_math.h>

typedef struct parameters{
	double Xo;
	double Ko;
	double SIGMA2_K;
	double SIGMA2_C;
	double R;
	double MU;
	double SIGMA2_MU;
}PAR;

/****** Function Prototypes ******/
/* Model */

double Nhat(double x, double Ko, double SIGMA2_K);
double C(double x, double y, double SIGMA2_C);

/* Jump moment equations */
double alpha1(double phi, void *params);
double alpha1_prime(double phi, void *params);

/* Useful functions */
double round( double x);               // rounding

/*Main function */
void canonical
(
	PAR parameters,
	double MAX_TIME,
	double ENSEMBLES,
	double Dt,
	double dt
);

void ratescalc
(
	gsl_rng *rng,
	double X, double Y,
	void *params,
	double *prob, double *wait_time
);

int func (double t, const double y[], double f[], void *params);
void gslode(void *params, double MAX_TIME, FILE *theory);
void euler(void *params, double MAX_TIME, FILE *theory);

max_stars_repo_name: kbroman/fluctuationDomains

max_issues_repo_name: kbroman/fluctuationDomains

max_forks_repo_name: kbroman/fluctuationDomains


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

/* parsec things */
#include "parsec.h"

/* system and io */
#include <stdlib.h>
#include <stdio.h>
#include <string.h>
/* Plasma and math libs */
#include <math.h>
//#include <cblas.h>
//#include <lapacke.h>
//#include <core_blas.h>

#include "parsec/profiling.h"
#include "parsec/parsec_internal.h"
#include "parsec/utils/debug.h"
//#include "dplasma.h"

#ifdef __cplusplus
extern "C" {
#endif

/* these are globals in common.c */
extern const char *PARSEC_SCHED_NAME[];
extern int unix_timestamp;
extern char cwd[];

/* Update PASTE_CODE_PROGRESS_KERNEL below if you change this list */
enum iparam_t {
  IPARAM_RANK,         /* Rank                              */
  IPARAM_NNODES,       /* Number of nodes                   */
  IPARAM_NCORES,       /* Number of cores                   */
  IPARAM_NGPUS,        /* Number of GPUs                    */
  IPARAM_P,            /* Rows in the process grid          */
  IPARAM_Q,            /* Columns in the process grid       */
  IPARAM_M,            /* Number of rows of the matrix      */
  IPARAM_N,            /* Number of columns of the matrix   */
  IPARAM_K,            /* RHS or K                          */
  IPARAM_LDA,          /* Leading dimension of A            */
  IPARAM_LDB,          /* Leading dimension of B            */
  IPARAM_LDC,          /* Leading dimension of C            */
  IPARAM_IB,           /* Inner-blocking size               */
  IPARAM_NB,           /* Number of columns in a tile       */
  IPARAM_MB,           /* Number of rows in a tile          */
  IPARAM_SNB,          /* Number of columns in a super-tile */
  IPARAM_SMB,          /* Number of rows in a super-tile    */
  IPARAM_HMB,          /* Small MB for recursive hdags */
  IPARAM_HNB,          /* Small NB for recursive hdags */
  IPARAM_CHECK,        /* Checking activated or not         */
  IPARAM_VERBOSE,      /* How much noise do we want?        */
  IPARAM_SCHEDULER,    /* User-selected scheduler */
  IPARAM_SIZEOF
};

#define PARSEC_SCHEDULER_DEFAULT 0
#define PARSEC_SCHEDULER_LFQ 1
#define PARSEC_SCHEDULER_LTQ 2
#define PARSEC_SCHEDULER_AP  3
#define PARSEC_SCHEDULER_LHQ 4
#define PARSEC_SCHEDULER_GD  5
#define PARSEC_SCHEDULER_PBQ 6
#define PARSEC_SCHEDULER_IP  7
#define PARSEC_SCHEDULER_RND 8

void iparam_default_gemm(int* iparam);
void iparam_default_ibnbmb(int* iparam, int ib, int nb, int mb);

#define PASTE_CODE_IPARAM_LOCALS(iparam)                                \
    rank  = iparam[IPARAM_RANK];                                    \
    nodes = iparam[IPARAM_NNODES];                                  \
    cores = iparam[IPARAM_NCORES];                                  \
    gpus  = iparam[IPARAM_NGPUS];                                   \
    P     = iparam[IPARAM_P];                                       \
    Q     = iparam[IPARAM_Q];                                       \
    M     = iparam[IPARAM_M];                                       \
    N     = iparam[IPARAM_N];                                       \
    K     = iparam[IPARAM_K];                                       \
    NRHS  = K;                                                      \
    LDA   = max(M, iparam[IPARAM_LDA]);                             \
    LDB   = max(N, iparam[IPARAM_LDB]);                             \
    LDC   = max(K, iparam[IPARAM_LDC]);                             \
    IB    = iparam[IPARAM_IB];                                      \
    MB    = iparam[IPARAM_MB];                                      \
    NB    = iparam[IPARAM_NB];                                      \
    SMB   = iparam[IPARAM_SMB];                                     \
    SNB   = iparam[IPARAM_SNB];                                     \
    HMB   = iparam[IPARAM_HMB];                                     \
    HNB   = iparam[IPARAM_HNB];                                     \
    MT    = (M%MB==0) ? (M/MB) : (M/MB+1);                          \
    NT    = (N%NB==0) ? (N/NB) : (N/NB+1);                          \
    KT    = (K%MB==0) ? (K/MB) : (K/MB+1);                          \
    check = iparam[IPARAM_CHECK];                                   \
    loud  = iparam[IPARAM_VERBOSE];                                 \
    scheduler = iparam[IPARAM_SCHEDULER];                           \
    (void)rank;(void)nodes;(void)cores;(void)gpus;(void)P;(void)Q;(void)M;(void)N;(void)K;(void)NRHS; \
    (void)LDA;(void)LDB;(void)LDC;(void)IB;(void)MB;(void)NB;(void)MT;(void)NT;(void)KT; \
    (void)SMB;(void)SNB;(void)HMB;(void)HNB;(void)check;(void)loud; \
    (void)scheduler;

/*******************************
 * globals values
 *******************************/

#if defined(PARSEC_HAVE_MPI)
extern MPI_Datatype SYNCHRO;
#endif  /* PARSEC_HAVE_MPI */

void print_usage(void);
void print_arguments(int* iparam);

parsec_context_t *setup_parsec(int argc, char* argv[], int *iparam);
void cleanup_parsec(parsec_context_t* parsec, int *iparam);

/**
 * No macro with the name max or min is acceptable as there is
 * no way to correctly define them without borderline effects.
 */
#undef max
#undef min
static inline int max(int a, int b) { return a > b ? a : b; }
static inline int min(int a, int b) { return a < b ? a : b; }

/* Paste code to allocate a matrix in desc if cond_init is true */
#define PASTE_CODE_ALLOCATE_MATRIX(DC, COND, TYPE, INIT_PARAMS)      \
    if(COND) {                                                          \
        DC.mat = parsec_data_allocate((size_t)DC.super.nb_local_tiles * \
                                        (size_t)DC.super.bsiz *      \
                                        (size_t)parsec_datadist_getsizeoftype(DC.super.mtype)); \
        parsec_data_collection_set_key((parsec_data_collection_t*)&DC, #DC);          \
    }

#ifdef __cplusplus
}
#endif

#endif /* _TESTSCOMMON_H */


text:
/*
 * VGMTrans (c) 2002-2019
 * Licensed under the zlib license,
 * refer to the included LICENSE.txt file
 */
#pragma once

#include <cstdint>
#include <vector>
#include <map>
#include <string>
#include <cassert>
#include <climits>
#include <gsl-lite.hpp>

class RawFile;

constexpr auto PSF_TAG_SIG = "[TAG]";
constexpr auto PSF_TAG_SIG_LEN = 5;
constexpr auto PSF_STRIP_BUF_SIZE = 4096;

class PSFFile2 {
   public:
    PSFFile2(const RawFile &file);
    ~PSFFile2() = default;

    uint8_t version() const noexcept { return m_version; }
    const std::map<std::string, std::string> &tags() const noexcept { return m_tags; }
    const std::vector<char> &exe() const noexcept { return m_exe_data; }
    const std::vector<char> &reservedSection() const noexcept { return m_reserved_data; }

    template <typename T>
    T getExe(size_t ind) const {
        assert(ind + sizeof(T) < m_exe_data.size());

        T value = 0;
        for (size_t i = 0; i < sizeof(T); i++) {
            value |= (m_exe_data[ind + i] << (i * CHAR_BIT));
        }

        return value;
    }

    template <typename T>
    T getRes(size_t ind) const {
        assert(ind + sizeof(T) < m_reserved_data.size());

        T value = 0;
        for (size_t i = 0; i < sizeof(T); i++) {
            value |= (m_reserved_data[ind + i] << (i * CHAR_BIT));
        }

        return value;
    }

   private:
    uint8_t m_version;
    uint32_t m_exe_CRC;
    std::vector<char> m_exe_data;
    std::vector<char> m_reserved_data;
    std::map<std::string, std::string> m_tags;

    void parseTags(gsl::span<const char> data);
};

max_stars_repo_name: sykhro/vgmtrans-qt

max_issues_repo_name: sykhro/vgmtrans-qt

max_forks_repo_name: sykhro/vgmtrans-qt


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

#include <config.h>
#include <stdlib.h>
#include <string.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_ieee_utils.h>

static int
lookup_string (const char * p, int * precision, int * rounding,
               int * exception_mask) ;

int
gsl_ieee_read_mode_string (const char * description,
                           int * precision,
                           int * rounding,
                           int * exception_mask)
{
  char * start ;
  char * end;
  char * p;

  int precision_count = 0 ;
  int rounding_count = 0 ;
  int exception_count = 0 ;

  start = (char *) malloc(strlen(description) + 1) ;

  if (start == 0)
    {
      GSL_ERROR ("no memory to parse mode string", GSL_ENOMEM) ;
    }

  strcpy (start, description) ;

  p = start ;

  *precision = 0 ;
  *rounding = 0 ;
  *exception_mask = 0 ;

  do {
    int status ;
    int new_precision, new_rounding, new_exception ;

    end = strchr (p,',') ;

    if (end)
      {
        *end = '\0' ;
        do
          {
            end++ ;  /* skip over trailing whitespace */
          }
        while (*end == ' ' || *end == ',') ;
      }

    new_precision = 0 ;
    new_rounding = 0 ;
    new_exception = 0 ;

    status = lookup_string (p, &new_precision, &new_rounding, &new_exception) ;

    if (status)
      GSL_ERROR ("unrecognized GSL_IEEE_MODE string.\nValid settings are:\n\n"
                 "  single-precision double-precision extended-precision\n"
                 "  round-to-nearest round-down round-up round-to-zero\n"
                 "  mask-invalid mask-denormalized mask-division-by-zero\n"
                 "  mask-overflow mask-underflow mask-all\n"
                 "  trap-common trap-inexact\n"
                 "\n"
                 "separated by commas. "
                 "(e.g. GSL_IEEE_MODE=\"round-down,mask-underflow\")",
                 GSL_EINVAL) ;

    if (new_precision)
      {
        *precision = new_precision ;
        precision_count ++ ;
        if (precision_count > 1)
          GSL_ERROR ("attempted to set IEEE precision twice", GSL_EINVAL) ;
      }

    if (new_rounding)
      {
        *rounding = new_rounding ;
        rounding_count ++ ;
        if (rounding_count > 1)
          GSL_ERROR ("attempted to set IEEE rounding mode twice", GSL_EINVAL) ;
      }

    if (new_exception)
      {
        *exception_mask |= new_exception ;
        exception_count ++ ;
      }

    p = end ;

  } while (end && *p != '\0') ;

  free(start) ;

  return GSL_SUCCESS ;
}

static int
lookup_string (const char * p, int * precision, int * rounding,
               int * exception_mask)
{
  if (strcmp(p,"single-precision") == 0)
    {
      *precision = GSL_IEEE_SINGLE_PRECISION ;
    }
  else if (strcmp(p,"double-precision") == 0)
    {
      *precision = GSL_IEEE_DOUBLE_PRECISION ;
    }
  else if (strcmp(p,"extended-precision") == 0)
    {
      *precision = GSL_IEEE_EXTENDED_PRECISION ;
    }
  else if (strcmp(p,"round-to-nearest") == 0)
    {
      *rounding = GSL_IEEE_ROUND_TO_NEAREST ;
    }
  else if (strcmp(p,"round-down") == 0)
    {
      *rounding = GSL_IEEE_ROUND_DOWN ;
    }
  else if (strcmp(p,"round-up") == 0)
    {
      *rounding = GSL_IEEE_ROUND_UP ;
    }
  else if (strcmp(p,"round-to-zero") == 0)
    {
      *rounding = GSL_IEEE_ROUND_TO_ZERO ;
    }
  else if (strcmp(p,"mask-all") == 0)
    {
      *exception_mask = GSL_IEEE_MASK_ALL ;
    }
  else if (strcmp(p,"mask-invalid") == 0)
    {
      *exception_mask = GSL_IEEE_MASK_INVALID ;
    }
  else if (strcmp(p,"mask-denormalized") == 0)
    {
      *exception_mask = GSL_IEEE_MASK_DENORMALIZED ;
    }
  else if (strcmp(p,"mask-division-by-zero") == 0)
    {
      *exception_mask = GSL_IEEE_MASK_DIVISION_BY_ZERO ;
    }
  else if (strcmp(p,"mask-overflow") == 0)
    {
      *exception_mask = GSL_IEEE_MASK_OVERFLOW ;
    }
  else if (strcmp(p,"mask-underflow") == 0)
    {
      *exception_mask = GSL_IEEE_MASK_UNDERFLOW ;
    }
  else if (strcmp(p,"trap-inexact") == 0)
    {
      *exception_mask = GSL_IEEE_TRAP_INEXACT ;
    }
  else if (strcmp(p,"trap-common") == 0)
    {
      return 0 ;
    }
  else
    {
      return 1 ;
    }

  return 0 ;
}

max_stars_repo_name: ruslankuzmin/julia


text:
/* fft/gsl_fft_complex.h
 *
 * Copyright (C) 1996, 1997, 1998, 1999, 2000 Brian Gough
 *
 * This program is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation; either version 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., 675 Mass Ave, Cambridge, MA 02139, USA.
 */

#ifndef __GSL_FFT_COMPLEX_H__
#define __GSL_FFT_COMPLEX_H__

#include <stddef.h>

#include <gsl/gsl_math.h>
#include <gsl/gsl_complex.h>
#include <gsl/gsl_fft.h>

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

__BEGIN_DECLS

/*  Power of 2 routines  */

int gsl_fft_complex_radix2_forward (gsl_complex_packed_array data,
                                    const size_t stride,
                                    const size_t n);

int gsl_fft_complex_radix2_backward (gsl_complex_packed_array data,
                                     const size_t stride,
                                     const size_t n);

int gsl_fft_complex_radix2_inverse (gsl_complex_packed_array data,
                                    const size_t stride,
                                    const size_t n);

int gsl_fft_complex_radix2_transform (gsl_complex_packed_array data,
                                      const size_t stride,
                                      const size_t n,
                                      const gsl_fft_direction sign);

int gsl_fft_complex_radix2_dif_forward (gsl_complex_packed_array data,
                                        const size_t stride,
                                        const size_t n);

int gsl_fft_complex_radix2_dif_backward (gsl_complex_packed_array data,
                                         const size_t stride,
                                         const size_t n);

int gsl_fft_complex_radix2_dif_inverse (gsl_complex_packed_array data,
                                        const size_t stride,
                                        const size_t n);

int gsl_fft_complex_radix2_dif_transform (gsl_complex_packed_array data,
                                          const size_t stride,
                                          const size_t n,
                                          const gsl_fft_direction sign);

int gsl_fft_complex_bitreverse_order (gsl_complex_packed_array data,
                                      size_t stride,
                                      size_t n,
                                      size_t n_bits);

/*  Mixed Radix general-N routines  */

typedef struct
  {
    size_t n;
    size_t nf;
    size_t factor[64];
    gsl_complex *twiddle[64];
    gsl_complex *trig;
  }
gsl_fft_complex_wavetable;

typedef struct
{
  size_t n;
  double *scratch;
}
gsl_fft_complex_workspace;

gsl_fft_complex_wavetable *gsl_fft_complex_wavetable_alloc (size_t n);

void gsl_fft_complex_wavetable_free (gsl_fft_complex_wavetable * wavetable);

gsl_fft_complex_workspace *gsl_fft_complex_workspace_alloc (size_t n);

void gsl_fft_complex_workspace_free (gsl_fft_complex_workspace * workspace);

int gsl_fft_complex_memcpy (gsl_fft_complex_wavetable * dest,
                            gsl_fft_complex_wavetable * src);

int gsl_fft_complex_forward (gsl_complex_packed_array data,
                             const size_t stride,
                             const size_t n,
                             const gsl_fft_complex_wavetable * wavetable,
                             gsl_fft_complex_workspace * work);

int gsl_fft_complex_backward (gsl_complex_packed_array data,
                              const size_t stride,
                              const size_t n,
                              const gsl_fft_complex_wavetable * wavetable,
                              gsl_fft_complex_workspace * work);

int gsl_fft_complex_inverse (gsl_complex_packed_array data,
                             const size_t stride,
                             const size_t n,
                             const gsl_fft_complex_wavetable * wavetable,
                             gsl_fft_complex_workspace * work);

int gsl_fft_complex_transform (gsl_complex_packed_array data,
                               const size_t stride, const size_t n,
                               const gsl_fft_complex_wavetable * wavetable,
                               gsl_fft_complex_workspace * work,
                               const gsl_fft_direction sign);

__END_DECLS

#endif /* __GSL_FFT_COMPLEX_H__ */


text:
#ifndef _AVA_COMMON_SUPPORT_GEN_STAT_H_
#define _AVA_COMMON_SUPPORT_GEN_STAT_H_

#include <fmt/format.h>
#include <stdint.h>

#include <gsl/gsl>

#include "common/extensions/cmd_batching.h"
#include "common/support/io.h"
#include "guestlib/guest_thread.h"
#include "time_util.h"

namespace ava {
namespace support {

inline void stats_end(const char *func_name, int64_t begin_ts) {
  auto end_ts = ava::GetMonotonicNanoTimestamp();
  fmt::memory_buffer output;
  fmt::format_to(output, "GuestlibStat {}, {}\n", func_name, gsl::narrow_cast<int32_t>(end_ts - begin_ts));
  ava::guest_write_stats(output.data(), output.size());
}

}  // namespace support
}  // namespace ava

#endif  // _AVA_COMMON_SUPPORT_GEN_STAT_H_

max_stars_repo_name: photoszzt/ava

max_issues_repo_name: photoszzt/ava

max_forks_repo_name: photoszzt/ava