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

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

#include "test_funcs.h"

gsl_multimin_function_fdf rosenbrock =
{&rosenbrock_f,
 &rosenbrock_df,
 &rosenbrock_fdf,
 2, 0};

gsl_multimin_function rosenbrock_fmin =
{&rosenbrock_f,
 2, 0};

void rosenbrock_initpt (gsl_vector * x)
{
  gsl_vector_set (x, 0, -1.2);
  gsl_vector_set (x, 1, 1.0);
}

double rosenbrock_f (const gsl_vector * x, void *params)
{
  double u = gsl_vector_get(x,0);
  double v = gsl_vector_get(x,1);
  double a = u - 1;
  double b = u * u - v;
  fcount++;
  return a * a + 10 * b * b;
}

void rosenbrock_df (const gsl_vector * x, void *params, gsl_vector * df)
{
  double u = gsl_vector_get(x,0);
  double v = gsl_vector_get(x,1);
  double b = u * u - v;
  gcount++;
  gsl_vector_set(df,0,2 * (u - 1) + 40 * u * b);
  gsl_vector_set(df,1,-20 * b);
}

void rosenbrock_fdf (const gsl_vector * x, void *params, double * f,
                     gsl_vector * df)
{
  double u = gsl_vector_get(x,0);
  double v = gsl_vector_get(x,1);
  double a = u - 1;
  double b = u * u - v;
  gcount++;
  *f = a * a + 10 * b * b;
  gsl_vector_set(df,0,2 * (u - 1) + 40 * u * b);
  gsl_vector_set(df,1,-20 * b);
}

gsl_multimin_function_fdf roth =
{&roth_f,
 &roth_df,
 &roth_fdf,
 2, 0};

gsl_multimin_function roth_fmin =
{&roth_f,
 2, 0};

void roth_initpt (gsl_vector * x)
{
  gsl_vector_set (x, 0, 4.5);
  gsl_vector_set (x, 1, 3.5);
}

double roth_f (const gsl_vector * x, void *params)
{
  double u = gsl_vector_get(x,0);
  double v = gsl_vector_get(x,1);
  double a = -13.0 + u + ((5.0 - v)*v - 2.0)*v;
  double b = -29.0 + u + ((v + 1.0)*v - 14.0)*v;
  fcount++;
  return a * a + b * b;
}

void roth_df (const gsl_vector * x, void *params, gsl_vector * df)
{
  double u = gsl_vector_get(x,0);
  double v = gsl_vector_get(x,1);
  double a = -13.0 + u + ((5.0 - v)*v - 2.0)*v;
  double b = -29.0 + u + ((v + 1.0)*v - 14.0)*v;
  double c = -2 + v * (10 - 3 * v);
  double d = -14 + v * (2 + 3 * v);
  gcount++;
  gsl_vector_set(df,0,2 * a + 2 * b);
  gsl_vector_set(df,1,2 * a * c + 2 * b * d);
}

void roth_fdf (const gsl_vector * x, void *params, double * f,
               gsl_vector * df)
{
  *f = roth_f (x,params);
  roth_df(x,params,df);
}

gsl_multimin_function_fdf wood =
{&wood_f,
 &wood_df,
 &wood_fdf,
 4, 0};

gsl_multimin_function wood_fmin =
{&wood_f,
 4, 0};

void
wood_initpt (gsl_vector * x)
{
  gsl_vector_set (x, 0, -3.0);
  gsl_vector_set (x, 1, -1.0);
  gsl_vector_set (x, 2, -3.0);
  gsl_vector_set (x, 3, -1.0);
}

double wood_f (const gsl_vector * x, void *params)
{
  double u1 = gsl_vector_get(x,0);
  double u2 = gsl_vector_get(x,1);
  double u3 = gsl_vector_get(x,2);
  double u4 = gsl_vector_get(x,3);

  double t1 = u1 * u1 - u2;
  double t2 = u3 * u3 - u4;
  fcount++;
  return 100 * t1 * t1 + (1 - u1) * (1 - u1)
    + 90 * t2 * t2 + (1 - u3) * (1 - u3)
    + 10.1 * ( (1 - u2) * (1 - u2) + (1 - u4) * (1 - u4) )
    + 19.8 * (1 - u2) * (1 - u4);
}

void wood_df (const gsl_vector * x, void *params, gsl_vector * df)
{
  double u1 = gsl_vector_get(x,0);
  double u2 = gsl_vector_get(x,1);
  double u3 = gsl_vector_get(x,2);
  double u4 = gsl_vector_get(x,3);

  double t1 = u1 * u1 - u2;
  double t2 = u3 * u3 - u4;
  gcount++;
  gsl_vector_set(df,0, 400 * u1 * t1 - 2 * (1 - u1) );
  gsl_vector_set(df,1, -200 * t1 - 20.2 * (1 - u2) - 19.8 * (1 - u4) );
  gsl_vector_set(df,2, 360 * u3 * t2 - 2 * (1 - u3) );
  gsl_vector_set(df,3, -180 * t2 - 20.2 * (1 - u4) - 19.8 * (1 - u2) );

}

void wood_fdf (const gsl_vector * x, void *params, double * f, gsl_vector * df)
{
  wood_df(x,params,df);
  *f=wood_f(x,params);
}

gsl_multimin_function_fdf Nrosenbrock =
{&rosenbrock_f,
 &Nrosenbrock_df,
 &Nrosenbrock_fdf,
 2, 0};

void Nrosenbrock_df (const gsl_vector * x, void *params, gsl_vector * df)
{
  gsl_multimin_function F ;
  F.f = rosenbrock_f;
  F.params = params;
  F.n = x->size;
  gsl_multimin_diff (&F, x, df);
}

void Nrosenbrock_fdf (const gsl_vector * x, void *params, double * f,
                     gsl_vector * df)
{
  *f = rosenbrock_f (x, params);
  Nrosenbrock_df (x, params, df);
}

gsl_multimin_function_fdf Nroth =
{&roth_f,
 &Nroth_df,
 &Nroth_fdf,
 2, 0};

void Nroth_df (const gsl_vector * x, void *params, gsl_vector * df)
{
  gsl_multimin_function F ;
  F.f = roth_f;
  F.params = params;
  F.n = x->size;
  gsl_multimin_diff (&F, x, df);
}

void Nroth_fdf (const gsl_vector * x, void *params, double * f,
                     gsl_vector * df)
{
  *f = roth_f (x, params);
  Nroth_df (x, params, df);
}

gsl_multimin_function_fdf Nwood =
{&wood_f,
 &Nwood_df,
 &Nwood_fdf,
 4, 0};

void Nwood_df (const gsl_vector * x, void *params, gsl_vector * df)
{
  gsl_multimin_function F ;
  F.f = wood_f;
  F.params = params;
  F.n = x->size;
  gsl_multimin_diff (&F, x, df);
}

void Nwood_fdf (const gsl_vector * x, void *params, double * f,
                     gsl_vector * df)
{
  *f = wood_f (x, params);
  Nwood_df (x, params, df);
}

max_issues_repo_name: manggoguy/parsec-modified

max_forks_repo_name: manggoguy/parsec-modified


text:
#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include <gsl/gsl_sf_legendre.h>
#include <gsl/gsl_integration.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_interp.h>
#include <gsl/gsl_spline.h>
#include <gsl/gsl_linalg.h>
#include <gsl/gsl_matrix.h>

/************************************************************************************************************************************************************************/
/***********\\Global variable\\******************************************************************************************************************************************/

size_t limit = 100;								// Limit for the convergence of the integrale
size_t size = 100;
double 	kmin,   								// minimum k in the power spectrum file (define in main function)
			 	kmax, 									// minimum k in the power spectrum file (define in main function)
			 	lk, 									    // lk = kmax-kmin lenth for the integral calculation (define in main function)
			 	q_v[3],									// unit vector q_i, q_j, q_k
			 	sig8;									// Value of sigma 8 (define in main function)

struct my_f_params {double a; double b; double c; double d;};
gsl_interp_accel *acc[25];
gsl_spline *spline[25];						//Use for interpolation

/************************************************************************************************************************************************************************/
/*******************\\ Gauss-Legendre integral quadrature\\***********************************************************************************************************/

static const double x[] = {
    1.56289844215430828714e-02,    4.68716824215916316162e-02,
    7.80685828134366366918e-02,    1.09189203580061115002e-01,
    1.40203137236113973212e-01,    1.71080080538603274883e-01,
    2.01789864095735997236e-01,    2.32302481844973969643e-01,
    2.62588120371503479163e-01,    2.92617188038471964730e-01,
    3.22360343900529151720e-01,    3.51788526372421720979e-01,
    3.80872981624629956772e-01,    4.09585291678301542532e-01,
    4.37897402172031513100e-01,    4.65781649773358042251e-01,
    4.93210789208190933576e-01,    5.20158019881763056670e-01,
    5.46597012065094167460e-01,    5.72501932621381191292e-01,
    5.97847470247178721259e-01,    6.22608860203707771585e-01,
    6.46761908514129279840e-01,    6.70283015603141015784e-01,
    6.93149199355801965946e-01,    7.15338117573056446485e-01,
    7.36828089802020705530e-01,    7.57598118519707176062e-01,
    7.77627909649495475605e-01,    7.96897892390314476375e-01,
    8.15389238339176254384e-01,    8.33083879888400823522e-01,
    8.49964527879591284320e-01,    8.66014688497164623416e-01,
    8.81218679385018415547e-01,    8.95561644970726986709e-01,
    9.09029570982529690453e-01,    9.21609298145333952679e-01,
    9.33288535043079545942e-01,    9.44055870136255977955e-01,
    9.53900782925491742847e-01,    9.62813654255815527284e-01,
    9.70785775763706331929e-01,    9.77809358486918288561e-01,
    9.83877540706057015509e-01,    9.88984395242991747997e-01,
    9.93124937037443459632e-01,    9.96295134733125149166e-01,
    9.98491950639595818382e-01,    9.99713726773441233703e-01
};

static const double A[] = {
    3.12554234538633569472e-02,    3.12248842548493577326e-02,
    3.11638356962099067834e-02,    3.10723374275665165874e-02,
    3.09504788504909882337e-02,    3.07983790311525904274e-02,
    3.06161865839804484966e-02,    3.04040795264548200160e-02,
    3.01622651051691449196e-02,    2.98909795933328309169e-02,
    2.95904880599126425122e-02,    2.92610841106382766198e-02,
    2.89030896011252031353e-02,    2.85168543223950979908e-02,
    2.81027556591011733175e-02,    2.76611982207923882944e-02,
    2.71926134465768801373e-02,    2.66974591835709626611e-02,
    2.61762192395456763420e-02,    2.56294029102081160751e-02,
    2.50575444815795897034e-02,    2.44612027079570527207e-02,
    2.38409602659682059633e-02,    2.31974231852541216230e-02,
    2.25312202563362727021e-02,    2.18430024162473863146e-02,
    2.11334421125276415432e-02,    2.04032326462094327666e-02,
    1.96530874944353058650e-02,    1.88837396133749045537e-02,
    1.80959407221281166640e-02,    1.72904605683235824399e-02,
    1.64680861761452126430e-02,    1.56296210775460027242e-02,
    1.47758845274413017686e-02,    1.39077107037187726882e-02,
    1.30259478929715422855e-02,    1.21314576629794974079e-02,
    1.12251140231859771176e-02,    1.03078025748689695861e-02,
    9.38041965369445795116e-03,    8.44387146966897140266e-03,
    7.49907325546471157895e-03,    6.54694845084532276405e-03,
    5.58842800386551515727e-03,    4.62445006342211935096e-03,
    3.65596120132637518238e-03,    2.68392537155348241939e-03,
    1.70939265351810523958e-03,    7.34634490505671730396e-04
};

#define NUM_OF_POSITIVE_ZEROS  sizeof(x) / sizeof(double)
#define NUM_OF_ZEROS           NUM_OF_POSITIVE_ZEROS+NUM_OF_POSITIVE_ZEROS

double
  gl_int(double a, double b, double (*f)(double, void *), void *prms)
{
   double integral = 0.0;
   double c = 0.5 * (b - a);
   double d = 0.5 * (b + a);
   double dum;
   const double *px = &x[NUM_OF_POSITIVE_ZEROS - 1];
   const double *pA = &A[NUM_OF_POSITIVE_ZEROS - 1];

   for (; px >= x; pA--, px--) {
      dum = c * *px;
      integral += *pA * ( (*f)(d - dum,prms) + (*f)(d + dum,prms) );
   }

   return c * integral;
}

void
  gl_int2(double a, double b, void (*f)(double, double, double*,double*), void *prms,double *int1,double *int2)
{
   double integral1 = 0.0;
   double integral2 = 0.0;
   double c = 0.5 * (b - a);
   double d = 0.5 * (b + a);
   double dum;
   const double *px = &x[NUM_OF_POSITIVE_ZEROS - 1];
   const double *pA = &A[NUM_OF_POSITIVE_ZEROS - 1];
   double outi[2], outs[2];
   double R=* (double *)prms;

   for (; px >= x; pA--, px--) {
      dum = c * *px;
      (*f)(d - dum,R,&outi[0],&outi[1]);
      (*f)(d + dum,R,&outs[0],&outs[1]);
      integral1 += *pA * (outi[0] + outs[0]);
      integral2 += *pA * (outi[1] + outs[1]);
   }

   *int1 = c * integral1;
   *int2 = c * integral2;
}

void
  gl_int5(double a, double b, void (*f)(double, double, double*,double*,double*,double*,double*), void *prms, double *int1,double *int2, double *int3, double *int4, double *int5)
{
    double integral1 = 0.0;
    double integral2 = 0.0;
    double integral3 = 0.0;
  	double integral4 = 0.0;
   	double integral5 = 0.0;
    double c = 0.5 * (b - a);
    double d = 0.5 * (b + a);
    double dum;
    double outi[5],outs[5];
    const double *px = &x[NUM_OF_POSITIVE_ZEROS - 1];
    const double *pA = &A[NUM_OF_POSITIVE_ZEROS - 1];
	double R=* (double *)prms;

	for (; px >= x; pA--, px--) {
 		dum = c * *px;
    	(*f)(d - dum,R,&outi[0],&outi[1],&outi[2],&outi[3],&outi[4]);
    	(*f)(d + dum,R,&outs[0],&outs[1],&outs[2],&outs[3],&outs[4]);
      	integral1 += *pA * (outi[0] + outs[0]);
      	integral2 += *pA * (outi[1] + outs[1]);
        integral3 += *pA * (outi[2] + outs[2]);
      	integral4 += *pA * (outi[3] + outs[3]);
		integral5 += *pA * (outi[4] + outs[4]);
	}
   	*int1 = c * integral1;
   	*int2 = c * integral2;
   	*int3 = c * integral3;
   	*int4 = c * integral4;
   	*int5 = c * integral5;
}

