text:
/*******************************************/
/* libScalar:  -lm                         */
/*             -lIO                        */
/*                                         */
/* Scalar Arithmetic Operations and Macros */
/*                                         */
/* Nov 11 1999 Ju Li <liju99@mit.edu>      */
/*******************************************/

#ifndef _Scalar_h
#define _Scalar_h

#include <math.h>
#include <time.h>
#include <IO.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_spline.h>

/** Mathematical Constants **/

#define HALF      (0.5)
#define THIRD     (1./3.)
#define FOURTH    (0.25)
#define FIFTH     (0.2)
#define SIXTH     (1./6.)
#define SEVENTH   (1./7.)
#define EIGHTH    (0.125)
#define NINTH     (1./9.)
#define TENTH     (1./10.)
#define ELEVENTH  (1./11.)
#define TWELFTH   (1./12.)

/* extended precision to 60 digits: see Notes/mathconst.mws */
#define SQRT2 (1.41421356237309504880168872420969807856967187537694807317668)
#define SQRT3 (1.73205080756887729352744634150587236694280525381038062805581)
#define SQRT5 (2.23606797749978969640917366873127623544061835961152572427090)
#define SQRT6 (2.44948974278317809819728407470589139196594748065667012843269)
#define SQRT7 (2.64575131106459059050161575363926042571025918308245018036833)
#define SQRT8 (2.82842712474619009760337744841939615713934375075389614635336)
#define CBRT2 (1.25992104989487316476721060727822835057025146470150798008198)
#define CBRT3 (1.44224957030740838232163831078010958839186925349935057754642)
#define CBRT4 (1.58740105196819947475170563927230826039149332789985300980829)
#define CBRT5 (1.70997594667669698935310887254386010986805511054305492438286)
#define CBRT6 (1.81712059283213965889121175632726050242821046314121967148133)
#define CBRT7 (1.91293118277238910119911683954876028286243905034587576621065)
#define CBRT9 (2.08008382305190411453005682435788538633780534037326210969759)
/* (-1+sqrt(5.))/2 */
#define GOLDEN_RATIO \
  (.618033988749894848204586834365638117720309179805762862135450)
#define GOLDEN_RATIO_RECIPROCAL (GOLDEN_RATIO + 1)

#ifndef PI
#define PI    (3.14159265358979323846264338327950288419716939937510582097494)
#endif
#define RADIAN_TO_DEGREE(r) ((r)/PI*180.)
#define DEGREE_TO_RADIAN(d) ((d)/180.*PI)

#define SIN(x) sin(DEGREE_TO_RADIAN(x))
#define COS(x) cos(DEGREE_TO_RADIAN(x))
#define TAN(x) tan(DEGREE_TO_RADIAN(x))
#define ASIN(x) RADIAN_TO_DEGREE(asin(x))
#define ACOS(x) RADIAN_TO_DEGREE(acos(x))
#define ATAN(x) RADIAN_TO_DEGREE(atan(x))

#ifndef EPS
/* IEEE 754 floating-point arithmetic relative error bound */
#define EPS (2.220446e-16)
#endif

/* Deemed reasonably large number of floating point operations */
#ifndef REASONABLE_FLOPS
#define REASONABLE_FLOPS (1000)
#endif
/* on variables to estimate their likely numerical relative error */

/* When the condition number gets bad, for instance if subtracting */
/* two close numbers, we should override the above to larger value */

#ifndef TINY /* as compared to unity */
#define TINY (REASONABLE_FLOPS*EPS)
#endif

#ifndef HUGE /* as compared to unity */
#define HUGE (1./TINY)
#endif

#ifndef HUGE_INT
#define HUGE_INT (INT_MAX/2)
#endif

/* Close to the floating-point overflow */
/* IEEE 754 32-bit floating-point standard: 1.175494e-38 to 3.402823e38 */
#define SINGLE_PRECISION_INFINITY (1e37)
/* IEEE 754 64-bit floating-point standard: 2.225074e-308 to 1.797693e308 */
#define DOUBLE_PRECISION_INFINITY (1e307)

#ifndef SMALL /* as compared to unity */
#define SMALL (sqrt(EPS))
#endif

#ifndef LARGE /* as compared to unity */
 #define LARGE (1./SMALL)
#endif

/* a wise choice for estimating quadratic curvature */
#ifndef DELTA /* as compared to unity */
#define DELTA (sqrt(sqrt(EPS)))
#endif

/* very large exponent for later purposes of defining inf-norm, etc. */
#ifndef EXPONENT_INFINITY
#define EXPONENT_INFINITY (16.)
#endif

#define EQ(x,y)  ( (x) == (y) )
#define GE(x,y)  ( (x) >= (y) )
#define GT(x,y)  ( (x) >  (y) )
#define LE(x,y)  ( (x) <= (y) )
#define LT(x,y)  ( (x) <  (y) )
/* x = [lower_bound, upper_bound] */
#define XIN(x,lower_bound,upper_bound) (GE(x,lower_bound) && LE(x,upper_bound))
#define XINW(x,width) XIN(x,0.,width)
/* more likely to overflow than underflow, I think */
#define XOU(x,lower_bound,upper_bound) (GT(x,upper_bound) || LT(x,lower_bound))
#define XOUW(x,width) XOU(x,0.,width)
/* with no assumption about which is bigger, bound1 or bound2 */
#define INSIDE(x,bound1,bound2) (XIN(x,bound1,bound2) || XIN(x,bound2,bound1))
#define OUSIDE(x,bound1,bound2) (XOU(x,bound1,bound2) && XOU(x,bound2,bound1))
/* i = [lower_bound, upper_bound) */
#define IN(i,lower_bound,upper_bound) (GE(i,lower_bound) && LT(i,upper_bound))
#define INW(i,n) IN(i,0,n)  /* i = 0..n-1, C convention */
#define OU(i,lower_bound,upper_bound) (GE(i,upper_bound) || LT(i,lower_bound))
#define OUW(i,n) OU(i,0,n)

#define ABS(x)  ((x)>0?(x):-(x))
#define ABS_LE(x,positive_bound) XIN(x,-(positive_bound),(positive_bound))
#define ABS_GT(x,positive_bound) XOU(x,-(positive_bound),(positive_bound))
#define ABS_LT(x,positive_bound) \
  (GT(x,-(positive_bound)) && LT(x,(positive_bound)))
#define ABS_GE(x,positive_bound) \
  (GE(x,(positive_bound)) || LE(x,-(positive_bound)))

#ifndef MIN
#define MIN(x,y)      ((x)<(y)?(x):(y))
#endif
#ifndef MAX
#define MAX(x,y)      ((x)>(y)?(x):(y))
#endif
#define SWAP(x,y,tmp) ((tmp)=(x), (x)=(y), (y)=(tmp))
#define SWITCH(p,q,state0,state1) CAN ( if ((p) == (state0)) \
  { (p)=(state1); (q)=(state0); } else { (p)=(state0); (q)=(state1); } )
#define SWITCHSTATE(p,q,state) SWITCH(p,q,state,(state)+1)
#define ENSURE(min,max,tmp) { if (min>max) SWAP(min,max,tmp); }
#define FIND3(x0,x1,x2,op,xm) CAN( xm = x0; \
  if (xm op x1) xm = x1; if (xm op x2) xm = x2 )
#define FIND4(x0,x1,x2,x3,op,xm) CAN( xm = x0; \
  if (xm op x1) xm = x1; if (xm op x2) xm = x2; if (xm op x3) xm = x3 )
#define BEPOSITIVE(x)       CAN( if ((x) < 0) x = -(x) )
#define SQUARE(i)     ((i)*(i))
#define CUBE(i)       ((i)*(i)*(i))
#define QUAD(i)       ((i)*(i)*(i)*(i))
#define CEIL(i,j)     (((i)<=0)?0:((i)-1)/(j)+1)
#define SEIL(i,j)     ((size_t)CEIL(i,j))
#define ROUNDUP_TO(i,j)         (CEIL(i,j)*(j))
#define SEPARATION(x,y)         (fabs((x)-(y)))
#define SMALLSEPARATION(x,y,z)  (fabs((x)-(y))<z)

#define SPHERE_VOLUME(radius) (4.*PI/3.*CUBE(radius))
#define safe_avg(sum,count) ((double)(sum)/(((count)!=0)?(count):1))

#define NONEED_TRIM(x)  ( ((x) >= 0.) && ((x) <  1.) )
#define NEED_TRIM(x)    ( ((x) <  0.) || ((x) >= 1.) )
#define Trim(x) { if ((x)<0) do (x)++; while((x)<0); else \
  if ((x)>=1) do (x)--; while ((x)>=1); }
/* faster than TRIM(x) when x is around 1 */
#define TRIM(x) ((x)-floor((double)(x)))  /* x's image on [0,1) */
/* faster than Trim(x) when x is very large */

#define NONEED_IMAGE(x)  ( ((x) >= -0.5) && ((x) <  0.5) )
#define NEED_IMAGE(x)    ( ((x) <  -0.5) || ((x) >= 0.5) )
#define Image(x) { if ((x)<-0.5) do (x)++; while((x)<-0.5); else \
  if ((x)>=0.5) do (x)--; while ((x)>=0.5); }
/* faster than IMAGE(x) when x is around 1 */
#define IMAGE(x)  (TRIM(x+0.5)-0.5) /* x's image on [-0.5,0.5) */
/* faster than Image(x) when x is very large */

#define DISTANCE2(x,y,z) (SQUARE(x)+SQUARE(y)+SQUARE(z))
#define DISTANCE(x,y,z)   sqrt(DISTANCE2(x,y,z))
#define DISTANCED2(x,y)  (SQUARE(x)+SQUARE(y))
#define DISTANCE2D(x,y)   sqrt(DISTANCED2(x,y))

/* b^2 = c^2 - a^2:  the answer b will be nonnegative */
#define SafePythagorasComplement(c,a)                   \
    ( (c)>fabs(a) ? sqrt(SQUARE(c)-SQUARE(a)) : 0 )

/* x is floating-point */
#define Floor(x) (((x)==(int)(x))?(int)(x):(((x)>0)?((int)(x)):((int)(x)-1)))
#define FLOOR(x,i) { (i)=(int)(x); if ((i)>(x)) (i)--; }

#if defined(_alpha)
/* Compaq C warns about negative shifts which cannot actually happen */
#define SAFE_RSHIFT(a,n)  ( ((n)>=0) ? ((a)>>ABS(n)) : ((a)<<ABS(n)) )
#define SAFE_LSHIFT(a,n)  ( ((n)>=0) ? ((a)<<ABS(n)) : ((a)>>ABS(n)) )
#else
#define SAFE_RSHIFT(a,n)  ( ((n)>=0) ? ((a)>>(n)) : ((a)<<(-(n))) )
#define SAFE_LSHIFT(a,n)  ( ((n)>=0) ? ((a)<<(n)) : ((a)>>(-(n))) )
#endif

/* x, target can now be of any numeric type */
#define ERR(x,target) Ferr((double)(x), (double)(target))

#define ISTINY(a) (fabs(a)<TINY)
#define NOTTINY(a) (fabs(a)>=TINY)
#define ISSMALL(a) (fabs(a)<SMALL)
#define NOTSMALL(a) (fabs(a)>=SMALL)
#define ISDELTA(a) (fabs(a)<DELTA)
#define NOTDELTA(a) (fabs(a)>=DELTA)
#define NEAR(x,y)    ISTINY((x)-(y))
#define NOTNEAR(x,y) NOTTINY((x)-(y))

/* whether x and y are equal to floating point accuracy */
#define EQUAL(x,y) ((fabs((double)(x-y))<TINY) || (fabs(ERR(x,y))<TINY))

/* If a number is tiny, set it to zero. USE WITH CARE! */
#define SET_ZERO_IF_TINY(a) { if ISTINY(a) (a)=0.; }

/* Genesis of all randomness */
#define DEFAULT_RANDOM_SEED (1975425) /* for marking and reproducing results */

/* randomizing subsequent Frandom() sequence with seed */
#define Randomize(seed) srandom((unsigned int)(seed))
/* if never called, Frandom() automatically uses seed 1 */
#define TimeRandomize() Randomize(time((time_t *)NULL))

#ifdef _Linux
#define Frandom() ((double)(random()+0.5)/(1.+RAND_MAX))
#define FRANDOM() ((double)(random()+0.5)/(1.+RAND_MAX)-0.5)
#else
#define Frandom() ((double)(random()+0.5)/pow(2.,31.))
#define FRANDOM() ((double)(random()+0.5)/pow(2.,31.)-0.5)
#endif
/* Genesis of all randomness */

#define HALF_DECISION() (Frandom() > 0.5)

/* randomly pick an integer in min..max (set {min,min+1... max}) */
#define Fran(min,max) (INT(min)+(int)floor((INT(max)-INT(min)+1)*Frandom()))
/* randomly pick an integer in min..max-1 (set {min,min+1... max-1}) */
#define fran(min,max) (INT(min)+(int)floor((INT(max)-INT(min))*Frandom()))

/* return normally distributed random number with <x>=E, <(x-E)^2>=sigma2 */
#define Frandnorm(E,sigma2) \
(sqrt(-2.*log(Frandom()))*sin(2*acos(-1.)*Frandom())*sqrt((double)sigma2)+E)
/* The above is not most efficient. For array assignment, use Vfrandnorm() */

/* return normally distributed random number with <x>=0, <x^2>=1 */
#define FRANDNORM() (sqrt(-2.*log(Frandom()))*sin(2*acos(-1.)*Frandom()))

