added [mtx_svd] for singular value decomposition with the GNU scientific library support. config.h.in and configure.ac correspondingly modified.

svn path=/trunk/externals/iem/iemmatrix/; revision=10493

4 files changed, 217 insertions, 1 deletions

diff --git a/doc/mtx_svd-help.pd b/doc/mtx_svd-help.pd
new file mode 100644
index 0000000..7fd0f3a
--- /dev/null
+++ b/doc/mtx_svd-help.pd
@@ -0,0 +1,23 @@
+#N canvas 482 188 717 423 12;
+#X obj 93 216 mtx_svd;
+#X obj 93 108 mtx_rand;
+#X obj 309 275 mtx_transpose;
+#X obj 201 273 mtx_diag;
+#X obj 93 131 t a a;
+#X text 64 24 [mtx_svd] singular value decomposition using gsl;
+#X text 442 338 Franz Zotter \, 2009;
+#X text 353 153 A = U * S * V^T;
+#X obj 93 298 mtx_print U;
+#X obj 201 299 mtx_print S;
+#X obj 309 298 mtx_print V^T;
+#X obj 169 155 mtx_print A;
+#X msg 93 85 4 3;
+#X connect 0 0 8 0;
+#X connect 0 1 3 0;
+#X connect 0 2 2 0;
+#X connect 1 0 4 0;
+#X connect 2 0 10 0;
+#X connect 3 0 9 0;
+#X connect 4 0 0 0;
+#X connect 4 1 11 0;
+#X connect 12 0 1 0;
diff --git a/src/config.h.in b/src/config.h.in
index f4a52f9..ce9d7bf 100644
--- a/src/config.h.in
+++ b/src/config.h.in
@@ -7,5 +7,8 @@
 /* do we have libsndfile installed? */
 #undef HAVE_SNDFILE_H
 
+/* do we have libgsl installed? */
+#undef HAVE_LIBGSL
+
 
 #endif /* CONFIG_H_ */
diff --git a/src/configure.ac b/src/configure.ac
index ca98ab4..7128c61 100644
--- a/src/configure.ac
+++ b/src/configure.ac
@@ -72,6 +72,11 @@ AC_CHECK_LIB(fftw3, fftw_destroy_plan)
 dnl for soundfile reading (and writing, if we do that...)
 AC_CHECK_LIB(sndfile, sf_close)
 
+dnl for gnu scientific library -lgsl:
+dnl AC_CHECK_LIB([m], [cos])
+AC_CHECK_LIB([gslcblas], [cblas_dgemm])
+AC_CHECK_LIB([gsl], [gsl_blas_dgemm])
+
 
 if test "x$with_pd" != "x"; then
  if test -d "${with_pd}/src"; then
@@ -99,7 +104,7 @@ AC_CHECK_LIB(pd, nullfn)
 
 dnl Checks for header files.
 AC_HEADER_STDC
-AC_CHECK_HEADERS(stdlib.h stdio.h string.h math.h time.h sys/time.h fftw3.h sndfile.h)
+AC_CHECK_HEADERS(stdlib.h stdio.h string.h math.h time.h sys/time.h fftw3.h sndfile.h gsl/gsl_linalg.h)
 
 dnl Checks for typedefs, structures, and compiler characteristics.
 AC_HEADER_TIME