void
  gl_int6(double a, double b, void (*f)(double, double, double []), void *prms,double result[], int n)
{
	int i;
   double integral[n];
   for (i=0; i<n; i++)
		integral[i] = 0;
	double c = 0.5 * (b - a);
    double d = 0.5 * (b + a);
    double dum;
    const double *px = &x[NUM_OF_POSITIVE_ZEROS - 1];
	const double *pA = &A[NUM_OF_POSITIVE_ZEROS - 1];
	double outi[n], outs[n];
	double R=* (double *)prms;

    for (; px >= x; pA--, px--) {
		dum = c * *px;
		(*f)(d - dum,R,outi);
		(*f)(d + dum,R,outs);
		for (i=0; i<n; i++)
			integral[i] += *pA * (outi[i] + outs[i]);
	   }
		for (i=0; i<n; i++)
			result[i] = c*integral[i];
}

void Gauss_Legendre_Zeros_100pts( double zeros[] ) {

   const double *px = &x[NUM_OF_POSITIVE_ZEROS - 1];
   double *pz = &zeros[NUM_OF_ZEROS - 1];

   for (; px >= x; px--)  {
      *(zeros++) = - *px;
      *(pz--) = *px;
   }
}

void Gauss_Legendre_Coefs_100pts( double coefs[] ) {

   const double *pA = &A[NUM_OF_POSITIVE_ZEROS - 1];
   double *pc = &coefs[NUM_OF_ZEROS - 1];

   for (; pA >= A; pA--)  {
      *(coefs++) =  *pA;
      *(pc--) = *pA;
   }
}

/************************************************************************************************************************************************************************/
/****************************\\gsl_cquad integration \\*****************************************************************************************************************/

double int_cquad (double func(double, void*), double alpha)
{
	double result, error;
    gsl_function F;
    F.function = func;
    F.params = α

    gsl_integration_cquad_workspace *w = gsl_integration_cquad_workspace_alloc(size);
    gsl_integration_cquad(&F, kmin, kmax, 0, 1e-9, w, &result, &error, NULL);
    gsl_integration_cquad_workspace_free(w);
    return result;
}

/************************************************************************************************************************************************************************/
/****************************\\Bessel function\\************************************************************************************************************************/

double func_bessel_0(double x)
{
	return sin(x)/x;
}

double func_bessel_1(double x)
{
	return sin(x)/(x*x)-cos(x)/x;
}

double func_bessel_2 (double x)
{
	return (3./(x*x)-1)*sin(x)/x-3*cos(x)/(x*x);
}

double func_bessel_3 (double x)
{
	return (15./pow(x,3)-6./(x))*sin(x)/x-(15./(x*x)-1)*cos(x)/x;
}

/************************************************************************************************************************************************************************/
/********************************\\ Interpolation Function \\***********************************************************************************************************/

double P_L (double k)
{
    if (k>=kmin && k<=kmax)  return gsl_spline_eval (spline[1], k, acc[1]);
    else return 0;
}

double R_1(double k)
{
    if (k>=kmin && k<=kmax)  return gsl_spline_eval (spline[2], k, acc[2]);
    else return 0;
}

double R_2(double k)
{
    if (k>=kmin && k<=kmax)  return gsl_spline_eval (spline[3], k, acc[3]);
    else return 0;
}

double Q_1(double k)
{
    if (k>=kmin && k<=kmax)  return gsl_spline_eval (spline[4], k, acc[4]);
    else return 0;
}

double Q_2(double k)
{
    if (k>=kmin && k<=kmax)  return gsl_spline_eval (spline[5], k, acc[5]);
    else return 0;
}

double Q_5(double k)
{
    if (k>=kmin && k<=kmax)  return gsl_spline_eval (spline[6], k, acc[6]);
    else return 0;
}

double Q_8(double k)
{
    if (k>=kmin && k<=kmax)  return gsl_spline_eval (spline[7], k, acc[7]);
    else return 0;
}

double iXi_L(double q)
{
    if (q>0 && q<=1999)  return gsl_spline_eval (spline[8], q, acc[8]);
    else return 0;
}

double iU_1(double q)
{
    if (q>0 && q<=1999)  return gsl_spline_eval (spline[9], q, acc[9]);
    else return 0;
}

double iU_3(double q)
{
    if (q>0 && q<=1999)  return gsl_spline_eval (spline[10], q, acc[10]);
    else return 0;
}

double iU_11(double q)
{
    if (q>0 && q<=1999)  return gsl_spline_eval (spline[11], q, acc[11]);
    else return 0;
}

double iU_20(double q)
{
    if (q>0 && q<=1999)  return gsl_spline_eval (spline[12], q, acc[12]);
    else return 0;
}

double iX_11(double q)
{
    if (q>0 && q<=1999)  return gsl_spline_eval (spline[13], q, acc[13]);
    else return 0;
}

double iX_13(double q)
{
    if (q>0 && q<=1999)  return gsl_spline_eval (spline[14], q, acc[14]);
    else return 0;
}

double iX_22(double q)
{
    if (q>0 && q<=1999)  return gsl_spline_eval (spline[15], q, acc[15]);
    else return 0;
}

double iX_10_12(double q)
{
    if (q>0 && q<=1999)  return gsl_spline_eval (spline[16], q, acc[16]);
    else return 0;
}

double iY_11(double q)
{
    if (q>0 && q<=1999)  return gsl_spline_eval (spline[17], q, acc[17]);
    else return 0;
}

double iY_13(double q)
{
    if (q>0 && q<=1999)  return gsl_spline_eval (spline[18], q, acc[18]);
    else return 0;
}

double iY_22(double q)
{
    if (q>0 && q<=1999)  return gsl_spline_eval (spline[19], q, acc[19]);
    else return 0;
}

double iY_10_12(double q)
{
    if (q>0 && q<=1999)  return gsl_spline_eval (spline[20], q, acc[20]);
    else return 0;
}

double iW_V1(double q)
{
    if (q>0 && q<=1999)  return gsl_spline_eval (spline[21], q, acc[21]);
    else return 0;
}

double iW_V3(double q)
{
    if (q>0 && q<=1999)  return gsl_spline_eval (spline[22], q, acc[22]);
    else return 0;
}

double iW_T(double q)
{
    if (q>0 && q<=1999)  return gsl_spline_eval (spline[23], q, acc[23]);
    else return 0;
}

double iXi_R(double r)
{
	if (r>0 && r<=129)  return gsl_spline_eval (spline[24], r, acc[24]);
	else return 0;
}

/************************************************************************************************************************************************************************/
/*****************\\delta kronecker\\***********************************************************************************************************************************/

int delta_K (int i, int j)
{
	if (i == j) return 1;
	else return 0;
}

/************************************************************************************************************************************************************************/
/***************************************************\\Sigma 8\\*********************************************************************************************************/

double fS8 (double k, void * params)
{
  double x = k*8;
  return P_L(k)*k*k*(sin(x)-x*cos(x))*(sin(x)-x*cos(x))/pow(x,6);
}

double S8 (void)
{
    double result, error;

    gsl_function F;
    F.function = &fS8;

    gsl_integration_cquad_workspace *w = gsl_integration_cquad_workspace_alloc(size);
    gsl_integration_cquad(&F, 0, 130, 0, 1e-12, w, &result, &error, NULL);
    gsl_integration_cquad_workspace_free(w);

    return 9./(2.*M_PI*M_PI)*result;
}
/*
double P_L (double k)
{
	return Pm(k)/sig8;
}
*/
/************************************************************************************************************************************************************************/
/****************************************\\ R_n FUNCTIONS \\*********************************************************************************************************/

//	R_1_tilde function
double fRt1 (double x, void *p)
{
	double r = *(double *) p;
	return (r*r*(1.-x*x)*(1.-x*x))/(1.+r*r-2.*x*r);
}

double Rt_1 (double r)
{
	return gl_int(-1, 1, &fRt1, &r);
}

//	R_2_tilde function
double fRt2 (double x, void *p)
{
	double r = *(double *) p;
	return ((1.-x*x)*r*x*(1.-x*r))/(1.+r*r-2.*r*x);
}

double Rt_2 (double r)
{
	return gl_int(-1, 1, &fRt2, &r);
}

// R1 function
double fR_1 (double r, void *p)
{
	double k = *(double *) p;
	return P_L(k*r)*Rt_1(r);
}

double R1(double k)
{
	double result, error;
    gsl_function F;
    F.function = &fR_1;
    F.params = &k;

    gsl_integration_workspace *w = gsl_integration_workspace_alloc(200);
    gsl_integration_qag(&F, 0, kmax/k, 0, 1e-3, 200, 6, w, &result, &error);			// Kmax/k car k*r<=! kmax sinon Pk = 0
    gsl_integration_workspace_free(w);
	return (k*k*k)/(4.*M_PI*M_PI)*P_L(k)*result;
}

// R2 function
double fR_2 (double r, void *p)
{
	double k = *(double *) p;
	return P_L(k*r)*Rt_2(r);
}

double R2(double k)
{
	double result, error;
    gsl_function F;
    F.function = &fR_2;
    F.params = &k;

    gsl_integration_workspace *w = gsl_integration_workspace_alloc(200);
    gsl_integration_qag(&F, 0, kmax/k, 0, 5e-3, 200, 6, w, &result, &error);
    gsl_integration_workspace_free(w);
	return (k*k*k)/(4.*M_PI*M_PI)*P_L(k)*result;//gl_int(0, 800, &fR_2, &k);
}

/*****************Compute R1 and R2 in text files for interpolation*************************************************************************************************/

void calc_R1 ()
{
	printf("getting R function 1/2 ...\n");
	FILE* f;
	double k;
	f = fopen("data/func_R1.dat","w+");
	for (k=kmin; k<1; k+=0.001)
		fprintf(f, "%le %le\n", k, R1(k));
	fclose(f);
	f = fopen("data/func_R1.dat","a+");
	for (k=1; k<(kmax+9.9e-3); k+=0.01)
		fprintf(f, "%le %le\n", k, R1(k));
	fclose(f);
	return;
}
void calc_R2 ()
{
	printf("getting R function 2/2 ...\n");
	FILE* f;
	double k;
	f = fopen("data/func_R2.dat","w+");
	for (k=kmin; k<1; k+=0.001)
		fprintf(f, "%le %le\n", k, R2(k));
	fclose(f);
	f = fopen("data/func_R2.dat","a+");
	for (k=1; k<(kmax+9.9e-3); k+=0.01)
		fprintf(f, "%le %le\n", k, R2(k));
	fclose(f);
	return;
}

/***********************************************************************************************************************************************************************/
/******************************\\FUNCTION Q\\***********************************************************************************************************************/

// Compute Qn(r,x)
double Qt_1 (double r, double x)
{
	return (r*r*(1.-x*x)*(1.-x*x))/((1.+r*r-2.*x*r)*(1.+r*r-2.*x*r));
}

double Qt_2 (double r, double x)
{
	return ((1.-x*x)*r*x*(1.-x*r))/((1.+r*r-2.*r*x)*(1.+r*r-2.*r*x));
}

double Qt_5 (double r, double x)
{
	return (r*x*(1.-x*x))/(1.+r*r-2.*x*r);
}

double Qt_8 (double r, double x)
{
	return (r*r*(1.-x*x))	/(1.+r*r-2.*x*r);
}

/***********\\Integration to compute the inner integral of Qn(k)\\**************************************************************************************************/

double lg_Q_n( double x, void *p)
{
	struct my_f_params * params = (struct my_f_params *)p;
	double k = (params->a);
    double r = (params->b);
    double n = (params->c);
    double y_p = 1.+r*r-2.*x*r;
    if  (y_p < 0)	return 0;
    else {
		if (n == 1)	return P_L(k*sqrt(y_p))*Qt_1(r,x);
		if (n == 2)	return P_L(k*sqrt(y_p))*Qt_2(r,x);
		if (n == 5)	return P_L(k*sqrt(y_p))*Qt_5(r,x);
		if (n == 8)	return P_L(k*sqrt(y_p))*Qt_8(r,x);
		else return -1;
	}

}
/***********\\Compute Qn(k)\\***************************************************************************************************************************************/

double fQ_n (double r, void *p)
{
	struct my_f_params * params = (struct my_f_params *)p;
	double k = (params->a);
    double n = (params->b);
	struct my_f_params par = {k, r, n};
	return P_L(k*r)*gl_int(-1, 1, &lg_Q_n, &par);
}

double Q_n (int n, double k)
{
	struct my_f_params p = {k, n};
	double result, error;
    gsl_function F;
    F.function = &fQ_n;
    F.params = &p;

    gsl_integration_workspace *w = gsl_integration_workspace_alloc(200);
    gsl_integration_qag(&F, 0, kmax/k, 0, 1e-3, 200, 6, w, &result, &error);
    gsl_integration_workspace_free(w);

	return (k*k*k)/(4.*M_PI*M_PI)*result;
}

/*************Compute Qn in text file for interpolation*************************************************************************************************************/
void calc_Qn ()
{
	FILE* f;
	int n;
	double k;
	for (n=1; n<9; n++){
		printf("getting Q function %d/8...\n",n);
		switch (n){
		case 1:
			f = fopen("data/func_Q1.dat","w+");
			for (k=kmin; k<1; k+=0.001)
				fprintf(f, "%le %le\n", k, Q_n(1,k));
			fclose(f);
			f = fopen("data/func_Q1.dat","a+");
			for (k=1; k<kmax+100; k+=0.1)
				fprintf(f, "%le %le\n", k, Q_n(1,k));
			fclose(f);
			break;
		case 2:
			f = fopen("data/func_Q2.dat","w+");
			for (k=kmin; k<1; k+=0.001)
				fprintf(f, "%le %le\n", k, Q_n(2,k));
			fclose(f);
			f = fopen("data/func_Q2.dat","a+");
			for (k=1; k<kmax+100; k+=0.1)
				fprintf(f, "%le %le\n", k, Q_n(2,k));
			fclose(f);
			break;
		case 5:
			f = fopen("data/func_Q5.dat","w+");
			for (k=kmin; k<1; k+=0.001)
				fprintf(f, "%le %le\n", k, Q_n(5,k));
			fclose(f);
			f = fopen("data/func_Q5.dat","a+");
			for (k=1; k<kmax+100; k+=0.1)
				fprintf(f, "%le %le\n", k, Q_n(5,k));
			fclose(f);
			break;
		case 8:
			f = fopen("data/func_Q8.dat","w+");
			for (k=kmin; k<1; k+=0.001)
				fprintf(f, "%le %le\n", k, Q_n(8,k));
			fclose(f);
			f = fopen("data/func_Q8.dat","a+");
			for (k=1; k<kmax+100; k+=0.1)
				fprintf(f, "%le %le\n", k, Q_n(8,k));
			fclose(f);
			break;
		default:
			break;
		}
	}
	return;
}

