aboutsummaryrefslogtreecommitdiff
diff options
context:
space:
mode:
authorFranz Zotter <fzotter@users.sourceforge.net>2009-06-09 20:01:55 +0000
committerFranz Zotter <fzotter@users.sourceforge.net>2009-06-09 20:01:55 +0000
commit678b221a9b3323cc66e13b31ae1f53fa1138f32f (patch)
tree5058527e095bd22eaf0a473988bca82cf5b800ff
parentecbfff1732bf4e2e555e648aa81b9f1ca4ac94c2 (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.pd17
-rw-r--r--src/iemmatrix_sources.c3
-rw-r--r--src/iemmatrix_sources.h3
-rw-r--r--src/mtx_qr.c198
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();
+}