split the objects from mtx_binops into several files: mtx_add, mtx_sub, mtx_mul, mtx_pow

mtx_binops is still there for glue functions
exposed some important operations on float-arrays via iemmatrix.h (mtx_doInvert, mtx_doTranspose, mtx_doMultiply)
as you can easily see, they all start with "mtx_do" and a capital letter;
they all return a pointer to (newly allocated) memory with the result

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

1 files changed, 63 insertions, 36 deletions

diff --git a/src/mtx_inverse.c b/src/mtx_inverse.c
index c0c4019..f994b39 100644
--- a/src/mtx_inverse.c
+++ b/src/mtx_inverse.c
@@ -16,33 +16,23 @@
 /* mtx_inverse */
 static t_class *mtx_inverse_class;
 
-static void mtx_inverse_matrix(t_matrix *x, t_symbol *s, int argc, t_atom *argv)
-{
-  /* maybe we should do this in double or long double ? */
-  int row=atom_getfloat(argv);
-  int col=atom_getfloat(argv+1);
-  int i, k, row2=row*row;
 
-  t_matrixfloat *original, *inverted;
-  t_matrixfloat *a1, *a2, *b1, *b2;  // dummy pointers
+t_matrixfloat* mtx_doInvert(t_matrixfloat*input, int rowcol){
+  /*
+   * row==col==rowclo
+   * input=t_matrixfloat[row*col]
+   * output=t_matrixfloat[row*col]
+   */
+  int i, k;
+  t_matrixfloat *a1, *b1, *a2, *b2;
 
-  int ok = 0;
+  int ok=0; /* error counter */
 
-  if(row*col+2>argc){
-    post("mtx_print : sparse matrices not yet supported : use \"mtx_check\"");
-    return;
-  }
-  if (row!=col){
-    post("mtx_inverse: only square matrices can be inverted");
-    return;
-  }
+  int col=rowcol, row=rowcol, row2=row*col;
+  t_matrixfloat *original=input;
 
-  // reserve memory for outputting afterwards
-  adjustsize(x, row, row);
-  // 1. get the 2 matrices : orig; invert (create as eye, but will be orig^(-1))
-  inverted = (t_matrixfloat *)getbytes(sizeof(t_matrixfloat)*row2);
-  // 1a extract values of A to float-buf
-  original=matrix2float(argv);
+  // 1a reserve space for the inverted matrix
+  t_matrixfloat *inverted = (t_matrixfloat *)getbytes(sizeof(t_matrixfloat)*row2);
   // 1b make an eye-shaped float-buf for B
   i=row2;
   b1=inverted;
@@ -58,7 +48,6 @@ static void mtx_inverse_matrix(t_matrix *x, t_symbol *s, int argc, t_atom *argv)
     t_matrixfloat i_diagel = diagel?1./diagel:0;
     if (!diagel)ok++;
 
-
     /* normalize current row (set the diagonal-element to 1 */
     a2=original+k*col;
     b2=inverted+k*col;
@@ -73,23 +62,61 @@ static void mtx_inverse_matrix(t_matrix *x, t_symbol *s, int argc, t_atom *argv)
     b2=inverted+k*row;
     for(i=0;i<row;i++)
       if (i-k) {
-	t_matrixfloat f=-*(original+i*row+k);
-	int j = row;
-	a1=original+i*row;
-	b1=inverted+i*row;
-	while (j--) {
-	  *(a1+j)+=f**(a2+j);
-	  *(b1+j)+=f**(b2+j);
-	}
+        t_matrixfloat f=-*(original+i*row+k);
+        int j = row;
+        a1=original+i*row;
+        b1=inverted+i*row;
+        while (j--) {
+          *(a1+j)+=f**(a2+j);
+          *(b1+j)+=f**(b2+j);
+        }
       }
   }
+
+  if (ok)post("mtx_inverse: couldn't really invert the matrix !!! %d error%c", ok, (ok-1)?'s':0);
+
+  return inverted;
+}
+
+static void mtx_inverse_matrix(t_matrix *x, t_symbol *s, int argc, t_atom *argv)
+{
+  /* maybe we should do this in double or long double ? */
+  int row=atom_getfloat(argv);
+  int col=atom_getfloat(argv+1);
+
+  t_matrixfloat *original, *inverted;
+
+  if(row*col+2>argc){
+    post("mtx_print : sparse matrices not yet supported : use \"mtx_check\"");
+    return;
+  }
+
+  // 1 extract values of A to float-buf
+  original=matrix2float(argv);
+
+  // reserve memory for outputting afterwards
+  adjustsize(x, col, row);
+
+  if (row!=col){
+    // we'll have to do the pseudo-inverse:
+    // P=A'*inv(A*A');
+    t_matrixfloat*transposed=mtx_doTranspose(original, row, col);
+    t_matrixfloat*invertee  =mtx_doMultiply(row, original, col, transposed, row);
+    inverted=mtx_doMultiply(col, transposed, row, mtx_doInvert(invertee, row), row);
+
+    freebytes(transposed, sizeof(t_matrixfloat)*col*row);
+    freebytes(invertee  , sizeof(t_matrixfloat)*row*row);
+    
+  } else {
+    inverted=mtx_doInvert(original, row);
+  }
+
   // 3. output the matrix
   // 3a convert the floatbuf to an atombuf;
   float2matrix(x->atombuffer, inverted);
   // 3b destroy the buffers
-  freebytes(original, sizeof(t_matrixfloat)*row2);
-
-  if (ok)post("mtx_inverse: couldn't really invert the matrix !!! %d error%c", ok, (ok-1)?'s':0);
+  freebytes(original, sizeof(t_matrixfloat)*row*col);
+  freebytes(inverted, sizeof(t_matrixfloat)*row*row);
 
   // 3c output the atombuf;
   matrix_bang(x);
@@ -107,7 +134,7 @@ static void *mtx_inverse_new(t_symbol *s, int argc, t_atom *argv)
 void mtx_inverse_setup(void)
 {
   mtx_inverse_class = class_new(gensym("mtx_inverse"), (t_newmethod)mtx_inverse_new, 
-				(t_method)matrix_free, sizeof(t_matrix), 0, A_GIMME, 0);
+                                (t_method)matrix_free, sizeof(t_matrix), 0, A_GIMME, 0);
   class_addbang  (mtx_inverse_class, matrix_bang);
   class_addmethod(mtx_inverse_class, (t_method)mtx_inverse_matrix, gensym("matrix"), A_GIMME, 0);
   class_sethelpsymbol(mtx_inverse_class, gensym("iemmatrix/mtx_inverse"));