/*******************************************************************************************************************************************************************/
/************************\\Compute Xi linear\\********************************************************************************************************************/

double fXi_L (double k, void * p)
{
	double q = *(double *) p;
    return P_L(k)*k*k*func_bessel_0(k*q);
}

double Xi_L (double q)
{
    return 1./(2.*M_PI*M_PI)*int_cquad(fXi_L, q);
}

/*******************************************************************************************************************************************************************/
/**************\\U function\\**************************************************************************************************************************************/
// U^1(q)
double fU_1 (double k, void *p)
{
	double q = *(double *) p;
	return -P_L(k)*k*func_bessel_1(k*q);
}

double U_1 (double q)
{
	return 1./(2.*M_PI*M_PI)*int_cquad(fU_1, q);
}

// U^3(q)
double fU_3 (double k, void *p)
{
	double q = *(double *) p;
	return -5./21.*R_1(k)*k*func_bessel_1(k*q);
}

double U_3 (double q)
{
	return 1./(2.*M_PI*M_PI)*int_cquad(fU_3, q);
}

// U_11(q)
double fU_11 (double k, void *p)
{
	double q = *(double *) p;
	return -6./7.*k*(R_1(k)+R_2(k))*func_bessel_1(k*q);
}

double U_11 (double q)
{
	return 1./(2.*M_PI*M_PI)*int_cquad(fU_11, q);
}

// U_20(q)
double fU_20 (double k, void *p)
{
	double q = *(double *) p;
	return -3./7.*k*Q_8(k)*func_bessel_1(k*q);
}

double U_20 (double q)
{
	return 1./(2.*M_PI*M_PI)*int_cquad(fU_20, q);
}

//Ui(q)
double get_U (double q)
{
	return iU_1(q)+iU_3(q);
}

/*******************************************************************************************************************************************************************/
/*******************\\X FUNCTIONS\\*****************************************************************************************************************************/
// X^11(q)
double fX_11 (double k, void *p)
{
	double q = *(double *) p;
	double x = k*q;
	return P_L(k)*(2./3.-2.*func_bessel_1(x)/x);
}

double X_11 (double q)
{
	return 1./(2.*M_PI*M_PI)*int_cquad(fX_11, q);
}

// X^22(q)
double fX_22 (double k, void *p)
{
		double q = *(double *) p;
		double x = k*q;
	return 9./98.*Q_1(k)*(2./3.-2.*func_bessel_1(x)/(x));
}

double X_22 (double q)
{
	return 1./(2.*M_PI*M_PI)*int_cquad(fX_22, q);
}

// X^13(q)
double fX_13 (double k, void *p)
{
	double q = *(double *) p;
	double x = k*q;
	return 5./21.*R_1(k)*(2./3.-2.*func_bessel_1(x)/x);
}

double X_13 (double q)
{
	return 1./(2.*M_PI*M_PI)*int_cquad(fX_13, q);
}

// X_10^12(q)
double fX_10_12 (double k, void *p)
{
	double q = *(double *) p;
	double x = k*q;
	return 1./14.*(2.*(R_1(k)-R_2(k))+3.*R_1(k)*sin(k*q)/(k*q)-3.*(3.*R_1(k)+4.*R_2(k)+2.*Q_5(k))*func_bessel_1(x)/x);
}

double X_10_12 (double q)
{
	return 1./(2.*M_PI*M_PI)*int_cquad(fX_10_12, q);
}

/*******************************************************************************************************************************************************************/
/*******************Y FUNCTIONS*****************************************************************************************************************************/
// Y^11(q)
double fY_11 (double k, void *p)
{
	double q = *(double *) p;
	double x = k*q;
	return P_L(k)*(-2.*sin(k*q)/(k*q)+6.*func_bessel_1(x)/x);
}

double Y_11 (double q)
{
	return 1./(2.*M_PI*M_PI)*int_cquad(fY_11, q);
}

// Y^22(q)
double fY_22 (double k, void *p)
{
	double q = *(double *) p;
	double x = k*q;
	return 9./98.*Q_1(k)*(-2.*sin(k*q)/(k*q)+6.*func_bessel_1(x)/x);
}

double Y_22 (double q)
{
	return 1./(2.*M_PI*M_PI)*int_cquad(fY_22, q);
}

// Y^13(q)
double fY_13 (double k, void *p)
{
	double q = *(double *) p;
	double x = k*q;
	return 5./21.*R_1(k)*(-2.*sin(k*q)/(k*q)+6.*func_bessel_1(x)/x);
}

double Y_13 (double q)
{
	return 1./(2.*M_PI*M_PI)*int_cquad(fY_13, q);
}

// Y^10_12(q)
double fY_10_12 (double k, void *p)
{
	double q = *(double *) p;
	double x = q*k;
	return -3./14.*(3.*R_1(k)+4.*R_2(k)+2.*Q_5(k))*(sin(k*q)/(k*q)-3.*func_bessel_1(x)/x);
}

double Y_10_12 (double q)
{
	return 1./(2.*M_PI*M_PI)*int_cquad(fY_10_12, q);
}

/************************************************************************************************************************************************************************/
/**************************\\V+S et T function\\************************************************************************************************************************/
// V^112_1(q) + S^112	(q)
double fV1_112 (double k, void*p)
{
	double q = *(double *) p;
	double x = k*q;
	return -3./7.*R_1(k)*func_bessel_1(x)/k + 3./(7.*k)*(2.*R_1(k)+4.*R_2(k)+Q_1(k)+2.*Q_2(k))*func_bessel_2(x)/x;
}								// V^112_1								+								S^112

double V1_112 (double q)
{
	return 1./(2.*M_PI*M_PI)*int_cquad(fV1_112, q);
}

// V^112_3(q)+ S^112	(q)
double fV3_112 (double k, void *p)
{
	double q = *(double *) p;
	double x = k*q;
	return -3./7.*Q_1(k)*func_bessel_1(x)/k + 3./(7.*k)*(2.*R_1(k)+4.*R_2(k)+Q_1(k)+2.*Q_2(k))*func_bessel_2(x)/x;
}								// V^112_3								+								S^112

double V3_112 (double q)
{
	return 1./(2.*M_PI*M_PI)*int_cquad(fV3_112, q);
}

// T^112(q)
double fT_112 (double k, void *p)
{
	double q = *(double *) p;
	double x = k*q;
	return -3./7.*(2.*R_1(k)+4.*R_2(k)+Q_1(k)+2.*Q_2(k))*func_bessel_3(x)/k;
}

double T_112 (double q)
{
	return 1./(2.*M_PI*M_PI)*int_cquad(fT_112, q);
}

/************************************************************************************************************************************************************************/
/***************************\\Compute U, X, Y, V and T (q) in text file for interpolation\\******************************************************************************/

void q_function ()
 {
 	FILE *fx, *fy, *fw, *fu, *fxi;
 	double q;

 	printf("getting Xi_L function...\n");
 	fxi = fopen("data/Xi_L_func.dat","w+");
	for (q=0; q<10; q+=.01)
		fprintf(fxi,"%lf %lf\n",q, Xi_L(q));
	fclose(fxi);
	fxi = fopen("data/Xi_L_func.dat","a+");
	for (q=11; q<2000; q+=1)
		fprintf(fxi,"%lf %lf\n",q, Xi_L(q));
	fclose(fxi);

 	printf("getting U function...\n");
 	fu = fopen("data/U_func.dat","w+");
	for (q=0; q<10; q+=.01)
		fprintf(fu,"%lf %lf %lf %lf %lf\n",q, U_1(q), U_3(q), U_11(q), U_20(q));
	fclose(fu);
	fu = fopen("data/U_func.dat","a+");
	for (q=11; q<2000; q+=1)
		fprintf(fu,"%lf %lf %lf %lf %lf\n",q, U_1(q), U_3(q), U_11(q), U_20(q));
	fclose(fu);

	printf("getting X function...\n");
	fx = fopen("data/X_func.dat","w+");
	for (q=0; q<10; q+=.01)
		fprintf(fx,"%lf %lf %lf %lf %lf\n",q, X_11(q), X_13(q), X_22(q), X_10_12(q));
	fclose(fx);
	fx = fopen("data/X_func.dat","a+");
	for (q=11; q<2000; q+=1)
		fprintf(fx,"%lf %lf %lf %lf %lf\n",q, X_11(q), X_13(q), X_22(q), X_10_12(q));
	fclose(fx);

	printf("getting Y function...\n");
	fy = fopen("data/Y_func.dat","w+");
	for (q=0; q<10; q+=.01)
		fprintf(fy,"%lf %lf %lf %lf %lf\n",q, Y_11(q), Y_13(q), Y_22(q), Y_10_12(q));
	fclose(fy);
	fy = fopen("data/Y_func.dat","a+");
	for (q=11; q<2000; q+=1)
		fprintf(fy,"%lf %lf %lf %lf %lf\n",q, Y_11(q), Y_13(q), Y_22(q), Y_10_12(q));
	fclose(fy);

	printf("getting W function...\n");
	fw = fopen("data/W_func.dat","w+");
	for (q=0; q<10; q+=.01)
		fprintf(fw,"%lf %lf %lf %lf\n",q, V1_112(q), V3_112(q), T_112(q));
	fclose(fw);
	fw = fopen("data/W_func.dat","a+");
	for (q=11; q<2000; q+=1)
		fprintf(fw,"%lf %lf %lf %lf\n",q, V1_112(q), V3_112(q), T_112(q));
	fclose(fw);
	return;
}

/************************************************************************************************************************************************************************/
/*********************\\Interpolation functions\\***********************************************************************************************************************/

void interpole(int n, char ficher[100],int val)
{
	double T_x[val];
	double T_y[val];
	FILE* f;
	int i = 0;
	f = fopen(ficher, "r");
	for(i=0; i < val; i++) fscanf(f, "%lf %lf\n", &T_x[i], &T_y[i]);
	acc[n] = gsl_interp_accel_alloc ();
    spline[n] = gsl_spline_alloc(gsl_interp_cspline, val);
    gsl_spline_init (spline[n], T_x, T_y, val);
    fclose(f);
}

void interp_qfunc()
{
	FILE *fu, *fx, *fy, *fw;
	int i;
	const int val = 2990;
	double X_q[val], X_x11[val], X_x13[val], X_x22[val], X_x1012[val], Y_q[val], Y_y11[val], Y_y13[val], Y_y22[val], Y_y1012[val], U_q[val], U_u1[val], U_u3[val], U_u11[val], U_u20[val], W_q[val], W_v1[val], W_v3[val], W_t[val];

	fu = fopen("data/U_func.dat", "r");
	for(i=0; i < val; i++) fscanf(fu, "%lf %lf %lf %lf %lf\n", &U_q[i], &U_u1[i], &U_u3[i], &U_u11[i], &U_u20[i]);
	acc[9] = gsl_interp_accel_alloc ();
    spline[9] = gsl_spline_alloc(gsl_interp_cspline, val);
    gsl_spline_init (spline[9], U_q, U_u1, val);
    acc[10] = gsl_interp_accel_alloc ();
    spline[10] = gsl_spline_alloc(gsl_interp_cspline, val);
    gsl_spline_init (spline[10], U_q, U_u3, val);
   	acc[11] = gsl_interp_accel_alloc ();
    spline[11] = gsl_spline_alloc(gsl_interp_cspline, val);
    gsl_spline_init (spline[11], U_q, U_u11, val);
    acc[12] = gsl_interp_accel_alloc ();
    spline[12] = gsl_spline_alloc(gsl_interp_cspline, val);
    gsl_spline_init (spline[12], U_q, U_u20, val);
    fclose(fu);

   	fx = fopen("data/X_func.dat", "r");
	for(i=0; i < val; i++) fscanf(fx, "%lf %lf %lf %lf %lf\n", &X_q[i], &X_x11[i], &X_x13[i], &X_x22[i], &X_x1012[i]);
	acc[13] = gsl_interp_accel_alloc ();
    spline[13] = gsl_spline_alloc(gsl_interp_cspline, val);
    gsl_spline_init (spline[13], X_q, X_x11, val);
    acc[14] = gsl_interp_accel_alloc ();
    spline[14] = gsl_spline_alloc(gsl_interp_cspline, val);
    gsl_spline_init (spline[14], X_q, X_x13, val);
   	acc[15] = gsl_interp_accel_alloc ();
    spline[15] = gsl_spline_alloc(gsl_interp_cspline, val);
    gsl_spline_init (spline[15], X_q, X_x22, val);
    acc[16] = gsl_interp_accel_alloc ();
    spline[16] = gsl_spline_alloc(gsl_interp_cspline, val);
    gsl_spline_init (spline[16], X_q, X_x1012, val);
    fclose(fx);

    fy = fopen("data/Y_func.dat", "r");
	for(i=0; i < val; i++) fscanf(fy, "%lf %lf %lf %lf %lf\n", &Y_q[i], &Y_y11[i], &Y_y13[i], &Y_y22[i], &Y_y1012[i]);
	acc[17] = gsl_interp_accel_alloc ();
    spline[17] = gsl_spline_alloc(gsl_interp_cspline, val);
    gsl_spline_init (spline[17], Y_q, Y_y11, val);
    acc[18] = gsl_interp_accel_alloc ();
    spline[18] = gsl_spline_alloc(gsl_interp_cspline, val);
    gsl_spline_init (spline[18], Y_q, Y_y13, val);
   	acc[19] = gsl_interp_accel_alloc ();
    spline[19] = gsl_spline_alloc(gsl_interp_cspline, val);
    gsl_spline_init (spline[19], Y_q, Y_y22, val);
    acc[20] = gsl_interp_accel_alloc ();
    spline[20] = gsl_spline_alloc(gsl_interp_cspline, val);
    gsl_spline_init (spline[20], Y_q, Y_y1012, val);
    fclose(fy);

    fw = fopen("data/W_func.dat", "r");
	for(i=0; i < val; i++) fscanf(fw, "%lf %lf %lf %lf\n", &W_q[i], &W_v1[i], &W_v3[i], &W_t[i]);
	acc[21] = gsl_interp_accel_alloc ();
    spline[21] = gsl_spline_alloc(gsl_interp_cspline, val);
    gsl_spline_init (spline[21], W_q, W_v1, val);
    acc[22] = gsl_interp_accel_alloc ();
    spline[22] = gsl_spline_alloc(gsl_interp_cspline, val);
    gsl_spline_init (spline[22], W_q, W_v3, val);
   	acc[23] = gsl_interp_accel_alloc ();
    spline[23] = gsl_spline_alloc(gsl_interp_cspline, val);
    gsl_spline_init (spline[23], W_q, W_t, val);
    fclose(fw);
}

