text:
/*****************************************************************************
*
* Rokko: Integrated Interface for libraries of eigenvalue decomposition
*
* Copyright (C) 2012-2020 by Rokko Developers https://github.com/t-sakashita/rokko
*
* Distributed under the Boost Software License, Version 1.0. (See accompanying
* file LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
*
*****************************************************************************/

#ifndef ROKKO_CSCALAPACK_H
#define ROKKO_CSCALAPACK_H

#include <rokko/config.h>
#include <lapacke.h>

#ifdef __cplusplus
extern "C" {
#endif

int cscalapack_descinit(int* desc, int m, int n, int mb, int nb, int irsrc, int icsrc,
                        int ictxt, int lld);

int cscalapack_indxg2p(int indxglob, int nb, int iproc, int isrcproc, int nprocs);

int cscalapack_numroc(int n, int nb, int iproc, int isrcproc, int nprocs);

double cscalapack_pdelget(char scope, char top, const double* A, int ia, int ja, const int* descA);

void cscalapack_pdelset(double* A, int ia, int ja, const int* descA, double alpha);

double cscalapack_pdlamch(int icnt, char cmch);

double cscalapack_pdlange(char norm, int m, int n, const double* A, const int* descA);

double cscalapack_pdlange_work(char norm, int m, int n, const double* A, const int* descA,
                               double* work);

void cscalapack_pdlaprnt(int m, int n, const double* A, int ia, int ja, const int* descA,
                         int irprnt, int icprnt, const char* cmatnm, int nout, double* work);

#define CSCALAPACK_PSYEV_DECL(NAMES, NAMEL, TYPE, TYPE_REAL) \
int cscalapack_ ## NAMES ## _work (char jobz, char uplo, int n, \
                                   TYPE* A, int ia, int ja, const int* descA, \
                                   TYPE_REAL* w, TYPE* Z, int iz, int jz, const int* descZ, \
                                   TYPE* work, int lwork); \
int cscalapack_ ## NAMES (char jobz, char uplo, int n, \
                          TYPE* A, int ia, int ja, const int* descA, \
                          TYPE_REAL* w, TYPE* Z, int iz, int jz, const int* descZ);

CSCALAPACK_PSYEV_DECL(pssyev, PSSYEV, float, float);
CSCALAPACK_PSYEV_DECL(pdsyev, PDSYEV, double, double);
CSCALAPACK_PSYEV_DECL(pcheev, PCHEEV, lapack_complex_float, float);
CSCALAPACK_PSYEV_DECL(pzheev, PZHEEV, lapack_complex_double, double);

#undef CSCALAPACK_PSYEV_DECL

#define CSCALAPACK_PSYEVD_DECL(NAMES, NAMEL, TYPE, TYPE_REAL) \
int cscalapack_ ## NAMES ## _work (char jobz, char uplo, int n, \
                                   TYPE* A, int ia, int ja, const int* descA, \
                                   TYPE_REAL* w, TYPE* Z, int iz, int jz, const int* descZ, \
                                   TYPE* work, int lwork, int* iwork, int liwork); \
int cscalapack_ ## NAMES (char jobz, char uplo, int n, \
                          TYPE* A, int ia, int ja, const int* descA, \
                          TYPE_REAL* w, TYPE* Z, int iz, int jz, const int* descZ);

CSCALAPACK_PSYEVD_DECL(pssyevd, PSSYEVD, float, float);
CSCALAPACK_PSYEVD_DECL(pdsyevd, PDSYEVD, double, double);
CSCALAPACK_PSYEVD_DECL(pcheevd, PCHEEVD, lapack_complex_float, float);
CSCALAPACK_PSYEVD_DECL(pzheevd, PZHEEVD, lapack_complex_double, double);

#undef CSCALAPACK_PSYEVD_DECL

#define CSCALAPACK_PSYEVX_DECL(NAMES, NAMEL, TYPE, TYPE_REAL) \
int cscalapack_ ## NAMES ## _work(char jobz, char range, char uplo, int n, \
                                  TYPE* A, int ia, int ja, const int* descA, \
                                  TYPE_REAL vl, TYPE_REAL vu, int il, int iu, \
                                  TYPE_REAL abstol, int* m, int* nZ, TYPE_REAL* w, TYPE_REAL orfac, \
                                  TYPE* Z, int iz, int jz, const int* descZ, \
                                  TYPE* work, int lwork, int* iwork, int liwork, \
                                  int* ifail, int* iclustr, TYPE_REAL* gap); \
int cscalapack_ ## NAMES (char jobz, char range, char uplo, int n, \
                          TYPE* A, int ia, int ja, const int* descA, \
                          TYPE_REAL vl, TYPE_REAL vu, int il, int iu, \
                          TYPE_REAL abstol, int* m, int* nZ, TYPE_REAL* w, TYPE_REAL orfac, \
                          TYPE* Z, int iz, int jz, const int* descZ, \
                          int* ifail, int* iclustr, TYPE_REAL* gap);

CSCALAPACK_PSYEVX_DECL(pssyevx, PSSYEVX, float, float);
CSCALAPACK_PSYEVX_DECL(pdsyevx, PDSYEVX, double, double);
CSCALAPACK_PSYEVX_DECL(pcheevx, PCHEEVX, lapack_complex_float, float);
CSCALAPACK_PSYEVX_DECL(pzheevx, PZHEEVX, lapack_complex_double, double);

#undef CSCALAPACK_PSYEVX_DECL

#ifdef ROKKO_HAVE_SCALAPACK_PDSYEVR
#define CSCALAPACK_PSYEVR_DECL(NAMES, NAMEL, TYPE, TYPE_REAL) \
int cscalapack_## NAMES ##_work(char jobz, char range, char uplo, int n, \
                                TYPE* A, int ia, int ja, const int* descA, \
                                TYPE_REAL vl, TYPE_REAL vu, int il, int iu, \
                                int* m, int* nz, \
                                TYPE_REAL* w, TYPE* Z, int iz, int jz, const int* descZ, \
                                TYPE* work, int lwork, int* iwork, int liwork); \
int cscalapack_## NAMES (char jobz, char range, char uplo, int n, \
                         TYPE* A, int ia, int ja, const int* descA, \
                         TYPE_REAL vl, TYPE_REAL vu, int il, int iu, \
                         int* m, int* nz, \
                         TYPE_REAL* w, TYPE* Z, int iz, int jz, const int* descZ);

CSCALAPACK_PSYEVR_DECL(pssyevr, PSSYEVR, float, float);
CSCALAPACK_PSYEVR_DECL(pdsyevr, PDSYEVR, double, double);
CSCALAPACK_PSYEVR_DECL(pcheevr, PCHEEVR, lapack_complex_float, float);
CSCALAPACK_PSYEVR_DECL(pzheevr, PZHEEVR, lapack_complex_double, double);

#undef CSCALAPACK_PSYEVR_DECL
#endif

float cscalapack_pselget(char scope, char top, const float* A, int ia, int ja, const int* descA);

void cscalapack_pselset(float* A, int ia, int ja, const int* descA, float alpha);

float cscalapack_pslamch(int icnt, char cmch);

float cscalapack_pslange(char norm, int m, int n, const float* A, const int* descA);

float cscalapack_pslange_work(char norm, int m, int n, const float* A, const int* descA,
                              float* work);

void cscalapack_pslaprnt(int m, int n, const float* A, int ia, int ja, const int* descA,
                         int irprnt, int icprnt, const char* cmatnm, int nout, float* work);

#define CSCALAPACK_PSTEBZ_DECL(NAMES, NAMEL, TYPE) \
int cscalapack_## NAMES ##_work(int ictxt, char range, char order, int n, \
                                TYPE vl, TYPE vu, int il, int iu, \
                                TYPE abstol, const TYPE* d, const TYPE* e, int* m, int* nsplit, \
                                TYPE* w, int* iblock, int* isplit, \
                                TYPE* work, int lwork, int* iwork, int liwork); \
int cscalapack_## NAMES (int ictxt, char range, char order, int n, \
                         TYPE vl, TYPE vu, int il, int iu, \
                         TYPE abstol, const TYPE* d, const TYPE* e, int* m, int* nsplit, \
                         TYPE* w, int* iblock, int* isplit);

CSCALAPACK_PSTEBZ_DECL(psstebz, PSSTEBZ, float);
CSCALAPACK_PSTEBZ_DECL(pdstebz, PDSTEBZ, double);

#undef CSCALAPACK_PSTEBZ_DECL

#define CSCALAPACK_PSTEIN_DECL(NAMES, NAMEL, TYPE, TYPE_REAL) \
int cscalapack_## NAMES ##_work(int n, const TYPE_REAL* d, const TYPE_REAL* e, int m, \
                                TYPE_REAL* w, const int* iblock, const int* isplit, TYPE_REAL orfac, \
                                TYPE* Z, const int* iZ, const int* jZ, const int* descZ, \
                                TYPE* work, int lwork, int* iwork, int liwork, \
                                int* ifail, int* iclustr, TYPE_REAL* gap); \
int cscalapack_## NAMES (int n, const TYPE_REAL* d, const TYPE_REAL* e, int m, \
                         TYPE_REAL* w, const int* iblock, const int* isplit, TYPE_REAL orfac, \
                         TYPE* Z, const int* iZ, const int* jZ, const int* descZ, \
                         int* ifail, int* iclustr, TYPE_REAL* gap);

CSCALAPACK_PSTEIN_DECL(psstein, PSSTEIN, float, float);
CSCALAPACK_PSTEIN_DECL(pdstein, PDSTEIN, double, double);
CSCALAPACK_PSTEIN_DECL(pcstein, PCSTEIN, lapack_complex_float, float);
CSCALAPACK_PSTEIN_DECL(pzstein, PZSTEIN, lapack_complex_double, double);

#undef CSCALAPACK_PSTEIN_DECL

#define CSCALAPACK_PSTEDC_DECL(NAMES, NAMEL, TYPE) \
int cscalapack_## NAMES ##_work(char compz, int n, TYPE* d, TYPE* e, \
                                TYPE* Q, int iq, int jq, const int* descQ, \
                                TYPE* work, int lwork, int* iwork, int liwork); \
int cscalapack_## NAMES (char compz, int n, TYPE* d, TYPE* e, \
                         TYPE* Q, int iq, int jq, const int* descQ);

CSCALAPACK_PSTEDC_DECL(psstedc, PSSTEDC, float);
CSCALAPACK_PSTEDC_DECL(pdstedc, PDSTEDC, double);

#undef CSCALAPACK_PSTEDC_DECL

#ifdef __cplusplus
}
#endif

#endif // ROKKO_CSCALAPACK_H

max_stars_repo_name: t-sakashita/rokko

max_issues_repo_name: t-sakashita/rokko

max_forks_repo_name: t-sakashita/rokko


text:
#pragma once

#include <gsl.h>
#include <memory>

namespace crimson {

class ISolverSetupManager;

/*! \brief   A solver setup service interface. */
class ISolverSetupService
{
public:

    /*!
     * \brief   Gets solver setup manager.
     */
    virtual gsl::not_null<ISolverSetupManager*> getSolverSetupManager() const = 0;
    virtual ~ISolverSetupService() {}
};

}

ext: h

max_stars_repo_name: carthurs/CRIMSONGUI

max_issues_repo_name: carthurs/CRIMSONGUI

max_forks_repo_name: carthurs/CRIMSONGUI


text:
// 02.05.2021 - Colloid in a fluctuating scalar field
// Davide Venturelli & Benjamin Walter
// Last update 08.11.21

/* COMMENTS:
- Stochastic Runge-Kutta II for colloid evolution, Euler-Maruyama for field evolution (can be enhanced, but the price is O(N) at least. We could even think of anisotropic resolution, i.e. better around the colloid).
- No boundary conditions on the colloid displacement; they only get enforced when locating the nearest site.
- Space is measured in units of the lattice spacing.
- I am saving X(t) in colloid_pos and Y(t) in colloid_msd. Printing only X(t) at the moment.
- Finite size colloid, Gaussian interaction potential.
- Zero set in correspondence of trap center when printing out colloid data.
- Random inizialitazion of the field in Fourier space
- There is no limit on the size of R (other than physical, L/2)
- Compile as $ gcc -o colloid colloid_RF.c -lm -lgsl -lfftw3 -O3 -Wall -Wextra -O3
*/

// LIBRARIES, TYPES, DEFINITIONS

#include <stdio.h>
#include <math.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>
#include <gsl/gsl_sf.h>
#include <time.h>
#include <unistd.h>
#include <fftw3.h>

#define ALLOC(p,n) (p)=malloc( (n) * sizeof(*(p))); if( (p) == NULL){printf("Allocation of '%s' failed. Terminate. \n", #p); exit(2); }
#define CALLOC(p,n) (p)=calloc( (n) , sizeof(*(p))); if( (p) == NULL){printf("Allocation of '%s' failed. Terminate. \n", #p); exit(2); }
//#define long unsigned long						 					// This is to make long hold twice as much
#define DD printf("# Debug: line %d \n",__LINE__);
#define DEBUG (1)
#define EPS 0.0000000001												// Used to compare non-integer variables to 0

// Data types
typedef struct{
	double mass;														// Mass of the field
	double lambda;														// Field-colloid coupling strength
	double quartic_u;													// Self-interaction coupling strength
	double temperature;													// Temperature of the common bath
	double relativeD;													// Ratio of colloid to field mobility (it's basically \nu of the particle)
	double trap_strength;												// Stiffness of the harmonic trap
	int rng_seed;														// Seed for random number generator
	int system_size;													// Side of the DIM-dimensional lattice
	double delta_t;														// Time-discretization step
	long n_timestep;													// Number of timesteps
	int mc_runs;														// Number of Monte Carlo iterations
	double R;															// Size of colloid in lattice units
	double X0;															// Initial displacement of the colloid
} parameter;

typedef struct{															// This structure bundles all observables
	double** field_average;												// Saves the measured field-average <phi[i]> for certain subset of i's (eg along an axis) AND at all writing times
	double** field_correlation;											// Saves the measured field-correlator < phi[0] phi[i]> for certain subset of i's (eg along an axis) AND at all writing times
	double* colloid_pos;												// Saves the measured displacement of the colloid at all writing times
	double* colloid_msd;												// Saves the measured mean square displacement of the colloid at all writing times
	double write_time_delta;											// This is the gap between writing times (at the moment linear, maybe later exponential writing times?)
	int write_count;
//	double* colloid_fpt_distribution;									// Saves the first-passage time distribution of the colloid
} observables;

