aboutsummaryrefslogtreecommitdiff
path: root/src/mtx_eig.c
diff options
context:
space:
mode:
Diffstat (limited to 'src/mtx_eig.c')
-rw-r--r--src/mtx_eig.c227
1 files changed, 227 insertions, 0 deletions
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 <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();
+}