/* return exponentially distributed random number with <x>=tau */
#define Frandexp(tau) (-log(Frandom())*(tau))

/* return exponentially distributed random number with <x>=1 */
#define FRANDEXP() (-log(Frandom()))

/* if a is an integer factor of b */
#define isfactor(a,b) (((int)(b))%((int)(a))==0)
#define isnotfactor(a,b) (((int)(b))%((int)(a))!=0)

/* f(x) := x * Heaviside(x) */
#define ZERO_RAMP(x) (((x)>0)?(x):0)

/** Scalar.c: **/

/* greatest common divisor of a and b */
long gcd (long a, long b);
#define GCD(a,b)  gcd((long)(a),(long)(b))
/* least common multiple */
#define LCM(a,b) (long)(a)*(b)/gcd((long)(a),(long)(b))
/* check if a is a factor of b */
#define ISFACTOR(a,b) (((long)(b))%((long)(a))==0)
#define ISEVEN(a) (((unsigned int)(a)&1)==0)
#define ISODD(a)  (((unsigned int)(a)&1))

/* find k,l such that k * x + l * y = gcd(x,y) */
void diophantine2 (int x, int y, int *k, int *l);

/* set Bitmap "b", starting from the bit "offset", every "period" */
/* bits to the lowest "period" bits of "value", for "n" periods.  */
Bmap *Bperiod (Bmap *b, int offset, int n, int period, int value);

/* relative error of x compared to target */
#define rerr(x,target) ((x)/(target)-1)
double Ferr(double x, double target);

/* calculates |_ a/b _| with mathematically correct floor */
int floorDiv (int a, int b);

/* calculates |^ a/b ^| with mathamatically correct ceiling */
int ceilDiv (int a,int b);

/* return (n!) = GAMMA(n+1) in double precision */
double factorial (int n);

/* Take average on structures whose first element is the counter  */
#define avg_what(type,s)  (sizeof(type)/sizeof(double)),((double *)(s))
#define Avg_what(type,s)  (sizeof(type)/sizeof(double)),((double *)(&(s)))
void avg_clear (int n, double *s);
int avg_add (int n, double *s, ...);
int avg_done (int n, double *s);
int avg_recover (int n, double *s);
/* to avoid alignment trouble please define the counter as double */

/* Generate a string of n-1 random base64 characters ended by EOS. */
char *RandomBase64String (int n, char *a);

/* (-1)%3=-1, but positive_remainder(-1,3)=2 */
int positive_remainder(int a, int b);

/* quadrature.c: */

/* Integration of one-dimensional function: */

/* Gaussian quadrature with W(x)=1, N=10, NR pp.148 */
double qgaus (double (*func)(double), double a, double b);

/* Romberg quadrature: */
#define QROMB_JMAX               20
#define QROMB_K                  5
#define QROMB_DEFAULT_TOLERANCE  1e-14
extern double qromb_tolerance;
double qromb (double (*func)(double), double a, double b);

/* newtonian.c: */

/************************************************************/
/* K-evaluation Integration of Newtonian Dynamics:          */
/*                                                          */
/* Given x(0),dotx(0) and routine ddotx(x) which is derived */
/* from a potential, we want to get x(h),dotx(h) where h is */
/* called a full step, during which K calls to ddotx(x) are */
/* made. Let g = h / K. For a scheme to be competitive, it  */
/* should give better results at t_f = Lh >> h compared to  */
/* other schemes with equal g, for some g/accuracy choices. */
/*                                                          */
/* Algorithms with so-called FSAL property can combine two  */
/* consecutive h-steps to save one evaluation. We provide   */
/* _start, _fsal, _stop routines where _start and _fsal     */
/* need K evaluations and _stop may need one. For complete- */
/* ness, we also give _nofsal which needs K+1 evaluations.  */
/************************************************************/

/* k=0..K-1 (some _nofsal/_stop call k=K), n, x[n], acc[n]=workspace */
typedef void (*Newtonian_ddotx) (int, int, double *, double *);

/***************************/
/* Symplectic Algorithm 1: */
/***************************/
/* First x, then dotx evolution. No great need to use FSAL savings. */

/* This algorithm is not the best for heat current computation or */
/* total energy checksum because when ddotx() is called at k=0, x */
/* is already shifted. Ruth83, Schlier98_6a, Tsitouras99.         */

/****************************************************/
/* Computer "Experiments" on Classical Fluids. I.   */
/* Thermodynamical Properties of Lennard-Jones      */
/* Molecules, L. Verlet, Phys. Rev. 159, 98 (1967). */
/* A Canonical Integration Technique, R.D. Ruth,    */
/* IEEE Trans. Nucl. Sci. NS-30, 2669 (1983).       */
/****************************************************/
/* evaluations per h-step:  1   */
/* global truncation error: h^2 */

/* First x, then dotx evolution. No great need to use FSAL savings. */
void symplectic_Ruth83_nofsal
(double h, int n, double *x, double *dotx, double *acc,
 Newtonian_ddotx ddotx);

/* FSAL sequence: _start, _fsal, _stop */
void symplectic_Ruth83_start
(double h, int n, double *x, double *dotx, double *acc,
 Newtonian_ddotx ddotx);
void symplectic_Ruth83_fsal
(double h, int n, double *x, double *dotx, double *acc,
 Newtonian_ddotx ddotx);
void symplectic_Ruth83_stop
(double h, int n, double *x, double *dotx, double *acc,
 Newtonian_ddotx ddotx);

/**************************************************************/
/* Ch. Schlier and A. Seiter, J. Phys. Chem. 102 (1998) 9399. */
/**************************************************************/
/* evaluations per h-step:  9   */
/* global truncation error: h^6 */

/* First x, then dotx evolution. No great need to use FSAL savings. */
void symplectic_Schlier98_6a_nofsal
(double h, int n, double *x, double *dotx, double *acc,
 Newtonian_ddotx ddotx);

/* FSAL sequence: _start, _fsal, _stop */
void symplectic_Schlier98_6a_start
(double h, int n, double *x, double *dotx, double *acc,
 Newtonian_ddotx ddotx);
void symplectic_Schlier98_6a_fsal
(double h, int n, double *x, double *dotx, double *acc,
 Newtonian_ddotx ddotx);
void symplectic_Schlier98_6a_stop
(double h, int n, double *x, double *dotx, double *acc,
 Newtonian_ddotx ddotx);

/************************************************************/
/* Haruo Yoshida, Phys. Lett. A 150, 262 (1990).            */
/* Celestial Mechanics and Dynamical Astronomy 74 (1999)    */
/* 223-230, Ch. Tsitouras,                                  */
/* www.math.ntua.gr/people/tsitoura/publications.html       */
/************************************************************/
/* evaluations per h-step:  33   */
/* global truncation error: h^10 */

/* First x, then dotx evolution. No great need to use FSAL savings. */
void symplectic_Tsitouras99_nofsal
(double h, int n, double *x, double *dotx, double *acc,
 Newtonian_ddotx ddotx);

/* FSAL sequence: _start, _fsal, _stop */
void symplectic_Tsitouras99_start
(double h, int n, double *x, double *dotx, double *acc,
 Newtonian_ddotx ddotx);
void symplectic_Tsitouras99_fsal
(double h, int n, double *x, double *dotx, double *acc,
 Newtonian_ddotx ddotx);
void symplectic_Tsitouras99_stop
(double h, int n, double *x, double *dotx, double *acc,
 Newtonian_ddotx ddotx);

/***************************/
/* Symplectic Algorithm 2: */
/***************************/
/* First dotx, then x evolution. You do need FSAL savings. */

/* This algorithm is more appropriate for heat current computation */
/* or total energy checksum when ddotx() of _start is called at k  */
/* = 0, if _stop was called at the last step, because exact p(nh), */
/* q(nh) would then be passed to ddotx(). Calvo93, Schlier00_6b,   */
/* Schlier00_8b, Schlier00_8c. */

/************************************************************/
/* The Development of Variable-Step Symplectic Integrators, */
/* With Application to the Two-Body Problem, M.P. Calvo,    */
/* J.M. Sanz-Serna, SIAM J. Sci. Comput. 14, 936-952 (1993) */
/* Numerical Hamiltonian Problems, J.M. Sanz-Serna,         */
/* M.P. Calvo, Chapman & Hall, London (1994).               */
/************************************************************/
/* evaluations per h-step:  4   */
/* global truncation error: h^4 */

/* First dotx, then x evolution. No FSAL savings yet. */
void symplectic_Calvo93_nofsal
(double h, int n, double *x, double *dotx, double *acc,
 Newtonian_ddotx ddotx);

/* FSAL sequence: _start, _fsal, _stop */
void symplectic_Calvo93_start
(double h, int n, double *x, double *dotx, double *acc,
 Newtonian_ddotx ddotx);
void symplectic_Calvo93_fsal
(double h, int n, double *x, double *dotx, double *acc,
 Newtonian_ddotx ddotx);
void symplectic_Calvo93_stop
(double h, int n, double *x, double *dotx, double *acc,
 Newtonian_ddotx ddotx);

/*******************************************************/
/* High-Order Symplectic Integration: An Assessment,   */
/* Ch. Schlier, A. Seiter, Comp. Phys. Comm. (2000)    */
/* phya8.physik.uni-freiburg.de/abt/papers/papers.html */
/*******************************************************/
/* evaluations per h-step:  8   */
/* global truncation error: h^6 */

/* First dotx, then x evolution. No FSAL savings yet. */
void symplectic_Schlier00_6b_nofsal
(double h, int n, double *x, double *dotx, double *acc,
 Newtonian_ddotx ddotx);

/* FSAL sequence: _start, _fsal, _stop */
void symplectic_Schlier00_6b_start
(double h, int n, double *x, double *dotx, double *acc,
 Newtonian_ddotx ddotx);
void symplectic_Schlier00_6b_fsal
(double h, int n, double *x, double *dotx, double *acc,
 Newtonian_ddotx ddotx);
void symplectic_Schlier00_6b_stop
(double h, int n, double *x, double *dotx, double *acc,
 Newtonian_ddotx ddotx);

/*******************************************************/
/* Ch. Schlier and A. Seiter, Comp. Phys. Comm. (2000) */
/*******************************************************/
/* evaluations per h-step:  17  */
/* global truncation error: h^8 */

/* First dotx, then x evolution. No FSAL savings yet.  */
void symplectic_Schlier00_8b_nofsal
(double h, int n, double *x, double *dotx, double *acc,
 Newtonian_ddotx ddotx);

/* FSAL sequence: _start, _fsal, _stop */
void symplectic_Schlier00_8b_start
(double h, int n, double *x, double *dotx, double *acc,
 Newtonian_ddotx ddotx);
void symplectic_Schlier00_8b_fsal
(double h, int n, double *x, double *dotx, double *acc,
 Newtonian_ddotx ddotx);
void symplectic_Schlier00_8b_stop
(double h, int n, double *x, double *dotx, double *acc,
 Newtonian_ddotx ddotx);

/*******************************************************/
/* Ch. Schlier and A. Seiter, Comp. Phys. Comm. (2000) */
/*******************************************************/
/* evaluations per h-step:  17  */
/* global truncation error: h^8 */

/* First dotx, then x evolution. No FSAL savings yet.  */
void symplectic_Schlier00_8c_nofsal
(double h, int n, double *x, double *dotx, double *acc,
 Newtonian_ddotx ddotx);

/* FSAL sequence: _start, _fsal, _stop */
void symplectic_Schlier00_8c_start
(double h, int n, double *x, double *dotx, double *acc,
 Newtonian_ddotx ddotx);
void symplectic_Schlier00_8c_fsal
(double h, int n, double *x, double *dotx, double *acc,
 Newtonian_ddotx ddotx);
void symplectic_Schlier00_8c_stop
(double h, int n, double *x, double *dotx, double *acc,
 Newtonian_ddotx ddotx);

/**********************************************************************/
/* Routers to various routines by K, the number of evaluations per h. */
/* To make fair comparison, one should use equal g for different K's. */
/* It gives "exact" x(wh) and dotx(wh), w=1..L, after each full step. */
/**********************************************************************/
void symplectic_nofsal ( int K, double g, int n, double *x, double *dotx,
                         double *acc, Newtonian_ddotx ddotx );

/**********************************************************************/
/* Routers to various routines by K, the number of evaluations per h. */
/* To make fair comparison, one should use equal g for different K's. */
/* To integrate t from 0 to Lh, one should call 1 x symplectic_start, */
/* (L-1) x symplectic_fsal, and 1 x symplectic_stop. One cannot have  */
/* "exact" x(wh) and dotx(wh) except when w=0,L, but capturing k=0    */
/* calls to ddotx(k,n,x,acc) at w=1..L-1 gives approximate answers.   */
/**********************************************************************/
void symplectic_start ( int K, double g, int n, double *x, double *dotx,
                        double *acc, Newtonian_ddotx ddotx );
void symplectic_fsal ( int K, double g, int n, double *x, double *dotx,
                       double *acc, Newtonian_ddotx ddotx );
void symplectic_stop ( int K, double g, int n, double *x, double *dotx,
                       double *acc, Newtonian_ddotx ddotx );

/*******************************************************/
/* Well known non-symplectic integrators & Benchmarks: */
/*******************************************************/

/************************************************************/
/* classic 4th-order Runge-Kutta method: Numerical Recipes  */
/* in C: The Art of Scientific Computing, William H. Press, */
/* Saul A. Teukolsky, William T. Vetterling, Brian P.       */
/* Flannery, Cambridge University Press, Cambridge (1992).  */
/* http://mathworld.wolfram.com/Runge-KuttaMethod.html      */
/************************************************************/
/* evaluations per h-step:  4   */
/* global truncation error: h^4 */