// Global variables
int DIM = 1;															// Physical dimensions
int MOD = 2;															// 0 for Model A, 2 for Model B
int top;																// Actual number of sites within 5R around a single site
double gauss_norm;														// Normalization constant of Gaussian functions
const gsl_rng_type *T;													// GSL RNG (it turns out Ziggurat is faster than Box-Muller)
gsl_rng *seed;

// FUNCTION PROTOTYPES

void default_parameters(parameter*);
void initialise(long***, parameter*, long***);
void initialise_observable(observables*, parameter*);
void wipe(double**, double**, parameter*);
void prethermalize(double**, parameter*);
void field_prepare(double**, parameter*);
void evolveB(double**, double**, long**, parameter*, observables*, long**);
void evolveA(double**, double**, long**, parameter*, observables*, long**);
void laplacian(double**, double*, long**, long);
void laplacian_of_cube(double**, double*, long**, long);
void generate_noise_field(double**, long, parameter*);
void gradient_field(double**, double*, long**, long);
void phi_evolveB(double**, double*, double*, double*, double*, long, parameter*, double*, long**);
void phi_evolveA(double**, double*, double*, long, parameter*, double*, long**);
void measure(double**, double**, long, parameter*, observables*);
void print_observables(observables*, parameter*);
void print_trajectory(observables*, parameter*);
void print_params(parameter*);
void printhelp(parameter*);
void print_source(void);
void neighborhood(long***, int);
int ind2coord(int, int, int);
int vec2ind(int*, int);
void ind2vecROW(int**, int, int, int);
int closest_site(double*, int);
double distance(double*, int, int);
unsigned modulo(int, unsigned);
double modular(double, double, double);
int factorial(int);
int double_factorial(int);
int intpow(int, int);
double floatpow(double, int);
void print_lookup(long **, int, int);
double gaussian(double, double);
double laplacian_gaussian(double, double);
void evolve_quenched(double**, double**, parameter*, observables*, long**);

// MAIN BODY

int main(int argc, char *argv[]){
	setlinebuf(stdout);

	// INPUT PARAMETERS
	parameter params;
	default_parameters(&params);

	opterr = 0;
	int c = 0;
	// ./colloid -r 2 -L 23 (order doesn't count)
    while( (c = getopt (argc, argv, "L:r:l:u:T:d:k:S:t:N:M:n:m:R:X:hP") ) != -1){
    	switch(c){
			case 'L':
				params.system_size = atoi(optarg);
				break;
			case 'r':
				params.mass = atof(optarg);
				break;
			case 'l':
				params.lambda = atof(optarg);
				break;
			case 'u':
				params.quartic_u = atof(optarg);
				break;
			case 'T':
				params.temperature = atof(optarg);
				break;
			case 'n':
				params.relativeD = atof(optarg);
				break;
			case 'k':
				params.trap_strength = atof(optarg);
				break;
			case 'S':
				params.rng_seed = atoi(optarg);
				break;
			case 't':
				params.delta_t = atof(optarg);
				break;
			case 'N':
				params.n_timestep = atol(optarg);
				break;
			case 'M':
				params.mc_runs = atoi(optarg);
				break;
			case 'd':
				DIM = atoi(optarg);
				break;
			case 'm':
				MOD = atoi(optarg);
				break;
			case 'R':
				params.R = atof(optarg);
				break;
			case 'X':
				params.X0 = atof(optarg);
				break;
			case 'P':
				print_source();
				exit(2);
			case 'h':
				printhelp(&params);
				exit(2);
			default:
				printhelp(&params);
                exit(EXIT_FAILURE);
       }
	}

	// VARIABLES, ALLOCATION
	int i, j;
	long L = params.system_size;
	long n_sites = intpow(L, DIM);
	int max_within_5R;
	double prefactor[3] = {2,M_PI,4.0/3*M_PI};
	if(DIM<4) max_within_5R = (int) ceil( prefactor[DIM-1]*floatpow(ceil(5*params.R), DIM) ) + 2;
	else max_within_5R=10000;											// Overkill, maybe I'll generalize it one day

	double* phi; 														// Field on lattice NxNxN
	double* y_colloid; 													// Colloid position (DIM real numbers)
	long** neighbours; 													// i x j - table with j neighbours of site i
	long** mosaic;														// Table with neighbours of site i with distance < 4R

	ALLOC(phi, n_sites);
	CALLOC(y_colloid, DIM);
	ALLOC(neighbours, n_sites);
	ALLOC(mosaic, n_sites);

	for(i = 0; i < n_sites; i++){
		neighbours[i]=calloc( 2 * DIM , sizeof(long)); 					// At each index there are 2*D neighbours
		if( neighbours[i] == NULL){printf("Allocation of neighbour list failed. Terminate. \n"); exit(2);}
		mosaic[i]=calloc( max_within_5R , sizeof(long)); 				// At each index, all the sites within 5R
		if( mosaic[i] == NULL){printf("Allocation of mosaic list failed. Terminate. \n"); exit(2);}
	}

	// INITIALIZATION
	initialise(&neighbours, &params, &mosaic);							// Random function, nearest neighbours list, mosaic
	observables obvs;													// Creates a pointer to an observables structure
	initialise_observable(&obvs, &params);								// Initialise observables
	print_params(&params);   											// Print header with all parameters
	gauss_norm = 1/(pow(2*M_PI,DIM*0.5) * floatpow(params.R,DIM));		// Normalization for Gaussian functions

	printf("# Max within 5R: %d\n", max_within_5R);
	printf("# Actual within 5R: %d\n", top);
	//print_lookup(neighbours, n_sites, 2*DIM);							// Prints nearest-neighbours list
	//print_lookup(mosaic, n_sites, top);								// Prints cell-list

	// MC ITERATION
	int mc_counter, flag;
	printf("\n# MC ITERATION BEGINS\n");

	for(mc_counter = 0; mc_counter < params.mc_runs; mc_counter++){
		//wipe(&phi, &y_colloid, &params); 								// Reset field and colloid to initial conditions
		prethermalize(&phi, &params);									// Pre-thermalization cycle
		//field_prepare(&phi, &params);
		for(j=0; j<DIM; j++) y_colloid[j] = L/2;						// Colloid initially in the middle of the trap
		y_colloid[0] = (L/2)+params.X0;									// Add initial displacement in one direction

		// Numerical integration of the dynamics

		if(MOD==0){														// Model A
			evolveA(&phi, &y_colloid, neighbours, &params, &obvs, mosaic);
		} else if(MOD==2){												// Model B
			evolveB(&phi, &y_colloid, neighbours, &params, &obvs, mosaic);
		} else if(MOD==1){												// Evolution in a quenched potential
			evolve_quenched(&phi, &y_colloid, &params, &obvs, mosaic);
		} else {
			printhelp(&params);
            exit(EXIT_FAILURE);
		}

		flag = params.mc_runs/10;										// Prints completing percentage
		if( flag !=0 && ((mc_counter +1)%flag) == 0 ) printf("# MC PROGRESS %d%%\n", (mc_counter +1)/flag*10);
	}

	//print_observables(&obvs, &params);
	print_trajectory(&obvs, &params);

	// Free memory and exit
	free(phi);
	free(y_colloid);
	for(i=0; i<n_sites; i++) free(neighbours[i]);
	for(i=0; i<n_sites; i++) free(mosaic[i]);
	free(obvs.colloid_pos);												// Add any other observable you are storing
	return 0;
}

// DEFINITION OF FUNCTIONS

// Default parameters
void default_parameters(parameter* params){
	params->mass = 0.0;
	params->lambda = 0.25;
	params->quartic_u = 0.0;
	params->temperature = 0.001;
	params->relativeD = 1.0;
	params->trap_strength = 0.1;
	params->rng_seed = -1; 												// if seed is -1 (not given by user), it will be picked randomly
	params->system_size = 128;
	params->delta_t = 0.01;
	params->n_timestep = 100000;
	params->mc_runs = 10000;
	params->R = 1;
	params->X0 = 2;
}

// All that needs to be done once
void initialise(long*** neighbours, parameter* params, long*** mosaic){
	// i) GSL random number generator setup
	gsl_rng_env_setup();
    T = gsl_rng_default;
    seed = gsl_rng_alloc (T);											// This is probably not really a seed, but sticazzi
    time_t t;															// If no seed provided, draw a random one
    if(params->rng_seed==-1) params->rng_seed = (unsigned) time(&t) % 100000;
    gsl_rng_set(seed, params->rng_seed);

	// ii) What are each position's neighbours?
	long L = params->system_size;
	neighborhood(neighbours, L);

	// iii) Initialize the mosaic
	double dist;
	long n_sites = intpow(L, DIM);
	int i, j, d;
	for(i=0; i<n_sites; i++){
		top=0;															// Number of sites within 5R around a given site (global variable)
		for(j=0; j<n_sites; j++){
			dist=0;
			for(d=0; d<DIM; d++) dist += floatpow( modular( ind2coord(d,i,L), ind2coord(d,j,L), L) , 2);
			dist = sqrt(dist);
			if(dist < ceil(5*params->R) + EPS ){
				(*mosaic)[i][top]=j;
				top++;
			}
		}
	}
}

// Reset field and colloid to initial conditions
void wipe(double** phi, double** y_colloid, parameter* params)
{
	int n_sites, i;
	n_sites = intpow(params->system_size, DIM);
	for(i = 0; i < n_sites; i++){
		 //(*phi)[i] = gsl_ran_gaussian_ziggurat(seed,1.0);				// Infinite temperature state
		 (*phi)[i] = 0;													// Initially flat field - default after calloc()
		 //(*phi)[i] = (i % 2 ? 1 : -1);								// Staggered
	}

	// y - colloid (initially in the middle of the lattice, where the harmonic well stands)
	long L = params->system_size;
	for(i=0; i<DIM; i++) (*y_colloid)[i] = L/2;
}

// This initialises the observables structure later containing the measurements
void initialise_observable(observables* obvs, parameter* params)
{
	obvs->write_time_delta = 0.2; 										// This is in physical time units, so writing occurs every (write_time_delta / n_timestep) integration step

	// This counts how many writing events will occur in time (including t=0, thus + 1). If you choose exp distributed measurements, it still works (but it's overestimated).
	int writing_times = (int)(1 + (((params->n_timestep)*params->delta_t)/obvs->write_time_delta));

	//CALLOC(obvs->colloid_msd, writing_times); 						// Save MSD vs time
	CALLOC(obvs->colloid_pos, writing_times); 							// Save position vs time

	/*
	// FIELD MEASUREMENTS ONLY
	CALLOC(obvs->field_average, writing_times); 						// Prepare writing_times many arrays to store averages
	//CALLOC(obvs->field_correlation, writing_times); 					// Prepare writing_times many arrays to store correlations
	int i;
	for(i = 0; i < writing_times; i++)
	{
		//obvs->field_correlation[i] = calloc( params->system_size , sizeof(double) );
		//if( (obvs->field_correlation[i]) == NULL){printf("Allocation of '(obvs->field_correlation[%i])'  failed. Terminate. \n", i); exit(2);}
		obvs->field_average[i] = calloc( params->system_size , sizeof(double) );
		if( (obvs->field_average[i]) == NULL){
			printf("Allocation of '(obvs->field_average[%i])'  failed. Terminate. \n", i);
			exit(2);
		}
	}
	*/
}

// Perform measurements and print them out
void measure(double** phi, double** y_colloid, long tstep, parameter* params, observables* obvs){
	long i;
	printf("%g\t", tstep*params->delta_t);
	for(i = 0; i < params->system_size; i++){
		obvs->field_average[obvs->write_count][i] += (*phi)[i];
		obvs->field_correlation[obvs->write_count][i] += ((*phi)[0] * (*phi)[i]);
	}

	long n_sites = intpow(params->system_size , DIM);
	for(i = 0; i < DIM; i++){
		obvs->colloid_msd[obvs->write_count] += ((*y_colloid)[i] - ind2coord(i, n_sites/2 , params->system_size))*((*y_colloid)[i] - ind2coord(i, n_sites/2 , params->system_size));
	}
	//printf("\n");
}

// Creates list of nearest neighbours in DIM dimensions
void neighborhood(long*** list, int L){
	int k,d;
	int vec[DIM], neigh[DIM];

	for(k=0; k<intpow(L,DIM); k++){										// Cycles over lattice sites
		for(d=0; d<DIM; d++){											// Finds the DIM coordinates of the current site
			vec[d] = ind2coord(d,k,L);
			neigh[d] = vec[d];
		}
		for(d=0; d<DIM; d++){
			neigh[d] = (vec[d]+1)%L;									// Finds right neighbour
			(*list)[k][2*d] = vec2ind(neigh,L);							// Stores it in the neighbour list
			neigh[d] = (vec[d]+L-1)%L;									// Finds left neighbour
			(*list)[k][2*d+1] = vec2ind(neigh,L);						// Stores it in the neighbour list
			neigh[d] = vec[d];											// Restores local copy (prepares for next dimension)
		}
	}
}

