Permalink
Browse files

fall back to old ode if v2 not available

  • Loading branch information...
1 parent 354d599 commit d109c791aef6dbb6a15fb65efe5c42121728aed0 @albertoruiz committed Mar 21, 2012
Showing with 204 additions and 90 deletions.
  1. +18 −4 Config.hs
  2. +2 −1 hmatrix.cabal
  3. +2 −85 lib/Numeric/GSL/gsl-aux.c
  4. +182 −0 lib/Numeric/GSL/gsl-ode.c
View
@@ -79,6 +79,18 @@ gsl112 bInfo buildInfo =
++ (join $ map ("-L"++) $ extraLibDirs buildInfo) ++ " -lgsl -lgslcblas"
++ " > /dev/null 2> /dev/null"
+-- test for odeiv2
+gslodeiv2 bInfo buildInfo =
+ "echo \"#include <gsl/gsl_odeiv2.h>\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
@@ -133,9 +145,11 @@ config bInfo = do
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, [])
View
@@ -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: examples/deriv.hs
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
View
@@ -32,7 +32,6 @@
#include <gsl/gsl_complex_math.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>
-#include <gsl/gsl_odeiv2.h>
#include <gsl/gsl_multifit_nlin.h>
#include <string.h>
#include <stdio.h>
@@ -1331,89 +1330,7 @@ int random_vector(int seed, int code, RVEC(r)) {
#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<xin; j++) {
- solp[i*xin + j] = y[j];
- }
- }
-
- free(y);
- gsl_odeiv2_driver_free (d);
-
- return status;
-}
+#include "gsl-ode.c"
+//////////////////////////////////////////////////////
View
@@ -0,0 +1,182 @@
+
+#ifdef GSLODE1
+
+////////////////////////////// ODE V1 //////////////////////////////////////////
+
+#include <gsl/gsl_odeiv.h>
+
+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<xin; j++) {
+ solp[i*xin + j] = y[j];
+ }
+ }
+
+ free(y);
+ gsl_odeiv_evolve_free (e);
+ gsl_odeiv_control_free (c);
+ gsl_odeiv_step_free (s);
+ return 0;
+}
+
+#else
+
+///////////////////// ODE V2 ///////////////////////////////////////////////////
+
+#include <gsl/gsl_odeiv2.h>
+
+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<xin; j++) {
+ solp[i*xin + j] = y[j];
+ }
+ }
+
+ free(y);
+ gsl_odeiv2_driver_free (d);
+
+ return status;
+}
+
+#endif
+

0 comments on commit d109c79

Please sign in to comment.