/************************************************************************************************************************************************************************/
/*****************\\Matrix A \\ *****************************************************************************************************************************************/

				/*A_ij(q) = X(q)*delta_k(i,j) + Y(q)*q_v[i]*q_v[j]*/

double Aij_13 (int i, int j, double q)
{
	return iX_13(q)*delta_K(i,j)+iY_13(q)*q_v[i]*q_v[j];		// q_v[i] unit vector over i component
}

double Aij_22 (int i, int j, double q)
{
	return iX_22(q)*delta_K(i,j)+iY_22(q)*q_v[i]*q_v[j];
}

double Aij_11 (int i, int j, double q)
{
	return iX_11(q)*delta_K(i,j)+iY_11(q)*q_v[i]*q_v[j];
}

double Aij (int i, int j, double q)		// Compute the ij component of matrix A
{
	return Aij_11(i, j, q)+Aij_22(i, j, q)+2.*Aij_13(i, j, q);
}

double A10_ij (int i, int j, double q)
{
		return 2.*(iX_10_12(q)*delta_K(i,j)+iY_10_12(q)*q_v[i]*q_v[j]);
}

// Initialize A matrix
gsl_matrix *get_M (double q)
{
	int i, j;
	gsl_matrix *M = gsl_matrix_calloc (3,3);
	for (i=0; i<3; i++){
		for (j=0; j<3;j++)
			gsl_matrix_set(M, i, j, Aij(i,j,q));
	}
	return M;
}

/************************************************************************************************************************************************************************/
/******************\\ g_i, G_ij and Gamma_ijk functions\\**************************************************************************************************************/

double g_i (int i, int j, gsl_matrix *M_inv, double y[])
{
	return gsl_matrix_get(M_inv, i, j)*y[j];
}

double G_ij (int i, int j, gsl_matrix *M_inv, double g[])
{
	return gsl_matrix_get(M_inv, i, j)-g[i]*g[j];
}

double Gamma_ijk (int i, int j, int k, gsl_matrix *M_inv, double g[])
{
	return 	gsl_matrix_get(M_inv, i, j)*g[k]+ gsl_matrix_get(M_inv, k, i)*g[j] + gsl_matrix_get(M_inv, j, k)*g[i] - g[i]*g[j]*g[k];
}

/************************************************************************************************************************************************************************/
/******************\\Compute W^112_ijk(q)\\**************************************************************************************************************************/

double W_112 (int i, int j, int k, double q)
{
	return 	iW_V1(q)*delta_K(j,k)*q_v[i] +
					iW_V1(q)*delta_K(k,i)*q_v[j] +
					iW_V3(q)*delta_K(i,j)*q_v[k] +
					iW_T(q)*q_v[i]*q_v[j]*q_v[k];
}

/************************************************************************************************************************************************************************/
/*************\\Dot functions\\*****************************************************************************************************************************************/
//Udot_n
double fU_dot (double q)
{
	return iU_1(q) + 3.*iU_3(q);
}

//Adot_in
double fA_dot (int i, int n, double q)
{
	return Aij_11(i,n,q)+4.*Aij_13(i,n,q)+2.*Aij_22(i,n,q);
}

//Adot^10_in
double fA10_dot (int i, int n, double q)
{
	return 3.*(iX_10_12(q)*delta_K(i,n)+iY_10_12(q)*q_v[i]*q_v[n]);
}

//Wdot_ijn
double fW_dot (int i, int j, int n, double q)
{
	return 2.*W_112(i,j,n,q)+W_112(n,i,j,q)+W_112(j,n,i,q);
}

/************************************************************************************************************************************************************************/
/*********\\Double dot functions\\*************************************************************************************************************************************/

//A2dot_nm
double fA_2dot (int n, int m, double q)
{
	return Aij_11(n,m,q)+6.*Aij_13(n,m,q)+4.*Aij_22(n,m,q);
}

//A2dot^10_nm
double fA10_2dot (int n, int m, double q)
{
	return 4.*(iX_10_12(q)*delta_K(n,m)+iY_10_12(q)*q_v[n]*q_v[m]);	//A_nm = A_mn
}

// W2dot_inm
double fW_2dot (int i, int n, int m, double q)
{
	return 2.*W_112(i,n,m,q)+2.*W_112(m,i,n,q)+W_112(n,m,i,q);		//W^121_inm = W^112_min and W^211_inm = W^112_nmi
}

/************************************************************************************************************************************************************************/
/**************\\M_0 integral\\*****************************************************************************************************************************************/
 void M_0(double y, double R, double M0_fin[])
{
	int l, m,o;
    double mu, q_n, det_A, y_t[3], q[3], U[3], G[3][3], U11[3], U20[3], g[3], Gamma[3][3][3], W[3][3][3], W_fin[3][3][3], A10[3][3], Xi, s, fac, F_b[6], M0_tab[6];
    double dmu = 0.001;
    int signum = 0;

	for (l = 0; l<6; l++)				// Inizialise at each step
		M0_tab[l] = 0;

	for (mu=-1; mu <= 1; mu +=dmu){					// Decomp integrale dtheta en somme de riemann int=somme f(x_i)*delta_x
	    gsl_matrix *A = gsl_matrix_calloc (3, 3);
		gsl_matrix *A_inv = gsl_matrix_calloc (3, 3);		// Declaration of matrix A (3x3)

//		change variable y = q-r
		y_t[0] = y*sqrt(1-mu*mu);
		y_t[1] = 0;
		y_t[2] = y*mu;

// 	Compute vector q = y+r with r= (0,0,1) along the LOS
		q[0] = y_t[0];
		q[1] = y_t[1];
		q[2] = y_t[2]+R;
	    q_n = sqrt(q[0]*q[0]+q[1]*q[1]+q[2]*q[2]); 	// norm of q

// 	unit vector
	    for (l=0; l<3; l++)
	    	q_v[l] = q[l]/q_n;

/********Compute function ****************************************/
// Xi linear
	    Xi = iXi_L(q_n);

//  	U functions
	    for (l = 0; l<3 ; l++){
	   		U[l] = get_U(q_n)*q_v[l];
	   		U11[l] = iU_11(q_n)*q_v[l];
	   		U20[l] = iU_20(q_n)*q_v[l];
	    }

// 	Matix A
		gsl_permutation *p = gsl_permutation_alloc (3);
	  	A = get_M(q_n);															//Initialize matix A
	    gsl_linalg_LU_decomp(A, p, &signum);
	    det_A =	gsl_linalg_LU_det(A, signum);					// Determinant of A
		gsl_linalg_LU_invert(A, p, A_inv); 	    					// Inverse of matrix A

// 	g_i components
		for (l=0; l<3; l++){
			g[l] = 0;
	   		for (m=0; m<3; m++)
	    		g[l] += g_i(l, m, A_inv, y_t);
	    }

// 	G_ij components
	    for (l=0; l<3; l++){
	    	for (m=0; m<3; m++)
	    		G[l][m] = G_ij(l, m, A_inv, g);
	    }

// Gamma_ijk
		for (l=0; l<3; l++){
			for	(m=0; m<3; m++){
				for (o=0; o<3; o++)
					Gamma[l][m][o] = Gamma_ijk(l, m, o, A_inv, g);
			}
		}

// 	W_ijk
		for (l=0; l<3; l++){
			for	(m=0; m<3; m++){
				for (o=0; o<3; o++)
					W[l][m][o] = W_112(l, m, o, q_n);
			}
		}
		for (l=0; l<3; l++){
			for	(m=0; m<3; m++){
				for (o=0; o<3; o++)
					W_fin[l][m][o] = W[l][m][o]+ W[o][l][m]+ W[m][o][ l ];
			}
		}

// 	A^10_ij
		for (l=0; l<3; l++){
			for	(m=0; m<3; m++)
				A10[l][m] = A10_ij(l, m, q_n);
		}

//		Sum over all the bias component
		for(l=0; l<6; l++)
			F_b[l] = 0;

//		<F'>
		s=0;
		for (l=0; l<3; l++)
			s += U[l]*g[l];
		F_b[1] -= 2.*s;
		F_b[4] -= 2.*Xi*s; 	//  <F'F''>
		s=0;
		for (l=0; l<3; l++){
			for	(m=0; m<3; m++)
				s += A10[l][m]*G[l][m];
		}
		F_b[1] -= s;

//		<F'>²
		F_b[3] += Xi;
		s=0;
		for (l=0; l<3; l++)
			s += U11[l]*g[l];
		F_b[3] -= s;
		s=0;
		for (l=0; l<3; l++){
			for	(m=0; m<3; m++)
				s += U[l]*U[m]*G[l][m];
		}
		F_b[3] -= s;

//		<F''>
		s=0;
		for (l=0; l<3; l++)
			s += U20[l]*g[l];
		F_b[2] -= s;
		s=0;
		for (l=0; l<3; l++){
			for	(m=0; m<3; m++)
				s += U[l]*U[m]*G[l][m];
		}
		F_b[2] -= s;

//		<F''>²
		F_b[5] = 0.5*Xi*Xi;

// 	F0
		s=0;
		for (l=0; l<3; l++){
			for	(m=0; m<3; m++){
				for (o=0; o<3; o++)
					s += Gamma[l][m][o]*W_fin[l][m][o];
			}
		}
		F_b[0] += s/6.;
		F_b[0] += 1.;

// 	gaussian factor of M
	    s = 0;
		for (l=0; l<3; l++)
			s += g[l]*y_t[l];
		fac = exp(-0.5*s)/(sqrt(2.*M_PI)*2.*M_PI*sqrt(det_A));

//		Sum over all component
		for (l=0; l<6; l++)
			M0_tab[l] += F_b[l]*fac;

		gsl_permutation_free (p);
		gsl_matrix_free(A);
		gsl_matrix_free(A_inv);
	}
	s = y*y*dmu;
	for (l=0; l<6; l++)
			M0_fin[l] = M0_tab[l]*s;
}

/********Generate correlation function in real space Xi_R(r) ************************************************************************************************************/

void get_Xi_R ()
{
	printf("Calculating Xi_r...\n");
	int i, n =6;
	double r, M0[n];
	FILE *fi;
	fi = fopen("data/Xi_r_CLPT.dat","w+");
	for (r=1; r<=160; r+=1){
		for (i=0; i<n; i++)
			M0[i] = 0;
		gl_int6(0, 50, &M_0, &r, M0, n);
		for (i=0; i<n; i++)
			M0[i] *= 2.*M_PI;
		M0[0] -= 1;
		fprintf(fi,"%lf %lf %lf %lf %lf %lf %lf\n", r, M0[0], M0[1], M0[2], M0[3], M0[4], M0[5]);
	}
	fclose(fi);
	return ;
}
/************************************************************************************************************************************************************************/
/**************\\M_1 integral\\*****************************************************************************************************************************************/
 void M_1(double y, double R, double M1_fin[])
{
	const double rn[ 3 ] = { 0, 0, 1 }; // unit vector to project over the LOS
	int l, m,o;
    double mu, q_n, det_A, y_t[3], q[3], U[3], U_dot[3], G[3][3], U11[3], U20[3], g[3], W[3][3][3], W_dot[3][3][3], A_dot[3][3], A10_dot[3][3], Xi, s, fac, F_b[5], M1_tab[5];
	double dmu = 0.005;
    int signum = 0;					// Use for matrix inversion

	for (l = 0; l<5 ; l++) 			// Initialize at each step
		M1_tab[l] = 0;

	for (mu=-1; mu <= 1; mu +=dmu){					// Decomp integrale dtheta en somme de riemann int=somme f(x_i)*delta_x
	    gsl_matrix *A = gsl_matrix_calloc (3, 3);
		gsl_matrix *A_inv = gsl_matrix_calloc (3, 3);		// Declaration of matrix A (3x3)

//		change variable y = q-r
		y_t[0] = y*sqrt(1-mu*mu);
		y_t[1] = 0;
		y_t[2] = y*mu;

// 	Compute vector q = y+r with r= (0,0,1) along the LOS
		q[0] = y_t[0];
		q[1] = y_t[1];
		q[2] = y_t[2]+R;
	    q_n = sqrt(q[0]*q[0]+q[1]*q[1]+q[2]*q[2]); 	// norm of q

// 	unit vector
	    for (l=0; l<3; l++)
	    	q_v[l] = q[l]/q_n;

/********Compute function ****************************************/
// Xi linear
	    Xi = iXi_L(q_n);

//  	U functions
	    for (l = 0; l<3 ; l++){
	   		U_dot[l] = fU_dot(q_n)*q_v[l];
	   		U[l] = get_U(q_n)*q_v[l];
	   		U11[l] = 2.*iU_11(q_n)*q_v[l];			// Pq 2 ?????
	   		U20[l] = 2.*iU_20(q_n)*q_v[l];
	    }

// 	Matix A
		gsl_permutation *p = gsl_permutation_alloc (3);
	  	A = get_M(q_n);															//Initialize matix A
	    gsl_linalg_LU_decomp(A, p, &signum);
	    det_A =	gsl_linalg_LU_det(A, signum);					// Determinant of A
		gsl_linalg_LU_invert(A, p, A_inv); 	    					// Inverse of matrix A

// 	g_i components
		for (l=0; l<3; l++){
			g[l] = 0;
	   		for (m=0; m<3; m++)
	    		g[l] += g_i(l, m, A_inv, y_t);
	    }

// 	G_ij components
	    for (l=0; l<3; l++){
	    	for (m=0; m<3; m++)
	    		G[l][m] = G_ij(l, m, A_inv, g);
	    }

//		 Adot_in
		for (l=0; l<3; l++){
			for	(m=0; m<3; m++)
				A_dot[l][m] = fA_dot(l, m, q_n);
		}

// 	Adot^10_in
		for (l=0; l<3; l++){
			for	(m=0; m<3; m++)
				A10_dot[l][m] = fA10_dot(l, m, q_n);
		}
//		 W_ijn
		for (l=0; l<3; l++){
			for	(m=0; m<3; m++){
				for (o=0; o<3; o++)
					W[l][m][o] = W_112(l, m, o, q_n);
			}
		}
		for (l=0; l<3; l++){
			for	(m=0; m<3; m++){
				for (o=0; o<3; o++)
					W_dot[l][m][o] = 2.*W[l][m][o]+ W[o][l][m]+ W[m][o][ l ];
			}
		}

//		Sum over all the bias component
		for(l=0; l<5; l++)
			F_b[l] = 0;
// 	<F'> 	F_b[1]
		s = 0;
		for (l=0; l<3; l++)
			s += U_dot[l]*rn[l];
	    F_b[1] += 2.*s;
		s = 0;
		for (l=0; l<3; l++){
			for (o=0; o<3; o++)
				s += g[l]*A10_dot[l][o]*rn[o];
		}
	    F_b[1] -= 2.*s;
		s = 0;
		for (l=0; l<3; l++){
			for (m=0; m<3; m++){
				for (o=0; o<3; o++)
					s += G[l][m]*U[l]*A_dot[m][o]*rn[o];
			}
		}
	    F_b[1] -= 2.*s;

//		<F''>
		s = 0;
		for (l=0; l<3; l++)
			s += U20[l]*rn[l];
		F_b[2] += s;
		s = 0;
		for (l=0; l<3; l++){
			for	(m=0; m<3; m++)
				s += U[l]*g[l]*U_dot[m]*rn[m];
		}
		F_b[2] -= 2.*s;

//		<F'>²
		s=0;
		for (l=0; l<3; l++)
			s += U11[l]*rn[l];
	    F_b[3] += s;
		s=0;
		for (l=0; l<3; l++){
			for	(m=0; m<3; m++)
				s += U[l]*g[l]*U_dot[m]*rn[m];
		}
	    F_b[3] -= 2.*s;
		s=0;
		for (l=0; l<3; l++){
			for	(m=0; m<3; m++)
				s += g[l]*A_dot[l][m]*rn[m];
		}
	    F_b[3] -= Xi*s;

//		<F'F''>
		s=0;
		for (l=0; l<3; l++)
			s += U_dot[l]*rn[l];
		F_b[4] += 2.*Xi*s;

//		F0
		s=0;
		for (l=0; l<3; l++){
			for	(m=0; m<3; m++){
				for (o=0; o<3; o++)
					s += G[l][m]*W_dot[l][m][o]*rn[o];
			}
		}
		F_b[0] -= s/2.;
		s=0;
		for (l=0; l<3; l++){
			for	(m=0; m<3; m++)
				s += g[l]*A_dot[l][m]*rn[m];
		}
		F_b[0] -= s;
// 	gaussian factor of M
		s = 0;
		for (l=0; l<3; l++)
			s += g[l]*y_t[l];
		fac = exp(-0.5*s)/(sqrt(2.*M_PI)*2.*M_PI*sqrt(det_A));

//		Sum over all component
		for (l=0; l<5; l++)
			M1_tab[l] += F_b[l]*fac;

		gsl_matrix_free(A);
		gsl_matrix_free(A_inv);
		gsl_permutation_free (p);
	}
	s = y*y*dmu;
	for (l=0; l<5; l++)
			M1_fin[l] = M1_tab[l]*s;
}
/********Generate V_12(r) **********************************************************************************************************************************************/
/*To compare with Wang +13 we do not divise by 1+Xi_R*/
void get_V12 ()
{
	printf("Calculating V12...\n");
	int i, n=5;
	double r, result	[n];
	FILE *fi;
	fi = fopen("data/V_12_CLPT.dat","w+");
	for (r=1; r<=160; r+=1){
		gl_int6(0, 50, &M_1, &r, result, n);
		for (i=0; i<n; i++)
			result[i] *= 2.*M_PI;
		fprintf(fi,"%lf %lf %lf %lf %lf %lf\n", r, result[0], result[1], result[2], result[3], result[4]);
	}
	fclose(fi);
	return ;
}

