/*
eigen.c
Ruby/GSL: Ruby extension library for GSL (GNU Scientific Library)
(C) Copyright 2001-2006 by Yoshiki Tsunesada
Ruby/GSL is free software: you can redistribute it and/or modify it
under the terms of the GNU General Public License.
This library is distributed in the hope that it will be useful, but
WITHOUT ANY WARRANTY.
*/
#include "rb_gsl.h"
#include "rb_gsl_array.h"
#include "rb_gsl_eigen.h"
#include "rb_gsl_complex.h"
#include <gsl/gsl_complex.h>
#include <gsl/gsl_complex_math.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_eigen.h>
static VALUE cgsl_eigen_symm_workspace;
static VALUE cgsl_eigen_symmv_workspace;
static VALUE cgsl_eigen_herm_workspace;
static VALUE cgsl_eigen_hermv_workspace;
static VALUE cgsl_eigen_values;
static VALUE cgsl_eigen_vectors;
static VALUE cgsl_eigen_vector;
static VALUE cgsl_eigen_vector_complex;
static VALUE cgsl_eigen_herm_vectors;
#ifdef HAVE_EIGEN_FRANCIS
static VALUE cgsl_eigen_francis_workspace;
#endif
#ifdef GSL_1_9_LATER
static VALUE cgsl_eigen_nonsymm_workspace;
static VALUE cgsl_eigen_nonsymmv_workspace;
#endif
#ifdef HAVE_EIGEN_GEN
static VALUE cgensymm, mgensymm;
#endif
static VALUE rb_gsl_eigen_symm_alloc(VALUE klass, VALUE nn)
{
gsl_eigen_symm_workspace *w = NULL;
CHECK_FIXNUM(nn);
w = gsl_eigen_symm_alloc(FIX2INT(nn));
return Data_Wrap_Struct(klass, 0,
gsl_eigen_symm_free, w);
}
static VALUE rb_gsl_eigen_symmv_alloc(VALUE klass, VALUE nn)
{
gsl_eigen_symmv_workspace *w = NULL;
CHECK_FIXNUM(nn);
w = gsl_eigen_symmv_alloc(FIX2INT(nn));
return Data_Wrap_Struct(klass, 0, gsl_eigen_symmv_free, w);
}
static VALUE rb_gsl_eigen_herm_alloc(VALUE klass, VALUE nn)
{
gsl_eigen_herm_workspace *w = NULL;
CHECK_FIXNUM(nn);
w = gsl_eigen_herm_alloc(FIX2INT(nn));
return Data_Wrap_Struct(klass, 0,
gsl_eigen_herm_free, w);
}
static VALUE rb_gsl_eigen_hermv_alloc(VALUE klass, VALUE nn)
{
gsl_eigen_hermv_workspace *w = NULL;
CHECK_FIXNUM(nn);
w = gsl_eigen_hermv_alloc(FIX2INT(nn));
return Data_Wrap_Struct(klass, 0,
gsl_eigen_hermv_free, w);
}
#ifdef HAVE_NARRAY_H
static VALUE rb_gsl_eigen_symm_narray(int argc, VALUE *argv, VALUE obj);
#endif
static VALUE rb_gsl_eigen_symm(int argc, VALUE *argv, VALUE obj)
{
gsl_matrix *Atmp = NULL, *A = NULL;
gsl_eigen_symm_workspace *w = NULL;
gsl_vector *v = NULL;
int flagw = 0;
switch (TYPE(obj)) {
case T_MODULE:
case T_CLASS:
case T_OBJECT:
switch (argc) {
case 2:
#ifdef HAVE_NARRAY_H
if (NA_IsNArray(argv[0])) return rb_gsl_eigen_symm_narray(argc, argv, obj);
#endif
CHECK_MATRIX(argv[0]);
Data_Get_Struct(argv[0], gsl_matrix, Atmp);
if (CLASS_OF(argv[1]) != cgsl_eigen_symm_workspace)
rb_raise(rb_eTypeError,
"argv[1]: wrong argument type %s (Eigen::Symm::Workspace expected)",
rb_class2name(CLASS_OF(argv[1])));
Data_Get_Struct(argv[1], gsl_eigen_symm_workspace, w);
break;
case 1:
#ifdef HAVE_NARRAY_H
if (NA_IsNArray(argv[0])) return rb_gsl_eigen_symm_narray(argc, argv, obj);
#endif
CHECK_MATRIX(argv[0]);
Data_Get_Struct(argv[0], gsl_matrix, Atmp);
w = gsl_eigen_symm_alloc(Atmp->size1);
flagw = 1;
break;
default:
rb_raise(rb_eArgError, "wrong number of arguments (%d for 1 or 2)", argc);
break;
}
break;
default:
CHECK_MATRIX(obj);
Data_Get_Struct(obj, gsl_matrix, Atmp);
switch (argc) {
case 1:
if (CLASS_OF(argv[0]) != cgsl_eigen_symm_workspace)
rb_raise(rb_eTypeError,
"argv[0]: wrong argument type %s (Eigen::Symm::Workspace expected",
rb_class2name(CLASS_OF(argv[0])));
Data_Get_Struct(argv[0], gsl_eigen_symm_workspace, w);
break;
case 0:
w = gsl_eigen_symm_alloc(Atmp->size1);
flagw = 1;
break;
default:
rb_raise(rb_eArgError, "wrong number of arguments (%d for 0 or 1)", argc);
}
}
A = make_matrix_clone(Atmp);
v = gsl_vector_alloc(A->size1);
gsl_eigen_symm(A, v, w);
/* gsl_sort_vector(v);*/
gsl_matrix_free(A);
if (flagw == 1) gsl_eigen_symm_free(w);
return Data_Wrap_Struct(cgsl_eigen_values, 0, gsl_vector_free, v);
}
#ifdef HAVE_NARRAY_H
static VALUE rb_gsl_eigen_symm_narray(int argc, VALUE *argv, VALUE obj)
{
struct NARRAY *na;
VALUE nary;
gsl_matrix *A = NULL;
gsl_eigen_symm_workspace *w = NULL;
gsl_vector_view vv;
int shape[1];
int flagw = 0;
switch (argc) {
case 2:
if (!NA_IsNArray(argv[0]))
rb_raise(rb_eTypeError, "wrong argument type %s (NArray expected)",
rb_class2name(CLASS_OF(argv[0])));
GetNArray(argv[0], na);
if (na->rank < 2) rb_raise(rb_eRuntimeError, "rank >= 2 required");
if (na->shape[0] != na->shape[1])
rb_raise(rb_eRuntimeError, "square matrix required");
A = gsl_matrix_alloc(na->shape[1], na->shape[0]);
memcpy(A->data, (double*) na->ptr, sizeof(double)*A->size1*A->size2);
if (CLASS_OF(argv[1]) != cgsl_eigen_symm_workspace)
rb_raise(rb_eTypeError,
"argv[1]: wrong argument type %s (Eigen::Symm::Workspace expected",
rb_class2name(CLASS_OF(argv[1])));
Data_Get_Struct(argv[1], gsl_eigen_symm_workspace, w);
flagw = 0;
break;
case 1:
if (!NA_IsNArray(argv[0]))
rb_raise(rb_eTypeError, "wrong argument type %s (NArray expected)",
rb_class2name(CLASS_OF(argv[0])));
GetNArray(argv[0], na);
if (na->rank < 2) rb_raise(rb_eRuntimeError, "rank >= 2 required");
if (na->shape[0] != na->shape[1])
rb_raise(rb_eRuntimeError, "square matrix required");
A = gsl_matrix_alloc(na->shape[1], na->shape[0]);
memcpy(A->data, (double*) na->ptr, sizeof(double)*A->size1*A->size2);
w = gsl_eigen_symm_alloc(A->size1);
flagw = 1;
break;
default:
rb_raise(rb_eArgError, "matrix not given");
break;
}
shape[0] = A->size1;
nary = na_make_object(NA_DFLOAT, 1, shape, cNVector);
vv = gsl_vector_view_array(NA_PTR_TYPE(nary,double*), A->size1);
gsl_eigen_symm(A, &vv.vector, w);
/* gsl_sort_vector(v);*/
gsl_matrix_free(A);
if (flagw == 1) gsl_eigen_symm_free(w);
return nary;
}
#endif
#ifdef HAVE_NARRAY_H
static VALUE rb_gsl_eigen_symmv_narray(int argc, VALUE *argv, VALUE obj);
#endif
static VALUE rb_gsl_eigen_symmv(int argc, VALUE *argv, VALUE obj)
{
gsl_matrix *Atmp = NULL, *A = NULL, *em = NULL;
gsl_eigen_symmv_workspace *w = NULL;
gsl_vector *v = NULL;
int flagw = 0;
VALUE vval, vvec;
switch (TYPE(obj)) {
case T_MODULE:
case T_CLASS:
case T_OBJECT:
switch (argc) {
case 2:
#ifdef HAVE_NARRAY_H
if (NA_IsNArray(argv[0])) return rb_gsl_eigen_symmv_narray(argc, argv, obj);
#endif
CHECK_MATRIX(argv[0]);
Data_Get_Struct(argv[0], gsl_matrix, Atmp);
if (CLASS_OF(argv[1]) != cgsl_eigen_symmv_workspace)
rb_raise(rb_eTypeError,
"argv[1]: wrong argument type %s (Eigen::Symmv::Workspace expected)",
rb_class2name(CLASS_OF(argv[1])));
Data_Get_Struct(argv[1], gsl_eigen_symmv_workspace, w);
break;
case 1:
#ifdef HAVE_NARRAY_H
if (NA_IsNArray(argv[0])) return rb_gsl_eigen_symmv_narray(argc, argv, obj);
#endif
CHECK_MATRIX(argv[0]);
Data_Get_Struct(argv[0], gsl_matrix, Atmp);
w = gsl_eigen_symmv_alloc(Atmp->size1);
flagw = 1;
break;
default:
rb_raise(rb_eArgError, "wrong number of arguments (%d for 1 or 2)", argc);
}
break;
default:
CHECK_MATRIX(obj);
Data_Get_Struct(obj, gsl_matrix, Atmp);
switch (argc) {
case 1:
if (CLASS_OF(argv[0]) != cgsl_eigen_symmv_workspace)
rb_raise(rb_eTypeError,
"argv[0]: wrong argument type %s (Eigen::Symmv::Workspace expected)",
rb_class2name(CLASS_OF(argv[0])));
Data_Get_Struct(argv[0], gsl_eigen_symmv_workspace, w);
break;
case 0:
w = gsl_eigen_symmv_alloc(Atmp->size1);
flagw = 1;
break;
default:
rb_raise(rb_eArgError, "wrong number of arguments (%d for 0 or 1)", argc);
break;
}
}
A = make_matrix_clone(Atmp);
em = gsl_matrix_alloc(A->size1, A->size2);
v = gsl_vector_alloc(A->size1);
gsl_eigen_symmv(A, v, em, w);
/* gsl_eigen_symmv_sort(v, em, GSL_EIGEN_SORT_VAL_DESC);*/
gsl_matrix_free(A);
if (flagw == 1) gsl_eigen_symmv_free(w);
vval = Data_Wrap_Struct(cgsl_eigen_values, 0, gsl_vector_free, v);
vvec = Data_Wrap_Struct(cgsl_eigen_vectors, 0, gsl_matrix_free, em);
return rb_ary_new3(2, vval, vvec);
}
#ifdef HAVE_NARRAY_H
static VALUE rb_gsl_eigen_symmv_narray(int argc, VALUE *argv, VALUE obj)
{
struct NARRAY *na;
VALUE eval, evec;
gsl_matrix *A = NULL;
gsl_eigen_symmv_workspace *w = NULL;
gsl_matrix_view mv;
gsl_vector_view vv;
int shape1[1], shape2[2];
int flagw = 0;
switch (argc) {
case 2:
if (!NA_IsNArray(argv[0]))
rb_raise(rb_eTypeError, "wrong argument type %s (NArray expected)",
rb_class2name(CLASS_OF(argv[0])));
GetNArray(argv[0], na);
if (na->rank < 2) rb_raise(rb_eRuntimeError, "rank >= 2 required");
if (na->shape[0] != na->shape[1])
rb_raise(rb_eRuntimeError, "square matrix required");
A = gsl_matrix_alloc(na->shape[1], na->shape[0]);
memcpy(A->data, (double*) na->ptr, sizeof(double)*A->size1*A->size2);
if (CLASS_OF(argv[1]) != cgsl_eigen_symmv_workspace)
rb_raise(rb_eTypeError,
"argv[1]: wrong argument type %s (Eigen::Symm::Workspace expected",
rb_class2name(CLASS_OF(argv[1])));
Data_Get_Struct(argv[1], gsl_eigen_symmv_workspace, w);
flagw = 0;
break;
case 1:
if (!NA_IsNArray(argv[0]))
rb_raise(rb_eTypeError, "wrong argument type %s (NArray expected)",
rb_class2name(CLASS_OF(argv[0])));
GetNArray(argv[0], na);
if (na->rank < 2) rb_raise(rb_eRuntimeError, "rank >= 2 required");
if (na->shape[0] != na->shape[1])
rb_raise(rb_eRuntimeError, "square matrix required");
A = gsl_matrix_alloc(na->shape[1], na->shape[0]);
memcpy(A->data, (double*) na->ptr, sizeof(double)*A->size1*A->size2);
w = gsl_eigen_symmv_alloc(A->size1);
flagw = 1;
break;
default:
rb_raise(rb_eArgError, "matrix not given");
break;
}
shape1[0] = A->size1;
shape2[0] = A->size1;
shape2[1] = A->size1;
eval = na_make_object(NA_DFLOAT, 1, shape1, cNVector);
evec = na_make_object(NA_DFLOAT, 2, shape2, CLASS_OF(argv[0]));
vv = gsl_vector_view_array(NA_PTR_TYPE(eval,double*), A->size1);
mv = gsl_matrix_view_array(NA_PTR_TYPE(evec,double*), A->size1, A->size2);
gsl_eigen_symmv(A, &vv.vector, &mv.matrix, w);
/* gsl_sort_vector(v);*/
gsl_matrix_free(A);
if (flagw == 1) gsl_eigen_symmv_free(w);
return rb_ary_new3(2, eval, evec);
}
#endif
static VALUE rb_gsl_eigen_herm(int argc, VALUE *argv, VALUE obj)
{
gsl_matrix_complex *Atmp = NULL, *A = NULL;
gsl_eigen_herm_workspace *w = NULL;
gsl_vector *v = NULL;
int flagw = 0;
switch (TYPE(obj)) {
case T_MODULE:
case T_CLASS:
case T_OBJECT:
switch (argc) {
case 2:
CHECK_MATRIX_COMPLEX(argv[0]);
Data_Get_Struct(argv[0], gsl_matrix_complex, Atmp);
if (CLASS_OF(argv[1]) != cgsl_eigen_herm_workspace)
rb_raise(rb_eTypeError,
"argv[1]: wrong argument type %s (Eigen::Herm::Workspace expected)",
rb_class2name(CLASS_OF(argv[1])));
Data_Get_Struct(argv[1], gsl_eigen_herm_workspace, w);
break;
case 1:
CHECK_MATRIX_COMPLEX(argv[0]);
Data_Get_Struct(argv[0], gsl_matrix_complex, Atmp);
w = gsl_eigen_herm_alloc(Atmp->size1);
flagw = 1;
break;
default:
rb_raise(rb_eArgError, "wrong number of arguments (%d for 1 or 2)", argc);
}
break;
default:
CHECK_MATRIX_COMPLEX(obj);
Data_Get_Struct(obj, gsl_matrix_complex, Atmp);
switch (argc) {
case 1:
if (CLASS_OF(argv[0]) != cgsl_eigen_herm_workspace)
rb_raise(rb_eTypeError,
"argv[0]: wrong argument type %s (Eigen::Herm::Workspace expected)",
rb_class2name(CLASS_OF(argv[0])));
Data_Get_Struct(argv[0], gsl_eigen_herm_workspace, w);
break;
case 0:
w = gsl_eigen_herm_alloc(Atmp->size1);
flagw = 1;
break;
default:
rb_raise(rb_eArgError, "wrong number of arguments (%d for 0 or 1)", argc);
}
}
A = make_matrix_complex_clone(Atmp);
v = gsl_vector_alloc(A->size1);
gsl_eigen_herm(A, v, w);
/* gsl_sort_vector(v);*/
gsl_matrix_complex_free(A);
if (flagw == 1) gsl_eigen_herm_free(w);
return Data_Wrap_Struct(cgsl_eigen_values, 0, gsl_vector_free, v);
}
static VALUE rb_gsl_eigen_hermv(int argc, VALUE *argv, VALUE obj)
{
gsl_matrix_complex *Atmp = NULL, *A = NULL, *em = NULL;
gsl_eigen_hermv_workspace *w = NULL;
gsl_vector *v = NULL;
int flagw = 0;
VALUE vval, vvec;
switch (TYPE(obj)) {
case T_MODULE:
case T_CLASS:
case T_OBJECT:
switch (argc) {
case 2:
CHECK_MATRIX_COMPLEX(argv[0]);
Data_Get_Struct(argv[0], gsl_matrix_complex, Atmp);
if (CLASS_OF(argv[1]) != cgsl_eigen_hermv_workspace)
rb_raise(rb_eTypeError,
"argv[1]: wrong argument type %s (Eigen::Hermv::Workspace expected)",
rb_class2name(CLASS_OF(argv[1])));
Data_Get_Struct(argv[1], gsl_eigen_hermv_workspace, w);
break;
case 1:
CHECK_MATRIX_COMPLEX(argv[0]);
Data_Get_Struct(argv[0], gsl_matrix_complex, Atmp);
w = gsl_eigen_hermv_alloc(Atmp->size1);
flagw = 1;
break;
default:
rb_raise(rb_eArgError, "wrong number of arguments (%d for 1 or 2)", argc);
}
break;
default:
CHECK_MATRIX_COMPLEX(obj);
Data_Get_Struct(obj, gsl_matrix_complex, Atmp);
switch (argc) {
case 1:
if (CLASS_OF(argv[0]) != cgsl_eigen_hermv_workspace)
rb_raise(rb_eTypeError,
"argv[0]: wrong argument type %s (Eigen::Hermv::Workspace expected)",
rb_class2name(CLASS_OF(argv[0])));
Data_Get_Struct(argv[0], gsl_eigen_hermv_workspace, w);
break;
case 0:
w = gsl_eigen_hermv_alloc(Atmp->size1);
flagw = 1;
break;
default:
rb_raise(rb_eArgError, "wrong number of arguments (%d for 0 or 1)", argc);
}
}
A = make_matrix_complex_clone(Atmp);
em = gsl_matrix_complex_alloc(A->size1, A->size2);
v = gsl_vector_alloc(A->size1);
gsl_eigen_hermv(A, v, em, w);
/* gsl_eigen_hermv_sort(v, em, GSL_EIGEN_SORT_VAL_DESC);*/
gsl_matrix_complex_free(A);
if (flagw == 1) gsl_eigen_hermv_free(w);
vval = Data_Wrap_Struct(cgsl_eigen_values, 0, gsl_vector_free, v);
vvec = Data_Wrap_Struct(cgsl_eigen_herm_vectors, 0, gsl_matrix_complex_free, em);
return rb_ary_new3(2, vval, vvec);
}
static VALUE rb_gsl_eigen_vectors_unpack(VALUE obj)
{
gsl_matrix *m = NULL;
gsl_vector *v = NULL;
size_t i, j;
double val;
VALUE ary, tmp;
Data_Get_Struct(obj, gsl_matrix, m);
ary = rb_ary_new2(m->size1);
for (i = 0; i < m->size1; i++) {
v = gsl_vector_alloc(m->size2);
for (j = 0; j < m->size2; j++) {
val = gsl_matrix_get(m, j, i);
gsl_vector_set(v, j, val);
}
tmp = Data_Wrap_Struct(cgsl_eigen_vector, 0, gsl_vector_free, v);
rb_ary_store(ary, i, tmp);
}
return ary;
}
static VALUE rb_gsl_eigen_vectors_complex_unpack(VALUE obj)
{
gsl_matrix_complex *m = NULL;
gsl_vector_complex *v = NULL;
size_t i, j;
gsl_complex z;
VALUE ary, tmp;
Data_Get_Struct(obj, gsl_matrix_complex, m);
ary = rb_ary_new2(m->size1);
for (i = 0; i < m->size1; i++) {
v = gsl_vector_complex_alloc(m->size2);
for (j = 0; j < m->size2; j++) {
z= gsl_matrix_complex_get(m, j, i);
gsl_vector_complex_set(v, j, z);
}
tmp = Data_Wrap_Struct(cgsl_eigen_vector_complex, 0, gsl_vector_complex_free, v);
rb_ary_store(ary, i, tmp);
}
return ary;
}
static void rb_gsl_eigen_define_const(VALUE module)
{
rb_define_const(module, "SORT_VAL_ASC", INT2FIX(GSL_EIGEN_SORT_VAL_ASC));
rb_define_const(module, "SORT_VAL_DESC", INT2FIX(GSL_EIGEN_SORT_VAL_DESC));
rb_define_const(module, "SORT_ABS_ASC", INT2FIX(GSL_EIGEN_SORT_ABS_ASC));
rb_define_const(module, "SORT_ABS_DESC", INT2FIX(GSL_EIGEN_SORT_ABS_DESC));
rb_define_const(module, "VAL_ASC", INT2FIX(GSL_EIGEN_SORT_VAL_ASC));
rb_define_const(module, "VAL_DESC", INT2FIX(GSL_EIGEN_SORT_VAL_DESC));
rb_define_const(module, "ABS_ASC", INT2FIX(GSL_EIGEN_SORT_ABS_ASC));
rb_define_const(module, "ABS_DESC", INT2FIX(GSL_EIGEN_SORT_ABS_DESC));
}
static VALUE rb_gsl_eigen_symmv_sort(int argc, VALUE *argv, VALUE obj)
{
gsl_vector *v;
gsl_matrix *m;
gsl_eigen_sort_t type = GSL_EIGEN_SORT_VAL_DESC;
switch (argc) {
case 3:
CHECK_FIXNUM(argv[2]);
type = FIX2INT(argv[2]);
/* no break, do next */
case 2:
CHECK_VECTOR(argv[0]); CHECK_MATRIX(argv[1]);
break;
default:
rb_raise(rb_eArgError, "wrong number of arguments (%d for 2 or 3)", argc);
}
Data_Get_Struct(argv[0], gsl_vector, v);
Data_Get_Struct(argv[1], gsl_matrix, m);
return INT2FIX(gsl_eigen_symmv_sort(v, m, type));
}
static VALUE rb_gsl_eigen_hermv_sort(int argc, VALUE *argv, VALUE obj)
{
gsl_vector *v;
gsl_matrix_complex *m;
gsl_eigen_sort_t type = GSL_EIGEN_SORT_VAL_DESC;
switch (argc) {
case 3:
CHECK_FIXNUM(argv[2]);
type = FIX2INT(argv[2]);
/* no break, do next */
case 2:
CHECK_VECTOR(argv[0]);
CHECK_MATRIX_COMPLEX(argv[1]);
break;
default:
rb_raise(rb_eArgError, "wrong number of arguments (%d for 2 or 3)", argc);
}
Data_Get_Struct(argv[0], gsl_vector, v);
Data_Get_Struct(argv[1], gsl_matrix_complex, m);
return INT2FIX(gsl_eigen_hermv_sort(v, m, type));
}
#ifdef HAVE_EIGEN_FRANCIS
static VALUE rb_gsl_eigen_francis_alloc(VALUE klass)
{
gsl_eigen_francis_workspace *w = NULL;
w = gsl_eigen_francis_alloc();
return Data_Wrap_Struct(klass, 0, gsl_eigen_francis_free, w);
}
static VALUE rb_gsl_eigen_francis_T(int argc, VALUE *argv, VALUE obj)
{
gsl_eigen_francis_workspace *w = NULL;
int istart = 0;
if (CLASS_OF(obj) == cgsl_eigen_francis_workspace) {
Data_Get_Struct(obj, gsl_eigen_francis_workspace, w);
istart = 0;
} else {
if (argc != 1) rb_raise(rb_eArgError, "too few arguments (%d for 1)\n", argc);
Data_Get_Struct(argv[0], gsl_eigen_francis_workspace, w);
istart = 1;
}
gsl_eigen_francis_T(FIX2INT(argv[istart]), w);
return Qtrue;
}
#ifdef HAVE_NARRAY_H
static VALUE rb_gsl_eigen_francis_narray(int argc, VALUE *argv, VALUE obj);
#endif
static VALUE rb_gsl_eigen_francis(int argc, VALUE *argv, VALUE obj)
{
gsl_matrix *m, *mtmp;
gsl_vector_complex *v;
gsl_eigen_francis_workspace *w;
int vflag = 0, wflag = 0;
int istart = 0;
VALUE *argv2;
#ifdef HAVE_NARRAY_H
if (NA_IsNArray(obj)) return rb_gsl_eigen_francis_narray(argc, argv, obj);
if (argc >= 1 && NA_IsNArray(argv[0]))
return rb_gsl_eigen_francis_narray(argc, argv, obj);
#endif
if (MATRIX_P(obj)) {
Data_Get_Struct(obj, gsl_matrix, m);
argv2 = argv;
istart = 0;
} else {
if (argc < 1) rb_raise(rb_eArgError, "Wrong number of arguments.\n");
Data_Get_Struct(argv[0], gsl_matrix, m);
istart = 1;
argv2 = argv + 1;
}
switch (argc-istart) {
case 0:
v = gsl_vector_complex_alloc(m->size1);
w = gsl_eigen_francis_alloc();
vflag = 1;
wflag = 1;
break;
case 1:
if (CLASS_OF(argv2[0]) == cgsl_vector_complex) {
Data_Get_Struct(argv2[0], gsl_vector_complex, v);
w = gsl_eigen_francis_alloc();
wflag = 1;
} else if (CLASS_OF(argv2[0]) == cgsl_eigen_francis_workspace) {
v = gsl_vector_complex_alloc(m->size1);
vflag = 1;
Data_Get_Struct(argv2[0], gsl_eigen_francis_workspace, w);
} else {
rb_raise(rb_eArgError, "Wrong argument type.\n");
}
break;
case 2:
CHECK_VECTOR_COMPLEX(argv2[0]);
if (CLASS_OF(argv2[1]) != cgsl_eigen_francis_workspace) {
rb_raise(rb_eArgError, "argv[1] must be a GSL::Eigen::francis::Workspace.\n");
}
Data_Get_Struct(argv2[0], gsl_vector_complex, v);
Data_Get_Struct(argv2[1], gsl_eigen_francis_workspace, w);
break;
default:
rb_raise(rb_eArgError, "Wrong number of arguments (%d for 0-2).\n", argc);
}
mtmp = make_matrix_clone(m);
gsl_eigen_francis(mtmp, v, w);
gsl_matrix_free(mtmp);
if (wflag == 1) gsl_eigen_francis_free(w);
if (vflag == 1)
return Data_Wrap_Struct(cgsl_vector_complex, 0, gsl_vector_complex_free, v);
else
return argv2[0];
}
#ifdef HAVE_NARRAY_H
static VALUE rb_gsl_eigen_francis_narray(int argc, VALUE *argv, VALUE obj)
{
struct NARRAY *na;
VALUE nary;
gsl_matrix *A = NULL;
gsl_eigen_francis_workspace *w = NULL;
gsl_vector_complex_view vv;
int shape[1];
int flagw = 0;
switch (argc) {
case 2:
if (!NA_IsNArray(argv[0]))
rb_raise(rb_eTypeError, "wrong argument type %s (NArray expected)",
rb_class2name(CLASS_OF(argv[0])));
GetNArray(argv[0], na);
if (na->rank < 2) rb_raise(rb_eRuntimeError, "rank >= 2 required");
if (na->shape[0] != na->shape[1])
rb_raise(rb_eRuntimeError, "square matrix required");
A = gsl_matrix_alloc(na->shape[1], na->shape[0]);
memcpy(A->data, (double*) na->ptr, sizeof(double)*A->size1*A->size2);
if (CLASS_OF(argv[1]) != cgsl_eigen_francis_workspace)
rb_raise(rb_eTypeError,
"argv[1]: wrong argument type %s (Eigen::Symm::Workspace expected",
rb_class2name(CLASS_OF(argv[1])));
Data_Get_Struct(argv[1], gsl_eigen_francis_workspace, w);
flagw = 0;
break;
case 1:
if (!NA_IsNArray(argv[0]))
rb_raise(rb_eTypeError, "wrong argument type %s (NArray expected)",
rb_class2name(CLASS_OF(argv[0])));
GetNArray(argv[0], na);
if (na->rank < 2) rb_raise(rb_eRuntimeError, "rank >= 2 required");
if (na->shape[0] != na->shape[1])
rb_raise(rb_eRuntimeError, "square matrix required");
A = gsl_matrix_alloc(na->shape[1], na->shape[0]);
memcpy(A->data, (double*) na->ptr, sizeof(double)*A->size1*A->size2);
w = gsl_eigen_francis_alloc();
flagw = 1;
break;
default:
rb_raise(rb_eArgError, "matrix not given");
break;
}
shape[0] = A->size1;
nary = na_make_object(NA_DCOMPLEX, 1, shape, cNVector);
vv = gsl_vector_complex_view_array(NA_PTR_TYPE(nary,double*), A->size1);
gsl_eigen_francis(A, &vv.vector, w);
/* gsl_sort_vector(v);*/
gsl_matrix_free(A);
if (flagw == 1) gsl_eigen_francis_free(w);
return nary;
}
#endif
static VALUE rb_gsl_eigen_francis_Z(int argc, VALUE *argv, VALUE obj)
{
gsl_matrix *m, *mtmp, *Z;
gsl_vector_complex *v;
gsl_eigen_francis_workspace *w;
int vflag = 0, wflag = 0;
int istart = 0;
VALUE *argv2, vv, ZZ;
if (MATRIX_P(obj)) {
Data_Get_Struct(obj, gsl_matrix, m);
argv2 = argv;
istart = 0;
} else {
if (argc < 1) rb_raise(rb_eArgError, "Wrong number of arguments.\n");
Data_Get_Struct(argv[0], gsl_matrix, m);
istart = 1;
argv2 = argv + 1;
}
switch (argc-istart) {
case 0:
v = gsl_vector_complex_alloc(m->size1);
Z = gsl_matrix_alloc(m->size1, m->size2);
w = gsl_eigen_francis_alloc();
vflag = 1;
wflag = 1;
break;
case 1:
if (CLASS_OF(argv2[0]) == cgsl_eigen_francis_workspace) {
v = gsl_vector_complex_alloc(m->size1);
Z = gsl_matrix_alloc(m->size1, m->size2);
vflag = 1;
Data_Get_Struct(argv2[0], gsl_eigen_francis_workspace, w);
} else {
rb_raise(rb_eArgError, "Wrong argument type.\n");
}
break;
case 3:
CHECK_VECTOR_COMPLEX(argv2[0]);
CHECK_MATRIX(argv2[1]);
if (CLASS_OF(argv2[2]) != cgsl_eigen_francis_workspace) {
rb_raise(rb_eArgError, "argv[1] must be a GSL::Eigen::francis::Workspace.\n");
}
Data_Get_Struct(argv2[0], gsl_vector_complex, v);
Data_Get_Struct(argv2[1], gsl_matrix, Z);
Data_Get_Struct(argv2[2], gsl_eigen_francis_workspace, w);
break;
default:
rb_raise(rb_eArgError, "Wrong number of arguments (%d for 0-2).\n", argc);
}
mtmp = make_matrix_clone(m);
gsl_eigen_francis_Z(mtmp, v, Z, w);
gsl_matrix_free(mtmp);
if (wflag == 1) gsl_eigen_francis_free(w);
if (vflag == 1) {
vv = Data_Wrap_Struct(cgsl_vector_complex, 0, gsl_vector_complex_free, v);
ZZ = Data_Wrap_Struct(cgsl_matrix, 0, gsl_matrix_free, Z);
} else {
vv = argv2[0];
ZZ = argv2[1];
}
return rb_ary_new3(2, vv, ZZ);
}
#endif
/* hogehoge */
#ifdef GSL_1_9_LATER
static VALUE rb_gsl_eigen_nonsymm_alloc(VALUE klass, VALUE nn)
{
size_t n;
gsl_eigen_nonsymm_workspace *w = NULL;
n = (size_t) FIX2UINT(nn);
w = gsl_eigen_nonsymm_alloc(n);
return Data_Wrap_Struct(klass, 0, gsl_eigen_nonsymm_free, w);
}
static VALUE rb_gsl_eigen_nonsymm_params(int argc, VALUE *argv, VALUE obj)
{
gsl_eigen_nonsymm_workspace *w = NULL;
int istart = 0;
if (CLASS_OF(obj) == cgsl_eigen_nonsymm_workspace) {
Data_Get_Struct(obj, gsl_eigen_nonsymm_workspace, w);
istart = 0;
} else {
if (argc != 3) rb_raise(rb_eArgError, "too few arguments (%d for 3)\n", argc);
Data_Get_Struct(argv[2], gsl_eigen_nonsymm_workspace, w);
istart = 1;
}
switch (argc - istart) {
case 2:
gsl_eigen_nonsymm_params(FIX2INT(argv[0]), FIX2INT(argv[1]), w);
break;
default:
rb_raise(rb_eArgError, "Wrong number of arguments.\n");
}
return Qtrue;
}
#ifdef HAVE_NARRAY_H
static VALUE rb_gsl_eigen_nonsymm_narray(int argc, VALUE *argv, VALUE obj);
#endif
static VALUE rb_gsl_eigen_nonsymm(int argc, VALUE *argv, VALUE obj)
{
gsl_matrix *m, *mtmp;
gsl_vector_complex *v;
gsl_eigen_nonsymm_workspace *w;
int vflag = 0, wflag = 0;
int istart = 0;
VALUE *argv2;
#ifdef HAVE_NARRAY_H
if (NA_IsNArray(obj)) return rb_gsl_eigen_nonsymm_narray(argc, argv, obj);
if (argc >= 1 && NA_IsNArray(argv[0]))
return rb_gsl_eigen_nonsymm_narray(argc, argv, obj);
#endif
if (MATRIX_P(obj)) {
Data_Get_Struct(obj, gsl_matrix, m);
argv2 = argv;
istart = 0;
} else {
if (argc < 1) rb_raise(rb_eArgError, "Wrong number of arguments.\n");
Data_Get_Struct(argv[0], gsl_matrix, m);
istart = 1;
argv2 = argv + 1;
}
switch (argc-istart) {
case 0:
v = gsl_vector_complex_alloc(m->size1);
w = gsl_eigen_nonsymm_alloc(m->size1);
vflag = 1;
wflag = 1;
break;
case 1:
if (CLASS_OF(argv2[0]) == cgsl_vector_complex) {
Data_Get_Struct(argv2[0], gsl_vector_complex, v);
w = gsl_eigen_nonsymm_alloc(m->size1);
wflag = 1;
} else if (CLASS_OF(argv2[0]) == cgsl_eigen_nonsymm_workspace) {
v = gsl_vector_complex_alloc(m->size1);
vflag = 1;
Data_Get_Struct(argv2[0], gsl_eigen_nonsymm_workspace, w);
} else {
rb_raise(rb_eArgError, "Wrong argument type.\n");
}
break;
case 2:
CHECK_VECTOR_COMPLEX(argv2[0]);
if (CLASS_OF(argv2[1]) != cgsl_eigen_nonsymm_workspace) {
rb_raise(rb_eArgError, "argv[1] must be a GSL::Eigen::Nonsymm::Workspace.\n");
}
Data_Get_Struct(argv2[0], gsl_vector_complex, v);
Data_Get_Struct(argv2[1], gsl_eigen_nonsymm_workspace, w);
break;
default:
rb_raise(rb_eArgError, "Wrong number of arguments (%d for 0-2).\n", argc);
}
mtmp = make_matrix_clone(m);
gsl_eigen_nonsymm(mtmp, v, w);
gsl_matrix_free(mtmp);
if (wflag == 1) gsl_eigen_nonsymm_free(w);
if (vflag == 1)
return Data_Wrap_Struct(cgsl_vector_complex, 0, gsl_vector_complex_free, v);
else
return argv2[0];
}
#ifdef HAVE_NARRAY_H
static VALUE rb_gsl_eigen_nonsymm_narray(int argc, VALUE *argv, VALUE obj)
{
struct NARRAY *na;
VALUE nary;
gsl_matrix *A = NULL;
gsl_eigen_nonsymm_workspace *w = NULL;
gsl_vector_complex_view vv;
int shape[1];
int flagw = 0;
switch (argc) {
case 2:
if (!NA_IsNArray(argv[0]))
rb_raise(rb_eTypeError, "wrong argument type %s (NArray expected)",
rb_class2name(CLASS_OF(argv[0])));
GetNArray(argv[0], na);
if (na->rank < 2) rb_raise(rb_eRuntimeError, "rank >= 2 required");
if (na->shape[0] != na->shape[1])
rb_raise(rb_eRuntimeError, "square matrix required");
A = gsl_matrix_alloc(na->shape[1], na->shape[0]);
memcpy(A->data, (double*) na->ptr, sizeof(double)*A->size1*A->size2);
if (CLASS_OF(argv[1]) != cgsl_eigen_nonsymm_workspace)
rb_raise(rb_eTypeError,
"argv[1]: wrong argument type %s (Eigen::Symm::Workspace expected",
rb_class2name(CLASS_OF(argv[1])));
Data_Get_Struct(argv[1], gsl_eigen_nonsymm_workspace, w);
flagw = 0;
break;
case 1:
if (!NA_IsNArray(argv[0]))
rb_raise(rb_eTypeError, "wrong argument type %s (NArray expected)",
rb_class2name(CLASS_OF(argv[0])));
GetNArray(argv[0], na);
if (na->rank < 2) rb_raise(rb_eRuntimeError, "rank >= 2 required");
if (na->shape[0] != na->shape[1])
rb_raise(rb_eRuntimeError, "square matrix required");
A = gsl_matrix_alloc(na->shape[1], na->shape[0]);
memcpy(A->data, (double*) na->ptr, sizeof(double)*A->size1*A->size2);
w = gsl_eigen_nonsymm_alloc(A->size1);
flagw = 1;
break;
default:
rb_raise(rb_eArgError, "matrix not given");
break;
}
shape[0] = A->size1;
nary = na_make_object(NA_DCOMPLEX, 1, shape, cNVector);
vv = gsl_vector_complex_view_array(NA_PTR_TYPE(nary,double*), A->size1);
gsl_eigen_nonsymm(A, &vv.vector, w);
/* gsl_sort_vector(v);*/
gsl_matrix_free(A);
if (flagw == 1) gsl_eigen_nonsymm_free(w);
return nary;
}
#endif
static VALUE rb_gsl_eigen_nonsymm_Z(int argc, VALUE *argv, VALUE obj)
{
gsl_matrix *m, *mtmp, *Z;
gsl_vector_complex *v;
gsl_eigen_nonsymm_workspace *w;
int vflag = 0, wflag = 0;
int istart = 0;
VALUE *argv2, vv, ZZ;
if (MATRIX_P(obj)) {
Data_Get_Struct(obj, gsl_matrix, m);
argv2 = argv;
istart = 0;
} else {
if (argc < 1) rb_raise(rb_eArgError, "Wrong number of arguments.\n");
Data_Get_Struct(argv[0], gsl_matrix, m);
istart = 1;
argv2 = argv + 1;
}
switch (argc-istart) {
case 0:
v = gsl_vector_complex_alloc(m->size1);
Z = gsl_matrix_alloc(m->size1, m->size2);
w = gsl_eigen_nonsymm_alloc(m->size1);
vflag = 1;
wflag = 1;
break;
case 1:
if (CLASS_OF(argv2[0]) == cgsl_eigen_nonsymm_workspace) {
v = gsl_vector_complex_alloc(m->size1);
Z = gsl_matrix_alloc(m->size1, m->size2);
vflag = 1;
Data_Get_Struct(argv2[0], gsl_eigen_nonsymm_workspace, w);
} else {
rb_raise(rb_eArgError, "Wrong argument type.\n");
}
break;
case 3:
CHECK_VECTOR_COMPLEX(argv2[0]);
CHECK_MATRIX(argv2[1]);
if (CLASS_OF(argv2[2]) != cgsl_eigen_nonsymm_workspace) {
rb_raise(rb_eArgError, "argv[1] must be a GSL::Eigen::Nonsymm::Workspace.\n");
}
Data_Get_Struct(argv2[0], gsl_vector_complex, v);
Data_Get_Struct(argv2[1], gsl_matrix, Z);
Data_Get_Struct(argv2[2], gsl_eigen_nonsymm_workspace, w);
break;
default:
rb_raise(rb_eArgError, "Wrong number of arguments (%d for 0-2).\n", argc);
}
mtmp = make_matrix_clone(m);
gsl_eigen_nonsymm_Z(mtmp, v, Z, w);
gsl_matrix_free(mtmp);
if (wflag == 1) gsl_eigen_nonsymm_free(w);
if (vflag == 1) {
vv = Data_Wrap_Struct(cgsl_vector_complex, 0, gsl_vector_complex_free, v);
ZZ = Data_Wrap_Struct(cgsl_matrix, 0, gsl_matrix_free, Z);
} else {
vv = argv2[0];
ZZ = argv2[1];
}
return rb_ary_new3(2, vv, ZZ);
}
static VALUE rb_gsl_eigen_nonsymmv_alloc(VALUE klass, VALUE nn)
{
size_t n;
gsl_eigen_nonsymmv_workspace *w = NULL;
n = (size_t) FIX2UINT(nn);
w = gsl_eigen_nonsymmv_alloc(n);
return Data_Wrap_Struct(klass, 0, gsl_eigen_nonsymmv_free, w);
}
#ifdef HAVE_NARRAY_H
static VALUE rb_gsl_eigen_nonsymmv_narray(int argc, VALUE *argv, VALUE obj);
#endif
static VALUE rb_gsl_eigen_nonsymmv(int argc, VALUE *argv, VALUE obj)
{
gsl_matrix *m, *mtmp;
gsl_vector_complex *v;
gsl_matrix_complex *evec;
gsl_eigen_nonsymmv_workspace *w;
int vflag = 0, wflag = 0;
int istart = 0;
VALUE *argv2;
#ifdef HAVE_NARRAY_H
if (NA_IsNArray(obj)) return rb_gsl_eigen_nonsymmv_narray(argc, argv, obj);
if (argc >= 1 && NA_IsNArray(argv[0]))
return rb_gsl_eigen_nonsymmv_narray(argc, argv, obj);
#endif
if (MATRIX_P(obj)) {
Data_Get_Struct(obj, gsl_matrix, m);
argv2 = argv;
istart = 0;
} else {
if (argc < 1) rb_raise(rb_eArgError, "Wrong number of arguments.\n");
Data_Get_Struct(argv[0], gsl_matrix, m);
istart = 1;
argv2 = argv + 1;
}
switch (argc-istart) {
case 0:
v = gsl_vector_complex_alloc(m->size1);
evec = gsl_matrix_complex_alloc(m->size1, m->size2);
w = gsl_eigen_nonsymmv_alloc(m->size1);
vflag = 1;
wflag = 1;
break;
case 1:
if (CLASS_OF(argv2[0]) == cgsl_eigen_nonsymm_workspace) {
v = gsl_vector_complex_alloc(m->size1);
evec = gsl_matrix_complex_alloc(m->size1, m->size2);
vflag = 1;
Data_Get_Struct(argv2[0], gsl_eigen_nonsymmv_workspace, w);
} else {
rb_raise(rb_eArgError, "Wrong argument type.\n");
}
break;
case 2:
CHECK_VECTOR_COMPLEX(argv2[0]);
CHECK_MATRIX_COMPLEX(argv2[1]);
w = gsl_eigen_nonsymmv_alloc(m->size1);
wflag = 1;
break;
case 3:
CHECK_VECTOR_COMPLEX(argv2[0]);
CHECK_MATRIX_COMPLEX(argv2[1]);
if (CLASS_OF(argv2[2]) != cgsl_eigen_nonsymm_workspace) {
rb_raise(rb_eArgError, "argv[1] must be a GSL::Eigen::Nonsymm::Workspace.\n");
}
Data_Get_Struct(argv2[0], gsl_vector_complex, v);
Data_Get_Struct(argv2[1], gsl_matrix_complex, evec);
Data_Get_Struct(argv2[2], gsl_eigen_nonsymmv_workspace, w);
break;
default:
rb_raise(rb_eArgError, "Wrong number of arguments (%d for 0-3).\n", argc);
}
mtmp = make_matrix_clone(m);
gsl_eigen_nonsymmv(mtmp, v, evec, w);
gsl_matrix_free(mtmp);
if (wflag == 1) gsl_eigen_nonsymmv_free(w);
if (vflag == 1) {
return rb_ary_new3(2,
Data_Wrap_Struct(cgsl_vector_complex, 0, gsl_vector_complex_free, v),
Data_Wrap_Struct(cgsl_matrix_complex, 0, gsl_matrix_complex_free, evec));
} else {
return rb_ary_new3(2, argv2[0], argv2[1]);
}
}
#ifdef HAVE_NARRAY_H
static VALUE rb_gsl_eigen_nonsymmv_narray(int argc, VALUE *argv, VALUE obj)
{
struct NARRAY *na;
VALUE nary, nvec;
gsl_matrix *A = NULL;
gsl_eigen_nonsymmv_workspace *w = NULL;
gsl_vector_complex_view vv;
gsl_matrix_complex_view mm;
int shape[1], shape2[2];
int flagw = 0;
switch (argc) {
case 2:
if (!NA_IsNArray(argv[0]))
rb_raise(rb_eTypeError, "wrong argument type %s (NArray expected)",
rb_class2name(CLASS_OF(argv[0])));
GetNArray(argv[0], na);
if (na->rank < 2) rb_raise(rb_eRuntimeError, "rank >= 2 required");
if (na->shape[0] != na->shape[1])
rb_raise(rb_eRuntimeError, "square matrix required");
A = gsl_matrix_alloc(na->shape[1], na->shape[0]);
memcpy(A->data, (double*) na->ptr, sizeof(double)*A->size1*A->size2);
if (CLASS_OF(argv[1]) != cgsl_eigen_nonsymmv_workspace)
rb_raise(rb_eTypeError,
"argv[1]: wrong argument type %s (Eigen::Symm::Workspace expected",
rb_class2name(CLASS_OF(argv[1])));
Data_Get_Struct(argv[1], gsl_eigen_nonsymmv_workspace, w);
flagw = 0;
break;
case 1:
if (!NA_IsNArray(argv[0]))
rb_raise(rb_eTypeError, "wrong argument type %s (NArray expected)",
rb_class2name(CLASS_OF(argv[0])));
GetNArray(argv[0], na);
if (na->rank < 2) rb_raise(rb_eRuntimeError, "rank >= 2 required");
if (na->shape[0] != na->shape[1])
rb_raise(rb_eRuntimeError, "square matrix required");
A = gsl_matrix_alloc(na->shape[1], na->shape[0]);
memcpy(A->data, (double*) na->ptr, sizeof(double)*A->size1*A->size2);
w = gsl_eigen_nonsymmv_alloc(A->size1);
flagw = 1;
break;
default:
rb_raise(rb_eArgError, "matrix not given");
break;
}
shape[0] = A->size1; shape2[0] = A->size1; shape2[1] = A->size2;
nary = na_make_object(NA_DCOMPLEX, 1, shape, cNVector);
vv = gsl_vector_complex_view_array(NA_PTR_TYPE(nary,double*), A->size1);
nvec = na_make_object(NA_DCOMPLEX, 2, shape2, CLASS_OF(argv[0]));
mm = gsl_matrix_complex_view_array(NA_PTR_TYPE(nvec,double*), A->size1, A->size2);
gsl_eigen_nonsymmv(A, &vv.vector, &mm.matrix, w);
/* gsl_sort_vector(v);*/
gsl_matrix_free(A);
if (flagw == 1) gsl_eigen_nonsymmv_free(w);
return rb_ary_new3(2, nary, nvec);
}
#endif
static VALUE rb_gsl_eigen_nonsymmv_Z(int argc, VALUE *argv, VALUE obj)
{
gsl_matrix *m, *mtmp, *Z;
gsl_vector_complex *v;
gsl_matrix_complex *evec;
gsl_eigen_nonsymmv_workspace *w;
int vflag = 0, wflag = 0;
int istart = 0;
VALUE *argv2;
if (MATRIX_P(obj)) {
Data_Get_Struct(obj, gsl_matrix, m);
argv2 = argv;
istart = 0;
} else {
if (argc < 1) rb_raise(rb_eArgError, "Wrong number of arguments.\n");
Data_Get_Struct(argv[0], gsl_matrix, m);
istart = 1;
argv2 = argv + 1;
}
switch (argc-istart) {
case 0:
v = gsl_vector_complex_alloc(m->size1);
evec = gsl_matrix_complex_alloc(m->size1, m->size2);
Z = gsl_matrix_alloc(m->size1, m->size2);
w = gsl_eigen_nonsymmv_alloc(m->size1);
vflag = 1;
wflag = 1;
break;
case 1:
if (CLASS_OF(argv2[0]) == cgsl_eigen_nonsymm_workspace) {
v = gsl_vector_complex_alloc(m->size1);
evec = gsl_matrix_complex_alloc(m->size1, m->size2);
vflag = 1;
Data_Get_Struct(argv2[0], gsl_eigen_nonsymmv_workspace, w);
} else {
rb_raise(rb_eArgError, "Wrong argument type.\n");
}
break;
case 3:
CHECK_VECTOR_COMPLEX(argv2[0]);
CHECK_MATRIX_COMPLEX(argv2[1]);
CHECK_MATRIX(argv2[2]);
w = gsl_eigen_nonsymmv_alloc(m->size1);
wflag = 1;
break;
case 4:
CHECK_VECTOR_COMPLEX(argv2[0]);
CHECK_MATRIX_COMPLEX(argv2[1]);
CHECK_MATRIX(argv2[2]);
if (CLASS_OF(argv2[3]) != cgsl_eigen_nonsymm_workspace) {
rb_raise(rb_eArgError, "argv[1] must be a GSL::Eigen::Nonsymm::Workspace.\n");
}
Data_Get_Struct(argv2[0], gsl_vector_complex, v);
Data_Get_Struct(argv2[1], gsl_matrix_complex, evec);
Data_Get_Struct(argv2[1], gsl_matrix, Z);
Data_Get_Struct(argv2[3], gsl_eigen_nonsymmv_workspace, w);
break;
default:
rb_raise(rb_eArgError, "Wrong number of arguments (%d for 0-3).\n", argc);
}
mtmp = make_matrix_clone(m);
gsl_eigen_nonsymmv_Z(mtmp, v, evec, Z, w);
gsl_matrix_free(mtmp);
if (wflag == 1) gsl_eigen_nonsymmv_free(w);
if (vflag == 1) {
return rb_ary_new3(3,
Data_Wrap_Struct(cgsl_vector_complex, 0, gsl_vector_complex_free, v),
Data_Wrap_Struct(cgsl_matrix_complex, 0, gsl_matrix_complex_free, evec),
Data_Wrap_Struct(cgsl_matrix, 0, gsl_matrix_free, Z));
} else {
return rb_ary_new3(2, argv2[0], argv2[1], argv2[2]);
}
}
static VALUE rb_gsl_eigen_nonsymmv_sort(int argc, VALUE *argv, VALUE obj)
{
gsl_vector_complex *v;
gsl_matrix_complex *m;
gsl_eigen_sort_t type = GSL_EIGEN_SORT_ABS_DESC;
switch (argc) {
case 3:
CHECK_FIXNUM(argv[2]);
type = FIX2INT(argv[2]);
/* no break, do next */
case 2:
CHECK_VECTOR_COMPLEX(argv[0]);
CHECK_MATRIX_COMPLEX(argv[1]);
break;
default:
rb_raise(rb_eArgError, "wrong number of arguments (%d for 2 or 3)", argc);
}
Data_Get_Struct(argv[0], gsl_vector_complex, v);
Data_Get_Struct(argv[1], gsl_matrix_complex, m);
return INT2FIX(gsl_eigen_nonsymmv_sort(v, m, type));
}
#endif
#ifdef HAVE_EIGEN_GEN
static VALUE rb_gsl_eigen_gensymm_alloc(VALUE klass, VALUE nn)
{
gsl_eigen_gen_workspace *w = NULL;
CHECK_FIXNUM(nn);
w = gsl_eigen_gen_alloc(FIX2INT(nn));
return Data_Wrap_Struct(cgensymm, 0, gsl_eigen_gen_free, w);
}
static int check_argv_gensymm(int argc, VALUE *argv, VALUE obj, gsl_matrix **A, gsl_matrix **B,
gsl_vector **eval, gsl_eigen_gensymm_workspace **w);
static VALUE rb_gsl_eigen_gensymm(int argc, VALUE *argv, VALUE obj)
{
gsl_matrix *A = NULL, *B = NULL;
gsl_vector *eval = NULL;
gsl_eigen_gensymm_workspace *w = NULL;
int flag;
VALUE veval;
flag = check_argv_gensymm(argc, argv, obj, &A, &B, &eval, &w);
// gsl_eigen_gensymm(A, B, eval, w);
switch (flag) {
case 0:
veval = argv[2];
break;
case 1:
veval = Data_Wrap_Struct(cgsl_eigen_values, 0, gsl_vector_free, eval);
break;
case 2:
veval = argv[2];
gsl_eigen_gensymm_free(w);
break;
case 3:
veval = Data_Wrap_Struct(cgsl_eigen_values, 0, gsl_vector_free, eval);
gsl_eigen_gensymm_free(w);
break;
}
return veval;
}
static int check_argv_gensymm(int argc, VALUE *argv, VALUE obj, gsl_matrix **A, gsl_matrix **B,
gsl_vector **eval, gsl_eigen_gensymm_workspace **w)
{
int argc2 = argc;
int flag = 0;
if (CLASS_OF(obj) == cgensymm) {
Data_Get_Struct(obj, gsl_eigen_gensymm_workspace, *w);
} else {
if (rb_obj_is_kind_of(argv[argc-1], cgensymm)) {
Data_Get_Struct(argv[argc-1], gsl_eigen_gensymm_workspace, *w);
argc2 = argc-1;
} else {
/* workspace is not given */
}
}
switch (argc2) {
case 2:
CHECK_MATRIX(argv[0]); CHECK_MATRIX(argv[1]);
Data_Get_Struct(argv[0], gsl_matrix, *A);
Data_Get_Struct(argv[1], gsl_matrix, *B);
break;
case 3:
if (rb_obj_is_kind_of(argv[2], cgensymm)) {
Data_Get_Struct(argv[2], gsl_eigen_gensymm_workspace, *w);
} else {
CHECK_VECTOR(argv[2]);
Data_Get_Struct(argv[2], gsl_vector, *eval);
}
CHECK_MATRIX(argv[0]); CHECK_MATRIX(argv[1]);
Data_Get_Struct(argv[0], gsl_matrix, *A);
Data_Get_Struct(argv[1], gsl_matrix, *B);
break;
default:
rb_raise(rb_eArgError, "wrong number of arguments (%d for 2 or 3)", argc);
}
if (*eval == NULL) {
*eval = gsl_vector_alloc((*A)->size1);
flag += 1;
}
if (*w == NULL) {
*w = gsl_eigen_gensymm_alloc((*A)->size1);
flag += 2;
}
return flag;
}
#endif
void Init_gsl_eigen(VALUE module)
{
VALUE mgsl_eigen;
VALUE mgsl_eigen_symm;
VALUE mgsl_eigen_symmv;
VALUE mgsl_eigen_herm;
VALUE mgsl_eigen_hermv;
#ifdef HAVE_EIGEN_FRANCIS
VALUE mgsl_eigen_francis;
#endif
#ifdef GSL_1_9_LATER
VALUE mgsl_eigen_nonsymmv;
VALUE mgsl_eigen_nonsymm;
#endif
mgsl_eigen = rb_define_module_under(module, "Eigen");
mgsl_eigen_symm = rb_define_module_under(mgsl_eigen, "Symm");
mgsl_eigen_symmv = rb_define_module_under(mgsl_eigen, "Symmv");
mgsl_eigen_herm = rb_define_module_under(mgsl_eigen, "Herm");
mgsl_eigen_hermv = rb_define_module_under(mgsl_eigen, "Hermv");
cgsl_eigen_values = rb_define_class_under(mgsl_eigen, "EigenValues",
cgsl_vector);
cgsl_eigen_vectors = rb_define_class_under(mgsl_eigen, "EigenVectors",
cgsl_matrix);
cgsl_eigen_vector = rb_define_class_under(mgsl_eigen, "EigenVector",
cgsl_vector);
cgsl_eigen_herm_vectors = rb_define_class_under(mgsl_eigen_herm, "EigenVectors",
cgsl_matrix_complex);
cgsl_eigen_vector_complex = rb_define_class_under(mgsl_eigen_herm, "EigenVector",
cgsl_vector_complex);
cgsl_eigen_symm_workspace = rb_define_class_under(mgsl_eigen_symm,
"Workspace", cGSL_Object);
cgsl_eigen_symmv_workspace = rb_define_class_under(mgsl_eigen_symmv,
"Workspace", cGSL_Object);
cgsl_eigen_herm_workspace = rb_define_class_under(mgsl_eigen_herm,
"Workspace", cGSL_Object);
cgsl_eigen_hermv_workspace = rb_define_class_under(mgsl_eigen_hermv,
"Workspace", cGSL_Object);
rb_define_singleton_method(cgsl_eigen_symm_workspace, "alloc",
rb_gsl_eigen_symm_alloc, 1);
rb_define_singleton_method(cgsl_eigen_symmv_workspace, "alloc",
rb_gsl_eigen_symmv_alloc, 1);
rb_define_singleton_method(cgsl_eigen_herm_workspace, "alloc",
rb_gsl_eigen_herm_alloc, 1);
rb_define_singleton_method(cgsl_eigen_hermv_workspace, "alloc",
rb_gsl_eigen_hermv_alloc, 1);
rb_define_module_function(mgsl_eigen, "symm",
rb_gsl_eigen_symm, -1);
rb_define_module_function(mgsl_eigen, "symmv",
rb_gsl_eigen_symmv, -1);
rb_define_module_function(mgsl_eigen, "herm",
rb_gsl_eigen_herm, -1);
rb_define_module_function(mgsl_eigen, "herm",
rb_gsl_eigen_hermv, -1);
rb_define_module_function(module, "eigen_symm",
rb_gsl_eigen_symm, -1);
rb_define_module_function(module, "eigen_symmv",
rb_gsl_eigen_symmv, -1);
rb_define_module_function(module, "eigen_herm",
rb_gsl_eigen_herm, -1);
rb_define_module_function(module, "eigen_hermv",
rb_gsl_eigen_hermv, -1);
rb_define_method(cgsl_matrix, "eigen_symm", rb_gsl_eigen_symm, -1);
rb_define_method(cgsl_matrix, "eigen_symmv", rb_gsl_eigen_symmv, -1);
rb_define_method(cgsl_matrix_complex, "eigen_herm", rb_gsl_eigen_herm, -1);
rb_define_method(cgsl_matrix_complex, "eigen_hermv", rb_gsl_eigen_hermv, -1);
rb_define_method(cgsl_eigen_vectors, "unpack", rb_gsl_eigen_vectors_unpack, 0);
rb_define_method(cgsl_eigen_herm_vectors, "unpack", rb_gsl_eigen_vectors_complex_unpack, 0);
rb_gsl_eigen_define_const(mgsl_eigen);
rb_define_module_function(mgsl_eigen, "symmv_sort",
rb_gsl_eigen_symmv_sort, -1);
rb_define_module_function(mgsl_eigen, "hermv_sort",
rb_gsl_eigen_hermv_sort, -1);
rb_define_module_function(mgsl_eigen_symmv, "sort",
rb_gsl_eigen_symmv_sort, -1);
rb_define_module_function(mgsl_eigen_hermv, "sort",
rb_gsl_eigen_hermv_sort, -1);
#ifdef HAVE_EIGEN_FRANCIS
mgsl_eigen_francis = rb_define_module_under(mgsl_eigen, "francis");
cgsl_eigen_francis_workspace = rb_define_class_under(mgsl_eigen_francis,
"Workspace", cGSL_Object);
rb_define_singleton_method(cgsl_eigen_francis_workspace, "alloc",
rb_gsl_eigen_francis_alloc, 0);
rb_define_method(cgsl_matrix, "eigen_francis", rb_gsl_eigen_francis, -1);
rb_define_module_function(mgsl_eigen, "francis", rb_gsl_eigen_francis, -1);
rb_define_module_function(module, "eigen_francis", rb_gsl_eigen_francis, -1);
rb_define_method(cgsl_matrix, "eigen_francis_Z", rb_gsl_eigen_francis_Z, -1);
rb_define_module_function(mgsl_eigen, "francis_Z", rb_gsl_eigen_francis_Z, -1);
rb_define_module_function(module, "eigen_francis_Z", rb_gsl_eigen_francis_Z, -1);
rb_define_method(cgsl_eigen_francis_workspace, "T", rb_gsl_eigen_francis_T, -1);
rb_define_module_function(mgsl_eigen_francis, "T", rb_gsl_eigen_francis_T, -1);
#endif
#ifdef GSL_1_9_LATER
mgsl_eigen_nonsymm = rb_define_module_under(mgsl_eigen, "Nonsymm");
mgsl_eigen_nonsymmv = rb_define_module_under(mgsl_eigen, "Nonsymmv");
cgsl_eigen_nonsymm_workspace = rb_define_class_under(mgsl_eigen_nonsymm,
"Workspace", cGSL_Object);
rb_define_singleton_method(cgsl_eigen_nonsymm_workspace, "alloc",
rb_gsl_eigen_nonsymm_alloc, 1);
rb_define_method(cgsl_matrix, "eigen_nonsymm", rb_gsl_eigen_nonsymm, -1);
rb_define_module_function(mgsl_eigen, "nonsymm", rb_gsl_eigen_nonsymm, -1);
rb_define_module_function(module, "eigen_nonsymm", rb_gsl_eigen_nonsymm, -1);
rb_define_method(cgsl_matrix, "eigen_nonsymm_Z", rb_gsl_eigen_nonsymm_Z, -1);
rb_define_module_function(mgsl_eigen, "nonsymm_Z", rb_gsl_eigen_nonsymm_Z, -1);
rb_define_module_function(module, "eigen_nonsymm_Z", rb_gsl_eigen_nonsymm_Z, -1);
rb_define_method(cgsl_eigen_nonsymm_workspace, "params", rb_gsl_eigen_nonsymm_params, -1);
rb_define_module_function(mgsl_eigen_nonsymm, "params", rb_gsl_eigen_nonsymm_params, -1);
cgsl_eigen_nonsymmv_workspace = rb_define_class_under(mgsl_eigen_nonsymmv,
"Workspace", cGSL_Object);
rb_define_singleton_method(cgsl_eigen_nonsymmv_workspace, "alloc",
rb_gsl_eigen_nonsymmv_alloc, 1);
rb_define_method(cgsl_matrix, "eigen_nonsymmv", rb_gsl_eigen_nonsymmv, -1);
rb_define_module_function(mgsl_eigen, "nonsymmv", rb_gsl_eigen_nonsymmv, -1);
rb_define_module_function(module, "eigen_nonsymmv", rb_gsl_eigen_nonsymmv, -1);
rb_define_method(cgsl_matrix, "eigen", rb_gsl_eigen_nonsymmv, -1);
rb_define_alias(cgsl_matrix, "eig", "eigen");
rb_define_method(cgsl_matrix, "eigen_nonsymmv_Z", rb_gsl_eigen_nonsymmv_Z, -1);
rb_define_module_function(mgsl_eigen, "nonsymmv_Z", rb_gsl_eigen_nonsymmv_Z, -1);
rb_define_module_function(module, "eigen_nonsymmv_Z", rb_gsl_eigen_nonsymmv_Z, -1);
rb_define_module_function(mgsl_eigen, "nonsymmv_sort",
rb_gsl_eigen_nonsymmv_sort, -1);
rb_define_module_function(mgsl_eigen_nonsymmv, "sort",
rb_gsl_eigen_nonsymmv_sort, -1);
rb_define_module_function(module, "eigen_nonsymmv_sort",
rb_gsl_eigen_nonsymmv_sort, -1);
#endif
#ifdef HAVE_EIGEN_GEN
mgensymm = rb_define_module_under(mgsl_eigen, "Gensymm");
cgensymm = rb_define_class_under(mgensymm, "Workspace", cGSL_Object);
rb_define_singleton_method(mgensymm, "alloc", rb_gsl_eigen_gensymm_alloc, 1);
rb_define_method(cgensymm, "gensymm", rb_gsl_eigen_gensymm, -1);
rb_define_module_function(module, "eigen_gensymm", rb_gsl_eigen_gensymm, -1);
rb_define_module_function(mgsl_eigen, "gensymm", rb_gsl_eigen_gensymm, -1);
rb_define_module_function(mgensymm, "gensymm", rb_gsl_eigen_gensymm, -1);
#endif
}
syntax highlighted by Code2HTML, v. 0.9.1