From 197e88c3b56d28840217010a2871c6ea3a4dd1a4 Mon Sep 17 00:00:00 2001 From: Alberto Ruiz Date: Wed, 21 May 2014 10:30:55 +0200 Subject: update dependencies, move examples etc --- packages/gsl/src/Numeric/GSL/gsl-aux.c | 945 +++++++++++++++++++++++++++++++++ 1 file changed, 945 insertions(+) create mode 100644 packages/gsl/src/Numeric/GSL/gsl-aux.c (limited to 'packages/gsl/src/Numeric/GSL/gsl-aux.c') diff --git a/packages/gsl/src/Numeric/GSL/gsl-aux.c b/packages/gsl/src/Numeric/GSL/gsl-aux.c new file mode 100644 index 0000000..ffc5c20 --- /dev/null +++ b/packages/gsl/src/Numeric/GSL/gsl-aux.c @@ -0,0 +1,945 @@ +#include + +#define RVEC(A) int A##n, double*A##p +#define RMAT(A) int A##r, int A##c, double* A##p +#define KRVEC(A) int A##n, const double*A##p +#define KRMAT(A) int A##r, int A##c, const double* A##p + +#define CVEC(A) int A##n, gsl_complex*A##p +#define CMAT(A) int A##r, int A##c, gsl_complex* A##p +#define KCVEC(A) int A##n, const gsl_complex*A##p +#define KCMAT(A) int A##r, int A##c, const gsl_complex* A##p + +#define FVEC(A) int A##n, float*A##p +#define FMAT(A) int A##r, int A##c, float* A##p +#define KFVEC(A) int A##n, const float*A##p +#define KFMAT(A) int A##r, int A##c, const float* A##p + +#define QVEC(A) int A##n, gsl_complex_float*A##p +#define QMAT(A) int A##r, int A##c, gsl_complex_float* A##p +#define KQVEC(A) int A##n, const gsl_complex_float*A##p +#define KQMAT(A) int A##r, int A##c, const gsl_complex_float* A##p + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#define MACRO(B) do {B} while (0) +#define ERROR(CODE) MACRO(return CODE;) +#define REQUIRES(COND, CODE) MACRO(if(!(COND)) {ERROR(CODE);}) +#define OK return 0; + +#define MIN(A,B) ((A)<(B)?(A):(B)) +#define MAX(A,B) ((A)>(B)?(A):(B)) + +#ifdef DBG +#define DEBUGMSG(M) printf("*** calling aux C function: %s\n",M); +#else +#define DEBUGMSG(M) +#endif + +#define CHECK(RES,CODE) MACRO(if(RES) return CODE;) + +#ifdef DBG +#define DEBUGMAT(MSG,X) printf(MSG" = \n"); gsl_matrix_fprintf(stdout,X,"%f"); printf("\n"); +#else +#define DEBUGMAT(MSG,X) +#endif + +#ifdef DBG +#define DEBUGVEC(MSG,X) printf(MSG" = \n"); gsl_vector_fprintf(stdout,X,"%f"); printf("\n"); +#else +#define DEBUGVEC(MSG,X) +#endif + +#define DVVIEW(A) gsl_vector_view A = gsl_vector_view_array(A##p,A##n) +#define DMVIEW(A) gsl_matrix_view A = gsl_matrix_view_array(A##p,A##r,A##c) +#define CVVIEW(A) gsl_vector_complex_view A = gsl_vector_complex_view_array((double*)A##p,A##n) +#define CMVIEW(A) gsl_matrix_complex_view A = gsl_matrix_complex_view_array((double*)A##p,A##r,A##c) +#define KDVVIEW(A) gsl_vector_const_view A = gsl_vector_const_view_array(A##p,A##n) +#define KDMVIEW(A) gsl_matrix_const_view A = gsl_matrix_const_view_array(A##p,A##r,A##c) +#define KCVVIEW(A) gsl_vector_complex_const_view A = gsl_vector_complex_const_view_array((double*)A##p,A##n) +#define KCMVIEW(A) gsl_matrix_complex_const_view A = gsl_matrix_complex_const_view_array((double*)A##p,A##r,A##c) + +#define FVVIEW(A) gsl_vector_float_view A = gsl_vector_float_view_array(A##p,A##n) +#define FMVIEW(A) gsl_matrix_float_view A = gsl_matrix_float_view_array(A##p,A##r,A##c) +#define QVVIEW(A) gsl_vector_complex_float_view A = gsl_vector_float_complex_view_array((float*)A##p,A##n) +#define QMVIEW(A) gsl_matrix_complex_float_view A = gsl_matrix_float_complex_view_array((float*)A##p,A##r,A##c) +#define KFVVIEW(A) gsl_vector_float_const_view A = gsl_vector_float_const_view_array(A##p,A##n) +#define KFMVIEW(A) gsl_matrix_float_const_view A = gsl_matrix_float_const_view_array(A##p,A##r,A##c) +#define KQVVIEW(A) gsl_vector_complex_float_const_view A = gsl_vector_complex_float_const_view_array((float*)A##p,A##n) +#define KQMVIEW(A) gsl_matrix_complex_float_const_view A = gsl_matrix_complex_float_const_view_array((float*)A##p,A##r,A##c) + +#define V(a) (&a.vector) +#define M(a) (&a.matrix) + +#define GCVEC(A) int A##n, gsl_complex*A##p +#define KGCVEC(A) int A##n, const gsl_complex*A##p + +#define GQVEC(A) int A##n, gsl_complex_float*A##p +#define KGQVEC(A) int A##n, const gsl_complex_float*A##p + +#define BAD_SIZE 2000 +#define BAD_CODE 2001 +#define MEM 2002 +#define BAD_FILE 2003 + + +void no_abort_on_error() { + gsl_set_error_handler_off(); +} + + + +int fft(int code, KCVEC(X), CVEC(R)) { + REQUIRES(Xn == Rn,BAD_SIZE); + DEBUGMSG("fft"); + int s = Xn; + gsl_fft_complex_wavetable * wavetable = gsl_fft_complex_wavetable_alloc (s); + gsl_fft_complex_workspace * workspace = gsl_fft_complex_workspace_alloc (s); + gsl_vector_const_view X = gsl_vector_const_view_array((double*)Xp, 2*Xn); + gsl_vector_view R = gsl_vector_view_array((double*)Rp, 2*Rn); + gsl_blas_dcopy(&X.vector,&R.vector); + if(code==0) { + gsl_fft_complex_forward ((double*)Rp, 1, s, wavetable, workspace); + } else { + gsl_fft_complex_inverse ((double*)Rp, 1, s, wavetable, workspace); + } + gsl_fft_complex_wavetable_free (wavetable); + gsl_fft_complex_workspace_free (workspace); + OK +} + + +int deriv(int code, double f(double, void*), double x, double h, double * result, double * abserr) +{ + gsl_function F; + F.function = f; + F.params = 0; + + if(code==0) return gsl_deriv_central (&F, x, h, result, abserr); + + if(code==1) return gsl_deriv_forward (&F, x, h, result, abserr); + + if(code==2) return gsl_deriv_backward (&F, x, h, result, abserr); + + return 0; +} + + +int integrate_qng(double f(double, void*), double a, double b, double aprec, double prec, + double *result, double*error) { + DEBUGMSG("integrate_qng"); + gsl_function F; + F.function = f; + F.params = NULL; + size_t neval; + int res = gsl_integration_qng (&F, a,b, aprec, prec, result, error, &neval); + CHECK(res,res); + OK +} + +int integrate_qags(double f(double,void*), double a, double b, double aprec, double prec, int w, + double *result, double* error) { + DEBUGMSG("integrate_qags"); + gsl_integration_workspace * wk = gsl_integration_workspace_alloc (w); + gsl_function F; + F.function = f; + F.params = NULL; + int res = gsl_integration_qags (&F, a,b, aprec, prec, w,wk, result, error); + CHECK(res,res); + gsl_integration_workspace_free (wk); + OK +} + +int integrate_qagi(double f(double,void*), double aprec, double prec, int w, + double *result, double* error) { + DEBUGMSG("integrate_qagi"); + gsl_integration_workspace * wk = gsl_integration_workspace_alloc (w); + gsl_function F; + F.function = f; + F.params = NULL; + int res = gsl_integration_qagi (&F, aprec, prec, w,wk, result, error); + CHECK(res,res); + gsl_integration_workspace_free (wk); + OK +} + + +int integrate_qagiu(double f(double,void*), double a, double aprec, double prec, int w, + double *result, double* error) { + DEBUGMSG("integrate_qagiu"); + gsl_integration_workspace * wk = gsl_integration_workspace_alloc (w); + gsl_function F; + F.function = f; + F.params = NULL; + int res = gsl_integration_qagiu (&F, a, aprec, prec, w,wk, result, error); + CHECK(res,res); + gsl_integration_workspace_free (wk); + OK +} + + +int integrate_qagil(double f(double,void*), double b, double aprec, double prec, int w, + double *result, double* error) { + DEBUGMSG("integrate_qagil"); + gsl_integration_workspace * wk = gsl_integration_workspace_alloc (w); + gsl_function F; + F.function = f; + F.params = NULL; + int res = gsl_integration_qagil (&F, b, aprec, prec, w,wk, result, error); + CHECK(res,res); + gsl_integration_workspace_free (wk); + OK +} + +int integrate_cquad(double f(double,void*), double a, double b, double aprec, double prec, + int w, double *result, double* error, int *neval) { + DEBUGMSG("integrate_cquad"); + gsl_integration_cquad_workspace * wk = gsl_integration_cquad_workspace_alloc (w); + gsl_function F; + F.function = f; + F.params = NULL; + size_t * sneval = NULL; + int res = gsl_integration_cquad (&F, a, b, aprec, prec, wk, result, error, sneval); + *neval = *sneval; + CHECK(res,res); + gsl_integration_cquad_workspace_free (wk); + OK +} + + +int polySolve(KRVEC(a), CVEC(z)) { + DEBUGMSG("polySolve"); + REQUIRES(an>1,BAD_SIZE); + gsl_poly_complex_workspace * w = gsl_poly_complex_workspace_alloc (an); + int res = gsl_poly_complex_solve ((double*)ap, an, w, (double*)zp); + CHECK(res,res); + gsl_poly_complex_workspace_free (w); + OK; +} + +int vector_fscanf(char*filename, RVEC(a)) { + DEBUGMSG("gsl_vector_fscanf"); + DVVIEW(a); + FILE * f = fopen(filename,"r"); + CHECK(!f,BAD_FILE); + int res = gsl_vector_fscanf(f,V(a)); + CHECK(res,res); + fclose (f); + OK +} + +int vector_fprintf(char*filename, char*fmt, RVEC(a)) { + DEBUGMSG("gsl_vector_fprintf"); + DVVIEW(a); + FILE * f = fopen(filename,"w"); + CHECK(!f,BAD_FILE); + int res = gsl_vector_fprintf(f,V(a),fmt); + CHECK(res,res); + fclose (f); + OK +} + +int vector_fread(char*filename, RVEC(a)) { + DEBUGMSG("gsl_vector_fread"); + DVVIEW(a); + FILE * f = fopen(filename,"r"); + CHECK(!f,BAD_FILE); + int res = gsl_vector_fread(f,V(a)); + CHECK(res,res); + fclose (f); + OK +} + +int vector_fwrite(char*filename, RVEC(a)) { + DEBUGMSG("gsl_vector_fwrite"); + DVVIEW(a); + FILE * f = fopen(filename,"w"); + CHECK(!f,BAD_FILE); + int res = gsl_vector_fwrite(f,V(a)); + CHECK(res,res); + fclose (f); + OK +} + +int matrix_fprintf(char*filename, char*fmt, int ro, RMAT(m)) { + DEBUGMSG("matrix_fprintf"); + FILE * f = fopen(filename,"w"); + CHECK(!f,BAD_FILE); + int i,j,sr,sc; + if (ro==1) { sr = mc; sc = 1;} else { sr = 1; sc = mr;} + #define AT(M,r,c) (M##p[(r)*sr+(c)*sc]) + for (i=0; isize,sizeof(double)); + int k; + for(k=0;ksize;k++) { + p[k] = gsl_vector_get(x,k); + } + double res = f(x->size,p); + free(p); + return res; +} + +double only_f_aux_root(double x, void *pars); +int uniMinimize(int method, double f(double), + double epsrel, int maxit, double min, + double xl, double xu, RMAT(sol)) { + REQUIRES(solr == maxit && solc == 4,BAD_SIZE); + DEBUGMSG("minimize_only_f"); + gsl_function my_func; + my_func.function = only_f_aux_root; + my_func.params = f; + size_t iter = 0; + int status; + const gsl_min_fminimizer_type *T; + gsl_min_fminimizer *s; + // Starting point + switch(method) { + case 0 : {T = gsl_min_fminimizer_goldensection; break; } + case 1 : {T = gsl_min_fminimizer_brent; break; } + case 2 : {T = gsl_min_fminimizer_quad_golden; break; } + default: ERROR(BAD_CODE); + } + s = gsl_min_fminimizer_alloc (T); + gsl_min_fminimizer_set (s, &my_func, min, xl, xu); + do { + double current_min, current_lo, current_hi; + status = gsl_min_fminimizer_iterate (s); + current_min = gsl_min_fminimizer_x_minimum (s); + current_lo = gsl_min_fminimizer_x_lower (s); + current_hi = gsl_min_fminimizer_x_upper (s); + solp[iter*solc] = iter + 1; + solp[iter*solc+1] = current_min; + solp[iter*solc+2] = current_lo; + solp[iter*solc+3] = current_hi; + iter++; + if (status) /* check if solver is stuck */ + break; + + status = + gsl_min_test_interval (current_lo, current_hi, 0, epsrel); + } + while (status == GSL_CONTINUE && iter < maxit); + int i; + for (i=iter; ifval; + solp[iter*solc+2] = size; + + int k; + for(k=0;kx,k); + } + iter++; + if (status) break; + status = gsl_multimin_test_size (size, tolsize); + } while (status == GSL_CONTINUE && iter < maxit); + int i,j; + for (i=iter; isize,sizeof(double)); + int k; + for(k=0;ksize;k++) { + p[k] = gsl_vector_get(x,k); + } + double res = fdf->f(x->size,p); + free(p); + return res; +} + + +void df_aux_min(const gsl_vector * x, void * pars, gsl_vector * g) { + Tfdf * fdf = ((Tfdf*) pars); + double* p = (double*)calloc(x->size,sizeof(double)); + double* q = (double*)calloc(g->size,sizeof(double)); + int k; + for(k=0;ksize;k++) { + p[k] = gsl_vector_get(x,k); + } + + fdf->df(x->size,p,g->size,q); + + for(k=0;ksize;k++) { + gsl_vector_set(g,k,q[k]); + } + free(p); + free(q); +} + +void fdf_aux_min(const gsl_vector * x, void * pars, double * f, gsl_vector * g) { + *f = f_aux_min(x,pars); + df_aux_min(x,pars,g); +} + + +int minimizeD(int method, double f(int, double*), int df(int, double*, int, double*), + double initstep, double minimpar, double tolgrad, int maxit, + KRVEC(xi), RMAT(sol)) { + REQUIRES(solr == maxit && solc == 2+xin,BAD_SIZE); + DEBUGMSG("minimizeWithDeriv (conjugate_fr)"); + gsl_multimin_function_fdf my_func; + // extract function from pars + my_func.f = f_aux_min; + my_func.df = df_aux_min; + my_func.fdf = fdf_aux_min; + my_func.n = xin; + Tfdf stfdf; + stfdf.f = f; + stfdf.df = df; + my_func.params = &stfdf; + size_t iter = 0; + int status; + const gsl_multimin_fdfminimizer_type *T; + gsl_multimin_fdfminimizer *s = NULL; + // Starting point + KDVVIEW(xi); + // conjugate gradient fr + switch(method) { + case 0 : {T = gsl_multimin_fdfminimizer_conjugate_fr; break; } + case 1 : {T = gsl_multimin_fdfminimizer_conjugate_pr; break; } + case 2 : {T = gsl_multimin_fdfminimizer_vector_bfgs; break; } + case 3 : {T = gsl_multimin_fdfminimizer_vector_bfgs2; break; } + case 4 : {T = gsl_multimin_fdfminimizer_steepest_descent; break; } + default: ERROR(BAD_CODE); + } + s = gsl_multimin_fdfminimizer_alloc (T, my_func.n); + gsl_multimin_fdfminimizer_set (s, &my_func, V(xi), initstep, minimpar); + do { + status = gsl_multimin_fdfminimizer_iterate (s); + solp[iter*solc+0] = iter+1; + solp[iter*solc+1] = s->f; + int k; + for(k=0;kx,k); + } + iter++; + if (status) break; + status = gsl_multimin_test_gradient (s->gradient, tolgrad); + } while (status == GSL_CONTINUE && iter < maxit); + int i,j; + for (i=iter; if)(x); +} + +double jf_aux_uni(double x, void * pars) { + uniTfjf * fjf = ((uniTfjf*) pars); + return (fjf->jf)(x); +} + +void fjf_aux_uni(double x, void * pars, double * f, double * g) { + *f = f_aux_uni(x,pars); + *g = jf_aux_uni(x,pars); +} + +int rootj(int method, double f(double), + double df(double), + double epsrel, int maxit, + double x, RMAT(sol)) { + REQUIRES(solr == maxit && solc == 2,BAD_SIZE); + DEBUGMSG("root_fjf"); + gsl_function_fdf my_func; + // extract function from pars + my_func.f = f_aux_uni; + my_func.df = jf_aux_uni; + my_func.fdf = fjf_aux_uni; + uniTfjf stfjf; + stfjf.f = f; + stfjf.jf = df; + my_func.params = &stfjf; + size_t iter = 0; + int status; + const gsl_root_fdfsolver_type *T; + gsl_root_fdfsolver *s; + // Starting point + switch(method) { + case 0 : {T = gsl_root_fdfsolver_newton;; break; } + case 1 : {T = gsl_root_fdfsolver_secant; break; } + case 2 : {T = gsl_root_fdfsolver_steffenson; break; } + default: ERROR(BAD_CODE); + } + s = gsl_root_fdfsolver_alloc (T); + + gsl_root_fdfsolver_set (s, &my_func, x); + + do { + double x0; + status = gsl_root_fdfsolver_iterate (s); + x0 = x; + x = gsl_root_fdfsolver_root(s); + solp[iter*solc+0] = iter+1; + solp[iter*solc+1] = x; + + iter++; + if (status) /* check if solver is stuck */ + break; + + status = + gsl_root_test_delta (x, x0, 0, epsrel); + } + while (status == GSL_CONTINUE && iter < maxit); + + int i; + for (i=iter; isize,sizeof(double)); + double* q = (double*)calloc(y->size,sizeof(double)); + int k; + for(k=0;ksize;k++) { + p[k] = gsl_vector_get(x,k); + } + f(x->size,p,y->size,q); + for(k=0;ksize;k++) { + gsl_vector_set(y,k,q[k]); + } + free(p); + free(q); + return 0; //hmmm +} + +int multiroot(int method, void f(int, double*, int, double*), + double epsabs, int maxit, + KRVEC(xi), RMAT(sol)) { + REQUIRES(solr == maxit && solc == 1+2*xin,BAD_SIZE); + DEBUGMSG("root_only_f"); + gsl_multiroot_function my_func; + // extract function from pars + my_func.f = only_f_aux_multiroot; + my_func.n = xin; + my_func.params = f; + size_t iter = 0; + int status; + const gsl_multiroot_fsolver_type *T; + gsl_multiroot_fsolver *s; + // Starting point + KDVVIEW(xi); + switch(method) { + case 0 : {T = gsl_multiroot_fsolver_hybrids;; break; } + case 1 : {T = gsl_multiroot_fsolver_hybrid; break; } + case 2 : {T = gsl_multiroot_fsolver_dnewton; break; } + case 3 : {T = gsl_multiroot_fsolver_broyden; break; } + default: ERROR(BAD_CODE); + } + s = gsl_multiroot_fsolver_alloc (T, my_func.n); + gsl_multiroot_fsolver_set (s, &my_func, V(xi)); + + do { + status = gsl_multiroot_fsolver_iterate (s); + + solp[iter*solc+0] = iter+1; + + int k; + for(k=0;kx,k); + } + for(k=xin;k<2*xin;k++) { + solp[iter*solc+k+1] = gsl_vector_get(s->f,k-xin); + } + + iter++; + if (status) /* check if solver is stuck */ + break; + + status = + gsl_multiroot_test_residual (s->f, epsabs); + } + while (status == GSL_CONTINUE && iter < maxit); + + int i,j; + for (i=iter; isize,sizeof(double)); + double* q = (double*)calloc(y->size,sizeof(double)); + int k; + for(k=0;ksize;k++) { + p[k] = gsl_vector_get(x,k); + } + (fjf->f)(x->size,p,y->size,q); + for(k=0;ksize;k++) { + gsl_vector_set(y,k,q[k]); + } + free(p); + free(q); + return 0; +} + +int jf_aux(const gsl_vector * x, void * pars, gsl_matrix * jac) { + Tfjf * fjf = ((Tfjf*) pars); + double* p = (double*)calloc(x->size,sizeof(double)); + double* q = (double*)calloc((jac->size1)*(jac->size2),sizeof(double)); + int i,j,k; + for(k=0;ksize;k++) { + p[k] = gsl_vector_get(x,k); + } + + (fjf->jf)(x->size,p,jac->size1,jac->size2,q); + + k=0; + for(i=0;isize1;i++) { + for(j=0;jsize2;j++){ + gsl_matrix_set(jac,i,j,q[k++]); + } + } + free(p); + free(q); + return 0; +} + +int fjf_aux(const gsl_vector * x, void * pars, gsl_vector * f, gsl_matrix * g) { + f_aux(x,pars,f); + jf_aux(x,pars,g); + return 0; +} + +int multirootj(int method, int f(int, double*, int, double*), + int jac(int, double*, int, int, double*), + double epsabs, int maxit, + KRVEC(xi), RMAT(sol)) { + REQUIRES(solr == maxit && solc == 1+2*xin,BAD_SIZE); + DEBUGMSG("root_fjf"); + gsl_multiroot_function_fdf my_func; + // extract function from pars + my_func.f = f_aux; + my_func.df = jf_aux; + my_func.fdf = fjf_aux; + my_func.n = xin; + Tfjf stfjf; + stfjf.f = f; + stfjf.jf = jac; + my_func.params = &stfjf; + size_t iter = 0; + int status; + const gsl_multiroot_fdfsolver_type *T; + gsl_multiroot_fdfsolver *s; + // Starting point + KDVVIEW(xi); + switch(method) { + case 0 : {T = gsl_multiroot_fdfsolver_hybridsj;; break; } + case 1 : {T = gsl_multiroot_fdfsolver_hybridj; break; } + case 2 : {T = gsl_multiroot_fdfsolver_newton; break; } + case 3 : {T = gsl_multiroot_fdfsolver_gnewton; break; } + default: ERROR(BAD_CODE); + } + s = gsl_multiroot_fdfsolver_alloc (T, my_func.n); + + gsl_multiroot_fdfsolver_set (s, &my_func, V(xi)); + + do { + status = gsl_multiroot_fdfsolver_iterate (s); + + solp[iter*solc+0] = iter+1; + + int k; + for(k=0;kx,k); + } + for(k=xin;k<2*xin;k++) { + solp[iter*solc+k+1] = gsl_vector_get(s->f,k-xin); + } + + iter++; + if (status) /* check if solver is stuck */ + break; + + status = + gsl_multiroot_test_residual (s->f, epsabs); + } + while (status == GSL_CONTINUE && iter < maxit); + + int i,j; + for (i=iter; if); + + int k; + for(k=0;kx,k); + } + + iter++; + if (status) /* check if solver is stuck */ + break; + + status = gsl_multifit_test_delta (s->dx, s->x, epsabs, epsrel); + } + while (status == GSL_CONTINUE && iter < maxit); + + int i,j; + for (i=iter; iJ, 0.0, M(cov)); + + gsl_multifit_fdfsolver_free (s); + OK +} + + +////////////////////////////////////////////////////// + + +#define RAN(C,F) case C: { for(k=0;k