/************************************************************************************************************************************************************************/
/************************************************************************************************************************************************************************/
//Calcul de M2
void M_2 (double y, double R, double M2_fin[])
{
	const double rn[ 3 ] = { 0, 0, 1 }; // unit vector to project on n and m
	int l, m,o, n;
    double mu, q_n, det_A, y_t[3], q[3], U[3], U_dot[3], G[3][3], g[3], W[3][3][3], W_2dot[3][3][3], A_dot[3][3], A_2dot[3][3], A10_2dot[3][3], Xi, s[3][3], f, fac, M2_tab[8];;
    double	F_par[4], 			//Bias factor for sigma_parallel
    				F_per[4];				//Bias factor for sigma_perpendicular
    double dmu = 0.005;
    int signum = 0; 				// use for matrix inversion

	for (l=0; l<8; l++)				//Initialize at each step
		M2_tab[l] = 0;

	for (mu=-1; mu <= 1; mu +=dmu){					// Decomp integrale dtheta en somme de riemann int=somme f(x_i)*delta_x
	    gsl_matrix *A = gsl_matrix_calloc (3, 3);
		gsl_matrix *A_inv = gsl_matrix_calloc (3, 3);		// Declaration of matrix A (3x3)

//		change variable y = q-r
		y_t[0] = y*sqrt(1-mu*mu);
		y_t[1] = 0;
		y_t[2] = y*mu;

// 	Compute vector q = y+r with r= (0,0,1) along the LOS
		q[0] = y_t[0];
		q[1] = y_t[1];
		q[2] = y_t[2]+R;
	    q_n = sqrt(q[0]*q[0]+q[1]*q[1]+q[2]*q[2]); 	// norm of q

// 	unit vector
	    for (l=0; l<3; l++)
	    	q_v[l] = q[l]/q_n;

/********Compute function ****************************************/
// Xi linear
	    Xi = iXi_L(q_n);

// 	U functions
	    for (l = 0; l<3 ; l++){
	    	U[l] = get_U(q_n)*q_v[l];
	   		U_dot[l] = fU_dot(q_n)*q_v[l];
	    }

// 	Matix A
		gsl_permutation *p = gsl_permutation_alloc (3);
	  	A = get_M(q_n);															//Initialize matix A
	    gsl_linalg_LU_decomp(A, p, &signum);
	    det_A =	gsl_linalg_LU_det(A, signum);					// Determinant of A
		gsl_linalg_LU_invert(A, p, A_inv); 	    					// Inverse of matrix A

// 	g_i components
		for (l=0; l<3; l++){
			g[l] = 0;
	   		for (m=0; m<3; m++)
	    		g[l] += g_i(l, m, A_inv, y_t);
	    }

// 	G_ij components
	    for (l=0; l<3; l++){
	    	for (m=0; m<3; m++)
	    		G[l][m] = G_ij(l, m, A_inv, g);
	    }

// 	Adot_nm
		for (l=0; l<3; l++){
			for	(m=0; m<3; m++)
				A_dot[l][m] = fA_dot(l, m, q_n);
		}

//		A2dot_nm
		for (l=0; l<3; l++){
			for	(m=0; m<3; m++)
				A_2dot[l][m] = fA_2dot(l, m, q_n);
		}

//		A2dot^10_nm
		for (l=0; l<3; l++){
			for	(m=0; m<3; m++)
				A10_2dot[l][m] = fA10_2dot(l, m, q_n);
		}
	// W_inm
		for (l=0; l<3; l++){
			for	(m=0; m<3; m++){
				for (o=0; o<3; o++)
					W[l][m][o] = W_112(l, m, o, q_n);
			}
		}
		for (l=0; l<3; l++){
			for	(m=0; m<3; m++){
				for (o=0; o<3; o++)
					W_2dot[l][m][o] = 2.*W[l][m][o]+ W[o][l][m]+ 2.*W[m][o][ l ];
			}
		}

//		Sum over all the bias component
		for (l=0; l<4; l++){			// Initialize at each step
			F_par[l]=0;
			F_per[l]=0;
		}

// 	<F'>
		for (l=0; l<3; l++){
			for (m=0; m<3; m++)
				s[l][m] = A10_2dot[l][m];
		}
		for (l=0; l<3; l++){
			for (m=0; m<3; m++){
				for (o=0; o<3; o++){
					s[m][o] -= (A_dot[l][m]*g[l]*U_dot[o]+A_dot[l][o]*g[l]*U_dot[m]);
					s[m][o] -= U[l]*g[l]*A_2dot[m][o];
				}
			}
		}
		for (l=0; l<3; l++){
			for (m=0; m<3; m++){
				F_par[1] += 2.*s[l][m]*rn[l]*rn[m];
				F_per[1] += 2.*s[l][m]*delta_K(l,m);
			}
		}
//		<F''>
		for (l=0; l<3; l++){
			for	(m=0; m<3; m++)
				s[l][m] = U_dot[l]*U_dot[m];
		}
		for (l=0; l<3; l++){
			for (m=0; m<3; m++){
				F_par[2] += 2.*s[l][m]*rn[l]*rn[m];
				F_per[2] += 2.*s[l][m]*delta_K(l,m);
			}
		}

//		<F'>²
		for (l=0; l<3; l++){
			for	(m=0; m<3; m++){
				s[l][m] = Xi*A_2dot[l][m];
				s[l][m] += 2.*U_dot[l]*U_dot[m];
			}
		}
		for (l=0; l<3; l++){
			for (m=0; m<3; m++){
				F_par[3] += s[l][m]*rn[l]*rn[m];
				F_per[3] += s[l][m]*delta_K(l,m);
			}
		}

//		F0
		for (l=0; l<3; l++){
			for (m=0; m<3; m++)
				s[l][m] = A_2dot[l][m];
		}
		for (l=0; l<3; l++){
			for (m=0; m<3; m++){
				for(n=0; n<3; n++){
					for (o=0; o<3; o++)
						s[n][o] -= A_dot[l][n]*A_dot[m][o]*G[l][m];
				}
			}
		}
		for (l=0; l<3; l++){
			for (m=0; m<3; m++){
				for (o=0; o<3; o++)
					s[m][o] -= W_2dot[l][m][o]*g[l];
			}
		}
		for (l=0; l<3; l++){
			for (m=0; m<3; m++){
				F_per[0] += s[l][m]*delta_K(l,m);
				F_par[0] += s[l][m]*rn[l]*rn[m];
			}
		}

// 	gaussian factor of M
		f = 0;
		for (l=0; l<3; l++)
			f += g[l]*y_t[l];
		fac = exp(-0.5*f)/(sqrt(2.*M_PI)*2.*M_PI*sqrt(fabs(det_A)));

//		Sum over all component
		for (l=0; l<4; l++){
			M2_tab[l] += F_par[l]*fac;
			M2_tab[l+4] += F_per[l]*fac;
		}
		gsl_permutation_free (p);
	}
	f = y*y*dmu;
	for (l=0; l<4; l++){
		M2_fin[l] = M2_tab[l]*f;
		M2_fin[l+4] = 0.5*f*(M2_tab[l+4]-M2_tab[l]);
	}
}

/********Generate Sigma_12(r) perp and par*******************************************************************************************************************************/
//To compare with Wang +13 we do not divise by 1+Xi_R
void get_Sigma ()
{
	int i, n=8;
	double r, result[n];
	FILE *fi;
	printf("Calculating Sigma...\n");
	fi = fopen("data/Sigma_12_CLPT.dat","w+");
	for (r=1; r<=160; r+=1){
		gl_int6(0, 50, &M_2, &r, result, n);
		for (i=0; i<n; i++)
			result[i] *= 2.*M_PI;
		fprintf(fi,"%lf %lf %lf %lf %lf %lf %lf %lf %lf\n", r, result[0], result[1], result[2], result[3], result[4], result[5], result[6], result[7]);
	}
	fclose(fi);
	return ;
}

/************************************************************************************************************************************************************************/
/************\\Count the number of ligne in a file\\********************************************************************************************************************/

int compte(FILE *fichier)
{
	int c;
	int nLignes = 0;
	int c2 = '\0';
	while((c=fgetc(fichier)) != EOF){	if(c=='\n')	nLignes++;	c2 = c;}
	if(c2 != '\n') nLignes++;
	return nLignes;
}

/************************************************************************************************************************************************************************/
/**********\\MAIN FUNCTION\\***************************************************************************************************************************************/
int
main (int argc, char *argv[])
{
	int i;
	double nLignes;
	FILE *ps;
// Count number of lines
	ps = fopen(argv[1], "r");
	if (ps!=0){
		nLignes = compte(ps);
		}
	else {
		printf("Error openning power spectrum file.\n");
		return 1;
	}

// Read kmin and kmax
	rewind(ps);
	kmin=1e30;
	kmax=0;
	for(i=0; i < nLignes; i++) {
		double val;
		fscanf(ps, "%lf %*f\n", &val);
		if (val<kmin) kmin=val;
		if (val>kmax) kmax=val;
	}
	fclose(ps);
	lk = kmax-kmin;

	interpole(1,argv[1], nLignes);

	calc_R1();
	calc_R2();
	calc_Qn();

	interpole(2,"data/func_R1.dat", 10901);
	interpole(3,"data/func_R2.dat", 10901);
	interpole(4,"data/func_Q1.dat", 2991);
	interpole(5,"data/func_Q2.dat", 2991);
	interpole(6,"data/func_Q5.dat", 2991);
	interpole(7,"data/func_Q8.dat", 2991);

	q_function();

	interpole(8,"data/Xi_L_func.dat", 2900);
	interp_qfunc();

	get_Xi_R();
	get_V12();
	get_Sigma();

    for (i=0; i<25; i++){
    	gsl_spline_free (spline[i]);
   		gsl_interp_accel_free (acc[i]);}

	return 0;
}

ext: c


text:
#pragma once

#include "datastructures/PointBuffer.h"
#include "math/Vector3.h"
#include "proj_api.h"

#include <containers/Range.h>

#include <string>

#include <gsl/gsl>

class Point;
struct AABB;

/// <summary>
/// Spatial reference systems that the SRSTransformHelper can transform into
/// </summary>
enum class TargetSRS
{
  /// <summary>
  /// WGS84 specified with latitude, longitude (radians) and height over
  /// ellipsoid (meters)
  /// </summary>
  WGS84,
  /// <summary>
  /// Cesium world coordinate system (WGS84 geocentric)
  /// </summary>
  CesiumWorld
};

/// <summary>
/// Utility class for transforming positions between spatial reference systems
/// </summary>
struct SRSTransformHelper
{
  virtual ~SRSTransformHelper();

  /// <summary>
  /// Transforms a range of positions into the target SRS. Transformation
  /// modifies the source points
  /// </summary>
  virtual void transformPositionsTo(TargetSRS targetSRS,
                                    gsl::span<Vector3<double>> positions) const = 0;

