gsl-shell.git - gsl-shell

diff --git a/Makefile b/Makefile
index 6036ea6e..2eb3e33c 100644
--- a/Makefile
+++ b/Makefile

@@ -51,7 +51,8 @@ C_SRC_FILES = gs-types.c lcomplex.c matrix.c matrix_arith.c nlinfit_helper.c \

integ.c ode_solver.c ode.c random.c randist.c \

pdf.c cdf.c sf.c fmultimin.c gradcheck.c fdfmultimin.c \

multimin.c eigen-systems.c mlinear.c bspline.c interp.c \

- cmatrix.c cnlinfit.c code.c fft.c lua-graph.c lua-gsl.c

+ cmatrix.c cnlinfit.c code.c fft.c lua-graph.c lu_decomp_real.c \

+ lu_decomp_complex.c lu_decomp.c lua-gsl.c

LUA_BASE_DIRS =

LUA_BASE_FILES = igsl.lua base.lua integ.lua csv.lua

diff --git a/gs-types.c b/gs-types.c
index f34c92d0..3ea85c4c 100644
--- a/gs-types.c
+++ b/gs-types.c

@@ -24,6 +24,9 @@ static int gs_type_string (lua_State *L);

#define GS_FMULTIMIN_NAME_DEF "GSL.fmmin"

#define GS_BSPLINE_NAME_DEF "GSL.bspline"

#define GS_INTERP_NAME_DEF "GSL.interp"

+#define GS_LU_DECOMP_NAME_DEF "GSL.LUdec"

+#define GS_CLU_DECOMP_NAME_DEF "GSL.cLUdec"

+#define GS_QR_DECOMP_NAME_DEF "GSL.QRdec"

#ifdef AGG_PLOT_ENABLED

#define GS_DRAW_PLOT_NAME_DEF "GSL.plot"

#define GS_DRAW_SCALABLE_NAME_DEF NULL

