diff --git a/Config.hs b/Config.hs --- a/Config.hs +++ b/Config.hs @@ -79,7 +79,19 @@ ++ (join $ map ("-L"++) $ extraLibDirs buildInfo) ++ " -lgsl -lgslcblas" ++ " > /dev/null 2> /dev/null" +-- test for odeiv2 +gslodeiv2 bInfo buildInfo = + "echo \"#include \nint main(){return 0;}\"" + ++" > " ++ (buildDir bInfo) ++ "/dummy.c; gcc " + ++ (buildDir bInfo) ++ "/dummy.c " + ++ (join $ ccOptions buildInfo) ++ " " + ++ (join $ cppOptions buildInfo) ++ " " + ++ (join $ map ("-I"++) $ includeDirs buildInfo) + ++" -o " ++ (buildDir bInfo) ++ "/dummy " + ++ (join $ map ("-L"++) $ extraLibDirs buildInfo) ++ " -lgsl -lgslcblas" + ++ " > /dev/null 2> /dev/null" + checkCommand c = (ExitSuccess ==) `fmap` system c -- test different configurations until the first one works @@ -133,9 +145,11 @@ return (Just emptyBuildInfo { buildable = False }, []) Just ops -> do putStrLn $ " OK " ++ ops - g <- checkCommand $ gsl112 bInfo buildInfo - let hbi = if g - then emptyBuildInfo { extraLibs = words ops} - else emptyBuildInfo { extraLibs = words ops, ccOptions = ["-DGSL110"]} + g1 <- checkCommand $ gsl112 bInfo buildInfo + let op1 = if g1 then "" else "-DGSL110" + g2 <- checkCommand $ gslodeiv2 bInfo buildInfo + let op2 = if g2 then "" else "-DGSLODE1" + opts = filter (not.null) [op1,op2] + let hbi = emptyBuildInfo { extraLibs = words ops, ccOptions = opts } return (Just hbi, []) diff --git a/INSTALL.md b/INSTALL.md --- a/INSTALL.md +++ b/INSTALL.md @@ -33,6 +33,7 @@ We use this [GSL binary](http://www.miscdebris.net/blog/2009/04/20/mingw-345-binaries-of-gnu-scientific-library-112-for-use-with-mingw-and-visual-c/), and blas/lapack dlls built with g77 (contributed by Gilberto Camara). All required files are in [gsl-lapack-windows.zip][winpack]. +(Due to [issue 21](https://github.com/albertoruiz/hmatrix/issues/21) we need hmatrix-0.13.1.0.) 1) Install the Haskell Platform (tested on 2011.2.0.1) @@ -44,11 +45,11 @@ 3.a) In a msys shell the installation should be fully automatic: - $ cabal install hmatrix --extra-lib-dir=${GSL} --extra-include-dir=${GSL} + $ cabal install hmatrix-0.13.1.0 --extra-lib-dir=${GSL} --extra-include-dir=${GSL} 3.b) Alternatively, in a normal windows cmd: - > cabal unpack hmatrix + > cabal unpack hmatrix-0.13.1.0 Edit hmatrix.cabal, in line 28 change build-type to "Simple", and then diff --git a/THANKS.md b/THANKS.md --- a/THANKS.md +++ b/THANKS.md @@ -113,3 +113,5 @@ - Jeffrey Burdges reported a glpk link problem on OS/X +- Jian Zhang reported the Windows installation problem due to new ODE interface. + diff --git a/hmatrix.cabal b/hmatrix.cabal --- a/hmatrix.cabal +++ b/hmatrix.cabal @@ -1,5 +1,5 @@ Name: hmatrix -Version: 0.14.0.0 +Version: 0.14.0.1 License: GPL License-file: LICENSE Author: Alberto Ruiz @@ -56,6 +56,7 @@ extra-source-files: lib/Numeric/LinearAlgebra/LAPACK/lapack-aux.h, lib/Numeric/LinearAlgebra/LAPACK/clapack.h + lib/Numeric/GSL/gsl-ode.c flag dd description: svd = zgesdd diff --git a/lib/Numeric/GSL/gsl-aux.c b/lib/Numeric/GSL/gsl-aux.c --- a/lib/Numeric/GSL/gsl-aux.c +++ b/lib/Numeric/GSL/gsl-aux.c @@ -32,7 +32,6 @@ #include #include #include -#include #include #include #include @@ -1331,89 +1330,7 @@ #undef RAN ////////////////////////////////////////////////////// -// ODE -typedef struct {int n; int (*f)(double,int, const double*, int, double *); int (*j)(double,int, const double*, int, int, double*);} Tode; - -int odefunc (double t, const double y[], double f[], void *params) { - Tode * P = (Tode*) params; - (P->f)(t,P->n,y,P->n,f); - return GSL_SUCCESS; -} - -int odejac (double t, const double y[], double *dfdy, double dfdt[], void *params) { - Tode * P = ((Tode*) params); - (P->j)(t,P->n,y,P->n,P->n,dfdy); - int j; - for (j=0; j< P->n; j++) - dfdt[j] = 0.0; - return GSL_SUCCESS; -} - - -int ode(int method, double h, double eps_abs, double eps_rel, - int f(double, int, const double*, int, double*), - int jac(double, int, const double*, int, int, double*), - KRVEC(xi), KRVEC(ts), RMAT(sol)) { - - const gsl_odeiv2_step_type * T; - - switch(method) { - case 0 : {T = gsl_odeiv2_step_rk2; break; } - case 1 : {T = gsl_odeiv2_step_rk4; break; } - case 2 : {T = gsl_odeiv2_step_rkf45; break; } - case 3 : {T = gsl_odeiv2_step_rkck; break; } - case 4 : {T = gsl_odeiv2_step_rk8pd; break; } - case 5 : {T = gsl_odeiv2_step_rk2imp; break; } - case 6 : {T = gsl_odeiv2_step_rk4imp; break; } - case 7 : {T = gsl_odeiv2_step_bsimp; break; } - case 8 : {T = gsl_odeiv2_step_rk1imp; break; } - case 9 : {T = gsl_odeiv2_step_msadams; break; } - case 10: {T = gsl_odeiv2_step_msbdf; break; } - default: ERROR(BAD_CODE); - } - - Tode P; - P.f = f; - P.j = jac; - P.n = xin; - - gsl_odeiv2_system sys = {odefunc, odejac, xin, &P}; - - gsl_odeiv2_driver * d = - gsl_odeiv2_driver_alloc_y_new (&sys, T, h, eps_abs, eps_rel); - - double t = tsp[0]; - - double* y = (double*)calloc(xin,sizeof(double)); - int i,j; - int status; - for(i=0; i< xin; i++) { - y[i] = xip[i]; - solp[i] = xip[i]; - } - - for (i = 1; i < tsn ; i++) - { - double ti = tsp[i]; - - status = gsl_odeiv2_driver_apply (d, &t, ti, y); - - if (status != GSL_SUCCESS) { - printf ("error in ode, return value=%d\n", status); - break; - } - -// printf ("%.5e %.5e %.5e\n", t, y[0], y[1]); - - for(j=0; j + +typedef struct {int n; int (*f)(double,int, const double*, int, double *); int (*j)(double,int, const double*, int, int, double*);} Tode; + +int odefunc (double t, const double y[], double f[], void *params) { + Tode * P = (Tode*) params; + (P->f)(t,P->n,y,P->n,f); + return GSL_SUCCESS; +} + +int odejac (double t, const double y[], double *dfdy, double dfdt[], void *params) { + Tode * P = ((Tode*) params); + (P->j)(t,P->n,y,P->n,P->n,dfdy); + int j; + for (j=0; j< P->n; j++) + dfdt[j] = 0.0; + return GSL_SUCCESS; +} + + +int ode(int method, double h, double eps_abs, double eps_rel, + int f(double, int, const double*, int, double*), + int jac(double, int, const double*, int, int, double*), + KRVEC(xi), KRVEC(ts), RMAT(sol)) { + + const gsl_odeiv_step_type * T; + + switch(method) { + case 0 : {T = gsl_odeiv_step_rk2; break; } + case 1 : {T = gsl_odeiv_step_rk4; break; } + case 2 : {T = gsl_odeiv_step_rkf45; break; } + case 3 : {T = gsl_odeiv_step_rkck; break; } + case 4 : {T = gsl_odeiv_step_rk8pd; break; } + case 5 : {T = gsl_odeiv_step_rk2imp; break; } + case 6 : {T = gsl_odeiv_step_rk4imp; break; } + case 7 : {T = gsl_odeiv_step_bsimp; break; } + case 8 : { printf("Sorry: ODE rk1imp not available in this GSL version\n"); exit(0); } + case 9 : { printf("Sorry: ODE msadams not available in this GSL version\n"); exit(0); } + case 10: { printf("Sorry: ODE msbdf not available in this GSL version\n"); exit(0); } + default: ERROR(BAD_CODE); + } + + gsl_odeiv_step * s = gsl_odeiv_step_alloc (T, xin); + gsl_odeiv_control * c = gsl_odeiv_control_y_new (eps_abs, eps_rel); + gsl_odeiv_evolve * e = gsl_odeiv_evolve_alloc (xin); + + Tode P; + P.f = f; + P.j = jac; + P.n = xin; + + gsl_odeiv_system sys = {odefunc, odejac, xin, &P}; + + double t = tsp[0]; + + double* y = (double*)calloc(xin,sizeof(double)); + int i,j; + for(i=0; i< xin; i++) { + y[i] = xip[i]; + solp[i] = xip[i]; + } + + for (i = 1; i < tsn ; i++) + { + double ti = tsp[i]; + while (t < ti) + { + gsl_odeiv_evolve_apply (e, c, s, + &sys, + &t, ti, &h, + y); + // if (h < hmin) h = hmin; + } + for(j=0; j + +typedef struct {int n; int (*f)(double,int, const double*, int, double *); int (*j)(double,int, const double*, int, int, double*);} Tode; + +int odefunc (double t, const double y[], double f[], void *params) { + Tode * P = (Tode*) params; + (P->f)(t,P->n,y,P->n,f); + return GSL_SUCCESS; +} + +int odejac (double t, const double y[], double *dfdy, double dfdt[], void *params) { + Tode * P = ((Tode*) params); + (P->j)(t,P->n,y,P->n,P->n,dfdy); + int j; + for (j=0; j< P->n; j++) + dfdt[j] = 0.0; + return GSL_SUCCESS; +} + + +int ode(int method, double h, double eps_abs, double eps_rel, + int f(double, int, const double*, int, double*), + int jac(double, int, const double*, int, int, double*), + KRVEC(xi), KRVEC(ts), RMAT(sol)) { + + const gsl_odeiv2_step_type * T; + + switch(method) { + case 0 : {T = gsl_odeiv2_step_rk2; break; } + case 1 : {T = gsl_odeiv2_step_rk4; break; } + case 2 : {T = gsl_odeiv2_step_rkf45; break; } + case 3 : {T = gsl_odeiv2_step_rkck; break; } + case 4 : {T = gsl_odeiv2_step_rk8pd; break; } + case 5 : {T = gsl_odeiv2_step_rk2imp; break; } + case 6 : {T = gsl_odeiv2_step_rk4imp; break; } + case 7 : {T = gsl_odeiv2_step_bsimp; break; } + case 8 : {T = gsl_odeiv2_step_rk1imp; break; } + case 9 : {T = gsl_odeiv2_step_msadams; break; } + case 10: {T = gsl_odeiv2_step_msbdf; break; } + default: ERROR(BAD_CODE); + } + + Tode P; + P.f = f; + P.j = jac; + P.n = xin; + + gsl_odeiv2_system sys = {odefunc, odejac, xin, &P}; + + gsl_odeiv2_driver * d = + gsl_odeiv2_driver_alloc_y_new (&sys, T, h, eps_abs, eps_rel); + + double t = tsp[0]; + + double* y = (double*)calloc(xin,sizeof(double)); + int i,j; + int status; + for(i=0; i< xin; i++) { + y[i] = xip[i]; + solp[i] = xip[i]; + } + + for (i = 1; i < tsn ; i++) + { + double ti = tsp[i]; + + status = gsl_odeiv2_driver_apply (d, &t, ti, y); + + if (status != GSL_SUCCESS) { + printf ("error in ode, return value=%d\n", status); + break; + } + +// printf ("%.5e %.5e %.5e\n", t, y[0], y[1]); + + for(j=0; j