From bacae787abd380690126e6adf82ab0cd19f5abec Mon Sep 17 00:00:00 2001 From: Franz Zotter Date: Sat, 10 Jan 2009 18:33:39 +0000 Subject: added [mtx_eig] to compute eigenvalues and eigenvectors using GSL. Something seems to be wrong with the eigenvectors: result differs from GNU octave result. svn path=/trunk/externals/iem/iemmatrix/; revision=10502 --- src/mtx_eig.c | 227 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 227 insertions(+) create mode 100644 src/mtx_eig.c (limited to 'src/mtx_eig.c') diff --git a/src/mtx_eig.c b/src/mtx_eig.c new file mode 100644 index 0000000..04dc068 --- /dev/null +++ b/src/mtx_eig.c @@ -0,0 +1,227 @@ +/* + * iemmatrix + * + * objects for manipulating simple matrices + * mostly refering to matlab/octave matrix functions + * this functions depends on the GNU scientific library + * + * Copyright (c) 2009, Franz Zotter + * IEM, Graz, Austria + * + * For information on usage and redistribution, and for a DISCLAIMER OF ALL + * WARRANTIES, see the file, "LICENSE.txt," in this distribution. + * + */ + +#include "iemmatrix.h" +#include + +#ifdef HAVE_LIBGSL +#include +#endif + +static t_class *mtx_eig_class; +enum WithEigenVectors {WITHEVS=1, WITHOUTEVS=0}; +typedef struct _MTXEig_ MTXEig; +struct _MTXEig_ +{ + t_object x_obj; +#ifdef HAVE_LIBGSL + gsl_matrix *a; + gsl_matrix_complex *q; + gsl_vector_complex *l; + gsl_eigen_nonsymm_workspace *w; + gsl_eigen_nonsymmv_workspace *wv; +#endif + t_outlet *list_q_out_re; + t_outlet *list_q_out_im; + t_outlet *list_l_out_re; + t_outlet *list_l_out_im; + t_atom *list_q_re; + t_atom *list_q_im; + t_atom *list_l_re; + t_atom *list_l_im; + int size; + enum WithEigenVectors withevs; +}; + +#ifdef HAVE_LIBGSL +static void allocMTXqlw (MTXEig *x) +{ + x->a=(gsl_matrix*)gsl_matrix_alloc(x->size,x->size); + x->l=(gsl_vector_complex*)gsl_vector_complex_alloc(x->size); + + switch (x->withevs) { + case WITHEVS: + x->wv=(gsl_eigen_nonsymmv_workspace*)gsl_eigen_nonsymmv_alloc(x->size); + x->q=(gsl_matrix_complex*)gsl_matrix_complex_alloc(x->size,x->size); + break; + case WITHOUTEVS: + x->w=(gsl_eigen_nonsymm_workspace*)gsl_eigen_nonsymm_alloc(x->size); + } + + x->list_q_re=(t_atom*)calloc(sizeof(t_atom),x->size*x->size+2); + x->list_q_im=(t_atom*)calloc(sizeof(t_atom),x->size*x->size+2); + x->list_l_re=(t_atom*)calloc(sizeof(t_atom),x->size); + x->list_l_im=(t_atom*)calloc(sizeof(t_atom),x->size); +} + +static void deleteMTXqlw (MTXEig *x) +{ + if (x->list_q_re!=0) + free(x->list_q_re); + if (x->list_q_im!=0) + free(x->list_q_im); + if (x->list_l_re!=0) + free(x->list_l_re); + if (x->list_l_im!=0) + free(x->list_l_im); + + x->list_q_re = 0; + x->list_q_im = 0; + x->list_l_re = 0; + x->list_l_im = 0; + + if (x->a!=0) + gsl_matrix_free(x->a); + if (x->q!=0) + gsl_matrix_complex_free(x->q); + if (x->l!=0) + gsl_vector_complex_free(x->l); + if (x->w!=0) + gsl_eigen_nonsymm_free(x->w); + if (x->wv!=0) + gsl_eigen_nonsymmv_free(x->wv); + + x->a = 0; + x->q = 0; + x->l = 0; + x->w = 0; + x->wv = 0; +} +#endif + +static void deleteMTXEig (MTXEig *x) +{ +#ifdef HAVE_LIBGSL + deleteMTXqlw(x); +#endif +} + +static void *newMTXEig (t_symbol *s, int argc, t_atom *argv) +{ + MTXEig *x = (MTXEig *) pd_new (mtx_eig_class); + + x->list_l_out_re = outlet_new (&x->x_obj, gensym("list")); + x->list_l_out_im = outlet_new (&x->x_obj, gensym("list")); + if (atom_getsymbol(argv)==gensym("v")) { + x->withevs=1; + x->list_q_out_re = outlet_new (&x->x_obj, gensym("matrix")); + x->list_q_out_im = outlet_new (&x->x_obj, gensym("matrix")); + } + + x->list_l_re = 0; + x->list_l_im = 0; + x->list_q_re = 0; + x->list_q_im = 0; +#ifdef HAVE_LIBGSL + x->a=0; + x->q=0; + x->l=0; + x->w=0; + x->wv=0; +#endif + + return ((void *) x); +} + +static void mTXEigBang (MTXEig *x) +{ + if (x->list_l_re) { + switch (x->withevs) { + case WITHEVS: + outlet_anything(x->list_q_out_im, gensym("matrix"), x->size*x->size+2, x->list_q_im); + outlet_anything(x->list_q_out_re, gensym("matrix"), x->size*x->size+2, x->list_q_re); + case WITHOUTEVS: + outlet_anything(x->list_l_out_im, gensym("list"), x->size, x->list_l_im); + outlet_anything(x->list_l_out_re, gensym("list"), x->size, x->list_l_re); + } + } +} + +static void mTXEigMatrix (MTXEig *x, t_symbol *s, + int argc, t_atom *argv) +{ + int rows = atom_getint (argv++); + int columns = atom_getint (argv++); + int size = rows * columns; + int in_size = argc-2; + int n,m; + float f; + +#ifdef HAVE_LIBGSL + gsl_complex c; + /* size check */ + if (!size) + post("mtx_eig: invalid dimensions"); + else if (in_sizesize=size; + + deleteMTXqlw(x); + allocMTXqlw(x); + + for (n=0;na->data[n]=(double) atom_getfloat(argv++); + + switch (x->withevs) { + case WITHOUTEVS: + gsl_eigen_nonsymm(x->a,x->l,x->w); + break; + case WITHEVS: + gsl_eigen_nonsymmv(x->a,x->l,x->q,x->wv); + SETFLOAT((x->list_q_re),(float) x->size); + SETFLOAT((x->list_q_im),(float) x->size); + SETFLOAT((x->list_q_re+1),(float) x->size); + SETFLOAT((x->list_q_im+1),(float) x->size); + for (n=0;nlist_q_im+2+n), (float) x->q->data[2*n+1]); + SETFLOAT((x->list_q_re+2+n), (float) x->q->data[2*n]); + } + break; + } + + for (n=0;nsize;n++) { + f=(float) GSL_VECTOR_IMAG(x->l, n); + SETFLOAT((x->list_l_im+n), f); + f=(float) GSL_VECTOR_REAL(x->l, n); + SETFLOAT((x->list_l_re+n), f); + } + + mTXEigBang(x); + } +#else + post("mtx_eig: implementation requires gsl"); +#endif + +} + +void mtx_eig_setup (void) +{ + mtx_eig_class = class_new + (gensym("mtx_eig"), + (t_newmethod) newMTXEig, + (t_method) deleteMTXEig, + sizeof (MTXEig), + CLASS_DEFAULT, A_GIMME, 0); + class_addbang (mtx_eig_class, (t_method) mTXEigBang); + class_addmethod (mtx_eig_class, (t_method) mTXEigMatrix, gensym("matrix"), A_GIMME,0); +} + +void iemtx_eig_setup(void){ + mtx_eig_setup(); +} -- cgit v1.2.1