// Pre-thermalization of the field
void prethermalize(double** phi, parameter* params){
	// Benjamin: 20/08/21, spectral thermalisation
	// \phi_k are complex normal random variables with <\phi_k> = 0, < |\phi_k|^2> = k_B T*(k^2+r)^-1, and \phi_k = \phi^*_{-k}
	// Using FFTW library, -lfftw3
	fftw_complex *phi_fourier;
	fftw_plan p;

	int n_sites = intpow(params->system_size, DIM);
	int i, j;

	phi_fourier = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * (n_sites/2 + 1)); // complex Hermitian array
	int *vec = malloc(DIM * sizeof(int));
	/* fill with random data */
	double inverse_variance, variance;
	double inv_L_square = (39.438/pow(((double) params->system_size), 2)); //(2\pi /L)^2
	double radius;
	double polar_angle;

	int centre = (n_sites/2);
	for(i = centre; i >=0 ; i--) // Fill half the reciprocal lattice without the zero mode first. The first entry is at the centre of the reciprocal lattice and therefore the zero mode (~ k_BT/r)
	{
		// What's the L2 norm of the reciprocal vector?
		inverse_variance = params->mass;
		ind2vecROW(&vec, i, DIM, params->system_size);
		for(j = 0; j < DIM; j++)
		{
			inverse_variance += inv_L_square*(vec[j]*vec[j]);
		}
		variance = (params->temperature/inverse_variance); // k_B T /(q^2 + r)

		if(inverse_variance > 0){
			radius = gsl_ran_gaussian_ziggurat(seed, sqrt(variance));
			polar_angle =  gsl_ran_flat(seed, 0, 6.2831); // pick random angle

			phi_fourier[centre-i][0] = radius*cos(polar_angle);
			phi_fourier[centre-i][1] = radius*sin(polar_angle);
		}
		else{ // Critical mode is force set to zero
			phi_fourier[centre-i][0] = 0.0;
			phi_fourier[centre-i][1] = 0.0;
		}
	}

	// _c2r transforms a half-complex array into real fft http://www.fftw.org/fftw3_doc/Real_002ddata-DFTs.html
	int lattice_dimensions[DIM];
	for(i = 0; i < DIM; i++){lattice_dimensions[i] = ((int) params->system_size);} //It's a cube

	double* phi_real;
 	phi_real = (double*)fftw_malloc(n_sites * sizeof(double));

	// http://www.fftw.org/fftw3_doc/Real_002ddata-DFTs.html
	//fftw_plan fftw_plan_dft_c2r(int rank, const int *n, fftw_complex *in, double *out, unsigned flags);
	p = fftw_plan_dft_c2r(DIM, lattice_dimensions, phi_fourier, phi_real,  FFTW_ESTIMATE);
	if(p == NULL){printf("fftw plan didn't work\n");}
	fftw_execute(p);

	// normalising inverse transform with V^(-d/2), ass fftw3 doesn't do that.
	double normalising_volume = (1.0/(((double) n_sites)));

	// Now transform rowmajor phi_real into col_major *phi -> matrix transpose
	/* Only for D > 1
	 * int row_maj_index, col_maj_index, dir;
	for(row_maj_index = 0; row_maj_index < n_sites; row_maj_index++)
	{
		for(dir = 0; dir < DIM; dir++)
		{
			vec[dir] = ((int) (row_maj_index / intpow(params->system_size, (DIM - dir -1))) % params->system_size);
		}
	       	col_maj_index = vec2ind(vec, params->system_size);
		if(row_maj_ind
		printf("RM %i CM %i | phi_real: %g\n", row_maj_index, col_maj_index, phi_real[row_maj_index]);
		(*phi)[col_maj_index] = phi_real[row_maj_index];
	}*/

	// For D = 1
	// Shift origin to centre
	for(i = 0; i <= centre; i++) (*phi)[i] = normalising_volume*phi_real[centre-i];
	for(i = centre+1; i < n_sites; i++) (*phi)[i] = normalising_volume*phi_real[n_sites+centre-i];

    fftw_free(phi_fourier);
    fftw_free(phi_real);
	fftw_destroy_plan(p);
}

// Preparation of the field in equilibrium at T=0 around X0
// 1d only at the moment. It can probably be generalized to finite T by adding noise later on.
void field_prepare(double** phi, parameter* params){
	// Spectral initialization using FFTW library, -lfftw3
	// \phi_k are complex normal random variables with <\phi_k> = 0, < |\phi_k|^2> = k_B T*(k^2+r)^-1, and \phi_k = \phi^*_{-k}
	fftw_complex *phi_fourier;
	fftw_plan p;

	int n_sites = intpow(params->system_size, DIM);
	int i, j;

	phi_fourier = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * (n_sites/2 + 1)); // complex Hermitian array
	int *vec = malloc(DIM * sizeof(int));
	/* fill with random data */
	//double inverse_variance, variance;
	double inverse_variance;
	double inv_L = 6.28318530718/params->system_size;									// 2pi/L
	double inv_L_square = (39.438/intpow(params->system_size, 2)); 						// (2\pi /L)^2
	double radius;
	double polar_angle;

	int centre = (n_sites/2);
	for(i = centre; i >=0 ; i--) // Fill half the reciprocal lattice without the zero mode first. The first entry is at the centre of the reciprocal lattice and therefore the zero mode (~ k_BT/r)
	{
		// What's the L2 norm of the reciprocal vector?
		inverse_variance = params->mass;
		ind2vecROW(&vec, i, DIM, params->system_size);
		for(j = 0; j < DIM; j++)
		{
			inverse_variance += inv_L_square*(vec[j]*vec[j]);
		}
		//variance = (params->temperature/inverse_variance); // k_B T /(q^2 + r)
		radius = params->lambda * exp(- floatpow(params->R,2)*(inverse_variance - params->mass)/2) /inverse_variance ; // l* exp(-q^2 R^2/2) / (q^2 + r)

		if(inverse_variance > 0){
			//radius = gsl_ran_gaussian_ziggurat(seed, sqrt(variance));
			//polar_angle =  gsl_ran_flat(seed, 0, 6.2831); 							// pick random angle
			polar_angle = - inv_L * vec[j] * (params->system_size/2+params->X0);		// Careful if you ever generalize to higher dimensions

			phi_fourier[centre-i][0] = radius*cos(polar_angle);
			phi_fourier[centre-i][1] = radius*sin(polar_angle);
		}
		else{ // Critical mode is force set to zero
			phi_fourier[centre-i][0] = 0.0;
			phi_fourier[centre-i][1] = 0.0;
		}
	}

	// _c2r transforms a half-complex array into real fft http://www.fftw.org/fftw3_doc/Real_002ddata-DFTs.html
	int lattice_dimensions[DIM];
	for(i = 0; i < DIM; i++){lattice_dimensions[i] = ((int) params->system_size);} //It's a cube

	double* phi_real;
 	phi_real = (double*)fftw_malloc(n_sites * sizeof(double));

	// http://www.fftw.org/fftw3_doc/Real_002ddata-DFTs.html
	//fftw_plan fftw_plan_dft_c2r(int rank, const int *n, fftw_complex *in, double *out, unsigned flags);
	p = fftw_plan_dft_c2r(DIM, lattice_dimensions, phi_fourier, phi_real,  FFTW_ESTIMATE);
	if(p == NULL){printf("fftw plan didn't work\n");}
	fftw_execute(p);

	// normalising inverse transform with V^(-d/2), ass fftw3 doesn't do that.
	double normalising_volume = (1.0/(((double) n_sites)));

	// Now transform rowmajor phi_real into col_major *phi -> matrix transpose
	/* Only for D > 1
	 * int row_maj_index, col_maj_index, dir;
	for(row_maj_index = 0; row_maj_index < n_sites; row_maj_index++)
	{
		for(dir = 0; dir < DIM; dir++)
		{
			vec[dir] = ((int) (row_maj_index / intpow(params->system_size, (DIM - dir -1))) % params->system_size);
		}
	       	col_maj_index = vec2ind(vec, params->system_size);
		if(row_maj_ind
		printf("RM %i CM %i | phi_real: %g\n", row_maj_index, col_maj_index, phi_real[row_maj_index]);
		(*phi)[col_maj_index] = phi_real[row_maj_index];
	}*/

	// For D = 1
	// Shift origin to centre
	for(i = 0; i <= centre; i++) (*phi)[i] = normalising_volume*phi_real[centre-i];
	for(i = centre+1; i < n_sites; i++) (*phi)[i] = normalising_volume*phi_real[n_sites+centre-i];

    fftw_free(phi_fourier);
    fftw_free(phi_real);
	fftw_destroy_plan(p);
}

// Evolution - colloid only (quenched potential)
void evolve_quenched(double** phi, double** y_colloid, parameter* params, observables* obvs, long** mosaic){
	// Preparing variables for field evolution (Euler-Maruyama)
	long tstep, n_timestep;
	n_timestep = params->n_timestep;

	// Preparing variables for colloid evolution (Stochastic Runge-Kutta II)
	int i, j, y_site;
	double grad, noise, F2, r, weight;
	double w[DIM], F1[DIM], Y[DIM];
	long L = params->system_size;
	double noise_intensity = sqrt(2.0 * params->temperature * params->relativeD * params->delta_t);
	for(j=0; j<DIM; j++) w[j] = L/2;												// Finds position of the harmonic well

	// Preparing measurement process
	int write_time_delta_step = (int) ((obvs->write_time_delta)/(params->delta_t));
	int next_writing_step = 0;
	obvs->write_count = 0;

	for(tstep = 0; tstep < n_timestep; tstep++){									// Time evolution

		// MEASUREMENT
		if(tstep == next_writing_step){
			obvs->colloid_pos[obvs->write_count] += (*y_colloid)[0];				// Save colloid position
			//obvs->colloid_msd[obvs->write_count] += (*y_colloid)[1];				// Save colloid position of 2nd coordinate
			next_writing_step += write_time_delta_step;								// Linearly distributed
			obvs->write_count++;
		}

		// EVOLUTION

		// i) Create local copy of the colloid variable
		for(i=0; i<DIM; i++) Y[i] = (*y_colloid)[i];

		// ii) Prediction step for the colloid
		y_site = closest_site(Y,L);													// Finds the lattice site closest to Y

		for(i=0; i<DIM; i++){ 														// Evolves each of the components separately
			// Compute gradient of the field under the colloid
			grad = 0;
			for(j=0; j<top; j++){
				r = distance(Y, mosaic[y_site][j], L);
				weight = gaussian(r,params->R);
				grad += 0.5 * ( (*phi)[mosaic[y_site][j]+intpow(L,i)] - (*phi)[mosaic[y_site][j]-intpow(L,i)] ) * weight;
			}

			// Compute temporary position
			noise = noise_intensity * gsl_ran_gaussian_ziggurat(seed, 1.0);
			F1[i] =  params->relativeD * (params->lambda*grad -params->trap_strength*( Y[i]-w[i] ));
			(*y_colloid)[i] += params->delta_t * F1[i] + noise;
		}

		// iv) Correction step for the colloid (with the evolved field)
		y_site = closest_site(*y_colloid,L);										// Compute new site of the colloid

		for(i=0; i< DIM; i++){														// Evolve each of the components separately
			// Compute new gradient
			grad = 0;
			for(j=0; j<top; j++){
				r = distance(*y_colloid, mosaic[y_site][j], L);
				weight = gaussian(r,params->R);
				grad += 0.5 * ( (*phi)[mosaic[y_site][j]+intpow(L,i)] - (*phi)[mosaic[y_site][j]-intpow(L,i)] ) * weight;
			}

			// Compute corrected contribution
			F2 = params->relativeD * ( params->lambda*grad -params->trap_strength*( (*y_colloid)[i]-w[i] ) );

			// Sum the two contributions: y_n+1 = y_n + 1/2(F1+F2)*dt + noise
			noise = noise_intensity * gsl_ran_gaussian_ziggurat(seed, 1.0);
			(*y_colloid)[i] = Y[i] + 0.5*(F1[i]+F2)*params->delta_t + noise;
		}
	}
}

// Model B evolution
void evolveB(double** phi, double** y_colloid, long** neighbours, parameter* params, observables* obvs, long** mosaic){
	// Preparing variables for field evolution (Euler-Maruyama)
	long tstep, n_sites, n_timestep;
	n_sites = intpow(params->system_size, DIM);
	n_timestep = params->n_timestep;

	// TODO it's a terrible idea to allocate memory here, because this function is called M times
	double* laplacian_phi;
	ALLOC(laplacian_phi, n_sites);
	double* laplacian_square_phi;
	ALLOC(laplacian_square_phi, n_sites);
	double* laplacian_phi_cubed;
	CALLOC(laplacian_phi_cubed, n_sites);
	double* noise_field;
	CALLOC(noise_field, DIM * n_sites);
	double* noise_gradient;
	CALLOC(noise_gradient, n_sites);

	// Preparing variables for colloid evolution (Stochastic Runge-Kutta II)
	int i, j, y_site;
	double grad, noise, F2, r, weight;
	double w[DIM], F1[DIM], Y[DIM];
	long L = params->system_size;
	double noise_intensity = sqrt(2.0 * params->temperature * params->relativeD * params->delta_t);
	for(j=0; j<DIM; j++) w[j] = L/2;												// Finds position of the harmonic well

	// Preparing measurement process
	int write_time_delta_step = (int) ((obvs->write_time_delta)/(params->delta_t));
	int next_writing_step = 0;
	obvs->write_count = 0;

	for(tstep = 0; tstep < n_timestep; tstep++){									// Time evolution

		// MEASUREMENT

		if(tstep == next_writing_step){
			//measure(phi, y_colloid, tstep, params, obvs); 						// Evaluate all sorts of correlators etc. here
			//for(i = 0; i < params->system_size; i++) obvs->field_average[obvs->write_count][i] += (*phi)[i]; // Save field average
			obvs->colloid_pos[obvs->write_count] += (*y_colloid)[0];				// Save colloid position
			//obvs->colloid_msd[obvs->write_count] += (*y_colloid)[1];				// Save colloid position of 2nd coordinate
			next_writing_step += write_time_delta_step;								// Linearly distributed
			obvs->write_count++;
		}

		// EVOLUTION

		// i) Create local copy of the colloid variable
		for(i=0; i<DIM; i++) Y[i] = (*y_colloid)[i];

		// ii) Prediction step for the colloid
		y_site = closest_site(Y,L);													// Finds the lattice site closest to Y

		for(i=0; i<DIM; i++){ 														// Evolves each of the components separately
			// Compute gradient of the field under the colloid
			grad = 0;
			for(j=0; j<top; j++){
				r = distance(Y, mosaic[y_site][j], L);
				weight = gaussian(r,params->R);
				grad += 0.5 * ( (*phi)[mosaic[y_site][j]+intpow(L,i)] - (*phi)[mosaic[y_site][j]-intpow(L,i)] ) * weight;
			}

			// Compute temporary position
			noise = noise_intensity * gsl_ran_gaussian_ziggurat(seed, 1.0);
			F1[i] =  params->relativeD * (params->lambda*grad -params->trap_strength*( Y[i]-w[i] ));
			(*y_colloid)[i] += params->delta_t * F1[i] + noise;
		}

		// iii) Evolve the field with the local copy of the colloid position
		laplacian(&laplacian_phi, *phi, neighbours, n_sites);						// Write Laplacian of phi into laplacian_phi
		laplacian(&laplacian_square_phi, laplacian_phi, neighbours, n_sites);		// Write (D^2)^2 phi into laplacian_square_phi
		if(params->quartic_u > EPS){
			laplacian_of_cube(&laplacian_phi_cubed, *phi, neighbours, n_sites);		// Write D2 [phi(x)^3] into laplacian_phi_cubed
		}
		if(params->temperature > EPS){
			generate_noise_field(&noise_field, DIM*n_sites, params);				// Fill \vec{Lambda} with randomness
			gradient_field(&noise_gradient, noise_field, neighbours, n_sites);		// Compute gradient noise term
		}

		// Add together to new step d/dt phi = -a * D2 phi - b D4 phi - u D2 (phi^3) + D * noise (D is nabla)
		phi_evolveB(phi, laplacian_phi, laplacian_square_phi, laplacian_phi_cubed, noise_gradient, n_sites, params, Y, mosaic);

		// iv) Correction step for the colloid (with the evolved field)
		y_site = closest_site(*y_colloid,L);										// Compute new site of the colloid

		for(i=0; i< DIM; i++){														// Evolve each of the components separately
			// Compute new gradient
			grad = 0;
			for(j=0; j<top; j++){
				r = distance(*y_colloid, mosaic[y_site][j], L);
				weight = gaussian(r,params->R);
				grad += 0.5 * ( (*phi)[mosaic[y_site][j]+intpow(L,i)] - (*phi)[mosaic[y_site][j]-intpow(L,i)] ) * weight;
			}

			// Compute corrected contribution
			F2 = params->relativeD * ( params->lambda*grad -params->trap_strength*( (*y_colloid)[i]-w[i] ) );

			// Sum the two contributions: y_n+1 = y_n + 1/2(F1+F2)*dt + noise
			noise = noise_intensity * gsl_ran_gaussian_ziggurat(seed, 1.0);
			(*y_colloid)[i] = Y[i] + 0.5*(F1[i]+F2)*params->delta_t + noise;
		}
	}

	// Free memory
	free(laplacian_phi);
	free(laplacian_square_phi);
	free(laplacian_phi_cubed);
	free(noise_field);
	free(noise_gradient);
}