/* suited for heat current computation & total energy checksum at k=0 */

/* allocate 8*n internal workspace. No integration is performed yet. */
void rk4_prep ( double h, int n, double *x, double *dotx, double *acc,
                Newtonian_ddotx ddotx );
/* "acc" would not be used, so it can be passed as NULL. */
void rk4 (double h, int n, double *x, double *dotx, double *acc,
          Newtonian_ddotx ddotx);
/* free 8*n internal workspace */
void rk4_free ( double h, int n, double *x, double *dotx, double *acc,
                Newtonian_ddotx ddotx );

/**********************************************************/
/* Gear predictor-corrector algorithm: Numerical Initial  */
/* Value Problems in Ordinary Differential Equation,      */
/* C.W. Gear, Englewood Cliffs N.J., Prentice-Hall, 1971. */
/* itp.nat.uni-magdeburg.de/~schinner/dip/node22.html     */
/**********************************************************/
/* k-value: keep k n-arrays: x^0,x^1,..,x^{k-1}:          */
/* initialized by "exact" time derivatives.               */
/* Before corrector, x0 has local truncation error h^k;   */
/* after corrector, x0 has local truncation error h^(k+1) */
/* therefore the global truncation error of x0 is h^k.    */
/* The global truncation error of dotx is h^{k-1}, etc.   */
/**********************************************************/

#define GEAR_TEST_H  50

/* 4-value Gear predictor-corrector algorithm */
/* evaluations per h-step:  1   */
/* local  truncation error: h^5 */
/* global truncation error: h^4 */

/* Initialize 4-value Gear integrator and allocate 2*n  */
/* internal workspace. No integration is performed yet. */
void gear4_prep ( double h, int n, double *x, double *dotx, double *acc,
                  Newtonian_ddotx ddotx );
void gear4 ( double h, int n, double *x, double *dotx, double *acc,
             Newtonian_ddotx ddotx );
/* free 2*n internal workspace */
void gear4_free ( double h, int n, double *x, double *dotx, double *acc,
                  Newtonian_ddotx ddotx );

/******************************************************/
/* 5-value Gear predictor-corrector algorithm:        */
/* Computer Simulation of Liquids, M.P. Allen and     */
/* D.J. Tildesley, Clarendon Press, (Oxford 1987).    */
/* ftp://ftp.dl.ac.uk/ccp5/ALLEN_TILDESLEY/F.02       */
/* itp.nat.uni-magdeburg.de/~schinner/dip/node22.html */
/******************************************************/
/* evaluations per h-step:  1   */
/* local  truncation error: h^6 */
/* global truncation error: h^5 */

/* Initialize 5-value Gear integrator and allocate 3*n  */
/* internal workspace. No integration is performed yet. */
void gear5_prep ( double h, int n, double *x, double *dotx, double *acc,
                  Newtonian_ddotx ddotx );
void gear5 ( double h, int n, double *x, double *dotx, double *acc,
             Newtonian_ddotx ddotx );
/* free 3*n internal workspace */
void gear5_free ( double h, int n, double *x, double *dotx, double *acc,
                  Newtonian_ddotx ddotx );

/* 6-value Gear predictor-corrector algorithm */
/* evaluations per h-step:  1   */
/* local  truncation error: h^7 */
/* global truncation error: h^6 */

/* Initialize 6-value Gear integrator and allocate 4*n  */
/* internal workspace. No integration is performed yet. */
void gear6_prep ( double h, int n, double *x, double *dotx, double *acc,
                  Newtonian_ddotx ddotx );
void gear6 ( double h, int n, double *x, double *dotx, double *acc,
             Newtonian_ddotx ddotx );
/* free 4*n internal workspace */
void gear6_free ( double h, int n, double *x, double *dotx, double *acc,
                  Newtonian_ddotx ddotx );

/* 7-value Gear predictor-corrector algorithm */
/* evaluations per h-step:  1   */
/* local  truncation error: h^8 */
/* global truncation error: h^7 */

/* Initialize 7-value Gear integrator and allocate 5*n  */
/* internal workspace. No integration is performed yet. */
void gear7_prep ( double h, int n, double *x, double *dotx, double *acc,
                  Newtonian_ddotx ddotx );
void gear7 ( double h, int n, double *x, double *dotx, double *acc,
             Newtonian_ddotx ddotx );
/* free 5*n internal workspace */
void gear7_free ( double h, int n, double *x, double *dotx, double *acc,
                  Newtonian_ddotx ddotx );

/* 8-value Gear predictor-corrector algorithm */
/* evaluations per h-step:  1   */
/* local  truncation error: h^9 */
/* global truncation error: h^8 */

/* Initialize 8-value Gear integrator and allocate 6*n  */
/* internal workspace. No integration is performed yet. */
void gear8_prep ( double h, int n, double *x, double *dotx, double *acc,
                  Newtonian_ddotx ddotx );
void gear8 ( double h, int n, double *x, double *dotx, double *acc,
             Newtonian_ddotx ddotx );
/* free 6*n internal workspace */
void gear8_free ( double h, int n, double *x, double *dotx, double *acc,
                  Newtonian_ddotx ddotx );

/* symplectic algorithm 1 */
#define NEWTONIAN_RUTH83          0
#define NEWTONIAN_SCHLIER98_6A   (NEWTONIAN_RUTH83 + 1)
#define NEWTONIAN_TSITOURAS99    (NEWTONIAN_SCHLIER98_6A + 1)
#define NEWTONIAN_SYMPLECTIC_ALGORITHM_1_MAX  (NEWTONIAN_TSITOURAS99 + 1)
/* symplectic algorithm 2 */
#define NEWTONIAN_CALVO93        NEWTONIAN_SYMPLECTIC_ALGORITHM_1_MAX
#define NEWTONIAN_SCHLIER00_6B   (NEWTONIAN_CALVO93 + 1)
#define NEWTONIAN_SCHLIER00_8B   (NEWTONIAN_SCHLIER00_6B + 1)
#define NEWTONIAN_SCHLIER00_8C   (NEWTONIAN_SCHLIER00_8B + 1)
#define NEWTONIAN_SYMPLECTIC_ALGORITHM_2_MAX  (NEWTONIAN_SCHLIER00_8C + 1)
/* non-symplectic algorithms */
#define NEWTONIAN_RK4            NEWTONIAN_SYMPLECTIC_ALGORITHM_2_MAX
#define NEWTONIAN_GEAR4          (NEWTONIAN_RK4 + 1)
#define NEWTONIAN_GEAR5          (NEWTONIAN_GEAR4 + 1)
#define NEWTONIAN_GEAR6          (NEWTONIAN_GEAR5 + 1)
#define NEWTONIAN_GEAR7          (NEWTONIAN_GEAR6 + 1)
#define NEWTONIAN_GEAR8          (NEWTONIAN_GEAR7 + 1)
#define NEWTONIAN_MAX            (NEWTONIAN_GEAR8 + 1)

typedef void (*Integrator) ( double, int, double *, double *, double *,
                             Newtonian_ddotx );
typedef struct
{
    char *name;
    bool is_symplectic;
    int K;
    int order;
    Integrator a,b,c;
} Newtonian;

extern Newtonian newtonian [NEWTONIAN_MAX];

/***************************************************************/
/* prep  +  [ firststep, L-1 nextstep, halt ] repeat  +  free  */
/*            |                                                */
/* k=0 property calculation for algorithm 2 and non-symplectic */
/***************************************************************/

/* Allocate and initialize states if necessary, but do NOT integrate. */
void newtonian_integrator_prep
( int id, double g, int n, double *x, double *dotx,
  double *acc, Newtonian_ddotx ddotx );

/* Get the integrator to run its first full step of a session. */
void newtonian_integrator_firststep
( int id, double g, int n, double *x, double *dotx,
  double *acc, Newtonian_ddotx ddotx );

/* Run other full steps of a session. */
void newtonian_integrator_nextstep
( int id, double g, int n, double *x, double *dotx,
  double *acc, Newtonian_ddotx ddotx );

/* Halt the session and do necessary mop-ups. To */
/* start a new session, call _firststep again.   */
void newtonian_integrator_halt
( int id, double g, int n, double *x, double *dotx,
  double *acc, Newtonian_ddotx ddotx );

/* Free workspace if any. No more sessions. To restart, call _prep. */
void newtonian_integrator_free
( int id, double g, int n, double *x, double *dotx,
  double *acc, Newtonian_ddotx ddotx );

/** spline.c: **/

#define SPLINE_CUBIC_NATURAL   0
#define SPLINE_CUBIC_PERIODIC  1
#define SPLINE_AKIMA_NATURAL   2
#define SPLINE_AKIMA_PERIODIC  3

typedef struct
{
    gsl_spline *spline;
    gsl_interp_accel *acc;
    double xmin, xmax, ymax, ymin;
} Spline;

/* Allocate and initialize Spline structure based on */
/* lookup table x[0..N-1],y[0..N-1] and spline type. */
Spline *spline_initialize (int N, double *x, double *y, int spline_type);

/* Evaluate spline value. */
double spline_value (Spline *sp, double x);

/* Evaluate first-order derivative. */
double spline_deriv (Spline *sp, double x);

/* Evaluate second-order derivative. */
double spline_deriv2 (Spline *sp, double x);

/* Free all data structure. */
void spline_finalize (Spline *sp);

/** histogram.c: **/

typedef struct
{
    double xmin, xmax;   /* domain of interest: {xmin,xmax} */
    int bins;            /* # of collection bins in the domain of interest */
    double xbin;         /* width of one bin */
    int totalcount;      /* counts of total submission */
    double domaincount;  /* counts that fall into the domain of interest */
    double *bincount;    /* count in each bin */
    double *x_sample;    /* center x of each bin */
    double *density_profile;  /* density profile which integrates to 1 */
} Histogram;

/* Given finite xmin < xmax and bins >= 1, reallocate */
/* data structure in h and reset all counters.        */
void histogram_reinitialize
(double xmin, double xmax, int bins, Histogram *h);

/* Update histogram by data array x[0..n-1]; return domaincount */
double histogram_intercept (int n, double *x, Histogram *h);

/* Return the interception rate = domaincount / totalcount */
double histogram_interception_rate (Histogram *h);

/* Update the density profile; return the interception rate */
double histogram_update_density_profile (Histogram *h);

/* Save histogram in matlab script "fname" as  */
/* symbol "token" and notify I/O stream "info" */
double histogram_save_as_matlab
(Histogram *h, FILE *info, char *token, char *fname);
#define histogram_save_as_MATLAB(h,info,token) histogram_save_as_matlab\
  (h,info,token,str2(token,"_h.m"))

/* Clear all counters as if it is just the beginning */
void histogram_reset_all_counters (Histogram *h);

/* Free the histogram data structure and reset all counters */
void histogram_free (Histogram *h);

#endif  /* _Scalar_h */

max_stars_repo_name: jameskermode/AtomEye

max_issues_repo_name: jameskermode/AtomEye

max_forks_repo_name: jameskermode/AtomEye


text:
/* linalg/test_lu.c
 *
 * Copyright (C) 2019 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 <gsl/gsl_test.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_ieee_utils.h>
#include <gsl/gsl_permute_vector.h>
#include <gsl/gsl_blas.h>
#include <gsl/gsl_linalg.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_permutation.h>

static int
test_LU_decomp_eps(const gsl_matrix * m, const double eps, const char * desc)
{
  int s = 0;
  const size_t M = m->size1;
  const size_t N = m->size2;

  gsl_matrix * A = gsl_matrix_alloc(M, N);
  gsl_matrix * PLU = gsl_matrix_calloc(M, N);
  gsl_permutation * p = gsl_permutation_alloc(M);
  size_t i, j;
  int signum;

  gsl_matrix_memcpy(A, m);
  gsl_linalg_LU_decomp(A, p, &signum);

  if (M >= N)
    {
      gsl_matrix_view U = gsl_matrix_submatrix(A, 0, 0, N, N);

      /* copy L factor into PLU */
      for (j = 0; j < GSL_MIN(M, N); ++j)
        {
          gsl_vector_view v1 = gsl_matrix_subcolumn(A, j, j, M - j);
          gsl_vector_view v2 = gsl_matrix_subcolumn(PLU, j, j, M - j);

          gsl_vector_memcpy(&v2.vector, &v1.vector);
          gsl_matrix_set(PLU, j, j, 1.0);
        }

      /* PLU := L * U */
      gsl_blas_dtrmm(CblasRight, CblasUpper, CblasNoTrans, CblasNonUnit, 1.0, &U.matrix, PLU);
    }
  else
    {
      gsl_matrix_view L = gsl_matrix_submatrix(A, 0, 0, M, M);

      /* copy U factor into PLU */
      for (i = 0; i < GSL_MIN(M, N); ++i)
        {
          gsl_vector_view v1 = gsl_matrix_subrow(A, i, i, N - i);
          gsl_vector_view v2 = gsl_matrix_subrow(PLU, i, i, N - i);
          gsl_vector_memcpy(&v2.vector, &v1.vector);
        }

      /* PLU := L * U */
      gsl_blas_dtrmm(CblasLeft, CblasLower, CblasNoTrans, CblasUnit, 1.0, &L.matrix, PLU);
    }

  /* PLU := P * L * U */
  for (j = 0; j < N; ++j)
    {
      gsl_vector_view v = gsl_matrix_column(PLU, j);
      gsl_permute_vector_inverse(p, &v.vector);
    }

  /* now test m = PLU */
  for (i = 0; i < M; ++i)
    {
      for (j = 0; j < N; ++j)
        {
          double aij = gsl_matrix_get(PLU, i, j);
          double mij = gsl_matrix_get(m, i, j);

          gsl_test_rel(aij, mij, eps,
                       "%s: (%3lu,%3lu)[%lu,%lu]: %22.18g   %22.18g\n",
                       desc, M, N, i, j, aij, mij);
        }
    }

  gsl_matrix_free(A);
  gsl_matrix_free(PLU);
  gsl_permutation_free(p);

  return s;
}