diff --git a/src/mtx_svd.c b/src/mtx_svd.c
new file mode 100644
index 0000000..8e784ad
--- /dev/null
+++ b/src/mtx_svd.c
@@ -0,0 +1,185 @@
+/*
+ *  iemmatrix
+ *
+ *  objects for manipulating simple matrices
+ *  mostly refering to matlab/octave matrix functions
+ *
+ * Copyright (c) 2008, 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_linalg.h>
+#endif
+
+static t_class *mtx_svd_class;
+
+typedef struct _MTXSvd_ MTXSvd;
+struct _MTXSvd_
+{
+  t_object x_obj;
+#ifdef HAVE_LIBGSL
+  gsl_matrix *u;
+  gsl_vector *s;
+  gsl_matrix *v;
+  gsl_vector *w;
+#endif
+  t_outlet *list_u_out;
+  t_outlet *list_s_out;
+  t_outlet *list_v_out;
+  t_atom *list_u;
+  t_atom *list_s;
+  t_atom *list_v;
+  int rows;
+  int columns;
+};
+
+#ifdef HAVE_LIBGSL
+static void allocMTXusvw (MTXSvd *x) 
+{
+     x->u=(gsl_matrix*)gsl_matrix_alloc(x->rows,x->columns);
+     x->s=(gsl_vector*)gsl_vector_alloc(x->columns);
+     x->v=(gsl_matrix*)gsl_matrix_alloc(x->columns,x->columns);
+     x->w=(gsl_vector*)gsl_vector_alloc(x->columns);
+
+     x->list_u=(t_atom*)calloc(sizeof(t_atom),x->rows*x->columns+2);
+     x->list_s=(t_atom*)calloc(sizeof(t_atom),x->columns);
+     x->list_v=(t_atom*)calloc(sizeof(t_atom),x->columns*x->columns+2);
+}
+
+static void deleteMTXusvw (MTXSvd *x) 
+{
+   if (x->list_u!=0)
+      free(x->list_u);
+   if (x->list_s!=0)
+      free(x->list_s);
+   if (x->list_v!=0)
+      free(x->list_v);
+
+   x->list_u = x->list_s = x->list_v = 0;
+
+   if (x->u!=0)
+      gsl_matrix_free(x->u);
+   if (x->s!=0)
+      gsl_vector_free(x->s);
+   if (x->v!=0)
+      gsl_matrix_free(x->v);
+   if (x->w!=0)
+      gsl_vector_free(x->w);
+
+   x->u = 0;
+   x->s = 0;
+   x->v = 0;
+   x->w = 0;
+}
+#endif
+
+static void deleteMTXSvd (MTXSvd *x) 
+{
+#ifdef HAVE_LIBGSL
+   deleteMTXSvd(x);
+#endif
+}
+
+static void *newMTXSvd (t_symbol *s, int argc, t_atom *argv)
+{
+  MTXSvd *x = (MTXSvd *) pd_new (mtx_svd_class);
+  x->list_u_out = outlet_new (&x->x_obj, gensym("matrix"));
+  x->list_s_out = outlet_new (&x->x_obj, gensym("list"));
+  x->list_v_out = outlet_new (&x->x_obj, gensym("matrix"));
+  x->list_u = 0; 
+  x->list_s = 0;
+  x->list_v = 0;
+#ifdef HAVE_LIBGSL
+  x->u=0;
+  x->s=0;
+  x->v=0;
+  x->w=0;
+#endif
+  
+  return ((void *) x);
+} 
+
+static void mTXSvdBang (MTXSvd *x)
+{
+  if (x->list_u) {
+    outlet_anything(x->list_v_out, gensym("matrix"), x->columns*x->columns+2, x->list_v);
+    outlet_anything(x->list_s_out, gensym("list"), x->columns, x->list_s);
+    outlet_anything(x->list_u_out, gensym("matrix"), x->rows*x->columns+2, x->list_u);
+  }
+}
+
+static void mTXSvdMatrix (MTXSvd *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;
+
+
+#ifdef HAVE_LIBGSL
+  /* size check */
+  if (!size) 
+    post("mtx_svd: invalid dimensions");
+  else if (in_size<size) 
+    post("mtx_svd: sparse matrix not yet supported: use \"mtx_check\"");
+  else if (rows<columns)
+     post("mtx_svd: gsl_linalg_SVD_decomp does not support M<N");
+  else {
+     x->rows=rows;
+     x->columns=columns;
+
+    deleteMTXusvw(x);
+    allocMTXusvw(x);
+
+    for (n=0;n<in_size;n++) 
+       x->u->data[n]=(double) atom_getfloat(argv++);
+    
+    gsl_linalg_SV_decomp(x->u,x->v,x->s,x->w);
+
+    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]);
+
+    for (n=0;n<x->columns;n++)
+       SETFLOAT((x->list_s+n), 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]);
+
+    mTXSvdBang(x);
+  }
+#else
+    post("mtx_svd: implementation requires gsl");
+#endif
+
+}
+
+void mtx_svd_setup (void)
+{
+  mtx_svd_class = class_new 
+    (gensym("mtx_svd"),
+     (t_newmethod) newMTXSvd,
+     (t_method) deleteMTXSvd,
+     sizeof (MTXSvd),
+     CLASS_DEFAULT, A_GIMME, 0);
+  class_addbang (mtx_svd_class, (t_method) mTXSvdBang);
+  class_addmethod (mtx_svd_class, (t_method) mTXSvdMatrix, gensym("matrix"), A_GIMME,0);
+}
+
+void iemtx_svd_setup(void){
+  mtx_svd_setup();
+}