| author | Francesco Abbate <francesco.bbt@gmail.com> | 2009年09月06日 20:38:01 +0000 |
|---|---|---|
| committer | Francesco Abbate <francesco.bbt@gmail.com> | 2009年09月06日 20:38:01 +0000 |
| commit | 2909c35d22f3b04cf6fd32fe79d1897b495f0d55 (patch) | |
| tree | a3a7a6d6fff4b5b9f6a447e1b7c8e45ea91f0441 | |
| parent | 8b5bf5e26a4b6f0a7460cb48d87403a032da0dd4 (diff) | |
| download | gsl-shell-2909c35d22f3b04cf6fd32fe79d1897b495f0d55.tar.gz | |
| -rw-r--r-- | Makefile | 4 | ||||
| -rw-r--r-- | cmatrix.c | 18 | ||||
| -rw-r--r-- | cmatrix.h | 21 | ||||
| -rw-r--r-- | cnlinfit.c | 2 | ||||
| -rw-r--r-- | errors.c | 73 | ||||
| -rw-r--r-- | igsl.lua | 14 | ||||
| -rw-r--r-- | integ.c | 1 | ||||
| -rw-r--r-- | lua-gsl.c | 7 | ||||
| -rw-r--r-- | lua-utils.c | 66 | ||||
| -rw-r--r-- | lua-utils.h | 16 | ||||
| -rw-r--r-- | matrix.c | 19 | ||||
| -rw-r--r-- | matrix.h | 20 | ||||
| -rw-r--r-- | matrix_headers_source.h | 31 | ||||
| -rw-r--r-- | matrix_helper_source.c | 68 | ||||
| -rw-r--r-- | nlinfit.c | 2 | ||||
| -rw-r--r-- | nlinfit_helper.c (renamed from solver-impl.c) | 2 | ||||
| -rw-r--r-- | nlinfit_helper.h (renamed from solver-impl.h) | 0 | ||||
| -rw-r--r-- | nlinfit_source.c | 75 | ||||
| -rw-r--r-- | ode.c | 37 | ||||
| -rw-r--r-- | ode.h (renamed from errors.h) | 58 | ||||
| -rw-r--r-- | ode_decls_source.c | 56 | ||||
| -rw-r--r-- | ode_solver.c | 98 | ||||
| -rw-r--r-- | ode_solver.h | 50 | ||||
| -rw-r--r-- | ode_source.c | 292 | ||||
| -rw-r--r-- | template_matrix_off.h | 7 | ||||
| -rw-r--r-- | template_matrix_on.h | 22 | ||||
| -rw-r--r-- | tests/ode-test.lua | 55 |
@@ -52,9 +52,9 @@ endif COMPILE = $(CC) $(CFLAGS) $(DEFS) $(INCLUDES)
-LUAGSL_SRC_FILES = math-types.c errors.c matrix.c cmatrix.c solver-impl.c \
+LUAGSL_SRC_FILES = math-types.c matrix.c cmatrix.c nlinfit_helper.c \
fdfsolver.c nlinfit.c cnlinfit.c lua-utils.c linalg.c \
- integ.c lua-gsl.c
+ integ.c ode_solver.c ode.c lua-gsl.c
SUBDIRS =
@@ -21,6 +21,8 @@ #include <lua.h> #include <lauxlib.h> #include <assert.h> +#include <string.h> +#include <gsl/gsl_complex_math.h> #include <gsl/gsl_matrix.h> #include <gsl/gsl_blas.h> #include <gsl/gsl_permutation.h> @@ -29,11 +31,26 @@ #include "lua-utils.h" #include "cmatrix.h" +static lua_Complex +value_retrieve_complex (gsl_complex v) +{ + return GSL_REAL(v) + GSL_IMAG(v) * I; +} + +static gsl_complex +value_assign_complex (lua_Complex v) +{ + gsl_complex z; + GSL_SET_COMPLEX(&z, creal(v), cimag(v)); + return z; +} + #define BASE_GSL_COMPLEX #include "template_matrix_on.h" #include "matrix_decls_source.c" #include "matrix_source.c" +#include "matrix_helper_source.c" #define OPER_ADD #include "template_matrix_oper_on.h" @@ -59,4 +76,5 @@ #include "template_matrix_oper_off.h" #undef OPER_DIV +#include "template_matrix_off.h" #undef BASE_GSL_COMPLEX @@ -22,20 +22,15 @@ #define CMATRIX_COMPLEX_H #include <lua.h> -#include <lauxlib.h> -#include "defs.h" - -extern void matrix_complex_register (lua_State *L); -extern gsl_matrix_complex * matrix_complex_push (lua_State *L, - int n1, int n2); -extern gsl_matrix_complex * matrix_complex_push_raw (lua_State *L, - int n1, int n2); -extern void matrix_complex_push_view (lua_State *L, - gsl_matrix_complex *m); +#include <gsl/gsl_matrix.h> -extern gsl_matrix_complex * matrix_complex_check (lua_State *L, int index); +#include "defs.h" +#include "math-types.h" -extern gsl_matrix_complex_view * \ - matrix_complex_check_view(lua_State *L, int idx); +#define BASE_GSL_COMPLEX +#include "template_matrix_on.h" +#include "matrix_headers_source.h" +#include "template_matrix_off.h" +#undef BASE_GSL_COMPLEX #endif diff --git a/cnlinfit.c b/cnlinfit.c index 1b01cdb2..2bea2f21 100644 --- a/cnlinfit.c +++ b/cnlinfit.c @@ -28,7 +28,7 @@ #include "matrix.h"
#include "cmatrix.h"
#include "fdfsolver.h"
-#include "solver-impl.h"
+#include "nlinfit_helper.h"
#include "lua-utils.h"
#include "math-types.h"
diff --git a/errors.c b/errors.c deleted file mode 100644 index d49dd82f..00000000 --- a/errors.c +++ /dev/null @@ -1,73 +0,0 @@ -
-/* errors.c
- *
- * Copyright (C) 2009 Francesco Abbate
- *
- * 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 3 of the License, or (at
- * your option) any later version.
- *
- * This program is distributed in the hope that it will be useful, but
- * WITHOUT ANY WARRANTY; without even the implied warranty of
- * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
- * General Public License for more details.
- *
- * You should have received a copy of the GNU General Public License
- * along with this program; if not, write to the Free Software
- * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
- */
-
-#include <lua.h>
-#include <lauxlib.h>
-#include <gsl/gsl_errno.h>
-
-#include "errors.h"
-
-static lua_State *sys_wide_lua_state = NULL;
-
-/*
-struct error_reg {
- int error_code;
- const char *id;
- const char *long_msg;
-};
-
-static const struct error_reg error_codes[] = {
- {GSL_EMAXITER, "maxiter"},
- {GSL_EROUND, "failed"},
- {GSL_ESING, "singular"},
- {GSL_EDIVERGE, "diverge"},
- {0, NULL}
-};
-
-static const char *
-integ_get_error_msg (int code)
-{
- const char *default_msg = "unknown";
- const struct error_reg *p;
- for (p = error_codes; p->msg; p++)
- {
- if (p->error_code == code)
- return p->msg;
- }
- return default_msg;
-}
-*/
-
-void
-my_default_handler (const char * reason,
- const char * file,
- int line,
- int gsl_errno)
-{
- luaL_error (sys_wide_lua_state, "GSL error: %s", reason);
-}
-
-void
-set_gsl_error_handler (lua_State *L)
-{
- gsl_error_handler_t *old_handler;
- sys_wide_lua_state = L;
- old_handler = gsl_set_error_handler (& my_default_handler);
-}
@@ -146,6 +146,20 @@ function matrix_columns (m, istart, iend) return mr end +function set(d, s) + local r, c = s:dims() + for i=0,r-1 do + for j=0,c-1 do + d:set(i,j,s:get(i,j)) + end + end +end + +function null(m) + local r, c = m:dims() + for i=0,r-1 do for j=0,c-1 do m:set(i,j,0) end end +end + function chop(m, eps) eps = eps and eps or 1e-4 m = m:copy() @@ -20,7 +20,6 @@ #include <lua.h> #include <lauxlib.h> -#include <gsl/gsl_errno.h> #include <gsl/gsl_integration.h> #include "integ.h" @@ -21,20 +21,21 @@ #include <lua.h> #include <lauxlib.h> -#include "errors.h" #include "nlinfit.h" #include "cnlinfit.h" #include "matrix.h" #include "cmatrix.h" #include "linalg.h" #include "integ.h" +#include "ode.h" +#include "ode_solver.h" static const struct luaL_Reg gsl_methods_dummy[] = {{NULL, NULL}}; int luaopen_gsl (lua_State *L) { - set_gsl_error_handler (L); + gsl_set_error_handler_off (); luaL_register (L, "gsl", gsl_methods_dummy); @@ -44,6 +45,8 @@ luaopen_gsl (lua_State *L) matrix_complex_register (L); linalg_register (L); integ_register (L); + ode_solver_register (L); + ode_register (L); return 1; } diff --git a/lua-utils.c b/lua-utils.c index 6ab90ac9..dc06e7d4 100644 --- a/lua-utils.c +++ b/lua-utils.c @@ -80,3 +80,69 @@ mlua_null_cache (lua_State *L, int index) lua_pushnil (L);
lua_setfield (L, index, CACHE_FIELD_NAME);
}
+
+int
+mlua_index_with_properties (lua_State *L,
+ const struct luaL_Reg *properties,
+ const struct luaL_Reg *methods)
+{
+ char const * key;
+ const struct luaL_Reg *reg;
+
+ key = lua_tostring (L, 2);
+ if (key == NULL)
+ return 0;
+
+ reg = mlua_find_method (properties, key);
+ if (reg)
+ {
+ return mlua_get_property (L, reg, true);
+ }
+
+ reg = mlua_find_method (methods, key);
+ if (reg)
+ {
+ lua_pushcfunction (L, reg->func);
+ return 1;
+ }
+
+ return 0;
+}
+
+void
+mlua_check_field_type (lua_State *L, int index, const char *key, int type,
+ const char *error_msg)
+{
+ lua_getfield (L, index, key);
+ if (lua_type (L, -1) != type)
+ {
+ if (error_msg)
+ luaL_error (L, "field \"%s\", ", key, error_msg);
+ else
+ luaL_error (L, "field \"%s\" should be an %s", key,
+ lua_typename (L, type));
+ }
+ lua_pop (L, 1);
+}
+
+lua_Number
+mlua_named_optnumber (lua_State *L, int index, const char *key,
+ lua_Number default_value)
+{
+ lua_Number r;
+ lua_getfield (L, index, key);
+ r = luaL_optnumber (L, -1, default_value);
+ lua_pop (L, 1);
+ return r;
+}
+
+const char *
+mlua_named_optstring (lua_State *L, int index, const char *key,
+ const char * default_value)
+{
+ const char * r;
+ lua_getfield (L, index, key);
+ r = luaL_optstring (L, -1, default_value);
+ lua_pop (L, 1);
+ return r;
+}
diff --git a/lua-utils.h b/lua-utils.h index 70bb7571..2c775c3d 100644 --- a/lua-utils.h +++ b/lua-utils.h @@ -33,4 +33,20 @@ mlua_find_method (const struct luaL_Reg *p, const char *key); extern void
mlua_null_cache (lua_State *L, int index);
+extern void
+mlua_check_field_type (lua_State *L, int index, const char *key, int type,
+ const char *error_msg);
+
+extern int
+mlua_index_with_properties (lua_State *L,
+ const struct luaL_Reg *properties,
+ const struct luaL_Reg *methods);
+extern const char *
+mlua_named_optstring (lua_State *L, int index, const char *key,
+ const char * default_value);
+
+extern lua_Number
+mlua_named_optnumber (lua_State *L, int index, const char *key,
+ lua_Number default_value);
+
#endif
@@ -21,6 +21,7 @@ #include <lua.h> #include <lauxlib.h> #include <assert.h> +#include <string.h> #include <gsl/gsl_matrix.h> #include <gsl/gsl_blas.h> #include <gsl/gsl_permutation.h> @@ -29,11 +30,28 @@ #include "matrix.h" #include "lua-utils.h" +/* danger: three mathematicians in the name of this function :-) */ +void +matrix_jacob_copy_cauchy_riemann (gsl_matrix *jreal, gsl_matrix_complex *jcmpl, + size_t n) +{ + size_t k; + for (k = 0; k < n; k++) + { + gsl_complex z = gsl_matrix_complex_get (jcmpl, k, k); + gsl_matrix_set (jreal, 2*k, 2*k, GSL_REAL(z)); + gsl_matrix_set (jreal, 2*k, 2*k+1, - GSL_IMAG(z)); + gsl_matrix_set (jreal, 2*k+1, 2*k, GSL_IMAG(z)); + gsl_matrix_set (jreal, 2*k+1, 2*k+1, GSL_REAL(z)); + } +} + #define BASE_DOUBLE #include "template_matrix_on.h" #include "matrix_decls_source.c" #include "matrix_source.c" +#include "matrix_helper_source.c" #define OPER_ADD #include "template_matrix_oper_on.h" @@ -59,4 +77,5 @@ #include "template_matrix_oper_off.h" #undef OPER_DIV +#include "template_matrix_off.h" #undef BASE_DOUBLE @@ -25,20 +25,16 @@ #include <gsl/gsl_matrix.h> #include "defs.h" - #include "math-types.h" -extern void matrix_register (lua_State *L); - -extern gsl_matrix * matrix_push (lua_State *L, int n1, int n2); -extern gsl_matrix * matrix_push_raw (lua_State *L, int n1, int n2); - -extern void matrix_push_view (lua_State *L, - gsl_matrix *m); - -extern gsl_matrix * matrix_check (lua_State *L, int index); - -extern gsl_matrix_view * matrix_check_view (lua_State *L, int idx); +extern void +matrix_jacob_copy_cauchy_riemann (gsl_matrix *jreal, gsl_matrix_complex *jcmpl, + size_t n); +#define BASE_DOUBLE +#include "template_matrix_on.h" +#include "matrix_headers_source.h" +#include "template_matrix_off.h" +#undef BASE_DOUBLE #endif diff --git a/matrix_headers_source.h b/matrix_headers_source.h new file mode 100644 index 00000000..65272519 --- /dev/null +++ b/matrix_headers_source.h @@ -0,0 +1,31 @@ +
+extern void FUNCTION (matrix, register) (lua_State *L);
+
+extern TYPE (gsl_matrix) * FUNCTION (matrix, push) (lua_State *L,
+ int n1, int n2);
+
+extern TYPE (gsl_matrix) * FUNCTION (matrix, push_raw) (lua_State *L,
+ int n1, int n2);
+
+extern void FUNCTION (matrix, push_view) (lua_State *L,
+ TYPE (gsl_matrix) *m);
+
+extern TYPE (gsl_matrix) * FUNCTION (matrix, check) (lua_State *L,
+ int index);
+
+extern VIEW (gsl_matrix) * FUNCTION (matrix, check_view) (lua_State *L,
+ int idx);
+
+
+/* matrix helper functions */
+extern void
+FUNCTION (matrix, set_view_and_push) (lua_State *L, int index, double *data,
+ size_t n1, size_t n2, const double *src);
+
+extern void
+FUNCTION (matrix, jacob_copy_real_to_cmpl) (double *dest_cmpl, double *src_real,
+ size_t n, size_t p, size_t mult);
+
+extern void
+FUNCTION (matrix, jacob_copy_cmpl_to_real) (double *dest_real, double *src_cmpl,
+ size_t n, size_t p, size_t mult);
diff --git a/matrix_helper_source.c b/matrix_helper_source.c new file mode 100644 index 00000000..5feb61fb --- /dev/null +++ b/matrix_helper_source.c @@ -0,0 +1,68 @@ +
+void
+FUNCTION (matrix, set_view_and_push) (lua_State *L, int index, double *data,
+ size_t n1, size_t n2, const double *src)
+{
+ VIEW (gsl_matrix) *view = FUNCTION (matrix, check_view) (L, index);
+ *view = FUNCTION (gsl_matrix, view_array) (data, n1, n2);
+ if (src)
+ memcpy (data, src, n1 * n2 * sizeof(BASE));
+ lua_pushvalue (L, index);
+}
+
+static void
+TYPE (copy_jacobian_raw) (double *cmpl, double *real, size_t n, size_t p,
+ size_t mult, bool inverse)
+{
+ gsl_matrix_view rview, cview;
+ gsl_vector_view vview;
+ size_t nu;
+
+ for (nu = 0; nu < mult; nu++)
+ {
+ rview = gsl_matrix_view_array_with_tda (real + nu*p, n, p, p * mult);
+ vview = gsl_vector_view_array_with_stride (cmpl + nu, mult, n * p);
+ cview = gsl_matrix_view_vector (& vview.vector, n, p);
+
+ if (inverse)
+ gsl_matrix_memcpy (& rview.matrix, & cview.matrix);
+ else
+ gsl_matrix_memcpy (& cview.matrix, & rview.matrix);
+ }
+
+ /*
+ gsl_vector_view dview, sview;
+ double *cp, *rp;
+ size_t k, nu;
+
+ for (nu = 0; nu < multiplicity; nu++)
+ {
+ cp = cmpl + nu;
+ rp = real + p*nu;
+ for (k = 0; k < p; k++, rp += 1, cp += multiplicity)
+ {
+ dview = gsl_vector_view_array_with_stride (cp, multiplicity * p, n);
+ sview = gsl_vector_view_array_with_stride (rp, multiplicity * p, n);
+ if (inverse)
+ gsl_vector_memcpy (& sview.vector, & dview.vector);
+ else
+ gsl_vector_memcpy (& dview.vector, & sview.vector);
+ }
+ }
+ */
+}
+
+void
+FUNCTION (matrix, jacob_copy_real_to_cmpl) (double *dest_cmpl, double *src_real,
+ size_t n, size_t p, size_t mult)
+{
+ TYPE (copy_jacobian_raw) (dest_cmpl, src_real, n, p, mult, false);
+}
+
+void
+FUNCTION (matrix, jacob_copy_cmpl_to_real) (double *dest_real, double *src_cmpl,
+ size_t n, size_t p, size_t mult)
+{
+ TYPE (copy_jacobian_raw) (src_cmpl, dest_real, n, p, mult, true);
+}
+
@@ -27,7 +27,7 @@ #include "matrix.h"
#include "fdfsolver.h"
-#include "solver-impl.h"
+#include "nlinfit_helper.h"
#include "lua-utils.h"
#include "math-types.h"
diff --git a/solver-impl.c b/nlinfit_helper.c index b1522d90..a39def54 100644 --- a/solver-impl.c +++ b/nlinfit_helper.c @@ -23,7 +23,7 @@ #include <gsl/gsl_vector.h>
#include <gsl/gsl_matrix.h>
-#include "solver-impl.h"
+#include "nlinfit_helper.h"
#include "matrix.h"
size_t
diff --git a/solver-impl.h b/nlinfit_helper.h index 8ee7af78..8ee7af78 100644 --- a/solver-impl.h +++ b/nlinfit_helper.h diff --git a/nlinfit_source.c b/nlinfit_source.c index 1e72b736..78a24ed6 100644 --- a/nlinfit_source.c +++ b/nlinfit_source.c @@ -139,41 +139,6 @@ FUNCTION (solver, check) (lua_State *L, int index) return sext;
}
-static void
-FUNCTION (set_matrix, view_and_push) (lua_State *L, int index,
- double *data, size_t n1, size_t n2)
-{
- VIEW (gsl_matrix) *view = FUNCTION (matrix, check_view) (L, index);
- *view = FUNCTION (gsl_matrix, view_array) (data, n1, n2);
- lua_pushvalue (L, index);
-}
-
-#if MULTIPLICITY >= 2
-static void
-copy_jacobian (double *cmpl, double *real, size_t n, size_t p,
- size_t multiplicity, bool inverse)
-{
- gsl_vector_view dview, sview;
- double *cp, *rp;
- size_t k, nu;
-
- for (nu = 0; nu < multiplicity; nu++)
- {
- cp = cmpl + nu;
- rp = real + p*nu;
- for (k = 0; k < p; k++, rp += 1, cp += multiplicity)
- {
- dview = gsl_vector_view_array_with_stride (cp, multiplicity * p, n);
- sview = gsl_vector_view_array_with_stride (rp, multiplicity * p, n);
- if (inverse)
- gsl_vector_memcpy (& sview.vector, & dview.vector);
- else
- gsl_vector_memcpy (& dview.vector, & sview.vector);
- }
- }
-}
-#endif
-
int
FUNCTION (solver, fdf_hook) (const gsl_vector * x, void * _params,
gsl_vector * f, gsl_matrix * J)
@@ -193,23 +158,25 @@ FUNCTION (solver, fdf_hook) (const gsl_vector * x, void * _params, lua_pushvalue (L, 3);
if (f)
- FUNCTION (set_matrix, view_and_push) (L, 4, f->data, n, 1);
+ FUNCTION (matrix, set_view_and_push) (L, 4, f->data, n, 1, NULL);
else
lua_pushnil (L);
if (J)
{
double *jptr = (MULTIPLICITY >= 2 ? params->j_raw->data : J->data);
- FUNCTION (set_matrix, view_and_push) (L, 5, jptr, n, p);
+ FUNCTION (matrix, set_view_and_push) (L, 5, jptr, n, p, NULL);
}
lua_call (L, nargs, 0);
#if MULTIPLICITY >= 2
if (J)
- copy_jacobian (params->j_raw->data, J->data, n, p, MULTIPLICITY, true);
+ {
+ double *dst = J->data, *src = params->j_raw->data;
+ FUNCTION (matrix, jacob_copy_cmpl_to_real) (dst, src, n, p, MULTIPLICITY);
+ }
#endif
-
return GSL_SUCCESS;
}
@@ -364,11 +331,8 @@ FUNCTION (solver, get_jacob) (lua_State *L) m = FUNCTION (matrix, push) (L, n, p);
-#if MULTIPLICITY >= 2
- copy_jacobian (m->data, src->data, n, p, MULTIPLICITY, false);
-#else
- gsl_matrix_memcpy (m, src);
-#endif
+ FUNCTION (matrix, jacob_copy_real_to_cmpl) (m->data, src->data, n, p,
+ MULTIPLICITY);
return 1;
}
@@ -376,25 +340,6 @@ FUNCTION (solver, get_jacob) (lua_State *L) int
FUNCTION (solver, index) (lua_State *L)
{
- char const * key;
- const struct luaL_Reg *reg;
-
- key = lua_tostring (L, 2);
- if (key == NULL)
- return 0;
-
- reg = mlua_find_method (FUNCTION (solver, properties), key);
- if (reg)
- {
- return mlua_get_property (L, reg, true);
- }
-
- reg = mlua_find_method (FUNCTION (solver, methods), key);
- if (reg)
- {
- lua_pushcfunction (L, reg->func);
- return 1;
- }
-
- return 0;
+ return mlua_index_with_properties (L, FUNCTION (solver, properties),
+ FUNCTION (solver, methods));
}
@@ -0,0 +1,37 @@ +
+/* ode.c
+ *
+ * Copyright (C) 2009 Francesco Abbate
+ *
+ * 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 3 of the License, or (at
+ * your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful, but
+ * WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+ * General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program; if not, write to the Free Software
+ * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
+ */
+
+#include <lua.h>
+#include <lauxlib.h>
+#include <assert.h>
+
+#include "ode.h"
+#include "ode_solver.h"
+#include "matrix.h"
+#include "lua-utils.h"
+
+#define BASE_DOUBLE
+#include "template_matrix_on.h"
+
+#include "ode_decls_source.c"
+#include "ode_source.c"
+
+#include "template_matrix_off.h"
+#undef BASE_DOUBLE
@@ -1,29 +1,29 @@ -
-/* errors.h
- *
- * Copyright (C) 2009 Francesco Abbate
- *
- * 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 3 of the License, or (at
- * your option) any later version.
- *
- * This program is distributed in the hope that it will be useful, but
- * WITHOUT ANY WARRANTY; without even the implied warranty of
- * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
- * General Public License for more details.
- *
- * You should have received a copy of the GNU General Public License
- * along with this program; if not, write to the Free Software
- * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
- */
-
-#ifndef ERRORS_H
-#define ERRORS_H
-
-#include <lua.h>
-#include "defs.h"
-
-extern void set_gsl_error_handler (lua_State *L);
-
-#endif
+ +/* ode.h + * + * Copyright (C) 2009 Francesco Abbate + * + * 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 3 of the License, or (at + * your option) any later version. + * + * This program is distributed in the hope that it will be useful, but + * WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + * General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program; if not, write to the Free Software + * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA. + */ + +#ifndef ODE_H +#define ODE_H + +#include <lua.h> +#include "defs.h" + +extern void ode_register (lua_State *L); + +#endif diff --git a/ode_decls_source.c b/ode_decls_source.c new file mode 100644 index 00000000..aa44d330 --- /dev/null +++ b/ode_decls_source.c @@ -0,0 +1,56 @@ +
+typedef struct {
+ lua_State *L;
+ TYPE (gsl_matrix) *y;
+ TYPE (gsl_matrix) *ybuff;
+#if MULTIPLICITY >= 2
+ TYPE (gsl_matrix) *J;
+#endif
+ double t, h;
+ size_t n; /* ODE system dimension */
+} TYPE (ode_params);
+
+#define ODE_DEFAULT_EPS_ABS 1e-4
+#define ODE_DEFAULT_EPS_REL 0.0
+#define ODE_DEFAULT_METHOD "rk8pd"
+
+static int FUNCTION (ode, evolve) (lua_State *L);
+static int FUNCTION (ode, set) (lua_State *L);
+static int FUNCTION (ode, new) (lua_State *L);
+static int FUNCTION (ode, index) (lua_State *L);
+
+static int FUNCTION (ode, get_t) (lua_State *L);
+static int FUNCTION (ode, get_y) (lua_State *L);
+
+static struct ode_solver * FUNCTION (ode, check) (lua_State *L, int index);
+
+
+static TYPE (ode_params) * FUNCTION (ode_params, push) (lua_State *L,
+ size_t dim, double h);
+static int FUNCTION (ode_params, free) (lua_State *L);
+static TYPE (ode_params) * FUNCTION (ode_params, check) (lua_State *L, int idx);
+
+static const struct luaL_Reg FUNCTION (ode, methods)[] = {
+ {"set", FUNCTION (ode, set)},
+ {"evolve", FUNCTION (ode, evolve)},
+ {NULL, NULL}
+};
+
+static const struct luaL_Reg FUNCTION (ode, properties)[] = {
+ {"t", FUNCTION (ode, get_t)},
+ {"y", FUNCTION (ode, get_y)},
+ {NULL, NULL}
+};
+
+static const struct luaL_Reg FUNCTION (ode, functions)[] = {
+ {PREFIX "ode", FUNCTION (ode, new)},
+ {NULL, NULL}
+};
+
+static const struct luaL_Reg FUNCTION (ode_params, methods)[] = {
+ {"__gc", FUNCTION (ode_params, free)},
+ {NULL, NULL}
+};
+
+static char const * const TYPE (name_ode) = "GSL." PREFIX "ode";
+static char const * const TYPE (name_ode_params) = "GSL." PREFIX "odep";
diff --git a/ode_solver.c b/ode_solver.c new file mode 100644 index 00000000..d81975c4 --- /dev/null +++ b/ode_solver.c @@ -0,0 +1,98 @@ +
+#include <lua.h>
+#include <lauxlib.h>
+#include <string.h>
+
+#include "ode_solver.h"
+
+struct method_entry {
+ const char *name;
+ const gsl_odeiv_step_type **method_type;
+};
+
+static char const * const ode_solver_type_name = "GSL.ode_solver";
+
+#define ODEIV_METHOD(t) #t, & gsl_odeiv_step_ ## t
+
+static struct method_entry methods_table[] = {
+ {ODEIV_METHOD (rk2)},
+ {ODEIV_METHOD (rk4)},
+ {ODEIV_METHOD (rkf45)},
+ {ODEIV_METHOD (rkck)},
+ {ODEIV_METHOD (rk8pd)},
+ {ODEIV_METHOD (rk2imp)},
+ {ODEIV_METHOD (rk4imp)},
+ {ODEIV_METHOD (bsimp)},
+ {ODEIV_METHOD (gear1)},
+ {ODEIV_METHOD (gear2)},
+ {NULL, NULL}
+};
+#undef ODEIV_METHOD
+
+static int
+ode_solver_dealloc (lua_State *L);
+
+static const struct luaL_Reg ode_solver_methods[] = {
+ {"__gc", ode_solver_dealloc},
+ {NULL, NULL}
+};
+
+struct ode_solver *
+ode_solver_push_new (lua_State *L, const gsl_odeiv_step_type *type,
+ size_t dim, double eps_abs, double eps_rel)
+{
+ struct ode_solver *s;
+
+ s = lua_newuserdata (L, sizeof (struct ode_solver));
+
+ s->step = gsl_odeiv_step_alloc (type, dim);
+ s->ctrl = gsl_odeiv_control_y_new (eps_abs, eps_rel);
+ s->evol = gsl_odeiv_evolve_alloc (dim);
+
+ s->dimension = dim;
+
+ luaL_getmetatable (L, ode_solver_type_name);
+ lua_setmetatable (L, -2);
+
+ return s;
+}
+
+
+struct ode_solver *
+check_ode_solver (lua_State *L, int index)
+{
+ return luaL_checkudata (L, index, ode_solver_type_name);
+}
+
+int
+ode_solver_dealloc (lua_State *L)
+{
+ struct ode_solver *s = check_ode_solver (L, 1);
+
+ gsl_odeiv_evolve_free (s->evol);
+ gsl_odeiv_control_free (s->ctrl);
+ gsl_odeiv_step_free (s->step);
+
+ return 0;
+}
+
+const gsl_odeiv_step_type *
+method_lookup (const char *method, const gsl_odeiv_step_type *default_type)
+{
+ const struct method_entry *p;
+ for (p = methods_table; p->name; p++)
+ {
+ if (strcmp (p->name, method) == 0)
+ return *(p->method_type);
+ }
+ return default_type;
+}
+
+void
+ode_solver_register (lua_State *L)
+{
+ /* ode solver declaration */
+ luaL_newmetatable (L, ode_solver_type_name);
+ luaL_register (L, NULL, ode_solver_methods);
+ lua_pop (L, 1);
+}
diff --git a/ode_solver.h b/ode_solver.h new file mode 100644 index 00000000..dcd9cd50 --- /dev/null +++ b/ode_solver.h @@ -0,0 +1,50 @@ +
+/* ode-solver.h
+ *
+ * Copyright (C) 2009 Francesco Abbate
+ *
+ * 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 3 of the License, or (at
+ * your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful, but
+ * WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+ * General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program; if not, write to the Free Software
+ * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
+ */
+
+#ifndef ODE_SOLVER_H
+#define ODE_SOLVER_H
+
+#include <lua.h>
+#include <gsl/gsl_vector.h>
+#include <gsl/gsl_odeiv.h>
+#include "defs.h"
+
+struct ode_solver {
+ gsl_odeiv_step * step;
+ gsl_odeiv_control * ctrl;
+ gsl_odeiv_evolve * evol;
+
+ size_t dimension;
+};
+
+extern struct ode_solver *
+ode_solver_push_new (lua_State *L, const gsl_odeiv_step_type *type,
+ size_t dim, double eps_abs, double eps_rel);
+
+extern struct ode_solver *
+check_ode_solver (lua_State *L, int index);
+
+extern void
+ode_solver_register (lua_State *L);
+
+extern const gsl_odeiv_step_type *
+method_lookup (const char *method, const gsl_odeiv_step_type *default_type);
+
+#endif
diff --git a/ode_source.c b/ode_source.c new file mode 100644 index 00000000..477c340c --- /dev/null +++ b/ode_source.c @@ -0,0 +1,292 @@ +
+/* ode_source.c
+ *
+ * Copyright (C) 2009 Francesco Abbate
+ *
+ * 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 3 of the License, or (at
+ * your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful, but
+ * WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+ * General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program; if not, write to the Free Software
+ * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
+ */
+
+TYPE (ode_params) *
+FUNCTION (ode_params, push) (lua_State *L, size_t dim, double h)
+{
+ TYPE (ode_params) *p = lua_newuserdata (L, sizeof(TYPE (ode_params)));
+
+ p->L = L;
+ p->n = dim;
+ p->h = h;
+ p->y = FUNCTION (gsl_matrix, alloc) (dim, 1);
+ p->ybuff = FUNCTION (gsl_matrix, alloc) (dim, 1);
+#if MULTIPLICITY >= 2
+ p->J = FUNCTION (gsl_matrix, alloc) (dim, dim);
+#endif
+
+ luaL_getmetatable (L, TYPE (name_ode_params));
+ lua_setmetatable (L, -2);
+
+ return p;
+}
+
+int
+FUNCTION (ode_params, free) (lua_State *L)
+{
+ TYPE (ode_params) *p = FUNCTION (ode_params, check) (L, 1);
+ FUNCTION (gsl_matrix, free) (p->y);
+ FUNCTION (gsl_matrix, free) (p->ybuff);
+#if MULTIPLICITY >= 2
+ FUNCTION (gsl_matrix, free) (p->J);
+#endif
+ return 0;
+}
+
+TYPE (ode_params) *
+FUNCTION (ode_params, check) (lua_State *L, int index)
+{
+ return luaL_checkudata (L, index, TYPE (name_ode_params));
+}
+
+static TYPE (ode_params) *
+FUNCTION (ode_params, get) (lua_State *L, int index)
+{
+ TYPE (ode_params) *p;
+ lua_getfield (L, index, "params");
+ p = luaL_checkudata (L, -1, TYPE (name_ode_params));
+ lua_pop (L, 1);
+ return p;
+}
+
+static int
+FUNCTION (ode, hook_f) (double t, const double y[], double f[],
+ void *params)
+{
+ TYPE (ode_params) *p = params;
+ const int args_index = 2;
+ lua_State *L = p->L;
+
+ /* push the ode system function */
+ lua_getfield (L, 1, "f");
+
+ lua_pushnumber (L, t);
+
+ FUNCTION (matrix, set_view_and_push) (L, args_index, p->ybuff->data,
+ p->n, 1, y);
+ FUNCTION (matrix, set_view_and_push) (L, args_index+1, f, p->n, 1, NULL);
+
+ lua_call (L, 3, 0);
+
+ return GSL_SUCCESS;
+}
+
+static int
+FUNCTION (ode, hook_jacob) (double t, const double y[], double *dfdy,
+ double dfdt[], void *params)
+{
+ TYPE (ode_params) *p = params;
+ const int args_index = 2;
+ lua_State *L = p->L;
+ size_t n = p->n;
+#if MULTIPLICITY >= 2
+ double *jacob = p->J->data;
+#else
+ double *jacob = dfdy;
+#endif
+
+ /* push the ode system function */
+ lua_getfield (L, 1, "df");
+
+ lua_pushnumber (L, t);
+
+ FUNCTION (matrix, set_view_and_push) (L, args_index, p->ybuff->data,
+ n, 1, y);
+ FUNCTION (matrix, set_view_and_push) (L, args_index+1, jacob, n, n, NULL);
+ FUNCTION (matrix, set_view_and_push) (L, args_index+2, dfdt, n, 1, NULL);
+
+ lua_call (L, 4, 0);
+
+#if MULTIPLICITY == 2
+ {
+ gsl_matrix_view dfdy_view = gsl_matrix_view_array (dfdy, 2*n, 2*n);
+ matrix_jacob_copy_cauchy_riemann (& dfdy_view.matrix, jacob, n);
+ }
+#elif MULTIPLICITY > 2
+#error MULTIPLICITY > 2 not supported
+#endif
+
+ return GSL_SUCCESS;
+}
+
+int
+FUNCTION (ode, new) (lua_State *L)
+{
+ const gsl_odeiv_step_type * T;
+ struct ode_solver *s;
+ double eps_abs, eps_rel;
+ const char *method;
+ int n;
+
+ luaL_checktype (L, 1, LUA_TTABLE);
+
+ lua_getfield (L, 1, "n");
+ if (! lua_isinteger (L, -1))
+ luaL_error (L, "missing dimension 'n' of the ODE system");
+ n = lua_tointeger (L, -1);
+ lua_pop (L, 1);
+
+ mlua_check_field_type (L, 1, "f", LUA_TFUNCTION, NULL);
+ mlua_check_field_type (L, 1, "df", LUA_TFUNCTION,
+ "Jacobian function expected");
+
+ eps_abs = mlua_named_optnumber (L, 1, "eps_abs", ODE_DEFAULT_EPS_ABS);
+ eps_rel = mlua_named_optnumber (L, 1, "eps_rel", ODE_DEFAULT_EPS_REL);
+ method = mlua_named_optstring (L, 1, "method", ODE_DEFAULT_METHOD);
+
+ T = method_lookup (method, gsl_odeiv_step_rk8pd);
+
+ s = ode_solver_push_new (L, T, n, eps_abs, eps_rel);
+
+ lua_setfield (L, 1, "solver");
+
+ luaL_getmetatable (L, TYPE (name_ode));
+ lua_setmetatable (L, -2);
+
+ return 1;
+}
+
+struct ode_solver *
+FUNCTION (ode, check) (lua_State *L, int index)
+{
+ if (lua_getmetatable(L, index))
+ {
+ struct ode_solver *s;
+
+ lua_getfield(L, LUA_REGISTRYINDEX, TYPE (name_ode));
+ if (! lua_rawequal(L, -1, -2))
+ luaL_typerror (L, index, "ODE solver");
+ lua_pop (L, 2);
+
+ lua_getfield (L, 1, "solver");
+ s = check_ode_solver (L, -1);
+ lua_pop (L, 1);
+
+ assert (s != NULL);
+
+ return s;
+ }
+
+ lua_pop (L, 1);
+ return NULL;
+}
+
+int
+FUNCTION (ode, set) (lua_State *L)
+{
+ struct ode_solver *s = FUNCTION (ode, check) (L, 1);
+ TYPE (ode_params) *p;
+ TYPE (gsl_matrix) *y;
+ double t;
+
+ t = luaL_checknumber (L, 2);
+ y = FUNCTION (matrix, check) (L, 3);
+
+ p = FUNCTION (ode_params, push) (L, s->dimension, 1e-6);
+
+ p->t = t;
+ FUNCTION (gsl_matrix, memcpy) (p->y, y);
+
+ lua_setfield (L, 1, "params");
+
+ mlua_null_cache (L, 1);
+
+ return 0;
+}
+
+int
+FUNCTION (ode, evolve) (lua_State *L)
+{
+ struct ode_solver *s = FUNCTION (ode, check) (L, 1);
+ TYPE (ode_params) *p = FUNCTION (ode_params, get) (L, 1);
+ gsl_odeiv_system system[1];
+ double t1;
+ int status;
+
+ t1 = luaL_checknumber (L, 2);
+
+ p->h = luaL_optnumber (L, 3, p->h);
+
+ lua_settop (L, 1);
+
+ FUNCTION (matrix, push_view) (L, NULL);
+ FUNCTION (matrix, push_view) (L, NULL);
+ FUNCTION (matrix, push_view) (L, NULL);
+
+ system->function = & FUNCTION (ode, hook_f);
+ system->jacobian = & FUNCTION (ode, hook_jacob);
+ system->dimension = s->dimension;
+ system->params = p;
+
+ status = gsl_odeiv_evolve_apply (s->evol, s->ctrl, s->step,
+ system, & p->t, t1,
+ & p->h, p->y->data);
+
+ if (status != GSL_SUCCESS)
+ luaL_error (L, "error in ODE evolve: %s", gsl_strerror (status));
+
+ mlua_null_cache (L, 1);
+
+ return 0;
+}
+
+int
+FUNCTION (ode, get_t) (lua_State *L)
+{
+ TYPE (ode_params) *p = FUNCTION (ode_params, get) (L, 1);
+ lua_pushnumber (L, p->t);
+ return 1;
+}
+
+
+int
+FUNCTION (ode, get_y) (lua_State *L)
+{
+ TYPE (ode_params) *p = FUNCTION (ode_params, get) (L, 1);
+ TYPE (gsl_matrix) *y = FUNCTION (matrix, push) (L, p->n, 1);
+ FUNCTION (gsl_matrix, memcpy) (y, p->y);
+ return 1;
+}
+
+int
+FUNCTION (ode, index) (lua_State *L)
+{
+ return mlua_index_with_properties (L,
+ FUNCTION (ode, properties),
+ FUNCTION (ode, methods));
+}
+
+void
+FUNCTION (ode, register) (lua_State *L)
+{
+ luaL_newmetatable (L, TYPE (name_ode_params));
+ lua_pushvalue (L, -1);
+ lua_setfield (L, -2, "__index");
+ luaL_register (L, NULL, FUNCTION (ode_params, methods));
+ lua_pop (L, 1);
+
+ luaL_newmetatable (L, TYPE (name_ode));
+ lua_pushcfunction (L, FUNCTION (ode, index));
+ lua_setfield (L, -2, "__index");
+ luaL_register (L, NULL, FUNCTION (ode, methods));
+ lua_pop (L, 1);
+
+ luaL_register (L, NULL, FUNCTION (ode, functions));
+}
diff --git a/template_matrix_off.h b/template_matrix_off.h index b42c3d93..dd39ed46 100644 --- a/template_matrix_off.h +++ b/template_matrix_off.h @@ -17,3 +17,10 @@ #undef LUA_TYPE #undef LUA_FUNCTION +#undef LUAL_FUNCTION +#undef BLAS_FUNCTION + +#ifdef value_retrieve +#undef value_retrieve +#undef value_assign +#endif diff --git a/template_matrix_on.h b/template_matrix_on.h index 04168485..8560e590 100644 --- a/template_matrix_on.h +++ b/template_matrix_on.h @@ -32,20 +32,6 @@ #define PREFIX "c" #define BASE_TYPE TYPE_COMPLEX -static lua_Complex -value_retrieve_complex (gsl_complex v) -{ - return GSL_REAL(v) + GSL_IMAG(v) * I; -} - -static gsl_complex -value_assign_complex (lua_Complex v) -{ - gsl_complex z; - GSL_SET_COMPLEX(&z, creal(v), cimag(v)); - return z; -} - #elif defined(BASE_DOUBLE) #define BASE double #define SHORT @@ -57,8 +43,8 @@ value_assign_complex (lua_Complex v) #define PREFIX "" #define BASE_TYPE TYPE_REAL -#define value_retrieve(x) x -#define value_assign(x) x +#define value_retrieve(x) (x) +#define value_assign(x) (x) #else #error unknown BASE_ directive @@ -70,6 +56,8 @@ value_assign_complex (lua_Complex v) #define CONCAT3(a,b,c) CONCAT3x(a,b,c) #define CONCAT4x(a,b,c,d) a ## _ ## b ## _ ## c ## _ ## d #define CONCAT4(a,b,c,d) CONCAT4x(a,b,c,d) +#define MYCATx(a,b) a ## b +#define MYCAT(a,b) MYCATx(a,b) #if defined(BASE_DOUBLE) #define FUNCTION(dir,name) CONCAT2(dir,name) @@ -86,8 +74,6 @@ value_assign_complex (lua_Complex v) #endif #define LUA_TYPE CONCAT2(lua,LUA_SHORTM) -#define MYCATx(a,b) a ## b -#define MYCAT(a,b) MYCATx(a,b) #define LUA_FUNCTION(oper) CONCAT2(lua,MYCAT(oper,LUA_SHORT)) #define LUAL_FUNCTION(oper) CONCAT2(luaL,MYCAT(oper,LUA_SHORT)) #define BLAS_FUNCTION(name) CONCAT2(gsl_blas,MYCAT(BLAS_ID,name)) diff --git a/tests/ode-test.lua b/tests/ode-test.lua new file mode 100644 index 00000000..42c567b5 --- /dev/null +++ b/tests/ode-test.lua @@ -0,0 +1,55 @@ +
+require 'igsl'
+
+function test_ode1()
+ local k, w = -0.3, 4
+ local m = tmatrix {{k, -w}, {w, k}}
+
+ local function cexpf(t, y, f)
+ set(f, gsl.mul(m, y))
+ end
+
+ local function cexpdf(t, y, dfdy, dfdt)
+ set(dfdy, m)
+ null(dfdt)
+ end
+
+ return gsl.ode {f= cexpf, df= cexpdf, n= 2, method= 'bsimp', eps_abs= 1e-6}
+end
+
+function test_ode2()
+ local mu = 10.0
+
+ local function lorhf(t, y, f)
+ f:set(0,0, y[1])
+ f:set(1,0, -y[0] - mu*y[1]*(y[0]*y[0] - 1))
+ end
+
+ local function lorhdf(t, y, dfdy, dfdt)
+ dfdy:set(0,0, 0.0)
+ dfdy:set(0,1, 1.0)
+ dfdy:set(1,0, -2.0*mu*y[0]*y[1] - 1.0)
+ dfdy:set(1,1, -mu*(y[0]*y[0] - 1.0))
+
+ null(dfdt)
+ end
+
+ return gsl.ode {f= lorhf, df= lorhf, n= 2}
+end
+
+function vai(s, tf, step)
+ local r = {}
+ repeat
+ s:evolve(tf, step)
+ r[#r+1] = string.format('%g, %g, %g', s.t, s.y[0], s.y[1])
+ io.write("\r" .. tostring(s.t))
+ until s.t >= tf
+ return table.concat(r, '\n')
+end
+
+function ivai(s, tf)
+ repeat
+ s:evolve(tf)
+ print(string.format('%g, %g, %g', s.t, s.y[0], s.y[1]))
+ until s.t >= tf
+end
|