Revision: 10493 http://pure-data.svn.sourceforge.net/pure-data/?rev=10493&view=rev Author: fzotter Date: 2009-01-09 10:02:23 +0000 (Fri, 09 Jan 2009)
Log Message: ----------- added [mtx_svd] for singular value decomposition with the GNU scientific library support. config.h.in and configure.ac correspondingly modified.
Modified Paths: -------------- trunk/externals/iem/iemmatrix/src/config.h.in trunk/externals/iem/iemmatrix/src/configure.ac
Added Paths: ----------- trunk/externals/iem/iemmatrix/doc/mtx_svd-help.pd trunk/externals/iem/iemmatrix/src/mtx_svd.c
Added: trunk/externals/iem/iemmatrix/doc/mtx_svd-help.pd =================================================================== --- trunk/externals/iem/iemmatrix/doc/mtx_svd-help.pd (rev 0) +++ trunk/externals/iem/iemmatrix/doc/mtx_svd-help.pd 2009-01-09 10:02:23 UTC (rev 10493) @@ -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;
Modified: trunk/externals/iem/iemmatrix/src/config.h.in =================================================================== --- trunk/externals/iem/iemmatrix/src/config.h.in 2009-01-09 05:10:16 UTC (rev 10492) +++ trunk/externals/iem/iemmatrix/src/config.h.in 2009-01-09 10:02:23 UTC (rev 10493) @@ -7,5 +7,8 @@ /* do we have libsndfile installed? */ #undef HAVE_SNDFILE_H
+/* do we have libgsl installed? */ +#undef HAVE_LIBGSL
+ #endif /* CONFIG_H_ */
Modified: trunk/externals/iem/iemmatrix/src/configure.ac =================================================================== --- trunk/externals/iem/iemmatrix/src/configure.ac 2009-01-09 05:10:16 UTC (rev 10492) +++ trunk/externals/iem/iemmatrix/src/configure.ac 2009-01-09 10:02:23 UTC (rev 10493) @@ -72,7 +72,12 @@ 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 INCLUDES="-I${with_pd}/src ${INCLUDES}" @@ -99,7 +104,7 @@
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
Added: trunk/externals/iem/iemmatrix/src/mtx_svd.c =================================================================== --- trunk/externals/iem/iemmatrix/src/mtx_svd.c (rev 0) +++ trunk/externals/iem/iemmatrix/src/mtx_svd.c 2009-01-09 10:02:23 UTC (rev 10493) @@ -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(); +}
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.