author | francesco <francesco.bbt@gmail.com> | 2011年04月24日 12:44:25 +0200 |
---|---|---|
committer | francesco <francesco.bbt@gmail.com> | 2011年04月24日 12:44:25 +0200 |
commit | 90a9dc178acca62749df1393eacffd264350cae8 (patch) | |
tree | ec2e1518ec60fbd14393a11b6a2d42404737664b | |
parent | 2b1733183834a93f4cfcff03de15b372bc69b4ac (diff) | |
download | gsl-shell-90a9dc178acca62749df1393eacffd264350cae8.tar.gz |
-rw-r--r-- | Makefile | 13 | ||||
-rw-r--r-- | bspline.c | 145 | ||||
-rw-r--r-- | bspline.h | 8 | ||||
-rw-r--r-- | c-lmtest.c | 185 | ||||
-rw-r--r-- | clftest.lua | 2 | ||||
-rw-r--r-- | clinfit.lua | 2 | ||||
-rw-r--r-- | cmatrix.c | 59 | ||||
-rw-r--r-- | cmatrix.h | 35 | ||||
-rw-r--r-- | cmatrix.lua | 547 | ||||
-rw-r--r-- | cnlinfit.c | 40 | ||||
-rw-r--r-- | cnlinfit.h | 26 | ||||
-rw-r--r-- | code.c | 39 | ||||
-rw-r--r-- | code.h | 29 | ||||
-rw-r--r-- | debug-support.c | 23 | ||||
-rw-r--r-- | debug-support.h | 12 | ||||
-rw-r--r-- | eigen-systems.c | 435 | ||||
-rw-r--r-- | eigen-systems.h | 6 | ||||
-rw-r--r-- | examples/nlinfit-jit.lua | 57 | ||||
-rw-r--r-- | examples/nlinfit.lua | 212 | ||||
-rw-r--r-- | examples/ode-jit.lua | 8 | ||||
-rw-r--r-- | fdfmultimin.c | 411 | ||||
-rw-r--r-- | fft.c | 511 | ||||
-rw-r--r-- | fft.h | 7 | ||||
-rw-r--r-- | fmultimin.c | 340 | ||||
-rw-r--r-- | gradcheck.c | 135 | ||||
-rw-r--r-- | gs-types.c | 32 | ||||
-rw-r--r-- | gs-types.h | 19 | ||||
-rw-r--r-- | gsl-shell.c | 552 | ||||
-rw-r--r-- | integ.c | 450 | ||||
-rw-r--r-- | integ.h | 9 | ||||
-rw-r--r-- | interp.c | 206 | ||||
-rw-r--r-- | interp.h | 8 | ||||
-rw-r--r-- | lcomplex.c | 289 | ||||
-rw-r--r-- | lcomplex.h | 26 | ||||
-rw-r--r-- | linalg.c | 69 | ||||
-rw-r--r-- | linalg.h | 26 | ||||
-rw-r--r-- | lu_decomp.c | 51 | ||||
-rw-r--r-- | lu_decomp.h | 9 | ||||
-rw-r--r-- | lu_decomp_complex.c | 17 | ||||
-rw-r--r-- | lu_decomp_imp.h | 15 | ||||
-rw-r--r-- | lu_decomp_real.c | 15 | ||||
-rw-r--r-- | lu_decomp_source.c | 189 | ||||
-rw-r--r-- | lua-gsl.c | 43 | ||||
-rw-r--r-- | lua-utils.c | 62 | ||||
-rw-r--r-- | lua-utils.h | 15 | ||||
-rw-r--r-- | matrix-init.lua | 629 | ||||
-rw-r--r-- | matrix.c | 62 | ||||
-rw-r--r-- | matrix.h | 39 | ||||
-rw-r--r-- | matrix_arith.c | 676 | ||||
-rw-r--r-- | matrix_arith.h | 18 | ||||
-rw-r--r-- | matrix_decls_source.c | 54 | ||||
-rw-r--r-- | matrix_headers_source.h | 46 | ||||
-rw-r--r-- | matrix_helper_source.c | 77 | ||||
-rw-r--r-- | matrix_op_source.c | 168 | ||||
-rw-r--r-- | matrix_source.c | 407 | ||||
-rw-r--r-- | misc.lua | 55 | ||||
-rw-r--r-- | mlinear.c | 70 | ||||
-rw-r--r-- | mlinear.h | 8 | ||||
-rw-r--r-- | multimin.c | 46 | ||||
-rw-r--r-- | multimin.h | 40 | ||||
-rw-r--r-- | nlinfit.c | 39 | ||||
-rw-r--r-- | nlinfit.h | 26 | ||||
-rw-r--r-- | nlinfit_decls_source.c | 55 | ||||
-rw-r--r-- | nlinfit_helper.c | 52 | ||||
-rw-r--r-- | nlinfit_helper.h | 35 | ||||
-rw-r--r-- | nlinfit_source.c | 363 | ||||
-rw-r--r-- | num/lmfit.lua.in | 2 | ||||
-rw-r--r-- | ode.c | 38 | ||||
-rw-r--r-- | ode.h | 29 | ||||
-rw-r--r-- | ode_decls_source.c | 52 | ||||
-rw-r--r-- | ode_solver.c | 71 | ||||
-rw-r--r-- | ode_solver.h | 59 | ||||
-rw-r--r-- | ode_source.c | 303 | ||||
-rw-r--r-- | qr_decomp.c | 171 | ||||
-rw-r--r-- | qr_decomp.h | 8 | ||||
-rw-r--r-- | sf.c | 1 | ||||
-rw-r--r-- | sf_gener.c | 2 | ||||
-rw-r--r-- | sf_implement.h | 33 | ||||
-rw-r--r-- | template_matrix_off.h | 30 | ||||
-rw-r--r-- | template_matrix_on.h | 85 | ||||
-rw-r--r-- | template_matrix_oper_off.h | 13 | ||||
-rw-r--r-- | template_matrix_oper_on.h | 61 |
@@ -46,21 +46,12 @@ endif SUBDIRS = $(LUADIR) -C_SRC_FILES = gs-types.c lcomplex.c matrix.c matrix_arith.c nlinfit_helper.c \ - nlinfit.c lua-utils.c linalg.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 lu_decomp_real.c \ - lu_decomp_complex.c lu_decomp.c qr_decomp.c lua-gsl.c +C_SRC_FILES = gs-types.c lua-utils.c random.c randist.c \ + pdf.c cdf.c sf.c lua-graph.c lua-gsl.c LUA_BASE_DIRS = LUA_BASE_FILES = igsl.lua base.lua integ.lua csv.lua -ifeq ($(strip $(DEBUG)), yes) - C_SRC_FILES += debug-support.c -endif - ifeq ($(LUADIR), luajit2) LUAGSL_LIBS = $(LUADIR)/src/libluajit.a C_SRC_FILES += gsl-shell-jit.c diff --git a/bspline.c b/bspline.c deleted file mode 100644 index 585d4075..00000000 --- a/bspline.c +++ /dev/null @@ -1,145 +0,0 @@ - -#include <lua.h> -#include <lauxlib.h> -#include <assert.h> -#include <gsl/gsl_vector.h> -#include <gsl/gsl_bspline.h> - -#include "bspline.h" -#include "gs-types.h" -#include "matrix.h" - -struct bspline { - gsl_bspline_workspace *ws; - size_t k; - size_t nbreak; -}; - -static int bspline_new (lua_State *L); -static int bspline_free (lua_State *L); -static int bspline_eval (lua_State *L); -static int bspline_model (lua_State *L); - -static const struct luaL_Reg bspline_methods[] = { - {"__gc", bspline_free}, - {"model", bspline_model}, - {"eval", bspline_eval}, - {NULL, NULL} -}; - -static const struct luaL_Reg bspline_functions[] = { - {"bspline", bspline_new}, - {NULL, NULL} -}; - -int -bspline_new (lua_State *L) -{ - gsl_matrix *breaks = NULL; - gsl_vector_view brk; - double a, b; - int nbreak; - struct bspline *bs; - int nargs = 3; - int k = 4; - - if (lua_isnumber (L, 1)) - { - a = lua_tonumber (L, 1); - b = luaL_checknumber (L, 2); - nbreak = luaL_checkinteger (L, 3); - } - else - { - breaks = matrix_check (L, 1); - brk = gsl_matrix_column (breaks, 0); - nbreak = brk.vector.size; - nargs = 1; - } - - if (nbreak <= 2) - return luaL_error (L, "number of knots should be >= 2"); - - if (lua_gettop (L) >= nargs+1 && !lua_isnil (L, nargs+1)) - { - k = luaL_checkinteger (L, nargs+1); - if (k <= 0) - return luaL_error (L, "bspline order should be > 0"); - } - - bs = lua_newuserdata (L, sizeof(struct bspline)); - bs->ws = gsl_bspline_alloc ((size_t) k, (size_t) nbreak); - bs->k = k; - bs->nbreak = (size_t) nbreak; - - if (breaks) - gsl_bspline_knots (&brk.vector, bs->ws); - else - gsl_bspline_knots_uniform (a, b, bs->ws); - - gs_set_metatable (L, GS_BSPLINE); - - return 1; -} - -static struct bspline * -bspline_check (lua_State *L, int index) -{ - return gs_check_userdata (L, index, GS_BSPLINE); -} - -int -bspline_free (lua_State *L) -{ - struct bspline *bs = bspline_check (L, 1); - gsl_bspline_free (bs->ws); - return 0; -} - -int -bspline_model (lua_State *L) -{ - struct bspline *bs = bspline_check (L, 1); - gsl_matrix *x = matrix_check (L, 2); - size_t j, n = x->size1; - gsl_matrix *M; - - M = matrix_push_raw (L, n, bs->nbreak + bs->k - 2); - - for (j = 0; j < n; j++) - { - double xj = gsl_matrix_get (x, j, 0); - gsl_vector_view B = gsl_matrix_row (M, j); - gsl_bspline_eval(xj, &B.vector, bs->ws); - } - - return 1; -} - -int -bspline_eval (lua_State *L) -{ - struct bspline *bs = bspline_check (L, 1); - double x = luaL_checknumber (L, 2); - gsl_vector_view Bv; - gsl_matrix *B; - - B = matrix_push_raw (L, bs->nbreak + bs->k - 2, 1); - - Bv = gsl_matrix_column (B, 0); - gsl_bspline_eval(x, &Bv.vector, bs->ws); - - return 1; -} - -void -bspline_register (lua_State *L) -{ - luaL_newmetatable (L, GS_METATABLE(GS_BSPLINE)); - lua_pushvalue (L, -1); - lua_setfield (L, -2, "__index"); - luaL_register (L, NULL, bspline_methods); - lua_pop (L, 1); - - luaL_register (L, NULL, bspline_functions); -} diff --git a/bspline.h b/bspline.h deleted file mode 100644 index 5f8be894..00000000 --- a/bspline.h +++ /dev/null @@ -1,8 +0,0 @@ -#ifndef BSPLINE_SYSTEM_H -#define BSPLINE_SYSTEM_H - -#include <lua.h> - -extern void bspline_register (lua_State *L); - -#endif diff --git a/c-lmtest.c b/c-lmtest.c deleted file mode 100644 index 488dc0a1..00000000 --- a/c-lmtest.c +++ /dev/null @@ -1,185 +0,0 @@ - -#include <stdlib.h> -#include <stdio.h> -#include <gsl/gsl_rng.h> -#include <gsl/gsl_randist.h> -#include <gsl/gsl_vector.h> -#include <gsl/gsl_blas.h> -#include <gsl/gsl_multifit_nlin.h> - -/* expfit.c -- model functions for exponential + background */ - -struct data { - size_t n; - double * y; - double * sigma; -}; - -int -expb_f (const gsl_vector * x, void *data, - gsl_vector * f) -{ - size_t n = ((struct data *)data)->n; - double *y = ((struct data *)data)->y; - double *sigma = ((struct data *) data)->sigma; - - double A = gsl_vector_get (x, 0); - double lambda = gsl_vector_get (x, 1); - double b = gsl_vector_get (x, 2); - - size_t i; - - for (i = 0; i < n; i++) - { - /* Model Yi = A * exp(-lambda * i) + b */ - double t = i; - double Yi = A * exp (-lambda * t) + b; - gsl_vector_set (f, i, (Yi - y[i])/sigma[i]); - } - - return GSL_SUCCESS; -} - -int -expb_df (const gsl_vector * x, void *data, - gsl_matrix * J) -{ - size_t n = ((struct data *)data)->n; - double *sigma = ((struct data *) data)->sigma; - - double A = gsl_vector_get (x, 0); - double lambda = gsl_vector_get (x, 1); - - size_t i; - - for (i = 0; i < n; i++) - { - /* Jacobian matrix J(i,j) = dfi / dxj, */ - /* where fi = (Yi - yi)/sigma[i], */ - /* Yi = A * exp(-lambda * i) + b */ - /* and the xj are the parameters (A,lambda,b) */ - double t = i; - double s = sigma[i]; - double e = exp(-lambda * t); - gsl_matrix_set (J, i, 0, e/s); - gsl_matrix_set (J, i, 1, -t * A * e/s); - gsl_matrix_set (J, i, 2, 1/s); - } - return GSL_SUCCESS; -} - -int -expb_fdf (const gsl_vector * x, void *data, - gsl_vector * f, gsl_matrix * J) -{ - expb_f (x, data, f); - expb_df (x, data, J); - - return GSL_SUCCESS; -} - -#define N 40 - -void print_state (size_t iter, gsl_multifit_fdfsolver * s); - -int -main (void) -{ - const gsl_multifit_fdfsolver_type *T; - gsl_multifit_fdfsolver *s; - int status; - unsigned int i, iter = 0; - const size_t n = N; - const size_t p = 3; - - gsl_matrix *covar = gsl_matrix_alloc (p, p); - double y[N], sigma[N]; - struct data d = { n, y, sigma}; - gsl_multifit_function_fdf f; - double x_init[3] = { 1.0, 0.0, 0.0 }; - gsl_vector_view x = gsl_vector_view_array (x_init, p); - const gsl_rng_type * type; - gsl_rng * r; - - gsl_rng_env_setup(); - - type = gsl_rng_default; - r = gsl_rng_alloc (type); - - f.f = &expb_f; - f.df = &expb_df; - f.fdf = &expb_fdf; - f.n = n; - f.p = p; - f.params = &d; - - /* This is the data to be fitted */ - - for (i = 0; i < n; i++) - { - double t = i; - y[i] = 1.0 + 5 * exp (-0.1 * t) - + gsl_ran_gaussian (r, 0.1); - sigma[i] = 0.1; - printf ("data: %u %g %g\n", i, y[i], sigma[i]); - }; - - T = gsl_multifit_fdfsolver_lmsder; - s = gsl_multifit_fdfsolver_alloc (T, n, p); - gsl_multifit_fdfsolver_set (s, &f, &x.vector); - - print_state (iter, s); - - do - { - iter++; - status = gsl_multifit_fdfsolver_iterate (s); - - printf ("status = %s\n", gsl_strerror (status)); - - print_state (iter, s); - - if (status) - break; - - status = gsl_multifit_test_delta (s->dx, s->x, - 1e-4, 1e-4); - } - while (status == GSL_CONTINUE && iter < 500); - - gsl_multifit_covar (s->J, 0.0, covar); - -#define FIT(i) gsl_vector_get(s->x, i) -#define ERR(i) sqrt(gsl_matrix_get(covar,i,i)) - - { - double chi = gsl_blas_dnrm2(s->f); - double dof = n - p; - double c = GSL_MAX_DBL(1, chi / sqrt(dof)); - - printf("chisq/dof = %g\n", pow(chi, 2.0) / dof); - - printf ("A = %.5f +/- %.5f\n", FIT(0), c*ERR(0)); - printf ("lambda = %.5f +/- %.5f\n", FIT(1), c*ERR(1)); - printf ("b = %.5f +/- %.5f\n", FIT(2), c*ERR(2)); - } - - printf ("status = %s\n", gsl_strerror (status)); - - gsl_multifit_fdfsolver_free (s); - gsl_matrix_free (covar); - gsl_rng_free (r); - return 0; -} - -void -print_state (size_t iter, gsl_multifit_fdfsolver * s) -{ - printf ("iter: %3u x = % 15.8f % 15.8f % 15.8f " - "|f(x)| = %g\n", - iter, - gsl_vector_get (s->x, 0), - gsl_vector_get (s->x, 1), - gsl_vector_get (s->x, 2), - gsl_blas_dnrm2 (s->f)); -} diff --git a/clftest.lua b/clftest.lua index 51beefc8..7993891e 100644 --- a/clftest.lua +++ b/clftest.lua @@ -1,8 +1,6 @@ use 'gsl' -require 'cmatrix' - local linfit = require 'clinfit' local exp = math.exp diff --git a/clinfit.lua b/clinfit.lua index 32c31945..23d86360 100644 --- a/clinfit.lua +++ b/clinfit.lua @@ -1,6 +1,4 @@ -require 'cmatrix' - local ffi = require 'ffi' local cgsl = require 'cgsl' diff --git a/cmatrix.c b/cmatrix.c deleted file mode 100644 index a4cf42a5..00000000 --- a/cmatrix.c +++ /dev/null @@ -1,59 +0,0 @@ - -/* cmatrix.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 <string.h> -#include <gsl/gsl_complex_math.h> -#include <gsl/gsl_matrix.h> -#include <gsl/gsl_blas.h> -#include <gsl/gsl_permutation.h> -#include <gsl/gsl_linalg.h> - -#include "gs-types.h" -#include "lcomplex.h" -#include "lua-utils.h" -#include "cmatrix.h" -#include "matrix_arith.h" - -static Complex -value_retrieve_complex (gsl_complex v) -{ - return GSL_REAL(v) + GSL_IMAG(v) * I; -} - -static gsl_complex -value_assign_complex (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" - -#include "template_matrix_off.h" -#undef BASE_GSL_COMPLEX diff --git a/cmatrix.h b/cmatrix.h deleted file mode 100644 index 2da15e15..00000000 --- a/cmatrix.h +++ /dev/null @@ -1,35 +0,0 @@ - -/* cmatrix.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 CMATRIX_COMPLEX_H -#define CMATRIX_COMPLEX_H - -#include <lua.h> -#include <gsl/gsl_matrix.h> - -#include "defs.h" - -#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/cmatrix.lua b/cmatrix.lua deleted file mode 100644 index b3d6bf4e..00000000 --- a/cmatrix.lua +++ /dev/null @@ -1,547 +0,0 @@ - -local ffi = require 'ffi' -local cgsl = require 'cgsl' - -local sqrt, abs = math.sqrt, math.abs -local fmt = string.format - -local gsl_matrix = ffi.typeof('gsl_matrix') -local gsl_matrix_complex = ffi.typeof('gsl_matrix_complex') -local gsl_complex = ffi.typeof('complex') - -local gsl_check = require 'gsl-check' - -local function isreal(x) return type(x) == 'number' end - -local function cartesian(x) - if isreal(x) then - return x, 0 - else - return x[0], x[1] - end -end - -local function complex_conj(a) - local x, y = cartesian(z) - return gsl_complex(x, -y) -end - -local function complex_real(z) - local x = cartesian(z) - return x -end - -local function complex_imag(z) - local x, y = cartesian(z) - return y -end - -function matrix_alloc(n1, n2) - local block = cgsl.gsl_block_alloc(n1 * n2) - local m = gsl_matrix(n1, n2, n2, block.data, block, 1) - return m -end - -function matrix_calloc(n1, n2) - local block = cgsl.gsl_block_complex_alloc(n1 * n2) - local m = gsl_matrix_complex(n1, n2, n2, block.data, block, 1) - return m -end - -function matrix_new(n1, n2, f) - local m - if f then - m = matrix_alloc(n1, n2) - for i=0, n1-1 do - for j=0, n2-1 do - m.data[i*n2+j] = f(i, j) - end - end - else - local block = cgsl.gsl_block_calloc(n1 * n2) - m = gsl_matrix(n1, n2, n2, block.data, block, 1) - end - return m -end - -function matrix_cnew(n1, n2, f) - local m - if f then - m = matrix_calloc(n1, n2) - for i=0, n1-1 do - for j=0, n2-1 do - local z = f(i, j) - m.data[2*i*n2+2*j ] = z[0] - m.data[2*i*n2+2*j+1] = z[1] - end - end - else - local block = cgsl.gsl_block_complex_calloc(n1 * n2) - m = gsl_matrix_complex(n1, n2, n2, block.data, block, 1) - end - return m -end - -local function matrix_dim(m) - return m.size1, m.size2 -end - -local function itostr(im, signed) - local sign = im < 0 and '-' or (signed and '+' or '') - if im == 0 then return '' else - return sign .. (abs(im) == 1 and 'i' or fmt('%gi', abs(im))) - end -end - -local function recttostr(x, y, eps) - if abs(x) < eps then x = 0 end - if abs(y) < eps then y = 0 end - if x ~= 0 then - return (y == 0 and fmt('%g', x) or fmt('%g%s', x, itostr(y, true))) - else - return (y == 0 and '0' or itostr(y)) - end -end - -local function concat_pad(t, pad) - local sep = ' ' - local row - for i, s in ipairs(t) do - local x = string.rep(' ', pad - #s) .. s - row = row and (row .. sep .. x) or x - end - return row -end - -local function matrix_tostring_gen(sel) - return function(m) - local n1, n2 = m.size1, m.size2 - local sq = 0 - for i=0, n1-1 do - for j=0, n2-1 do - local x, y = sel(m, i, j) - sq = sq + x*x + y*y - end - end - local eps = sqrt(sq) * 1e-10 - - lsrow = {} - local lmax = 0 - for i=0, n1-1 do - local row = {} - for j=0, n2-1 do - local x, y = sel(m, i, j) - local s = recttostr(x, y, eps) - if #s > lmax then lmax = #s end - row[j+1] = s - end - lsrow[i+1] = row - end - - local ss = {} - for i=0, n1-1 do - ss[i+1] = '[ ' .. concat_pad(lsrow[i+1], lmax) .. ' ]' - end - - return table.concat(ss, '\n') - end -end - -local function matrix_copy(a) - local n1, n2 = a.size1, a.size2 - local b = matrix_alloc(n1, n2) - cgsl.gsl_matrix_memcpy(b, a) - return b -end - -local function matrix_complex_copy(a) - local n1, n2 = a.size1, a.size2 - local b = matrix_calloc(n1, n2) - cgsl.gsl_matrix_complex_memcpy(b, a) - return b -end - -local function matrix_col(m, j) - local r = matrix_alloc(m.size1, 1) - local tda = m.tda - for i = 0, m.size1 - 1 do - r.data[i] = m.data[i * tda + j] - end - return r -end - -local function matrix_row(m, i) - local r = matrix_alloc(1, m.size2) - local tda = m.tda - for j = 0, m.size2 - 1 do - r.data[j] = m.data[i * tda + j] - end - return r -end - -local function matrix_vect_def(t) - return matrix_new(#t, 1, |i| t[i+1]) -end - -local function get_typeid(a) - if isreal(a) then return true, true - elseif ffi.istype(gsl_complex, a) then return false, true - elseif ffi.istype(gsl_matrix, a) then return true, false - elseif ffi.istype(gsl_matrix_complex, a) then return false, false end -end - -local function mat_op_gen(n1, n2, opa, a, opb, b, oper) - local c = matrix_alloc(n1, n2) - for i = 0, n1-1 do - for j = 0, n2-1 do - local ar = opa(a,i,j) - local br = opb(b,i,j) - c.data[i*n2+j] = oper(ar, br) - end - end - return c -end - -local function mat_comp_op_gen(n1, n2, opa, a, opb, b, oper) - local c = matrix_calloc(n1, n2) - for i = 0, n1-1 do - for j = 0, n2-1 do - local ar, ai = opa(a,i,j) - local br, bi = opb(b,i,j) - local zr, zi = oper(ar, br, ai, bi) - c.data[2*i*n2+2*j ] = zr - c.data[2*i*n2+2*j+1] = zi - end - end - return c -end - -local function real_get(x) return x, 0 end -local function complex_get(z) return z[0], z[1] end -local function mat_real_get(m,i,j) return m.data[i*m.tda+j], 0 end - -local function mat_complex_get(m,i,j) - local idx = 2*i*m.tda+2*j - return m.data[idx], m.data[idx+1] -end - -local function selector(r, s) - if s then - return r and real_get or complex_get - else - return r and mat_real_get or mat_complex_get - end -end - -local function mat_complex_of_real(m) - local n1, n2 = m.size1, m.size2 - local mc = matrix_calloc(n1, n2) - for i=0, n1-1 do - for j=0, n2-1 do - mc.data[2*i*n2+2*j ] = m.data[i*n2+j] - mc.data[2*i*n2+2*j+1] = 0 - end - end - return mc -end - -local function opadd(ar, br, ai, bi) - if ai then return ar+br, ai+bi else return ar+br end -end - -local function opsub(ar, br, ai, bi) - if ai then return ar-br, ai-bi else return ar-br end -end - -local function opmul(ar, br, ai, bi) - if ai then return ar*br-ai*bi, ar*bi+ai*br else return ar*br end -end - -local function opdiv(ar, br, ai, bi) - if ai then - local d = br^2 + bi^2 - return (ar*br + ai*bi)/d, (-ar*bi + ai*br)/d - else - return ar/br - end -end - -local function vector_op(scalar_op, element_wise, no_inverse) - return function(a, b) - local ra, sa = get_typeid(a) - local rb, sb = get_typeid(b) - if not sb and no_inverse then - error 'invalid operation on matrix' - end - if sa and sb then - local ar, ai = cartesian(a) - local br, bi = cartesian(b) - local zr, zi = scalar_op(ar, br, ai, bi) - return gsl_complex(zr, zi) - elseif element_wise or sa or sb then - local sela, selb = selector(ra, sa), selector(rb, sb) - local n1 = (sa and b.size1 or a.size1) - local n2 = (sa and b.size2 or a.size2) - if ra and rb then - return mat_op_gen(n1, n2, sela, a, selb, b, scalar_op) - else - return mat_comp_op_gen(n1, n2, sela, a, selb, b, scalar_op) - end - else - if ra and rb then - local n1, n2 = a.size1, b.size2 - local c = matrix_new(n1, n2) - local NT = cgsl.CblasNoTrans - gsl_check(cgsl.gsl_blas_dgemm(NT, NT, 1, a, b, 1, c)) - return c - else - if ra then a = mat_complex_of_real(a) end - if rb then b = mat_complex_of_real(b) end - local n1, n2 = a.size1, b.size2 - local c = matrix_cnew(n1, n2) - local NT = cgsl.CblasNoTrans - gsl_check(cgsl.gsl_blas_zgemm(NT, NT, 1, a, b, 1, c)) - return c - end - end - end -end - -complex = { - new = gsl_complex, - conj = complex_conj, - real = complex_real, - imag = complex_imag, -} - -local generic_add = vector_op(opadd, true) -local generic_sub = vector_op(opsub, true) -local generic_mul = vector_op(opmul, false) -local generic_div = vector_op(opdiv, true, true) - -local complex_mt = { - - __add = generic_add, - __sub = generic_sub, - __mul = generic_mul, - __div = generic_div, - - __pow = function(z,n) - if isreal(n) then - return cgsl.gsl_complex_pow_real (z, n) - else - if isreal(z) then z = gsl_complex(z,0) end - return cgsl.gsl_complex_pow (z, n) - end - end, -} - -ffi.metatype(gsl_complex, complex_mt) - -matrix = { - new = matrix_new, - cnew = matrix_cnew, - alloc = matrix_alloc, - calloc = matrix_calloc, - copy = matrix_copy, - dim = matrix_dim, - vec = matrix_vect_def, - col = matrix_col, - row = matrix_row, - get = cgsl.gsl_matrix_get, - set = cgsl.gsl_matrix_set, -} - -local matrix_methods = { - col = matrix_col, - row = matrix_row, - get = cgsl.gsl_matrix_get, - set = cgsl.gsl_matrix_set, -} - -local matrix_mt = { - - __gc = function(m) if m.owner then cgsl.gsl_block_free(m.block) end end, - - __add = generic_add, - __sub = generic_sub, - __mul = generic_mul, - __div = generic_div, - - __index = function(m, k) - if type(k) == 'number' then - if m.size2 == 1 then - return cgsl.gsl_matrix_get(m, k, 0) - else - return matrix_row(m, k) - end - end - return matrix_methods[k] - end, - - - __newindex = function(m, k, v) - if type(k) == 'number' then - if m.size2 == 1 then - cgsl.gsl_matrix_set(m, k, 0, v) - else - local row = cgsl.gsl_matrix_submatrix(m, k, 0, 1, m.size2) - gsl_check(cgsl.gsl_matrix_memcpy(row, v)) - end - else - error 'cannot set a matrix field' - end - end, - - __tostring = matrix_tostring_gen(mat_real_get), -} - -ffi.metatype(gsl_matrix, matrix_mt) - -local matrix_complex_mt = { - - __gc = function(m) if m.owner then cgsl.gsl_block_free(m.block) end end, - - __add = generic_add, - __sub = generic_sub, - __mul = generic_mul, - __div = generic_div, - - __index = function(m, k) - if type(k) == 'number' then - if m.size2 == 1 then - return cgsl.gsl_matrix_complex_get(m, k, 0) - else - return matrix_row(m, k) - end - end - return matrix_complex_methods[k] - end, - - __tostring = matrix_tostring_gen(mat_complex_get), -} - -ffi.metatype(gsl_matrix_complex, matrix_complex_mt) - -local function cwrap(name) - local fc = cgsl['gsl_complex_' .. name] - local fr = math[name] - return function(x) - return isreal(x) and fr(x) or fc(x) - end -end - -gsl_function_list = { - 'sqrt', 'exp', 'log', 'log10', - 'sin', 'cos', 'sec', 'csc', 'tan', 'cot', - 'arcsin', 'arccos', 'arcsec', 'arccsc', 'arctan', 'arccot', - 'sinh', 'cosh', 'sech', 'csch', 'tanh', 'coth', - 'arcsinh', 'arccosh', 'arcsech', 'arccsch', 'arctanh', 'arccoth' -} - -for _, name in ipairs(gsl_function_list) do - complex[name] = cwrap(name) -end - -local function matrix_def(t) - local r, c = #t, #t[1] - local real = true - for i, row in ipairs(t) do - for j, x in ipairs(row) do - if not isreal(x) then - real = false - break - end - end - if not real then break end - end - if real then - local m = matrix_alloc(r, c) - for i= 0, r-1 do - local row = t[i+1] - for j = 0, c-1 do - m.data[i*c+j] = row[j+1] - end - end - return m - else - local m = matrix_calloc(r, c) - for i= 0, r-1 do - local row = t[i+1] - for j = 0, c-1 do - local x, y = cartesian(row[j+1]) - m.data[2*i*c+2*j ] = x - m.data[2*i*c+2*j+1] = y - end - end - return m - end -end - -local signum = ffi.new('int[1]') - -function matrix_inv(m) - local n = m.size1 - local lu = matrix_copy(m) - local p = ffi.gc(cgsl.gsl_permutation_alloc(n), cgsl.gsl_permutation_free) - gsl_check(cgsl.gsl_linalg_LU_decomp(lu, p, signum)) - local mi = matrix_alloc(n, n) - gsl_check(cgsl.gsl_linalg_LU_invert(lu, p, mi)) - return mi -end - -function matrix_solve(m, b) - local n = m.size1 - local lu = matrix_copy(m) - local p = ffi.gc(cgsl.gsl_permutation_alloc(n), cgsl.gsl_permutation_free) - gsl_check(cgsl.gsl_linalg_LU_decomp(lu, p, signum)) - local x = matrix_alloc(n, 1) - local xv = cgsl.gsl_matrix_column(x, 0) - local bv = cgsl.gsl_matrix_column(b, 0) - gsl_check(cgsl.gsl_linalg_LU_solve(lu, p, bv, xv)) - return x -end - -function matrix_complex_inv(m) - local n = m.size1 - local lu = matrix_complex_copy(m) - local p = ffi.gc(cgsl.gsl_permutation_alloc(n), cgsl.gsl_permutation_free) - gsl_check(cgsl.gsl_linalg_complex_LU_decomp(lu, p, signum)) - local mi = matrix_calloc(n, n) - gsl_check(cgsl.gsl_linalg_complex_LU_invert(lu, p, mi)) - return mi -end - -function matrix_complex_solve(m, b) - local n = m.size1 - local lu = matrix_complex_copy(m) - local p = ffi.gc(cgsl.gsl_permutation_alloc(n), cgsl.gsl_permutation_free) - gsl_check(cgsl.gsl_linalg_complex_LU_decomp(lu, p, signum)) - local x = matrix_calloc(n, 1) - local xv = cgsl.gsl_matrix_complex_column(x, 0) - local bv = cgsl.gsl_matrix_complex_column(b, 0) - gsl_check(cgsl.gsl_linalg_complex_LU_solve(lu, p, bv, xv)) - return x -end - -matrix.inv = function(m) - if ffi.istype(gsl_matrix, m) then - return matrix_inv(m) - else - return matrix_complex_inv(m) - end - end - -matrix.solve = function(m, b) - local mr = ffi.istype(gsl_matrix, m) - local br = ffi.istype(gsl_matrix, b) - if mr and br then - return matrix_solve(m, b) - else - if mr then m = mat_complex_of_real(m) end - if br then b = mat_complex_of_real(b) end - return matrix_complex_solve(m, b) - end - end - -matrix.def = matrix_def diff --git a/cnlinfit.c b/cnlinfit.c deleted file mode 100644 index ee76be3b..00000000 --- a/cnlinfit.c +++ /dev/null @@ -1,40 +0,0 @@ - -/* cnlinfit.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 <gsl/gsl_vector.h> -#include <gsl/gsl_matrix.h> -#include <gsl/gsl_multifit_nlin.h> - -#include "gs-types.h" -#include "matrix.h" -#include "cmatrix.h" -#include "nlinfit_helper.h" -#include "lua-utils.h" - -#define BASE_GSL_COMPLEX -#include "template_matrix_on.h" -#include "nlinfit_decls_source.c" -#include "nlinfit_source.c" -#include "template_matrix_off.h" - -#undef BASE_GSL_COMPLEX diff --git a/cnlinfit.h b/cnlinfit.h deleted file mode 100644 index 6643fea4..00000000 --- a/cnlinfit.h +++ /dev/null @@ -1,26 +0,0 @@ - -/* cnlinfit.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 CNLINFIT_H -#define CNLINFIT_H - -extern void solver_complex_register (lua_State *L); - -#endif diff --git a/code.c b/code.c deleted file mode 100644 index 471274d9..00000000 --- a/code.c +++ /dev/null @@ -1,39 +0,0 @@ - -/* code.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 <math.h> - -#include "ode.h" -#include "ode_solver.h" -#include "matrix.h" -#include "cmatrix.h" -#include "lua-utils.h" - -#define BASE_GSL_COMPLEX -#include "template_matrix_on.h" - -#include "ode_decls_source.c" -#include "ode_source.c" - -#include "template_matrix_off.h" -#undef BASE_DOUBLE diff --git a/code.h b/code.h deleted file mode 100644 index 25386091..00000000 --- a/code.h +++ /dev/null @@ -1,29 +0,0 @@ - -/* code.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 CODE_H -#define CODE_H - -#include <lua.h> -#include "defs.h" - -extern void ode_complex_register (lua_State *L); - -#endif diff --git a/debug-support.c b/debug-support.c deleted file mode 100644 index 78645c12..00000000 --- a/debug-support.c +++ /dev/null @@ -1,23 +0,0 @@ -#include "debug-support.h" - -#ifdef WIN32 - -void __attribute__((__stdcall__)) Sleep(long); - -void -msleep(int msec) -{ - Sleep (msec); -} -#else - -#include <unistd.h> - -void -msleep(int msec) -{ - unsigned long us = (msec >= 0 ? msec * 1000 : 0); - usleep (us); -} - -#endif diff --git a/debug-support.h b/debug-support.h deleted file mode 100644 index ce55e5b4..00000000 --- a/debug-support.h +++ /dev/null @@ -1,12 +0,0 @@ -#ifndef DEBUG_SUPPORT_H -#define DEBUG_SUPPORT_H - -#include "defs.h" - -__BEGIN_DECLS - -extern void msleep(int msec); - -__END_DECLS - -#endif diff --git a/eigen-systems.c b/eigen-systems.c deleted file mode 100644 index f7e798d5..00000000 --- a/eigen-systems.c +++ /dev/null @@ -1,435 +0,0 @@ - -#include <lua.h> -#include <lauxlib.h> -#include <assert.h> -#include <gsl/gsl_eigen.h> -#include <gsl/gsl_vector.h> -#include <gsl/gsl_matrix.h> -#include <gsl/gsl_complex_math.h> -#include <gsl/gsl_linalg.h> - -#include "matrix.h" -#include "cmatrix.h" -#include "eigen-systems.h" - -struct eigen_symm_cache { - size_t n; - gsl_eigen_symm_workspace *ws; - gsl_eigen_symmv_workspace *vws; - gsl_vector *diag; -}; - -struct eigen_herm_cache { - size_t n; - gsl_eigen_herm_workspace *ws; - gsl_eigen_hermv_workspace *vws; - gsl_vector_complex *diag; -}; - -struct eigen_nonsymm_cache { - size_t n; - gsl_eigen_nonsymm_workspace *ws; - gsl_eigen_nonsymmv_workspace *vws; - gsl_matrix *mcopy; -}; - -struct eigen_cache { - struct eigen_symm_cache symm[1]; - struct eigen_herm_cache herm[1]; - struct eigen_nonsymm_cache nonsymm[1]; -}; - -#define EIGEN_CACHE_MT_NAME "GSL.eigcache" - -static int eigen_symm (lua_State *L); -static int eigen_symmv (lua_State *L); -static int eigen_herm (lua_State *L); -static int eigen_hermv (lua_State *L); -static int eigen_nonsymm (lua_State *L); -static int eigen_nonsymmv (lua_State *L); -static int schur_decomp (lua_State *L); -static int eigen_cache_free (lua_State *L); - -static void eigen_push_cache (lua_State *L); - -static struct eigen_symm_cache * eigen_cache_symm_set (lua_State *L, size_t n); -static struct eigen_herm_cache * eigen_cache_herm_set (lua_State *L, size_t n); -static struct eigen_nonsymm_cache * eigen_cache_nonsymm_set (lua_State *L, size_t n); - -static const struct luaL_Reg eigen_cache_methods[] = { - {"__gc", eigen_cache_free}, - {NULL, NULL} -}; - -static const struct luaL_Reg eigen_functions[] = { - {"eigs", eigen_symm}, - {"eigsv", eigen_symmv}, - {"eigh", eigen_herm}, - {"eighv", eigen_hermv}, - {"eigns", eigen_nonsymm}, - {"eignsv", eigen_nonsymmv}, - {"schur", schur_decomp}, - {NULL, NULL} -}; - -static void -free_symm_cache (struct eigen_symm_cache *symm) -{ - gsl_eigen_symm_free (symm->ws ); - gsl_eigen_symmv_free (symm->vws); - gsl_vector_free (symm->diag); - symm->n = 0; -} - -static void -free_herm_cache (struct eigen_herm_cache *herm) -{ - gsl_eigen_herm_free (herm->ws ); - gsl_eigen_hermv_free (herm->vws); - gsl_vector_complex_free (herm->diag); - herm->n = 0; -} - -static void -free_nonsymm_cache (struct eigen_nonsymm_cache *nonsymm) -{ - gsl_eigen_nonsymm_free (nonsymm->ws ); - gsl_eigen_nonsymmv_free (nonsymm->vws); - gsl_matrix_free (nonsymm->mcopy); - nonsymm->n = 0; -} - -int -eigen_cache_free (lua_State *L) -{ - struct eigen_cache *cache = luaL_checkudata (L, 1, EIGEN_CACHE_MT_NAME); - - if (cache->symm->ws) - free_symm_cache (cache->symm); - - if (cache->herm->ws) - free_herm_cache (cache->herm); - - if (cache->nonsymm->ws) - free_nonsymm_cache (cache->nonsymm); - - return 0; -} - -void -eigen_push_cache (lua_State *L) -{ - struct eigen_cache *cache; - - cache = lua_newuserdata (L, sizeof(struct eigen_cache)); - - luaL_getmetatable (L, EIGEN_CACHE_MT_NAME); - lua_setmetatable (L, -2); - - cache->symm->ws = NULL; - cache->symm->n = 0; - - cache->herm->ws = NULL; - cache->herm->n = 0; - - cache->nonsymm->ws = NULL; - cache->nonsymm->n = 0; -} - -struct eigen_symm_cache * -eigen_cache_symm_set (lua_State *L, size_t n) -{ - struct eigen_cache *cache; - - lua_getfield(L, LUA_ENVIRONINDEX, "cache"); - cache = lua_touserdata (L, -1); - lua_pop (L, 1); - - assert (cache != NULL); - - if (!cache->symm->ws || cache->symm->n != n) - { - if (cache->symm->ws) - free_symm_cache (cache->symm); - - cache->symm->ws = gsl_eigen_symm_alloc (n); - cache->symm->vws = gsl_eigen_symmv_alloc (n); - cache->symm->diag = gsl_vector_alloc (n); - cache->symm->n = n; - } - - return cache->symm; -} - -struct eigen_herm_cache * -eigen_cache_herm_set (lua_State *L, size_t n) -{ - struct eigen_cache *cache; - - lua_getfield(L, LUA_ENVIRONINDEX, "cache"); - cache = lua_touserdata (L, -1); - lua_pop (L, 1); - - assert (cache != NULL); - - if (!cache->herm->ws || cache->herm->n != n) - { - if (cache->herm->ws) - free_herm_cache (cache->herm); - - cache->herm->ws = gsl_eigen_herm_alloc (n); - cache->herm->vws = gsl_eigen_hermv_alloc (n); - cache->herm->diag = gsl_vector_complex_alloc (n); - cache->herm->n = n; - } - - return cache->herm; -} - -struct eigen_nonsymm_cache * -eigen_cache_nonsymm_set (lua_State *L, size_t n) -{ - struct eigen_cache *cache; - - lua_getfield(L, LUA_ENVIRONINDEX, "cache"); - cache = lua_touserdata (L, -1); - lua_pop (L, 1); - - assert (cache != NULL); - - if (!cache->nonsymm->ws || cache->nonsymm->n != n) - { - if (cache->nonsymm->ws) - free_nonsymm_cache (cache->nonsymm); - - cache->nonsymm->ws = gsl_eigen_nonsymm_alloc (n); - cache->nonsymm->vws = gsl_eigen_nonsymmv_alloc (n); - cache->nonsymm->mcopy = gsl_matrix_alloc (n, n); - cache->nonsymm->n = n; - } - - return cache->nonsymm; -} - -int -eigen_symm_raw (lua_State *L, int compute_evecs) -{ - gsl_matrix *m = matrix_check (L, 1); - size_t i, j, n = m->size1; - struct eigen_symm_cache *cache; - gsl_matrix *eval, *evec; - gsl_vector *diag; - gsl_vector_view eview; - - if (m->size1 != m->size2) - return luaL_typerror (L, 1, "real symmetric matrix"); - - cache = eigen_cache_symm_set (L, n); - - eval = matrix_push_raw (L, n, 1); - eview = gsl_matrix_column (eval, 0); - - if (compute_evecs) - evec = matrix_push_raw (L, n, n); - - diag = cache->diag; - - for (j = 0; j < n; j++) - gsl_vector_set (diag, j, gsl_matrix_get (m, j, j)); - - if (compute_evecs) - gsl_eigen_symmv (m, &eview.vector, evec, cache->vws); - else - gsl_eigen_symm (m, &eview.vector, cache->ws); - - for (j = 0; j < n; j++) - gsl_matrix_set (m, j, j, gsl_vector_get (diag, j)); - - for (i = 0; i < n; i++) - for (j = 0; j < i; j++) - gsl_matrix_set (m, i, j, gsl_matrix_get (m, j, i)); - - return (compute_evecs ? 2 : 1); -} - -int -eigen_symm (lua_State *L) -{ - return eigen_symm_raw (L, 0); -} - -int -eigen_symmv (lua_State *L) -{ - return eigen_symm_raw (L, 1); -} - -static int -eigen_herm_raw (lua_State *L, int compute_evecs) -{ - gsl_matrix_complex *m = matrix_complex_check (L, 1); - size_t i, j, n = m->size1; - struct eigen_herm_cache *cache; - gsl_matrix *eval; - gsl_matrix_complex *evec; - gsl_vector_complex *diag; - gsl_vector_view eview; - - if (m->size1 != m->size2) - return luaL_typerror (L, 1, "real hermitian matrix"); - - cache = eigen_cache_herm_set (L, n); - - eval = matrix_push_raw (L, n, 1); - eview = gsl_matrix_column (eval, 0); - - if (compute_evecs) - evec = matrix_complex_push_raw (L, n, n); - - diag = cache->diag; - - for (j = 0; j < n; j++) - gsl_vector_complex_set (diag, j, gsl_matrix_complex_get (m, j, j)); - - if (compute_evecs) - gsl_eigen_hermv (m, &eview.vector, evec, cache->vws); - else - gsl_eigen_herm (m, &eview.vector, cache->ws); - - for (j = 0; j < n; j++) - gsl_matrix_complex_set (m, j, j, gsl_vector_complex_get (diag, j)); - - for (i = 0; i < n; i++) - for (j = 0; j < i; j++) - { - gsl_complex v = gsl_matrix_complex_get (m, j, i); - gsl_complex vc = gsl_complex_conjugate (v); - gsl_matrix_complex_set (m, i, j, vc); - } - - return (compute_evecs ? 2 : 1); -} - -int -eigen_herm (lua_State *L) -{ - return eigen_herm_raw (L, 0); -} - -int -eigen_hermv (lua_State *L) -{ - return eigen_herm_raw (L, 1); -} - -static int -eigen_nonsymm_raw (lua_State *L, int compute_evec) -{ - gsl_matrix *m = matrix_check (L, 1), *mcopy; - size_t n = m->size1; - struct eigen_nonsymm_cache *cache; - gsl_matrix_complex *eval, *evec; - gsl_vector_complex_view eview; - int retnval = 1; - int status; - - if (m->size1 != m->size2) - return luaL_typerror (L, 1, "real matrix"); - - cache = eigen_cache_nonsymm_set (L, n); - - eval = matrix_complex_push_raw (L, n, 1); - eview = gsl_matrix_complex_column (eval, 0); - - mcopy = cache->mcopy; - - gsl_matrix_memcpy (mcopy, m); - - if (compute_evec) - { - evec = matrix_complex_push_raw (L, n, n); - retnval ++; - status = gsl_eigen_nonsymmv (mcopy, &eview.vector, evec, cache->vws); - } - else - { - gsl_eigen_nonsymm_params (0, 0, cache->ws); - status = gsl_eigen_nonsymm (mcopy, &eview.vector, cache->ws); - } - - - if (status) - { - return luaL_error (L, "error during non-symmetric eigenvalues" - " determination: %s", - gsl_strerror (status)); - } - - return retnval; -} - -int -eigen_nonsymm (lua_State *L) -{ - return eigen_nonsymm_raw (L, 0); -} - -int -eigen_nonsymmv (lua_State *L) -{ - return eigen_nonsymm_raw (L, 1); -} - -int -schur_decomp (lua_State *L) -{ - gsl_matrix *m = matrix_check (L, 1); - size_t n = m->size1; - struct eigen_nonsymm_cache *cache; - gsl_matrix_complex *eval, *evec; - gsl_vector_complex_view eview; - gsl_matrix *t, *z; - int status; - - if (m->size1 != m->size2) - return luaL_typerror (L, 1, "real matrix"); - - cache = eigen_cache_nonsymm_set (L, n); - - eval = matrix_complex_push_raw (L, n, 1); - eview = gsl_matrix_complex_column (eval, 0); - - evec = matrix_complex_push_raw (L, n, n); - - t = matrix_push_raw (L, n, n); - gsl_matrix_memcpy (t, m); - - z = matrix_push_raw (L, n, n); - - status = gsl_eigen_nonsymmv_Z (t, &eview.vector, evec, z, cache->vws); - - if (status) - { - return luaL_error (L, "error during Schur decomposition: %s", - gsl_strerror (status)); - } - - gsl_linalg_hessenberg_set_zero (t); - - return 2; -} - -void -eigen_register (lua_State *L) -{ - luaL_newmetatable (L, EIGEN_CACHE_MT_NAME); - luaL_register (L, NULL, eigen_cache_methods); - lua_pop (L, 1); - - lua_newtable (L); - eigen_push_cache (L); - lua_setfield (L, -2, "cache"); - lua_replace (L, LUA_ENVIRONINDEX); - - luaL_register (L, NULL, eigen_functions); -} diff --git a/eigen-systems.h b/eigen-systems.h deleted file mode 100644 index 04c6433a..00000000 --- a/eigen-systems.h +++ /dev/null @@ -1,6 +0,0 @@ -#ifndef EIGEN_SYSTEM_H -#define EIGEN_SYSTEM_H - -extern void eigen_register (lua_State *L); - -#endif diff --git a/examples/nlinfit-jit.lua b/examples/nlinfit-jit.lua deleted file mode 100644 index b434f069..00000000 --- a/examples/nlinfit-jit.lua +++ /dev/null @@ -1,57 +0,0 @@ - -require 'cmatrix' - -use 'math' - -function demo1() - local n = 40 - - local yrf, sigrf - - local fdf = function(x, f, J) - for i=0, n-1 do - local A, lambda, b = x[0], x[1], x[2] - local t, y, sig = i, yrf[i], sigrf[i] - local e = exp(- lambda * t) - if f then f[i] = (A*e+b - y)/sig end - if J then - J:set(i, 0, e / sig) - J:set(i, 1, - t * A * e / sig) - J:set(i, 2, 1 / sig) - end - end - end - - local model = function(x, t) - local A, lambda, b = x[0], x[1], x[2] - return A * exp(- lambda * t) + b - end - - local xref = matrix.vec {5, 0.1, 1} - - local r = gsl.rng('mt19937') - r:set(0) - - yrf = matrix.new(n, 1, function(i) return model(xref, i) + gsl.rnd.gaussian(r, 0.1) end) - sigrf = matrix.new(n, 1, function() return 0.1 end) - - local s = gsl.nlinfit {n= n, p= 3} - - s:set(fdf, matrix.vec {1, 0, 0}) - print(s.x, s.chisq) - - for i=1, 10 do - s:iterate() - print('ITER=', i, ': ', s.x, s.chisq) - if s:test(0, 1e-8) then break end - end - - local p = graph.plot('Non-linear fit example') - local pts = graph.ipath(gsl.sequence(function(i) return i, yrf[i] end, 0, n-1)) - local fitln = graph.fxline(function(t) return model(s.x, t) end, 0, n-1) - p:addline(pts, 'blue', {{'marker', size=4}}) - p:addline(fitln) - p:show() -end - -echo 'demo1() - Simple non-linear fit example' diff --git a/examples/nlinfit.lua b/examples/nlinfit.lua index 8b412111..09c32f0f 100644 --- a/examples/nlinfit.lua +++ b/examples/nlinfit.lua @@ -1,195 +1,55 @@ - -- Non-linear Fit Examples / nlinfit.lua - -- - -- Copyright (C) 2009, 2010 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. +use 'math' function demo1() - local n, I = 50, complex.I - local p = {a= (-1+4*I) * 4, phi= 0.23, A= 0.55} - local y = cnew(n, 1, |i,j| p.A * complex.exp(p.a * (i-1)/n + I * p.phi)) + local n = 40 - local function cexpf(x, f, J) - for k=1, n do - local t, y = (k-1)/n, y[k] - local A, a, phi = x[1], x[2] + I * x[3], x[4] - local e = complex.exp(a * t + I * phi) - if f then f[k] = A * e - y end - if J then - J:set(k, 1, e) - J:set(k, 2, t * A * e) - J:set(k, 3, I * t * A * e) - J:set(k, 4, I * A * e) - end - end - end - - local function print_state(s) - print ("x: ", tr(s.p)) - print ("chi square: ", (hc(s.f) * s.f)[1]) - end - - s = cnlfsolver {fdf= cexpf, n= n, p0= vector {2.1, -2.8, 18, 0}} - repeat - print_state (s) - local status = s:iterate() - until status ~= 'continue' - print_state (s) -end - -function demo2() - local n = 50 - local px = vector {1.55, -1.1, 12.5} - local p0 = vector {2.5, -1.5, 5.3} - local xs = |i| (i-1)/n - local r = rng() + local yrf, sigrf - local fmodel = function(p, t, J) - local e, s = exp(p[2] * t), sin(p[3] * t) + local fdf = function(x, f, J) + for i=0, n-1 do + local A, lambda, b = x[0], x[1], x[2] + local t, y, sig = i, yrf[i], sigrf[i] + local e = exp(- lambda * t) + if f then f[i] = (A*e+b - y)/sig end if J then - J:set(1,1, e * s) - J:set(1,2, t * p[1] * e * s) - J:set(1,3, t * p[1] * e * cos(p[3] * t)) + J:set(i, 0, e / sig) + J:set(i, 1, - t * A * e / sig) + J:set(i, 2, 1 / sig) end - return p[1] * e * s end + end - local y = new(n, 1, |i,j| fmodel(px, xs(i)) * (1 + rnd.gaussian(r, 0.1))) - local x = new(n, 1, |i,j| xs(i)) - - local function expf(x, f, J) - for k=1, n do - local ym = fmodel(x, xs(k), J and J:row(k)) - if f then f[k] = ym - y[k] end - end - end - - local pl = plot('Non-linear fit / A * exp(a t) sin(w t)') - pl:add(xyline(x, y), 'blue', {{'stroke'}, {'marker', size= 5, mark="triangle"}}) - - local function print_state(s) - print ("x: ", tr(s.p)) - print ("chi square: ", prod(s.f, s.f)[1]) - end - - s = nlfsolver {fdf= expf, n= n, p0= p0} - - pl:addline(fxline(|x| fmodel(s.p, x), 0, xs(n)), 'red', {{'dash', 7, 3, 3, 3}}) - - repeat - print_state (s) - local status = s:iterate() - until status ~= 'continue' - print_state (s) - - pl:addline(fxline(|x| fmodel(s.p, x), 0, xs(n)), 'red') - pl:show() - - return pl -end + local model = function(x, t) + local A, lambda, b = x[0], x[1], x[2] + return A * exp(- lambda * t) + b + end -function demo2bis() - local n = 50 - local px = vector {1.55, -1.1, 12.5} - local p0 = vector {2.5, -1.5, 5.3} - local xs = |i| (i-1)/n - local r = rng() + local xref = matrix.vec {5, 0.1, 1} - local fmodel = function(p, t, J) - local e, s = exp(p[2] * t), sin(p[3] * t) - if J then - J:set(1,1, e * s) - J:set(1,2, t * p[1] * e * s) - J:set(1,3, t * p[1] * e * cos(p[3] * t)) - end - return p[1] * e * s - end + local r = gsl.rng('mt19937') + r:set(0) - local y = new(n, 1, |i,j| fmodel(px, xs(i)) * (1 + rnd.gaussian(r, 0.1))) - local x = new(n, 1, |i,j| xs(i)) + yrf = matrix.new(n, 1, function(i) return model(xref, i) + gsl.rnd.gaussian(r, 0.1) end) + sigrf = matrix.new(n, 1, function() return 0.1 end) - local function expf(x, f, J) - for k=1, n do - local ym = fmodel(x, xs(k), J and J:row(k)) - if f then f[k] = ym - y[k] end - end - end + local s = gsl.nlinfit {n= n, p= 3} - pl = plot('Non-linear fit / A * exp(a t) sin(w t)') - pl:addline(xyline(x, y), 'blue', {{'marker', size= 5}}) - pl:show() - pl.sync = false - pl:pushlayer() + s:set(fdf, matrix.vec {1, 0, 0}) + print(s.x, s.chisq) - local function print_state(s) - print ("x: ", tr(s.p)) - print ("chi square: ", prod(s.f, s.f)[1]) + for i=1, 10 do + s:iterate() + print('ITER=', i, ': ', s.x, s.chisq) + if s:test(0, 1e-8) then break end end - s = nlfsolver {fdf= expf, n= n, p0= p0} - - repeat - print_state (s) - pl:clear() - pl:addline(fxline(|x| fmodel(s.p, x), 0, xs(n))) - pl:flush() - io.read('*l') - local status = s:iterate() - until status ~= 'continue' - print_state (s) - - return pl -end - -function demo3() - -- This demo does the same things of demo2 but using the - -- higher level function 'nlinfit' - local px = vector {1.55, -1.1, 12.5} - local p0 = vector {2.5, -1.5, 5.3} - local n = 50 - local xs = |i| (i-1)/n - local r = rng() - - local f = function(p, x, J) - local e, s = exp(p[2] * x), sin(p[3] * x) - if J then - J:set(1,1, e * s) - J:set(1,2, x * p[1] * e * s) - J:set(1,3, x * p[1] * e * cos(p[3] * x)) - end - return p[1] * e * s - end - - local y = new(n, 1, |i,j| f(px, xs(i)) * (1 + rnd.gaussian(r, 0.1))) - local x = new(n, 1, |i,j| xs(i)) - - local fit, pr = nlinfit(f, x, y, p0) - - print('Fit result:', tr(pr)) - - pl = plot('Non-linear fit / A * exp(a t) sin(w t)') - pl:addline(xyline(x, y), 'blue', {{'marker', size= 4}}) - - pl:addline(fxline(|x| f(p0, x), 0, xs(n)), 'red', {{'dash', 7, 3, 3, 3}}) - pl:addline(fxline(fit, 0, xs(n)), 'red') - pl:show() - - return pl + local p = graph.plot('Non-linear fit example') + local pts = graph.ipath(gsl.sequence(function(i) return i, yrf[i] end, 0, n-1)) + local fitln = graph.fxline(function(t) return model(s.x, t) end, 0, n-1) + p:addline(pts, 'blue', {{'marker', size=4}}) + p:addline(fitln) + p:show() end -echo 'demo1() - examples on non-linear fit of complex data' -echo 'demo2() - examples on non-linear fit of real data and plots' -echo 'demo3() - the same of demo2() using a slightly different approach' +echo 'demo1() - Simple non-linear fit example' diff --git a/examples/ode-jit.lua b/examples/ode-jit.lua index 1cc3cd37..48fc5081 100644 --- a/examples/ode-jit.lua +++ b/examples/ode-jit.lua @@ -36,7 +36,7 @@ local function xyodeplot(f, t0, t1, x0, y0, h0, tsmp) while s.t < t do evol(s, f, t) end - ln:line_to(s.y[1], s.y[2]) + ln:line_to(s.y[0], s.y[1]) end local p = graph.plot('ODE integration example') @@ -57,14 +57,14 @@ local function modeplot(s, f, t0, y0, t1, tsmp) evol(s, f, t) end for k=1, n do - t[k]:line_to(s.t, s.y[k]) + t[k]:line_to(s.t, s.y[k-1]) end end else while s.t < t1 do evol(s, f, t1) for k=1, n do - t[k]:line_to(s.t, s.y[k]) + t[k]:line_to(s.t, s.y[k-1]) end end end @@ -115,7 +115,7 @@ function demo3() while s.t < t do s:evolve(f, t) end - ln:line_to(s.t, s.y[1]) + ln:line_to(s.t, s.y[0]) end local p = graph.plot() p:addline(ln) diff --git a/fdfmultimin.c b/fdfmultimin.c deleted file mode 100644 index afd96030..00000000 --- a/fdfmultimin.c +++ /dev/null @@ -1,411 +0,0 @@ - -/* fdfmultimin.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 <string.h> -#include <gsl/gsl_multimin.h> -#include <gsl/gsl_vector.h> -#include <gsl/gsl_matrix.h> - -#include "gs-types.h" -#include "multimin.h" -#include "random.h" -#include "lua-utils.h" -#include "matrix.h" - -struct params { - lua_State *L; - size_t n; - int type_error; - int args_index; -}; - -struct fdfmultimin { - gsl_multimin_fdfminimizer *s; - gsl_multimin_function_fdf fdf[1]; - struct params p[1]; -}; - -enum fenvidx { - FENV_FUNCTION = 1, - FENV_X = 2, - FENV_GRAD = 3 -}; - -static int fdfmultimin_free (lua_State *L); -static int fdfmultimin_set (lua_State *L); -static int fdfmultimin_run (lua_State *L); -static int fdfmultimin_step (lua_State *L); -static int fdfmultimin_index (lua_State *L); - -static int fdfmultimin_get_x (lua_State *L); -static int fdfmultimin_get_value (lua_State *L); -static int fdfmultimin_get_gradient (lua_State *L); - -static double fdfmultimin_f_hook (const gsl_vector * x, void * p); -static void fdfmultimin_df_hook (const gsl_vector * x, void * p, gsl_vector *G); -static void fdfmultimin_fdf_hook (const gsl_vector * x, void * p, double *f, - gsl_vector *G); - -static const struct luaL_Reg fdfmultimin_properties[] = { - {"x", fdfmultimin_get_x}, - {"value", fdfmultimin_get_value}, - {"gradient", fdfmultimin_get_gradient}, - {NULL, NULL} -}; - -static const struct luaL_Reg fdfmultimin_methods[] = { - {"step", fdfmultimin_step}, - {"run", fdfmultimin_run}, - {"set", fdfmultimin_set}, - {NULL, NULL} -}; - -const struct luaL_Reg fdfmultimin_metatable[] = { - {"__gc", fdfmultimin_free}, - {"__index", fdfmultimin_index}, - {NULL, NULL} -}; - -static const gsl_multimin_fdfminimizer_type * -fdfmin_algo_lookup (const char *name) -{ - if (strcmp (name, "conjugate_fr") == 0) - return gsl_multimin_fdfminimizer_conjugate_fr; - else if (strcmp (name, "conjugate_pr") == 0) - return gsl_multimin_fdfminimizer_conjugate_pr; - else if (strcmp (name, "bfgs") == 0) - return gsl_multimin_fdfminimizer_vector_bfgs2; - else if (strcmp (name, "steepest_descent") == 0) - return gsl_multimin_fdfminimizer_steepest_descent; - return NULL; -} - -static const char * const fdfmin_algo_list = "conjugate_fr, conjugate_pr, bfgs, steepest_descent"; - -#define MULTIMIN_MAX_ITER 20 - -struct fdfmultimin * -push_new_fdf_multimin (lua_State *L, int findex, size_t n, - const gsl_multimin_fdfminimizer_type *T) -{ - struct fdfmultimin *m = lua_newuserdata (L, sizeof (struct fdfmultimin)); - m->s = gsl_multimin_fdfminimizer_alloc (T, n); - - if (m->s == NULL) - luaL_error (L, OUT_OF_MEMORY_MSG); - - luaL_getmetatable (L, GS_METATABLE(GS_FDFMULTIMIN)); - lua_setmetatable (L, -2); - - lua_newtable (L); - - lua_pushvalue (L, findex); - lua_rawseti (L, -2, FENV_FUNCTION); - - matrix_push_raw (L, n, 1); - lua_rawseti (L, -2, FENV_X); - - matrix_push_view (L, NULL); - lua_rawseti (L, -2, FENV_GRAD); - - lua_setfenv (L, -2); - - return m; -} - -static struct fdfmultimin * -check_fdf_multimin (lua_State *L, int index) -{ - return luaL_checkudata (L, index, GS_METATABLE(GS_FDFMULTIMIN)); -} - -static struct fdfmultimin * -check_init_fdf_multimin (lua_State *L, int index) -{ - struct fdfmultimin *m = luaL_checkudata (L, index, GS_METATABLE(GS_FDFMULTIMIN)); - if (m->fdf->params == NULL) - luaL_error (L, "minimizer is not initialised"); - return m; -} - -int -fdfmultimin_new (lua_State *L) -{ - const gsl_multimin_fdfminimizer_type *T; - struct fdfmultimin *m; - int ni = luaL_checkinteger (L, 2); - int findex = 1; - size_t n; - - if (ni <= 0) - luaL_error (L, "argument #1 should be a positive integer"); - else - n = (size_t) ni; - - if (!lua_isfunction (L, findex)) - luaL_typerror (L, findex, "function"); - - if (lua_gettop (L) > 2) - { - const char *req_algo = luaL_checkstring (L, 3); - T = fdfmin_algo_lookup (req_algo); - if (T == NULL) - return luaL_error (L, "unknown algorithm. Known algorithm are: %s", - fdfmin_algo_list); - } - else - { - T = gsl_multimin_fdfminimizer_conjugate_fr; - } - - m = push_new_fdf_multimin (L, findex, n, T); - - m->fdf->n = n; - m->fdf->f = fdfmultimin_f_hook; - m->fdf->df = fdfmultimin_df_hook; - m->fdf->fdf = fdfmultimin_fdf_hook; - m->fdf->params = NULL; - - return 1; -} - -int -fdfmultimin_free (lua_State *L) -{ - struct fdfmultimin *m = check_fdf_multimin (L, 1); - lua_getfenv (L, 1); - lua_rawgeti (L, 2, FENV_GRAD); - matrix_null_view (L, -1); - gsl_multimin_fdfminimizer_free (m->s); - return 0; -} - -void -fdfmultimin_fdf_hook (const gsl_vector * x, void * _p, double *f, gsl_vector *G) -{ - struct params *p = _p; - lua_State *L = p->L; - size_t n = p->n, k; - gsl_matrix *lx; - int aindex = p->args_index; - int nargs = 1; - - lua_pushvalue (L, aindex); - - lx = matrix_check (L, aindex+1); - for (k = 0; k < n; k++) - lx->data[k] = x->data[k]; - lua_pushvalue (L, aindex+1); - - if (G) - { - matrix_set_view_and_push (L, aindex+2, G->data, n, 1, NULL); - nargs ++; - } - - lua_call (L, nargs, 1); - - if (lua_isnumber (L, -1)) - *f = lua_tonumber (L, -1); - else - { - *f = GSL_NAN; - p->type_error = 1; - } - - lua_pop (L, 1); -} - -double -fdfmultimin_f_hook (const gsl_vector * x, void * params) -{ - double f; - fdfmultimin_fdf_hook (x, params, &f, NULL); - return f; -} - -void -fdfmultimin_df_hook (const gsl_vector * x, void * params, gsl_vector *G) -{ - double f; - fdfmultimin_fdf_hook (x, params, &f, G); -} - -static void -fcall_args_prepare (lua_State *L, struct params *p, int index) -{ - int bindex; - lua_getfenv (L, index); - - bindex = lua_gettop (L); - lua_rawgeti (L, bindex, FENV_FUNCTION); - lua_rawgeti (L, bindex, FENV_X); - lua_rawgeti (L, bindex, FENV_GRAD); - lua_remove (L, bindex); - - p->args_index = bindex; -} - -int -fdfmultimin_set (lua_State *L) -{ - struct fdfmultimin *m = check_fdf_multimin (L, 1); - gsl_matrix *x0m = matrix_check (L, 2); - double step_size = gs_check_number (L, 3, FP_CHECK_NORMAL); - double tol = 0.1; - gsl_vector_view x0v; - size_t n = m->fdf->n; - int status; - - if (lua_gettop (L) > 3) - tol = gs_check_number (L, 4, FP_CHECK_NORMAL); - - lua_pushcfunction (L, gradient_auto_check); - mlua_fenv_get (L, 1, FENV_FUNCTION); - lua_pushvalue (L, 2); - lua_pushnumber (L, step_size); - - lua_call (L, 3, 0); - - lua_settop (L, 2); /* get rid of the step_size */ - - m->fdf->params = m->p; - - m->p->L = L; - m->p->n = n; - m->p->type_error = 0; - - x0v = gsl_matrix_column (x0m, 0); - - mlua_null_cache (L, 1); - fcall_args_prepare (L, m->p, 1); - - status = gsl_multimin_fdfminimizer_set (m->s, m->fdf, &x0v.vector, step_size, tol); - - if (status != GSL_SUCCESS) - return luaL_error (L, "minimizer:set %s", gsl_strerror (status)); - - return 0; -} - -int -fdfmultimin_step (lua_State *L) -{ - struct fdfmultimin *m = check_fdf_multimin (L, 1); - int status; - - mlua_null_cache (L, 1); - fcall_args_prepare (L, m->p, 1); - - status = gsl_multimin_fdfminimizer_iterate (m->s); - - if (status != GSL_SUCCESS) - return luaL_error (L, "minimizer:step %s", gsl_strerror (status)); - - status = gsl_multimin_test_gradient (m->s->gradient, 1e-4); - - if (status == GSL_CONTINUE) - { - lua_pushstring (L, "continue"); - return 1; - } - else if (status == GSL_SUCCESS) - { - lua_pushstring (L, "success"); - return 1; - } - - return luaL_error (L, "minimizer:step %s", gsl_strerror (status)); -} - -int -fdfmultimin_run (lua_State *L) -{ - struct fdfmultimin *m = check_fdf_multimin (L, 1); - size_t iter = 0; - int status; - - mlua_null_cache (L, 1); - fcall_args_prepare (L, m->p, 1); - - do { - iter++; - - status = gsl_multimin_fdfminimizer_iterate (m->s); - - if (status) - break; - - status = gsl_multimin_test_gradient (m->s->gradient, 1e-4); - - if (status == GSL_SUCCESS) - return 0; - } - while (status == GSL_CONTINUE && iter < MULTIMIN_MAX_ITER); - - if (m->p->type_error) - return luaL_error (L, "function should return a real number"); - - return luaL_error (L, "minimizer:run %s", gsl_strerror (status)); -} - -int -fdfmultimin_get_gradient (lua_State *L) -{ - struct fdfmultimin *m = check_init_fdf_multimin (L, 1); - gsl_matrix *gm = matrix_push (L, m->p->n, 1); - gsl_vector *g = gsl_multimin_fdfminimizer_gradient (m->s); - size_t k; - for (k = 0; k < m->p->n; k++) - gm->data[k] = g->data[k]; - return 1; -} - -int -fdfmultimin_get_x (lua_State *L) -{ - struct fdfmultimin *m = check_init_fdf_multimin (L, 1); - gsl_matrix *xm = matrix_push (L, m->p->n, 1); - gsl_vector *x = gsl_multimin_fdfminimizer_x (m->s); - size_t k; - for (k = 0; k < m->p->n; k++) - xm->data[k] = x->data[k]; - return 1; -} - -int -fdfmultimin_get_value (lua_State *L) -{ - struct fdfmultimin *m = check_init_fdf_multimin (L, 1); - double v = gsl_multimin_fdfminimizer_minimum (m->s); - lua_pushnumber (L, v); - return 1; -} - -int -fdfmultimin_index (lua_State *L) -{ - return mlua_index_with_properties (L, - fdfmultimin_properties, - fdfmultimin_methods, - true); -} @@ -1,511 +0,0 @@ - -/* fft.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 <gsl/gsl_fft_real.h> -#include <gsl/gsl_fft_complex.h> -#include <gsl/gsl_fft_halfcomplex.h> - -#include "gs-types.h" -#include "lcomplex.h" -#include "matrix.h" -#include "cmatrix.h" -#include "fft.h" -#include "lua-utils.h" - -#define FFT_CACHE_MT_NAME "GSL.fftcache" - -struct fft_cache { - size_t size; - gsl_fft_real_workspace *ws; - gsl_fft_real_wavetable *rwt; - gsl_fft_halfcomplex_wavetable *hcwt; - - size_t csize; - gsl_fft_complex_workspace *cws; - gsl_fft_complex_wavetable *cwt; -}; - -struct fft_hc_sel { - int (*length)(gsl_matrix *); - int (*get_index)(size_t, int, int *, int *, int *); - void (*transform)(lua_State *, gsl_matrix *); -}; - -static gsl_matrix * fft_hc_check (lua_State *L, int index, - struct fft_hc_sel ** selptr); -static struct fft_cache * check_fft_cache_dim (lua_State *L, size_t n, - bool want_complex); - -static int fft_hc_length (lua_State *L); -static int fft_hc_get (lua_State *L); -static int fft_hc_set (lua_State *L); -static int fft_hc_free (lua_State *L); -static int fft_hc_index (lua_State *L); -static int fft_real (lua_State *L); -static int fft_complex (lua_State *L); -static int fft_real_inverse (lua_State *L); - -static int fft_cache_free (lua_State *L); - - -static int fft_hc_mixed_radix_length (gsl_matrix *v); -static int fft_hc_mixed_radix_get_index (size_t n, int index, - int *rindex, int *cindex, int *csign); -static void fft_hc_mixed_radix_transform (lua_State *L, gsl_matrix *hc); -static int fft_hc_radix2_length (gsl_matrix *v); -static int fft_hc_radix2_get_index (size_t n, int index, - int *rindex, int *cindex, int *csign); -static void fft_hc_radix2_transform (lua_State *L, gsl_matrix *hc); - -static struct fft_hc_sel fft_hc_radix2_sel[1] = {{ - .length = fft_hc_radix2_length, - .get_index = fft_hc_radix2_get_index, - .transform = fft_hc_radix2_transform, - }}; - -static struct fft_hc_sel fft_hc_mixed_radix_sel[1] = {{ - .length = fft_hc_mixed_radix_length, - .get_index = fft_hc_mixed_radix_get_index, - .transform = fft_hc_mixed_radix_transform, - }}; - -static const struct luaL_Reg fft_hc_metatable[] = { - {"__gc", fft_hc_free}, - {"__index", fft_hc_index}, - {NULL, NULL} -}; - -static const struct luaL_Reg fft_hc_methods[] = { - {"get", fft_hc_get}, - {"set", fft_hc_set}, - {NULL, NULL} -}; - -static const struct luaL_Reg fft_hc_properties[] = { - {"length", fft_hc_length}, - {NULL, NULL} -}; - -static const struct luaL_Reg fft_cache_methods[] = { - {"__gc", fft_cache_free}, - {NULL, NULL} -}; - -static const struct luaL_Reg fft_functions[] = { - {"fft", fft_real}, - {"fft_inv", fft_real_inverse}, - {"cfft", fft_complex}, - {NULL, NULL} -}; - -static int -is_twopower (size_t n) -{ - for (; n > 0; n = n/2) - { - int r = n % 2; - if (r && n > 1) - return 0; - } - return 1; -} - -int -fft_hc_radix2_length (gsl_matrix *v) -{ - return v->size1 / 2; -}; - -int -fft_hc_radix2_get_index (size_t _n, int i, int *rindex, int *cindex, int *csign) -{ - const int n = (int) _n; - - if (i < -n/2+1 || i >= n/2+1) - return 1; - - if (i > 0) - { - *rindex = i; - *cindex = n-i; - *csign = (i == n/2 ? 0 : 1); - } - else if (i < 0) - { - *rindex = -i; - *cindex = n+i; - *csign = -1; - } - else - { - *rindex = 0; - *csign = 0; - } - - return 0; -}; - -void -fft_hc_radix2_transform (lua_State *L, gsl_matrix *hc) -{ - gsl_fft_halfcomplex_radix2_inverse (hc->data, hc->tda, hc->size1); -} - -int -fft_hc_mixed_radix_length (gsl_matrix *v) -{ - return v->size1 / 2; -}; - -int -fft_hc_mixed_radix_get_index (size_t _n, int i, - int *rindex, int *cindex, int *csign) -{ - int n = (int) _n; - int is = (n % 2 == 0 ? -n/2 + 1 : -(n-1)/2); - - if (i < is || i >= is + n) - return 1; - - if (i > 0) - { - *rindex = 2*i-1; - *cindex = 2*i; - *csign = (n % 2 == 0 && i == n/2 ? 0 : 1); - } - else if (i < 0) - { - *rindex = -2*i-1; - *cindex = -2*i; - *csign = -1; - } - else - { - *rindex = 0; - *csign = 0; - } - - return 0; -}; - -void -fft_hc_mixed_radix_transform (lua_State *L, gsl_matrix *hc) -{ - const size_t n = hc->size1; - struct fft_cache *cache = check_fft_cache_dim (L, n, false); - gsl_fft_halfcomplex_transform (hc->data, hc->tda, n, cache->hcwt, cache->ws); - gsl_matrix_scale (hc, 1/(double)n); -} - -gsl_matrix * -fft_hc_check (lua_State *L, int index, struct fft_hc_sel ** selptr) -{ - gsl_matrix *p; - - p = gs_is_userdata (L, index, GS_HALFCMPL_R2); - - if (p != NULL) - { - if (selptr) - *selptr = fft_hc_radix2_sel; - return p; - } - - p = gs_is_userdata (L, index, GS_HALFCMPL_MR); - - if (p != NULL) - { - if (selptr) - *selptr = fft_hc_mixed_radix_sel; - return p; - } - else - { - const char *msg = lua_pushfstring(L, "%s or %s", - type_qualified_name (GS_HALFCMPL_R2), - type_qualified_name (GS_HALFCMPL_MR)); - gs_type_error (L, index, msg); - } - - return 0; -} - -int -fft_hc_length (lua_State *L) -{ - struct fft_hc_sel *sel; - gsl_matrix *hc = fft_hc_check (L, 1, &sel); - lua_pushnumber (L, sel->length (hc)); - return 1; -} - -int -fft_hc_get (lua_State *L) -{ - struct fft_hc_sel *sel; - gsl_matrix *hc = fft_hc_check (L, 1, &sel); - int hcindex = lua_tonumber (L, 2); - size_t n = hc->size1; - lua_Complex r; - int rindex, cindex; - int csign; - int status; - - status = sel->get_index (n, hcindex, &rindex, &cindex, &csign); - if (status) - return luaL_error (L, "index out of range"); - - if (csign != 0) - r = gsl_matrix_get (hc, rindex, 0) + csign * gsl_matrix_get (hc, cindex, 0) * I; - else - r = gsl_matrix_get (hc, rindex, 0); - - lua_pushcomplex (L, r); - return 1; -} - -int -fft_hc_set (lua_State *L) -{ - struct fft_hc_sel *sel; - gsl_matrix *hc = fft_hc_check (L, 1, &sel); - int hcindex = lua_tonumber (L, 2); - lua_Complex val = lua_tocomplex (L, 3); - size_t n = hc->size1; - int rindex, cindex; - int csign; - int status; - - status = sel->get_index (n, hcindex, &rindex, &cindex, &csign); - if (status) - return luaL_error (L, "index out of range"); - - if (csign != 0) - { - gsl_matrix_set (hc, rindex, 0, creal(val)); - gsl_matrix_set (hc, cindex, 0, csign * cimag(val)); - } - else - { - if (cimag(val) != 0) - return luaL_error (L, "imaginary part should be 0 for this term"); - gsl_matrix_set (hc, rindex, 0, creal(val)); - } - - return 0; -} - -struct fft_cache * -check_fft_cache_dim (lua_State *L, size_t n, bool want_complex) -{ - struct fft_cache *cache; - - lua_getfield(L, LUA_ENVIRONINDEX, "cache"); - cache = lua_touserdata (L, -1); - lua_pop (L, 1); - - assert (cache != NULL); - - if (want_complex) - { - if (cache->cws && cache->csize == n) - return cache; - - if (cache->cws) - { - gsl_fft_complex_workspace_free (cache->cws); - gsl_fft_complex_wavetable_free (cache->cwt); - } - - cache->cws = gsl_fft_complex_workspace_alloc (n); - cache->cwt = gsl_fft_complex_wavetable_alloc (n); - cache->csize = n; - } - else - { - if (cache->ws && cache->size == n) - return cache; - - if (cache->ws) - { - gsl_fft_real_workspace_free (cache->ws); - gsl_fft_real_wavetable_free (cache->rwt); - gsl_fft_halfcomplex_wavetable_free (cache->hcwt); - } - - cache->ws = gsl_fft_real_workspace_alloc (n); - cache->rwt = gsl_fft_real_wavetable_alloc (n); - cache->hcwt = gsl_fft_halfcomplex_wavetable_alloc (n); - cache->size = n; - } - - return cache; -} - -int -fft_cache_free (lua_State *L) -{ - struct fft_cache *cache = luaL_checkudata (L, 1, FFT_CACHE_MT_NAME); - - if (cache->ws) - { - gsl_fft_real_workspace_free (cache->ws); - gsl_fft_real_wavetable_free (cache->rwt); - gsl_fft_halfcomplex_wavetable_free (cache->hcwt); - cache->size = 0; - } - - if (cache->cws) - { - gsl_fft_complex_workspace_free (cache->cws); - gsl_fft_complex_wavetable_free (cache->cwt); - cache->csize = 0; - } - - return 0; -} - -int -fft_real (lua_State *L) -{ - gsl_matrix *v = matrix_check (L, 1); - size_t n = v->size1; - - if (v->size2 != 1) - luaL_error (L, "single column matrix expected"); - if (lua_gettop (L) > 1) - luaL_error (L, "single argument expected"); - - if (is_twopower (n)) - { - gsl_fft_real_radix2_transform (v->data, v->tda, n); - luaL_getmetatable (L, GS_METATABLE(GS_HALFCMPL_R2)); - lua_setmetatable (L, -2); - } - else - { - struct fft_cache *cache = check_fft_cache_dim (L, n, false); - gsl_fft_real_transform (v->data, v->tda, n, cache->rwt, cache->ws); - luaL_getmetatable (L, GS_METATABLE(GS_HALFCMPL_MR)); - lua_setmetatable (L, -2); - } - - return 0; -} - -int -fft_real_inverse (lua_State *L) -{ - struct fft_hc_sel *sel; - gsl_matrix *hc = fft_hc_check (L, 1, &sel); - sel->transform (L, hc); - luaL_getmetatable (L, GS_METATABLE(GS_MATRIX)); - lua_setmetatable (L, -2); - return 0; -} - -int -fft_hc_free (lua_State *L) -{ - gsl_matrix *m = fft_hc_check (L, 1, NULL); - assert (m->block); - gsl_block_free (m->block); - return 0; -} - -int -fft_hc_index (lua_State *L) -{ - return mlua_index_with_properties (L, fft_hc_properties, fft_hc_methods, false); -} - -int -fft_complex (lua_State *L) -{ - gsl_matrix_complex *v = matrix_complex_check (L, 1); - lua_Integer sign = luaL_optinteger (L, 2, -1); - size_t n = v->size1; - struct fft_cache *cache; - int csign; - - if (v->size2 != 1) - luaL_error (L, "single column matrix expected"); - - csign = (sign > 0 ? -1 : 1); - - cache = check_fft_cache_dim (L, n, true); - gsl_fft_complex_transform (v->data, v->tda, n, cache->cwt, cache->cws, csign); - - if (csign < 0) - { - gsl_complex ff = {{1/(double)n, 0}}; - gsl_matrix_complex_scale (v, ff); - } - - return 0; -} - -static void -fft_pushcache (lua_State *L) -{ - struct fft_cache *cache; - - cache = lua_newuserdata (L, sizeof(struct fft_cache)); - - luaL_getmetatable (L, FFT_CACHE_MT_NAME); - lua_setmetatable (L, -2); - - cache->ws = NULL; - cache->rwt = NULL; - cache->hcwt = NULL; - cache->size = 0; - - cache->cws = NULL; - cache->cwt = NULL; - cache->csize = 0; -} - - - -void -fft_register (lua_State *L) -{ - luaL_newmetatable (L, GS_METATABLE(GS_HALFCMPL_R2)); - luaL_register (L, NULL, fft_hc_metatable); - lua_setfield (L, -2, "FFT_hc_radix2"); - - luaL_newmetatable (L, GS_METATABLE(GS_HALFCMPL_MR)); - luaL_register (L, NULL, fft_hc_metatable); - lua_setfield (L, -2, "FFT_hc_mixed_radix"); - - luaL_newmetatable (L, FFT_CACHE_MT_NAME); - luaL_register (L, NULL, fft_cache_methods); - lua_pop (L, 1); - - lua_newtable (L); - fft_pushcache (L); - lua_setfield (L, -2, "cache"); - lua_replace (L, LUA_ENVIRONINDEX); - - luaL_register (L, NULL, fft_functions); -} @@ -1,7 +0,0 @@ - -#ifndef FFT_H -#define FFT_H - -extern void fft_register (lua_State *L); - -#endif diff --git a/fmultimin.c b/fmultimin.c deleted file mode 100644 index b06b672b..00000000 --- a/fmultimin.c +++ /dev/null @@ -1,340 +0,0 @@ - -/* fmultimin.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 <string.h> -#include <gsl/gsl_multimin.h> -#include <gsl/gsl_vector.h> -#include <gsl/gsl_matrix.h> -#include <gsl/gsl_deriv.h> -#include <gsl/gsl_rng.h> - -#include "gs-types.h" -#include "multimin.h" -#include "random.h" -#include "lua-utils.h" -#include "matrix.h" - -#include <lua.h> -#include <lauxlib.h> - -static int fmultimin_free (lua_State *L); -static int fmultimin_set (lua_State *L); -static int fmultimin_run (lua_State *L); -static int fmultimin_step (lua_State *L); -static int fmultimin_index (lua_State *L); - -static int fmultimin_get_x (lua_State *L); -static int fmultimin_get_value (lua_State *L); - -static double fmultimin_f_hook (const gsl_vector * x, void * p); - -struct params { - lua_State *L; - size_t n; - int type_error; - int args_index; -}; - -struct fmultimin { - gsl_multimin_fminimizer *s; - gsl_multimin_function f[1]; - struct params p[1]; - double size_tol; -}; - -enum fenvidx { - FENV_FUNCTION = 1, - FENV_X = 2, -}; - -static const struct luaL_Reg fmultimin_properties[] = { - {"x", fmultimin_get_x}, - {"value", fmultimin_get_value}, - {NULL, NULL} -}; - -static const struct luaL_Reg fmultimin_methods[] = { - {"step", fmultimin_step}, - {"run", fmultimin_run}, - {"set", fmultimin_set}, - {NULL, NULL} -}; - -const struct luaL_Reg fmultimin_metatable[] = { - {"__gc", fmultimin_free}, - {"__index", fmultimin_index}, - {NULL, NULL} -}; - -#define FMULTIMIN_MAX_ITER 20 - -struct fmultimin * -push_new_fmultimin (lua_State *L, int findex, size_t n, - const gsl_multimin_fminimizer_type *T) -{ - struct fmultimin *m = lua_newuserdata (L, sizeof (struct fmultimin)); - m->s = gsl_multimin_fminimizer_alloc (T, n); - - if (m->s == NULL) - luaL_error (L, OUT_OF_MEMORY_MSG); - - luaL_getmetatable (L, GS_METATABLE(GS_FMULTIMIN)); - lua_setmetatable (L, -2); - - lua_newtable (L); - - lua_pushvalue (L, findex); - lua_rawseti (L, -2, FENV_FUNCTION); - - matrix_push_raw (L, n, 1); - lua_rawseti (L, -2, FENV_X); - - lua_setfenv (L, -2); - - return m; -} - -static struct fmultimin * -check_fmultimin (lua_State *L, int index) -{ - return luaL_checkudata (L, index, GS_METATABLE(GS_FMULTIMIN)); -} - -static struct fmultimin * -check_init_fmultimin (lua_State *L, int index) -{ - struct fmultimin *m = luaL_checkudata (L, index, GS_METATABLE(GS_FMULTIMIN)); - if (m->f->params == NULL) - luaL_error (L, "minimizer is not initialised"); - return m; -} - -int -fmultimin_new (lua_State *L) -{ - struct fmultimin *m; - int ni = luaL_checkinteger (L, 2); - int findex = 1; - size_t n; - - if (ni <= 0) - return luaL_error (L, "argument #1 should be a positive integer"); - else - n = (size_t) ni; - - if (!lua_isfunction (L, findex)) - luaL_typerror (L, findex, "function"); - - m = push_new_fmultimin (L, findex, n, gsl_multimin_fminimizer_nmsimplex2); - - m->f->n = n; - m->f->f = fmultimin_f_hook; - m->f->params = NULL; - - return 1; -} - -int -fmultimin_free (lua_State *L) -{ - struct fmultimin *m = check_fmultimin (L, 1); - gsl_multimin_fminimizer_free (m->s); - return 0; -} - -double -fmultimin_f_hook (const gsl_vector * x, void * _p) -{ - struct params *p = _p; - lua_State *L = p->L; - size_t n = p->n, k; - gsl_matrix *lx; - int aindex = p->args_index; - - lua_pushvalue (L, aindex); - - lx = matrix_check (L, aindex+1); - for (k = 0; k < n; k++) - lx->data[k] = x->data[k]; - lua_pushvalue (L, aindex+1); - - lua_call (L, 1, 1); - - if (lua_isnumber (L, -1)) - return lua_tonumber (L, -1); - else - { - p->type_error = 1; - return GSL_NAN; - } - - lua_pop (L, 1); -} - -static void -fcall_args_prepare (lua_State *L, struct params *p, int index) -{ - int bindex; - lua_getfenv (L, index); - - bindex = lua_gettop (L); - lua_rawgeti (L, bindex, FENV_FUNCTION); - lua_rawgeti (L, bindex, FENV_X); - lua_remove (L, bindex); - - p->args_index = bindex; -} - -int -fmultimin_set (lua_State *L) -{ - struct fmultimin *m = check_fmultimin (L, 1); - gsl_matrix *x0m = matrix_check (L, 2); - gsl_matrix *stepm = matrix_check (L, 3); - double size_tol = gs_check_number (L, 4, FP_CHECK_NORMAL); - gsl_vector_view x0v, stepv; - size_t n = m->f->n; - int status; - - if (x0m->size2 > 1) - return gs_type_error (L, 2, "column vector"); - - if (stepm->size2 > 1) - return gs_type_error (L, 2, "column vector"); - - x0v = gsl_matrix_column (x0m, 0); - stepv = gsl_matrix_column (stepm, 0); - - m->f->params = m->p; - m->size_tol = size_tol; - - m->p->L = L; - m->p->n = n; - m->p->type_error = 0; - - lua_settop (L, 2); /* get rid of the step_size */ - - mlua_null_cache (L, 1); - fcall_args_prepare (L, m->p, 1); - - status = gsl_multimin_fminimizer_set (m->s, m->f, &x0v.vector, &stepv.vector); - - if (status != GSL_SUCCESS) - return luaL_error (L, "minimizer:set %s", gsl_strerror (status)); - - return 0; -} - -int -fmultimin_step (lua_State *L) -{ - struct fmultimin *m = check_init_fmultimin (L, 1); - double size; - int status; - - mlua_null_cache (L, 1); - fcall_args_prepare (L, m->p, 1); - - status = gsl_multimin_fminimizer_iterate (m->s); - - if (status != GSL_SUCCESS) - return luaL_error (L, "minimizer:step %s", gsl_strerror (status)); - - size = gsl_multimin_fminimizer_size (m->s); - status = gsl_multimin_test_size (size, m->size_tol); - - if (status == GSL_CONTINUE) - { - lua_pushstring (L, "continue"); - return 1; - } - else if (status == GSL_SUCCESS) - { - lua_pushstring (L, "success"); - return 1; - } - - return luaL_error (L, "minimizer:step %s", gsl_strerror (status)); -} - -int -fmultimin_run (lua_State *L) -{ - struct fmultimin *m = check_init_fmultimin (L, 1); - size_t iter = 0; - double size; - int status; - - mlua_null_cache (L, 1); - fcall_args_prepare (L, m->p, 1); - - do { - iter++; - - status = gsl_multimin_fminimizer_iterate (m->s); - - if (status) - break; - - size = gsl_multimin_fminimizer_size (m->s); - status = gsl_multimin_test_size (size, m->size_tol); - - if (status == GSL_SUCCESS) - return 0; - } - while (status == GSL_CONTINUE && iter < FMULTIMIN_MAX_ITER); - - if (m->p->type_error) - return luaL_error (L, "function should return a real number"); - - return luaL_error (L, "minimizer:run %s", gsl_strerror (status)); -} - -int -fmultimin_get_x (lua_State *L) -{ - struct fmultimin *m = check_init_fmultimin (L, 1); - gsl_matrix *xm = matrix_push (L, m->p->n, 1); - gsl_vector *x = gsl_multimin_fminimizer_x (m->s); - size_t k; - for (k = 0; k < m->p->n; k++) - xm->data[k] = x->data[k]; - return 1; -} - -int -fmultimin_get_value (lua_State *L) -{ - struct fmultimin *m = check_init_fmultimin (L, 1); - double v = gsl_multimin_fminimizer_minimum (m->s); - lua_pushnumber (L, v); - return 1; -} - -int -fmultimin_index (lua_State *L) -{ - return mlua_index_with_properties (L, - fmultimin_properties, - fmultimin_methods, - true); -} diff --git a/gradcheck.c b/gradcheck.c deleted file mode 100644 index 55d65769..00000000 --- a/gradcheck.c +++ /dev/null @@ -1,135 +0,0 @@ - -/* gradcheck.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_vector.h> -#include <gsl/gsl_matrix.h> -#include <gsl/gsl_deriv.h> -#include <gsl/gsl_rng.h> - -#include "multimin.h" -#include "matrix.h" -#include "random.h" - -struct grad_check { - lua_State *L; - int f_index, x_index; - double step; - size_t var_vector_index; - size_t n; -}; - -static double -f_check_hook (double x, void *_p) -{ - struct grad_check *p = _p; - lua_State *L = p->L; - gsl_matrix *xm = matrix_check (L, p->x_index); - gsl_matrix_set (xm, p->var_vector_index, 0, x); - lua_pushvalue (L, p->f_index); - lua_pushvalue (L, p->x_index); - lua_call (p->L, 1, 1); - - if (lua_isnumber (L, -1)) - { - double v = lua_tonumber (L, -1); - lua_pop (L, 1); - return v; - } - - return luaL_error (L, "function should return a real number"); -} - -void -fdfmultimin_check_gradient (lua_State *L, struct grad_check *gc) -{ - const double abs_err = 1e-5, rel_err = 1e-5; - gsl_matrix *x = matrix_check (L, gc->x_index); - gsl_matrix *gnum, *gfunc; - gsl_function F[1]; - size_t k; - - gnum = matrix_push_raw (L, gc->n, 1); - - F->function = & f_check_hook; - F->params = gc; - - for (k = 0; k < gc->n; k++) - { - double result, abserr; - double x_start = gsl_matrix_get (x, k, 0); - gc->var_vector_index = k; - gsl_deriv_central (F, x_start, gc->step / 50.0, &result, &abserr); - gsl_matrix_set (x, k, 0, x_start); - gsl_matrix_set (gnum, k, 0, result); - } - - gfunc = matrix_push (L, gc->n, 1); - - lua_pushvalue (L, gc->f_index); - lua_pushvalue (L, gc->x_index); - lua_pushvalue (L, -3); // push gfunc - - lua_call (L, 2, 0); - - for (k = 0; k < gc->n; k++) - { - double gn, gf; - gn = gsl_matrix_get (gnum, k, 0); - gf = gsl_matrix_get (gfunc, k, 0); - if (fabs(gn - gf) > abs_err && fabs(gn - gf) >= fabs(gn) * rel_err) - luaL_error (L, "component #%d of gradient is wrong, " - "should be %f but the function gives %f", k+1, gn, gf); - } - - lua_pop (L, 2); -} - -int -gradient_auto_check (lua_State *L) -{ - gsl_matrix *x = matrix_check (L, 2); - double step = luaL_checknumber (L, 3); - gsl_rng *r = push_rng (L, gsl_rng_default); - size_t k, count, nb_tests = 3, n = x->size1; - gsl_matrix *xrnd = matrix_push_raw (L, n, 1); - struct grad_check gc[1]; - - gc->L = L; - gc->x_index = lua_gettop (L); - gc->f_index = 1; - gc->step = step; - gc->n = n; - - for (count = 0; count < nb_tests; count++) - { - for (k = 0; k < n; k++) - { - double u = gsl_matrix_get (x, k, 0) + step * gsl_rng_uniform (r); - gsl_matrix_set (xrnd, k, 0, u); - } - - fdfmultimin_check_gradient (L, gc); - } - - return 0; -} diff --git a/gs-types.c b/gs-types.c index 548a8331..26dcf695 100644 --- a/gs-types.c +++ b/gs-types.c @@ -10,23 +10,7 @@ static int gs_type_string (lua_State *L); -#define GS_MATRIX_NAME_DEF "GSL.matrix" -#define GS_CMATRIX_NAME_DEF "GSL.cmatrix" -#define GS_COMPLEX_NAME_DEF "GSL.cmpl" #define GS_RNG_NAME_DEF "GSL.rng" -#define GS_NLINFIT_NAME_DEF "GSL.solver" -#define GS_CNLINFIT_NAME_DEF "GSL.csolver" -#define GS_ODESOLV_NAME_DEF "GSL.ode" -#define GS_CODESOLV_NAME_DEF "GSL.code" -#define GS_HALFCMPL_R2_NAME_DEF "GSL.ffthcr2" -#define GS_HALFCMPL_MR_NAME_DEF "GSL.ffthcmr" -#define GS_FDFMULTIMIN_NAME_DEF "GSL.fdfmmin" -#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 @@ -48,23 +32,7 @@ static int gs_type_string (lua_State *L); #define MY_EXPAND_DER(NM,DESCR,BASE) {MYCAT2(GS,NM), MYCAT3(GS,NM,NAME_DEF), DESCR, MYCAT2(GS,BASE)} const struct gs_type gs_type_table[] = { - MY_EXPAND(MATRIX, "real matrix"), - MY_EXPAND(CMATRIX, "complex matrix"), - MY_EXPAND(COMPLEX, "complex number"), MY_EXPAND(RNG, "random number generator"), - MY_EXPAND(NLINFIT, "real values non-linear solver"), - MY_EXPAND(CNLINFIT, "complex values non-linear solver"), - MY_EXPAND(ODESOLV, "real values ODE solver"), - MY_EXPAND(CODESOLV, "complex values ODE solver"), - MY_EXPAND(HALFCMPL_R2, "half complex array (radix2)"), - MY_EXPAND(HALFCMPL_MR, "half complex array (mixed radix)"), - MY_EXPAND(FDFMULTIMIN, "fdf multimin solver"), - 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 40e74071..1c104082 100644 --- a/gs-types.h +++ b/gs-types.h @@ -10,24 +10,7 @@ __BEGIN_DECLS enum gs_type_e { GS_NO_TYPE = -1, - GS_MATRIX = 0, /* needs to start from zero because it is used as a - table index */ - GS_CMATRIX, - GS_COMPLEX, - GS_RNG, - GS_NLINFIT, - GS_CNLINFIT, - GS_ODESOLV, - GS_CODESOLV, - GS_HALFCMPL_R2, - GS_HALFCMPL_MR, - GS_FDFMULTIMIN, - GS_FMULTIMIN, - GS_BSPLINE, - GS_INTERP, - GS_LU_DECOMP, - GS_CLU_DECOMP, - GS_QR_DECOMP, + GS_RNG = 0, #ifdef AGG_PLOT_ENABLED GS_DRAW_PLOT, GS_DRAW_SCALABLE, /* derived types are declared only after their base class */ diff --git a/gsl-shell.c b/gsl-shell.c deleted file mode 100644 index f3c31ba2..00000000 --- a/gsl-shell.c +++ /dev/null @@ -1,552 +0,0 @@ -/* - * GSL shell interactive interface to GSL library - * Based on LUA stand-alone interpreter - * - * 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. - -** Lua stand-alone interpreter -** Lua - An Extensible Extension Language -** Lua.org, PUC-Rio, Brazil (http://www.lua.org) - -****************************************************************************** -* Copyright (C) 1994-2008 Lua.org, PUC-Rio. All rights reserved. -* -* Permission is hereby granted, free of charge, to any person obtaining -* a copy of this software and associated documentation files (the -* "Software"), to deal in the Software without restriction, including -* without limitation the rights to use, copy, modify, merge, publish, -* distribute, sublicense, and/or sell copies of the Software, and to -* permit persons to whom the Software is furnished to do so, subject to -* the following conditions: -* -* The above copyright notice and this permission notice shall be -* included in all copies or substantial portions of the Software. -* -* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, -* EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF -* MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. -* IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY -* CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, -* TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE -* SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. -******************************************************************************/ - - -#include <signal.h> -#include <stdio.h> -#include <stdlib.h> -#include <string.h> -#include <pthread.h> -// #include <sys/resource.h> - -#define lua_c - -#include "lua.h" - -#include <lauxlib.h> -#include <lualib.h> -#include "llimits.h" -#include "gsl-shell.h" -#include "lua-gsl.h" -#include "lua-utils.h" - -#if defined(GSL_SHELL_DEBUG) && defined(AGG_PLOT_ENABLED) -#include "window.h" -#include "window_registry.h" -#include "debug-support.h" -#endif - -#define report error_report - -static lua_State *globalL = NULL; - -static const char *progname = LUA_PROGNAME; - -extern int luaopen_gsl (lua_State *L); - -static const luaL_Reg gshlibs[] = { - {"", luaopen_base}, - {LUA_LOADLIBNAME, luaopen_package}, - {LUA_TABLIBNAME, luaopen_table}, - {LUA_IOLIBNAME, luaopen_io}, - {LUA_OSLIBNAME, luaopen_os}, - {LUA_STRLIBNAME, luaopen_string}, - {LUA_DBLIBNAME, luaopen_debug}, -#ifdef LUA_STRICT - {LUA_MATHLIBNAME, luaopen_math}, - {MLUA_GSLLIBNAME, luaopen_gsl}, -#endif - {NULL, NULL} -}; - -pthread_mutex_t gsl_shell_mutex[1]; - -static void lstop (lua_State *L, lua_Debug *ar) { - (void)ar; /* unused arg. */ - lua_sethook(L, NULL, 0, 0); - luaL_error(L, "interrupted!"); -} - - -static void laction (int i) { - signal(i, SIG_DFL); /* if another SIGINT happens before lstop, - terminate process (default action) */ - lua_sethook(globalL, lstop, LUA_MASKCALL | LUA_MASKRET | LUA_MASKCOUNT, 1); -} - - -static void print_usage (void) { - fprintf(stderr, - "usage: %s [options] [script [args]].\n" - "Available options are:\n" - " -e stat execute string " LUA_QL("stat") "\n" - " -l name require library " LUA_QL("name") "\n" - " -i enter interactive mode after executing " LUA_QL("script") "\n" - " -v show version information\n" - " -- stop handling options\n" - " - execute stdin and stop handling options\n" - , - progname); - fflush(stderr); -} - -void -gsl_shell_openlibs (lua_State *L) -{ - const luaL_Reg *lib = gshlibs; - for (; lib->func; lib++) - { - lua_pushcfunction(L, lib->func); - lua_pushstring(L, lib->name); - lua_call(L, 1, 0); - } - -#ifndef LUA_STRICT - lua_pushvalue (L, LUA_GLOBALSINDEX); /* open math in global scope */ - lua_setglobal (L, LUA_MATHLIBNAME); - luaopen_math (L); - lua_pop (L, 1); - lua_pushnil (L); /* remove math */ - lua_setglobal (L, LUA_MATHLIBNAME); - - lua_pushvalue (L, LUA_GLOBALSINDEX); /* open gsl in global scope */ - lua_setglobal (L, MLUA_GSLLIBNAME); - luaopen_gsl (L); - lua_pop (L, 1); - lua_pushnil (L); /* remove gsl */ - lua_setglobal (L, MLUA_GSLLIBNAME); -#endif -} - -static void l_message (const char *pname, const char *msg) { - if (pname) fprintf(stderr, "%s: ", pname); - fprintf(stderr, "%s\n", msg); - fflush(stderr); -} - - -static int error_report (lua_State *L, int status) { - if (status && !lua_isnil(L, -1)) { - const char *msg = lua_tostring(L, -1); - if (msg == NULL) msg = "(error object is not a string)"; - l_message(progname, msg); - lua_pop(L, 1); - } - return status; -} - - -static int traceback (lua_State *L) { - if (!lua_isstring(L, 1)) /* 'message' not a string? */ - return 1; /* keep it intact */ - lua_getfield(L, LUA_GLOBALSINDEX, "debug"); - if (!lua_istable(L, -1)) { - lua_pop(L, 1); - return 1; - } - lua_getfield(L, -1, "traceback"); - if (!lua_isfunction(L, -1)) { - lua_pop(L, 2); - return 1; - } - lua_pushvalue(L, 1); /* pass error message */ - lua_pushinteger(L, 2); /* skip this function and traceback */ - lua_call(L, 2, 1); /* call debug.traceback */ - return 1; -} - - -static int docall (lua_State *L, int narg, int clear) { - int status; - int base = lua_gettop(L) - narg; /* function index */ - lua_pushcfunction(L, traceback); /* push traceback function */ - lua_insert(L, base); /* put it under chunk and args */ - signal(SIGINT, laction); - status = lua_pcall(L, narg, (clear ? 0 : LUA_MULTRET), base); - signal(SIGINT, SIG_DFL); - lua_remove(L, base); /* remove traceback function */ - /* force a complete garbage collection in case of errors */ - if (status != 0) lua_gc(L, LUA_GCCOLLECT, 0); - return status; -} - - -static void print_version (void) { - l_message(NULL, "GSL Shell, Copyright (C) 2009, 2010 Francesco Abbate"); - l_message(NULL, "GNU Scientific Library, Copyright (C) The GSL Team"); - l_message(NULL, LUA_RELEASE ", " LUA_COPYRIGHT); -} - - -static int getargs (lua_State *L, char **argv, int n) { - int narg; - int i; - int argc = 0; - while (argv[argc]) argc++; /* count total number of arguments */ - narg = argc - (n + 1); /* number of arguments to the script */ - luaL_checkstack(L, narg + 3, "too many arguments to script"); - for (i=n+1; i < argc; i++) - lua_pushstring(L, argv[i]); - lua_createtable(L, narg, n + 1); - for (i=0; i < argc; i++) { - lua_pushstring(L, argv[i]); - lua_rawseti(L, -2, i - n); - } - return narg; -} - - -static int dofile (lua_State *L, const char *name) { - int status = luaL_loadfile(L, name) || docall(L, 0, 1); - return report(L, status); -} - - -static int dostring (lua_State *L, const char *s, const char *name) { - int status = luaL_loadbuffer(L, s, strlen(s), name) || docall(L, 0, 1); - return report(L, status); -} - - -static int dolibrary (lua_State *L, const char *name) { - lua_getglobal(L, "require"); - lua_pushstring(L, name); - return report(L, docall(L, 1, 1)); -} - - -static const char *get_prompt (lua_State *L, int firstline) { - const char *p; - lua_getfield(L, LUA_GLOBALSINDEX, firstline ? "_PROMPT" : "_PROMPT2"); - p = lua_tostring(L, -1); - if (p == NULL) p = (firstline ? LUA_PROMPT : LUA_PROMPT2); - lua_pop(L, 1); /* remove global */ - return p; -} - - -static int incomplete (lua_State *L, int status) { - if (status == LUA_ERRSYNTAX) { - size_t lmsg; - const char *msg = lua_tolstring(L, -1, &lmsg); - const char *tp = msg + lmsg - (sizeof(LUA_QL("<eof>")) - 1); - if (strstr(msg, LUA_QL("<eof>")) == tp) { - lua_pop(L, 1); - return 1; - } - } - return 0; /* else... */ -} - - -static int pushline (lua_State *L, int firstline) { - char buffer[LUA_MAXINPUT]; - char *b = buffer; - size_t l; - const char *prmt = get_prompt(L, firstline); - int ok; - - GSL_SHELL_UNLOCK(); - ok = lua_readline(L, b, prmt); - GSL_SHELL_LOCK(); - - if (ok == 0) - { - return 0; /* no input */ - } - - l = strlen(b); - if (l > 0 && b[l-1] == '\n') /* line ends with newline? */ - b[l-1] = '0円'; /* remove it */ - lua_pushstring(L, b); - lua_freeline(L, b); - return 1; -} - - -static int loadline (lua_State *L) { - int status; - lua_settop(L, 0); - if (!pushline(L, 1)) - return -1; /* no input */ - - if (strcmp (lua_tostring(L, 1), "exit") == 0) - return -1; - - lua_pushfstring(L, "return %s", lua_tostring(L, 1)); - status = luaL_loadbuffer(L, lua_tostring(L, 2), lua_strlen(L, 2), "=stdin"); - if (status == 0) - { - lua_saveline(L, 1); - lua_replace (L, 1); - lua_settop (L, 1); - return 0; - } - - lua_settop (L, 1); - for (;;) { /* repeat until gets a complete line */ - status = luaL_loadbuffer(L, lua_tostring(L, 1), lua_strlen(L, 1), "=stdin"); - if (!incomplete(L, status)) break; /* cannot try to add lines? */ - if (!pushline(L, 0)) /* no more input? */ - return -1; - lua_pushliteral(L, "\n"); /* add a new line... */ - lua_insert(L, -2); /* ...between the two lines */ - lua_concat(L, 3); /* join them */ - } - lua_saveline(L, 1); - lua_remove(L, 1); /* remove line */ - return status; -} - -static void dotty (lua_State *L) { - const char *oldprogname = progname; - progname = NULL; - - GSL_SHELL_LOCK(); - - for (;;) - { - int status = loadline(L); - - if (status == -1) - break; - if (status == 0) - status = docall(L, 0, 0); - report(L, status); - if (status == 0 && lua_gettop(L) > 0) - { /* any result to print? */ - lua_getglobal(L, "print"); - lua_insert(L, 1); - if (lua_pcall(L, lua_gettop(L)-1, 0, 0) != 0) - { - const char *fstr = "error calling " LUA_QL("print") " (%s)"; - const char *emsg = lua_pushfstring(L, fstr, lua_tostring(L, -1)); - l_message(progname, emsg); - } - } - } - -#if defined(GSL_SHELL_DEBUG) && defined(AGG_PLOT_ENABLED) - window_index_apply_all (L, window_close); - - do - { - GSL_SHELL_UNLOCK(); - msleep(50); - GSL_SHELL_LOCK(); - } - while (window_index_count (L) > 0); -#endif - - lua_settop(L, 0); /* clear stack */ - - GSL_SHELL_UNLOCK(); - - fputs("\n", stdout); - fflush(stdout); - progname = oldprogname; -} - - -static int handle_script (lua_State *L, char **argv, int n) { - int status; - const char *fname; - int narg = getargs(L, argv, n); /* collect arguments */ - lua_setglobal(L, "arg"); - fname = argv[n]; - if (strcmp(fname, "-") == 0 && strcmp(argv[n-1], "--") != 0) - fname = NULL; /* stdin */ - status = luaL_loadfile(L, fname); - lua_insert(L, -(narg+1)); - if (status == 0) - status = docall(L, narg, 0); - else - lua_pop(L, narg); - return report(L, status); -} - - -/* check that argument has no extra characters at the end */ -#define notail(x) {if ((x)[2] != '0円') return -1;} - - -static int collectargs (char **argv, int *pi, int *pv, int *pe) { - int i; - for (i = 1; argv[i] != NULL; i++) { - if (argv[i][0] != '-') /* not an option? */ - return i; - switch (argv[i][1]) { /* option */ - case '-': - notail(argv[i]); - return (argv[i+1] != NULL ? i+1 : 0); - case '0円': - return i; - case 'i': - notail(argv[i]); - *pi = 1; /* go through */ - case 'v': - notail(argv[i]); - *pv = 1; - break; - case 'e': - *pe = 1; /* go through */ - case 'l': - if (argv[i][2] == '0円') { - i++; - if (argv[i] == NULL) return -1; - } - break; - default: return -1; /* invalid option */ - } - } - return 0; -} - - -static int runargs (lua_State *L, char **argv, int n) { - int i; - for (i = 1; i < n; i++) { - if (argv[i] == NULL) continue; - lua_assert(argv[i][0] == '-'); - switch (argv[i][1]) { /* option */ - case 'e': { - const char *chunk = argv[i] + 2; - if (*chunk == '0円') chunk = argv[++i]; - lua_assert(chunk != NULL); - if (dostring(L, chunk, "=(command line)") != 0) - return 1; - break; - } - case 'l': { - const char *filename = argv[i] + 2; - if (*filename == '0円') filename = argv[++i]; - lua_assert(filename != NULL); - if (dolibrary(L, filename)) - return 1; /* stop if file fails */ - break; - } - default: break; - } - } - return 0; -} - - -static int handle_luainit (lua_State *L) { - const char *init = getenv(LUA_INIT); - if (init == NULL) return 0; /* status OK */ - else if (init[0] == '@') - return dofile(L, init+1); - else - return dostring(L, init, "=" LUA_INIT); -} - - -struct Smain { - int argc; - char **argv; - int status; -}; - - -static int pmain (lua_State *L) { - struct Smain *s = (struct Smain *)lua_touserdata(L, 1); - char **argv = s->argv; - int script; - int has_i = 0, has_v = 0, has_e = 0; - - globalL = L; - if (argv[0] && argv[0][0]) progname = argv[0]; - lua_gc(L, LUA_GCSTOP, 0); /* stop collector during initialization */ - gsl_shell_openlibs (L); /* open libraries */ - lua_gc(L, LUA_GCRESTART, 0); - dolibrary (L, "gslext"); - s->status = handle_luainit(L); - if (s->status != 0) return 0; - script = collectargs(argv, &has_i, &has_v, &has_e); - if (script < 0) { /* invalid args? */ - print_usage(); - s->status = 1; - return 0; - } - if (has_v) print_version(); - s->status = runargs(L, argv, (script > 0) ? script : s->argc); - if (s->status != 0) return 0; - if (script) - s->status = handle_script(L, argv, script); - if (s->status != 0) return 0; - if (has_i) - dotty(L); - else if (script == 0 && !has_e && !has_v) { - if (lua_stdin_is_tty()) { - print_version(); - dotty(L); - } - else dofile(L, NULL); /* executes stdin as a file */ - } - return 0; -} - -int main (int argc, char **argv) { - int status; - struct Smain s; - /* - struct rlimit lmt[1]; - - lmt->rlim_cur = 32 * 1024 * 1024; - lmt->rlim_max = 32 * 1024 * 1024; - - setrlimit (RLIMIT_AS, lmt); - */ - - pthread_mutex_init (gsl_shell_mutex, NULL); - - lua_State *L = lua_open(); /* create state */ - if (L == NULL) { - l_message(argv[0], "cannot create state: not enough memory"); - return EXIT_FAILURE; - } - s.argc = argc; - s.argv = argv; - status = lua_cpcall(L, &pmain, &s); - report(L, status); - lua_close(L); - - pthread_mutex_destroy (gsl_shell_mutex); - - return (status || s.status) ? EXIT_FAILURE : EXIT_SUCCESS; -} diff --git a/integ.c b/integ.c deleted file mode 100644 index 1013c66a..00000000 --- a/integ.c +++ /dev/null @@ -1,450 +0,0 @@ - -/* integ.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 <string.h> -#include <gsl/gsl_integration.h> -#include <gsl/gsl_vector.h> - -#include "lua-utils.h" -#include "integ.h" - -#define INTEG_CACHE_MT_NAME "GSL.integ_cache" - -struct integ_ws_cache { - gsl_integration_workspace *w; - size_t limit; - - gsl_integration_workspace *cycle_ws; - size_t cycle_ws_size; - - gsl_integration_qaws_table *qaws_tb; - - gsl_integration_qawo_table *wf; - size_t blevel; - - gsl_vector *pts; -}; - -struct integ_params_tr { - lua_State *L; - double x_coeff; - double y_coeff; -}; - -static void integ_pushcache (lua_State *L); -static struct integ_ws_cache * integ_retrieve_cache (lua_State *L); - -static int integ_raw(lua_State *L); - -static const struct luaL_Reg integ_functions[] = { - {"integ_raw", integ_raw}, - {NULL, NULL} -}; - -static double -integ_hook_f (double x, void *params) -{ - lua_State *L = params; - double y; - lua_pushvalue (L, 1); - lua_pushnumber (L, x); - lua_call (L, 1, 1); - y = lua_tonumber (L, -1); - lua_pop (L, 1); - return y; -} - -static double -integ_hook_f_tr (double x, void * _params) -{ - struct integ_params_tr *params = _params; - lua_State *L = params->L; - double y; - - lua_pushvalue (L, 1); - lua_pushnumber (L, x * params->x_coeff); - lua_call (L, 1, 1); - y = params->y_coeff * lua_tonumber (L, -1); - lua_pop (L, 1); - return y; -} - -static int -integ_get_gauss_rule (int n) -{ - switch (n) - { - case 15: return GSL_INTEG_GAUSS15; - case 21: return GSL_INTEG_GAUSS21; - case 31: return GSL_INTEG_GAUSS31; - case 41: return GSL_INTEG_GAUSS41; - case 51: return GSL_INTEG_GAUSS15; - case 61: return GSL_INTEG_GAUSS15; - default: - return -1; - } -} - -static void -integ_cache_check_workspace (struct integ_ws_cache *cache, int limit) -{ - const int limit_min = 512; - - if (limit < limit_min) - limit = limit_min; - - if (cache->limit < limit) - { - if (cache->w) - gsl_integration_workspace_free (cache->w); - - cache->limit = limit; - cache->w = gsl_integration_workspace_alloc (cache->limit); - } -} - -static void -integ_cache_check_cycle_ws (struct integ_ws_cache *cache, int size) -{ - if (cache->cycle_ws_size < size) - { - if (cache->cycle_ws) - gsl_integration_workspace_free (cache->cycle_ws); - - cache->cycle_ws_size = size; - cache->cycle_ws = gsl_integration_workspace_alloc (size); - } -} - -static void -integ_cache_check_qawo_workspace (struct integ_ws_cache *cache, int blevel) -{ - if (! cache->wf) - { - cache->blevel = (blevel > 8 ? blevel : 8); - cache->wf = gsl_integration_qawo_table_alloc (1.0, 1.0, GSL_INTEG_SINE, - cache->blevel); - } - else if (cache->blevel < blevel) - { - gsl_integration_qawo_table_free (cache->wf); - cache->blevel = blevel; - cache->wf = gsl_integration_qawo_table_alloc (1.0, 1.0, GSL_INTEG_SINE, - cache->blevel); - } -} - - -static gsl_integration_qaws_table * -integ_cache_set_qaws_table (struct integ_ws_cache *cache, - double alpha, double beta, int mu, int nu) -{ - if (cache->qaws_tb != NULL) - gsl_integration_qaws_table_free (cache->qaws_tb); - - cache->qaws_tb = gsl_integration_qaws_table_alloc (alpha, beta, mu, nu); - - return cache->qaws_tb; -} - -static void -integ_cache_free_qaws_table (struct integ_ws_cache *cache) -{ - if (cache->qaws_tb != NULL) - { - gsl_integration_qaws_table_free (cache->qaws_tb); - cache->qaws_tb = NULL; - } -} - - -static void -integ_cache_check_pts (struct integ_ws_cache *cache, int npts) -{ - if (cache->pts == NULL) - { - size_t my_npts = (npts > 8 ? npts : 8); - cache->pts = gsl_vector_alloc (my_npts); - } - else if (npts > cache->pts->size) - { - gsl_vector_free (cache->pts); - cache->pts = gsl_vector_alloc (npts); - } -} - -struct integ_ws_cache * -integ_retrieve_cache (lua_State *L) -{ - struct integ_ws_cache *cache; - lua_getfield(L, LUA_ENVIRONINDEX, "cache"); - cache = lua_touserdata (L, -1); - lua_pop (L, 1); - return cache; -} - -int -integ_cache_free (lua_State *L) -{ - struct integ_ws_cache *cache = luaL_checkudata (L, 1, INTEG_CACHE_MT_NAME); - if (cache->w) - gsl_integration_workspace_free (cache->w); - if (cache->cycle_ws) - gsl_integration_workspace_free (cache->cycle_ws); - if (cache->wf) - gsl_integration_qawo_table_free (cache->wf); - if (cache->pts) - gsl_vector_free (cache->pts); - if (cache->qaws_tb) - gsl_integration_qaws_table_free (cache->qaws_tb); - return 0; -} - -int -integ_raw(lua_State *L) -{ - struct integ_ws_cache *cache = NULL; - gsl_function f[1]; - const char *inttype = lua_tostring (L, 2); - double epsabs, epsrel; - double a, b; - double result, error; - int rule_key; - int status; - int is_adaptive; - - if (inttype == NULL) - luaL_error (L, "invalid integral type"); - - f->function = & integ_hook_f; - f->params = L; - - epsabs = mlua_named_number (L, 3, "eps_abs"); - epsrel = mlua_named_number (L, 3, "eps_rel"); - - is_adaptive = (inttype[0] == 'a'); - - if (is_adaptive) - { - int limit = (int) mlua_named_number (L, 3, "limit"); - cache = integ_retrieve_cache (L); - integ_cache_check_workspace (cache, limit); - } - - a = mlua_named_optnumber (L, 3, "a", 0); - b = mlua_named_optnumber (L, 3, "b", 1); - - if (strcmp (inttype, "ng") == 0) - { - size_t neval; - status = gsl_integration_qng (f, a, b, epsabs, epsrel, - &result, &error, &neval); - } - else if (strcmp (inttype, "ag") == 0) - { - int rule = mlua_named_optnumber (L, 3, "rule", 21); - - rule_key = integ_get_gauss_rule (rule); - if (rule_key < 0) - luaL_error (L, "invalid integration rule"); - - status = gsl_integration_qag (f, a, b, epsabs, epsrel, cache->limit, - rule_key, - cache->w, &result, &error); - } - else if (strcmp (inttype, "ags") == 0) - { - status = gsl_integration_qags (f, a, b, epsabs, epsrel, cache->limit, - cache->w, &result, &error); - } - else if (strcmp (inttype, "agp") == 0) - { - int npts, j; - - lua_getfield (L, 3, "points"); - npts = lua_objlen (L, -1); - - integ_cache_check_pts (cache, npts); - for (j = 0; j < npts; j++) - { - lua_rawgeti (L, -1, j+1); - gsl_vector_set (cache->pts, j, lua_tonumber (L, -1)); - lua_pop (L, 1); - } - lua_pop (L, 1); - - status = gsl_integration_qagp (f, cache->pts->data, npts, epsabs, epsrel, - cache->limit, - cache->w, &result, &error); - } - else if (strcmp (inttype, "awo") == 0) - { - double omega = mlua_named_number (L, 4, "omega"); - const char *ftype = mlua_named_string (L, 4, "type"); - enum gsl_integration_qawo_enum sine; - const size_t n_bisect_max = 8; - - sine = (strcmp (ftype, "sin") == 0 ? GSL_INTEG_SINE : GSL_INTEG_COSINE); - - integ_cache_check_qawo_workspace (cache, n_bisect_max); - gsl_integration_qawo_table_set (cache->wf, omega, b-a, sine); - status = gsl_integration_qawo (f, a, epsabs, epsrel, cache->limit, - cache->w, cache->wf, &result, &error); - } - else if (strcmp (inttype, "awfu") == 0) - { - double omega = mlua_named_number (L, 4, "omega"); - const char *ftype = mlua_named_string (L, 4, "type"); - enum gsl_integration_qawo_enum sine; - const size_t n_bisect_max = 8; - const size_t cycle_ws_limit = 512; - - sine = (strcmp (ftype, "sin") == 0 ? GSL_INTEG_SINE : GSL_INTEG_COSINE); - - integ_cache_check_qawo_workspace (cache, n_bisect_max); - integ_cache_check_cycle_ws (cache, cycle_ws_limit); - - gsl_integration_qawo_table_set (cache->wf, omega, 1.0, sine); - - status = gsl_integration_qawf (f, a, epsabs, cache->limit, - cache->w, cache->cycle_ws, cache->wf, - &result, &error); - } - else if (strcmp (inttype, "awfl") == 0) - { - double omega = mlua_named_number (L, 4, "omega"); - const char *ftype = mlua_named_string (L, 4, "type"); - enum gsl_integration_qawo_enum sine; - const size_t n_bisect_max = 8; - const size_t cycle_ws_limit = 512; - struct integ_params_tr params[1] = {{L, -1.0, 1.0}}; - - if (strcmp (ftype, "sin") == 0) - { - sine = GSL_INTEG_SINE; - params[0].y_coeff = -1.0; - } - else - { - sine = GSL_INTEG_COSINE; - params[0].y_coeff = 1.0; - } - - f->function = & integ_hook_f_tr; - f->params = params; - - integ_cache_check_qawo_workspace (cache, n_bisect_max); - integ_cache_check_cycle_ws (cache, cycle_ws_limit); - - gsl_integration_qawo_table_set (cache->wf, omega, 1.0, sine); - - status = gsl_integration_qawf (f, -b, epsabs, cache->limit, - cache->w, cache->cycle_ws, cache->wf, - &result, &error); - } - else if (strcmp (inttype, "awc") == 0) - { - double c = mlua_named_number (L, 4, "singularity"); - status = gsl_integration_qawc (f, a, b, c, epsabs, epsrel, cache->limit, - cache->w, &result, &error); - } - else if (strcmp (inttype, "aws") == 0) - { - gsl_integration_qaws_table *t; - double alpha = mlua_named_optnumber (L, 4, "alpha", 0); - double beta = mlua_named_optnumber (L, 4, "beta", 0); - int mu = mlua_named_optnumber (L, 4, "mu", 0); - int nu = mlua_named_optnumber (L, 4, "nu", 0); - - t = integ_cache_set_qaws_table (cache, alpha, beta, mu, nu); - if (t == NULL) - luaL_error (L, "invalid algebraic-logarithmic coefficients"); - - status = gsl_integration_qaws (f, a, b, t, epsabs, epsrel, cache->limit, - cache->w, &result, &error); - - integ_cache_free_qaws_table (cache); - } - else if (strcmp (inttype, "agi") == 0) - { - status = gsl_integration_qagi (f, epsabs, epsrel, cache->limit, - cache->w, &result, &error); - } - else if (strcmp (inttype, "agil") == 0) - { - status = gsl_integration_qagil (f, b, epsabs, epsrel, cache->limit, - cache->w, &result, &error); - } - else if (strcmp (inttype, "agiu") == 0) - { - status = gsl_integration_qagiu (f, a, epsabs, epsrel, cache->limit, - cache->w, &result, &error); - } - else - { - luaL_error (L, "GSL shell internal error integ function"); - } - - if (status != GSL_SUCCESS) - luaL_error (L, "GSL error: %s", gsl_strerror (status)); - - lua_pushnumber (L, result); - lua_pushnumber (L, error); - - return 2; -} - -void -integ_pushcache (lua_State *L) -{ - struct integ_ws_cache *cache; - cache = lua_newuserdata (L, sizeof(struct integ_ws_cache)); - - luaL_getmetatable (L, INTEG_CACHE_MT_NAME); - lua_setmetatable (L, -2); - - cache->w = NULL; - cache->limit = 0; - cache->cycle_ws = NULL; - cache->cycle_ws_size = 0; - cache->qaws_tb = NULL; - cache->wf = NULL; - cache->blevel = 0; - cache->pts = NULL; -} - -void -integ_register (lua_State *L) -{ - luaL_newmetatable (L, INTEG_CACHE_MT_NAME); - lua_pushcfunction (L, integ_cache_free); - lua_setfield (L, -2, "__gc"); - lua_pop (L, 1); - - lua_newtable (L); - integ_pushcache (L); - lua_setfield (L, -2, "cache"); - lua_replace (L, LUA_ENVIRONINDEX); - - /* gsl module registration */ - luaL_register (L, NULL, integ_functions); -} diff --git a/integ.h b/integ.h deleted file mode 100644 index c5fbd58b..00000000 --- a/integ.h +++ /dev/null @@ -1,9 +0,0 @@ -#ifndef INTEG_H -#define INTEG_H - -#include <lua.h> -#include "defs.h" - -extern void integ_register (lua_State *L); - -#endif diff --git a/interp.c b/interp.c deleted file mode 100644 index 3510fc5e..00000000 --- a/interp.c +++ /dev/null @@ -1,206 +0,0 @@ -#include <lua.h> -#include <lauxlib.h> -#include <assert.h> -#include <string.h> -#include <gsl/gsl_interp.h> - -#include "interp.h" -#include "gs-types.h" -#include "lua-utils.h" -#include "matrix.h" - -struct interp { - gsl_interp *interp; - gsl_interp_accel *acc; - const double *xsrc, *ysrc; - int n; -}; - -enum fenv_pos { - FENV_X_VECTOR = 1, - FENV_Y_VECTOR = 2, -}; - -static int interp_new (lua_State *L); -static int interp_index (lua_State *L); -static int interp_free (lua_State *L); -static int interp_eval (lua_State *L); -static int interp_deriv (lua_State *L); -static int interp_deriv2 (lua_State *L); -static int interp_integ (lua_State *L); - -static const struct luaL_Reg interp_metatable[] = { - {"__gc", interp_free}, - {"__index", interp_index}, - {NULL, NULL} -}; - -static const struct luaL_Reg interp_methods[] = { - {"eval", interp_eval}, - {"deriv", interp_deriv}, - {"deriv2", interp_deriv2}, - {"integ", interp_integ}, - {NULL, NULL} -}; - -static const struct luaL_Reg interp_functions[] = { - {"interp", interp_new}, - {NULL, NULL} -}; - -static const gsl_interp_type * -interp_algo_lookup (const char *req_name) -{ - if (strcmp ("linear", req_name) == 0) - return gsl_interp_linear; - else if (strcmp ("polynomial", req_name) == 0) - return gsl_interp_polynomial; - else if (strcmp ("cspline", req_name) == 0) - return gsl_interp_cspline; - else if (strcmp ("cspline_periodic", req_name) == 0) - return gsl_interp_cspline_periodic; - else if (strcmp ("akima", req_name) == 0) - return gsl_interp_akima; - else if (strcmp ("akima_periodic", req_name) == 0) - return gsl_interp_akima_periodic; - - return NULL; -} - -static void -interp_set_ref (lua_State *L, int index_x, int index_y) -{ - /* INDEX_SET_ABS_2(L, index_x, index_y); */ - lua_newtable (L); - lua_pushvalue (L, index_x); - lua_rawseti (L, -2, FENV_X_VECTOR); - lua_pushvalue (L, index_y); - lua_rawseti (L, -2, FENV_Y_VECTOR); - lua_setfenv (L, -2); -} - -int -interp_new (lua_State *L) -{ - gsl_matrix *x = matrix_check (L, 1); - gsl_matrix *y = matrix_check (L, 2); - const gsl_interp_type *T; - size_t n = x->size1; - struct interp *obj; - bool need_copy; - - if (x->size2 != 1 || y->size2 != 1) - return luaL_error (L, "both arguments should be column matrix"); - - if (x->size1 != y->size1) - return luaL_error (L, "mismatch in argument's dimensions"); - - if (lua_gettop (L) > 2) - { - const char *interp_name = luaL_checkstring (L, 3); - T = interp_algo_lookup (interp_name); - if (T == NULL) - return luaL_error (L, "unknown algorithm type"); - } - else - { - T = gsl_interp_linear; - } - - obj = gs_new_object (sizeof(struct interp), L, GS_INTERP); - obj->interp = gsl_interp_alloc (T, n); - obj->acc = gsl_interp_accel_alloc (); - obj->n = n; - - need_copy = (x->tda != 1 || y->tda != 1); - - if (need_copy) - { - gsl_matrix *xsrc = matrix_push_raw (L, n, 1); - gsl_matrix *ysrc = matrix_push_raw (L, n, 1); - gsl_matrix_memcpy (xsrc, x); - gsl_matrix_memcpy (ysrc, y); - obj->xsrc = xsrc->data; - obj->ysrc = ysrc->data; - } - else - { - obj->xsrc = x->data; - obj->ysrc = y->data; - } - - gsl_interp_init (obj->interp, obj->xsrc, obj->ysrc, n); - - if (need_copy) - lua_pop (L, 2); - - interp_set_ref (L, 2, 3); - - return 1; -}; - -int -interp_free (lua_State *L) -{ - struct interp *obj = gs_check_userdata (L, 1, GS_INTERP); - gsl_interp_free (obj->interp); - gsl_interp_accel_free (obj->acc); - return 0; -} - -int -interp_eval (lua_State *L) -{ - struct interp *obj = gs_check_userdata (L, 1, GS_INTERP); - double x = gs_check_number (L, 2, true); - double y = gsl_interp_eval (obj->interp, obj->xsrc, obj->ysrc, x, obj->acc); - lua_pushnumber (L, y); - return 1; -}; - -int -interp_deriv (lua_State *L) -{ - struct interp *obj = gs_check_userdata (L, 1, GS_INTERP); - double x = gs_check_number (L, 2, true); - double d = gsl_interp_eval_deriv (obj->interp, obj->xsrc, obj->ysrc, x, obj->acc); - lua_pushnumber (L, d); - return 1; -}; - -int -interp_deriv2 (lua_State *L) -{ - struct interp *obj = gs_check_userdata (L, 1, GS_INTERP); - double x = gs_check_number (L, 2, true); - double d = gsl_interp_eval_deriv2 (obj->interp, obj->xsrc, obj->ysrc, x, obj->acc); - lua_pushnumber (L, d); - return 1; -}; - -int -interp_integ (lua_State *L) -{ - struct interp *obj = gs_check_userdata (L, 1, GS_INTERP); - double a = gs_check_number (L, 2, true); - double b = gs_check_number (L, 3, true); - double v = gsl_interp_eval_integ (obj->interp, obj->xsrc, obj->ysrc, a, b, obj->acc); - lua_pushnumber (L, v); - return 1; -}; - -int -interp_index (lua_State *L) -{ - return mlua_index_methods (L, interp_methods); -} - -void -interp_register (lua_State *L) -{ - luaL_newmetatable (L, GS_METATABLE(GS_INTERP)); - luaL_register (L, NULL, interp_metatable); - lua_pop (L, 1); - - luaL_register (L, NULL, interp_functions); -} diff --git a/interp.h b/interp.h deleted file mode 100644 index 07e0c79f..00000000 --- a/interp.h +++ /dev/null @@ -1,8 +0,0 @@ -#ifndef INTERP_SYSTEM_H -#define INTERP_SYSTEM_H - -#include <lua.h> - -extern void interp_register (lua_State *L); - -#endif diff --git a/lcomplex.c b/lcomplex.c deleted file mode 100644 index 4739bb17..00000000 --- a/lcomplex.c +++ /dev/null @@ -1,289 +0,0 @@ -/* -* lcomplex.c -* C99 complex nummbers for Lua -* Luiz Henrique de Figueiredo <lhf@tecgraf.puc-rio.br> -* 02 Nov 2009 23:15:43 -* This code is hereby placed in the public domain. -*/ - -#include <math.h> - -#include "lcomplex.h" -#include "matrix_arith.h" - -#include "lua.h" -#include "lauxlib.h" - -#include "gs-types.h" - -#define Z(i) Pget(L,i) -#define O(i) luaL_optnumber(L,i,0) - -#define cadd(z,w) ((z)+(w)) -#define csub(z,w) ((z)-(w)) -#define cmul(z,w) ((z)*(w)) -#define cdiv(z,w) ((z)/(w)) -#define cneg(z) (-(z)) -#define cconj conj - -static Complex Pget(lua_State *L, int i) -{ - switch (lua_type(L,i)) - { - case LUA_TNUMBER: - case LUA_TSTRING: - return luaL_checknumber(L,i); - default: - return *((Complex*) gs_check_userdata (L, i, GS_COMPLEX)); - } -} - -int lua_pushcomplex(lua_State *L, Complex z) -{ - Complex *p= gs_new_object (sizeof(Complex), L, GS_COMPLEX); - *p=z; - return 1; -} - -static int Leq(lua_State *L) /** __eq(z,w) */ -{ - lua_pushboolean(L,Z(1)==Z(2)); - return 1; -} - -static void -fmt_number (double x, char *buffer, size_t bufsize) -{ - size_t len = snprintf (buffer, bufsize, "%g", x); - if (len+1 > bufsize) - buffer[bufsize-1] = '0円'; -} - -#define FMTBUF_SIZE 32 - -static const char * -img_part (lua_State *L, double y, double eps, bool with_sign) -{ - double ay = (y >= 0.0 ? y : -y); - const char *ysign = (y >= 0.0 ? (with_sign ? "+" : "") : "-"); - char buf[FMTBUF_SIZE]; - - if (fabs(ay - 1.0) > eps) - { - fmt_number (ay, buf, FMTBUF_SIZE); - lua_pushfstring (L, "%s%si", ysign, buf); - } - else - { - lua_pushfstring (L, "%si", ysign); - } - - return lua_tostring (L, -1); -} - -static int Ltostring(lua_State *L) /** tostring(z) */ -{ - Complex z = Z(1); - double x = creal(z); - double y = cimag(z); - double eps; - char buf[FMTBUF_SIZE]; - - if (lua_isnoneornil (L, 2)) - { - double nn = sqrt(x*x + y*y); - eps = fmax(1.0e-8 * nn, 1.0e-16); - } - else - { - eps = luaL_checknumber (L, 2); - } - - lua_settop (L, 0); - - if (fabs(y) > eps) - { - if (fabs(x) > eps) - { - const char *img = img_part (L, y, eps, true); - fmt_number (x, buf, FMTBUF_SIZE); - lua_pushfstring (L, "%s%s", buf, img); - } - else - { - img_part (L, y, eps, false); - } - } - else - { - if (fabs(x) <= eps) - { - lua_pushliteral (L, "0"); - } - else - { - fmt_number (x, buf, FMTBUF_SIZE); - lua_pushstring (L, buf); - } - } - - return 1; -} - -int lua_iscomplex (lua_State *L, int i) -{ - if (lua_isnumber (L, i)) - return 1; - else - { - void *p = gs_is_userdata (L, i, GS_COMPLEX); - return (p != NULL); - } - return 0; -} - -Complex luaL_checkcomplex (lua_State *L, int i) -{ - return Pget(L, i); -}; - -Complex lua_tocomplex (lua_State *L, int i) -{ - if (lua_isnumber (L, i)) - { - double n = lua_tonumber (L, i); - return n; - } - else if (gs_is_userdata (L, i, GS_COMPLEX)) - { - Complex *p = lua_touserdata (L, i); - return *p; - } - - return 0; -} - -#define A(f,e) static int L##f(lua_State *L) { return lua_pushcomplex(L,e); } -#define B(f) A(f,c##f(Z(1),Z(2))) -#define F(f) A(f,c##f(Z(1))) -#define G(f) static int L##f(lua_State *L) { lua_pushnumber(L,c##f(Z(1))); return 1; } - -#define RFIMP(f,er,ez) static int L##f(lua_State *L) { \ - if (lua_isnumber (L, 1)) \ - lua_pushnumber(L,er(lua_tonumber(L,1))); \ - else { \ - Complex *z = gs_check_userdata (L, 1, GS_COMPLEX); \ - lua_pushcomplex (L, ez(*z)); \ - } \ - return 1; \ -} - -#define RF(f) RFIMP(f,f,c##f) - -int -Lsqrt(lua_State *L) -{ - Complex z; - - switch (lua_type(L, 1)) - { - case LUA_TNUMBER: - case LUA_TSTRING: - { - double x = luaL_checknumber(L, 1); - if (x >= 0) - { - lua_pushnumber (L, sqrt(x)); - return 1; - } - - z = x; - break; - } - default: - z = *((Complex*) gs_check_userdata (L, 1, GS_COMPLEX)); - } - - lua_pushcomplex (L, csqrt(z)); - return 1; -} - -A(new,O(1)+O(2)*I) /** new(x,y) */ -F(neg) /** __unm(z) */ -G(abs) /** abs(z) */ -RF(acos) /** acos(z) */ -RF(acosh) /** acosh(z) */ -G(arg) /** arg(z) */ -RF(asin) /** asin(z) */ -RF(asinh) /** asinh(z) */ -RF(atan) /** atan(z) */ -RF(atanh) /** atanh(z) */ -F(conj) /** conj(z) */ -RF(cos) /** cos(z) */ -RF(cosh) /** cosh(z) */ -RF(exp) /** exp(z) */ -G(imag) /** imag(z) */ -RF(log) /** log(z) */ -B(pow) /** pow(z,w) */ -F(proj) /** proj(z) */ -G(real) /** real(z) */ -RF(sin) /** sin(z) */ -RF(sinh) /** sinh(z) */ -RF(tan) /** tan(z) */ -RF(tanh) /** tanh(z) */ - -static const luaL_Reg lcomplex_methods[] = -{ - { "__add", matrix_op_add }, - { "__div", matrix_op_div }, - { "__eq", Leq }, - { "__mul", matrix_op_mul }, - { "__sub", matrix_op_sub }, - { "__unm", Lneg }, - { "__pow", Lpow }, - { "__tostring", Ltostring}, - { NULL, NULL } -}; - - -static const luaL_Reg lcomplex_functions[] = -{ - { "new", Lnew }, - { "tostring_eps", Ltostring}, - { "abs", Labs }, - { "acos", Lacos }, - { "acosh", Lacosh }, - { "arg", Larg }, - { "asin", Lasin }, - { "asinh", Lasinh }, - { "atan", Latan }, - { "atanh", Latanh }, - { "conj", Lconj }, - { "cos", Lcos }, - { "cosh", Lcosh }, - { "exp", Lexp }, - { "imag", Limag }, - { "log", Llog }, - { "pow", Lpow }, - { "proj", Lproj }, - { "real", Lreal }, - { "sin", Lsin }, - { "sinh", Lsinh }, - { "sqrt", Lsqrt }, - { "tan", Ltan }, - { "tanh", Ltanh }, - { NULL, NULL } -}; - -int luaopen_lcomplex (lua_State *L) -{ - luaL_newmetatable (L, GS_METATABLE(GS_COMPLEX)); - luaL_register (L, NULL, lcomplex_methods); - lua_pop (L, 1); - - luaL_register(L, "complex", lcomplex_functions); - - lua_pushcomplex(L, I); - lua_setfield(L, -2, "I"); - return 1; -} diff --git a/lcomplex.h b/lcomplex.h deleted file mode 100644 index 976f76eb..00000000 --- a/lcomplex.h +++ /dev/null @@ -1,26 +0,0 @@ -#ifndef LUA_COMPLEX_H -#define LUA_COMPLEX_H - -#include "defs.h" - -#ifdef __cplusplus -#error "cannot include C99 complex in C++" -#else -#include <complex.h> -#endif - -__BEGIN_DECLS -#include "lua.h" - -#define Complex double _Complex - -extern int lua_pushcomplex(lua_State *L, Complex z); -extern int lua_iscomplex (lua_State *L, int i); -extern Complex lua_tocomplex (lua_State *L, int i); -extern Complex luaL_checkcomplex (lua_State *L, int i); - -extern int luaopen_lcomplex (lua_State *L); - -__END_DECLS - -#endif diff --git a/linalg.c b/linalg.c deleted file mode 100644 index f364fa3d..00000000 --- a/linalg.c +++ /dev/null @@ -1,69 +0,0 @@ - -/* linalg.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_linalg.h> -#include <gsl/gsl_blas.h> -#include "matrix.h" - -static int linalg_svd (lua_State *L); - -static const struct luaL_Reg linalg_functions[] = { - {"svd", linalg_svd}, - {NULL, NULL} -}; - -int -linalg_svd (lua_State *L) -{ - const gsl_matrix *a = matrix_check (L, 1); - int sm = a->size1, sn = a->size2; - gsl_matrix *u, *v, *s; - gsl_vector *s_vec, *work; - int k; - - u = matrix_push_raw (L, sm, sn); - s = matrix_push_raw (L, sn, sn); - v = matrix_push_raw (L, sn, sn); - - s_vec = gsl_vector_alloc (sn); - work = gsl_vector_alloc (sn); - - gsl_matrix_memcpy (u, a); - gsl_linalg_SV_decomp (u, v, s_vec, work); - - for (k = 0; k < sn; k++) - { - double z = gsl_vector_get (s_vec, k); - gsl_matrix_set (s, k, k, z); - } - - gsl_vector_free (s_vec); - gsl_vector_free (work); - - return 3; -} - -void -linalg_register (lua_State *L) -{ - luaL_register (L, NULL, linalg_functions); -} diff --git a/linalg.h b/linalg.h deleted file mode 100644 index 88da4bca..00000000 --- a/linalg.h +++ /dev/null @@ -1,26 +0,0 @@ - -/* linalg.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 LINALG_H -#define LINALG_H - -extern void linalg_register (lua_State *L); - -#endif diff --git a/lu_decomp.c b/lu_decomp.c deleted file mode 100644 index 7af9990d..00000000 --- a/lu_decomp.c +++ /dev/null @@ -1,51 +0,0 @@ - -#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 deleted file mode 100644 index 96eec6a3..00000000 --- a/lu_decomp.h +++ /dev/null @@ -1,9 +0,0 @@ - -#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 deleted file mode 100644 index 32de3058..00000000 --- a/lu_decomp_complex.c +++ /dev/null @@ -1,17 +0,0 @@ - -#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 deleted file mode 100644 index e93eb40c..00000000 --- a/lu_decomp_imp.h +++ /dev/null @@ -1,15 +0,0 @@ - -#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 deleted file mode 100644 index 0fd1b6f1..00000000 --- a/lu_decomp_real.c +++ /dev/null @@ -1,15 +0,0 @@ - -#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 deleted file mode 100644 index 754a33f6..00000000 --- a/lu_decomp_source.c +++ /dev/null @@ -1,189 +0,0 @@ - -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); -} @@ -24,30 +24,11 @@ #include "lua-gsl.h" #include "gs-types.h" #include "lua-utils.h" -#include "nlinfit.h" -#include "cnlinfit.h" -#include "matrix.h" -#include "cmatrix.h" -#include "lcomplex.h" -#include "matrix_arith.h" -#include "linalg.h" -#include "integ.h" -#include "fft.h" -#include "ode_solver.h" -#include "ode.h" -#include "code.h" #include "random.h" #include "randist.h" #include "pdf.h" #include "cdf.h" #include "sf.h" -#include "multimin.h" -#include "eigen-systems.h" -#include "mlinear.h" -#include "bspline.h" -#include "interp.h" -#include "lu_decomp.h" -#include "qr_decomp.h" #ifdef AGG_PLOT_ENABLED #include "lua-graph.h" @@ -78,39 +59,15 @@ luaopen_gsl (lua_State *L) lua_pop (L, 1); #endif - luaopen_lcomplex (L); - lua_pop (L, 1); - - luaL_register (L, MLUA_MATRIXLIBNAME, matrix_functions); - matrix_register (L); - matrix_complex_register (L); - matrix_arith_register (L); - linalg_register (L); - lu_decomp_register (L); - qr_decomp_register (L); - lua_pop (L, 1); - luaL_register (L, MLUA_GSLLIBNAME, gsl_shell_functions); luaL_register (L, NULL, gs_type_functions); - solver_register (L); - integ_register (L); - ode_register (L); random_register (L); randist_register (L); pdf_register (L); cdf_register (L); sf_register (L); - multimin_register (L); - eigen_register (L); - mlinear_register (L); - bspline_register (L); - interp_register (L); - - fft_register (L); - ode_complex_register (L); - solver_complex_register (L); lua_pop (L, 1); diff --git a/lua-utils.c b/lua-utils.c index cb7b8bea..39a779dc 100644 --- a/lua-utils.c +++ b/lua-utils.c @@ -79,15 +79,6 @@ mlua_get_property (lua_State *L, const struct luaL_Reg *p, bool use_cache) return rval; } -void -mlua_null_cache (lua_State *L, int index) -{ - lua_getfenv (L, index); - lua_pushnil (L); - lua_setfield (L, -2, CACHE_FIELD_NAME); - lua_pop (L, 1); -} - int mlua_index_with_properties (lua_State *L, const struct luaL_Reg *properties, @@ -118,26 +109,6 @@ mlua_index_with_properties (lua_State *L, } int -mlua_index_methods (lua_State *L, 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 (methods, key); - if (reg) - { - lua_pushcfunction (L, reg->func); - return 1; - } - - return 0; -} - -int mlua_newindex_with_properties (lua_State *L, const struct luaL_Reg *properties) { char const * key; @@ -157,22 +128,6 @@ mlua_newindex_with_properties (lua_State *L, const struct luaL_Reg *properties) return luaL_error (L, "invalid property for %s object", full_type_name (L, 1)); } -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) @@ -218,20 +173,3 @@ mlua_named_string (lua_State *L, int index, const char *key) lua_pop (L, 1); return r; } - -void -mlua_fenv_set (lua_State *L, int index, int fenv_index) -{ - lua_getfenv (L, index); - lua_insert (L, -2); - lua_rawseti (L, -2, fenv_index); - lua_pop (L, 1); -} - -void -mlua_fenv_get (lua_State *L, int index, int fenv_index) -{ - lua_getfenv (L, index); - lua_rawgeti (L, -1, fenv_index); - lua_remove (L, -2); -} diff --git a/lua-utils.h b/lua-utils.h index a48e1a7a..16e21a98 100644 --- a/lua-utils.h +++ b/lua-utils.h @@ -36,15 +36,6 @@ mlua_get_property (lua_State *L, const struct luaL_Reg *p, bool use_cache); extern const struct luaL_Reg * 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_methods (lua_State *L, const struct luaL_Reg *methods); - extern int mlua_index_with_properties (lua_State *L, const struct luaL_Reg *properties, const struct luaL_Reg *methods, @@ -67,12 +58,6 @@ extern lua_Number mlua_named_optnumber (lua_State *L, int index, extern lua_Number mlua_named_number (lua_State *L, int index, const char *key); -extern void mlua_fenv_set (lua_State *L, int index, int fenv_index); -extern void mlua_fenv_get (lua_State *L, int index, int fenv_index); -extern void mlua_table_clear (lua_State *L, int index); - -extern void mlua_set_fenv_ref (lua_State *L, int refidx); - __END_DECLS #endif diff --git a/matrix-init.lua b/matrix-init.lua index af50d334..b3d6bf4e 100644 --- a/matrix-init.lua +++ b/matrix-init.lua @@ -1,198 +1,547 @@ - -- matrix.lua - -- - -- 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. - -- - -local cat, insert, fmt = table.concat, table.insert, string.format - -local sqrt, abs, tostring_eps = math.sqrt, math.abs, complex.tostring_eps - -local function matrix_f_set(m, f) - local r, c = matrix.dim(m) - local mset = m.set - for i = 1, r do - for j = 1, c do - local z = f(i, j) - mset(m, i, j, z) +local ffi = require 'ffi' +local cgsl = require 'cgsl' + +local sqrt, abs = math.sqrt, math.abs +local fmt = string.format + +local gsl_matrix = ffi.typeof('gsl_matrix') +local gsl_matrix_complex = ffi.typeof('gsl_matrix_complex') +local gsl_complex = ffi.typeof('complex') + +local gsl_check = require 'gsl-check' + +local function isreal(x) return type(x) == 'number' end + +local function cartesian(x) + if isreal(x) then + return x, 0 + else + return x[0], x[1] + end +end + +local function complex_conj(a) + local x, y = cartesian(z) + return gsl_complex(x, -y) +end + +local function complex_real(z) + local x = cartesian(z) + return x +end + +local function complex_imag(z) + local x, y = cartesian(z) + return y +end + +function matrix_alloc(n1, n2) + local block = cgsl.gsl_block_alloc(n1 * n2) + local m = gsl_matrix(n1, n2, n2, block.data, block, 1) + return m +end + +function matrix_calloc(n1, n2) + local block = cgsl.gsl_block_complex_alloc(n1 * n2) + local m = gsl_matrix_complex(n1, n2, n2, block.data, block, 1) + return m +end + +function matrix_new(n1, n2, f) + local m + if f then + m = matrix_alloc(n1, n2) + for i=0, n1-1 do + for j=0, n2-1 do + m.data[i*n2+j] = f(i, j) + end end + else + local block = cgsl.gsl_block_calloc(n1 * n2) + m = gsl_matrix(n1, n2, n2, block.data, block, 1) end return m end -function matrix.matrix_reduce(m, f, accu) - local r, c = matrix.dim(m) - local mget = m.get - for i = 1, r do - for j = 1, c do - accu = f(accu, mget(m, i, j)) +function matrix_cnew(n1, n2, f) + local m + if f then + m = matrix_calloc(n1, n2) + for i=0, n1-1 do + for j=0, n2-1 do + local z = f(i, j) + m.data[2*i*n2+2*j ] = z[0] + m.data[2*i*n2+2*j+1] = z[1] + end end + else + local block = cgsl.gsl_block_complex_calloc(n1 * n2) + m = gsl_matrix_complex(n1, n2, n2, block.data, block, 1) end - return accu + return m +end + +local function matrix_dim(m) + return m.size1, m.size2 end -local ctor_table = {['number'] = matrix.new, ['complex number'] = matrix.cnew} +local function itostr(im, signed) + local sign = im < 0 and '-' or (signed and '+' or '') + if im == 0 then return '' else + return sign .. (abs(im) == 1 and 'i' or fmt('%gi', abs(im))) + end +end -local function number_type(a, t) - if gsl.type(t) == 'number' then - if a ~= 'complex number' then a = 'number' end - elseif gsl.type(t) == 'complex number' then - a = 'complex number' +local function recttostr(x, y, eps) + if abs(x) < eps then x = 0 end + if abs(y) < eps then y = 0 end + if x ~= 0 then + return (y == 0 and fmt('%g', x) or fmt('%g%s', x, itostr(y, true))) else - error('expecting real or complex number, got :' .. t) + return (y == 0 and '0' or itostr(y)) end - return a - end - -local function matrix_from_table(t) - local tp - for i, ts in ipairs(t) do - if type(ts) ~= 'table' then error 'invalid matrix definition' end - for j, v in ipairs(ts) do - tp = number_type(tp, v) - end +end + +local function concat_pad(t, pad) + local sep = ' ' + local row + for i, s in ipairs(t) do + local x = string.rep(' ', pad - #s) .. s + row = row and (row .. sep .. x) or x end - local ctor = ctor_table[tp] - if not ctor then error 'empty list for matrix definition' end + return row +end - local r, c = #t, #t[1] - return matrix_f_set(ctor(r, c), function(i,j) return t[i][j] end) +local function matrix_tostring_gen(sel) + return function(m) + local n1, n2 = m.size1, m.size2 + local sq = 0 + for i=0, n1-1 do + for j=0, n2-1 do + local x, y = sel(m, i, j) + sq = sq + x*x + y*y + end + end + local eps = sqrt(sq) * 1e-10 + + lsrow = {} + local lmax = 0 + for i=0, n1-1 do + local row = {} + for j=0, n2-1 do + local x, y = sel(m, i, j) + local s = recttostr(x, y, eps) + if #s > lmax then lmax = #s end + row[j+1] = s + end + lsrow[i+1] = row + end + + local ss = {} + for i=0, n1-1 do + ss[i+1] = '[ ' .. concat_pad(lsrow[i+1], lmax) .. ' ]' + end + + return table.concat(ss, '\n') + end end -local function vector_from_table(t) - local tp - for i, v in ipairs(t) do tp = number_type(tp, v) end - local ctor = ctor_table[tp] - if not ctor then error 'empty list for vector definition' end +local function matrix_copy(a) + local n1, n2 = a.size1, a.size2 + local b = matrix_alloc(n1, n2) + cgsl.gsl_matrix_memcpy(b, a) + return b +end - local v = ctor (#t, 1) - for i, x in ipairs(t) do v:set(i,1, x) end - return v +local function matrix_complex_copy(a) + local n1, n2 = a.size1, a.size2 + local b = matrix_calloc(n1, n2) + cgsl.gsl_matrix_complex_memcpy(b, a) + return b end -matrix.vec = vector_from_table -matrix.def = matrix_from_table +local function matrix_col(m, j) + local r = matrix_alloc(m.size1, 1) + local tda = m.tda + for i = 0, m.size1 - 1 do + r.data[i] = m.data[i * tda + j] + end + return r +end -local function padstr(s, w) - return fmt('%s%s', string.rep(' ', w - #s), s) +local function matrix_row(m, i) + local r = matrix_alloc(1, m.size2) + local tda = m.tda + for j = 0, m.size2 - 1 do + r.data[j] = m.data[i * tda + j] + end + return r end -local function matrix_to_string(m) - local eps = m:norm() * 1e-8 - local fwidth = function(w, val) - local ln = # tostring_eps(val, eps) - return (ln > w and ln or w) - end - local width = matrix.matrix_reduce(m, fwidth, 0) - local r, c = matrix.dim(m) - local lines = {} - for i=1,r do - local ln = {} - for j=1,c do - insert(ln, padstr(tostring_eps(m:get(i,j), eps), width)) +local function matrix_vect_def(t) + return matrix_new(#t, 1, |i| t[i+1]) +end + +local function get_typeid(a) + if isreal(a) then return true, true + elseif ffi.istype(gsl_complex, a) then return false, true + elseif ffi.istype(gsl_matrix, a) then return true, false + elseif ffi.istype(gsl_matrix_complex, a) then return false, false end +end + +local function mat_op_gen(n1, n2, opa, a, opb, b, oper) + local c = matrix_alloc(n1, n2) + for i = 0, n1-1 do + for j = 0, n2-1 do + local ar = opa(a,i,j) + local br = opb(b,i,j) + c.data[i*n2+j] = oper(ar, br) end - insert(lines, fmt('[ %s ]', cat(ln, ' '))) end - return cat(lines, '\n') + return c end -local function csqr(z) - local r, i = complex.real(z), complex.imag(z) - return r*r + i*i +local function mat_comp_op_gen(n1, n2, opa, a, opb, b, oper) + local c = matrix_calloc(n1, n2) + for i = 0, n1-1 do + for j = 0, n2-1 do + local ar, ai = opa(a,i,j) + local br, bi = opb(b,i,j) + local zr, zi = oper(ar, br, ai, bi) + c.data[2*i*n2+2*j ] = zr + c.data[2*i*n2+2*j+1] = zi + end + end + return c +end + +local function real_get(x) return x, 0 end +local function complex_get(z) return z[0], z[1] end +local function mat_real_get(m,i,j) return m.data[i*m.tda+j], 0 end + +local function mat_complex_get(m,i,j) + local idx = 2*i*m.tda+2*j + return m.data[idx], m.data[idx+1] +end + +local function selector(r, s) + if s then + return r and real_get or complex_get + else + return r and mat_real_get or mat_complex_get + end end -function matrix.tr(m) - local r, c = matrix.dim(m) - return matrix.new(c, r, function(i,j) return m:get(j,i) end) +local function mat_complex_of_real(m) + local n1, n2 = m.size1, m.size2 + local mc = matrix_calloc(n1, n2) + for i=0, n1-1 do + for j=0, n2-1 do + mc.data[2*i*n2+2*j ] = m.data[i*n2+j] + mc.data[2*i*n2+2*j+1] = 0 + end + end + return mc end -function matrix.hc(m) - local r, c = matrix.dim(m) - return matrix.cnew(c, r, function(i,j) return complex.conj(m:get(j,i)) end) +local function opadd(ar, br, ai, bi) + if ai then return ar+br, ai+bi else return ar+br end end -function matrix.diag(v) - local n = matrix.dim(v) - return matrix.new(n, n, function(i,j) return i == j and v:get(i,1) or 0 end) +local function opsub(ar, br, ai, bi) + if ai then return ar-br, ai-bi else return ar-br end end -function matrix.unit(n) - return matrix.new(n, n, function(i,j) return i == j and 1 or 0 end) +local function opmul(ar, br, ai, bi) + if ai then return ar*br-ai*bi, ar*bi+ai*br else return ar*br end end -local function matrix_norm(m) - local r, c = matrix.dim(m) - local s = 0 - for i=1, r do - for j=1, c do - s = s + csqr(m:get(i,j)) - end +local function opdiv(ar, br, ai, bi) + if ai then + local d = br^2 + bi^2 + return (ar*br + ai*bi)/d, (-ar*bi + ai*br)/d + else + return ar/br end - return sqrt(s) end -local function matrix_column (m, c) - local r = matrix.dim(m) - return m:slice(1, c, r, 1) +local function vector_op(scalar_op, element_wise, no_inverse) + return function(a, b) + local ra, sa = get_typeid(a) + local rb, sb = get_typeid(b) + if not sb and no_inverse then + error 'invalid operation on matrix' + end + if sa and sb then + local ar, ai = cartesian(a) + local br, bi = cartesian(b) + local zr, zi = scalar_op(ar, br, ai, bi) + return gsl_complex(zr, zi) + elseif element_wise or sa or sb then + local sela, selb = selector(ra, sa), selector(rb, sb) + local n1 = (sa and b.size1 or a.size1) + local n2 = (sa and b.size2 or a.size2) + if ra and rb then + return mat_op_gen(n1, n2, sela, a, selb, b, scalar_op) + else + return mat_comp_op_gen(n1, n2, sela, a, selb, b, scalar_op) + end + else + if ra and rb then + local n1, n2 = a.size1, b.size2 + local c = matrix_new(n1, n2) + local NT = cgsl.CblasNoTrans + gsl_check(cgsl.gsl_blas_dgemm(NT, NT, 1, a, b, 1, c)) + return c + else + if ra then a = mat_complex_of_real(a) end + if rb then b = mat_complex_of_real(b) end + local n1, n2 = a.size1, b.size2 + local c = matrix_cnew(n1, n2) + local NT = cgsl.CblasNoTrans + gsl_check(cgsl.gsl_blas_zgemm(NT, NT, 1, a, b, 1, c)) + return c + end + end + end end -local function matrix_row (m, r) - local _, c = matrix.dim(m) - return m:slice(r, 1, 1, c) +complex = { + new = gsl_complex, + conj = complex_conj, + real = complex_real, + imag = complex_imag, +} + +local generic_add = vector_op(opadd, true) +local generic_sub = vector_op(opsub, true) +local generic_mul = vector_op(opmul, false) +local generic_div = vector_op(opdiv, true, true) + +local complex_mt = { + + __add = generic_add, + __sub = generic_sub, + __mul = generic_mul, + __div = generic_div, + + __pow = function(z,n) + if isreal(n) then + return cgsl.gsl_complex_pow_real (z, n) + else + if isreal(z) then z = gsl_complex(z,0) end + return cgsl.gsl_complex_pow (z, n) + end + end, +} + +ffi.metatype(gsl_complex, complex_mt) + +matrix = { + new = matrix_new, + cnew = matrix_cnew, + alloc = matrix_alloc, + calloc = matrix_calloc, + copy = matrix_copy, + dim = matrix_dim, + vec = matrix_vect_def, + col = matrix_col, + row = matrix_row, + get = cgsl.gsl_matrix_get, + set = cgsl.gsl_matrix_set, +} + +local matrix_methods = { + col = matrix_col, + row = matrix_row, + get = cgsl.gsl_matrix_get, + set = cgsl.gsl_matrix_set, +} + +local matrix_mt = { + + __gc = function(m) if m.owner then cgsl.gsl_block_free(m.block) end end, + + __add = generic_add, + __sub = generic_sub, + __mul = generic_mul, + __div = generic_div, + + __index = function(m, k) + if type(k) == 'number' then + if m.size2 == 1 then + return cgsl.gsl_matrix_get(m, k, 0) + else + return matrix_row(m, k) + end + end + return matrix_methods[k] + end, + + + __newindex = function(m, k, v) + if type(k) == 'number' then + if m.size2 == 1 then + cgsl.gsl_matrix_set(m, k, 0, v) + else + local row = cgsl.gsl_matrix_submatrix(m, k, 0, 1, m.size2) + gsl_check(cgsl.gsl_matrix_memcpy(row, v)) + end + else + error 'cannot set a matrix field' + end + end, + + __tostring = matrix_tostring_gen(mat_real_get), +} + +ffi.metatype(gsl_matrix, matrix_mt) + +local matrix_complex_mt = { + + __gc = function(m) if m.owner then cgsl.gsl_block_free(m.block) end end, + + __add = generic_add, + __sub = generic_sub, + __mul = generic_mul, + __div = generic_div, + + __index = function(m, k) + if type(k) == 'number' then + if m.size2 == 1 then + return cgsl.gsl_matrix_complex_get(m, k, 0) + else + return matrix_row(m, k) + end + end + return matrix_complex_methods[k] + end, + + __tostring = matrix_tostring_gen(mat_complex_get), +} + +ffi.metatype(gsl_matrix_complex, matrix_complex_mt) + +local function cwrap(name) + local fc = cgsl['gsl_complex_' .. name] + local fr = math[name] + return function(x) + return isreal(x) and fr(x) or fc(x) + end end -local function matrix_rows(m) - local r, c = matrix.dim(m) - return matrix.sequence(function(i) m:slice(i, 1, 1, c) end, r) +gsl_function_list = { + 'sqrt', 'exp', 'log', 'log10', + 'sin', 'cos', 'sec', 'csc', 'tan', 'cot', + 'arcsin', 'arccos', 'arcsec', 'arccsc', 'arctan', 'arccot', + 'sinh', 'cosh', 'sech', 'csch', 'tanh', 'coth', + 'arcsinh', 'arccosh', 'arcsech', 'arccsch', 'arctanh', 'arccoth' +} + +for _, name in ipairs(gsl_function_list) do + complex[name] = cwrap(name) end -function matrix.null(m) - local r, c = matrix.dim(m) - local mset = m.set - for i=1, r do - for j=1, c do - mset(m, i, j, 0) +local function matrix_def(t) + local r, c = #t, #t[1] + local real = true + for i, row in ipairs(t) do + for j, x in ipairs(row) do + if not isreal(x) then + real = false + break + end end + if not real then break end + end + if real then + local m = matrix_alloc(r, c) + for i= 0, r-1 do + local row = t[i+1] + for j = 0, c-1 do + m.data[i*c+j] = row[j+1] + end + end + return m + else + local m = matrix_calloc(r, c) + for i= 0, r-1 do + local row = t[i+1] + for j = 0, c-1 do + local x, y = cartesian(row[j+1]) + m.data[2*i*c+2*j ] = x + m.data[2*i*c+2*j+1] = y + end + end + return m end end -function matrix.fset(m, f) - matrix_f_set(m, f) +local signum = ffi.new('int[1]') + +function matrix_inv(m) + local n = m.size1 + local lu = matrix_copy(m) + local p = ffi.gc(cgsl.gsl_permutation_alloc(n), cgsl.gsl_permutation_free) + gsl_check(cgsl.gsl_linalg_LU_decomp(lu, p, signum)) + local mi = matrix_alloc(n, n) + gsl_check(cgsl.gsl_linalg_LU_invert(lu, p, mi)) + return mi end -local function add_matrix_method(s, m) - matrix.Matrix[s] = m - matrix.cMatrix[s] = m +function matrix_solve(m, b) + local n = m.size1 + local lu = matrix_copy(m) + local p = ffi.gc(cgsl.gsl_permutation_alloc(n), cgsl.gsl_permutation_free) + gsl_check(cgsl.gsl_linalg_LU_decomp(lu, p, signum)) + local x = matrix_alloc(n, 1) + local xv = cgsl.gsl_matrix_column(x, 0) + local bv = cgsl.gsl_matrix_column(b, 0) + gsl_check(cgsl.gsl_linalg_LU_solve(lu, p, bv, xv)) + return x end -local function add_matrix_meta_method(key, method) - local m, mt - m = matrix.new(1,1) - mt = getmetatable(m) - mt[key] = method +function matrix_complex_inv(m) + local n = m.size1 + local lu = matrix_complex_copy(m) + local p = ffi.gc(cgsl.gsl_permutation_alloc(n), cgsl.gsl_permutation_free) + gsl_check(cgsl.gsl_linalg_complex_LU_decomp(lu, p, signum)) + local mi = matrix_calloc(n, n) + gsl_check(cgsl.gsl_linalg_complex_LU_invert(lu, p, mi)) + return mi +end - m = matrix.cnew(1,1) - mt = getmetatable(m) - mt[key] = method +function matrix_complex_solve(m, b) + local n = m.size1 + local lu = matrix_complex_copy(m) + local p = ffi.gc(cgsl.gsl_permutation_alloc(n), cgsl.gsl_permutation_free) + gsl_check(cgsl.gsl_linalg_complex_LU_decomp(lu, p, signum)) + local x = matrix_calloc(n, 1) + local xv = cgsl.gsl_matrix_complex_column(x, 0) + local bv = cgsl.gsl_matrix_complex_column(b, 0) + gsl_check(cgsl.gsl_linalg_complex_LU_solve(lu, p, bv, xv)) + return x end -add_matrix_meta_method('__tostring', matrix_to_string) +matrix.inv = function(m) + if ffi.istype(gsl_matrix, m) then + return matrix_inv(m) + else + return matrix_complex_inv(m) + end + end + +matrix.solve = function(m, b) + local mr = ffi.istype(gsl_matrix, m) + local br = ffi.istype(gsl_matrix, b) + if mr and br then + return matrix_solve(m, b) + else + if mr then m = mat_complex_of_real(m) end + if br then b = mat_complex_of_real(b) end + return matrix_complex_solve(m, b) + end + end -add_matrix_method('norm', matrix_norm) -add_matrix_method('col', matrix_column) -add_matrix_method('row', matrix_row) -add_matrix_method('rows', matrix_rows) +matrix.def = matrix_def diff --git a/matrix.c b/matrix.c deleted file mode 100644 index 30333e82..00000000 --- a/matrix.c +++ /dev/null @@ -1,62 +0,0 @@ - -/* matrix.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 <string.h> -#include <gsl/gsl_matrix.h> -#include <gsl/gsl_blas.h> -#include <gsl/gsl_permutation.h> -#include <gsl/gsl_linalg.h> - -#include "gs-types.h" -#include "matrix.h" -#include "matrix_arith.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 i, j; - for (i = 0; i < n; i++) - { - for (j = 0; j < n; j++) - { - gsl_complex z = gsl_matrix_complex_get (jcmpl, i, j); - gsl_matrix_set (jreal, 2*i, 2*j, GSL_REAL(z)); - gsl_matrix_set (jreal, 2*i, 2*j+1, - GSL_IMAG(z)); - gsl_matrix_set (jreal, 2*i+1, 2*j, GSL_IMAG(z)); - gsl_matrix_set (jreal, 2*i+1, 2*j+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" - -#include "template_matrix_off.h" -#undef BASE_DOUBLE diff --git a/matrix.h b/matrix.h deleted file mode 100644 index 617c4e1f..00000000 --- a/matrix.h +++ /dev/null @@ -1,39 +0,0 @@ - -/* matrix.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 MATRIX_H -#define MATRIX_H - -#include <lua.h> -#include <gsl/gsl_matrix.h> - -#include "defs.h" - -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_arith.c b/matrix_arith.c deleted file mode 100644 index b3e45e8e..00000000 --- a/matrix_arith.c +++ /dev/null @@ -1,676 +0,0 @@ - -/* matrix_arith.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 <math.h> -#include <gsl/gsl_matrix.h> -#include <gsl/gsl_blas.h> -#include <gsl/gsl_permutation.h> -#include <gsl/gsl_linalg.h> - -#include "lcomplex.h" -#include "gs-types.h" -#include "matrix.h" -#include "cmatrix.h" -#include "matrix_arith.h" -#include "lua-utils.h" - -struct pmatrix { - int tp; - union { - gsl_matrix *real; - gsl_matrix_complex *cmpl; - } m; -}; - -static int matrix_inv (lua_State *L); -static int matrix_solve (lua_State *L); -static int matrix_dim (lua_State *L); -static int matrix_copy (lua_State *L); -static int matrix_prod (lua_State *L); -static int matrix_dot (lua_State *L); -static int matrix_set (lua_State *L); -static int matrix_newindex (lua_State *L); - -static const struct luaL_Reg matrix_arith_functions[] = { - {"dim", matrix_dim}, - {"copy", matrix_copy}, - {"solve", matrix_solve}, - {"inv", matrix_inv}, - {"prod", matrix_prod}, - {"dot", matrix_dot}, - {"set", matrix_set}, - {NULL, NULL} -}; - -static char const * const genop_dim_err_msg = "matrices should have the same size in %s"; -static char const * const mm_dim_err_msg = "incompatible matrix dimensions in multiplication"; - -static void -check_matrix_mul_dim (lua_State *L, struct pmatrix *a, struct pmatrix *b, - bool atrans, bool btrans) -{ - size_t col1 = (a->tp == GS_MATRIX ? (atrans ? a->m.real->size1 : a->m.real->size2) : (atrans ? a->m.cmpl->size1 : a->m.cmpl->size2)); - size_t row2 = (b->tp == GS_MATRIX ? (btrans ? b->m.real->size2 : b->m.real->size1) : (btrans ? b->m.cmpl->size2 : b->m.cmpl->size1)); - - if (col1 != row2) - luaL_error (L, mm_dim_err_msg); -} - -static gsl_matrix_complex * -push_matrix_complex_of_real (lua_State *L, const gsl_matrix *a) -{ - size_t n1 = a->size1, n2 = a->size2; - gsl_matrix_complex *r = matrix_complex_push_raw (L, n1, n2); - size_t i; - for (i = 0; i < n1; i++) - { - double *rp0 = r->data + 2 * (r->tda * i); - double *rp1 = r->data + 2 * (r->tda * i + n2); - double *ap = a->data + a->tda * i; - double *rp; - for (rp = rp0; rp < rp1; rp += 2, ap += 1) - { - rp[0] = *ap; - rp[1] = 0.0; - } - } - - return r; -} - -static void -check_matrix_type (lua_State *L, int index, struct pmatrix *r) -{ - if (gs_is_userdata (L, index, GS_MATRIX)) - { - r->tp = GS_MATRIX; - r->m.real = lua_touserdata (L, index); - } - else if (gs_is_userdata (L, index, GS_CMATRIX)) - { - r->tp = GS_CMATRIX; - r->m.cmpl = lua_touserdata (L, index); - } - else - { - gs_type_error (L, index, "matrix"); - } -} - -static void -matrix_complex_promote (lua_State *L, int index, struct pmatrix *a) -{ - a->tp = GS_CMATRIX; - a->m.cmpl = push_matrix_complex_of_real (L, a->m.real); - lua_replace (L, index); -} - -#define OPER_ADD -#include "template_matrix_oper_on.h" -#include "matrix_op_source.c" -#include "template_matrix_oper_off.h" -#undef OPER_ADD - -#define OPER_SUB -#include "template_matrix_oper_on.h" -#include "matrix_op_source.c" -#include "template_matrix_oper_off.h" -#undef OPER_SUB - -#define OPER_MUL -#include "template_matrix_oper_on.h" -#include "matrix_op_source.c" -#include "template_matrix_oper_off.h" -#undef OPER_MUL - -#define OPER_DIV -#include "template_matrix_oper_on.h" -#include "matrix_op_source.c" -#include "template_matrix_oper_off.h" -#undef OPER_DIV - -int -matrix_op_mul (lua_State *L) -{ - bool a_is_scalar = lua_iscomplex (L, 1); - bool b_is_scalar = lua_iscomplex (L, 2); - - if (a_is_scalar && b_is_scalar) - { - Complex a = lua_tocomplex (L, 1), b = lua_tocomplex (L, 2); - lua_pushcomplex (L, a * b); - return 1; - } - - if (a_is_scalar) - { - return scalar_matrix_mul (L, 1, 2, true); - } - else if (b_is_scalar) - { - return scalar_matrix_mul (L, 2, 1, true); - } - - struct pmatrix pa, pb; - int rtp; - - check_matrix_type (L, 1, &pa); - check_matrix_type (L, 2, &pb); - - rtp = (pa.tp == GS_MATRIX && pb.tp == GS_MATRIX ? GS_MATRIX : GS_CMATRIX); - - if (pa.tp != rtp) - matrix_complex_promote (L, 1, &pa); - - if (pb.tp != rtp) - matrix_complex_promote (L, 1, &pb); - - if (rtp == GS_MATRIX) - { - const gsl_matrix *a = pa.m.real, *b = pb.m.real; - gsl_matrix *r = matrix_push (L, a->size1, b->size2); - - if (a->size2 != b->size1) - return luaL_error (L, mm_dim_err_msg); - - gsl_blas_dgemm (CblasNoTrans, CblasNoTrans, 1.0, a, b, 1.0, r); - } - else - { - const gsl_matrix_complex *a = pa.m.cmpl, *b = pb.m.cmpl; - gsl_matrix_complex *r = matrix_complex_push (L, a->size1, b->size2); - gsl_complex u = {{1.0, 0.0}}; - - if (a->size2 != b->size1) - return luaL_error (L, mm_dim_err_msg); - - gsl_blas_zgemm (CblasNoTrans, CblasNoTrans, u, a, b, u, r); - } - - return 1; -} - -int -matrix_op_div (lua_State *L) -{ - bool a_is_scalar = lua_iscomplex (L, 1); - bool b_is_scalar = lua_iscomplex (L, 2); - - if (a_is_scalar && b_is_scalar) - { - Complex a = lua_tocomplex (L, 1), b = lua_tocomplex (L, 2); - lua_pushcomplex (L, a / b); - return 1; - } - - if (b_is_scalar) - { - return scalar_matrix_div (L, 2, 1, false); - } - - return luaL_error (L, "cannot divide by a matrix"); -} - -int -matrix_unm (lua_State *L) -{ - struct pmatrix p; - - check_matrix_type (L, 1, &p); - - if (p.tp == GS_MATRIX) - { - const gsl_matrix *a = p.m.real; - size_t n1 = a->size1, n2 = a->size2; - gsl_matrix *r = matrix_push_raw (L, n1, n2); - size_t i; - - for (i = 0; i < n1; i++) - { - double *rp0 = r->data + (r->tda * i); - double *rp1 = r->data + (r->tda * i + n2); - double *ap = a->data + a->tda * i; - double *rp; - - for (rp = rp0; rp < rp1; rp++, ap++) - *rp = - (*ap); - } - } - else - { - const gsl_matrix_complex *a = p.m.cmpl; - size_t n1 = a->size1, n2 = a->size2; - gsl_matrix_complex *r = matrix_complex_push_raw (L, n1, n2); - size_t i; - - for (i = 0; i < n1; i++) - { - double *rp0 = r->data + 2* (r->tda * i); - double *rp1 = r->data + 2* (r->tda * i + n2); - double *ap = a->data + 2* (a->tda * i); - double *rp; - - for (rp = rp0; rp < rp1; rp += 2, ap += 2) - { - rp[0] = - ap[0]; - rp[1] = - ap[1]; - } - } - } - - return 1; -} - -int -matrix_inv (lua_State *L) -{ - struct pmatrix a; - check_matrix_type (L, 1, &a); - switch (a.tp) - { - case GS_MATRIX: - return matrix_inverse_raw (L, a.m.real); - case GS_CMATRIX: - return matrix_complex_inverse_raw (L, a.m.cmpl); - default: - /* */; - } - return 0; -} - -int -matrix_solve (lua_State *L) -{ - struct pmatrix a, b, r; - check_matrix_type (L, 1, &a); - check_matrix_type (L, 2, &b); - - r.tp = (a.tp == GS_MATRIX && b.tp == GS_MATRIX ? GS_MATRIX : GS_CMATRIX); - - if (a.tp != r.tp) - matrix_complex_promote (L, 1, &a); - - if (b.tp != r.tp) - matrix_complex_promote (L, 2, &b); - - switch (r.tp) - { - case GS_MATRIX: - return matrix_solve_raw (L, a.m.real, b.m.real); - case GS_CMATRIX: - return matrix_complex_solve_raw (L, a.m.cmpl, b.m.cmpl); - default: - /* */; - } - - return 0; -} - -int -matrix_dim (lua_State *L) -{ - struct pmatrix p; - check_matrix_type (L, 1, &p); - size_t n1, n2; - - if (p.tp == GS_MATRIX) - { - const gsl_matrix *a = lua_touserdata (L, 1); - n1 = a->size1; - n2 = a->size2; - } - else - { - const gsl_matrix_complex *a = lua_touserdata (L, 1); - n1 = a->size1; - n2 = a->size2; - } - - lua_pushinteger (L, n1); - lua_pushinteger (L, n2); - return 2; -} - -int -matrix_copy (lua_State *L) -{ - struct pmatrix p; - check_matrix_type (L, 1, &p); - - if (p.tp == GS_MATRIX) - { - const gsl_matrix *a = lua_touserdata (L, 1); - gsl_matrix *cp = matrix_push_raw (L, a->size1, a->size2); - gsl_matrix_memcpy (cp, a); - } - else - { - const gsl_matrix_complex *a = lua_touserdata (L, 1); - gsl_matrix_complex *cp = matrix_complex_push_raw (L, a->size1, a->size2); - gsl_matrix_complex_memcpy (cp, a); - } - - return 1; -} - -int -matrix_prod (lua_State *L) -{ - struct pmatrix a, b, r; - gsl_complex u = {{1.0, 0.0}}; - - check_matrix_type (L, 1, &a); - check_matrix_type (L, 2, &b); - - check_matrix_mul_dim (L, &a, &b, true, false); - - r.tp = (a.tp == GS_MATRIX && b.tp == GS_MATRIX ? GS_MATRIX : GS_CMATRIX); - - if (a.tp != r.tp) - matrix_complex_promote (L, 1, &a); - - if (b.tp != r.tp) - matrix_complex_promote (L, 2, &b); - - switch (r.tp) - { - case GS_MATRIX: - r.m.real = matrix_push (L, a.m.real->size2, b.m.real->size2); - gsl_blas_dgemm (CblasTrans, CblasNoTrans, 1.0, a.m.real, b.m.real, 1.0, r.m.real); - break; - case GS_CMATRIX: - r.m.cmpl = matrix_complex_push (L, a.m.cmpl->size2, b.m.cmpl->size2); - gsl_blas_zgemm (CblasConjTrans, CblasNoTrans, u, a.m.cmpl, b.m.cmpl, u, r.m.cmpl); - break; - default: - /* */; - } - - return 1; -} - -int -matrix_dot (lua_State *L) -{ - struct pmatrix a, b, r; - size_t sz; - - check_matrix_type (L, 1, &a); - check_matrix_type (L, 2, &b); - - sz = (a.tp == GS_MATRIX ? a.m.real->size2 : a.m.cmpl->size2); - if (sz != 1) - return gs_type_error (L, 1, "column matrix"); - - sz = (b.tp == GS_MATRIX ? b.m.real->size2 : b.m.cmpl->size2); - if (sz != 1) - return gs_type_error (L, 2, "column matrix"); - - check_matrix_mul_dim (L, &a, &b, true, false); - - r.tp = (a.tp == GS_MATRIX && b.tp == GS_MATRIX ? GS_MATRIX : GS_CMATRIX); - - if (a.tp != r.tp) - matrix_complex_promote (L, 1, &a); - - if (b.tp != r.tp) - matrix_complex_promote (L, 2, &b); - - switch (r.tp) - { - case GS_MATRIX: - { - double result; - gsl_vector_view av = gsl_matrix_column (a.m.real, 0); - gsl_vector_view bv = gsl_matrix_column (b.m.real, 0); - gsl_blas_ddot (&av.vector, &bv.vector, &result); - lua_pushnumber (L, result); - return 1; - } - case GS_CMATRIX: - { - gsl_complex result; - gsl_vector_complex_view av = gsl_matrix_complex_column (a.m.cmpl, 0); - gsl_vector_complex_view bv = gsl_matrix_complex_column (b.m.cmpl, 0); - Complex z; - - gsl_blas_zdotc (&av.vector, &bv.vector, &result); - z = result.dat[0] + _Complex_I * result.dat[1]; - lua_pushcomplex (L, z); - return 1; - } - default: - /* */; - } - - return 0; -} - -int -matrix_set (lua_State *L) -{ - struct pmatrix a, b; - int rtp; - - check_matrix_type (L, 1, &a); - check_matrix_type (L, 2, &b); - - rtp = (a.tp == GS_MATRIX && b.tp == GS_MATRIX ? GS_MATRIX : GS_CMATRIX); - - if (a.tp != rtp) - matrix_complex_promote (L, 1, &a); - - if (b.tp != rtp) - matrix_complex_promote (L, 2, &b); - - switch (rtp) - { - case GS_MATRIX: - { - gsl_matrix *dst = a.m.real, *src = b.m.real; - if (dst->size1 != src->size1 || dst->size2 != src->size2) - luaL_error (L, "matrix dimensions does not match"); - gsl_matrix_memcpy (dst, src); - break; - } - case GS_CMATRIX: - { - gsl_matrix_complex *dst = a.m.cmpl, *src = b.m.cmpl; - if (dst->size1 != src->size1 || dst->size2 != src->size2) - luaL_error (L, "matrix dimensions does not match"); - gsl_matrix_complex_memcpy (dst, src); - break; - } - default: - /* */; - } - - return 0; -} - -static void -matrix_set_row_matrix (lua_State *L, - struct pmatrix *lhs, int index, - struct pmatrix *rhs) -{ - if (rhs->tp == GS_CMATRIX && lhs->tp == GS_MATRIX) - gs_type_error (L, 3, "real matrix"); - - if (lhs->tp == GS_CMATRIX && rhs->tp == GS_MATRIX) - matrix_complex_promote (L, 3, rhs); - - if (lhs->tp == GS_MATRIX) - { - gsl_matrix *a = lhs->m.real, *b = rhs->m.real; - size_t ncol = a->size2; - - if (index >= a->size1 || index < 0) - { - luaL_error (L, INVALID_INDEX_MSG); - } - else - { - gsl_matrix_view ar = gsl_matrix_submatrix (a, index, 0, 1, ncol); - - if (b->size1 != 1 || a->size2 != b->size2) - luaL_error (L, "matrix dimensions does not match"); - - gsl_matrix_memcpy (&ar.matrix, b); - } - } - else - { - gsl_matrix_complex *a = lhs->m.cmpl, *b = rhs->m.cmpl; - size_t ncol = a->size2; - - if (index >= a->size1 || index < 0) - { - luaL_error (L, INVALID_INDEX_MSG); - } - else - { - gsl_matrix_complex_view ar = gsl_matrix_complex_submatrix (a, index, 0, 1, ncol); - - if (b->size1 != 1 || a->size2 != b->size2) - luaL_error (L, "matrix dimensions does not match"); - - gsl_matrix_complex_memcpy (&ar.matrix, b); - } - } -} - -static void -matrix_set_element (lua_State *L, struct pmatrix *lhs, int index, bool vdir1) -{ - size_t i1 = (vdir1 ? index : 0); - size_t i2 = (vdir1 ? 0 : index); - - if (lhs->tp == GS_MATRIX) - { - gsl_matrix *m = lhs->m.real; - size_t n = (vdir1 ? m->size1 : m->size2); - - if (index >= n || index < 0) - luaL_error (L, INVALID_INDEX_MSG); - - if (lua_isnumber (L, 3)) - { - double val = lua_tonumber (L, 3); - gsl_matrix_set (m, i1, i2, val); - } - else - { - gs_type_error (L, 3, "real number"); - } - } - else - { - gsl_matrix_complex *m = lhs->m.cmpl; - size_t n = (vdir1 ? m->size1 : m->size2); - - if (index >= n || index < 0) - luaL_error (L, INVALID_INDEX_MSG); - - if (lua_iscomplex (L, 3)) - { - Complex z = lua_tocomplex (L, 3); - gsl_complex gslz; - GSL_SET_COMPLEX(&gslz, creal(z), cimag(z)); - gsl_matrix_complex_set (m, i1, i2, gslz); - } - else - { - gs_type_error (L, 3, "complex number"); - } - } -} - -static int -matrix_set_row (lua_State *L, struct pmatrix *lhs, int index) -{ - bool vdir1, vdir2; - -#ifdef LUA_INDEX_CONVENTION - index -= 1; -#endif - - if (lhs->tp == GS_MATRIX) - { - gsl_matrix *m = lhs->m.real; - vdir1 = (m->size2 == 1); - vdir2 = (m->size1 == 1); - } - else - { - gsl_matrix_complex *m = lhs->m.cmpl; - vdir1 = (m->size2 == 1); - vdir2 = (m->size1 == 1); - } - - if (!vdir1 && !vdir2) - { - struct pmatrix rhs; - check_matrix_type (L, 3, &rhs); - matrix_set_row_matrix (L, lhs, index, &rhs); - } - else - { - matrix_set_element (L, lhs, index, vdir1); - } - - return 0; -} - -int -matrix_newindex (lua_State *L) -{ - struct pmatrix lhs; - - check_matrix_type (L, 1, &lhs); - - if (lua_isnumber (L, 2)) - { - int index = lua_tointeger (L, 2); - matrix_set_row (L, &lhs, index); - return 0; - } - - return luaL_error (L, "attempt to index matrix with non-integer value"); -} - -void -matrix_arith_register (lua_State *L) -{ - luaL_getmetatable (L, GS_METATABLE(GS_MATRIX)); - lua_getfield (L, -2, "Matrix"); - lua_pushcclosure (L, matrix_newindex, 1); - lua_setfield (L, -2, "__newindex"); - lua_pop (L, 1); - - luaL_getmetatable (L, GS_METATABLE(GS_CMATRIX)); - lua_getfield (L, -2, "cMatrix"); - lua_pushcclosure (L, matrix_newindex, 1); - lua_setfield (L, -2, "__newindex"); - lua_pop (L, 1); - - luaL_register (L, NULL, matrix_arith_functions); -} diff --git a/matrix_arith.h b/matrix_arith.h deleted file mode 100644 index 40825df7..00000000 --- a/matrix_arith.h +++ /dev/null @@ -1,18 +0,0 @@ - -#ifndef MATRIX_ARITH_H -#define MATRIX_ARITH_H - -#include <lua.h> - -extern int matrix_elemop_add (lua_State *L); -extern int matrix_elemop_sub (lua_State *L); -extern int matrix_op_mul (lua_State *L); -extern int matrix_op_div (lua_State *L); -extern int matrix_unm (lua_State *L); - -extern void matrix_arith_register (lua_State *L); - -#define matrix_op_add matrix_elemop_add -#define matrix_op_sub matrix_elemop_sub - -#endif diff --git a/matrix_decls_source.c b/matrix_decls_source.c deleted file mode 100644 index 19cdbd1a..00000000 --- a/matrix_decls_source.c +++ /dev/null @@ -1,54 +0,0 @@ - -/* matrix_decls_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. - */ - -#define NLINFIT_MAX_ITER 30 - -static int FUNCTION (matrix, index) (lua_State *L); -static int FUNCTION (matrix, len) (lua_State *L); -static int FUNCTION (matrix, get_elem) (lua_State *L); -static int FUNCTION (matrix, set_elem) (lua_State *L); -static int FUNCTION (matrix, free) (lua_State *L); -static int FUNCTION (matrix, new) (lua_State *L); -static int FUNCTION (matrix, slice) (lua_State *L); - -static void FUNCTION (matrix, set_ref) (lua_State *L, int index); - -static const struct luaL_Reg FUNCTION (matrix, meta_methods)[] = { - {"__add", matrix_op_add}, - {"__sub", matrix_op_sub}, - {"__mul", matrix_op_mul}, - {"__div", matrix_op_div}, - {"__unm", matrix_unm}, - {"__len", FUNCTION (matrix, len)}, - {"__gc", FUNCTION (matrix, free)}, - {NULL, NULL} -}; - -static const struct luaL_Reg FUNCTION (matrix, methods)[] = { - {"get", FUNCTION (matrix, get_elem)}, - {"set", FUNCTION (matrix, set_elem)}, - {"slice", FUNCTION (matrix, slice)}, - {NULL, NULL} -}; - -static const struct luaL_Reg FUNCTION (matrix, functions)[] = { - {PREFIX "new", FUNCTION (matrix, new)}, - {NULL, NULL} -}; diff --git a/matrix_headers_source.h b/matrix_headers_source.h deleted file mode 100644 index 0515d990..00000000 --- a/matrix_headers_source.h +++ /dev/null @@ -1,46 +0,0 @@ - -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 void FUNCTION (matrix, null_view) (lua_State *L, - int index); - -extern void FUNCTION (matrix, check_size) (lua_State *L, - TYPE (gsl_matrix) *m, - size_t n1, size_t n2); - -extern int FUNCTION (matrix, inverse_raw)(lua_State *L, - const TYPE (gsl_matrix) *a); - -extern int FUNCTION (matrix, solve_raw) (lua_State *L, - const TYPE (gsl_matrix) *a, - const TYPE (gsl_matrix) *b); - - -/* 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, set_view) (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 deleted file mode 100644 index 286b2d70..00000000 --- a/matrix_helper_source.c +++ /dev/null @@ -1,77 +0,0 @@ - -void -FUNCTION (matrix, set_view_and_push) (lua_State *L, int index, double *data, - size_t n1, size_t n2, const double *src) -{ - TYPE (gsl_matrix) *view = FUNCTION (matrix, check) (L, index); - - assert (view->owner == 0); - - view->data = data; - view->size1 = n1; - view->size2 = n2; - view->tda = n2; - view->block = NULL; - - if (src) - memcpy (data, src, n1 * n2 * sizeof(BASE)); - - lua_pushvalue (L, index); -} - -void -FUNCTION (matrix, set_view) (lua_State *L, int index, double *data, - size_t n1, size_t n2, const double *src) -{ - TYPE (gsl_matrix) *view = FUNCTION (matrix, check) (L, index); - - assert (view->owner == 0); - - view->data = data; - view->size1 = n1; - view->size2 = n2; - view->tda = n2; - view->block = NULL; - - if (src) - memcpy (data, src, n1 * n2 * sizeof(BASE)); -} - -static void -TYPE (copy_jacobian_raw) (double *cmpl, double *real, size_t n, size_t p, - size_t mult, bool inverse) -{ - gsl_vector_view dview, sview; - double *cp, *rp; - size_t k, nu; - - for (nu = 0; nu < mult; nu++) - { - cp = cmpl + nu; - rp = real + p*nu; - for (k = 0; k < p; k++, rp += 1, cp += mult) - { - dview = gsl_vector_view_array_with_stride (cp, mult * p, n); - sview = gsl_vector_view_array_with_stride (rp, mult * 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); -} - diff --git a/matrix_op_source.c b/matrix_op_source.c deleted file mode 100644 index 5ce008dc..00000000 --- a/matrix_op_source.c +++ /dev/null @@ -1,168 +0,0 @@ - -/* matrix_op_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. - */ - -static int -SCALAR_MAT_FUNCTION (scalar_matrix) (lua_State *L, int sidx, int midx, - bool direct) -{ - struct pmatrix p; - Complex sc = lua_tocomplex (L, sidx); - int rtp; - - check_matrix_type (L, midx, &p); - - rtp = (cimag(sc) == 0.0 && p.tp == GS_MATRIX ? GS_MATRIX : GS_CMATRIX); - - if (rtp == GS_MATRIX) - { - double s = creal(sc); - gsl_matrix *m = p.m.real; - size_t n1 = m->size1, n2 = m->size2; - gsl_matrix *r = matrix_push_raw (L, n1, n2); - size_t i; - - for (i = 0; i < n1; i++) - { - double *mp0 = m->data + (m->tda * i); - double *mp1 = m->data + (m->tda * i + n2); - double *rp = r->data + r->tda * i; - double *mp; - - for (mp = mp0; mp < mp1; mp++, rp++) - { - double a = (direct ? s : *mp), b = (direct ? *mp : s); - *rp = BASE_OPER(a, b); - } - } - } - else - { - Complex s = sc; - if (p.tp == GS_CMATRIX) - { - gsl_matrix_complex *m = p.m.cmpl; - size_t n1 = m->size1, n2 = m->size2; - gsl_matrix_complex *r = matrix_complex_push_raw (L, n1, n2); - size_t i; - - for (i = 0; i < n1; i++) - { - Complex *mp0 = (Complex *) (m->data + 2 * (m->tda * i)); - Complex *mp1 = (Complex *) (m->data + 2 * (m->tda * i + n2)); - Complex *rp = (Complex *) (r->data + 2 * (r->tda * i)); - Complex *mp; - - for (mp = mp0; mp < mp1; mp++, rp++) - { - Complex a = (direct ? s : *mp), b = (direct ? *mp : s); - *rp = BASE_OPER(a, b); - } - } - } - else - { - gsl_matrix *m = p.m.real; - size_t n1 = m->size1, n2 = m->size2; - gsl_matrix_complex *r = matrix_complex_push_raw (L, n1, n2); - size_t i; - - for (i = 0; i < n1; i++) - { - double *mp0 = m->data + (m->tda * i); - double *mp1 = m->data + (m->tda * i + n2); - Complex *rp = (Complex *) (r->data + 2 * (r->tda * i)); - double *mp; - - for (mp = mp0; mp < mp1; mp++, rp++) - { - Complex mc = *mp; - Complex a = (direct ? s : mc), b = (direct ? mc : s); - *rp = BASE_OPER(a, b); - } - } - } - } - - return 1; -} - -#if (OP_ELEM_DEF == 1) -int -OPER_FUNCTION(matrix_elemop) (lua_State *L) -{ - bool a_is_scalar = lua_iscomplex (L, 1); - bool b_is_scalar = lua_iscomplex (L, 2); - - if (a_is_scalar && b_is_scalar) - { - Complex a = lua_tocomplex (L, 1), b = lua_tocomplex (L, 2); - lua_pushcomplex (L, BASE_OPER(a, b)); - return 1; - } - - if (a_is_scalar) - { - return SCALAR_MAT_FUNCTION (scalar_matrix) (L, 1, 2, true); - } - else if (b_is_scalar) - { - return SCALAR_MAT_FUNCTION (scalar_matrix) (L, 2, 1, false); - } - - struct pmatrix pa, pb; - int rtp; - - check_matrix_type (L, 1, &pa); - check_matrix_type (L, 2, &pb); - - rtp = (pa.tp == GS_MATRIX && pb.tp == GS_MATRIX ? GS_MATRIX : GS_CMATRIX); - - if (pa.tp != rtp) - matrix_complex_promote (L, 1, &pa); - - if (pb.tp != rtp) - matrix_complex_promote (L, 1, &pb); - - if (rtp == GS_MATRIX) - { - const gsl_matrix *a = pa.m.real, *b = pb.m.real; - gsl_matrix *r = matrix_push_raw (L, a->size1, a->size2); - - if (a->size1 != b->size1 || a->size2 != b->size2) - return luaL_error (L, genop_dim_err_msg, OP_NAME); - - gsl_matrix_memcpy (r, a); - OPER_FUNCTION (gsl_matrix) (r, b); - } - else - { - const gsl_matrix_complex *a = pa.m.cmpl, *b = pb.m.cmpl; - gsl_matrix_complex *r = matrix_complex_push_raw (L, a->size1, a->size2); - - if (a->size1 != b->size1 || a->size2 != b->size2) - return luaL_error (L, genop_dim_err_msg, OP_NAME); - - gsl_matrix_complex_memcpy (r, a); - OPER_FUNCTION (gsl_matrix_complex) (r, b); - } - - return 1; -} -#endif diff --git a/matrix_source.c b/matrix_source.c deleted file mode 100644 index 43aeb8b1..00000000 --- a/matrix_source.c +++ /dev/null @@ -1,407 +0,0 @@ - -/* matrix_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. - */ - -static TYPE (gsl_matrix) * -TYPE (_push_matrix) (lua_State *L, int n1, int n2, bool clean) -{ - TYPE (gsl_matrix) *m = lua_newuserdata (L, sizeof(gsl_matrix)); - - if (clean) - m->block = FUNCTION (gsl_block, calloc) ((size_t) n1 * n2); - else - m->block = FUNCTION (gsl_block, alloc) ((size_t) n1 * n2); - - if (m->block == NULL) - luaL_error (L, OUT_OF_MEMORY_MSG); - - m->data = m->block->data; - m->size1 = n1; - m->size2 = n2; - m->tda = n2; - m->owner = 1; - - gs_set_metatable (L, GS_TYPE(MATRIX)); - - return m; -} - -TYPE (gsl_matrix) * -FUNCTION (matrix, push_raw) (lua_State *L, int n1, int n2) -{ - return TYPE (_push_matrix) (L, n1, n2, false); -} - -TYPE (gsl_matrix) * -FUNCTION (matrix, push) (lua_State *L, int n1, int n2) -{ - return TYPE (_push_matrix) (L, n1, n2, true); -} - -TYPE (gsl_matrix) * -FUNCTION (matrix, check) (lua_State *L, int index) -{ - return gs_check_userdata (L, index, GS_TYPE(MATRIX)); -} - -void -FUNCTION (matrix, check_size) (lua_State *L, TYPE (gsl_matrix) *m, - size_t n1, size_t n2) -{ - if (m->size1 != n1 || m->size2 != n2) - { - luaL_error (L, "expecting matrix of size %dx%d, got matrix %dx%d", - n1, n2, m->size1, m->size2); - } -} - -void -FUNCTION (matrix, push_view) (lua_State *L, TYPE (gsl_matrix) *m) -{ - TYPE (gsl_matrix) *mview; - - mview = lua_newuserdata (L, sizeof(TYPE (gsl_matrix))); - - if (m) - { - *mview = *m; - } - else - { - mview->data = NULL; - mview->block = NULL; - } - - mview->owner = 0; - - gs_set_metatable (L, GS_TYPE(MATRIX)); -} - -void -FUNCTION (matrix, null_view) (lua_State *L, int index) -{ - TYPE (gsl_matrix) *m = FUNCTION (matrix, check) (L, index); - assert (m->owner == 0); - m->data = NULL; -} - -void -FUNCTION (matrix, set_ref) (lua_State *L, int index) -{ - lua_newtable (L); - lua_pushvalue (L, index); - lua_rawseti (L, -2, 1); - lua_setfenv (L, -2); -} - -int -FUNCTION(matrix, slice) (lua_State *L) -{ - TYPE (gsl_matrix) *a = FUNCTION (matrix, check) (L, 1); - lua_Integer k1 = luaL_checkinteger (L, 2), k2 = luaL_checkinteger (L, 3); - lua_Integer n1 = luaL_checkinteger (L, 4), n2 = luaL_checkinteger (L, 5); - VIEW (gsl_matrix) view; - -#ifdef LUA_INDEX_CONVENTION - k1 -= 1; - k2 -= 1; -#endif - - if (k1 < 0 || k2 < 0 || n1 < 0 || n2 < 0) - luaL_error (L, INVALID_INDEX_MSG); - - if (k1 >= a->size1 || k2 >= a->size2 || - k1 + n1 > a->size1 || k2 + n2 > a->size2) - luaL_error (L, INVALID_INDEX_MSG); - - view = FUNCTION (gsl_matrix, submatrix) (a, k1, k2, n1, n2); - FUNCTION (matrix, push_view) (L, &view.matrix); - FUNCTION (matrix, set_ref) (L, 1); - - return 1; -} - -int -FUNCTION(matrix, get_elem) (lua_State *L) -{ - const TYPE (gsl_matrix) *m = FUNCTION (matrix, check) (L, 1); - lua_Integer r = luaL_checkinteger (L, 2); - lua_Integer c = luaL_checkinteger (L, 3); - LUA_TYPE v; - BASE gslval; - -#ifdef LUA_INDEX_CONVENTION - r -= 1; - c -= 1; -#endif - - if (r < 0 || r >= (int) m->size1 || c < 0 || c >= (int) m->size2) - return 0; - - gslval = FUNCTION (gsl_matrix, get) (m, (size_t) r, (size_t) c); - v = TYPE (value_retrieve) (gslval); - - LUA_FUNCTION (push) (L, v); - - return 1; -} - -int -FUNCTION(matrix, set_elem) (lua_State *L) -{ - TYPE (gsl_matrix) *m = FUNCTION (matrix, check) (L, 1); - lua_Integer r = luaL_checkinteger (L, 2); - lua_Integer c = luaL_checkinteger (L, 3); - LUA_TYPE v = LUAL_FUNCTION(check) (L, 4); - BASE gslval; - -#ifdef LUA_INDEX_CONVENTION - r -= 1; - c -= 1; -#endif - - luaL_argcheck (L, r >= 0 && r < (int) m->size1, 2, - "row index out of limits"); - luaL_argcheck (L, c >= 0 && c < (int) m->size2, 3, - "column index out of limits"); - - gslval = TYPE (value_assign) (v); - FUNCTION (gsl_matrix, set) (m, (size_t) r, (size_t) c, gslval); - - return 0; -} - -int -FUNCTION(matrix, free) (lua_State *L) -{ - TYPE (gsl_matrix) *m = FUNCTION (matrix, check) (L, 1); - if (m->owner) - { - assert (m->block); - FUNCTION (gsl_block, free) (m->block); - } - return 0; -}; - -int -FUNCTION(matrix, new) (lua_State *L) -{ - lua_Integer nr = luaL_checkinteger (L, 1); - lua_Integer nc = luaL_checkinteger (L, 2); - - luaL_argcheck (L, nr > 0, 1, "row number should be positive"); - luaL_argcheck (L, nc > 0, 2, "column number should be positive"); - - if (lua_gettop (L) == 3) - { - if (lua_isfunction (L, 3)) - { - TYPE (gsl_matrix) *m; - size_t i, j; - - m = FUNCTION (matrix, push) (L, (size_t) nr, (size_t) nc); - - for (i = 0; i < nr; i++) - { - for (j = 0; j < nc; j++) - { -#ifdef LUA_INDEX_CONVENTION - size_t ig = i+1, jg = j+1; -#else - size_t ig = i, jg = j; -#endif - LUA_TYPE z; - BASE gslz; - - lua_pushvalue (L, 3); - lua_pushnumber (L, ig); - lua_pushnumber (L, jg); - lua_call (L, 2, 1); - - if (! LUA_FUNCTION(is) (L, 5)) - return luaL_error (L, "expecting real or complex number" \ - " in matrix initialization"); - z = LUA_FUNCTION(to) (L, 5); - gslz = TYPE (value_assign) (z); - FUNCTION (gsl_matrix, set) (m, i, j, gslz); - - lua_pop (L, 1); - } - } - return 1; - } - } - - FUNCTION (matrix, push) (L, (size_t) nr, (size_t) nc); - return 1; -} - -int -FUNCTION(matrix, inverse_raw) (lua_State *L, const TYPE (gsl_matrix) *a) -{ - TYPE (gsl_matrix) *lu, *inverse; - gsl_permutation *p; - size_t n = a->size1; - int status, sign; - - if (a->size2 != n) - gs_type_error (L, 1, "square matrix"); - - p = gsl_permutation_alloc (n); - - lu = FUNCTION (gsl_matrix, alloc) (n, n); - FUNCTION (gsl_matrix, memcpy) (lu, a); - - inverse = FUNCTION (matrix, push_raw) (L, n, n); - - FUNCTION (gsl_linalg, LU_decomp) (lu, p, &sign); - status = FUNCTION (gsl_linalg, LU_invert) (lu, p, inverse); - - FUNCTION (gsl_matrix, free) (lu); - gsl_permutation_free (p); - - if (status) - { - return luaL_error (L, "error during matrix inversion: %s", - gsl_strerror (status)); - } - - return 1; -} -int -FUNCTION(matrix, solve_raw) (lua_State *L, - const TYPE (gsl_matrix) *a, - const TYPE (gsl_matrix) *b) -{ - TYPE (gsl_matrix) *x; - CONST_VIEW (gsl_vector) b_view = CONST_FUNCTION (gsl_matrix, column) (b, 0); - VIEW (gsl_vector) x_view; - TYPE (gsl_matrix) *lu; - gsl_permutation *p; - size_t n = a->size1; - int sign; - - if (a->size2 != n) - gs_type_error (L, 1, "square matrix"); - if (b->size2 != 1) - gs_type_error (L, 1, "vector"); - if (b->size1 != n) - luaL_error (L, "dimensions of vector does not match"); - - x = FUNCTION (matrix, push_raw) (L, n, 1); - x_view = FUNCTION (gsl_matrix, column) (x, 0); - - p = gsl_permutation_alloc (n); - - lu = FUNCTION (gsl_matrix, alloc) (n, n); - FUNCTION (gsl_matrix, memcpy) (lu, a); - - FUNCTION (gsl_linalg, LU_decomp) (lu, p, &sign); - FUNCTION (gsl_linalg, LU_solve) (lu, p, & b_view.vector, & x_view.vector); - - FUNCTION (gsl_matrix, free) (lu); - gsl_permutation_free (p); - - return 1; -} - -int -FUNCTION(matrix, len) (lua_State *L) -{ - TYPE (gsl_matrix) *m = FUNCTION (matrix, check) (L, 1); - lua_pushinteger (L, m->size1); - return 1; -} - -int -FUNCTION(matrix, index) (lua_State *L) -{ - TYPE (gsl_matrix) *m = FUNCTION (matrix, check) (L, 1); - - if (lua_isnumber (L, 2)) - { - int index = lua_tointeger (L, 2); - bool vdir1, vdir2; - -#ifdef LUA_INDEX_CONVENTION - index -= 1; -#endif - - vdir1 = (m->size2 == 1); - vdir2 = (m->size1 == 1); - - if (!vdir1 && !vdir2) - { - VIEW (gsl_matrix) view; - - if (index >= m->size1 || index < 0) - return luaL_error (L, INVALID_INDEX_MSG); - - view = FUNCTION (gsl_matrix, submatrix) (m, index, 0, 1, m->size2); - FUNCTION (matrix, push_view) (L, &view.matrix); - FUNCTION (matrix, set_ref) (L, 1); - } - else - { - int n = (vdir1 ? m->size1 : m->size2); - BASE gslval; - LUA_TYPE v; - - if (index >= n || index < 0) - return luaL_error (L, INVALID_INDEX_MSG); - - if (vdir1) - gslval = FUNCTION (gsl_matrix, get) (m, (size_t) index, 0); - else - gslval = FUNCTION (gsl_matrix, get) (m, 0, (size_t) index); - - v = TYPE (value_retrieve) (gslval); - LUA_FUNCTION (push) (L, v); - } - - return 1; - } - - lua_pushvalue (L, lua_upvalueindex(1)); - lua_replace (L, 1); - lua_rawget (L, 1); - - return 1; -} - -/* register matrix methods in a table (module) gives in the stack */ -void -FUNCTION (matrix, register) (lua_State *L) -{ - luaL_newmetatable (L, GS_TYPENAME(MATRIX)); - lua_newtable (L); - luaL_register (L, NULL, FUNCTION (matrix, methods)); - - lua_pushvalue (L, -1); - lua_pushcclosure (L, FUNCTION (matrix, index), 1); - lua_setfield (L, -3, "__index"); - - lua_insert (L, -2); - luaL_register (L, NULL, FUNCTION (matrix, meta_methods)); - lua_pop (L, 1); - - lua_setfield (L, -2, PREFIX "Matrix"); - - luaL_register (L, NULL, FUNCTION (matrix, functions)); -} @@ -1,6 +1,4 @@ -require 'cmatrix' - local template = require 'template' function gsl.ode(spec) @@ -38,8 +36,8 @@ function gsl.ode(spec) if k == 't' then return s._state.t end if k == 'y' then if not s._sync then - for k=1, s.dim do - s._y[k] = s._state.y[k-1] + for k=0, s.dim-1 do + s._y[k] = s._state.y[k] end s._sync = true end @@ -90,52 +88,3 @@ function gsl.nlinfit(spec) return s end - -function gsl.ode_iter(s, t0, y0, t1, h, tstep) - s:set(t0, y0, h) - return function() - local t, y = s.t, s.y - if t < t1 then - s:evolve(t1, tstep) - return t, y - end - end -end - -local function hc_reduce(hc, f, accu) - local n = hc.length - for i=0, n do accu = f(accu, hc:get(i)) end - return accu -end - -local function hc_print(hc) - local eps = 1e-8 * hc_reduce(hc, function(p,z) return p + csqr(z) end, 0) - local f = function(p, z) - insert(p, fmt('%6i: %s', #p, tostring_eps(z, eps))) - return p - end - return cat(hc_reduce(hc, f, {}), '\n') -end - -gsl.FFT_hc_mixed_radix.__tostring = hc_print -gsl.FFT_hc_radix2.__tostring = hc_print - -function gsl.linmodel(f, x) - local p, n = #f(x[1]), matrix.dim(x) - local A = matrix.new(n, p) - for k=1,n do - local y = f(x[k]) - for j=1,p do A:set(k, j, y[j]) end - end - return A -end - -function gsl.linfit(gener, x, y, w) - local X = gsl.linmodel(gener, x) - local c, cov = gsl.mlinear(X, y, w) - local f = function(xe) - local xs = matrix.vec(gener(xe)) - return matrix.prod(xs, c)[1] - end - return f, c -end diff --git a/mlinear.c b/mlinear.c deleted file mode 100644 index fe489181..00000000 --- a/mlinear.c +++ /dev/null @@ -1,70 +0,0 @@ - -#include <lua.h> -#include <lauxlib.h> -#include <assert.h> -#include <gsl/gsl_vector.h> -#include <gsl/gsl_matrix.h> -#include <gsl/gsl_multifit.h> - -#include "gs-types.h" -#include "matrix.h" -#include "cmatrix.h" -#include "mlinear.h" - -static int mlinear_fit (lua_State *L); - -static const struct luaL_Reg mlinear_functions[] = { - {"mlinear", mlinear_fit}, - {NULL, NULL} -}; - -int -mlinear_fit (lua_State *L) -{ - gsl_multifit_linear_workspace *ws; - gsl_matrix *X = matrix_check (L, 1); - gsl_matrix *y = matrix_check (L, 2); - gsl_matrix *w = NULL; - gsl_vector_const_view yv = gsl_matrix_const_column (y, 0); - size_t n = X->size1, p = X->size2; - gsl_matrix *c, *cov; - gsl_vector_view cv; - double chisq; - int status; - - if (lua_gettop(L) > 2 && !lua_isnil (L, 3)) - { - w = matrix_check (L, 3); - matrix_check_size (L, w, n, 1); - } - - matrix_check_size (L, y, n, 1); - - ws = gsl_multifit_linear_alloc (n, p); - - c = matrix_push_raw (L, p, 1); - cov = matrix_push_raw (L, p, p); - - cv = gsl_matrix_column (c, 0); - - if (w == NULL) - status = gsl_multifit_linear (X, &yv.vector, &cv.vector, cov, &chisq, ws); - else - { - gsl_vector_const_view wv = gsl_matrix_const_column (w, 0); - status = gsl_multifit_wlinear (X, &wv.vector, &yv.vector, - &cv.vector, cov, &chisq, ws); - } - - gsl_multifit_linear_free (ws); - - gs_gsl_errorcheck (L, "multilinear fit", status); - - return 2; -} - -void -mlinear_register (lua_State *L) -{ - luaL_register (L, NULL, mlinear_functions); -} diff --git a/mlinear.h b/mlinear.h deleted file mode 100644 index ec37609b..00000000 --- a/mlinear.h +++ /dev/null @@ -1,8 +0,0 @@ -#ifndef MLINEAR_SYSTEM_H -#define MLINEAR_SYSTEM_H - -#include <lua.h> - -extern void mlinear_register (lua_State *L); - -#endif diff --git a/multimin.c b/multimin.c deleted file mode 100644 index 966f8135..00000000 --- a/multimin.c +++ /dev/null @@ -1,46 +0,0 @@ - -/* multimin.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 <math.h> - -#include <lua.h> -#include "gs-types.h" -#include "multimin.h" - -static const struct luaL_Reg multimin_functions[] = { - {"fdfmultimin", fdfmultimin_new}, - {"fmultimin", fmultimin_new}, - {"gradcheck", gradient_auto_check}, - {NULL, NULL} -}; - -void -multimin_register (lua_State *L) -{ - luaL_newmetatable (L, GS_METATABLE(GS_FMULTIMIN)); - luaL_register (L, NULL, fmultimin_metatable); - lua_pop (L, 1); - - luaL_newmetatable (L, GS_METATABLE(GS_FDFMULTIMIN)); - luaL_register (L, NULL, fdfmultimin_metatable); - lua_pop (L, 1); - - luaL_register (L, NULL, multimin_functions); -} diff --git a/multimin.h b/multimin.h deleted file mode 100644 index 9b35f8ae..00000000 --- a/multimin.h +++ /dev/null @@ -1,40 +0,0 @@ - -/* multimin.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 MULTIMIN_H -#define MULTIMIN_H - -#include <lua.h> -#include <lauxlib.h> - -#include <gsl/gsl_vector.h> - -#include "defs.h" - -extern void multimin_register (lua_State *L); -extern int fdfmultimin_new (lua_State *L); -extern int fmultimin_new (lua_State *L); -extern int gradient_auto_check (lua_State *L); - - -extern const struct luaL_Reg fmultimin_metatable[]; -extern const struct luaL_Reg fdfmultimin_metatable[]; - -#endif diff --git a/nlinfit.c b/nlinfit.c deleted file mode 100644 index e9d3fb4d..00000000 --- a/nlinfit.c +++ /dev/null @@ -1,39 +0,0 @@ - -/* nlinfit.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 <gsl/gsl_vector.h> -#include <gsl/gsl_matrix.h> -#include <gsl/gsl_multifit_nlin.h> - -#include "gs-types.h" -#include "matrix.h" -#include "nlinfit_helper.h" -#include "lua-utils.h" - -#define BASE_DOUBLE -#include "template_matrix_on.h" -#include "nlinfit_decls_source.c" -#include "nlinfit_source.c" -#include "template_matrix_off.h" - -#undef BASE_DOUBLE diff --git a/nlinfit.h b/nlinfit.h deleted file mode 100644 index d7e9e3ea..00000000 --- a/nlinfit.h +++ /dev/null @@ -1,26 +0,0 @@ - -/* nlinfit.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 NLINFIT_H -#define NLINFIT_H - -extern void solver_register (lua_State *L); - -#endif diff --git a/nlinfit_decls_source.c b/nlinfit_decls_source.c deleted file mode 100644 index f0c70c75..00000000 --- a/nlinfit_decls_source.c +++ /dev/null @@ -1,55 +0,0 @@ - -/* nlinfit_decls_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. - */ - -/* declaration of lua function for solver methods */ -static int FUNCTION (solver, new) (lua_State *); -static int FUNCTION (solver, free) (lua_State *); -static int FUNCTION (solver, run) (lua_State *); -static int FUNCTION (solver, iterate) (lua_State *); -static int FUNCTION (solver, index) (lua_State *); -static int FUNCTION (solver, covar) (lua_State *); -static int FUNCTION (solver, get_p) (lua_State *); -static int FUNCTION (solver, get_f) (lua_State *); -static int FUNCTION (solver, get_jacob) (lua_State *); - -static const struct luaL_Reg FUNCTION (solver, metatable)[] = { - {"__gc", FUNCTION (solver, free)}, - {"__index", FUNCTION (solver, index)}, - {NULL, NULL} -}; - -static const struct luaL_Reg FUNCTION (solver, methods)[] = { - {"iterate", FUNCTION (solver, iterate)}, - {"run", FUNCTION (solver, run)}, - {NULL, NULL} -}; - -static const struct luaL_Reg FUNCTION (solver, properties)[] = { - {"covar", FUNCTION (solver, covar)}, - {"p", FUNCTION (solver, get_p)}, - {"f", FUNCTION (solver, get_f)}, - {"J", FUNCTION (solver, get_jacob)}, - {NULL, NULL} -}; - -static const struct luaL_Reg FUNCTION (solver, functions)[] = { - {PREFIX "nlfsolver", FUNCTION (solver, new)}, - {NULL, NULL} -}; diff --git a/nlinfit_helper.c b/nlinfit_helper.c deleted file mode 100644 index ef041805..00000000 --- a/nlinfit_helper.c +++ /dev/null @@ -1,52 +0,0 @@ - -/* solver-impl.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_vector.h> -#include <gsl/gsl_matrix.h> - -#include "nlinfit_helper.h" -#include "matrix.h" - -size_t -check_positive_arg (lua_State *L, const char *name, const char *fullname) -{ - int nb; - lua_getfield (L, 1, name); - nb = luaL_checkinteger (L, -1); - if (nb <= 0) - luaL_error (L, "the number of %s should be a number > 0", fullname); - lua_pop (L, 1); - return (size_t) nb; -} - -void -solver_get_p0 (lua_State *L, gsl_vector_view *p0, size_t *p) -{ - gsl_matrix *m; - lua_getfield (L, 1, "p0"); - m = matrix_check (L, -1); - if (m->size2 != 1) - luaL_error (L, "p0 should be a column matrix"); - *p0 = gsl_matrix_column (m, 0); - *p = m->size1; - lua_pop (L, 1); -} diff --git a/nlinfit_helper.h b/nlinfit_helper.h deleted file mode 100644 index e031a979..00000000 --- a/nlinfit_helper.h +++ /dev/null @@ -1,35 +0,0 @@ - -/* solver-impl.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. - */ - -#ifndef SOLVER_IMPL_H -#define SOLVER_IMPL_H - -#include "defs.h" - -#include <lua.h> -#include <gsl/gsl_vector.h> - -extern size_t -check_positive_arg (lua_State *L, const char *name, const char *fullname); - -extern void -solver_get_p0 (lua_State *L, gsl_vector_view *p0, size_t *p); - -#endif diff --git a/nlinfit_source.c b/nlinfit_source.c deleted file mode 100644 index 33cb029e..00000000 --- a/nlinfit_source.c +++ /dev/null @@ -1,363 +0,0 @@ - -/* nlinfit_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. - */ - -struct params { - lua_State *L; - size_t n; -}; - -struct solver { - gsl_multifit_fdfsolver *base; - gsl_multifit_function_fdf fdf[1]; - struct params p[1]; -}; - -static int -FUNCTION (solver, fdf_hook) (const gsl_vector * x, void * params, - gsl_vector * f, gsl_matrix * J); - -static int -FUNCTION (solver, f_hook) (const gsl_vector * x, void * params, - gsl_vector * f); - -static int -FUNCTION (solver, df_hook) (const gsl_vector * x, void * params, - gsl_matrix * J); - -#define NLINFIT_MAX_ITER 20 - -struct solver * -FUNCTION (solver, push) (lua_State *L, size_t n, size_t p, - int findex) -{ - gsl_multifit_fdfsolver_type const * const T = gsl_multifit_fdfsolver_lmsder; - struct solver *s = gs_new_object (sizeof (struct solver), L, GS_TYPE(NLINFIT)); - - s->base = gsl_multifit_fdfsolver_alloc (T, n, p); - - if (s->base == NULL) - luaL_error (L, OUT_OF_MEMORY_MSG); - - lua_newtable (L); - - lua_pushvalue (L, findex); - lua_rawseti (L, -2, 1); // 1 <- fdf lua function - - matrix_push_raw (L, p, 1); - lua_rawseti (L, -2, 2); // 2 <- x - - FUNCTION (matrix, push_view) (L, NULL); - lua_rawseti (L, -2, 3); // 3 <- f - - FUNCTION (matrix, push_view) (L, NULL); - lua_rawseti (L, -2, 4); // 4 <- f - - /* if multiplicity is 1 this is not really needed */ - matrix_push_raw (L, n, p); - lua_rawseti (L, -2, 5); // 5 <- Jbuffer - - lua_setfenv (L, -2); - - return s; -} - -static void -fcall_args_prepare (lua_State *L, int nargs) -{ - int j, index; - lua_getfenv (L, -1); - index = lua_gettop (L); - for (j = 1; j <= nargs; j++) - lua_rawgeti (L, index, j); - lua_remove (L, index); -} - -int -FUNCTION (solver, new) (lua_State *L) -{ - struct solver *s; - gsl_vector_view p0; - size_t n, nreal, p; - int findex; - - luaL_checktype (L, 1, LUA_TTABLE); - - n = check_positive_arg (L, "n", "observations"); - nreal = n * MULTIPLICITY; - - solver_get_p0 (L, &p0, &p); - - lua_settop (L, 1); - lua_getfield (L, 1, "fdf"); - findex = 2; - - s = FUNCTION (solver, push) (L, nreal, p, findex); - lua_replace (L, 1); - lua_settop (L, 1); - - s->p->L = L; - s->p->n = n; - - s->fdf->f = & FUNCTION (solver, f_hook); - s->fdf->df = & FUNCTION (solver, df_hook); - s->fdf->fdf = & FUNCTION (solver, fdf_hook); - s->fdf->n = nreal; - s->fdf->p = p; - s->fdf->params = s->p; - - fcall_args_prepare (L, 5); - - gsl_multifit_fdfsolver_set (s->base, s->fdf, & p0.vector); - - lua_pop (L, 5); - - return 1; -} - -static struct solver * -FUNCTION (solver, check) (lua_State *L, int index) -{ - return gs_check_userdata (L, index, GS_TYPE(NLINFIT)); -} - -int -FUNCTION (solver, free) (lua_State *L) -{ - struct solver *s = FUNCTION (solver, check) (L, 1); - lua_getfenv (L, 1); - lua_rawgeti (L, -1, 3); - FUNCTION (matrix, null_view) (L, -1); - lua_pop (L, 1); - lua_rawgeti (L, -1, 4); - FUNCTION (matrix, null_view) (L, -1); - gsl_multifit_fdfsolver_free (s->base); - return 0; -} - -int -FUNCTION (solver, fdf_hook) (const gsl_vector * x, void * _params, - gsl_vector * f, gsl_matrix * J) -{ - struct params *params = _params; - lua_State *L = params->L; - size_t n = params->n, p = x->size; - size_t nargs = 2; - gsl_vector_view vv; - gsl_matrix *xm; - gsl_matrix *jraw; - - lua_pushvalue (L, 2); - - xm = matrix_check (L, 3); - vv = gsl_matrix_column (xm, 0); - gsl_vector_memcpy (&vv.vector, x); - lua_pushvalue (L, 3); - - if (f) - FUNCTION (matrix, set_view_and_push) (L, 4, f->data, n, 1, NULL); - else - lua_pushnil (L); - - if (J) - { - double *jptr; - if (MULTIPLICITY >= 2) - { - jraw = matrix_check (L, 6); - jptr = jraw->data; - } - else - jptr = J->data; - - FUNCTION (matrix, set_view_and_push) (L, 5, jptr, n, p, NULL); - nargs ++; - } - - lua_call (L, nargs, 0); - -#if MULTIPLICITY >= 2 - if (J) - { - double *dst = J->data, *src = jraw->data; - FUNCTION (matrix, jacob_copy_cmpl_to_real) (dst, src, n, p, MULTIPLICITY); - } -#endif - - return GSL_SUCCESS; -} - -int -FUNCTION (solver, f_hook) (const gsl_vector * x, void * params, gsl_vector * f) -{ - return FUNCTION (solver, fdf_hook) (x, params, f, NULL); -} - -int -FUNCTION (solver, df_hook) (const gsl_vector * x, void * params, gsl_matrix * J) -{ - return FUNCTION (solver, fdf_hook) (x, params, NULL, J); -} - -int -FUNCTION (solver, iterate) (lua_State *L) -{ - struct solver *s = FUNCTION (solver, check) (L, 1); - int status; - - mlua_null_cache (L, 1); - lua_settop (L, 1); - - fcall_args_prepare (L, 5); - - status = gsl_multifit_fdfsolver_iterate (s->base); - - lua_pop (L, 5); - - if (status) - { - return luaL_error (L, "error during non-linear fit: %s", - gsl_strerror (status)); - } - - status = gsl_multifit_test_delta (s->base->dx, s->base->x, 1e-4, 1e-4); - - if (status == GSL_CONTINUE) - lua_pushstring (L, "continue"); - else - lua_pushstring (L, "success"); - - return 1; -} - -int -FUNCTION (solver, run) (lua_State *L) -{ - struct solver *s; - lua_Integer max_iter, iter = 0; - int iter_status, fit_status; - - s = FUNCTION (solver, check) (L, 1); - - max_iter = (lua_isnumber (L, 2) ? lua_tointeger (L, 2) : NLINFIT_MAX_ITER); - - mlua_null_cache (L, 1); - lua_settop (L, 1); - - fcall_args_prepare (L, 5); - - do - { - iter ++; - - iter_status = gsl_multifit_fdfsolver_iterate (s->base); - - if (iter_status) - { - return luaL_error (L, "error during non-linear fit: %s", - gsl_strerror (iter_status)); - } - - fit_status = gsl_multifit_test_delta (s->base->dx, s->base->x, 1e-4, 1e-4); - } - while (fit_status == GSL_CONTINUE && iter < max_iter); - - lua_pop (L, 5); - - return 0; -} - -int -FUNCTION (solver, covar) (lua_State *L) -{ - struct solver *s; - gsl_matrix *covar; - size_t p; - - s = FUNCTION (solver, check) (L, 1); - p = MULTIPLICITY * s->fdf->p; - covar = matrix_push (L, p, p); - gsl_multifit_covar (s->base->J, 0.0, covar); - return 1; -} - -int -FUNCTION (solver, get_p) (lua_State *L) -{ - struct solver *s = FUNCTION (solver, check) (L, 1); - gsl_vector *src = s->base->x; - size_t P = s->fdf->p; - gsl_matrix *p; - - p = matrix_push (L, P, 1); - gsl_matrix_set_col (p, 0, src); - - return 1; -} - -int -FUNCTION (solver, get_f) (lua_State *L) -{ - struct solver *s = FUNCTION (solver, check) (L, 1); - size_t n = s->fdf->n / MULTIPLICITY; - gsl_vector *src = s->base->f; - TYPE (gsl_matrix) *f; - VIEW (gsl_matrix) fview; - - f = FUNCTION (matrix, push) (L, n, 1); - fview = FUNCTION (gsl_matrix, view_array) (src->data, n, 1); - FUNCTION (gsl_matrix, memcpy) (f, & fview.matrix); - return 1; -} - -int -FUNCTION (solver, get_jacob) (lua_State *L) -{ - struct solver *s = FUNCTION (solver, check) (L, 1); - size_t n = s->fdf->n / MULTIPLICITY, p = s->fdf->p; - gsl_matrix *src = s->base->J; - TYPE (gsl_matrix) *m; - - m = FUNCTION (matrix, push) (L, n, p); - - FUNCTION (matrix, jacob_copy_real_to_cmpl) (m->data, src->data, n, p, - MULTIPLICITY); - - return 1; -} - -int -FUNCTION (solver, index) (lua_State *L) -{ - return mlua_index_with_properties (L, - FUNCTION (solver, properties), - FUNCTION (solver, methods), true); -} - - -void -FUNCTION (solver, register) (lua_State *L) -{ - luaL_newmetatable (L, GS_TYPENAME(NLINFIT)); - luaL_register (L, NULL, FUNCTION (solver, metatable)); - lua_setfield (L, -2, PREFIX "Solver"); - - /* gsl module registration */ - luaL_register (L, NULL, FUNCTION (solver, functions)); -} diff --git a/num/lmfit.lua.in b/num/lmfit.lua.in index 0cd80cf1..b7a997d8 100644 --- a/num/lmfit.lua.in +++ b/num/lmfit.lua.in @@ -1,6 +1,4 @@ -require 'cmatrix' - local abs, min, max, sqrt = math.abs, math.min, math.max, math.sqrt local ffi = require 'ffi' @@ -1,38 +0,0 @@ - -/* 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 <math.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 +0,0 @@ - -/* 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 deleted file mode 100644 index 6ad57ffa..00000000 --- a/ode_decls_source.c +++ /dev/null @@ -1,52 +0,0 @@ - -#define ODE_DEFAULT_EPS_ABS 1e-4 -#define ODE_DEFAULT_EPS_REL 0.0 -/* #define ODE_DEFAULT_STEP 0.1 */ -#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, free) (lua_State *L); - -static int FUNCTION (ode, get_t) (lua_State *L); -static int FUNCTION (ode, get_y) (lua_State *L); - -static struct solver * FUNCTION (ode, check) (lua_State *L, int index); - -enum fenv_e { - FENV_Y = 0, - FENV_Y_BUFFER, - FENV_DYDT, - FENV_DFDY, - FENV_DFDY_BUFFER, - FENV_DFDT, - FENV_F, - FENV_DF, -}; - -static const struct luaL_Reg FUNCTION (ode, metatable)[] = { - {"__gc", FUNCTION (ode, free)}, - {"__index", FUNCTION (ode, index)}, - {NULL, NULL} -}; - -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 struct solver_type TYPE (ode_solver_type)[1] = {{"GSL." PREFIX "ode"}}; diff --git a/ode_solver.c b/ode_solver.c deleted file mode 100644 index 603d1b33..00000000 --- a/ode_solver.c +++ /dev/null @@ -1,71 +0,0 @@ - -#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; - bool needs_jacobian; -}; - -#define ODEIV_METHOD(t) #t, & gsl_odeiv_step_ ## t - -static struct method_entry methods_table[] = { - {ODEIV_METHOD (rk2), false}, - {ODEIV_METHOD (rk4), false}, - {ODEIV_METHOD (rkf45), false}, - {ODEIV_METHOD (rkck), false}, - {ODEIV_METHOD (rk8pd), false}, - {ODEIV_METHOD (rk2imp), false}, - {ODEIV_METHOD (rk4imp), false}, - {ODEIV_METHOD (bsimp), true}, - {ODEIV_METHOD (gear1), false}, - {ODEIV_METHOD (gear2), false}, - {NULL, NULL} -}; -#undef ODEIV_METHOD - -struct solver * -ode_solver_push_new (lua_State *L, const gsl_odeiv_step_type *type, - size_t dim, double eps_abs, double eps_rel, - struct solver_type *st, bool use_jacob) -{ - struct solver *s; - - s = lua_newuserdata (L, sizeof (struct solver)); - - s->step = gsl_odeiv_step_alloc (type, dim); - - if (use_jacob) - s->ctrl = gsl_odeiv_control_standard_new (eps_abs, eps_rel, 1.0, 1.0); - else - s->ctrl = gsl_odeiv_control_y_new (eps_abs, eps_rel); - - s->evol = gsl_odeiv_evolve_alloc (dim); - - s->dimension = dim; - - luaL_getmetatable (L, st->metatable_name); - lua_setmetatable (L, -2); - - return s; -} - -const gsl_odeiv_step_type * -method_lookup (const char *method, const gsl_odeiv_step_type *default_type, - bool *needs_jacobian) -{ - const struct method_entry *p; - for (p = methods_table; p->name; p++) - { - if (strcmp (p->name, method) == 0) - { - *needs_jacobian = p->needs_jacobian; - return *(p->method_type); - } - } - return default_type; -} diff --git a/ode_solver.h b/ode_solver.h deleted file mode 100644 index b0fd7934..00000000 --- a/ode_solver.h +++ /dev/null @@ -1,59 +0,0 @@ - -/* 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 params { - lua_State *L; - size_t n; -}; - -struct solver { - gsl_odeiv_step * step; - gsl_odeiv_control * ctrl; - gsl_odeiv_evolve * evol; - - gsl_odeiv_system system[1]; - struct params params[1]; - - double t, h; - size_t dimension; -}; - -struct solver_type { - const char * const metatable_name; -}; - -extern struct solver * -ode_solver_push_new (lua_State *L, const gsl_odeiv_step_type *type, - size_t dim, double eps_abs, double eps_rel, - struct solver_type *st, bool use_jacob); - -extern const gsl_odeiv_step_type * -method_lookup (const char *method, const gsl_odeiv_step_type *default_type, - bool *needs_jacobian); - -#endif diff --git a/ode_source.c b/ode_source.c deleted file mode 100644 index 2b1012b0..00000000 --- a/ode_source.c +++ /dev/null @@ -1,303 +0,0 @@ - -/* 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. - */ - -void -FUNCTION (linear, array_copy) (TYPE (gsl_matrix) *dest, const double *src, size_t n) -{ - FUNCTION (gsl_matrix, const_view) srcv = FUNCTION (gsl_matrix, const_view_array) (src, n, 1); - FUNCTION (gsl_matrix, memcpy) (dest, &srcv.matrix); -} - -void -FUNCTION (ode, fenv_setup) (lua_State *L, int n) -{ - lua_newtable (L); - - FUNCTION (matrix, push_raw) (L, n, 1); - lua_rawseti (L, -2, FENV_Y); - - FUNCTION (matrix, push_raw) (L, n, 1); - lua_rawseti (L, -2, FENV_Y_BUFFER); - - FUNCTION (matrix, push_view) (L, NULL); - lua_rawseti (L, -2, FENV_DYDT); - - FUNCTION (matrix, push_view) (L, NULL); - lua_rawseti (L, -2, FENV_DFDY); - - FUNCTION (matrix, push_view) (L, NULL); - lua_rawseti (L, -2, FENV_DFDT); - - /* if multiplicity is 1 this is not really needed */ - FUNCTION (matrix, push_raw) (L, n, n); - lua_rawseti (L, -2, FENV_DFDY_BUFFER); - - lua_setfenv (L, -2); -} - -static int -FUNCTION (ode, hook_f) (double t, const double y[], double f[], void *params) -{ - struct params *p = params; - TYPE (gsl_matrix) *ym; - lua_State *L = p->L; - size_t n = p->n; - - /* push the ode system f function */ - lua_rawgeti (L, 2, FENV_F); - - lua_pushnumber (L, t); - - lua_rawgeti (L, 2, FENV_Y_BUFFER); - ym = FUNCTION (matrix, check) (L, -1); - FUNCTION (linear, array_copy) (ym, y, n); - - lua_rawgeti (L, 2, FENV_DYDT); - FUNCTION (matrix, set_view) (L, -1, f, 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) -{ - struct params *p = params; - lua_State *L = p->L; - size_t n = p->n; - TYPE (gsl_matrix) *ym; - double *jacob; - -#if MULTIPLICITY >= 2 - { - lua_rawgeti (L, 2, FENV_DFDY_BUFFER); - TYPE (gsl_matrix) *jm = FUNCTION (matrix, check) (L, -1); - jacob = jm->data; - lua_pop (L, 1); - } -#else - jacob = dfdy; -#endif - - /* push the ode system df function */ - lua_rawgeti (L, 2, FENV_DF); - - lua_pushnumber (L, t); - - lua_rawgeti (L, 2, FENV_Y_BUFFER); - ym = FUNCTION (matrix, check) (L, -1); - FUNCTION (linear, array_copy) (ym, y, n); - - lua_rawgeti (L, 2, FENV_DFDY); - FUNCTION (matrix, set_view) (L, -1, jacob, n, n, NULL); - - lua_rawgeti (L, 2, FENV_DFDT); - FUNCTION (matrix, set_view) (L, -1, dfdt, n, 1, NULL); - - lua_call (L, 4, 0); - -#if MULTIPLICITY == 2 - { - gsl_matrix_view dest = gsl_matrix_view_array (dfdy, 2*n, 2*n); - gsl_matrix_complex_view src = gsl_matrix_complex_view_array (jacob, n, n); - matrix_jacob_copy_cauchy_riemann (&dest.matrix, &src.matrix, n); - } -#endif - - return GSL_SUCCESS; -} - -int -FUNCTION (ode, new) (lua_State *L) -{ - const gsl_odeiv_step_type * T; - struct solver *s; - double eps_abs, eps_rel; - const char *method; - bool needs_jacobian; - int n; - - luaL_checktype (L, 1, LUA_TTABLE); - - lua_getfield (L, 1, "n"); - if (! lua_isnumber (L, -1)) - luaL_error (L, "missing dimension 'n' of the ODE system"); - n = lua_tointeger (L, -1); - lua_pop (L, 1); - - 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, & needs_jacobian); - - s = ode_solver_push_new (L, T, MULTIPLICITY * n, eps_abs, eps_rel, - TYPE (ode_solver_type), needs_jacobian); - - FUNCTION (ode, fenv_setup) (L, n); - - lua_getfield (L, 1, "f"); - if (lua_type (L, -1) != LUA_TFUNCTION) - return luaL_error (L, "field \"f\" should be a function"); - - mlua_fenv_set (L, -2, FENV_F); - - if (needs_jacobian) - { - lua_getfield (L, 1, "df"); - if (lua_type (L, -1) != LUA_TFUNCTION) - return luaL_error (L, "field \"df\" should be a function (jacobian)"); - mlua_fenv_set (L, -2, FENV_DF); - } - - s->h = -1; - - s->params->L = L; - s->params->n = n; - - s->system->function = & FUNCTION (ode, hook_f); - s->system->jacobian = & FUNCTION (ode, hook_jacob); - s->system->dimension = s->dimension; - s->system->params = s->params; - - return 1; -} - -struct solver * -FUNCTION (ode, check) (lua_State *L, int index) -{ - if (lua_getmetatable(L, index)) - { - struct solver_type *tp = TYPE (ode_solver_type); - struct solver *s; - - lua_getfield(L, LUA_REGISTRYINDEX, tp->metatable_name); - if (! lua_rawequal(L, -1, -2)) - luaL_typerror (L, index, "ODE solver"); - lua_pop (L, 2); - - s = lua_touserdata (L, index); - - return s; - } - - lua_pop (L, 1); - return NULL; -} - -int -FUNCTION (ode, free) (lua_State *L) -{ - struct solver *s = FUNCTION (ode, check) (L, 1); - - gsl_odeiv_evolve_free (s->evol); - gsl_odeiv_control_free (s->ctrl); - gsl_odeiv_step_free (s->step); - - return 0; -} - -int -FUNCTION (ode, set) (lua_State *L) -{ - struct solver *s = FUNCTION (ode, check) (L, 1); - TYPE (gsl_matrix) *y, *yode; - - s->t = luaL_checknumber (L, 2); - y = FUNCTION (matrix, check) (L, 3); - s->h = luaL_checknumber (L, 4); - - mlua_fenv_get (L, 1, FENV_Y); - yode = FUNCTION (matrix, check) (L, -1); - FUNCTION (gsl_matrix, memcpy) (yode, y); - - return 0; -} - -int -FUNCTION (ode, evolve) (lua_State *L) -{ - struct solver *s = FUNCTION (ode, check) (L, 1); - TYPE (gsl_matrix) *y; - int status; - double t1; - - t1 = luaL_checknumber (L, 2); - s->h = luaL_optnumber (L, 3, s->h); - - lua_settop (L, 1); - lua_getfenv (L, 1); - - lua_rawgeti (L, 2, FENV_Y); - y = FUNCTION (matrix, check) (L, -1); - lua_pop (L, 1); - - status = gsl_odeiv_evolve_apply (s->evol, s->ctrl, s->step, - s->system, & s->t, t1, & s->h, y->data); - - if (status != GSL_SUCCESS) - luaL_error (L, "error in ODE evolve: %s", gsl_strerror (status)); - - if (isnan (s->evol->yerr[0])) - luaL_error (L, "failed to converge, step too big"); - - return 0; -} - -int -FUNCTION (ode, get_t) (lua_State *L) -{ - struct solver *s = FUNCTION (ode, check) (L, 1); - lua_pushnumber (L, s->t); - return 1; -} - -int -FUNCTION (ode, get_y) (lua_State *L) -{ - FUNCTION (ode, check) (L, 1); - mlua_fenv_get (L, 1, FENV_Y); - return 1; -} - -int -FUNCTION (ode, index) (lua_State *L) -{ - return mlua_index_with_properties (L, FUNCTION (ode, properties), - FUNCTION (ode, methods), false); -} - -void -FUNCTION (ode, register) (lua_State *L) -{ - struct solver_type *st = TYPE (ode_solver_type); - - luaL_newmetatable (L, st->metatable_name); - luaL_register (L, NULL, FUNCTION (ode, metatable)); - lua_pop (L, 1); - - luaL_register (L, NULL, FUNCTION (ode, functions)); -} - -#if MULTIPLICITY != 2 && MULTIPLICITY != 1 -#error MULTIPLICITY not supported -#endif diff --git a/qr_decomp.c b/qr_decomp.c deleted file mode 100644 index c13035f3..00000000 --- a/qr_decomp.c +++ /dev/null @@ -1,171 +0,0 @@ - -#include <lua.h> -#include <lauxlib.h> - -#include <gsl/gsl_linalg.h> - -#include "gs-types.h" -#include "matrix.h" -#include "lua-utils.h" - -typedef struct { - gsl_matrix *m; - gsl_vector *tau; -} qr_decomp; - -static int qr_decomp_new (lua_State *L); -static int qr_decomp_free (lua_State *L); -static int qr_decomp_index (lua_State *L); -static int qr_decomp_solve (lua_State *L); -static int qr_decomp_lssolve (lua_State *L); -static int qr_decomp_unpack (lua_State *L); - -const struct luaL_Reg qr_decomp_metatable[] = { - {"__gc", qr_decomp_free}, - {"__index", qr_decomp_index}, - {NULL, NULL} -}; - -static const struct luaL_Reg qr_decomp_methods[] = { - {"solve", qr_decomp_solve}, - {"lssolve", qr_decomp_lssolve}, - {"unpack", qr_decomp_unpack}, - {NULL, NULL} -}; - - -static const struct luaL_Reg qr_decomp_functions[] = { - {"QR", qr_decomp_new}, - {NULL, NULL} -}; - -int -qr_decomp_new (lua_State *L) -{ - int status; - gsl_matrix *a = matrix_check (L, 1); - size_t m = a->size1, n = a->size2; - qr_decomp *qr = gs_new_object (sizeof(qr_decomp), L, GS_QR_DECOMP); - size_t tn = (m <= n ? m : n); - - qr->m = gsl_matrix_alloc (m, n); - gsl_matrix_memcpy (qr->m, a); - - qr->tau = gsl_vector_alloc (tn); - if (qr->tau == NULL) - return luaL_error (L, "out of memory"); - - status = gsl_linalg_QR_decomp (qr->m, qr->tau); - - if (status != GSL_SUCCESS) - { - return luaL_error (L, "error during QR decomposition: %s", - gsl_strerror (status)); - } - - lua_newtable (L); - lua_setfenv (L, -2); - - return 1; -} - -int -qr_decomp_free (lua_State *L) -{ - qr_decomp *qr = gs_check_userdata (L, 1, GS_QR_DECOMP); - gsl_vector_free (qr->tau); - gsl_matrix_free (qr->m); - return 0; -} - -int -qr_decomp_unpack (lua_State *L) -{ - qr_decomp *qr = gs_check_userdata (L, 1, GS_QR_DECOMP); - size_t m = qr->m->size1, n = qr->m->size2; - gsl_matrix *q = matrix_push_raw (L, m, m); - gsl_matrix *r = matrix_push_raw (L, m, n); - - gsl_linalg_QR_unpack (qr->m, qr->tau, q, r); - - return 2; -} - -int -qr_decomp_solve (lua_State *L) -{ - qr_decomp *qr = gs_check_userdata (L, 1, GS_QR_DECOMP); - size_t m = qr->m->size1, n = qr->m->size2; - gsl_matrix *b = matrix_check (L, 2); - gsl_vector_view b1 = gsl_matrix_column (b, 0); - - if (m != n) - return luaL_error (L, "\"solve\" requires a square matrix"); - - if (b->size1 != n || b->size2 != 1) - return luaL_error (L, "matrix dimensions mismatch"); - - { - int status; - gsl_matrix *x = matrix_push_raw (L, n, 1); - gsl_vector_view x1 = gsl_matrix_column (x, 0); - status = gsl_linalg_QR_solve (qr->m, qr->tau, &b1.vector, &x1.vector); - if (status != GSL_SUCCESS) - { - return luaL_error (L, "cannot solve the system: %s", - gsl_strerror (status)); - } - } - - return 1; -} - -int -qr_decomp_lssolve (lua_State *L) -{ - qr_decomp *qr = gs_check_userdata (L, 1, GS_QR_DECOMP); - size_t m = qr->m->size1, n = qr->m->size2; - gsl_matrix *b = matrix_check (L, 2); - gsl_vector_view b1 = gsl_matrix_column (b, 0); - - if (m <= n) - return luaL_error (L, "too few rows"); - - if (b->size1 != m || b->size2 != 1) - return luaL_error (L, "matrix dimensions mismatch"); - - { - int status; - gsl_matrix *x = matrix_push_raw (L, n, 1); - gsl_matrix *r = matrix_push_raw (L, m, 1); - gsl_vector_view x1 = gsl_matrix_column (x, 0); - gsl_vector_view r1 = gsl_matrix_column (r, 0); - - status = gsl_linalg_QR_lssolve (qr->m, qr->tau, &b1.vector, &x1.vector, - &r1.vector); - - if (status != GSL_SUCCESS) - { - return luaL_error (L, "cannot solve the system: %s", - gsl_strerror (status)); - } - } - - return 2; -} - -int -qr_decomp_index (lua_State *L) -{ - return mlua_index_methods (L, qr_decomp_methods); -} - -void -qr_decomp_register (lua_State *L) -{ - luaL_newmetatable (L, GS_METATABLE(GS_QR_DECOMP)); - luaL_register (L, NULL, qr_decomp_metatable); - lua_pop (L, 1); - - luaL_register (L, NULL, qr_decomp_functions); -} diff --git a/qr_decomp.h b/qr_decomp.h deleted file mode 100644 index 79ea8a84..00000000 --- a/qr_decomp.h +++ /dev/null @@ -1,8 +0,0 @@ -#ifndef QR_DECOMP_H -#define QR_DECOMP_H - -#include <lua.h> - -extern void qr_decomp_register (lua_State *L); - -#endif @@ -26,7 +26,6 @@ #include <lauxlib.h> #include "defs.h" -#include "lcomplex.h" static gsl_mode_t gsl_mode_from_string (const char *s) diff --git a/sf_gener.c b/sf_gener.c index 3438c5b3..6e301988 100644 --- a/sf_gener.c +++ b/sf_gener.c @@ -22,7 +22,6 @@ GSH_SF_CUSTOM(debye) GSH_SF_CUSTOM(fermi_dirac) GSH_SF(D, dilog) -GSH_SF_CUSTOM(cdilog) GSH_SF(D, erf) @@ -56,6 +55,5 @@ GSH_SF_COMP(DD, hyperg, 0F1) GSH_SF_CUSTOM(hyperg1F1) GSH_SF_CUSTOM(hypergU) GSH_SF_COMP(DDDD, hyperg, 2F1) -GSH_SF_CUSTOM(hyperg2F1conj) GSH_SF_CUSTOM(zeta) diff --git a/sf_implement.h b/sf_implement.h index 65030f06..34e12824 100644 --- a/sf_implement.h +++ b/sf_implement.h @@ -313,22 +313,6 @@ int GSH_LUA_NAME(hypergU) (lua_State *L) return push_gsl_result (L, &res); } -int GSH_LUA_NAME(hyperg2F1conj) (lua_State *L) -{ - Complex a = luaL_checkcomplex(L, 1); - double c = luaL_checknumber(L, 2); - double x = luaL_checknumber (L, 3); - gsl_sf_result res; - int status; - - status = gsl_sf_hyperg_2F1_conj_e (creal(a), cimag(a), c, x, &res); - - if (status != GSL_SUCCESS) - return luaL_error (L, "hyperg2F1conj: %s", gsl_strerror (status)); - - return push_gsl_result (L, &res); -} - int GSH_LUA_NAME(laguerre) (lua_State *L) { int i = luaL_checkinteger (L, 1); @@ -517,20 +501,3 @@ int GSH_LUA_NAME(legendreQ) (lua_State *L) return push_gsl_result (L, &res); } - -int GSH_LUA_NAME(cdilog) (lua_State *L) -{ - Complex z = luaL_checkcomplex (L, 1); - gsl_sf_result rr, ri; - double zr = creal(z), zi = cimag(z); - double r = sqrt(zr*zr+zi*zi), th = atan2(cimag(z), creal(z)); - int status; - - status = gsl_sf_complex_dilog_e (r, th, &rr, &ri); - - if (status != GSL_SUCCESS) - return luaL_error (L, "cdilog: %s", gsl_strerror (status)); - - lua_pushcomplex (L, rr.val + I * ri.val); - return 1; -} diff --git a/template_matrix_off.h b/template_matrix_off.h deleted file mode 100644 index 1fddac16..00000000 --- a/template_matrix_off.h +++ /dev/null @@ -1,30 +0,0 @@ - -#undef BASE -#undef SHORT -#undef MULTIPLICITY -#undef LUA_SHORT -#undef LUA_SHORTM -#undef ONE -#undef BLAS_ID -#undef PREFIX -#undef GSPREFIX - -#undef FUNCTION -#undef CONST_FUNCTION -#undef TYPE -#undef VIEW -#undef CONST_VIEW -#undef GS_TYPE -#undef GS_TYPENAME - -#undef LUA_TYPE -#undef LUA_FUNCTION -#undef LUAL_FUNCTION -#undef BLAS_FUNCTION - -#ifdef value_retrieve -#undef value_retrieve -#undef value_assign -#endif - -#define Complex double _Complex diff --git a/template_matrix_on.h b/template_matrix_on.h deleted file mode 100644 index 26892559..00000000 --- a/template_matrix_on.h +++ /dev/null @@ -1,85 +0,0 @@ - -/* template_matrix_on.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. - */ - -#define lua_Complex double _Complex - -#if defined(BASE_GSL_COMPLEX) -#undef complex -#undef Complex -#define BASE gsl_complex -#define SHORT complex -#define MULTIPLICITY 2 -#define ONE {{1.0, 0.0}} -#define LUA_SHORT complex -#define LUA_SHORTM Complex -#define BLAS_ID z -#define PREFIX "c" -#define GSPREFIX C - -#elif defined(BASE_DOUBLE) -#define BASE double -#define SHORT -#define MULTIPLICITY 1 -#define ONE 1.0 -#define LUA_SHORT number -#define LUA_SHORTM Number -#define BLAS_ID d -#define PREFIX "" -#define GSPREFIX - -#define value_retrieve(x) (x) -#define value_assign(x) (x) - -#else -#error unknown BASE_ directive -#endif - -#define CONCAT2x(a,b) a ## _ ## b -#define CONCAT2(a,b) CONCAT2x(a,b) -#define CONCAT3x(a,b,c) a ## _ ## b ## _ ## c -#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) -#define CONST_FUNCTION(dir,name) CONCAT3(dir,const,name) -#define TYPE(dir) dir -#define VIEW(dir) CONCAT2(dir,view) -#define CONST_VIEW(dir) CONCAT3(dir,const,view) -#define GS_TYPE(tp) CONCAT2(GS,tp) -#define GS_TYPENAME(tp) GS_METATABLE(CONCAT2(GS,tp)) -#else -#define FUNCTION(a,c) CONCAT3(a,SHORT,c) -#define CONST_FUNCTION(a,c) CONCAT4(a,SHORT,const,c) -#define TYPE(dir) CONCAT2(dir,SHORT) -#define VIEW(dir) CONCAT3(dir,SHORT,view) -#define CONST_VIEW(dir) CONCAT4(dir,SHORT,const,view) -#define GS_TYPE(tp) CONCAT2(GS,MYCAT(GSPREFIX,tp)) -#define GS_TYPENAME(tp) GS_METATABLE(GS_TYPE(tp)) -#endif - -#define LUA_TYPE CONCAT2(lua,LUA_SHORTM) -#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/template_matrix_oper_off.h b/template_matrix_oper_off.h deleted file mode 100644 index 3d7a7f3f..00000000 --- a/template_matrix_oper_off.h +++ /dev/null @@ -1,13 +0,0 @@ - -#undef SCALAR_OP -#undef OP_NAME -#undef OPER -#undef OPER_ELEM -#undef OP_ELEM_DEF -#undef BASE_OPER - -#undef CONCAT2x -#undef CONCAT2 -#undef OPER_FUNCTION -#undef SCALAR_MAT_FUNCTION - diff --git a/template_matrix_oper_on.h b/template_matrix_oper_on.h deleted file mode 100644 index a105f0ec..00000000 --- a/template_matrix_oper_on.h +++ /dev/null @@ -1,61 +0,0 @@ - -/* template_matrix_oper_on.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. - */ - -#if defined(OPER_ADD) -#define OPER add -#define OPER_ELEM add -#define OP_ELEM_DEF 1 -#define SCALAR_OP add_constant -#define OP_NAME "addition" -#define BASE_OPER(a,b) ((a) + (b)) - -#elif defined(OPER_MUL) -#define OPER mul -#define OPER_ELEM mul_elements -#define OP_ELEM_DEF 0 -#define SCALAR_OP scale -#define OP_NAME "multiplication" -#define BASE_OPER(a,b) ((a) * (b)) - -#elif defined(OPER_SUB) -#define OPER sub -#define OPER_ELEM sub -#define OP_ELEM_DEF 1 -#define SCALAR_OP add_constant -#define OP_NAME "subtraction" -#define BASE_OPER(a,b) ((a) - (b)) - -#elif defined(OPER_DIV) -#define OPER div -#define OPER_ELEM div_elements -#define OP_ELEM_DEF 0 -#define SCALAR_OP scale -#define OP_NAME "division" -#define BASE_OPER(a,b) ((a) / (b)) - -#else -#error matrix operation directive unknown -#endif - -#define CONCAT2x(a,b) a ## _ ## b -#define CONCAT2(a,b) CONCAT2x(a,b) - -#define OPER_FUNCTION(base) CONCAT2(base,OPER_ELEM) -#define SCALAR_MAT_FUNCTION(base) CONCAT2(base,OPER) |