static int
test_LU_decomp(gsl_rng * r)
{
  int s = 0;
  size_t n;

  /* test square matrices */
  for (n = 1; n <= 50; ++n)
    {
      gsl_matrix * m = gsl_matrix_alloc(n, n);

      create_random_matrix(m, r);
      test_LU_decomp_eps(m, 4096.0 * n * GSL_DBL_EPSILON, "LU_decomp random");

      create_hilbert_matrix2(m);
      test_LU_decomp_eps(m, 256.0 * n * GSL_DBL_EPSILON, "LU_decomp hilbert");

      gsl_matrix_free(m);
    }

  {
    gsl_matrix * m = gsl_matrix_alloc(100, 50);
    create_random_matrix(m, r);
    test_LU_decomp_eps(m, 256.0 * n * GSL_DBL_EPSILON, "LU_decomp rect1");
    gsl_matrix_free(m);
  }

  {
    gsl_matrix * m = gsl_matrix_alloc(50, 100);
    create_random_matrix(m, r);
    test_LU_decomp_eps(m, 256.0 * n * GSL_DBL_EPSILON, "LU_decomp rect2");
    gsl_matrix_free(m);
  }

  {
    gsl_matrix * m = gsl_matrix_alloc(80, 100);
    create_random_matrix(m, r);
    test_LU_decomp_eps(m, 256.0 * n * GSL_DBL_EPSILON, "LU_decomp rect3");
    gsl_matrix_free(m);
  }

  return s;
}

static int
test_LU_solve_eps(const gsl_matrix * m, const gsl_vector * rhs, const gsl_vector * sol, const double eps, const char * desc)
{
  int s = 0;
  const size_t N = m->size1;
  int signum;
  size_t i;

  gsl_permutation * perm = gsl_permutation_alloc(N);
  gsl_matrix * lu  = gsl_matrix_alloc(N, N);
  gsl_vector * x = gsl_vector_alloc(N);
  gsl_vector * residual = gsl_vector_alloc(N);

  gsl_matrix_memcpy(lu, m);

  s += gsl_linalg_LU_decomp(lu, perm, &signum);
  s += gsl_linalg_LU_solve(lu, perm, rhs, x);

  for (i = 0; i < N; i++)
    {
      double xi = gsl_vector_get(x, i);
      double yi = gsl_vector_get(sol, i);

      gsl_test_rel(xi, yi, eps, "%s: %3lu[%lu]: %22.18g    %22.18g\n",
                   desc, N, i, xi, yi);
    }

  s += gsl_linalg_LU_refine(m, lu, perm, rhs, x, residual);

  for (i = 0; i < N; i++)
    {
      double xi = gsl_vector_get(x, i);
      double yi = gsl_vector_get(sol, i);

      gsl_test_rel(xi, yi, eps, "%s: improved %3lu[%lu]: %22.18g    %22.18g\n",
                   desc, N, i, xi, yi);
    }

  gsl_vector_free(residual);
  gsl_vector_free(x);
  gsl_matrix_free(lu);
  gsl_permutation_free(perm);

  return s;
}

static int
test_LU_solve(gsl_rng * r)
{
  int s = 0;
  size_t n;

  for (n = 1; n <= 50; ++n)
    {
      gsl_matrix * m = gsl_matrix_alloc(n, n);
      gsl_vector * rhs = gsl_vector_alloc(n);
      gsl_vector * sol = gsl_vector_alloc(n);

      create_random_matrix(m, r);
      create_random_vector(sol, r);
      gsl_blas_dgemv(CblasNoTrans, 1.0, m, sol, 0.0, rhs);
      test_LU_solve_eps(m, rhs, sol, 1.0e5 * n * GSL_DBL_EPSILON, "LU_solve random");

      if (n <= 4)
        {
          create_hilbert_matrix2(m);
          gsl_blas_dsymv(CblasLower, 1.0, m, sol, 0.0, rhs);
          test_LU_solve_eps(m, rhs, sol, 4096.0 * n * GSL_DBL_EPSILON, "LU_solve hilbert");
        }

      gsl_matrix_free(m);
      gsl_vector_free(rhs);
      gsl_vector_free(sol);
    }

  return s;
}

int
test_LU_invert_eps(const gsl_matrix * m, const double eps, const char *desc)
{
  int s = 0;
  size_t i, j, N = m->size1;

  gsl_matrix * lu  = gsl_matrix_alloc(N, N);
  gsl_matrix * c  = gsl_matrix_alloc(N, N);
  gsl_permutation * perm = gsl_permutation_alloc(N);
  int signum;

  gsl_matrix_memcpy(lu, m);

  s += gsl_linalg_LU_decomp(lu, perm, &signum);
  s += gsl_linalg_LU_invx(lu, perm);

  /* c = m m^{-1} */
  gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1.0, m, lu, 0.0, c);

  /* c should be the identity matrix */

  for (i = 0; i < N; ++i)
    {
      for (j = 0; j < N; ++j)
        {
          double cij = gsl_matrix_get(c, i, j);
          double expected = (i == j) ? 1.0 : 0.0;

          gsl_test_rel(cij, expected, eps, "%s (%3lu,%3lu)[%lu,%lu]: %22.18g   %22.18g\n",
                       desc, N, N, i, j, cij, expected);
        }
    }

  gsl_matrix_free(lu);
  gsl_matrix_free(c);
  gsl_permutation_free(perm);

  return s;
}

static int
test_LU_invert(gsl_rng * r)
{
  int s = 0;
  size_t n;

  for (n = 1; n <= 50; ++n)
    {
      gsl_matrix * m = gsl_matrix_alloc(n, n);

      create_random_matrix(m, r);
      test_LU_invert_eps(m, 1.0e5 * n * GSL_DBL_EPSILON, "LU_invert random");

      gsl_matrix_free(m);
    }

  return s;
}

ext: c

max_stars_repo_name: ielomariala/Hex-Game

max_issues_repo_name: ielomariala/Hex-Game

max_forks_repo_name: karanbirsandhu/nu-sense


text:
/**
 *
 * @file testing_zposv.c
 *
 *  PLASMA testing routines
 *  PLASMA is a software package provided by Univ. of Tennessee,
 *  Univ. of California Berkeley and Univ. of Colorado Denver
 *
 * @version 2.6.0
 * @author Bilel Hadri, Hatem Ltaief
 * @date 2010-11-15
 * @precisions normal z -> c d s
 *
 **/
#include <stdlib.h>
#include <stdio.h>
#include <string.h>
#include <math.h>

#include <plasma.h>
#include <cblas.h>
#include <lapacke.h>
#include <core_blas.h>
#include "testing_zmain.h"

enum blas_order_type {
            blas_rowmajor = 101,
            blas_colmajor = 102 };

enum blas_cmach_type {
            blas_base      = 151,
            blas_t         = 152,
            blas_rnd       = 153,
            blas_ieee      = 154,
            blas_emin      = 155,
            blas_emax      = 156,
            blas_eps       = 157,
            blas_prec      = 158,
            blas_underflow = 159,
            blas_overflow  = 160,
            blas_sfmin     = 161};

enum blas_norm_type {
            blas_one_norm       = 171,
            blas_real_one_norm  = 172,
            blas_two_norm       = 173,
            blas_frobenius_norm = 174,
            blas_inf_norm       = 175,
            blas_real_inf_norm  = 176,
            blas_max_norm       = 177,
            blas_real_max_norm  = 178 };

static void
BLAS_error(char *rname, int err, int val, int x) {
  fprintf( stderr, "%s %d %d %d\n", rname, err, val, x );
  abort();
}

static
void
BLAS_zge_norm(enum blas_order_type order, enum blas_norm_type norm,
  int m, int n, const PLASMA_Complex64_t *a, int lda, double *res) {
  int i, j; float anorm, v;
  char rname[] = "BLAS_zge_norm";

  if (order != blas_colmajor) BLAS_error( rname, -1, order, 0 );

  if (norm == blas_frobenius_norm) {
    anorm = 0.0f;
    for (j = n; j; --j) {
      for (i = m; i; --i) {
        v = a[0];
        anorm += v * v;
        a++;
      }
      a += lda - m;
    }
    anorm = sqrt( anorm );
  } else if (norm == blas_inf_norm) {
    anorm = 0.0f;
    for (i = 0; i < m; ++i) {
      v = 0.0f;
      for (j = 0; j < n; ++j) {
        v += cabs( a[i + j * lda] );
      }
      if (v > anorm)
        anorm = v;
    }
  } else {
    BLAS_error( rname, -2, norm, 0 );
    return;
  }

  if (res) *res = anorm;
}

static
double
BLAS_dpow_di(double x, int n) {
  double rv = 1.0;

  if (n < 0) {
    n = -n;
    x = 1.0 / x;
  }

  for (; n; n >>= 1, x *= x) {
    if (n & 1)
      rv *= x;
  }

  return rv;
}

static
double
BLAS_dfpinfo(enum blas_cmach_type cmach) {
  double eps = 1.0, r = 1.0, o = 1.0, b = 2.0;
  int t = 53, l = 1024, m = -1021;
  char rname[] = "BLAS_dfpinfo";

  if ((sizeof eps) == sizeof(float)) {
    t = 24;
    l = 128;
    m = -125;
  } else {
    t = 53;
    l = 1024;
    m = -1021;
  }

  /* for (i = 0; i < t; ++i) eps *= half; */
  eps = BLAS_dpow_di( b, -t );
  /* for (i = 0; i >= m; --i) r *= half; */
  r = BLAS_dpow_di( b, m-1 );

  o -= eps;
  /* for (i = 0; i < l; ++i) o *= b; */
  o = (o * BLAS_dpow_di( b, l-1 )) * b;

  switch (cmach) {
    case blas_eps: return eps;
    case blas_sfmin: return r;
    default:
      BLAS_error( rname, -1, cmach, 0 );
      break;
  }
  return 0.0;
}

static int check_factorization(int, PLASMA_Complex64_t*, PLASMA_Complex64_t*, int, int , double);
static int check_solution(int, int, PLASMA_Complex64_t*, int, PLASMA_Complex64_t*, PLASMA_Complex64_t*, int, double);
static int check_estimator(PLASMA_enum, int, PLASMA_Complex64_t *, int,
                           PLASMA_Complex64_t *, double, double,
                           double);

