Revision: 10502 http://pure-data.svn.sourceforge.net/pure-data/?rev=10502&view=rev Author: fzotter Date: 2009-01-10 18:33:39 +0000 (Sat, 10 Jan 2009)
Log Message: ----------- added [mtx_eig] to compute eigenvalues and eigenvectors using GSL. Something seems to be wrong with the eigenvectors: result differs from GNU octave result.
Modified Paths: -------------- trunk/externals/iem/iemmatrix/src/mtx_svd.c
Added Paths: ----------- trunk/externals/iem/iemmatrix/doc/mtx_eig-help.pd trunk/externals/iem/iemmatrix/src/mtx_eig.c
Added: trunk/externals/iem/iemmatrix/doc/mtx_eig-help.pd =================================================================== --- trunk/externals/iem/iemmatrix/doc/mtx_eig-help.pd (rev 0) +++ trunk/externals/iem/iemmatrix/doc/mtx_eig-help.pd 2009-01-10 18:33:39 UTC (rev 10502) @@ -0,0 +1,35 @@ +#N canvas 0 0 710 531 10; +#X obj 82 93 mtx_eig; +#X obj 82 66 mtx_rand; +#X obj 82 146 print l_re; +#X obj 128 124 print l_im; +#X msg 82 38 3; +#X obj 221 61 mtx_rand; +#X msg 221 33 3; +#X obj 230 119 mtx_eig v; +#X obj 230 192 print l_re; +#X obj 250 171 print l_im; +#X obj 314 193 mtx_print Q_re; +#X obj 335 171 mtx_print Q_im; +#X obj 230 98 t a a; +#X obj 375 123 mtx_print A; +#X msg 350 56 matrix 2 2 1 1 -1 2; +#X text 61 -33 [mtx_eig] computes (complex valued) eigenvalues , and +if desired eigenvectors of a matrix.; +#X text 425 21 requires GSL; +#X text 212 276 Franz Zotter , 2009; +#X text 98 238 Attention: eigenvectors seem to be wrong. Don't know +why!; +#X connect 0 0 2 0; +#X connect 0 1 3 0; +#X connect 1 0 0 0; +#X connect 4 0 1 0; +#X connect 5 0 12 0; +#X connect 6 0 5 0; +#X connect 7 0 8 0; +#X connect 7 1 9 0; +#X connect 7 2 10 0; +#X connect 7 3 11 0; +#X connect 12 0 7 0; +#X connect 12 1 13 0; +#X connect 14 0 12 0;
Added: trunk/externals/iem/iemmatrix/src/mtx_eig.c =================================================================== --- trunk/externals/iem/iemmatrix/src/mtx_eig.c (rev 0) +++ trunk/externals/iem/iemmatrix/src/mtx_eig.c 2009-01-10 18:33:39 UTC (rev 10502) @@ -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 <stdlib.h> + +#ifdef HAVE_LIBGSL +#include <gsl/gsl_eigen.h> +#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_size<size) + post("mtx_eig: sparse matrix not yet supported: use "mtx_check""); + else if (rows!=columns) + post("mtx_eig: Eigendecomposition works for square matrices only!"); + else { + size=rows; + x->size=size; + + deleteMTXqlw(x); + allocMTXqlw(x); + + for (n=0;n<in_size;n++) + x->a->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;n<in_size;n++) { + SETFLOAT((x->list_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;n<x->size;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(); +}
Modified: trunk/externals/iem/iemmatrix/src/mtx_svd.c =================================================================== --- trunk/externals/iem/iemmatrix/src/mtx_svd.c 2009-01-10 17:23:46 UTC (rev 10501) +++ trunk/externals/iem/iemmatrix/src/mtx_svd.c 2009-01-10 18:33:39 UTC (rev 10502) @@ -3,8 +3,9 @@ * * objects for manipulating simple matrices * mostly refering to matlab/octave matrix functions + * this functions depends on the GNU scientific library * - * Copyright (c) 2008, Franz Zotter + * Copyright (c) 2009, Franz Zotter * IEM, Graz, Austria * * For information on usage and redistribution, and for a DISCLAIMER OF ALL @@ -149,16 +150,16 @@ SETFLOAT((x->list_u),(float) x->rows); SETFLOAT((x->list_u+1),(float) x->columns); for (n=0;n<in_size;n++) - SETFLOAT((x->list_u+2+n), x->u->data[n]); + SETFLOAT((x->list_u+2+n), (float) x->u->data[n]);
for (n=0;n<x->columns;n++) - SETFLOAT((x->list_s+n), x->s->data[n]); + SETFLOAT((x->list_s+n),(float) x->s->data[n]);
SETFLOAT((x->list_v),(float) x->columns); SETFLOAT((x->list_v+1),(float) x->columns); in_size=x->columns*x->columns; for (n=0;n<in_size;n++) - SETFLOAT((x->list_v+n+2), x->v->data[n]); + SETFLOAT((x->list_v+n+2), (float) x->v->data[n]);
mTXSvdBang(x); }
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.