// Model A evolution
void evolveA(double** phi, double** y_colloid, long** neighbours, parameter* params, observables* obvs, long** mosaic){
	// Preparing variables for field evolution (Euler-Maruyama)
	long tstep, n_sites, n_timestep;
	n_sites = intpow(params->system_size, DIM);
	n_timestep = params->n_timestep;

	double* laplacian_phi;
	ALLOC(laplacian_phi, n_sites);
	double* noise_field;
	CALLOC(noise_field, n_sites);

	int i, j, y_site;
	double grad, noise, F2, r, weight;
	double w[DIM], F1[DIM], Y[DIM];
	long L = params->system_size;
	double noise_intensity = sqrt(2.0 * params->temperature * params->relativeD * params->delta_t);

	for(i=0; i<DIM; i++) w[i] = L/2;												// Finds position of the harmonic well

	// Preparing measurement process
	int write_time_delta_step = (int) ((obvs->write_time_delta)/(params->delta_t));
	int next_writing_step = 0;														// Linearly distributed writing times
	//int next_writing_step = write_time_delta_step;								// Exponentially distributed writing times TODO
	obvs->write_count = 0;

	for(tstep = 0; tstep < n_timestep; tstep++){									// Time evolution

		// MEASUREMENT

		if(tstep == next_writing_step){
			//measure(phi, y_colloid, tstep, params, obvs); 						// Evaluate all sorts of correlators etc. here
			obvs->colloid_pos[obvs->write_count] += (*y_colloid)[0];				// Save colloid position
			//for(i = 0; i < params->system_size; i++) obvs->field_average[obvs->write_count][i] += (*phi)[i]; // Save field average
			//obvs->colloid_msd[obvs->write_count] += (*y_colloid)[1];				// Save colloid position of 2nd coordinate
			next_writing_step += write_time_delta_step;								// Linearly distributed
			//next_writing_step = (int)(next_writing_step * exp( write_time_delta_step));	// Exponentially distributed TODO
			obvs->write_count++;
		}

		// EVOLUTION

		// i) Create local copy of the colloid variable
		for(i=0; i<DIM; i++) Y[i] = (*y_colloid)[i];

		// ii) Prediction step for the colloid
		y_site = closest_site(Y,L);													// Finds the lattice site closest to Y

		for(i=0; i<DIM; i++){ 														// Evolve each of the components separately
			// Compute gradient of the field under the colloid
			grad = 0;
			for(j=0; j<top; j++){
				r = distance(Y, mosaic[y_site][j], L);
				weight = gaussian(r,params->R);
				grad += 0.5 * ( (*phi)[mosaic[y_site][j]+intpow(L,i)] - (*phi)[mosaic[y_site][j]-intpow(L,i)] ) * weight;
			}

			// Compute temporary position
			noise = noise_intensity * gsl_ran_gaussian_ziggurat(seed, 1.0);
			F1[i] =  params->relativeD * (params->lambda*grad -params->trap_strength*( Y[i]-w[i] ));
			(*y_colloid)[i] += params->delta_t * F1[i] + noise;
		}

		// iii) Evolve the field with the local copy of the colloid position
		laplacian(&laplacian_phi, *phi, neighbours, n_sites);						// Write Laplacian of phi into laplacian_phi
		if(params->temperature > EPS){
			generate_noise_field(&noise_field, n_sites, params);					// Fill Lambda with randomness
		}

		// Add together to new step d/dt phi = -r * phi + \nabla^2 phi - u * (phi^3) + l * V(x-Y) + noise
		phi_evolveA(phi, laplacian_phi, noise_field, n_sites, params, Y, mosaic);

		// iv) Correction step for the colloid (with the evolved field)
		y_site = closest_site(*y_colloid,L);										// Compute new site of the colloid

		for(i=0; i< DIM; i++){														// Evolve each of the components separately
			// Compute new gradient
			grad = 0;
			for(j=0; j<top; j++){
				r = distance(*y_colloid, mosaic[y_site][j], L);
				weight = gaussian(r,params->R);
				grad += 0.5 * ( (*phi)[mosaic[y_site][j]+intpow(L,i)] - (*phi)[mosaic[y_site][j]-intpow(L,i)] ) * weight;
			}

			// Compute corrected contribution
			F2 = params->relativeD * (params->lambda*grad -params->trap_strength*( (*y_colloid)[i]-w[i] ));

			// Sum the two contributions: y_n+1 = y_n + 1/2(F1+F2)*dt + noise
			noise = noise_intensity * gsl_ran_gaussian_ziggurat(seed, 1.0);
			(*y_colloid)[i] = Y[i] + 0.5*(F1[i]+F2)*params->delta_t + noise;
		}
	}

	// Free memory
	free(laplacian_phi);
	free(noise_field);
}

// Field evolution - Model B
void phi_evolveB(double** phi, double* laplacian_phi, double* laplacian_square_phi, double* laplacian_phi_cubed, double* noise_gradient, long n_sites, parameter* params, double* Y, long** mosaic){
	long i;
	long L = params->system_size;
	double delta_t = params->delta_t;
	double r, weight;

	for(i = 0; i < n_sites; i++){
		(*phi)[i] += (delta_t*( -laplacian_square_phi[i] + params->mass*laplacian_phi[i] + params->quartic_u * laplacian_phi_cubed[i]) + noise_gradient[i]);
	}

	// Interaction with the colloid
	int y_site = closest_site(Y,L);													// Finds the lattice site closest to Y, the colloid
	for(i=0; i<top; i++){
		r = distance(Y, mosaic[y_site][i], L);
		weight = laplacian_gaussian(r,params->R);
		(*phi)[mosaic[y_site][i]] -= delta_t * params->lambda * weight;
	}
}

// Field evolution - Model A
void phi_evolveA(double** phi, double* laplacian_phi, double* noise_field, long n_sites, parameter* params, double* Y, long** mosaic){
	long i;
	double delta_t = params->delta_t;
	long L = params->system_size;
	double r, weight;

	for(i = 0; i < n_sites; i++){													// Notice noise_field contains delta_t in its variance
		(*phi)[i] += delta_t*(- params->mass*(*phi)[i] + laplacian_phi[i] - params->quartic_u * floatpow((*phi)[i],3) ) + noise_field[i];
	}

	// Interaction with the colloid
	int y_site = closest_site(Y,L);													// Finds the lattice site closest to Y, the colloid
	for(i=0; i<top; i++){
		r = distance(Y, mosaic[y_site][i], L);
		weight = gaussian(r,params->R);
		(*phi)[mosaic[y_site][i]] += delta_t * params->lambda * weight;
	}
}

// Returns Laplacian as calculated from cubic neighbour cells in DIM dimensions
void laplacian(double** laplacian,  double* field, long** neighbours, long n_sites){
	double buffer;
	long pos, i;
	for(pos = 0; pos < n_sites; pos++){
		buffer = 0;
		for(i = 0; i < 2*DIM; i++) {buffer += field[neighbours[pos][i]]; }
		buffer -= (2*DIM*field[pos]);
		(*laplacian)[pos] = buffer;
	}
}

// Returns the Laplacian of field^3
void laplacian_of_cube(double** laplacian,  double* field, long** neighbours, long n_sites){
	double buffer;
	long pos, i;
	for(pos = 0; pos < n_sites; pos++){
		buffer = 0;
		for(i = 0; i < 2*DIM; i++){
			buffer += pow( (field[neighbours[pos][i]]), 3);
		}
		buffer -= (2*DIM*field[pos]*field[pos]*field[pos]);
		(*laplacian)[pos] = buffer;
	}
}

// This function generates a completely uncorrelated random field on a line
void generate_noise_field(double** noise_field, long length, parameter* params){
	long i;
	double noise_intensity = sqrt(2.0 * params->temperature * params->delta_t);
	for(i = 0; i < length; i++){
		(*noise_field)[i] = noise_intensity * gsl_ran_gaussian_ziggurat(seed, 1.0);
	}
}

// Computes the gradient of a (noisy) field
void gradient_field(double** grad_noise, double* noise, long** neighbours, long n_sites){
	long i;
	int j; // neighbour of i
	double buffer;
	for(i = 0; i < n_sites; i++){
		buffer = 0;
		for(j = 0; j < DIM; j++){
			buffer += noise[neighbours[i][2*j]];
			buffer -= noise[neighbours[i][2*j+1]];
		}
		(*grad_noise)[i] = 0.5*buffer;
	}
}

// Returns the value of a Gaussian of variance R evaluated at r
double gaussian(double r, double R){
	return gauss_norm * exp(-r*r/(2*R*R));
}

// Returns the value of the laplacian of a Gaussian of variance R evaluated at r
double laplacian_gaussian(double r, double R){
	return gauss_norm * exp(-r*r/(2*R*R)) * (r*r - DIM*R*R) / floatpow(R,4);
}

// Translate from a list index (k) to d-th lattice index, and viceversa.
inline int ind2coord(int d, int k, int L) {return ( (int)( k/intpow(L,d) ) )%L ;}
int vec2ind(int *vec, int L){
	int i, res=0;
	for(i=0; i<DIM; i++) res += vec[i] * intpow(L,i);
	return res;
}

void ind2vecROW(int** vec, int ind, int dim, int L)
{
	// Places origin in the middle, row-major format.
	int dir;
	int offset = (int) ((L-1)/2.0); // So if L = 7 then the middle is at ___3___
	for(dir = 0; dir < dim; dir++)
	{
		(*vec)[dir] = (((int) (ind / intpow(L, (dim - dir -1))) % L) - offset) ; // row format, last index is fastest moving
	}
}

// Finds index of closest lattice site to the vector "vec"
// I'm using a special modulo function to avoid getting negative return values, eg -3 % 10 = -3, but modulo(-3,10) = 7.
int closest_site(double *vec, int L){
	int i, site=0;
	for(i=0; i<DIM; i++){
		site += (int)(modulo( round(vec[i])  , L) * intpow(L,i));
	}
	return site;
}

// According to https://stackoverflow.com/questions/14997165/fastest-way-to-get-a-positive-modulo-in-c-c this is still fast
unsigned modulo( int value, unsigned m) {
    int mod = value % (int)m;
    if (mod < 0) {
        mod += m;
    }
    return mod;
}

// Return the distance mod(L)
double modular(double a, double b, double L){
    return fabs(L/2 - fmod(3*L/2 + a - b, L));
}

// Returns the distance between the position vector and a given site (in a lattice of side L)
double distance(double* pos, int site, int L){
	int d;
	double dist=0;
	//for(d=0; d<DIM; d++) dist += floatpow( pos[d] - ind2coord(d,site,L) ,2);
	for(d=0; d<DIM; d++) dist += floatpow( modular(pos[d], ind2coord(d,site,L), L) , 2);
	dist = sqrt(dist);
	return dist;
}

// Factorial of a number
int factorial(int n){
    if (n == 0) return 1;
    return n * factorial(n - 1);
}

// Double factorial of a number
int double_factorial(int n){
    if (n <= 1) return 1;
    return n * factorial(n - 2);
}

// The usual pow(a,b)=exp(log(a) * b) is slow AF
int intpow(int a, int b){
	int i, res=1;
	for(i=0; i<b; i++) res *= a;
	return res;
}

// The usual pow(a,b)=exp(log(a) * b) is slow AF
double floatpow(double a, int b){
	int i;
	double res=1;
	for(i=0; i<b; i++) res *= a;
	return res;
}

