diff options
author | Franz Zotter <fzotter@users.sourceforge.net> | 2009-06-09 20:01:55 +0000 |
---|---|---|
committer | Franz Zotter <fzotter@users.sourceforge.net> | 2009-06-09 20:01:55 +0000 |
commit | 678b221a9b3323cc66e13b31ae1f53fa1138f32f (patch) | |
tree | 5058527e095bd22eaf0a473988bca82cf5b800ff | |
parent | ecbfff1732bf4e2e555e648aa81b9f1ca4ac94c2 (diff) |
adding [mtx_qr] using gsl support; performs QR-decomposition, but Q is not implemented yet.
svn path=/trunk/externals/iem/iemmatrix/; revision=11727
-rw-r--r-- | doc/mtx_qr-help.pd | 17 | ||||
-rw-r--r-- | src/iemmatrix_sources.c | 3 | ||||
-rw-r--r-- | src/iemmatrix_sources.h | 3 | ||||
-rw-r--r-- | src/mtx_qr.c | 198 |
4 files changed, 219 insertions, 2 deletions
diff --git a/doc/mtx_qr-help.pd b/doc/mtx_qr-help.pd new file mode 100644 index 0000000..eb88ea1 --- /dev/null +++ b/doc/mtx_qr-help.pd @@ -0,0 +1,17 @@ +#N canvas 0 0 670 367 10; +#X obj 93 216 mtx_qr; +#X obj 93 108 mtx_rand; +#X obj 93 131 t a a; +#X text 442 338 Franz Zotter \, 2009; +#X obj 169 155 mtx_print A; +#X msg 93 85 4 3; +#X obj 201 298 mtx_print R; +#X obj 93 298 mtx_print Q; +#X text 353 153 A = Q * R; +#X text 64 24 [mtx_qr] qr decomposition using gsl; +#X connect 0 0 7 0; +#X connect 0 1 6 0; +#X connect 1 0 2 0; +#X connect 2 0 0 0; +#X connect 2 1 4 0; +#X connect 5 0 1 0; diff --git a/src/iemmatrix_sources.c b/src/iemmatrix_sources.c index b65ca8e..a246ec1 100644 --- a/src/iemmatrix_sources.c +++ b/src/iemmatrix_sources.c @@ -57,8 +57,8 @@ void iemmatrix_sources_setup(void) iemtx_mean_setup(); /* mtx_mean.c */ iemtx_min2_setup(); /* mtx_min2.c */ iemtx_minmax_setup(); /* mtx_minmax.c */ - iemtx_mul_setup(); /* mtx_mul.c */ iemtx_mul__setup(); /* mtx_mul~.c */ + iemtx_mul_setup(); /* mtx_mul.c */ iemtx_neq_setup(); /* mtx_neq.c */ iemtx_not_setup(); /* mtx_not.c */ iemtx_ones_setup(); /* mtx_ones.c */ @@ -69,6 +69,7 @@ void iemmatrix_sources_setup(void) iemtx_powtodb_setup(); /* mtx_powtodb.c */ iemtx_print_setup(); /* mtx_print.c */ iemtx_prod_setup(); /* mtx_prod.c */ + iemtx_qr_setup(); /* mtx_qr.c */ iemtx_rand_setup(); /* mtx_rand.c */ iemtx_repmat_setup(); /* mtx_repmat.c */ iemtx_resize_setup(); /* mtx_resize.c */ diff --git a/src/iemmatrix_sources.h b/src/iemmatrix_sources.h index 3ccd923..cda8fc5 100644 --- a/src/iemmatrix_sources.h +++ b/src/iemmatrix_sources.h @@ -55,8 +55,8 @@ void iemtx_max2_setup(void); /* mtx_max2.c */ void iemtx_mean_setup(void); /* mtx_mean.c */ void iemtx_min2_setup(void); /* mtx_min2.c */ void iemtx_minmax_setup(void); /* mtx_minmax.c */ -void iemtx_mul_setup(void); /* mtx_mul.c */ void iemtx_mul__setup(void); /* mtx_mul~.c */ +void iemtx_mul_setup(void); /* mtx_mul.c */ void iemtx_neq_setup(void); /* mtx_neq.c */ void iemtx_not_setup(void); /* mtx_not.c */ void iemtx_ones_setup(void); /* mtx_ones.c */ @@ -67,6 +67,7 @@ void iemtx_pow_setup(void); /* mtx_pow.c */ void iemtx_powtodb_setup(void); /* mtx_powtodb.c */ void iemtx_print_setup(void); /* mtx_print.c */ void iemtx_prod_setup(void); /* mtx_prod.c */ +void iemtx_qr_setup(void); /* mtx_qr.c */ void iemtx_rand_setup(void); /* mtx_rand.c */ void iemtx_repmat_setup(void); /* mtx_repmat.c */ void iemtx_resize_setup(void); /* mtx_resize.c */ diff --git a/src/mtx_qr.c b/src/mtx_qr.c new file mode 100644 index 0000000..ef09d0f --- /dev/null +++ b/src/mtx_qr.c @@ -0,0 +1,198 @@ +/* + * 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_linalg.h> +#endif + +static t_class *mtx_qr_class; + +typedef struct _MTXQr_ MTXQr; +struct _MTXQr_ +{ + t_object x_obj; +#ifdef HAVE_LIBGSL + gsl_matrix *a; + gsl_vector *tau; +#endif + t_outlet *list_q_out; + t_outlet *list_r_out; + t_atom *list_q; + t_atom *list_r; + int rows; + int columns; +}; + +#ifdef HAVE_LIBGSL +static void allocMTXqr (MTXQr *x) +{ + x->a=(gsl_matrix*)gsl_matrix_alloc(x->rows,x->columns); + x->tau=(gsl_vector*)gsl_vector_alloc( + ((x->columns<x->rows)?x->columns:x->rows)); + + x->list_q=(t_atom*)calloc(sizeof(t_atom),x->rows*x->rows+2); + x->list_r=(t_atom*)calloc(sizeof(t_atom),x->rows*x->columns+2); +} + +static void deleteMTXqr (MTXQr *x) +{ + if (x->list_q!=0) + free(x->list_q); + if (x->list_r!=0) + free(x->list_r); + + x->list_q = x->list_r = 0; + + if (x->a!=0) + gsl_matrix_free(x->a); + if (x->tau!=0) + gsl_vector_free(x->tau); + + x->a = 0; + x->tau = 0; +} +#endif + +static void deleteMTXQr (MTXQr *x) +{ +#ifdef HAVE_LIBGSL + deleteMTXqr(x); +#endif +} + +static void *newMTXQr (t_symbol *s, int argc, t_atom *argv) +{ + MTXQr *x = (MTXQr *) pd_new (mtx_qr_class); + x->list_q_out = outlet_new (&x->x_obj, gensym("matrix")); + x->list_r_out = outlet_new (&x->x_obj, gensym("matrix")); + x->list_q = 0; + x->list_r = 0; +#ifdef HAVE_LIBGSL + x->a=0; + x->tau=0; +#endif + + return ((void *) x); +} + +static void mTXQrBang (MTXQr *x) +{ + if (x->list_q) { + outlet_anything(x->list_r_out, gensym("matrix"), x->rows*x->columns+2, x->list_r); + post("mtx_qr: implementation outputs only R currently! Q has to be implemented..."); + // outlet_anything(x->list_q_out, gensym("matrix"), x->rows*x->rows+2, x->list_q); + } +} + +static void mTXQrMatrix (MTXQr *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 m,n; + + +#ifdef HAVE_LIBGSL + /* size check */ + if (!size) + post("mtx_qr: invalid dimensions"); + else if (in_size<size) + post("mtx_qr: sparse matrix not yet supported: use \"mtx_check\""); + else if (rows<columns) + post("mtx_qr: gsl_linalg_SVD_decomp does not support M<N"); + else { + x->rows=rows; + x->columns=columns; + + deleteMTXqr(x); + allocMTXqr(x); + + for (n=0;n<in_size;n++) + x->a->data[n]=(double) atom_getfloat(argv++); + + gsl_linalg_QR_decomp(x->a,x->tau); + + + SETFLOAT((x->list_r),(float) x->rows); + SETFLOAT((x->list_r+1),(float) x->columns); + for (n=0,in_size=0;n<x->rows;n++) { + for (m=0;m<n;m++) { + SETFLOAT((x->list_r+2+in_size), 0); + in_size++; + } + for (;m<x->columns;m++) { + SETFLOAT((x->list_r+2+in_size), (float) x->a->data[in_size]); + in_size++; + } + } + + SETFLOAT((x->list_q),(float) x->rows); + SETFLOAT((x->list_q+1),(float) x->rows); + +// TODO: Housholder transformations have to be decoded from +// +// x->tau and lower triangular part of x->a. +// +// with L=min(rows,columns) +// +// Matrix multiplications have to be done: +// +// Q=QL QL-1 ... Q1 +// +// using the matrix factors +// +// Qi = I - tau[i] vi^T vi +// +// employing the dyadic vector product of +// +// [ 0 ] +// [ : ] +// [ 0 ] +// vi=[A[i+1,i]] +// [A[i+1,i]] +// [ : ] +// [ A[L,i] ] +// +// on itself (out of x->a), +// and the scalar tau[i] in the vector x->tau. + + mTXQrBang(x); + } +#else + post("mtx_qr: implementation requires gsl"); +#endif + +} + +void mtx_qr_setup (void) +{ + mtx_qr_class = class_new + (gensym("mtx_qr"), + (t_newmethod) newMTXQr, + (t_method) deleteMTXQr, + sizeof (MTXQr), + CLASS_DEFAULT, A_GIMME, 0); + class_addbang (mtx_qr_class, (t_method) mTXQrBang); + class_addmethod (mtx_qr_class, (t_method) mTXQrMatrix, gensym("matrix"), A_GIMME,0); +} + +void iemtx_qr_setup(void){ + mtx_qr_setup(); +} |