  /// <summary>
  /// Transforms a range of position from the given SRS back into source
  /// coordinate space
  /// </summary>
  virtual void transformPositionToSourceFrom(TargetSRS currentSRS,
                                             gsl::span<Vector3<double>> positions) const = 0;

  virtual void transformPointsTo(TargetSRS targetSRS,
                                 gsl::span<PointBuffer::PointReference> points) const = 0;

  virtual void transformPointsTo(TargetSRS targetSRS,
                                 util::Range<PointBuffer::PointIterator> points) const = 0;

  /// <summary>
  /// Transforms a range of axis-aligned bounding boxes into the target SRS.
  /// Transformation happens in-place. The axis-aligned property will be
  /// preserved after the transformation, which can result in an increase in
  /// volume of the AABBs
  /// </summary>
  virtual void transformAABBsTo(TargetSRS targetSRS, gsl::span<AABB> aabbs) const = 0;
};

/// <summary>
/// Identity transformation, i.e. a transformation that does nothing
/// </summary>
struct IdentityTransform : SRSTransformHelper
{
  virtual ~IdentityTransform();

  void transformPositionsTo(TargetSRS targetSRS,
                            gsl::span<Vector3<double>> positions) const override;
  void transformPositionToSourceFrom(TargetSRS currentSRS,
                                     gsl::span<Vector3<double>> positions) const override;
  void transformPointsTo(TargetSRS targetSRS,
                         gsl::span<PointBuffer::PointReference> points) const override;
  void transformPointsTo(TargetSRS targetSRS,
                         util::Range<PointBuffer::PointIterator> points) const override;
  void transformAABBsTo(TargetSRS targetSRS, gsl::span<AABB> aabbs) const override;
};

/// <summary>
/// Transformation based on proj4 coordinate system conversion
/// </summary>
struct Proj4Transform : SRSTransformHelper
{
  explicit Proj4Transform(const std::string& sourceTransformation);
  virtual ~Proj4Transform();

  void transformPositionsTo(TargetSRS targetSRS,
                            gsl::span<Vector3<double>> positions) const override;
  void transformPositionToSourceFrom(TargetSRS currentSRS,
                                     gsl::span<Vector3<double>> positions) const override;
  void transformPointsTo(TargetSRS targetSRS,
                         gsl::span<PointBuffer::PointReference> points) const override;
  void transformPointsTo(TargetSRS targetSRS,
                         util::Range<PointBuffer::PointIterator> points) const override;
  void transformAABBsTo(TargetSRS targetSRS, gsl::span<AABB> aabbs) const override;

private:
  projPJ getTargetTransformation(TargetSRS targetSRS) const;

  projPJ _sourceTransformation;
  projPJ _wgs84;
  projPJ _cesiumWorld;
};

/// <summary>
/// Finds the smallest point of the given points and subtracts it from all
/// points. This sets the origin of the points to the smallest point. The
/// smallest point prior to subtraction is returned
/// </summary>
Vector3<double>
setOriginToSmallestPoint(std::vector<Vector3<double>>& points);

template<typename Iter>
Vector3<double>
setOriginToSmallestPoint(Iter points_begin, Iter points_end)
{
  const auto dblMax = std::numeric_limits<double>::max();
  const auto smallestPoint =
    std::accumulate(points_begin,
                    points_end,
                    Vector3<double>{ dblMax, dblMax, dblMax },
                    [](const auto& accum, const auto& point) {
                      return Vector3<double>::minByAxis(accum, point.position());
                    });

  std::for_each(
    points_begin, points_end, [smallestPoint](auto& point) { point.position() -= smallestPoint; });

  return smallestPoint;
}

ext: h

max_stars_repo_name: igd-geo/schwarzwald

max_issues_repo_name: igd-geo/schwarzwald

max_forks_repo_name: igd-geo/schwarzwald


text:
/*!
 *  Copyright (c) 2014 by Contributors
 * \file base.h
 * \brief definitions of base types, operators, macros functions
 *
 * \author Bing Xu, Tianqi Chen
 */
#ifndef MSHADOW_BASE_H_
#define MSHADOW_BASE_H_
#ifdef _MSC_VER
#define _CRT_SECURE_NO_WARNINGS
#define _CRT_SECURE_NO_DEPRECATE
#define NOMINMAX
#endif
#include <cmath>
#include <cstdio>
#include <cfloat>
#include <climits>
#include <algorithm>
// macro defintiions
/*!
 * \brief if this macro is define to be 1,
 * mshadow should compile without any of other libs
 */
#ifndef MSHADOW_STAND_ALONE
#define MSHADOW_STAND_ALONE 0
#endif
/*! \brief whether do padding during allocation */
#ifndef MSHADOW_ALLOC_PAD
#define MSHADOW_ALLOC_PAD true
#endif
/*!
 * \brief
 *  x dimension of data must be bigger pad_size * ratio to be alloced padded memory,
 *  otherwise use tide allocation
 *  for example, if pad_ratio=2, GPU memory alignement size is 32,
 *  then we will only allocate padded memory if x dimension > 64
 *  set it to 0 then we will always allocate padded memory
 */
#ifndef MSHADOW_MIN_PAD_RATIO
  #define MSHADOW_MIN_PAD_RATIO 2
#endif

#if MSHADOW_STAND_ALONE
  #define MSHADOW_USE_CBLAS 0
  #define MSHADOW_USE_MKL   0
  #define MSHADOW_USE_CUDA  0
#endif

/*!
 * \brief force user to use GPU stream during computation
 *  error will be shot when default stream NULL is used
 */
#ifndef MSHADOW_FORCE_STREAM
#define MSHADOW_FORCE_STREAM 0
#endif

/*! \brief use CBLAS for CBLAS */
#ifndef MSHADOW_USE_CBLAS
  #define MSHADOW_USE_CBLAS 0
#endif
/*! \brief use MKL for BLAS */
#ifndef MSHADOW_USE_MKL
  #define MSHADOW_USE_MKL   1
#endif
/*!
 * \brief use CUDA support, must ensure that the cuda include path is correct,
 * or directly compile using nvcc
 */
#ifndef MSHADOW_USE_CUDA
  #define MSHADOW_USE_CUDA   1
#endif

/*!
 * \brief seems CUDAARCH is deprecated in future NVCC
 * set this to 1 if you want to use CUDA version smaller than 2.0
 */
#ifndef MSHADOW_OLD_CUDA
#define MSHADOW_OLD_CUDA 0
#endif

/*! \brief whether use SSE */
#ifndef MSHADOW_USE_SSE
  #define MSHADOW_USE_SSE 1
#endif
/*! \brief whether use NVML to get dynamic info */
#ifndef MSHADOW_USE_NVML
  #define MSHADOW_USE_NVML 0
#endif
// SSE is conflict with cudacc
#ifdef __CUDACC__
  #undef MSHADOW_USE_SSE
  #define MSHADOW_USE_SSE 0
#endif

#if MSHADOW_USE_CBLAS
extern "C" {
    #include <cblas.h>
}
#elif MSHADOW_USE_MKL
  #include <mkl.h>
  #include <mkl_cblas.h>
  #include <mkl_vsl.h>
  #include <mkl_vsl_functions.h>
#endif

#if MSHADOW_USE_CUDA
  #include <cublas.h>
  #include <curand.h>
#endif

#if MSHADOW_USE_NVML
  #include <nvml.h>
#endif
// --------------------------------
// MSHADOW_XINLINE is used for inlining template code for both CUDA and CPU code
#ifdef MSHADOW_XINLINE
  #error "MSHADOW_XINLINE must not be defined"
#endif
#ifdef _MSC_VER
#define MSHADOW_FORCE_INLINE __forceinline
#pragma warning( disable : 4068 )
#else
#define MSHADOW_FORCE_INLINE inline __attribute__((always_inline))
#endif
#ifdef __CUDACC__
  #define MSHADOW_XINLINE MSHADOW_FORCE_INLINE __device__ __host__
#else
  #define MSHADOW_XINLINE MSHADOW_FORCE_INLINE
#endif
/*! \brief cpu force inline */
#define MSHADOW_CINLINE MSHADOW_FORCE_INLINE

#if defined(__GXX_EXPERIMENTAL_CXX0X) ||\
    defined(__GXX_EXPERIMENTAL_CXX0X__) || __cplusplus >= 201103L
  #define MSHADOW_CONSTEXPR constexpr
#else
  #define MSHADOW_CONSTEXPR const
#endif

/*!
 * \brief default data type for tensor string
 *  in code release, change it to default_real_t
 *  during development, change it to empty string so that missing
 *  template arguments can be detected
 */
#ifndef MSHADOW_DEFAULT_DTYPE
#define MSHADOW_DEFAULT_DTYPE = default_real_t
//#define MSHADOW_DEFAULT_DTYPE
#endif

/*! \brief namespace for mshadow */
namespace mshadow {
/*! \brief buffer size for each random number generator */
const unsigned kRandBufferSize = 1000000;
/*! \brief pi  */
const float kPi = 3.1415926f;
/*! \brief type that will be used for index */
typedef unsigned index_t;
/*! \brief float point type that will be used in default by mshadow */
typedef float default_real_t;

/*! \brief namespace for operators */
namespace op {
// binary operator
/*! \brief mul operator */
struct mul{
  /*! \brief map a, b to result using defined operation */
  template<typename DType>
  MSHADOW_XINLINE static DType Map(DType a, DType b) {
    return a * b;
  }
};
/*! \brief plus operator */
struct plus {
  /*! \brief map a, b to result using defined operation */
  template<typename DType>
  MSHADOW_XINLINE static DType Map(DType a, DType b) {
    return a + b;
  }
};
/*! \brief minus operator */
struct minus {
  /*! \brief map a, b to result using defined operation */
  template<typename DType>
  MSHADOW_XINLINE static DType Map(DType a, DType b) {
    return a - b;
  }
};
/*! \brief divide operator */
struct div {
  /*! \brief map a, b to result using defined operation */
  template<typename DType>
  MSHADOW_XINLINE static DType Map(DType a, DType b) {
    return a / b;
  }
};
/*! \brief get rhs */
struct right {
  /*! \brief map a, b to result using defined operation */
  template<typename DType>
  MSHADOW_XINLINE static DType Map(DType a, DType b) {
    return b;
  }
};
// unary operator/ function: example
// these operators can be defined by user,
// in the same style as binary and unary operator
// to use, simply write F<op::identity>( src )
/*! \brief identity function that maps a real number to it self */
struct identity{
  /*! \brief map a to result using defined operation */
  template<typename DType>
  MSHADOW_XINLINE static DType Map(DType a) {
    return a;
  }
};
}  // namespace op
/*! \brief namespace for savers */
namespace sv {
/*! \brief save to saver: = */
struct saveto {
  /*! \brief save b to a using save method */
  template<typename DType>
  MSHADOW_XINLINE static void Save(DType &a, DType b) {
    a = b;
  }
  /*! \brief helper constant to use BLAS, alpha */
  inline static default_real_t AlphaBLAS(void) { return 1.0f; }
  /*! \brief helper constant to use BLAS, beta */
  inline static default_real_t BetaBLAS(void) { return 0.0f; }
  /*! \brief corresponding binary operator type */
  typedef op::right OPType;
};
/*! \brief save to saver: += */
struct plusto {
  /*! \brief save b to a using save method */
  template<typename DType>
  MSHADOW_XINLINE static void Save(DType &a, DType b) {
    a += b;
  }
  /*! \brief helper constant to use BLAS, alpha */
  inline static default_real_t AlphaBLAS(void) { return 1.0f; }
  /*! \brief helper constant to use BLAS, beta */
  inline static default_real_t BetaBLAS(void) { return 1.0f; }
  /*! \brief corresponding binary operator type */
  typedef op::plus OPType;
};
/*! \brief minus to saver: -= */
struct minusto {
  /*! \brief save b to a using save method */
  template<typename DType>
  MSHADOW_XINLINE static void Save(DType &a, DType b) {
    a -= b;
  }
  /*! \brief helper constant to use BLAS, alpha */
  inline static default_real_t AlphaBLAS(void) { return -1.0f; }
  /*! \brief helper constant to use BLAS, beta */
  inline static default_real_t BetaBLAS(void) { return 1.0f; }
  /*! \brief corresponding binary operator type */
  typedef op::minus OPType;
};
/*! \brief multiply to saver: *= */
struct multo {
  /*! \brief save b to a using save method */
  template<typename DType>
  MSHADOW_XINLINE static void Save(DType &a, DType b) {
    a *= b;
  }
  /*! \brief corresponding binary operator type */
  typedef op::mul OPType;
};
/*! \brief divide to saver: /= */
struct divto {
  /*! \brief save b to a using save method */
  template<typename DType>
  MSHADOW_XINLINE static void Save(DType& a, DType b) {
    a /= b;
  }
  /*! \brief corresponding binary operator type */
  typedef op::div OPType;
};
}  // namespace sv
/*! \brief namespace for potential reducer operations */
namespace red {
namespace limits {
/*!
 * \brief minimum value of certain types
 * \tparam DType data type
 */
template<typename DType>
MSHADOW_XINLINE DType MinValue(void);
/*! \brief minimum value of float */
template<>
MSHADOW_XINLINE float MinValue<float>(void) {
  return -FLT_MAX;
}
/*! \brief minimum value of double */
template<>
MSHADOW_XINLINE double MinValue<double>(void) {
  return -DBL_MAX;
}
/*! \brief minimum value of int */
template<>
MSHADOW_XINLINE int MinValue<int>(void) {
  return INT_MIN;
}
}  // namespace limits

/*! \brief sum reducer */
struct sum {
  /*! \brief do reduction into dst */
  template<typename DType>
  MSHADOW_XINLINE static void Reduce(volatile DType& dst,  volatile DType src) {
    dst += src;
  }
  /*!
   *\brief calculate gradient of redres with respect to redsrc,
   * redres: reduced result, redsrc: one of reduction element
   */
  template<typename DType>
  MSHADOW_XINLINE static DType PartialGrad(DType redres, DType redsrc) {
    return 1;
  }
  /*!
   *\brief set the initial value during reduction
   */
  template<typename DType>
  MSHADOW_XINLINE static void SetInitValue(DType &initv) {
    initv = 0;
  }
};
/*! \brief maximum reducer */
struct maximum {
  /*! \brief do reduction into dst */
  template<typename DType>
  MSHADOW_XINLINE static void Reduce(volatile DType& dst,  volatile DType src) {
    using namespace std;
    dst = max(dst, src);
  }
  /*!
   * \brief calculate gradient of redres with respect to redsrc,
   * redres: reduced result, redsrc: one of reduction element
   */
  template<typename DType>
  MSHADOW_XINLINE static DType PartialGrad(DType redres, DType redsrc) {
    return redres == redsrc ? 1: 0;
  }
  /*!
   *\brief set the initial value during reduction
   */
  template<typename DType>
  MSHADOW_XINLINE static void SetInitValue(DType &initv) {
    initv = limits::MinValue<DType>();
  }
};
}  // namespace red
}  // namespace mshadow
#endif  // MSHADOW_BASE_H_