// Prints observables including field
void print_observables(observables* obvs, parameter* params){
	int i,j;
	int system_size = params->system_size;
	int write_count = obvs->write_count;
	double weight = 1/((double) params->mc_runs);

	for(i = 0; i < write_count; i ++){
		printf("#FIELDAVG %g", i * (obvs->write_time_delta));
		for(j = 0; j < system_size; j++){
			printf("\t%g", weight*(obvs->field_average[i][j]));
		}
		printf("\n");
	}

	/*
	for(i = 0; i < write_count; i ++){
		printf("#FIELDCORR %g", i * (obvs->write_time_delta));
		for(j = 0; j < system_size; j++){
			printf("\t%g", weight*(obvs->field_correlation[i][j]- (obvs->field_average[i][0] * obvs->field_average[i][j]) ));
		}
		printf("\n");
	}

	// Output MSD of colloid
	for(i = 0; i < write_count; i++){
		printf("# COLLOIDMSD %g", i * (obvs->write_time_delta));
		printf("\t%.3f",weight*(obvs->colloid_msd[i]));
		printf("\n");
	}
	*/

	// Output trajectory of colloid
	for(i = 0; i < write_count; i++){
		printf("#COLLOID_X %g", i * (obvs->write_time_delta));
		//printf("%.4f", i * (obvs->write_time_delta));
		printf("\t%.15f",weight*obvs->colloid_pos[i]- params->system_size/2);
		//printf("\t%.12f",weight*(obvs->colloid_msd[i]));
		printf("\n");
	}
}

// Prints colloid-related observables only
void print_trajectory(observables* obvs, parameter* params)
{
	int i;
	int write_count = obvs->write_count;
	double weight = 1/((double) params->mc_runs);

	// Output trajectory of colloid
	for(i = 0; i < write_count; i++){
		//printf("#COLLOID_X %g", i * (obvs->write_time_delta));
		printf("%.4f", i * (obvs->write_time_delta));
		printf("\t%.15f", weight*obvs->colloid_pos[i]- params->system_size/2);
		//printf("\t%.12f",weight*(obvs->colloid_msd[i]));
		printf("\n");
	}
}

// Prints out simulation parameters
void print_params(parameter* params){
printf("# Parameters\n\
# MASS %g\n\
# LAMBDA (field-colloid-coupling) %g\n\
# U (quartic coupling) %g\n\
# TEMPERATURE %g\n\
# RELATIVE MOTILITY D %g\n\
# TRAP STRENGTH %g\n\
# RNG SEED %u\n\
# L (System Size) %i\n\
# DIM %i\n\
# RADIUS %g\n\
# INITIAL DISPLACEMENT %g\n\
# DELTA T %g\n\
# TIMESTEPS %lu\n\
# MONTE CARLO RUNS %i\n",
params->mass, params->lambda, params->quartic_u, params->temperature, params->relativeD, params->trap_strength, params->rng_seed, params->system_size, DIM, params->R, params->X0, params->delta_t, params->n_timestep, params->mc_runs);

	if(MOD==0) printf("# MODEL A\n");
	if(MOD==2) printf("# MODEL B\n");
}

// Prints out the instructions
void printhelp(parameter* params){
	printf("# Colloid in Gaussian Field\n# '%s' built %s\n\
# Use with options flags\n\
# -L Length of d-dimensional lattice (default %i)\n\
# -r Mass of Gaussian field (default %g)\n\
# -l Coupling strength between colloid and field (default %g)\n\
# -u Quartic coupling strength (default %g)\n\
# -T Temperature of bath (default %g)\n\
# -n Relative motility colloid/field (default %g)\n\
# -k Strength of harmonic trap (default %g)\n\
# -S Seed for RNG (default random)\n\
# -t Integration timestep (default %g)\n\
# -d Dimension (default %i)\n\
# -N Number of timesteps (default %lu)\n\
# -M Number of Monte Carlo samples (default %i)\n\
# -R Radius of the colloid (default %g)\n\
# -X Initial displacement of the colloid (default %g)\n\
# -m Field dynamics (model A -> 0, model B -> 2, quenched -> 1, default %i)\n\
# -h To see this helpscreen\n\
# -P Output source code\n",__FILE__,__DATE__, params->system_size, params->mass, params->lambda, params->quartic_u, params->temperature, params->relativeD, params->trap_strength, params->delta_t, DIM, params->n_timestep, params->mc_runs, params->R, params->X0, MOD);
}

// Prints out the nearest-neighbour list
// For neighbour list, size1=n_sites and size2=2*DIM. For mosaic, size1=n_sites and size2=top.
void print_lookup(long ** list, int size1, int size2){
	int i,j;
	for(i=0; i<size1; i++){
		printf("# Around site %d : \t",i);
		for(j=0; j<size2; j++) printf("%li\t", list[i][j]);
		printf("\n");
	}
}

// Prints out this whole code
void print_source(){
    printf("/* Source Code %s, created %s */\n",__FILE__,__DATE__);
    FILE *fp;
    int c;

    fp = fopen(__FILE__,"r");

    do {
     c = getc(fp);
 	 putchar(c);
 	}
 	while(c != EOF);
    fclose(fp);
}

ext: c


text:
/**
 * @file postfilter.h

 * @brief implementation of post-filters for a microphone array.

   The following post-filters are implemented:
   [1] Zelinski post-filter
   [2] APAB post-filter
   [3] McCowan's post-filter
   [4] Lefkimmiatis's post-filter

   The correspondig references are:
   [1] C.Claude Marro et al. "Analysis of noise reduction and dereverberation techniques based on microphone arrays with postfiltering", IEEE Trans. ASP, vol. 6, pp 240-259, May 1998.
   [2] M.Brandstein, "Microphone Arrays", Springer, ISBN 3-540-41953-5, pp.39-60.
   [3] Iain A. Mccowan et al., "Microphone array post-filter based on noise field coherence", IEEE Trans. SAP, vol. 11, pp. Nov. 709--716, 2003.
   [4] Stamatios Lefkimmiatis et al., "A generalized estimation approach for linear and nonlinear microphone array post-filters",  Speech Communication, 2007.

  * @author Kenichi Kumatani
*/
#ifndef POSTFILTER_H
#define POSTFILTER_H

#include <stdio.h>
#include <assert.h>
#include <float.h>

#include <gsl/gsl_block.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_complex.h>
#include <gsl/gsl_complex_math.h>
#include <gsl/gsl_fft_complex.h>
#include <common/refcount.h>
#include "common/jexception.h"

#include "stream/stream.h"
#include "postfilter/spectralsubtraction.h"
#include "postfilter/binauralprocessing.h"
#include "beamformer/spectralinfoarray.h"
#include "beamformer/beamformer.h"

typedef enum {
  TYPE_ZELINSKI1_REAL = 0x01,
  TYPE_ZELINSKI1_ABS  = 0x02,
  TYPE_APAB = 0x04,
  TYPE_ZELINSKI2 = 0x08,
  NO_USE_POST_FILTER = 0x00
} PostfilterType;

void ZelinskiFilter(gsl_vector_complex **arrayManifold,
		    SnapShotArrayPtr     snapShotArray,
		    bool halfBandShift,
		    gsl_vector_complex *beamformedSignal,
		    gsl_vector_complex **prevCSDs,
		    gsl_vector_complex *pfweights,
		    double alpha, int Ropt );

void ApabFilter( gsl_vector_complex **arrayManifold,
		 SnapShotArrayPtr     snapShotArray,
		 int fftLen, int nChan, bool halfBandShift,
		 gsl_vector_complex *beamformedSignal,
		 int channelX );

/**
   @class Zelinski post-filtering

   @brief filter beamformer's outputs under the assumption that noise signals between sensors are uncorrelated.
   @usage
   1. construct an object,
   2. set the array snapshot and array manifold vectors with
      either setBeamformer() or setSnapShotArray() and setArrayManifold(),
   3. process data at each frame by caling next().
*/
class ZelinskiPostFilter: public VectorComplexFeatureStream {
public:
  ZelinskiPostFilter( VectorComplexFeatureStreamPtr &output, unsigned fftLen, double alpha=0.6, int type=2, int minFrames=0, const String& nm = "ZelinskPostFilter" );
  ~ZelinskiPostFilter();

  virtual const gsl_vector_complex* next(int frame_no = -5);
  virtual void reset();

  void set_beamformer(SubbandDSPtr &beamformer);
  void set_snapshot_array(SnapShotArrayPtr &snapShotArray);
  void set_array_manifold_vector(unsigned fbinX, gsl_vector_complex *arrayManifoldVector, bool halfBandShift, unsigned NC = 1);
  const gsl_vector_complex* postfilter_weights(){
    if( NULL == bf_weights_ )
      return NULL;
    return(bf_weights_->wp1());
  }

#ifdef ENABLE_LEGACY_BTK_API
  void setBeamformer(SubbandDSPtr &beamformer){ set_beamformer(beamformer); }
  void setSnapShotArray(SnapShotArrayPtr &snapShotArray){ set_snapshot_array(snapShotArray); }
  void setArrayManifoldVector(unsigned fbinX, gsl_vector_complex *arrayManifoldVector, bool halfBandShift, unsigned NC = 1){ set_array_manifold_vector(fbinX, arrayManifoldVector, halfBandShift, NC); }
  const gsl_vector_complex* getPostFilterWeights(){ return postfilter_weights(); }
#endif

protected:
  unsigned                      fftLen_;
  VectorComplexFeatureStreamPtr samp_; /* output of the beamformer */
  PostfilterType                type_; /* the type of the Zelinski-postfilters */
  double                        alpha_; /* forgetting factor */
  int                           min_frames_;
  SubbandDSPtr                  bf_ptr_; /* */
  BeamformerWeights*            bf_weights_;
  bool                          has_bf_ptr_; /* true if bf_ptr_ is set with setBeamformer() */
  SnapShotArrayPtr              snapshot_array_; /* multi-channel input */
};

typedef Inherit<ZelinskiPostFilter, VectorComplexFeatureStreamPtr> ZelinskiPostFilterPtr;

/**
   @class McCowan post-filtering

   @brief process the beamformer's outputs with McCowan's post-filtering
   @usage
   1. construct an object,
   2. compute the noise coherence matrix through setDiffuseNoiseModel( micPositions, ssampleRate)
   3. set the array snapshot and array manifold vectors with
      either setBeamformer() or setSnapShotArray() and setArrayManifold(),
   4. process data at each frame by caling next().
*/
class McCowanPostFilter: public ZelinskiPostFilter {
public:
  McCowanPostFilter( VectorComplexFeatureStreamPtr &output, unsigned fftLen, double alpha=0.6, int type=2, int minFrames=0, float threshold=0.99, const String& nm = "McCowanPostFilter" );
  ~McCowanPostFilter();

  virtual const gsl_vector_complex* next(int frame_no = -5);
  virtual void reset();
  /* micPositions[][x,y,z] */

  const gsl_matrix_complex *noise_spatial_spectral_matrix( unsigned fbinX );
  bool set_noise_spatial_spectral_matrix( unsigned fbinX, gsl_matrix_complex* Rnn );
  bool set_diffuse_noise_model( const gsl_matrix* micPositions, double sampleRate, double sspeed = 343740.0 );
  void set_all_diagonal_loading(float diagonalWeight);
  void set_diagonal_looading(unsigned fbinX, float diagonalWeight);
  void divide_all_nondiagonal_elements(float mu);
  void divide_nondiagonal_elements(unsigned fbinX, float mu);

#ifdef ENABLE_LEGACY_BTK_API
  const gsl_matrix_complex *getNoiseSpatialSpectralMatrix(unsigned fbinX){ return noise_spatial_spectral_matrix(fbinX); }
  bool setNoiseSpatialSpectralMatrix(unsigned fbinX, gsl_matrix_complex* Rnn){ return set_noise_spatial_spectral_matrix(fbinX, Rnn); }
  bool setDiffuseNoiseModel(const gsl_matrix* micPositions, double sampleRate, double sspeed = 343740.0){ return set_diffuse_noise_model( micPositions, sampleRate, sspeed); }
  void setAllLevelsOfDiagonalLoading(float diagonalWeight){ set_all_diagonal_loading(diagonalWeight); }
  void setLevelOfDiagonalLoading(unsigned fbinX, float diagonalWeight){ set_diagonal_looading(fbinX, diagonalWeight); }
  void divideAllNonDiagonalElements( float mu );
  void divideNonDiagonalElements( unsigned fbinX, float mu );
#endif

protected:
  double estimate_average_clean_PSD_( unsigned fbinX, gsl_vector_complex* currCSDf );
  virtual void post_filtering_();

protected:
  gsl_matrix_complex**                           R_; /* Noise spatial spectral matrices */
  float*                                         diagonal_weights_;
  float                                          threshold_of_Rij_; /* to avoid the indeterminate solution*/
  gsl_vector_complex*                            time_aligned_signal_f_; /* workspace */
  bool                                           invR_computed_;
};

typedef Inherit<McCowanPostFilter, ZelinskiPostFilterPtr> McCowanPostFilterPtr;

/**
   @class Lefkimmiatis post-filtering
   @brief compute a Winer filter under the the diffuse noise field assumption
   @usage
   1. construct an object,
   2. compute the noise coherence matrix through setDiffuseNoiseModel( micPositions, ssampleRate)
   3. set the array snapshot and array manifold vectors with
      either setBeamformer() or setSnapShotArray() and setArrayManifold(),
   4. process data at each frame by caling next().
*/
class LefkimmiatisPostFilter: public McCowanPostFilter {
public:
  LefkimmiatisPostFilter( VectorComplexFeatureStreamPtr &output, unsigned fftLen, double minSV=1.0E-8, unsigned fbinX1=0, double alpha=0.6, int type=2, int minFrames=0, float threshold=0.99, const String& nm = "LefkimmiatisPostFilte" );
  ~LefkimmiatisPostFilter();

  virtual const gsl_vector_complex* next(int frame_no = -5);
  virtual void reset();
  /* micPositions[][x,y,z] */

  void calc_inverse_noise_spatial_spectral_matrix();

#ifdef ENABLE_LEGACY_BTK_API
  void calcInverseNoiseSpatialSpectralMatrix(){ calc_inverse_noise_spatial_spectral_matrix(); }
#endif

protected:
  double estimate_average_noise_PSD_( unsigned fbinX, gsl_vector_complex* currCSDf );
  virtual void post_filtering_();

private:
  gsl_complex calcLambda( unsigned fbinX );

  gsl_matrix_complex** invR_;
  gsl_vector_complex*  tmpH_;
  double              minSV_;
  unsigned            fbinX1_;
};

typedef Inherit<LefkimmiatisPostFilter, McCowanPostFilterPtr> LefkimmiatisPostFilterPtr;

/**
   @class high pass filter
*/
class HighPassFilter: public VectorComplexFeatureStream {
public:
  HighPassFilter( VectorComplexFeatureStreamPtr &output, float cutOffFreq, int sampleRate, const String& nm = "HighPassFilter" );
  virtual const gsl_vector_complex* next(int frame_no = -5);
  virtual void reset();

private:
  VectorComplexFeatureStreamPtr src_;
  unsigned cutoff_fbinX_;
};