int testing_zposv(int argc, char **argv)
{

    /* Check for number of arguments*/
    if (argc != 4){
        USAGE("POSV", "N LDA NRHS LDB",
              "   - N    : the size of the matrix\n"
              "   - LDA  : leading dimension of the matrix A\n"
              "   - NRHS : number of RHS\n"
              "   - LDB  : leading dimension of the RHS B\n");
        return -1;
    }

    int N     = atoi(argv[0]);
    int LDA   = atoi(argv[1]);
    int NRHS  = atoi(argv[2]);
    int LDB   = atoi(argv[3]);
    double eps;
    int info_solution, info_factorization;
    int u, trans1, trans2;

    PLASMA_Complex64_t *A1   = (PLASMA_Complex64_t *)malloc(LDA*N*sizeof(PLASMA_Complex64_t));
    PLASMA_Complex64_t *A2   = (PLASMA_Complex64_t *)malloc(LDA*N*sizeof(PLASMA_Complex64_t));
    PLASMA_Complex64_t *B1   = (PLASMA_Complex64_t *)malloc(LDB*NRHS*sizeof(PLASMA_Complex64_t));
    PLASMA_Complex64_t *B2   = (PLASMA_Complex64_t *)malloc(LDB*NRHS*sizeof(PLASMA_Complex64_t));

    /* Check if unable to allocate memory */
    if ((!A1)||(!A2)||(!B1)||(!B2)){
        printf("Out of Memory \n ");
        return -2;
    }

    eps = BLAS_dfpinfo( blas_eps );

    for(u=0; u<2; u++) {

        trans1 = uplo[u] == PlasmaUpper ? PlasmaConjTrans : PlasmaNoTrans;
        trans2 = uplo[u] == PlasmaUpper ? PlasmaNoTrans : PlasmaConjTrans;

        /*-------------------------------------------------------------
         *  TESTING ZPOSV
         */

        /* Initialize A1 and A2 for Symmetric Positif Matrix */
        PLASMA_zplghe( (double)N, N, A1, LDA, 51 );
        PLASMA_zlacpy( PlasmaUpperLower, N, N, A1, LDA, A2, LDA );

        /* Initialize B1 and B2 */
        PLASMA_zplrnt( N, NRHS, B1, LDB, 371 );
        PLASMA_zlacpy( PlasmaUpperLower, N, NRHS, B1, LDB, B2, LDB );

        printf("\n");
        printf("------ TESTS FOR PLASMA ZPOSV ROUTINE -------  \n");
        printf("            Size of the Matrix %d by %d\n", N, N);
        printf("\n");
        printf(" The matrix A is randomly generated for each test.\n");
        printf("============\n");
        printf(" The relative machine precision (eps) is to be %e \n", eps);
        printf(" Computational tests pass if scaled residuals are less than 60.\n");

        /* PLASMA ZPOSV */
        PLASMA_zposv(uplo[u], N, NRHS, A2, LDA, B2, LDB);

        /* Check the factorization and the solution */
        info_factorization = check_factorization( N, A1, A2, LDA, uplo[u], eps);
        info_solution = check_solution(N, NRHS, A1, LDA, B1, B2, LDB, eps);

        if ( (info_solution == 0) && (info_factorization == 0) ) {
            printf("***************************************************\n");
            printf(" ---- TESTING ZPOSV(%s) ...................... PASSED !\n", uplostr[u]);
            printf("***************************************************\n");
        }
        else {
            printf("***************************************************\n");
            printf(" - TESTING ZPOSV(%s) ... FAILED !\n", uplostr[u]);
            printf("***************************************************\n");
        }

        /*-------------------------------------------------------------
         *  TESTING ZPOTRF + ZPOTRS
         */

        /* Initialize A1 and A2 for Symmetric Positif Matrix */
        PLASMA_zplghe( (double)N, N, A1, LDA, 51 );
        PLASMA_zlacpy( PlasmaUpperLower, N, N, A1, LDA, A2, LDA );

        /* Initialize B1 and B2 */
        PLASMA_zplrnt( N, NRHS, B1, LDB, 371 );
        PLASMA_zlacpy( PlasmaUpperLower, N, NRHS, B1, LDB, B2, LDB );

        /* Plasma routines */
        PLASMA_zpotrf(uplo[u], N, A2, LDA);
        PLASMA_zpotrs(uplo[u], N, NRHS, A2, LDA, B2, LDB);

        printf("\n");
        printf("------ TESTS FOR PLASMA ZPOTRF + ZPOTRS ROUTINE -------  \n");
        printf("            Size of the Matrix %d by %d\n", N, N);
        printf("\n");
        printf(" The matrix A is randomly generated for each test.\n");
        printf("============\n");
        printf(" The relative machine precision (eps) is to be %e \n", eps);
        printf(" Computational tests pass if scaled residuals are less than 60.\n");

        /* Check the factorization and the solution */
        info_factorization = check_factorization( N, A1, A2, LDA, uplo[u], eps);
        info_solution = check_solution(N, NRHS, A1, LDA, B1, B2, LDB, eps);

        if ((info_solution == 0)&(info_factorization == 0)){
            printf("***************************************************\n");
            printf(" ---- TESTING ZPOTRF + ZPOTRS (%s)............ PASSED !\n", uplostr[u]);
            printf("***************************************************\n");
        }
        else{
            printf("****************************************************\n");
            printf(" - TESTING ZPOTRF + ZPOTRS (%s)... FAILED !\n", uplostr[u]);
            printf("****************************************************\n");
        }

        /*-------------------------------------------------------------
         *  TESTING ZPOTRF + ZPTRSM + ZTRSM
         */

        /* Initialize A1 and A2 for Symmetric Positif Matrix */
        PLASMA_zplghe( (double)N, N, A1, LDA, 51 );
        PLASMA_zlacpy( PlasmaUpperLower, N, N, A1, LDA, A2, LDA );

        /* Initialize B1 and B2 */
        PLASMA_zplrnt( N, NRHS, B1, LDB, 371 );
        PLASMA_zlacpy( PlasmaUpperLower, N, NRHS, B1, LDB, B2, LDB );

        /* PLASMA routines */
        PLASMA_zpotrf(uplo[u], N, A2, LDA);
        PLASMA_ztrsm(PlasmaLeft, uplo[u], trans1, PlasmaNonUnit,
                     N, NRHS, 1.0, A2, LDA, B2, LDB);
        PLASMA_ztrsm(PlasmaLeft, uplo[u], trans2, PlasmaNonUnit,
                     N, NRHS, 1.0, A2, LDA, B2, LDB);

        printf("\n");
        printf("------ TESTS FOR PLASMA ZPOTRF + ZTRSM + ZTRSM  ROUTINE -------  \n");
        printf("            Size of the Matrix %d by %d\n", N, N);
        printf("\n");
        printf(" The matrix A is randomly generated for each test.\n");
        printf("============\n");
        printf(" The relative machine precision (eps) is to be %e \n", eps);
        printf(" Computational tests pass if scaled residuals are less than 60.\n");

        /* Check the factorization and the solution */
        info_factorization = check_factorization( N, A1, A2, LDA, uplo[u], eps);
        info_solution = check_solution(N, NRHS, A1, LDA, B1, B2, LDB, eps);

        if ((info_solution == 0)&(info_factorization == 0)){
            printf("***************************************************\n");
            printf(" ---- TESTING ZPOTRF + ZTRSM + ZTRSM (%s)..... PASSED !\n", uplostr[u]);
            printf("***************************************************\n");
        }
        else{
            printf("***************************************************\n");
            printf(" - TESTING ZPOTRF + ZTRSM + ZTRSM (%s)... FAILED !\n", uplostr[u]);
            printf("***************************************************\n");
        }

        /*-------------------------------------------------------------
         *  TESTING ZPOCON on the last call
         */
        {
            double Anorm = PLASMA_zlanhe( PlasmaOneNorm, uplo[u], N, A1, LDA );
            double Acond;

            info_solution = PLASMA_zpocon(uplo[u], N, A2, LDA, Anorm, &Acond);
            if ( info_solution == 0 ) {
                info_solution = check_estimator(uplo[u], N, A1, LDA, A2, Anorm, Acond, eps);
            } else {
                printf(" PLASMA_zpocon returned info = %d\n", info_solution );
            }
            if ((info_solution == 0)){
                printf("***************************************************\n");
                printf(" ---- TESTING ZPOTRF + ZPOCON (%s) ........... PASSED !\n", uplostr[u]);
                printf("***************************************************\n");
            }
            else{
                printf("**************************************************\n");
                printf(" - TESTING ZPOTRF + ZPOCON (%s) ... FAILED !\n", uplostr[u]);
                printf("**************************************************\n");
            }
        }
    }

    free(A1); free(A2); free(B1); free(B2);

    return 0;
}

/*------------------------------------------------------------------------
 *  Check the factorization of the matrix A2
 */
static int check_factorization(int N, PLASMA_Complex64_t *A1, PLASMA_Complex64_t *A2, int LDA, int uplo, double eps)
{
    double Anorm, Rnorm;
    PLASMA_Complex64_t alpha;
    int info_factorization;
    int i,j;

    PLASMA_Complex64_t *Residual = (PLASMA_Complex64_t *)malloc(N*N*sizeof(PLASMA_Complex64_t));
    PLASMA_Complex64_t *L1       = (PLASMA_Complex64_t *)malloc(N*N*sizeof(PLASMA_Complex64_t));
    PLASMA_Complex64_t *L2       = (PLASMA_Complex64_t *)malloc(N*N*sizeof(PLASMA_Complex64_t));
    double *work              = (double *)malloc(N*sizeof(double));

    memset((void*)L1, 0, N*N*sizeof(PLASMA_Complex64_t));
    memset((void*)L2, 0, N*N*sizeof(PLASMA_Complex64_t));

    alpha= 1.0;

    LAPACKE_zlacpy_work(LAPACK_COL_MAJOR,' ', N, N, A1, LDA, Residual, N);

    /* Dealing with L'L or U'U  */
    if (uplo == PlasmaUpper){
        LAPACKE_zlacpy_work(LAPACK_COL_MAJOR,'u', N, N, A2, LDA, L1, N);
        LAPACKE_zlacpy_work(LAPACK_COL_MAJOR,'u', N, N, A2, LDA, L2, N);
        cblas_ztrmm(CblasColMajor, CblasLeft, CblasUpper, CblasConjTrans, CblasNonUnit, N, N, CBLAS_SADDR(alpha), L1, N, L2, N);
    }
    else{
        LAPACKE_zlacpy_work(LAPACK_COL_MAJOR,'l', N, N, A2, LDA, L1, N);
        LAPACKE_zlacpy_work(LAPACK_COL_MAJOR,'l', N, N, A2, LDA, L2, N);
        cblas_ztrmm(CblasColMajor, CblasRight, CblasLower, CblasConjTrans, CblasNonUnit, N, N, CBLAS_SADDR(alpha), L1, N, L2, N);
    }

    /* Compute the Residual || A -L'L|| */
    for (i = 0; i < N; i++)
        for (j = 0; j < N; j++)
           Residual[j*N+i] = L2[j*N+i] - Residual[j*N+i];

    BLAS_zge_norm( blas_colmajor, blas_inf_norm, N, N, Residual, N, &Rnorm );
    BLAS_zge_norm( blas_colmajor, blas_inf_norm, N, N, A1, LDA, &Anorm );

    printf("============\n");
    printf("Checking the Cholesky Factorization \n");
    printf("-- ||L'L-A||_oo/(||A||_oo.N.eps) = %e \n",Rnorm/(Anorm*N*eps));

    if ( isnan(Rnorm/(Anorm*N*eps)) || isinf(Rnorm/(Anorm*N*eps)) || (Rnorm/(Anorm*N*eps) > 60.0) ){
        printf("-- Factorization is suspicious ! \n");
        info_factorization = 1;
    }
    else{
        printf("-- Factorization is CORRECT ! \n");
        info_factorization = 0;
    }

    free(Residual); free(L1); free(L2); free(work);

    return info_factorization;
}

/*------------------------------------------------------------------------
 *  Check the accuracy of the solution of the linear system
 */

static int check_solution(int N, int NRHS, PLASMA_Complex64_t *A1, int LDA, PLASMA_Complex64_t *B1, PLASMA_Complex64_t *B2, int LDB, double eps )
{
    int info_solution;
    double Rnorm, Anorm, Xnorm, Bnorm, result;
    PLASMA_Complex64_t alpha, beta;
    double *work = (double *)malloc(N*sizeof(double));

    alpha = 1.0;
    beta  = -1.0;

    BLAS_zge_norm( blas_colmajor, blas_inf_norm, N, NRHS, B2, LDB, &Xnorm );
    BLAS_zge_norm( blas_colmajor, blas_inf_norm, N, N,    A1, LDA, &Anorm );
    BLAS_zge_norm( blas_colmajor, blas_inf_norm, N, NRHS, B1, LDB, &Bnorm );

    cblas_zgemm(CblasColMajor, CblasNoTrans, CblasNoTrans, N, NRHS, N, CBLAS_SADDR(alpha), A1, LDA, B2, LDB, CBLAS_SADDR(beta), B1, LDB);
    BLAS_zge_norm( blas_colmajor, blas_inf_norm, N, NRHS, B1, LDB, &Rnorm );

    if (getenv("PLASMA_TESTING_VERBOSE"))
      printf( "||A||_oo=%f\n||X||_oo=%f\n||B||_oo=%f\n||A X - B||_oo=%e\n", Anorm, Xnorm, Bnorm, Rnorm );

    result = Rnorm / ( (Anorm*Xnorm+Bnorm)*N*eps ) ;
    printf("============\n");
    printf("Checking the Residual of the solution \n");
    printf("-- ||Ax-B||_oo/((||A||_oo||x||_oo+||B||_oo).N.eps) = %e \n", result);

    if (  isnan(Xnorm) || isinf(Xnorm) || isnan(result) || isinf(result) || (result > 60.0) ) {
        printf("-- The solution is suspicious ! \n");
        info_solution = 1;
     }
    else{
        printf("-- The solution is CORRECT ! \n");
        info_solution = 0;
    }

    free(work);

    return info_solution;
}

/*------------------------------------------------------------------------
 *  Check the accuracy of the condition estimator
 */
static int check_estimator(PLASMA_enum uplo, int N, PLASMA_Complex64_t *A1, int LDA,
                           PLASMA_Complex64_t *A2, double Anorm, double Acond,
                           double eps)
{
    int info_solution;
    double result, Acond_lapack;
    double invcond, invcond_lapack;

    info_solution = LAPACKE_zpocon(LAPACK_COL_MAJOR, lapack_const(uplo), N, A2, LDA, Anorm, &Acond_lapack);

    if ( info_solution != 0 ) {
        printf(" PLASMA_zgecon returned info = %d\n", info_solution );
        return info_solution;
    }

    invcond_lapack = 1. / ( Acond_lapack );
    invcond        = 1. / ( Acond );

    printf("============\n");
    printf("Checking the condition number \n");
    printf("-- Acond_plasma = %e, Acond_lapack = %e \n"
           "-- Ainvcond_plasma = %e, Ainvcond_lapack = %e \n",
           Acond, Acond_lapack, invcond, invcond_lapack );

    result = fabs( Acond_lapack - Acond ) / eps;
    if ( result > 60. ) {
        printf("-- The solution is suspicious ! \n");
        info_solution = 1;
     }
    else{
        printf("-- The solution is CORRECT ! \n");
        info_solution = 0;
    }
    return info_solution;
}


text:
/* =======================================
 * Code for 8.06 Term Paper (Spring 2015)
 *
 * Author : To Chin Yu (ytc@mit.edu)
 * Purpose: 1D lattice simulation
 * V(x) = alpha*x^2 + lambda*x^4
 * =======================================
 */
#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#define sq(x) pow((x),2.0)
#define cb(x) pow((x),3.0)
#define qd(x) pow((x),4.0)

// GSL Random Number Generator (https://www.gnu.org/software/gsl/)
#include <gsl/gsl_rng.h>
const gsl_rng_type * T; gsl_rng * r;

// Parameters
#define N_lat   5000  // Lattice size
#define a       0.1   // Lattice spacing
#define N_therm 500   // Thermalization
#define N_meas  1000  // Number of configuration to measure
#define N_skip  100   // Decorrelation
#define alpha   0     // Strength of potential (harmonic)
#define lambda  1     // Strength of potential (anharmonic)
#define epsilon 0.2   // Size of fluctuation