@@ -59,6 +62,9 @@ const struct gs_type gs_type_table[] = {

MY_EXPAND(FMULTIMIN, "f multimin solver"),

MY_EXPAND(BSPLINE, "B-spline"),

MY_EXPAND(INTERP, "Interpolation object"),

+ MY_EXPAND(LU_DECOMP, "real matrix LU decomposition"),

+ MY_EXPAND(CLU_DECOMP, "complex matrix LU decomposition"),

+ MY_EXPAND(QR_DECOMP, "QR decomposition"),

#ifdef AGG_PLOT_ENABLED

MY_EXPAND(DRAW_PLOT, "plot"),

MY_EXPAND(DRAW_SCALABLE, "graphical object"),

diff --git a/gs-types.h b/gs-types.h
index 75220653..40e74071 100644
--- a/gs-types.h
+++ b/gs-types.h

@@ -25,6 +25,9 @@ enum gs_type_e {

GS_FMULTIMIN,

GS_BSPLINE,

GS_INTERP,

+ GS_LU_DECOMP,

+ GS_CLU_DECOMP,

+ GS_QR_DECOMP,

#ifdef AGG_PLOT_ENABLED

GS_DRAW_PLOT,

GS_DRAW_SCALABLE, /* derived types are declared only after their base class */

diff --git a/lu_decomp.c b/lu_decomp.c
new file mode 100644
index 00000000..7af9990d
--- /dev/null
+++ b/lu_decomp.c

@@ -0,0 +1,51 @@

+

+#include <lua.h>

+#include <lauxlib.h>

+

+#include "lu_decomp.h"

+#include "lu_decomp_imp.h"

+#include "gs-types.h"

+#include "matrix.h"

+#include "cmatrix.h"

+

+static int lu_decomp (lua_State *L);

+

+const struct luaL_Reg lu_decomp_functions[] = {

+ {"LU", lu_decomp},

+ {NULL, NULL}

+};

+

+int

+lu_decomp (lua_State *L)

+{

+ if (gs_is_userdata (L, 1, GS_MATRIX))

+ {

+ gsl_matrix *m = lua_touserdata (L, 1);

+ if (m->size1 != m->size2)

+ return gs_type_error (L, 1, "square matrix");

+ return lu_decomp_raw (L, m->size1, m);

+ }

+ else if (gs_is_userdata (L, 1, GS_CMATRIX))

+ {

+ gsl_matrix_complex *m = lua_touserdata (L, 1);

+ if (m->size1 != m->size2)

+ return gs_type_error (L, 1, "square matrix");

+ return lu_decomp_complex_raw (L, m->size1, m);

+ }

+

+ return gs_type_error (L, 1, "matrix");

+}

+

+void

+lu_decomp_register (lua_State *L)

+{

+ luaL_newmetatable (L, GS_METATABLE(GS_LU_DECOMP));

+ luaL_register (L, NULL, lu_decomp_metatable);

+ lua_pop (L, 1);

+

+ luaL_newmetatable (L, GS_METATABLE(GS_CLU_DECOMP));

+ luaL_register (L, NULL, lu_decomp_complex_metatable);

+ lua_pop (L, 1);

+

+ luaL_register (L, NULL, lu_decomp_functions);

+}

diff --git a/lu_decomp.h b/lu_decomp.h
new file mode 100644
index 00000000..96eec6a3
--- /dev/null
+++ b/lu_decomp.h

@@ -0,0 +1,9 @@

+

+#ifndef LU_DECOMP_H

+#define LU_DECOMP_H

+

+#include <lua.h>

+

+extern void lu_decomp_register (lua_State *L);

+

+#endif

diff --git a/lu_decomp_complex.c b/lu_decomp_complex.c
new file mode 100644
index 00000000..32de3058
--- /dev/null
+++ b/lu_decomp_complex.c

@@ -0,0 +1,17 @@

+

+#include <lua.h>

+#include <lauxlib.h>

+

+#include <gsl/gsl_linalg.h>

+

+#include "gs-types.h"

+#include "matrix.h"

+#include "cmatrix.h"

+#include "lcomplex.h"

+#include "lua-utils.h"

+#include "lu_decomp_imp.h"

+

+#define BASE_GSL_COMPLEX

+#include "template_matrix_on.h"

+#include "lu_decomp_source.c"

+#include "template_matrix_off.h"

diff --git a/lu_decomp_imp.h b/lu_decomp_imp.h
new file mode 100644
index 00000000..e93eb40c
--- /dev/null
+++ b/lu_decomp_imp.h

@@ -0,0 +1,15 @@

+

+#ifndef LU_DECOMP_IMP_H

+#define LU_DECOMP_IMP_H

+

+#include <gsl/gsl_matrix.h>

+

+#include <lua.h>

+

+extern const struct luaL_Reg lu_decomp_metatable[];

+extern const struct luaL_Reg lu_decomp_complex_metatable[];

+

+extern int lu_decomp_raw (lua_State *L, size_t n, gsl_matrix *m);

+extern int lu_decomp_complex_raw (lua_State *L, size_t n, gsl_matrix_complex *m);

+

+#endif

diff --git a/lu_decomp_real.c b/lu_decomp_real.c
new file mode 100644
index 00000000..0fd1b6f1
--- /dev/null
+++ b/lu_decomp_real.c

@@ -0,0 +1,15 @@

+

+#include <lua.h>

+#include <lauxlib.h>

+

+#include <gsl/gsl_linalg.h>

+

+#include "gs-types.h"

+#include "matrix.h"

+#include "lua-utils.h"

+#include "lu_decomp_imp.h"

+

+#define BASE_DOUBLE

+#include "template_matrix_on.h"

+#include "lu_decomp_source.c"

+#include "template_matrix_off.h"

diff --git a/lu_decomp_source.c b/lu_decomp_source.c
new file mode 100644
index 00000000..754a33f6
--- /dev/null
+++ b/lu_decomp_source.c

@@ -0,0 +1,189 @@

+

+typedef struct {

+ TYPE(gsl_matrix) *m;

+ gsl_permutation *p;

+ int signum;

+} TYPE(lu_decomp);

+

+static int FUNCTION(lu_decomp, free) (lua_State *L);

+static int FUNCTION(lu_decomp, index) (lua_State *L);

+static int FUNCTION(lu_decomp, solve) (lua_State *L);

+static int FUNCTION(lu_decomp, det) (lua_State *L);

+static int FUNCTION(lu_decomp, L) (lua_State *L);

+static int FUNCTION(lu_decomp, U) (lua_State *L);

+static int FUNCTION(lu_decomp, P) (lua_State *L);

+static int FUNCTION(lu_decomp, signum) (lua_State *L);

+

+const struct luaL_Reg FUNCTION (lu_decomp, metatable)[] = {

+ {"__gc", FUNCTION(lu_decomp, free)},

+ {"__index", FUNCTION(lu_decomp, index)},

+ {NULL, NULL}

+};

+

+static const struct luaL_Reg FUNCTION (lu_decomp, methods)[] = {

+ {"solve", FUNCTION(lu_decomp, solve)},

+ {"det", FUNCTION(lu_decomp, det)},

+ {NULL, NULL}

+};

+

+static const struct luaL_Reg FUNCTION(lu_decomp, properties)[] = {

+ {"L", FUNCTION(lu_decomp, L)},

+ {"U", FUNCTION(lu_decomp, U)},

+ {"P", FUNCTION(lu_decomp, P)},

+ {"signum", FUNCTION(lu_decomp, signum)},

+ {NULL, NULL}

+};

+

+int

+FUNCTION(lu_decomp, raw) (lua_State *L, size_t n, TYPE(gsl_matrix) *m)

+{

+ int status;

+ TYPE(lu_decomp) *lu = gs_new_object (sizeof(TYPE(lu_decomp)), L,

+ GS_TYPE(LU_DECOMP));

+

+ lu->m = FUNCTION(gsl_matrix, alloc) (n, n);

+ FUNCTION(gsl_matrix, memcpy) (lu->m, m);

+

+ lu->p = gsl_permutation_alloc (n);

+ if (lu->p == NULL)

+ return luaL_error (L, "out of memory");

+

+ status = FUNCTION(gsl_linalg, LU_decomp) (lu->m, lu->p, &lu->signum);

+

+ if (status != GSL_SUCCESS)

+ {

+ return luaL_error (L, "error during LU decomposition: %s",

+ gsl_strerror (status));

+ }

+

+ lua_newtable (L);

+ lua_setfenv (L, -2);

+

+ return 1;

+}

+

+int

+FUNCTION(lu_decomp, free) (lua_State *L)

+{

+ TYPE(lu_decomp) *lu = gs_check_userdata (L, 1, GS_TYPE(LU_DECOMP));

+ gsl_permutation_free (lu->p);

+ FUNCTION(gsl_matrix, free) (lu->m);

+ return 0;

+}

+

+int

+FUNCTION(lu_decomp, L) (lua_State *L)

+{

+ TYPE(lu_decomp) *lu = gs_check_userdata (L, 1, GS_TYPE(LU_DECOMP));

+ size_t n = lu->m->size1;

+ TYPE(gsl_matrix) *r = FUNCTION (matrix, push_raw) (L, n, n);

+ LUA_TYPE *rp = (LUA_TYPE *) r->data, *mp = (LUA_TYPE *) lu->m->data;

+ int i;

+

+ for (i = 0; i < n; i++)

+ {

+ LUA_TYPE *lm1 = mp + i, *lm2 = mp + n;

+ for (; mp < lm1; rp++, mp++)

+ *rp = *mp;

+

+ *(rp++) = (LUA_TYPE) 1;

+ mp++;

+

+ for (; mp < lm2; rp++, mp++)

+ *rp = (LUA_TYPE) 0;

+ }

+

+ return 1;

+}

+

+

+int

+FUNCTION(lu_decomp, U) (lua_State *L)

+{

+ TYPE(lu_decomp) *lu = gs_check_userdata (L, 1, GS_TYPE(LU_DECOMP));

+ size_t n = lu->m->size1;

+ TYPE(gsl_matrix) *r = FUNCTION (matrix, push_raw) (L, n, n);

+ LUA_TYPE *rp = (LUA_TYPE *) r->data, *mp = (LUA_TYPE *) lu->m->data;

+ int i;

+

+ for (i = 0; i < n; i++)

+ {

+ LUA_TYPE *lm1 = mp + i, *lm2 = mp + n;

+ for (; mp < lm1; rp++, mp++)

+ *rp = (LUA_TYPE) 0;

+

+ *(rp++) = *(mp++);

+

+ for (; mp < lm2; rp++, mp++)

+ *rp = *mp;

+ }

+

+ return 1;

+}

+

+

+int

+FUNCTION(lu_decomp, P) (lua_State *L)

+{

+ TYPE(lu_decomp) *lu = gs_check_userdata (L, 1, GS_TYPE(LU_DECOMP));

+ size_t n = lu->m->size1;

+ gsl_matrix *pm = matrix_push (L, n, n);

+ double *ptr = pm->data;

+ gsl_permutation *p = lu->p;

+ int j;

+

+ for (j = 0; j < n; j++)

+ {

+ int pj = p->data[j];

+ ptr[j * n + pj] = 1.0;

+ }

+

+ return 1;

+}

+

+int

+FUNCTION(lu_decomp, signum) (lua_State *L)

+{

+ TYPE(lu_decomp) *lu = gs_check_userdata (L, 1, GS_TYPE(LU_DECOMP));

+ lua_pushnumber (L, lu->signum);

+ return 1;

+}

+

+int

+FUNCTION(lu_decomp, solve) (lua_State *L)

+{

+ TYPE(lu_decomp) *lu = gs_check_userdata (L, 1, GS_TYPE(LU_DECOMP));

+ TYPE(gsl_matrix) *b = FUNCTION(matrix, check) (L, 2);

+ VIEW(gsl_vector) b1 = FUNCTION(gsl_matrix, column) (b, 0);

+ size_t n = lu->m->size1;

+

+ if (b->size1 != n || b->size2 != 1)

+ return gs_type_error (L, 2, "column matrix");

+

+ {

+ TYPE(gsl_matrix) *x = FUNCTION(matrix, push_raw) (L, n, 1);

+ VIEW(gsl_vector) x1 = FUNCTION(gsl_matrix, column) (x, 0);

+ FUNCTION (gsl_linalg, LU_solve) (lu->m, lu->p, &b1.vector, &x1.vector);

+ }

+

+ return 1;

+}

+

+int

+FUNCTION(lu_decomp, det) (lua_State *L)

+{

+ TYPE(lu_decomp) *lu = gs_check_userdata (L, 1, GS_TYPE(LU_DECOMP));

+ BASE d = FUNCTION (gsl_linalg, LU_det) (lu->m, lu->signum);

+ LUA_TYPE *zp = (LUA_TYPE *) &d;

+ LUA_FUNCTION(push) (L, *zp);

+ return 1;

+}

+

+int

+FUNCTION(lu_decomp, index) (lua_State *L)

+{

+ return mlua_index_with_properties (L,

+ FUNCTION(lu_decomp, properties),

+ FUNCTION(lu_decomp, methods),

+ true);

+}

diff --git a/lua-gsl.c b/lua-gsl.c
index 8d117a4b..236feef8 100644
--- a/lua-gsl.c
+++ b/lua-gsl.c

@@ -46,6 +46,7 @@

#include "mlinear.h"

#include "bspline.h"

#include "interp.h"

+#include "lu_decomp.h"

#ifdef AGG_PLOT_ENABLED

#include "lua-graph.h"

@@ -102,6 +103,7 @@ luaopen_gsl (lua_State *L)

mlinear_register (L);

bspline_register (L);

interp_register (L);

+ lu_decomp_register (L);

fft_register (L);

matrix_complex_register (L);