max_stars_repo_name: Kingnus/mshadow

max_issues_repo_name: Kingnus/mshadow

max_forks_repo_name: Kingnus/mshadow


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

#include <starneig_config.h>
#include <starneig/configuration.h>
#include "cpu.h"
#include "../common/common.h"
#include "../common/sanity.h"
#include "../common/tiles.h"
#include "../common/trace.h"
#include <starpu.h>
#include <starpu_scheduler.h>
#include <hwloc.h>
#include <cblas.h>
#include <omp.h>

void starneig_hessenberg_cpu_prepare_column(
    void *buffers[], void *cl_args)
{
    // LAPACK subroutine that generates a real elementary reflector H
    extern void dlarfg_(int const *, double *, double *, int const *, double *);

    int i; // the index of the currect column inside the panel
    struct range_packing_info v_pi;
    starpu_codelet_unpack_args(cl_args, &i, &v_pi);

    int k = 0;

    double *Y = NULL; int ldY = 0;
    if (0 < i) {
        Y = (double *) STARPU_MATRIX_GET_PTR(buffers[k]);
        ldY = STARPU_MATRIX_GET_LD(buffers[k]);
        k++;
    }

    double *V = (double *) STARPU_MATRIX_GET_PTR(buffers[k]);
    int ldV = STARPU_MATRIX_GET_LD(buffers[k]);
    int m = STARPU_MATRIX_GET_NX(buffers[k]);
    int nb = STARPU_MATRIX_GET_NY(buffers[k]);
    k++;

    double *T = (double *) STARPU_MATRIX_GET_PTR(buffers[k]);
    int ldT = STARPU_MATRIX_GET_LD(buffers[k]);
    k++;

    double *P = (double *) STARPU_MATRIX_GET_PTR(buffers[k]);
    int ldP = STARPU_MATRIX_GET_LD(buffers[k]);
    k++;

    // an intemediate vector interface for the trailing matrix operation
    struct starpu_vector_interface **v_i =
        (struct starpu_vector_interface **)buffers + k;
    k += v_pi.handles;

    // current column
    double *p = P+i*ldP;

    //
    // update the current column
    //

    if (0 < i) {

        // A <- A - Y * V' (update column from the right)
        cblas_dgemv(CblasColMajor, CblasNoTrans, m, i,
            -1.0, Y, ldY, V+i-1, ldV, 1.0, p, 1);

        //
        // update column from the left
        //

        // we use the last column of T as a work space
        double *w = T+(nb-1)*ldT;

        // w <- V1' * b1 (upper part)
        cblas_dcopy(i, p, 1, w, 1);
        cblas_dtrmv(
            CblasColMajor, CblasLower, CblasTrans, CblasUnit, i, V, ldV, w, 1);

        // w <- w + V2' * b2 (lower part)
        cblas_dgemv(CblasColMajor, CblasTrans, m-i, i,
            1.0, V+i, ldV, p+i, 1, 1.0, w, 1);

        // w <- T' * w
        cblas_dtrmv(
            CblasColMajor, CblasUpper, CblasTrans, CblasNonUnit, i,
            T, ldT, w, 1);

        // b2 <- b2 - V2 * w
        cblas_dgemv(CblasColMajor, CblasNoTrans, m-i, i,
            -1.0, V+i, ldV, w, 1, 1.0, p+i, 1);

        // b1 <- b1 - V1 * w
        cblas_dtrmv(
            CblasColMajor, CblasLower, CblasNoTrans, CblasUnit, i,
            V, ldV, w, 1);
        cblas_daxpy(i, -1.0, w, 1, p, 1);
    }

    //
    // compute the current unit vector
    //

    int height = m-i;
    double tau, *v = V+i*ldV+i;
    memcpy(v, p+i, height*sizeof(double));
    dlarfg_(&height, p+i, v+1, (const int[]){1}, &tau);
    v[0] = 1.0;

    //
    // copy the current unit vector to the intemediate vector interface
    //

    starneig_join_range(&v_pi, v_i, v, 1);

    //
    // set elements below the subdiagonal to zero
    //

    for (int j = i+1; j < m; j++)
        p[j] = 0.0;

    //
    // store tau for future use
    //

    T[i*ldT+i] = tau;
}

void starneig_hessenberg_cpu_compute_column(
    void *buffers[], void *cl_args)
{
    struct packing_info A_pi;
    struct range_packing_info v_pi, y_pi;
    starpu_codelet_unpack_args(cl_args, &A_pi, &v_pi, &y_pi);

    STARNEIG_EVENT_BEGIN(&A_pi, starneig_event_red);

    int k = 0;

    // involved trailing matrix tiles
    struct starpu_matrix_interface **A_i =
        (struct starpu_matrix_interface **)buffers + k;
    k += A_pi.handles;

    // intemediate vector interface for the trailing matrix operation
    struct starpu_vector_interface **v_i =
        (struct starpu_vector_interface **)buffers + k;
    k += v_pi.handles;

    // intemediate vector interface from the trailing matrix operation
    struct starpu_vector_interface **y_i =
        (struct starpu_vector_interface **)buffers + k;
    k += y_pi.handles;

    int t_rows = (A_pi.rend-1) / A_pi.bm + 1 - (A_pi.rbegin-1) / A_pi.bm;
    int t_cols = (A_pi.cend-1) / A_pi.bn + 1 - (A_pi.cbegin-1) / A_pi.bn;

    //
    // loop oper tile rows
    //

    for (int i = 0; i < t_rows; i++) {

        double *y = (double *) STARPU_VECTOR_GET_PTR(y_i[i]);

        int rbegin = MAX(     0,  A_pi.rbegin - i * A_pi.bm);
        int rend =   MIN(A_pi.bm, A_pi.rend   - i * A_pi.bm);

        //
        // loop over tile columns
        //

        for (int j = 0; j < t_cols; j++) {

            double *A = (double *) STARPU_MATRIX_GET_PTR(A_i[j*t_rows+i]);
            int ldA = STARPU_MATRIX_GET_LD(A_i[j*t_rows+i]);

            double *v = (double *) STARPU_VECTOR_GET_PTR(v_i[j]);

            int cbegin = MAX(      0, A_pi.cbegin - j * A_pi.bn);
            int cend =   MIN(A_pi.bn, A_pi.cend   - j * A_pi.bn);

            cblas_dgemv(
                CblasColMajor, CblasNoTrans, rend-rbegin, cend-cbegin,
                1.0, A+cbegin*ldA+rbegin, ldA, v+cbegin, 1, 1.0, y+rbegin, 1);
        }
    }

    STARNEIG_EVENT_END();
}

void starneig_hessenberg_cpu_finish_column(
    void *buffers[], void *cl_args)
{
    struct range_packing_info y_pi;
    int i;
    starpu_codelet_unpack_args(cl_args, &i, &y_pi);

    int k = 0;

    double *V = (double *) STARPU_MATRIX_GET_PTR(buffers[k]);
    int m = STARPU_MATRIX_GET_NX(buffers[k]);
    int ldV = STARPU_MATRIX_GET_LD(buffers[k]);
    k++;

    double *T = (double *) STARPU_MATRIX_GET_PTR(buffers[k]);
    int ldT = STARPU_MATRIX_GET_LD(buffers[k]);
    k++;

    double *Y = (double *) STARPU_MATRIX_GET_PTR(buffers[k]);
    int ldY = STARPU_MATRIX_GET_LD(buffers[k]);
    k++;

    // intemediate vector interface from the trailing matrix operation
    struct starpu_vector_interface **y_i =
        (struct starpu_vector_interface **)buffers + k;
    k += y_pi.handles;

    double tau = T[i*ldT+i];
    double *v = V+i*ldV+i;

    //
    // finish Y update
    //

    starneig_join_range(&y_pi, y_i, Y+i*ldY, 0);

    // w <- V' * v (shared result)
    cblas_dgemv(CblasColMajor, CblasTrans, m-i, i,
        1.0, V+i, ldV, v, 1, 0.0, T+i*ldT, 1);

    // Y(:,i) <- Y(:,i) - Y * w
    cblas_dgemv(CblasColMajor, CblasNoTrans, m, i,
        -1.0, Y, ldY, T+i*ldT, 1, 1.0, Y+i*ldY, 1);

    cblas_dscal(m, tau, Y+i*ldY, 1);

    //
    // update T
    //

    // w <- tau * w
    cblas_dscal(i, -tau, T+i*ldT, 1);

    // T(0:i,i) = T * w
    cblas_dtrmv(
        CblasColMajor, CblasUpper, CblasNoTrans, CblasNonUnit, i,
        T, ldT, T+i*ldT, 1);

    T[i*ldT+i] = tau;
}

void starneig_hessenberg_cpu_update_trail_right(
    void *buffers[], void *cl_args)
{
    struct packing_info A_pi;
    int nb, roffset, coffset;
    starpu_codelet_unpack_args(cl_args, &A_pi, &nb, &roffset, &coffset);

    STARNEIG_EVENT_BEGIN(&A_pi, starneig_event_blue);

    int m = A_pi.rend - A_pi.rbegin;
    int n = A_pi.cend - A_pi.cbegin;

    double *V = (double *) STARPU_MATRIX_GET_PTR(buffers[0]);
    int ldV = STARPU_MATRIX_GET_LD(buffers[0]);

    double *Y = (double *) STARPU_MATRIX_GET_PTR(buffers[1]);
    int ldY = STARPU_MATRIX_GET_LD(buffers[1]);

    double *A = (double *) STARPU_MATRIX_GET_PTR(buffers[2]);
    int ldA = STARPU_MATRIX_GET_LD(buffers[2]);

    struct starpu_matrix_interface **A_i =
        (struct starpu_matrix_interface **)buffers + 3;

    // join tiles
    starneig_join_window(&A_pi, ldA, A_i, A, 0);

    // A <- Y V^T
    cblas_dgemm(CblasColMajor, CblasNoTrans, CblasTrans,
        m, n, nb, -1.0, Y+roffset, ldY, V+coffset+nb-1, ldV, 1.0, A, ldA);

    // split tiles
    starneig_join_window(&A_pi, ldA, A_i, A, 1);

    STARNEIG_EVENT_END();
}

void starneig_hessenberg_cpu_update_left_a(
    void *buffers[], void *cl_args)
{
    struct packing_info A_pi, W_pi;
    int nb, offset;
    starpu_codelet_unpack_args(cl_args, &A_pi, &W_pi, &nb, &offset);

    STARNEIG_EVENT_BEGIN(&A_pi, starneig_event_green);

    int m = A_pi.rend - A_pi.rbegin;
    int n = A_pi.cend - A_pi.cbegin;

    int k = 0;

    double *V = (double *) STARPU_MATRIX_GET_PTR(buffers[k]);
    int ldV = STARPU_MATRIX_GET_LD(buffers[k]);
    k++;

    double *T = (double *) STARPU_MATRIX_GET_PTR(buffers[k]);
    int ldT = STARPU_MATRIX_GET_LD(buffers[k]);
    k++;

    double *A = (double *) STARPU_MATRIX_GET_PTR(buffers[k]);
    int ldA = STARPU_MATRIX_GET_LD(buffers[k]);
    k++;

    double *W = (double *) STARPU_MATRIX_GET_PTR(buffers[k]);
    int ldW = STARPU_MATRIX_GET_LD(buffers[k]);
    k++;

    double *P = (double *) STARPU_MATRIX_GET_PTR(buffers[k]);
    int ldP = STARPU_MATRIX_GET_LD(buffers[k]);
    k++;

    struct starpu_matrix_interface **A_i =
        (struct starpu_matrix_interface **)buffers + k;
    k += A_pi.handles;

    struct starpu_matrix_interface **W_i =
        (struct starpu_matrix_interface **)buffers + k;
    k += W_pi.handles;

    // join A tiles
    starneig_join_window(&A_pi, ldA, A_i, A, 0);

    // join W tiles
    starneig_join_window(&W_pi, ldW, W_i, W, 0);

    // P <- A^T * V
    cblas_dgemm(
        CblasColMajor, CblasTrans, CblasNoTrans, n, nb, m,
        1.0, A, ldA, V+offset, ldV, 0.0, P, ldP);

    // P <- P * T
    cblas_dtrmm(
        CblasColMajor, CblasRight, CblasUpper, CblasNoTrans,
        CblasNonUnit, n, nb, 1.0, T, ldT, P, ldP);

    // W <- W + P
    for (int j = 0; j < nb; j++)
        cblas_daxpy(n, 1.0, P+j*ldP, 1, W+j*ldW, 1);

    // split W tiles
    starneig_join_window(&W_pi, ldW, W_i, W, 1);

    STARNEIG_EVENT_END();
}

void starneig_hessenberg_cpu_update_left_b(
    void *buffers[], void *cl_args)
{
    struct packing_info A_pi, W_pi;
    int nb, offset;
    starpu_codelet_unpack_args(cl_args, &A_pi, &W_pi, &nb, &offset);

    STARNEIG_EVENT_BEGIN(&packing_info, starneig_event_green);

    int m = A_pi.rend - A_pi.rbegin;
    int n = A_pi.cend - A_pi.cbegin;

    int k = 0;

    double *V = (double *) STARPU_MATRIX_GET_PTR(buffers[k]);
    int ldV = STARPU_MATRIX_GET_LD(buffers[k]);
    k++;

    double *W = (double *) STARPU_MATRIX_GET_PTR(buffers[k]);
    int ldW = STARPU_MATRIX_GET_LD(buffers[k]);
    k++;

    double *A = (double *) STARPU_MATRIX_GET_PTR(buffers[k]);
    int ldA = STARPU_MATRIX_GET_LD(buffers[k]);
    k++;

    struct starpu_matrix_interface **W_i =
        (struct starpu_matrix_interface **)buffers + k;
    k += W_pi.handles;

    struct starpu_matrix_interface **A_i =
        (struct starpu_matrix_interface **)buffers + k;
    k += A_pi.handles;

    // join A tiles
    starneig_join_window(&A_pi, ldA, A_i, A, 0);

    // join W tiles
    starneig_join_window(&W_pi, ldW, W_i, W, 0);

    //  A <- A - V * W^T
    cblas_dgemm(
        CblasColMajor, CblasNoTrans, CblasTrans, m, n,
        nb, -1.0, V+offset, ldV, W, ldW, 1.0, A, ldA);

    // split A tiles
    starneig_join_window(&A_pi, ldA, A_i, A, 1);

    STARNEIG_EVENT_END();
}