#define N_bin 1000
#define dx 0.01
#define x_range 5.0 // x in [-range, range]

// Globals
double x[N_lat];      // Lattice
int acc,tot;          // Acceptance rate = acc/tot
double E0;            // Ground state energy
double C[N_lat];      // Correlator
double psi[N_bin];    // Wavefunction

// Output Files
FILE *fE0, *fCorRatio, *fpsi;

// Prototypes
void setup();
void cleanup();
void update();

int main(void){
  setup();

  // Thermalize
  for(int n=0;n<N_therm;n++){
    update();

    // Plot E0 during thermalization
    if(0){
      E0=0;
      for(int i=0;i<N_lat;i++){
        E0 += ( 2*x[i]*lambda*cb(x[i]) + lambda*qd(x[i]) )/N_lat;
      }
      fprintf(fE0,"%f\n", E0);
    }
  }

  // Do updates and measurements
  for(int n=0;n<N_meas;n++){
    // Decorrelate
    for(int m=0;m<N_skip;m++){ update(); }
    // Measurements:
    // 1. E0

    // Plot E0
    E0=0;
    for(int i=0;i<N_lat;i++){
      E0 += ( 2*x[i]*lambda*cb(x[i]) + lambda*qd(x[i]) )/N_lat;
    }
    fprintf(fE0,"%f\n", E0);

    // 2. Correlator
    for(int i=0;i<N_lat;i++){
      for(int T=0;T<N_lat;T++){
        C[T] += x[i]*x[(i+T)%N_lat] / (N_lat*N_meas);
      }
    }

    // 3. Wavefunction
    for(int i=0;i<N_lat;i++){
      if(abs(x[i])<x_range){
        psi[(int)floor((x[i]+x_range)/dx)]+=1;
      }
    }

  }

  // Plot correlator ratio
  for(int T=0;T<N_lat;T++){
    fprintf(fCorRatio,"%f \n", C[T]/C[(T+1)%N_lat]); // for large T, this gives a(E1-E0)
  }

  // Plot wavefunction
  for(int i=0;i<N_bin;i++){
    fprintf(fpsi,"%f, %f\n", i*dx-x_range, sqrt(psi[i]/(dx*N_meas*N_lat)));
  }

  printf("Acceptance rate: %f \n", (double)acc/tot);

  cleanup();
  return 0;
}

// Calculate difference in action
double dS(int i, double x_old){
  return (x[i]*(x[i]-x[(i+1)%N_lat]-x[(i-1+N_lat)%N_lat])/a
          + alpha*a*sq(x[i]) + lambda*a*qd(x[i]))
          - (x_old*(x_old-x[(i+1)%N_lat]-x[(i-1+N_lat)%N_lat])/a
               + alpha*a*sq(x_old) + lambda*a*qd(x_old));
}

void update(){
  double x_old;
  double d; // Difference in action
  for(int i=0;i<N_lat;i++){
    x_old=x[i];
    x[i]=x[i] + epsilon * (gsl_rng_uniform(r)-0.5)*2;
    acc++; tot++;
    d=dS(i,x_old);
    if(d>0){
      if(exp(-d)<gsl_rng_uniform(r)){ x[i]=x_old; acc--; } // Revert to old value
    }
  }
}

void setup(){
  // Cold start
  for(int i=0;i<N_lat;i++){ x[i]=0; C[i]=0; }

  // Initialize variables
  acc = 0;
  tot = 0;
  E0 = 0;
  for(int i=0;i<N_bin;i++){ psi[i]=0; }

  // Create Files
  fE0=fopen("E0.txt", "w+");
  fCorRatio=fopen("CorRatio.txt", "w+");
  fpsi=fopen("psi.txt", "w+");

  // Setup random environment
  gsl_rng_env_setup();
  T = gsl_rng_mt19937;
  r = gsl_rng_alloc (T);
}

void cleanup(){
  gsl_rng_free (r);

  // Close Files
  fclose(fE0);
  fclose(fCorRatio);
  fclose(fpsi);
}

max_stars_repo_name: uchuutamashi/lqcd

max_issues_repo_name: uchuutamashi/lqcd

max_forks_repo_name: uchuutamashi/lqcd


text:
 // MIT License

 // Copyright (c) [2017] [Vinay Yuvashankar]

 // Permission is hereby granted, free of charge, to any person obtaining a copy
 // of this software and associated documentation files (the "Software"), to deal
 // in the Software without restriction, including without limitation the rights
 // to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
 // copies of the Software, and to permit persons to whom the Software is
 // furnished to do so, subject to the following conditions:

 // The above copyright notice and this permission notice shall be included in all
 // copies or substantial portions of the Software.

 // THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
 // IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
 // FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
 // AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
 // LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
 // OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
 // SOFTWARE.

#include <assert.h>
#include "wavelet.h"
#include <inttypes.h>
#include <gsl/gsl_statistics.h>
#include <omp.h>

int main(int argc, char const *argv[])
{
    if (argc < 2)
    {
        printf("Not enough input arguments please provide a BDF file\n");
        return(-1);
    }

    double t = omp_get_wtime();

    //Variable Declarations
    int handle,
        numberOfChannels,
        openFlag,
        channel,
        filteredTriggerNumber,
        readFlag,
        writeFlag,
        samplesToRead,
        numberOfTriggers,
        read_channel,
        min_i,
        max_i,
        J; //needed by Wavelet

    int *rawStatus,
        *buffer,
        *value,
        *filteredBuffer;

    double  sampleFrequency,
            dt;

    int64_t numberOfRecords;

    double  *tempBuffer,
            *data,
            *result, //The Final Result
            *scales,
            *frequency; //The Corrosponding Frequencies

    char    graph_title[256],
            file_name[256];

    int64_t * triggerList;

    struct edf_hdr_struct edfHeader;

    // Debug File ... just in case.
    // FILE *debug_file = fopen("Main_debug.log", "w");
    // assert(debug_file != NULL);

    //Opening and reading the file into the edfHeader input the first argument is the input file.
    openFlag = OpenFile(argv[2], &edfHeader);
    read_channel = atoi(argv[1]);
    assert( openFlag == 0 );

    //Get File Information
    handle           = edfHeader.handle;
    sampleFrequency  = (( double ) edfHeader.signalparam[1].smp_in_datarecord /
                        ( double ) edfHeader.datarecord_duration              ) * EDFLIB_TIME_DIMENSION;
    numberOfChannels = edfHeader.edfsignals;
    numberOfRecords  = edfHeader.signalparam[numberOfChannels - 1].smp_in_file;
    channel          = numberOfChannels - 1; //The status channel.
    samplesToRead    = (PRE_EVENT_TIME + POST_EVENT_TIME) * sampleFrequency;

    //Generate Constants for the CWT
    dt               = 1.0/sampleFrequency;
    min_i            = floor( ( log2( (W_0) / (dt * 2 * M_PI * sampleFrequency/2.0) ) )/D_J);
    max_i            = floor( ( log2( (W_0) / (dt * 2 * M_PI * MIN_FREQUENCY) ) )      /D_J);
    assert(min_i > 0); assert(max_i > 0);

    J     = (int) max_i - min_i;

    /*Allocate Necessary Memory*/
    rawStatus   =    (int*)     malloc( numberOfRecords  * sizeof   (int)     );
    triggerList =    (int64_t*) malloc( MAXIMUM_TRIGGERS * sizeof   (int64_t) );
    tempBuffer  =    (double*)  malloc( samplesToRead    * sizeof   (double)  );
    assert(rawStatus != NULL); assert(triggerList!= NULL); assert(tempBuffer!= NULL);

    //Read the status Signal --> Output to rawStatus
    readFlag = edfread_digital_samples(handle, channel, numberOfRecords, rawStatus);
    assert (readFlag != 1);

    //Find and parse the file and find the triggers -->  Output to TriggerList
    numberOfTriggers = FindTriggers(rawStatus, numberOfRecords, triggerList);
    assert (numberOfTriggers != -1);

    //Allocating Read Buffer
    buffer         = (int *) malloc( numberOfTriggers * sizeof(int) );
    filteredBuffer = (int *) malloc( numberOfTriggers * sizeof(int) );
    assert (buffer!= NULL); assert(filteredBuffer != NULL);

    scales         = (double*) malloc ( J * sizeof(double) );
    frequency      = (double*) malloc ( J * sizeof(double) );
    assert(scales != NULL); assert(frequency != NULL);

    //Wavelet Memory Allocations
    data           = (double*) malloc( samplesToRead * filteredTriggerNumber * sizeof(double));
    result         = (double*) malloc( J * samplesToRead * sizeof(double));
    assert(result != NULL); assert(data != NULL);

    /*Main Program Begins*/
    GenerateScalesAndFrequency(min_i, max_i, dt, scales, frequency);

    //Read the status channel when there is a trigger and put it in the buffer
    value = (int*) malloc( sizeof(int) );
    for (int i = 0; i < numberOfTriggers; ++i)
    {
        edfseek(handle, channel, triggerList[i], EDFSEEK_SET);

        readFlag = edfread_digital_samples(handle, channel, 1, value);
        assert(readFlag != -1);
        buffer[i] = *value;
    }

    //Filter the Triggers to what you want.
    filteredTriggerNumber = FilterTriggers(1, 2, numberOfTriggers, triggerList,
                                            buffer, filteredBuffer);
    printf("Number of Filtered Triggers Found: %d\n", filteredTriggerNumber);

    printf("Beginning Wavelet Analysis\n");
    for (int i = 0; i < filteredTriggerNumber; ++i)
    {
        edfseek(handle, 0, filteredBuffer[i], EDFSEEK_SET);

        readFlag = edfread_physical_samples(handle, read_channel, samplesToRead, tempBuffer);
        // TestCases( tempBuffer, 5, 16.0 , sampleFrequency, samplesToRead);

        // Preform a Z-Score on the read data.
        CleanData(tempBuffer, samplesToRead);

        for (int j = 0; j < samplesToRead; ++j)
        {
            data[i * samplesToRead + j] = tempBuffer[j];
        }
    }

    ERSP (data, scales, sampleFrequency, samplesToRead,
                J, filteredTriggerNumber, 1,
                result);

    printf("Writing to File\n");
    sprintf(graph_title, "Channel %d", read_channel);
    sprintf(file_name,   "Channel_%d.log", read_channel);

    writeFlag = WriteFile(result, frequency, J, samplesToRead, sampleFrequency,
                          file_name);
    WriteGnuplotScript(graph_title, file_name);

    assert(writeFlag != -1);

    printf("Freeing Memory and closing files\n");
    edfclose_file(handle);

    free(data);   free(result);
    free(scales); free(frequency);
    free(buffer); free(filteredBuffer);
    free(rawStatus); free(triggerList); free(tempBuffer);

    t = omp_get_wtime() - t;

    printf("Execution Time = %f\n", t);

    return 0;
}

max_stars_repo_name: yuvashankar/Research

max_issues_repo_name: yuvashankar/Research

max_forks_repo_name: yuvashankar/Research


text:
#include <math.h>
#include <stdlib.h>
#include <gsl/gsl_linalg.h>
#include "userFunc.h"

