From 8721acc0c0329ef6e63107c44e1bb61fae26f015 Mon Sep 17 00:00:00 2001 From: Franz Zotter Date: Fri, 9 Jan 2009 10:02:23 +0000 Subject: 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 --- doc/mtx_svd-help.pd | 23 +++++++ src/config.h.in | 3 + src/configure.ac | 7 +- src/mtx_svd.c | 185 ++++++++++++++++++++++++++++++++++++++++++++++++++++ 4 files changed, 217 insertions(+), 1 deletion(-) create mode 100644 doc/mtx_svd-help.pd create mode 100644 src/mtx_svd.c 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 + +#ifdef HAVE_LIBGSL +#include +#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_sizerows=rows; + x->columns=columns; + + deleteMTXusvw(x); + allocMTXusvw(x); + + for (n=0;nu->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;nlist_u+2+n), x->u->data[n]); + + for (n=0;ncolumns;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;nlist_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(); +} -- cgit v1.2.1