typedef Inherit<HighPassFilter, VectorComplexFeatureStreamPtr> HighPassFilterPtr;

#endif


text:
/**
 * @file ssygvd.c
 *
 *  PLASMA computational routines
 *  Release Date: November, 15th 2009
 *  PLASMA is a software package provided by Univ. of Tennessee,
 *  Univ. of California Berkeley and Univ. of Colorado Denver
 *
 * @version 2.6.0
 * @author Azzam Haidar
 * @author Hatem Ltaief
 * @date 2010-11-15
 * @generated s Tue Jan  7 11:45:09 2014
 *
 **/
#include <lapacke.h>
#include "common.h"

/***************************************************************************//**
 *
 * @ingroup float
 *
 *  PLASMA_ssygvd - Computes all eigenvalues and, optionally,
 *  eigenvectors of a complex generalized Hermitian-definite
 *  eigenproblem of the form: A*x=(lambda)*B*x, A*Bx=(lambda)*x, or
 *  B*A*x=(lambda)*x.
 *  Here A and B are assumed to be Hermitian and B is also positive
 *  definite.
 *
 *******************************************************************************
 *
 * @param[in] itype
 *          Intended usage:
 *          = 1: A*x=(lambda)*B*x
 *          = 2: A*Bx=(lambda)*x
 *          = 3: B*A*x=(lambda)*x
 *
 * @param[in] jobz
 *          Intended usage:
 *          = PlasmaNoVec: computes eigenvalues only;
 *          = PlasmaVec: computes eigenvalues and eigenvectors.
 *
 * @param[in] uplo
 *          Specifies whether the matrix A is upper triangular or
 *          lower triangular:
 *          = PlasmaUpper: Upper triangle of A and B are stored;
 *          = PlasmaLower: Lower triangle of A and B are stored.
 *
 * @param[in] N
 *          The order of the matrix A. N >= 0.
 *
 * @param[in,out] A
 *          On entry, the symmetric (or Hermitian) matrix A.
 *          If uplo = PlasmaUpper, the leading N-by-N upper triangular
 *          part of A contains the upper triangular part of the matrix
 *          A, and the strictly lower triangular part of A is not
 *          referenced.
 *          If uplo = PlasmaLower, the leading N-by-N lower triangular
 *          part of A contains the lower triangular part of the matrix
 *          A, and the strictly upper triangular part of A is not
 *          referenced.
 *          On exit, if jobz = PlasmaVec, then if return value = 0, A
 *          contains the matrix Z of eigenvectors.
 *          The eigenvectors are normalized as follows:
 *          if ITYPE = 1 or 2, Z**T*B*Z = I;
 *          if ITYPE = 3,      Z**T*inv(B)*Z = I.
 *          If jobz = PlasmaNoVec, then on exit the lower triangle (if
 *          uplo = PlasmaLower) or the upper triangle (if uplo =
 *          PlasmaUpper) of A, including the diagonal, is destroyed.
 *
 * @param[in] LDA
 *          The leading dimension of the array A. LDA >= max(1,N).
 *
 * @param[in,out] B
 *          On entry, the symmetric (or Hermitian) positive definite
 *          matrix B.
 *          If uplo = PlasmaUpper, the leading N-by-N upper triangular
 *          part of B contains the upper triangular part of the matrix
 *          B, and the strictly lower triangular part of B is not
 *          referenced.
 *          If uplo = PlasmaLower, the leading N-by-N lower triangular
 *          part of B contains the lower triangular part of the matrix
 *          B, and the strictly upper triangular part of B is not
 *          referenced.
 *          On exit, if return value <= N, the part of B containing
 *          the matrix is overwritten by the triangular factor U or L
 *          from the Cholesky factorization B = U**T*U or B = L*L**T.
 *
 * @param[in] LDB
 *          The leading dimension of the array B. LDA >= max(1,N).
 *
 * @param[out] W
 *          On exit, if info = 0, the eigenvalues.
 *
 * @param[in, out] descT
 *          On entry, descriptor as return by PLASMA_Alloc_Workspace_ssygvd
 *          On exit, contains auxiliary factorization data.
 *
 * @param[out] Q
 *          On exit, if jobz = PlasmaVec and info = 0, the eigenvectors.
 *
 * @param[in] LDQ
 *          The leading dimension of Q.
 *
 *******************************************************************************
 *
 * @return
 *          \retval PLASMA_SUCCESS successful exit
 *          \retval <0 if -i, the i-th argument had an illegal value
 *          \retval <=N if INFO = i, plasma_ssygvd failed to converge; i
 *               off-diagonal elements of an intermediate tridiagonal
 *               form did not converge to zero.
 *          \retval >N if INFO = N + i, for 1 <= i <= N, then the leading
 *                     minor of order i of B is not positive definite.
 *                     The factorization of B could not be completed and
 *                     no eigenvalues or eigenvectors were computed.
 *
 *******************************************************************************
 *
 * @sa PLASMA_ssygvd_Tile
 * @sa PLASMA_ssygvd_Tile_Async
 * @sa PLASMA_chegvd
 * @sa PLASMA_dsygvd
 * @sa PLASMA_ssygvd
 *
 ******************************************************************************/
int PLASMA_ssygvd(PLASMA_enum itype, PLASMA_enum jobz, PLASMA_enum uplo, int N,
                 float *A, int LDA,
                 float *B, int LDB,
                 float *W,
                 PLASMA_desc *descT,
                 float *Q, int LDQ)
{
    int NB;
    int status;
    plasma_context_t *plasma;
    PLASMA_sequence *sequence = NULL;
    PLASMA_request request = PLASMA_REQUEST_INITIALIZER;
    PLASMA_desc descA, descB, descQ;

    plasma = plasma_context_self();
    if (plasma == NULL) {
        plasma_error("PLASMA_ssygvd", "PLASMA not initialized");
        return PLASMA_ERR_NOT_INITIALIZED;
    }

    /* Check input arguments */
    if (itype != 1 && itype != 2 && itype != 3) {
        plasma_error("PLASMA_ssygvd", "Illegal value of itype");
        return -1;
    }
    if (jobz != PlasmaNoVec && jobz != PlasmaVec) {
        plasma_error("PLASMA_ssygvd", "illegal value of jobz");
        return -2;
    }
    if (uplo != PlasmaLower && uplo!= PlasmaUpper) {
        plasma_error("PLASMA_ssygvd", "only PlasmaLower supported");
        return -3;
    }
    if (N < 0) {
        plasma_error("PLASMA_ssygvd", "illegal value of N");
        return -4;
    }
    if (LDA < max(1, N)) {
        plasma_error("PLASMA_ssygvd", "illegal value of LDA");
        return -6;
    }
    if (LDB < max(1, N)) {
        plasma_error("PLASMA_ssygvd", "illegal value of LDB");
        return -8;
    }
    if (LDQ < max(1, N)) {
        plasma_error("PLASMA_ssygvd", "illegal value of LDQ");
        return -12;
    }

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

    /* Tune NB & IB depending on N; Set NBNB */
    status = plasma_tune(PLASMA_FUNC_SSYGV, N, N, 0);
    if (status != PLASMA_SUCCESS) {
        plasma_error("PLASMA_ssygvd", "plasma_tune() failed");
        return status;
    }

    /* Set NT */
    NB = PLASMA_NB;

    plasma_sequence_create(plasma, &sequence);

    if ( PLASMA_TRANSLATION == PLASMA_OUTOFPLACE ) {
        plasma_sooplap2tile( descA, A, NB, NB, LDA, N, 0, 0, N, N, sequence, &request,
                             plasma_desc_mat_free(&(descA)) );
        plasma_sooplap2tile( descB, B, NB, NB, LDB, N, 0, 0, N, N, sequence, &request,
                             plasma_desc_mat_free(&(descA)); plasma_desc_mat_free(&(descB)) );
        /* No conversion for Q */
    } else {
        plasma_siplap2tile( descA, A, NB, NB, LDA, N, 0, 0, N, N,
                            sequence, &request);
        plasma_siplap2tile( descB, B, NB, NB, LDB, N, 0, 0, N, N,
                            sequence, &request);
    }

    /* Q is used only as an output */
    descQ = plasma_desc_init(
        PlasmaRealFloat, NB, NB, NB*NB,
        LDQ, N, 0, 0, N, N);
    descQ.mat = Q;

    /* Call the tile interface */
    PLASMA_ssygvd_Tile_Async(itype, jobz, uplo,
                            &descA, &descB, W,
                            descT, &descQ,
                            sequence, &request);

    plasma_siptile2lap( descQ, Q, NB, NB, LDQ, N,  sequence, &request);
    if ( PLASMA_TRANSLATION == PLASMA_OUTOFPLACE ) {
        plasma_sooptile2lap( descA, A, NB, NB, LDA, N, sequence, &request);
        plasma_sooptile2lap( descB, B, NB, NB, LDB, N, sequence, &request);
        plasma_dynamic_sync();
        plasma_desc_mat_free(&descA);
        plasma_desc_mat_free(&descB);
    } else {
        plasma_siptile2lap( descA, A, NB, NB, LDA, N, sequence, &request);
        plasma_siptile2lap( descB, B, NB, NB, LDB, N, sequence, &request);
        plasma_dynamic_sync();
    }

    status = sequence->status;
    plasma_sequence_destroy(plasma, sequence);
    return status;
}
/***************************************************************************//**
 *
 * @ingroup float_Tile
 *
 *  PLASMA_ssygvd_Tile - Computes all eigenvalues and, optionally,
 *  eigenvectors of a complex generalized Hermitian-definite
 *  eigenproblem of the form: A*x=(lambda)*B*x, A*Bx=(lambda)*x, or
 *  B*A*x=(lambda)*x.
 *  Here A and B are assumed to be Hermitian and B is also
 *  positive definite.
 *
 *  Tile equivalent of PLASMA_ssygvd().
 *  Operates on matrices stored by tiles.
 *  All matrices are passed through descriptors.
 *  All dimensions are taken from the descriptors.
 *
 *******************************************************************************
 *
 * @param[in] itype
 *          Intended usage:
 *          = 1: A*x=(lambda)*B*x
 *          = 2: A*Bx=(lambda)*x
 *          = 3: B*A*x=(lambda)*x
 *
 * @param[in] jobz
 *          Intended usage:
 *          = PlasmaNoVec: computes eigenvalues only;
 *          = PlasmaVec: computes eigenvalues and eigenvectors.
 *
 * @param[in] uplo
 *          Specifies whether the matrix A is upper triangular or
 *          lower triangular:
 *          = PlasmaUpper: Upper triangle of A and B are stored;
 *          = PlasmaLower: Lower triangle of A and B are stored.
 *
 * @param[in,out] A
 *          On entry, the symmetric (or Hermitian) matrix A.
 *          If uplo = PlasmaUpper, the leading N-by-N upper triangular
 *          part of A contains the upper triangular part of the matrix
 *          A, and the strictly lower triangular part of A is not
 *          referenced.
 *          If uplo = PlasmaLower, the leading N-by-N lower triangular
 *          part of A contains the lower triangular part of the matrix
 *          A, and the strictly upper triangular part of A is not
 *          referenced.
 *          On exit, if jobz = PlasmaVec, then if return value = 0, A
 *          contains the matrix Z of eigenvectors.
 *          The eigenvectors are normalized as follows:
 *          if ITYPE = 1 or 2, Z**T*B*Z = I;
 *          if ITYPE = 3,      Z**T*inv(B)*Z = I.
 *          If jobz = PlasmaNoVec, then on exit the lower triangle (if
 *          uplo = PlasmaLower) or the upper triangle (if uplo =
 *          PlasmaUpper) of A, including the diagonal, is destroyed.
 *
 * @param[in,out] B
 *          On entry, the symmetric (or Hermitian) positive definite
 *          matrix B.
 *          If uplo = PlasmaUpper, the leading N-by-N upper triangular
 *          part of B contains the upper triangular part of the matrix
 *          B, and the strictly lower triangular part of B is not
 *          referenced.
 *          If uplo = PlasmaLower, the leading N-by-N lower triangular
 *          part of B contains the lower triangular part of the matrix
 *          B, and the strictly upper triangular part of B is not
 *          referenced.
 *          On exit, if return value <= N, the part of B containing
 *          the matrix is overwritten by the triangular factor U or L
 *          from the Cholesky factorization B = U**T*U or B = L*L**T.
 *
 * @param[in,out] T
 *          On entry, descriptor as return by
 *          PLASMA_Alloc_Workspace_ssygvd
 *          On exit, contains auxiliary factorization data.
 *
 * @param[out] W
 *          On exit, if info = 0, the eigenvalues.
 *
 * @param[out] Q
 *          On exit, if jobz = PlasmaVec and info = 0, the eigenvectors.
 *
 *******************************************************************************
 *
 * @return
 *          \retval PLASMA_SUCCESS successful exit
 *          \retval <0 if -i, the i-th argument had an illegal value
 *          \retval <=N if INFO = i, plasma_ssygvd failed to converge; i
 *               off-diagonal elements of an intermediate tridiagonal
 *               form did not converge to zero.
 *          \retval >N if INFO = N + i, for 1 <= i <= N, then the leading
 *                     minor of order i of B is not positive definite.
 *                     The factorization of B could not be completed and
 *                     no eigenvalues or eigenvectors were computed.
 *
 *******************************************************************************
 *
 * @sa PLASMA_ssygvd
 * @sa PLASMA_ssygvd_Tile_Async
 * @sa PLASMA_chegvd_Tile
 * @sa PLASMA_dsygvd_Tile
 * @sa PLASMA_ssygvd_Tile
 *
 ******************************************************************************/
int PLASMA_ssygvd_Tile(PLASMA_enum itype, PLASMA_enum jobz, PLASMA_enum uplo,
                      PLASMA_desc *A, PLASMA_desc *B, float *W,
                      PLASMA_desc *T, PLASMA_desc *Q)
{
    plasma_context_t *plasma;
    PLASMA_sequence *sequence = NULL;
    PLASMA_request request = PLASMA_REQUEST_INITIALIZER;
    int status;

    plasma = plasma_context_self();
    if (plasma == NULL) {
        plasma_fatal_error("PLASMA_ssygvd_Tile", "PLASMA not initialized");
        return PLASMA_ERR_NOT_INITIALIZED;
    }
    plasma_sequence_create(plasma, &sequence);
    PLASMA_ssygvd_Tile_Async(itype, jobz, uplo, A, B, W, T, Q, sequence, &request);
    plasma_dynamic_sync();
    status = sequence->status;
    plasma_sequence_destroy(plasma, sequence);
    return status;
}

