[PD-cvs] SF.net SVN: pure-data:[10502] trunk/externals/iem/iemmatrix

fzotter at users.sourceforge.net fzotter at users.sourceforge.net
Sat Jan 10 19:33:39 CET 2009


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.




More information about the Pd-cvs mailing list