#define G        6.67384e-8    /* 2010 CODATA value in CGS units */
#define QCRIT    1.0           /* Q value at which GI turns on */
#define QDAMP    10.0          /* Rate of damping of alpha above Qcrit */
#define ALPHAMIN 1.0e-6        /* Minimum allowed alpha value */

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

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

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

  static int nr = 0;
  static double rmin = -1.0, rmax = -1.0;
  static double *x = NULL;
  static double *h1 = NULL;
  static double *h0 = NULL;
  static double *H = NULL;
  static double *Q = NULL;
  static gsl_vector *ud_g = NULL;
  static gsl_vector *ld_g = NULL;
  static gsl_vector *diag_g = NULL;
  static gsl_vector *rhs_g = NULL;
  static gsl_vector *tau_g = NULL;
  double eta = ((double *) params)[0];
  double chi = ((double *) params)[1];
  double tQ = ((double *) params)[2];
  double u, dlnx, sL, s, sR, dsdx, beta, dbetadx;
  double mdot;
  double dlnQdlnt;
  int i;

  /* Allocate memory if necessary */
  if (grd->nr != nr) {
    nr = grd->nr;
    rmin = -1.0;
    rmax = -1.0;
    /* If previously allocated, free */
    if (x != NULL) free(x);
    if (h1 != NULL) free(h1);
    if (h0 != NULL) free(h0);
    if (H != NULL) free(H);
    if (Q != NULL) free(Q);
    if (ud_g != NULL) gsl_vector_free(ud_g);
    if (ld_g != NULL) gsl_vector_free(ld_g);
    if (diag_g != NULL) gsl_vector_free(diag_g);
    if (rhs_g != NULL) gsl_vector_free(rhs_g);
    if (tau_g != NULL) gsl_vector_free(tau_g);
    /* Allocate new memory */
    if (!(x = calloc(grd->nr, sizeof(double)))) {
      fprintf(stderr, "Error: unable to allocate memory in userAlpha!\n");
      exit(1);
    }
    if (!(h1 = calloc(grd->nr, sizeof(double)))) {
      fprintf(stderr, "Error: unable to allocate memory in userAlpha!\n");
      exit(1);
    }
    if (!(h0 = calloc(grd->nr, sizeof(double)))) {
      fprintf(stderr, "Error: unable to allocate memory in userAlpha!\n");
      exit(1);
    }
    if (!(H = calloc(grd->nr, sizeof(double)))) {
      fprintf(stderr, "Error: unable to allocate memory in userAlpha!\n");
      exit(1);
    }
    if (!(Q = calloc(grd->nr, sizeof(double)))) {
      fprintf(stderr, "Error: unable to allocate memory in userAlpha!\n");
      exit(1);
    }
    if (!(ud_g = gsl_vector_alloc(grd->nr+1))) {
      fprintf(stderr, "Error: unable to allocate memory in userAlpha!\n");
      exit(1);
    }
    if (!(ld_g = gsl_vector_alloc(grd->nr+1))) {
      fprintf(stderr, "Error: unable to allocate memory in userAlpha!\n");
      exit(1);
    }
    if (!(diag_g = gsl_vector_alloc(grd->nr+2))) {
      fprintf(stderr, "Error: unable to allocate memory in userAlpha!\n");
      exit(1);
    }
    if (!(rhs_g = gsl_vector_alloc(grd->nr+2))) {
      fprintf(stderr, "Error: unable to allocate memory in userAlpha!\n");
      exit(1);
    }
    if (!(tau_g = gsl_vector_alloc(grd->nr+2))) {
      fprintf(stderr, "Error: unable to allocate memory in userAlpha!\n");
      exit(1);
    }
  }
  /* Initailize x if necessary */
  if ((rmin != grd->r_h[0]) || (rmax != grd->r_h[grd->nr+1])) {
    for (i=0; i<grd->nr; i++) x[i] = grd->r_g[i+1]/grd->r_h[grd->nr];
    rmin = grd->r_h[0];
    rmax = grd->r_h[grd->nr+1];
  }

  /* Compute coefficients of torque equation */
  dlnx = log(grd->r_g[1]/grd->r_g[0]);
  for (i=0; i<grd->nr; i++) {

    /* Define dimensionless velocity */
    u = grd->vphi_g[i+1] / grd->vphi_h[grd->nr];

    /* Note: compute derivatives of s using one-sided differences at
       edge cells, centered-limited differences elsewhere. */
    if (i==0) {
      s = sL = sqrt(pres[0]/col[0])/grd->vphi_h[grd->nr];
      sR = sqrt(pres[1]/col[1])/grd->vphi_h[grd->nr];
      dsdx = (sR - sL) / (x[i]*dlnx);
    } else if (i==grd->nr-1) {
      sL = sqrt(pres[grd->nr-2]/col[grd->nr-2])/grd->vphi_h[grd->nr];
      s = sR = sqrt(pres[grd->nr-1]/col[grd->nr-1])/grd->vphi_h[grd->nr];
      dsdx = (sR - sL) / (x[i]*dlnx);
    } else {
      sL = sqrt(pres[i+1]/col[i+1])/grd->vphi_h[grd->nr];
      s = sqrt(pres[i]/col[i])/grd->vphi_h[grd->nr];
      sR = sqrt(pres[i-1]/col[i-1])/grd->vphi_h[grd->nr];
      dsdx = ((sL-s)*(s-sR)) > 0 ? (sL-sR)/(2.0*x[i]*dlnx) : 0;
    }

    /* Compute derivative of beta using centered difference */
    beta = grd->beta_g[i+1];
    dbetadx = (grd->beta_h[i+1] - grd->beta_h[i]) / (x[i]*dlnx);

    /* Compute torque equation coefficients */
    h1[i] = -(5.0*(beta+1)*x[i]*dsdx + 2.0*s*(beta+SQR(beta)+x[i]*dbetadx)) /
      (2.0*(beta+1)*s*x[i]);
    h0[i] = (SQR(beta)-1) * SQR(u/(x[i]*s)) / 2.0;
    Q[i] = sqrt(2.0*(1+beta))*grd->vphi_g[i+1]/grd->r_g[i+1] *
      s*grd->vphi_h[grd->nr]/(M_PI*G*col[i]);
    dlnQdlnt = Q[i] < QCRIT ? -(exp(QCRIT/Q[i]) - exp(1.0))*u/x[i] : 0.0;
    H[i] = G*grd->r_h[grd->nr]*u*(1+beta)*col[i] *
      (2*M_PI*u*eta - 3*dlnQdlnt/tQ) / (2.0*SQR(grd->vphi_h[grd->nr])*chi);
  }

  /* Load matrix elements for tridiagonal solve */
  for (i=0; i<grd->nr; i++) {
    gsl_vector_set(ld_g, i, 1.0/SQR(x[i]*dlnx) +
		   1.0/(2.0*SQR(x[i])*dlnx) - h1[i]/(2.0*x[i]*dlnx));
    gsl_vector_set(diag_g, i+1, -2.0/SQR(x[i]*dlnx) + h0[i]);
    gsl_vector_set(ud_g, i+1, 1.0/SQR(x[i]*dlnx) -
		   1.0/(2.0*SQR(x[i])*dlnx) + h1[i]/(2.0*x[i]*dlnx));
    gsl_vector_set(rhs_g, i+1, H[i]);
  }

  /* Load boundary conditions into matrix */
  /* Inner BC: tau = -x in ghost cell */
  gsl_vector_set(diag_g, 0, 1);
  gsl_vector_set(ud_g, 0, 0);
  gsl_vector_set(rhs_g, 0, -grd->r_g[0]/grd->r_h[grd->nr]);
  /* Outer BC: tau' = -beta-1 */
  gsl_vector_set(ld_g, grd->nr, -1.0/dlnx);
  gsl_vector_set(diag_g, grd->nr+1, 1.0/dlnx);
  gsl_vector_set(rhs_g, grd->nr+1, -grd->beta_g[grd->nr+1]-1);

  /* Solve matrix equation to get dimensionless torque */
  gsl_linalg_solve_tridiag(diag_g, ud_g, ld_g, rhs_g, tau_g);

  /* Use torque to compute alpha */
  mdot = chi*SQR(grd->vphi_h[grd->nr])*grd->vphi_h[grd->nr] / G;
  for (i=0; i<grd->nr; i++) {
    alpha[i] = -gsl_vector_get(tau_g, i+1)*mdot*grd->vphi_h[grd->nr] *
      grd->r_h[grd->nr] /
      (2.0*M_PI*SQR(grd->r_g[i+1])*pres[i]);
    if (Q[i] > QCRIT) alpha[i] *= exp(-QDAMP*(Q[i]-QCRIT));
    if (alpha[i] < 0.0) alpha[i] = ALPHAMIN;
  }
}

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

void
userIntEnSrc(const double t, const double dt, const grid *grd,
	     const double *col, const double *pres, const double *eInt,
	     const double *gamma, const double *delta,
	     void *params,
	     double *intEnSrc) {
  /* Cooling rate = eta Sigma sigma^2 Omega = eta P vphi/r */
  int i;
  double eta = ((double *) params)[0];

  for (i=0; i<grd->nr; i++)
    intEnSrc[i] = -eta * pres[i] * grd->vphi_g[i+1] / grd->r_g[i+1];
}

void
userIBC(const double t, const double dt, const grid *grd,
	const double *col, const double *pres, const double *eInt,
	const double *gamma, const double *delta,
	const pres_bc_type ibc_pres, const enth_bc_type ibc_enth,
	void *params,
	double *ibc_pres_val, double *ibc_enth_val) {

  /* Set torque. If Q <= QCRIT in innermost cell, so that disk is
     gravitationally unstable at the center, dimensionless torque =
     -x, dimensional torque = -(r_in/R) Mdot vphi/R. If Q > QCRIT in
     innermost cell, scale down torque in the same way as in
     userAlpha. */
  double chi = ((double *) params)[1];
  double x, mdot, s, Q;

  s = sqrt(pres[0]/col[0])/grd->vphi_h[grd->nr];
  Q = sqrt(2.0*(1+grd->beta_g[1]))*grd->vphi_g[1]/grd->r_g[1] *
    s*grd->vphi_h[grd->nr]/(M_PI*G*col[0]);
  mdot = chi*SQR(grd->vphi_h[grd->nr])*grd->vphi_h[grd->nr] / G;
  x = grd->r_g[0] / grd->r_h[grd->nr];
  *ibc_pres_val = -x * mdot * grd->vphi_g[0] * grd->r_h[grd->nr];
  if (Q > QCRIT) *ibc_pres_val *= exp(-QDAMP*(Q-QCRIT));
  *ibc_enth_val = gamma[0]/(gamma[0]-1.0)*pres[0]/col[0];
}

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

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

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

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

max_stars_repo_name: franciscaconcha/amuse-vader

max_issues_repo_name: franciscaconcha/amuse-vader

max_forks_repo_name: franciscaconcha/amuse-vader


text:
#ifndef TREE_H
#define TREE_H
#include "node.h"
#include <string>
#include <vector>
#include <gsl/gsl_rng.h>
#include "data.h"

using namespace std;

namespace weakarg
{

extern gsl_rng * rng;

/**
    @brief Genealogy of the isolates under study
*/
class Tree
{
protected:
    Node * root;///<Root node of the genealogy
    std::vector<Node*> nodes;///<Vector of all the nodes in the genealogy
    int n;///<Number of isolates
    double ttotal;///<Sum of branch lengths
public:
    Tree(std::string newick,bool isFilename=true,bool forceages=true);///<Creates a tree from a Newick file
    Tree(int n);///< Creates a Coalescent tree by simulation
    Tree(Tree *intree);///< Creates a copy of a tree via newick (== slow)
    Tree(Data * data);///< Creates a UPGMA tree
    Tree(const Tree& t);///< Copy constructor
    void assign(const Tree& t);///< Copy a tree while allocating as little new storage as possible
    ~Tree();
    std::string newick(int p=6) const; ///<Returns a Newick description of the tree (precision options)
    std::string newickNoInternalLabels(int p=6) const; ///<Returns a Newick description of the tree (precision options)
    void makeFromNewick(std::string newick,bool forceages=false);///< Creates a tree from a newick string
    inline int getN() const
    {
        return n;
    }///<Returns the number of isolates
    inline Node* getNode(int i) const
    {
        return nodes[i];
    }///<Returns a node given its index
    inline Node* getRoot() const
    {
        return root;
    }///<Returns the root node
    inline double getDist(int i) const
    {
        return nodes[i]->getDist();
    }///<Returns the time of a node given its index
    double prior() const; ///<Returns log-prior of the tree
    inline double getTTotal() const
    {
        return ttotal;
    }///<Returns the sum of branch lengths
    void computeTTotal(); ///<Computes the sum of branch lengths
    int getPoint(double * dist, std::vector<int> * samplespace=NULL) const; ///<Returns a point chosen uniformly at random on the tree, or uniformally on a subspace of specified clonal edges
    std::vector<int> getAllChildren(int e);///<Returns a vector of the indices of the children of a given node
    std::vector<int> getAllSampledSeqs(int e);///<Returns a vector of all the observed sequences of the children of a given node

    void orderNodes(double dist=1);///<orders the list by age
    int orderNodes(int which, double dist);///<places a node in a new place in the list
    void swapFather(int a,int b);///<Swaps the father of two nodes
    void swapNode(int a, int b);///<Swaps two nodes in the list
    int getOldestReversedNode();///<Returns the oldest node that has been age reversed

    void testNodeAges() const; ///<Tests node ages
    double tavare() const;

    inline int otherChild(int par,int onechild){
	if(getNode(par)->getLeft()->getId()==onechild) return(getNode(par)->getRight()->getId());
	else if(getNode(par)->getRight()->getId()==onechild) return(getNode(par)->getLeft()->getId());
	else{cerr<<"Error in RecTree::otherChild: parent doesn't have a suitable child!"<<endl;throw("unknown child error");}
    };
    std::vector<int> getMinEdgeList(std::vector<int> seqs);
    int getNextGroup(std::vector<int> *seqs);
    bool isParentToOnly(int e, std::vector<int> seqs,std::vector<int> *which);

};

} // end namespace weakarg
#endif

ext: h


text:
/**
  * Copyright 2019 José Manuel Abuín Mosquera <josemanuel.abuin@usc.es>
  *
  * This file is part of Matrix Market Suite.
  *
  * Matrix Market Suite 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.
  *
  * Matrix Market Suite 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 Matrix Market Suite. If not, see <http://www.gnu.org/licenses/>.
  */

#include <cblas.h>
#include "ConjugateGradientSolver.h"

int ConjugateGradientSolver(unsigned long *II, unsigned long *J, double *A, unsigned long M, unsigned long N, unsigned long long nz, double *b, unsigned long M_Vector, unsigned long N_Vector, unsigned long long nz_vector, double *x, int iterationNumber) {

	//A*x=b

        double *Ap=(double *) malloc(nz_vector * sizeof(double));
        double *r=(double *) malloc(nz_vector * sizeof(double));
        double *p=(double *) malloc(nz_vector * sizeof(double));
	//double *x=(double *) calloc(nz_vector,sizeof(double));

	//r = b-A*x
	//If we take x=0 the init multiplication is avoided and r=b

	memcpy(r, b, N*sizeof(double));

	//p=r

	memcpy(p, r, N*sizeof(double));

	//rsold = r*r
	double rsold = cblas_ddot(N,r,1,r,1);

	int stop = 0;

	double alphaCG = 0.0;

	double rsnew = 0.0;
	unsigned long k = 0;

	unsigned long maxIterations = M*2;

	if(iterationNumber != 0 ){
		maxIterations = iterationNumber;
	}

	//int i  = 0;

	while(!stop){

		//Ap=A*p
		// for(i=0; i<M; i++){
		//	Ap[i] = 0.0;
		// }
		cblas_dgemv(CblasRowMajor, CblasNoTrans, M,N , 1.0, A, N, p, 1, 0.0, Ap, 1);

		//alphaCG=rsold/(p'*Ap)
		alphaCG = rsold/cblas_ddot(N,p,1,Ap,1);

		//x=x+alphaCG*p
		cblas_daxpy(N,alphaCG,p,1,x,1);

		//r=r-alphaCG*Ap
		cblas_daxpy(N,-alphaCG,Ap,1,r,1);

		//rsnew = r'*r
		rsnew = cblas_ddot(N,r,1,r,1);

		if((sqrt(rsnew)<=EPSILON)||(k == maxIterations)){
			stop = 1;
		}

		//p=r+rsnew/rsold*p
		cblas_dscal(N, rsnew/rsold, p, 1);
		cblas_daxpy(N,1.0,r,1,p,1);

		rsold = rsnew;

		k++;
	}

	//memcpy(b, x, N*sizeof(double));

        free(Ap);
        free(r);
        free(p);
	//free(x);

	fprintf(stderr, "[%s] Number of iterations %lu\n",__func__,k);

	return 1;
}