/***************************************************************************//**
 *
 * @ingroup float_Tile_Async
 *
 *  PLASMA_ssygvd_Tile - Computes all eigenvalues and, optionally,
 *  eigenvectors of a complex generalized Hermitian-definite
 *  eigenproblem of the form: A*x=(lambda)*B*x, A*Bx=(lambda)*x, or
 *  B*A*x=(lambda)*x.
 *  Here A and B are assumed to be Hermitian and B is also
 *  positive definite.
 *
 *  Non-blocking equivalent of PLASMA_ssygvd_Tile().
 *  May return before the computation is finished.
 *  Allows for pipelining of operations at runtime.
 *
 *******************************************************************************
 *
 * @param[in] sequence
 *          Identifies the sequence of function calls that this call belongs to
 *          (for completion checks and exception handling purposes).
 *
 * @param[out] request
 *          Identifies this function call (for exception handling purposes).
 *
 *******************************************************************************
 *
 * @sa PLASMA_ssygvd
 * @sa PLASMA_ssygvd_Tile
 * @sa PLASMA_chegvd_Tile_Async
 * @sa PLASMA_dsygvd_Tile_Async
 * @sa PLASMA_ssygvd_Tile_Async
 *
 ******************************************************************************/
int PLASMA_ssygvd_Tile_Async(PLASMA_enum itype, PLASMA_enum jobz, PLASMA_enum uplo,
                             PLASMA_desc *A, PLASMA_desc *B, float *W,
                             PLASMA_desc *T, PLASMA_desc *Q,
PLASMA_sequence *sequence, PLASMA_request *request)
{
    int status;
    PLASMA_desc descA;
    PLASMA_desc descB;
    PLASMA_desc descQ;
    PLASMA_enum trans;
    plasma_context_t *plasma;

    plasma = plasma_context_self();
    if (plasma == NULL) {
        plasma_fatal_error("PLASMA_ssygvd_Tile_Async", "PLASMA not initialized");
        return PLASMA_ERR_NOT_INITIALIZED;
    }
    if (sequence == NULL) {
        plasma_fatal_error("PLASMA_ssygvd_Tile_Async", "NULL sequence");
        return PLASMA_ERR_UNALLOCATED;
    }
    if (request == NULL) {
        plasma_fatal_error("PLASMA_ssygvd_Tile_Async", "NULL request");
        return PLASMA_ERR_UNALLOCATED;
    }
    /* Check sequence status */
    if (sequence->status == PLASMA_SUCCESS)
        request->status = PLASMA_SUCCESS;
    else
        return plasma_request_fail(sequence, request, PLASMA_ERR_SEQUENCE_FLUSHED);

    /* Check descriptors for correctness */
    if (plasma_desc_check(A) != PLASMA_SUCCESS) {
        plasma_error("PLASMA_ssygvd_Tile_Async", "invalid first descriptor");
        return plasma_request_fail(sequence, request, PLASMA_ERR_ILLEGAL_VALUE);
    } else {
        descA = *A;
    }
    if (plasma_desc_check(B) != PLASMA_SUCCESS) {
        plasma_error("PLASMA_ssygvd_Tile_Async", "invalid second descriptor");
        return plasma_request_fail(sequence, request, PLASMA_ERR_ILLEGAL_VALUE);
    } else {
        descB = *B;
    }
    if (plasma_desc_check(T) != PLASMA_SUCCESS) {
        plasma_error("PLASMA_ssygvd_Tile_Async", "invalid T descriptor");
        return plasma_request_fail(sequence, request, PLASMA_ERR_ILLEGAL_VALUE);
    }
    if (plasma_desc_check(Q) != PLASMA_SUCCESS) {
        plasma_error("PLASMA_ssygvd_Tile_Async", "invalid Q descriptor");
        return plasma_request_fail(sequence, request, PLASMA_ERR_ILLEGAL_VALUE);
    } else {
        descQ = *Q;
    }

    /* Check input arguments */
    if (itype != 1 && itype != 2 && itype != 3) {
        plasma_error("PLASMA_ssygvd_Tile_Async", "Illegal value of itype");
        return plasma_request_fail(sequence, request, PLASMA_ERR_ILLEGAL_VALUE);
    }
    if (jobz != PlasmaNoVec && jobz != PlasmaVec) {
        plasma_error("PLASMA_ssygvd_Tile_Async", "illegal value of jobz");
        return plasma_request_fail(sequence, request, PLASMA_ERR_ILLEGAL_VALUE);
    }
    if (uplo != PlasmaLower && uplo != PlasmaUpper) {
        plasma_error("PLASMA_ssygvd_Tile_Async", "illegal value of uplo");
        return plasma_request_fail(sequence, request, PLASMA_ERR_ILLEGAL_VALUE);
    }
    if (descA.nb != descA.mb) {
        plasma_error("PLASMA_ssygvd_Tile_Async", "only square tiles supported");
        return plasma_request_fail(sequence, request, PLASMA_ERR_ILLEGAL_VALUE);
    }
    if ( (jobz == PlasmaVec) && (descA.nb != descA.mb) ) {
        plasma_error("PLASMA_ssygv_Tile_Async", "only square tiles supported");
        return plasma_request_fail(sequence, request, PLASMA_ERR_ILLEGAL_VALUE);
    }

    /*
     * Form a Cholesky factorization of B.
     */
    plasma_parallel_call_4(plasma_pspotrf,
        PLASMA_enum,      uplo,
        PLASMA_desc,      descB,
        PLASMA_sequence*, sequence,
        PLASMA_request*,  request);

    status = sequence->status;
    if (status != 0){
       status = descA.n + status;
       return status;
    }

    /*
     * Transform problem to standard eigenvalue problem and solve.
     */
    plasma_dynamic_call_6(plasma_pssygst,
        PLASMA_enum,      itype,
        PLASMA_enum,      uplo,
        PLASMA_desc,      descA,
        PLASMA_desc,      descB,
        PLASMA_sequence*, sequence,
        PLASMA_request*,  request);

    /*
     * Call Eigensolver
     */
    status = PLASMA_ssyevd_Tile_Async( jobz, uplo,
                                       A, W, T,
                                       (float*)(Q->mat), Q->lm,
                                       sequence, request );

    /*
     * Backtransform eigenvectors to the original problem.
     */
    if( jobz == PlasmaVec ) {
        if ( status != PLASMA_SUCCESS ) {
            plasma_error("PLASMA_sheegvd", "ssygvd");
            return status;
        }

        /*
         * Convert back the lapack layout used in ssyev to tile layout
         */
        plasma_siplap2tile( descQ, (float*)Q->mat, descQ.mb, descQ.nb, descQ.lm, descQ.ln,
                            0, 0, descQ.m, descQ.n, sequence, request );

        if( (itype == 1) || (itype == 2) ) {
            /*
             * For A*x=(lambda)*B*x and A*B*x=(lambda)*x;
             * backtransform eigenvectors: x = inv(L)**T *y or inv(U)*y
             */
             if( uplo == PlasmaUpper )
                 trans = PlasmaNoTrans;
             else
                 trans = PlasmaTrans;

             plasma_parallel_call_9(plasma_pstrsm,
                 PLASMA_enum,        PlasmaLeft,
                 PLASMA_enum,        uplo,
                 PLASMA_enum,        trans,
                 PLASMA_enum,        PlasmaNonUnit,
                 float, 1.0,
                 PLASMA_desc,        descB,
                 PLASMA_desc,        descQ,
                 PLASMA_sequence*,   sequence,
                 PLASMA_request*,    request);
        }
        else if( itype == 3 ) {
            /*
             * For B*A*x=(lambda)*x;
             * backtransform eigenvectors: x = L*y or U**T *y
             */
             if( uplo == PlasmaUpper )
                 trans = PlasmaTrans;
             else
                 trans = PlasmaNoTrans;

             plasma_dynamic_call_9(plasma_pstrmm,
                 PLASMA_enum,        PlasmaLeft,
                 PLASMA_enum,        uplo,
                 PLASMA_enum,        trans,
                 PLASMA_enum,        PlasmaNonUnit,
                 float, 1.0,
                 PLASMA_desc,        descB,
                 PLASMA_desc,        descQ,
                 PLASMA_sequence*,   sequence,
                 PLASMA_request*,    request);
        }
    }

    return PLASMA_SUCCESS;
}


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

#include <config.h>
#include <math.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_multiroots.h>

#include "test_funcs.h"

/* For information on testing see the following paper,

   J.J More, B.S. Garbow, K.E. Hillstrom, "Testing Unconstrained
   Optimization Software", ACM Transactions on Mathematical Software,
   Vol 7, No 1, (1981) p 17-41

   */

/* Rosenbrock Function */

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

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

int
rosenbrock_f (const gsl_vector * x, void *params, gsl_vector * f)
{
  double x0 = gsl_vector_get (x, 0);
  double x1 = gsl_vector_get (x, 1);

  double y0 = 1 - x0;
  double y1 = 10 * (x1 - x0 * x0);

  gsl_vector_set (f, 0, y0);
  gsl_vector_set (f, 1, y1);

  params = 0;			/* avoid warning about unused parameters */

  return GSL_SUCCESS;
}

int
rosenbrock_df (const gsl_vector * x, void *params, gsl_matrix * df)
{
  double x0 = gsl_vector_get (x, 0);

  double df00 = -1;
  double df01 = 0;
  double df10 = -20 * x0;
  double df11 = 10;

  gsl_matrix_set (df, 0, 0, df00);
  gsl_matrix_set (df, 0, 1, df01);
  gsl_matrix_set (df, 1, 0, df10);
  gsl_matrix_set (df, 1, 1, df11);

  params = 0;			/* avoid warning about unused parameters */

  return GSL_SUCCESS;
}

int
rosenbrock_fdf (const gsl_vector * x, void *params,
		gsl_vector * f, gsl_matrix * df)
{
  rosenbrock_f (x, params, f);
  rosenbrock_df (x, params, df);

  return GSL_SUCCESS;
}

/* Freudenstein and Roth function */

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

void
roth_initpt (gsl_vector * x)
{
  gsl_vector_set (x, 0, 4.5);  /* changed from the value in the paper */
  gsl_vector_set (x, 1, 3.5);  /* otherwise the problem is too hard */
}

int
roth_f (const gsl_vector * x, void *params, gsl_vector * f)
{
  double x0 = gsl_vector_get (x, 0);
  double x1 = gsl_vector_get (x, 1);

  double y0 = -13.0 + x0 + ((5.0 - x1)*x1 - 2.0)*x1;
  double y1 = -29.0 + x0 + ((x1 + 1.0)*x1 - 14.0)*x1;

  gsl_vector_set (f, 0, y0);
  gsl_vector_set (f, 1, y1);

  params = 0;			/* avoid warning about unused parameters */

  return GSL_SUCCESS;
}

int
roth_df (const gsl_vector * x, void *params, gsl_matrix * df)
{
  double x1 = gsl_vector_get (x, 1);

  double df00 = 1;
  double df01 = -3 * x1 * x1 + 10 * x1 - 2;
  double df10 = 1;
  double df11 = 3 * x1 * x1 + 2 * x1 - 14;

  gsl_matrix_set (df, 0, 0, df00);
  gsl_matrix_set (df, 0, 1, df01);
  gsl_matrix_set (df, 1, 0, df10);
  gsl_matrix_set (df, 1, 1, df11);

  params = 0;			/* avoid warning about unused parameters */

  return GSL_SUCCESS;
}

int
roth_fdf (const gsl_vector * x, void *params,
		gsl_vector * f, gsl_matrix * df)
{
  roth_f (x, params, f);
  roth_df (x, params, df);

  return GSL_SUCCESS;
}

/* Powell badly scaled function */

gsl_multiroot_function_fdf powellscal =
{&powellscal_f,
 &powellscal_df,
 &powellscal_fdf,
 2, 0};

void
powellscal_initpt (gsl_vector * x)
{
  gsl_vector_set (x, 0, 0.0);
  gsl_vector_set (x, 1, 1.0);
}

int
powellscal_f (const gsl_vector * x, void *params, gsl_vector * f)
{
  double x0 = gsl_vector_get (x, 0);
  double x1 = gsl_vector_get (x, 1);

  double y0 = 10000.0 * x0 * x1 - 1.0;
  double y1 = exp (-x0) + exp (-x1) - 1.0001;

  gsl_vector_set (f, 0, y0);
  gsl_vector_set (f, 1, y1);

  params = 0;			/* avoid warning about unused parameters */

  return GSL_SUCCESS;
}

int
powellscal_df (const gsl_vector * x, void *params, gsl_matrix * df)
{
  double x0 = gsl_vector_get (x, 0);
  double x1 = gsl_vector_get (x, 1);

  double df00 = 10000.0 * x1, df01 = 10000.0 * x0;
  double df10 = -exp (-x0), df11 = -exp (-x1);

  gsl_matrix_set (df, 0, 0, df00);
  gsl_matrix_set (df, 0, 1, df01);
  gsl_matrix_set (df, 1, 0, df10);
  gsl_matrix_set (df, 1, 1, df11);

  params = 0;			/* avoid warning about unused parameters */

  return GSL_SUCCESS;
}

int
powellscal_fdf (const gsl_vector * x, void *params,
		  gsl_vector * f, gsl_matrix * df)
{
  powellscal_f (x, params, f);
  powellscal_df (x, params, df);

  return GSL_SUCCESS;
}

/* Brown badly scaled function */

gsl_multiroot_function_fdf brownscal =
{&brownscal_f,
 &brownscal_df,
 &brownscal_fdf,
 2, 0};

void
brownscal_initpt (gsl_vector * x)
{
  gsl_vector_set (x, 0, 1.0);
  gsl_vector_set (x, 1, 1.0);
}