void starneig_hessenberg_cpu_update_right_a(
    void *buffers[], void *cl_args)
{
    struct packing_info A_pi, W_pi;
    int nb, offset;
    starpu_codelet_unpack_args(cl_args, &A_pi, &W_pi, &nb, &offset);

    STARNEIG_EVENT_BEGIN(&packing_info, starneig_event_blue);

    int m = A_pi.rend - A_pi.rbegin;
    int n = A_pi.cend - A_pi.cbegin;

    int k = 0;

    double *V = (double *) STARPU_MATRIX_GET_PTR(buffers[k]);
    int ldV = STARPU_MATRIX_GET_LD(buffers[k]);
    k++;

    double *T = (double *) STARPU_MATRIX_GET_PTR(buffers[k]);
    int ldT = STARPU_MATRIX_GET_LD(buffers[k]);
    k++;

    double *A = (double *) STARPU_MATRIX_GET_PTR(buffers[k]);
    int ldA = STARPU_MATRIX_GET_LD(buffers[k]);
    k++;

    double *W = (double *) STARPU_MATRIX_GET_PTR(buffers[k]);
    int ldW = STARPU_MATRIX_GET_LD(buffers[k]);
    k++;

    double *P = (double *) STARPU_MATRIX_GET_PTR(buffers[k]);
    int ldP = STARPU_MATRIX_GET_LD(buffers[k]);
    k++;

    struct starpu_matrix_interface **A_i =
        (struct starpu_matrix_interface **)buffers + k;
    k += A_pi.handles;

    struct starpu_matrix_interface **W_i =
        (struct starpu_matrix_interface **)buffers + k;
    k += W_pi.handles;

    // join A tiles
    starneig_join_window(&A_pi, ldA, A_i, A, 0);

    // join W tiles
    starneig_join_window(&W_pi, ldW, W_i, W, 0);

    // P <- A * V
    cblas_dgemm(
        CblasColMajor, CblasNoTrans, CblasNoTrans, m, nb, n,
        1.0, A, ldA, V+offset, ldV, 0.0, P, ldP);

    // P <- P * T
    cblas_dtrmm(
        CblasColMajor, CblasRight, CblasUpper, CblasNoTrans,
        CblasNonUnit, m, nb, 1.0, T, ldT, P, ldP);

    // W <- W + P
    for (int j = 0; j < nb; j++)
        cblas_daxpy(m, 1.0, P+j*ldP, 1, W+j*ldW, 1);

    // split W tiles
    starneig_join_window(&W_pi, ldW, W_i, W, 1);

    STARNEIG_EVENT_END();
}

void starneig_hessenberg_cpu_update_right_b(
    void *buffers[], void *cl_args)
{
    struct packing_info A_pi, W_pi;
    int nb, offset;
    starpu_codelet_unpack_args(cl_args, &A_pi, &W_pi, &nb, &offset);

    STARNEIG_EVENT_BEGIN(&packing_info, starneig_event_blue);

    int m = A_pi.rend - A_pi.rbegin;
    int n = A_pi.cend - A_pi.cbegin;

    int k = 0;

    double *V = (double *) STARPU_MATRIX_GET_PTR(buffers[k]);
    int ldV = STARPU_MATRIX_GET_LD(buffers[k]);
    k++;

    double *W = (double *) STARPU_MATRIX_GET_PTR(buffers[k]);
    int ldW = STARPU_MATRIX_GET_LD(buffers[k]);
    k++;

    double *A = (double *) STARPU_MATRIX_GET_PTR(buffers[k]);
    int ldA = STARPU_MATRIX_GET_LD(buffers[k]);
    k++;

    struct starpu_matrix_interface **W_i =
        (struct starpu_matrix_interface **)buffers + k;
    k += W_pi.handles;

    struct starpu_matrix_interface **A_i =
        (struct starpu_matrix_interface **)buffers + k;
    k += A_pi.handles;

    // join A tiles
    starneig_join_window(&A_pi, ldA, A_i, A, 0);

    // join W tiles
    starneig_join_window(&W_pi, ldW, W_i, W, 0);

    //  A <- A - W * V
    cblas_dgemm(
        CblasColMajor, CblasNoTrans, CblasTrans, m, n,
        nb, -1.0, W, ldW, V+offset, ldV, 1.0, A, ldA);

    // split A tiles
    starneig_join_window(&A_pi, ldA, A_i, A, 1);

    STARNEIG_EVENT_END();
}

max_stars_repo_name: NLAFET/StarNEig

max_issues_repo_name: NLAFET/StarNEig

max_forks_repo_name: NLAFET/StarNEig


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

#include <ndlinear/gsl_multifit_ndlinear.h>

#include <gsl/gsl_math.h>
#include <gsl/gsl_sf_laguerre.h>
#include <gsl/gsl_sf_legendre.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_multifit.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>
#include <gsl/gsl_bspline.h>
#include <gsl/gsl_statistics.h>

/* dimension of fit */
#define N_DIM            3

/* number of basis functions for each variable */
#define N_SUM_R          10
#define N_SUM_THETA      10
#define N_SUM_PHI        9

#define R_MAX            3.0

double
psi_real_exact(int k, int l, int m, double r, double theta, double phi)
{
  double R, T, P;

  R = pow(r, (double) l) *
      exp(-r*r) *
      gsl_sf_laguerre_n(k, l + 0.5, 2 * r * r);
  T = gsl_sf_legendre_sphPlm(l, m, cos(theta));
  P = cos(m * phi);

  return (R * T * P);
}

/* basis functions for each variable */

int
basis_r(double r, double y[], void *params)
{
  gsl_bspline_workspace *bw = params;
  gsl_vector_view v = gsl_vector_view_array(y, N_SUM_R);
  int s;

  /* use B-splines for r dependence */
  s = gsl_bspline_eval(r, &v.vector, bw);

  return s;
}

int
basis_theta(double theta, double y[], void *params)
{
  size_t i;

  /* use Legendre polynomials for theta dependence */
  for (i = 0; i < N_SUM_THETA; ++i)
    y[i] = gsl_sf_legendre_Pl(i, cos(theta));

  return GSL_SUCCESS;
}

int
basis_phi(double phi, double y[], void *params)
{
  size_t i;

  /* use standard Fourier basis (sin/cos) for phi dependence */
  for (i = 0; i < N_SUM_PHI; ++i)
    {
      if ((i % 2) == 0)
        y[i] = cos((double)(i/2) * phi);
      else
        y[i] = sin((double)((i+1)/2) * phi);
    }

  return GSL_SUCCESS;
}

int
main(int argc, char *argv[])
{
  const size_t ndim = N_DIM; /* dimension of fit */
  const size_t ndata = 3000; /* number of data points to fit */
  size_t N[N_DIM];           /* upper bounds on model sums */
  int (*u[N_DIM])(double x, double y[], void *params);
  size_t i;                  /* looping */
  int k, l, m;               /* quantum numbers */
  gsl_rng *rng_p;
  gsl_bspline_workspace *bspline_p;
  gsl_multifit_linear_workspace *multifit_p;
  gsl_multifit_ndlinear_workspace *ndlinear_p;
  gsl_vector *data;          /* psi data */
  gsl_matrix *vars;          /* parameters corresponding to psi data */
  gsl_matrix *X;             /* matrix for least squares fit */
  gsl_vector *coeffs;        /* fit coefficients */
  gsl_matrix *cov;           /* covariance matrix */
  double chisq;              /* chi^2 */
  double Rsq;                /* R^2 */
  size_t ncoeffs;            /* total number of fit coefficients */

  gsl_rng_env_setup();

  k = 5;
  l = 4;
  m = 2;

  N[0] = N_SUM_R;
  N[1] = N_SUM_THETA;
  N[2] = N_SUM_PHI;

  u[0] = &basis_r;
  u[1] = &basis_theta;
  u[2] = &basis_phi;

  rng_p = gsl_rng_alloc(gsl_rng_default);
  bspline_p = gsl_bspline_alloc(4, N_SUM_R - 2);
  ndlinear_p = gsl_multifit_ndlinear_alloc(ndim, N, u, bspline_p);

  ncoeffs = gsl_multifit_ndlinear_ncoeffs(ndlinear_p);

  multifit_p = gsl_multifit_linear_alloc(ndata, ncoeffs);
  data = gsl_vector_alloc(ndata);
  vars = gsl_matrix_alloc(ndata, ndim);
  X = gsl_matrix_alloc(ndata, ncoeffs);
  coeffs = gsl_vector_alloc(ncoeffs);
  cov = gsl_matrix_alloc(ncoeffs, ncoeffs);

  gsl_bspline_knots_uniform(0.0, R_MAX, bspline_p);

  /* this is the data to be fitted */

  for (i = 0; i < ndata; ++i)
    {
      double r = gsl_rng_uniform(rng_p) * R_MAX;
      double theta = gsl_rng_uniform(rng_p) * M_PI;
      double phi = gsl_rng_uniform(rng_p) * 2.0 * M_PI;
      double psi = psi_real_exact(k, l, m, r, theta, phi);
      double dpsi = gsl_ran_gaussian(rng_p, 0.05 * psi);

      /* keep track of (r, theta, phi) points */
      gsl_matrix_set(vars, i, 0, r);
      gsl_matrix_set(vars, i, 1, theta);
      gsl_matrix_set(vars, i, 2, phi);

      /* fill in RHS data vector */
      gsl_vector_set(data, i, psi + dpsi);
    }

  /* construct the design matrix X */
  gsl_multifit_ndlinear_design(vars, X, ndlinear_p);

  /* now do the actual least squares fit */
  gsl_multifit_linear(X, data, coeffs, cov, &chisq, multifit_p);

  /* compute R^2 */
  Rsq = 1.0 - chisq / gsl_stats_tss(data->data, 1, data->size);

  fprintf(stderr, "chisq = %e, Rsq = %f\n", chisq, Rsq);

  /* now print out the model and the exact solution and compute rms error */
  {
    double eps_rms = 0.0;
    double volume = 0.0;
    double r, theta, phi;
    double dr = 0.05;
    double dtheta = 5.0 * M_PI / 180.0;
    double dphi = 5.0 * M_PI / 180.0;
    double x[N_DIM];
    gsl_vector_view xv = gsl_vector_view_array(x, N_DIM);
    double psi;
    double psi_model;
    double err;

    for (r = 0.01; r < R_MAX; r += dr)
      {
        for (theta = 0.0; theta < M_PI; theta += dtheta)
          {
            for (phi = 0.0; phi < 2.0 * M_PI; phi += dphi)
              {
                double dV = r * r * sin(theta) * dr * dtheta * dphi;

                x[0] = r;
                x[1] = theta;
                x[2] = phi;

                /* compute model value for this (r, theta, phi) */
                psi_model = gsl_multifit_ndlinear_calc(&xv.vector,
                                                       coeffs,
                                                       ndlinear_p);

                /* compute exact value for this (r, theta, phi) */
                psi = psi_real_exact(k, l, m, r, theta, phi);

                err = psi_model - psi;
                eps_rms += err * err * dV;
                volume += dV;

                if (phi == 0.0)
                  printf("%e %e %e %e\n", r, theta, psi, psi_model);
              }
          }
        printf("\n");
      }

    eps_rms /= volume;
    eps_rms = sqrt(eps_rms);
    fprintf(stderr, "rms error over all parameter space = %e\n", eps_rms);
  }

  gsl_rng_free(rng_p);
  gsl_bspline_free(bspline_p);
  gsl_multifit_ndlinear_free(ndlinear_p);
  gsl_multifit_linear_free(multifit_p);
  gsl_vector_free(data);
  gsl_matrix_free(vars);
  gsl_matrix_free(X);
  gsl_vector_free(coeffs);
  gsl_matrix_free(cov);

  return 0;
} /* main() */

max_stars_repo_name: rennhak/Keyposes

max_issues_repo_name: rennhak/Keyposes

max_forks_repo_name: rennhak/Keyposes


text:
#include <petsc.h>
#include <bpetsc_impl.h>
#include <assert.h>

#define LEFT_SUBSPACE Full
  #define RIGHT_SUBSPACE Full
    #include <buildmat_template.h>
  #undef RIGHT_SUBSPACE
#undef LEFT_SUBSPACE

int main(int argc, char const *argv[]) {

  PetscErrorCode ierr;
  PetscScalar* x_array;
  PetscScalar* b_array;
  shell_context ctx;
  data_Full subspace_data;

  /* test identity */
  ctx.nmasks = 1;
  ierr = PetscMalloc1(1, &(ctx.masks));CHKERRQ(ierr);
  ierr = PetscMalloc1(2, &(ctx.mask_offsets));CHKERRQ(ierr);
  ierr = PetscMalloc1(1, &(ctx.signs));CHKERRQ(ierr);
  ierr = PetscMalloc1(1, &(ctx.real_coeffs));CHKERRQ(ierr);

  ctx.masks[0] = 0;
  ctx.mask_offsets[0] = 0;
  ctx.mask_offsets[1] = 1;
  ctx.signs[0] = 0;
  ctx.real_coeffs[0] = 1.;

  subspace_data.L = -1;
  ctx.left_subspace_data = &subspace_data;
  ctx.right_subspace_data = &subspace_data;

  #define DIM 100

  ierr = PetscMalloc1(DIM, &x_array);CHKERRQ(ierr);
  ierr = PetscMalloc1(DIM, &b_array);CHKERRQ(ierr);

  for (PetscInt i = 0; i < DIM; ++i) {
    x_array[i] = i;
    b_array[i] = 0;
  }

  MatMult_CPU_kernel_Full_Full(x_array, b_array, &ctx, 0, DIM, 0, DIM);

  for (PetscInt i = 0; i < DIM; ++i) {
    assert(b_array[i] == i);
  }

  ierr = PetscFree(x_array);CHKERRQ(ierr);
  ierr = PetscFree(b_array);CHKERRQ(ierr);

  #undef DIM

  ierr = PetscFree(ctx.masks);CHKERRQ(ierr);
  ierr = PetscFree(ctx.mask_offsets);CHKERRQ(ierr);
  ierr = PetscFree(ctx.signs);CHKERRQ(ierr);
  ierr = PetscFree(ctx.real_coeffs);CHKERRQ(ierr);

  return 0;
}

max_stars_repo_name: GregDMeyer/dynamite

max_issues_repo_name: GregDMeyer/dynamite

max_forks_repo_name: GregDMeyer/dynamite