text:
/*
* -----------------------------------------------------------------
*  Binary Search Tree Library --- bst_lib.c
*  Version: 1.6180
*  Date: Mar 12, 2010
* -----------------------------------------------------------------
*  Programmer: Americo Barbosa da Cunha Junior
*              americo.cunhajr@gmail.com
* -----------------------------------------------------------------
*  Copyright (c) 2010 by Americo Barbosa da Cunha Junior
*
*  This program is free software: you can redistribute it and/or
*  modify it under the terms of the GNU General Public License as
*  published by the Free Software Foundation, either version 3 of
*  the License, or (at your option) any later version.
*
*  This program is distributed in the hope that it will be useful,
*  but WITHOUT ANY WARRANTY; without even the implied warranty of
*  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
*  GNU General Public License for more details.
*
*  A copy of the GNU General Public License is available in
*  LICENSE.txt or http://www.gnu.org/licenses/.
* -----------------------------------------------------------------
*  This is the implementation file of a library
*  to work with binary search trees.
* -----------------------------------------------------------------
*/

#include <stdlib.h>
#include <gsl/gsl_blas.h>
#include <gsl/gsl_errno.h>

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

/*
*------------------------------------------------------------
*   bst_leaf_alloc
*
*   This function allocates a struct leaf.
*
*   Input:
*   void
*
*   Output:
*   leaf - pointer to a struct leaf
*
*   last update: May 10, 2009
*------------------------------------------------------------
*/

bst_leaf *bst_leaf_alloc()
{
    bst_leaf *leaf = NULL;

    /* memory allocation for bst_leaf */
    leaf = (bst_leaf *) malloc(sizeof(bst_leaf));
    if ( leaf == NULL )
        return NULL;

    /* setting bst_leaf elements equal to NULL */
    leaf->phi  = NULL;
    leaf->Rphi = NULL;
    leaf->A    = NULL;
    leaf->L    = NULL;

    return leaf;
}
/*------------------------------------------------------------*/

/*
*------------------------------------------------------------
*   bst_node_alloc
*
*   This function allocates a struct node.
*
*   Input:
*   void
*
*   Output:
*   node - pointer to a struct node
*
*   last update: May 10, 2009
*------------------------------------------------------------
*/

bst_node* bst_node_alloc()
{
    bst_node *node = NULL;

    /* memory allocation for bst_node */
    node = (bst_node *) malloc(sizeof(bst_node));
    if ( node == NULL )
        return NULL;

    /* setting bst_node elements equal to NULL */
    node->v      = NULL;
    node->r_node = NULL;
    node->l_node = NULL;
    node->r_leaf = NULL;
    node->l_leaf = NULL;

    return node;
}
/*------------------------------------------------------------*/

/*
*------------------------------------------------------------
*   bst_leaf_free
*
*   This function frees the memory used by a struct leaf.
*
*   Input:
*   leaf - pointer to a struct leaf
*
*   Output:
*   void
*
*   last update: May 9, 2009
*------------------------------------------------------------
*/

void bst_leaf_free(void **leaf_bl)
{
    bst_leaf *leaf = NULL;

    /* checking if leaf_bl is NULL */
    if ( *leaf_bl == NULL )
        return;

    leaf = (bst_leaf *) (*leaf_bl);

    /* releasing bst_leaf elements */
    if( leaf->phi != NULL )
    {
        gsl_vector_free(leaf->phi);
        leaf->phi = NULL;
    }

    if( leaf->Rphi != NULL )
    {
        gsl_vector_free(leaf->Rphi);
        leaf->Rphi = NULL;
    }

    if( leaf->A != NULL )
    {
        gsl_matrix_free(leaf->A);
        leaf->A = NULL;
    }

    if( leaf->L != NULL )
    {
        gsl_matrix_free(leaf->L);
        leaf->L = NULL;
    }

    /* releasing allocated memory by bst_leaf */
    free(*leaf_bl);
    *leaf_bl = NULL;

    return;
}
/*------------------------------------------------------------*/

/*
*------------------------------------------------------------
*   bst_node_free
*
*   This function frees the memory used by a struct node.
*
*   Input:
*   node - pointer to a struct node
*
*   Output:
*   void
*
*   last update: Feb 19, 2010
*------------------------------------------------------------
*/

void bst_node_free(void **node_bl)
{
    bst_node *node = NULL;

    /* checking if node_bl is NULL */
    if ( *node_bl == NULL )
        return;

    node = (bst_node *) (*node_bl);

    /* releasing allocated memory by bst_node elements */
    if( node->l_leaf != NULL )
    {
        bst_leaf_free((void **)&(node->l_leaf));
        node->l_leaf = NULL;
    }

    if( node->r_leaf != NULL )
    {
        bst_leaf_free((void **)&(node->r_leaf));
        node->r_leaf = NULL;
    }

    if( node->l_node != NULL )
    {
        bst_node_free((void **)&(node->l_node));
        node->l_node = NULL;
    }

    if( node->r_node != NULL )
    {
        bst_node_free((void **)&(node->r_node));
        node->r_node = NULL;
    }

    if( node->v != NULL )
    {
        gsl_vector_free(node->v);
        node->v = NULL;
    }

    /* releasing allocated memory by bst_node */
    free(*node_bl);
    *node_bl = NULL;

    return;
}
/*------------------------------------------------------------*/

/*
*------------------------------------------------------------
*   bst_leaf_set
*
*   This function sets the bst_leaf elements equal to
*   the input elements and returns GSL_SUCCESS if
*   there is no problem during its execution.
*
*   Input:
*   phi  - composition
*   Rphi - reaction mapping
*   A    - mapping gradient matrix
*   L    - ellipsoid Cholesky matrix
*
*   Output:
*   leaf - pointer to a struct leaf
*   success or error
*
*   last update: Feb 22, 2009
*------------------------------------------------------------
*/

int bst_leaf_set(gsl_vector *phi,gsl_vector *Rphi,
                    gsl_matrix *A,gsl_matrix *L,bst_leaf *leaf)
{
    /* checking if leaf is NULL */
    if ( leaf == NULL )
        return GSL_EFAILED;

    /* setting phi equal to the input vector phi */
    if ( leaf->phi == NULL )
        leaf->phi = gsl_vector_calloc(phi->size);

    gsl_vector_memcpy(leaf->phi,phi);

    /* setting Rphi equal to the input vector Rphi */
    if ( leaf->Rphi == NULL )
        leaf->Rphi = gsl_vector_calloc(Rphi->size);

    gsl_vector_memcpy(leaf->Rphi,Rphi);

    /* setting A equal to the input matrix A */
    if ( leaf->A == NULL )
        leaf->A = gsl_matrix_calloc(A->size1,A->size2);

    gsl_matrix_memcpy(leaf->A,A);

    /* setting L equal to the input matrix L */
    if ( leaf->L == NULL )
        leaf->L = gsl_matrix_calloc(L->size1,L->size2);

    gsl_matrix_memcpy(leaf->L,L);

    return GSL_SUCCESS;
}
/*------------------------------------------------------------*/

/*
*------------------------------------------------------------
*   bst_cutplane
*
*   This function compute the vector ( v = phi_r - phi_l )
*   and the scalar ( a = (|phi_r|^2-|phi_l|^2)/2 )
*   which define a cutting plane used to search any
*   information in the binary search tree.
*
*   Input:
*   phi_l - left  composition
*   phi_r - right compositon
*
*   Output:
*   v    - cutting plane normal vector
*   a    - cutting plane scalar
*
*   last update: Feb 27, 2009
*------------------------------------------------------------
*/

double bst_cutplane(gsl_vector *phi_l,gsl_vector *phi_r,gsl_vector *v)
{
    double a1, a2;

    /* v := phi_r */
    gsl_vector_memcpy(v,phi_r);

    /* v := phi_r - phi_l */
    gsl_vector_sub(v,phi_l);

    /* a1 := |phi_r|^2 */
    gsl_blas_ddot(phi_r,phi_r,&a1);

    /* a2 := |phi_l|^2 */
    gsl_blas_ddot(phi_l,phi_l,&a2);

    return 0.5*(a1-a2);
}
/*------------------------------------------------------------*/

/*
*------------------------------------------------------------
*   bst_node_set
*
*   This function set the bst_node elements equal to
*   the input elements and define the vector and
*   scalar to do the search for a information.
*
*   Input:
*   l_leaf - left  leaf
*   r_leaf - right leaf
*
*   Output:
*   node - pointer to a struct node
*   success or error
*
*   last update: Feb 27, 2009
*------------------------------------------------------------
*/

int bst_node_set(bst_leaf *l_leaf,bst_leaf *r_leaf,bst_node *node)
{
    /* checking if node is NULL */
    if ( node == NULL )
        return GSL_EFAILED;

    /* memory allocation for v */
    if ( node->v == NULL )
    	node->v = gsl_vector_calloc(r_leaf->phi->size);

    /* setting v and a */
    node->a = bst_cutplane(l_leaf->phi,r_leaf->phi,node->v);

    /* setting leaf elements equal to the input leaves */
    node->l_leaf = l_leaf;
    node->r_leaf = r_leaf;

    /* setting node elements equal to NULL */
    node->l_node = NULL;
    node->r_node = NULL;

    return GSL_SUCCESS;
}
/*------------------------------------------------------------*/

/*
*------------------------------------------------------------
*   bst_node_add
*
*   This function adds a new node to a given side
*   of a old node and returns GSL_SUCCESS if
*   there is no problem during its execution.
*
*   Input:
*   side     - new node side
*   old_node - pointer to the old node
*   new_node - pointer to the new node
*
*   Output:
*   success or error
*
*   last update: Feb 22, 2009
*------------------------------------------------------------
*/

int bst_node_add(int side,bst_node *old_node,bst_node *new_node)
{
    if( side == RIGHT )
    {
        if( old_node->r_node != NULL )
            return GSL_EFAILED;
        else
        {
            /* adding a new right node */
            old_node->r_node = new_node;
            old_node->r_leaf = NULL;

            return GSL_SUCCESS;
        }
    }
    else if( side == LEFT )
    {
        if( old_node->l_node != NULL )
            return GSL_EFAILED;
        else
        {
            /* adding a new left node */
            old_node->l_node = new_node;
            old_node->l_leaf = NULL;

            return GSL_SUCCESS;
        }
    }
    else
        return GSL_EFAILED;
}
/*------------------------------------------------------------*/

/*
*------------------------------------------------------------
*   bst_search
*
*   This function searchs for a near composition phi0 in the
*   binary search tree.
*
*   Input:
*   root     - binary tree root
*   phi      - composition
*
*   Output:
*   end_leaf - leaf with the near composition
*   end_node - node with the near composition
*   side     - near compositon leaf side
*
*   last update: Feb 3, 2009
*------------------------------------------------------------
*/

int bst_search(bst_node *root,gsl_vector *phi,
                        bst_node **end_node,bst_leaf **end_leaf)
{
    double dot;

    if( root == NULL )
    {
        *end_leaf = NULL;
        *end_node = NULL;

        return GSL_EFAILED;
    }

    /* dot := v^T*phi */
    gsl_blas_ddot(root->v,phi,&dot);

    /* if dot > a then right node else left node */
    if( dot > root->a )
    {
        /* extern node */
        if( root->r_node == NULL )
        {
            *end_leaf = root->r_leaf;
            *end_node = root;

            return RIGHT;
        }
        else
            return bst_search(root->r_node,phi,end_node,end_leaf);
    }
    else
    {
        /* extern node */
        if( root->l_node == NULL )
        {
            *end_leaf = root->l_leaf;
            *end_node = root;

            return LEFT;
        }
        else
            return bst_search(root->l_node,phi,end_node,end_leaf);
    }
}
/*------------------------------------------------------------*/

/*
*------------------------------------------------------------
*   bst_height
*
*   This function computes the binary search tree height.
*
*   Input:
*   root - binary tree root
*
*   Output:
*   height - binary search tree height
*
*   last update: Nov 2, 2009
*------------------------------------------------------------
*/

int bst_height(bst_node *root)
{
    int l_tree_height = 0;
    int r_tree_height = 0;

    if ( root != NULL )
    {
        l_tree_height = bst_height(root->l_node);
        r_tree_height = bst_height(root->r_node);

        if ( l_tree_height > r_tree_height )
            return l_tree_height + 1;
        else
            return r_tree_height + 1;
    }
    else
        return 0;

}
/*------------------------------------------------------------*/

max_stars_repo_name: americocunhajr/CRFlowLib

max_issues_repo_name: americocunhajr/CRFlowLib

max_forks_repo_name: americocunhajr/CRFlowLib