int
brownscal_f (const gsl_vector * x, void *params, gsl_vector * f)
{
  double x0 = gsl_vector_get (x, 0);
  double x1 = gsl_vector_get (x, 1);

  double y0 = x0 - 1e6;
  double y1 = x0 * x1 - 2;

  gsl_vector_set (f, 0, y0);
  gsl_vector_set (f, 1, y1);

  params = 0;			/* avoid warning about unused parameters */

  return GSL_SUCCESS;
}

int
brownscal_df (const gsl_vector * x, void *params, gsl_matrix * df)
{
  double x0 = gsl_vector_get (x, 0);
  double x1 = gsl_vector_get (x, 1);

  double df00 = 1.0, df01 = 0.0;
  double df10 = x1, df11 = x0;

  gsl_matrix_set (df, 0, 0, df00);
  gsl_matrix_set (df, 0, 1, df01);
  gsl_matrix_set (df, 1, 0, df10);
  gsl_matrix_set (df, 1, 1, df11);

  params = 0;			/* avoid warning about unused parameters */

  return GSL_SUCCESS;
}

int
brownscal_fdf (const gsl_vector * x, void *params,
		  gsl_vector * f, gsl_matrix * df)
{
  brownscal_f (x, params, f);
  brownscal_df (x, params, df);

  return GSL_SUCCESS;
}

/* Powell Singular Function */

gsl_multiroot_function_fdf powellsing =
{&powellsing_f,
 &powellsing_df,
 &powellsing_fdf,
 4, 0};

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

int
powellsing_f (const gsl_vector * x, void *params, gsl_vector * f)
{
  double x0 = gsl_vector_get (x, 0);
  double x1 = gsl_vector_get (x, 1);
  double x2 = gsl_vector_get (x, 2);
  double x3 = gsl_vector_get (x, 3);

  double y0 = x0 + 10 * x1;
  double y1 = sqrt (5.0) * (x2 - x3);
  double y2 = pow (x1 - 2 * x2, 2.0);
  double y3 = sqrt (10.0) * pow (x0 - x3, 2.0);

  gsl_vector_set (f, 0, y0);
  gsl_vector_set (f, 1, y1);
  gsl_vector_set (f, 2, y2);
  gsl_vector_set (f, 3, y3);

  params = 0;			/* avoid warning about unused parameters */

  return GSL_SUCCESS;
}

int
powellsing_df (const gsl_vector * x, void *params, gsl_matrix * df)
{
  double x0 = gsl_vector_get (x, 0);
  double x1 = gsl_vector_get (x, 1);
  double x2 = gsl_vector_get (x, 2);
  double x3 = gsl_vector_get (x, 3);

  double df00 = 1, df01 = 10, df02 = 0, df03 = 0;
  double df10 = 0, df11 = 0, df12 = sqrt (5.0), df13 = -df12;
  double df20 = 0, df21 = 2 * (x1 - 2 * x2), df22 = -2 * df21, df23 = 0;
  double df30 = 2 * sqrt (10.0) * (x0 - x3), df31 = 0, df32 = 0, df33 = -df30;

  gsl_matrix_set (df, 0, 0, df00);
  gsl_matrix_set (df, 0, 1, df01);
  gsl_matrix_set (df, 0, 2, df02);
  gsl_matrix_set (df, 0, 3, df03);

  gsl_matrix_set (df, 1, 0, df10);
  gsl_matrix_set (df, 1, 1, df11);
  gsl_matrix_set (df, 1, 2, df12);
  gsl_matrix_set (df, 1, 3, df13);

  gsl_matrix_set (df, 2, 0, df20);
  gsl_matrix_set (df, 2, 1, df21);
  gsl_matrix_set (df, 2, 2, df22);
  gsl_matrix_set (df, 2, 3, df23);

  gsl_matrix_set (df, 3, 0, df30);
  gsl_matrix_set (df, 3, 1, df31);
  gsl_matrix_set (df, 3, 2, df32);
  gsl_matrix_set (df, 3, 3, df33);

  params = 0;			/* avoid warning about unused parameters */

  return GSL_SUCCESS;
}

int
powellsing_fdf (const gsl_vector * x, void *params,
		    gsl_vector * f, gsl_matrix * df)
{
  powellsing_f (x, params, f);
  powellsing_df (x, params, df);

  return GSL_SUCCESS;
}

/* Wood function */

gsl_multiroot_function_fdf wood =
{&wood_f,
 &wood_df,
 &wood_fdf,
 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);
}

int
wood_f (const gsl_vector * x, void *params, gsl_vector * f)
{
  double x0 = gsl_vector_get (x, 0);
  double x1 = gsl_vector_get (x, 1);
  double x2 = gsl_vector_get (x, 2);
  double x3 = gsl_vector_get (x, 3);

  double t1 = x1 - x0 * x0;
  double t2 = x3 - x2 * x2;

  double y0 = -200.0 * x0 * t1 - (1 - x0);
  double y1 = 200.0 * t1 + 20.2 * (x1 - 1) + 19.8 * (x3 - 1);
  double y2 = -180.0 * x2 * t2 - (1 - x2);
  double y3 = 180.0 * t2 + 20.2 * (x3 - 1) + 19.8 * (x1 - 1);

  gsl_vector_set (f, 0, y0);
  gsl_vector_set (f, 1, y1);
  gsl_vector_set (f, 2, y2);
  gsl_vector_set (f, 3, y3);

  params = 0;			/* avoid warning about unused parameters */

  return GSL_SUCCESS;
}

int
wood_df (const gsl_vector * x, void *params, gsl_matrix * df)
{
  double x0 = gsl_vector_get (x, 0);
  double x1 = gsl_vector_get (x, 1);
  double x2 = gsl_vector_get (x, 2);
  double x3 = gsl_vector_get (x, 3);

  double t1 = x1 - 3 * x0 * x0;
  double t2 = x3 - 3 * x2 * x2;

  double df00 = -200.0 * t1 + 1, df01 = -200.0 * x0, df02 = 0, df03 = 0;
  double df10 = -400.0*x0, df11 = 200.0 + 20.2, df12 = 0, df13 = 19.8;
  double df20 = 0, df21 = 0, df22 = -180.0 * t2 + 1, df23 = -180.0 * x2;
  double df30 = 0, df31 = 19.8, df32 = -2 * 180 * x2, df33 = 180.0 + 20.2;

  gsl_matrix_set (df, 0, 0, df00);
  gsl_matrix_set (df, 0, 1, df01);
  gsl_matrix_set (df, 0, 2, df02);
  gsl_matrix_set (df, 0, 3, df03);

  gsl_matrix_set (df, 1, 0, df10);
  gsl_matrix_set (df, 1, 1, df11);
  gsl_matrix_set (df, 1, 2, df12);
  gsl_matrix_set (df, 1, 3, df13);

  gsl_matrix_set (df, 2, 0, df20);
  gsl_matrix_set (df, 2, 1, df21);
  gsl_matrix_set (df, 2, 2, df22);
  gsl_matrix_set (df, 2, 3, df23);

  gsl_matrix_set (df, 3, 0, df30);
  gsl_matrix_set (df, 3, 1, df31);
  gsl_matrix_set (df, 3, 2, df32);
  gsl_matrix_set (df, 3, 3, df33);

  params = 0;			/* avoid warning about unused parameters */

  return GSL_SUCCESS;
}

int
wood_fdf (const gsl_vector * x, void *params,
		    gsl_vector * f, gsl_matrix * df)
{
  wood_f (x, params, f);
  wood_df (x, params, df);

  return GSL_SUCCESS;
}

/* Helical Valley Function */

gsl_multiroot_function_fdf helical =
{&helical_f,
 &helical_df,
 &helical_fdf,
 3, 0};

void
helical_initpt (gsl_vector * x)
{
  gsl_vector_set (x, 0, -1.0);
  gsl_vector_set (x, 1, 0.0);
  gsl_vector_set (x, 2, 0.0);
}

int
helical_f (const gsl_vector * x, void *params, gsl_vector * f)
{
  double x0 = gsl_vector_get (x, 0);
  double x1 = gsl_vector_get (x, 1);
  double x2 = gsl_vector_get (x, 2);

  double t1, t2;
  double y0, y1, y2;

  if (x0 > 0)
    {
      t1 = atan(x1/x0) / (2.0 * M_PI);
    }
  else if (x0 < 0)
    {
      t1 = 0.5 + atan(x1/x0) / (2.0 * M_PI);
    }
  else
    {
      t1 = 0.25 * (x1 > 0 ? +1 : -1);
    }

  t2 = sqrt(x0*x0 + x1*x1) ;

  y0 = 10 * (x2 - 10 * t1);
  y1 = 10 * (t2 - 1);
  y2 = x2 ;

  gsl_vector_set (f, 0, y0);
  gsl_vector_set (f, 1, y1);
  gsl_vector_set (f, 2, y2);

  params = 0;			/* avoid warning about unused parameters */

  return GSL_SUCCESS;
}

int
helical_df (const gsl_vector * x, void *params, gsl_matrix * df)
{
  double x0 = gsl_vector_get (x, 0);
  double x1 = gsl_vector_get (x, 1);

  double t = x0 * x0 + x1 * x1 ;
  double t1 = 2 * M_PI * t ;
  double t2 = sqrt(t) ;

  double df00 = 100*x1/t1, df01 = -100.0 * x0/t1, df02 = 10.0;
  double df10 = 10*x0/t2, df11 = 10*x1/t2, df12 = 0;
  double df20 = 0, df21 = 0, df22 = 1.0;

  gsl_matrix_set (df, 0, 0, df00);
  gsl_matrix_set (df, 0, 1, df01);
  gsl_matrix_set (df, 0, 2, df02);

  gsl_matrix_set (df, 1, 0, df10);
  gsl_matrix_set (df, 1, 1, df11);
  gsl_matrix_set (df, 1, 2, df12);

  gsl_matrix_set (df, 2, 0, df20);
  gsl_matrix_set (df, 2, 1, df21);
  gsl_matrix_set (df, 2, 2, df22);

  params = 0;			/* avoid warning about unused parameters */

  return GSL_SUCCESS;
}

int
helical_fdf (const gsl_vector * x, void *params,
		    gsl_vector * f, gsl_matrix * df)
{
  helical_f (x, params, f);
  helical_df (x, params, df);

  return GSL_SUCCESS;
}

/* Discrete Boundary Value Function */

#define N 10

gsl_multiroot_function_fdf dbv =
{&dbv_f,
 &dbv_df,
 &dbv_fdf,
 N, 0};

void
dbv_initpt (gsl_vector * x)
{
  size_t i;
  double h = 1.0 / (N + 1.0);

  for (i = 0; i < N; i++)
    {
      double t = (i + 1) * h;
      double z = t * (t - 1);
      gsl_vector_set (x, i, z);
    }
}

int
dbv_f (const gsl_vector * x, void *params, gsl_vector * f)
{
  size_t i;

  double h = 1.0 / (N + 1.0);

  for (i = 0; i < N; i++)
    {
      double z, ti = (i + 1) * h;
      double xi = 0, xim1 = 0, xip1 = 0;

      xi = gsl_vector_get (x, i);

      if (i > 0)
        xim1 = gsl_vector_get (x, i - 1);

      if (i < N - 1)
        xip1 = gsl_vector_get (x, i + 1);

      z = 2 * xi - xim1 - xip1 + h * h * pow(xi + ti + 1, 3.0) / 2.0;

      gsl_vector_set (f, i, z);

    }

  params = 0;			/* avoid warning about unused parameters */

  return GSL_SUCCESS;
}

int
dbv_df (const gsl_vector * x, void *params, gsl_matrix * df)
{
  size_t i, j;

  double h = 1.0 / (N + 1.0);

  for (i = 0; i < N; i++)
    for (j = 0; j < N; j++)
      gsl_matrix_set (df, i, j, 0.0);

  for (i = 0; i < N; i++)
    {
      double dz_dxi, ti = (i + 1) * h;

      double xi = gsl_vector_get (x, i);

      dz_dxi = 2.0 + (3.0 / 2.0) * h * h * pow(xi + ti + 1, 2.0) ;

      gsl_matrix_set (df, i, i, dz_dxi);

      if (i > 0)
        gsl_matrix_set (df, i, i-1, -1.0);

      if (i < N - 1)
        gsl_matrix_set (df, i, i+1, -1.0);

    }

  params = 0;			/* avoid warning about unused parameters */

  return GSL_SUCCESS;
}

int
dbv_fdf (const gsl_vector * x, void *params,
		    gsl_vector * f, gsl_matrix * df)
{
  dbv_f (x, params, f);
  dbv_df (x, params, df);

  return GSL_SUCCESS;
}

/* Trigonometric Function */

gsl_multiroot_function_fdf trig =
{&trig_f,
 &trig_df,
 &trig_fdf,
 N, 0};

void
trig_initpt (gsl_vector * x)
{
  size_t i;

  for (i = 0; i < N; i++)       /* choose an initial point which converges */
    {
      gsl_vector_set (x, i, 0.05);
    }
}

int
trig_f (const gsl_vector * x, void *params, gsl_vector * f)
{
  size_t i;
  double sum = 0;

  for (i = 0; i < N; i++)
    {
      sum += cos(gsl_vector_get(x,i));
    }

  for (i = 0; i < N; i++)
    {
      double xi = gsl_vector_get (x,i);
      double z = N - sum + (i + 1) * (1 - cos(xi)) - sin(xi);

      gsl_vector_set (f, i, z);
    }

  params = 0;			/* avoid warning about unused parameters */

  return GSL_SUCCESS;
}

int
trig_df (const gsl_vector * x, void *params, gsl_matrix * df)
{
  size_t i, j;

  for (i = 0; i < N; i++)
    {
      for (j = 0; j < N; j++)
        {
          double dz;
          double xi = gsl_vector_get(x, i);
          double xj = gsl_vector_get(x, j);

          if (j == i)
            dz = sin(xi) + (i + 1) * sin(xi) - cos(xi);
          else
            dz = sin(xj);

          gsl_matrix_set(df, i, j, dz);
        }
    }

  params = 0;			/* avoid warning about unused parameters */

  return GSL_SUCCESS;
}

int
trig_fdf (const gsl_vector * x, void *params,
		    gsl_vector * f, gsl_matrix * df)
{
  trig_f (x, params, f);
  trig_df (x, params, df);

  return GSL_SUCCESS;
}