small fix

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

4 files changed, 104 insertions, 97 deletions

diff --git a/src/iemmatrix_utility.c b/src/iemmatrix_utility.c
index 0781226..553e308 100644
--- a/src/iemmatrix_utility.c
+++ b/src/iemmatrix_utility.c
@@ -416,3 +416,107 @@ void matrix_free(t_matrix *x)
   x->atombuffer=0;
   x->col=x->row=0;
 }
+
+
+/* some math */
+
+/*  invert a square matrix (row=col=rowcol) */
+/* if "error" is non-NULL, it's content will be set to 0 if the matrix was invertable, else to non-0 */
+t_matrixfloat* mtx_doInvert(t_matrixfloat*input, int rowcol, int*err){
+  /*
+   * 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; /* error counter */
+
+  int col=rowcol, row=rowcol, row2=row*col;
+  t_matrixfloat *original=input;
+  t_matrixfloat *inverted = 0;
+
+  if(input==0)return 0;
+
+  /* 1a reserve space for the inverted matrix */
+  inverted=(t_matrixfloat *)getbytes(sizeof(t_matrixfloat)*row2);
+  if(inverted==0)return 0;
+
+  /* 1b make an eye-shaped float-buf for B */
+  i=row2;
+  b1=inverted;
+  while(i--)*b1++=0;
+  i=row;
+  b1=inverted;
+  while(i--)b1[i*(row+1)]=1;
+
+  /* 2. do the Gauss-Jordan */
+  for (k=0;k<row;k++) {
+    /* adjust current row */
+    t_matrixfloat diagel = original[k*(col+1)];
+    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;
+    i=row;
+    while(i--){
+      *a2++*=i_diagel;
+      *b2++*=i_diagel;
+    }
+
+    /* eliminate the k-th element in each row by adding the weighted normalized row */
+    a2=original+k*row;
+    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);
+        }
+      }
+  }
+  if(err!=0)*err=ok;
+
+  return inverted;
+}
+
+/*  transpose a matrix */
+t_matrixfloat*mtx_doTranspose(t_matrixfloat*transposee, int row, int col){
+  int r,c;
+  t_matrixfloat*transposed=0;
+  if(!transposee||!row||!col)return 0;
+  transposed=(t_matrixfloat*)getbytes(sizeof(t_matrixfloat)*row*col);
+  r=row;
+  while(r--){
+    c=col;
+    while(c--)
+      transposed[c*row+r]=transposee[r*col+c];
+  }
+  return transposed;
+}
+
+/*  multiply matrix A=[rowA*colA] with matrix B=[rowB*colB]; C=A*B; colA=rowB=colArowB */
+t_matrixfloat*mtx_doMultiply(int rowA, t_matrixfloat*A, int colArowB, t_matrixfloat*B, int colB){
+  t_matrixfloat*result=0;
+  int r, c, n;
+
+  if(!A || !B || !rowA || !colArowB || !colB)return 0;
+  result=(t_matrixfloat*)getbytes(sizeof(t_matrixfloat)*rowA*colB);
+
+  for(r=0; r<rowA; r++){
+    for(c=0; c<colB; c++){
+      t_matrixfloat sum=0.f;
+      for(n=0;n<colArowB; n++)
+        sum+=A[colArowB*r+n]*B[colB*n+c];
+      result[colB*r+c]=sum;
+    }
+  }
+  return result;
+}
diff --git a/src/mtx_inverse.c b/src/mtx_inverse.c
index 2555835..f8c1d6a 100644
--- a/src/mtx_inverse.c
+++ b/src/mtx_inverse.c
@@ -16,72 +16,6 @@
 /* mtx_inverse */
 static t_class *mtx_inverse_class;
 
-
-t_matrixfloat* mtx_doInvert(t_matrixfloat*input, int rowcol, int*err){
-  /*
-   * 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; /* error counter */
-
-  int col=rowcol, row=rowcol, row2=row*col;
-  t_matrixfloat *original=input;
-  t_matrixfloat *inverted = 0;
-
-  if(input==0)return 0;
-
-  /* 1a reserve space for the inverted matrix */
-  inverted=(t_matrixfloat *)getbytes(sizeof(t_matrixfloat)*row2);
-  if(inverted==0)return 0;
-
-  /* 1b make an eye-shaped float-buf for B */
-  i=row2;
-  b1=inverted;
-  while(i--)*b1++=0;
-  i=row;
-  b1=inverted;
-  while(i--)b1[i*(row+1)]=1;
-
-  /* 2. do the Gauss-Jordan */
-  for (k=0;k<row;k++) {
-    /* adjust current row */
-    t_matrixfloat diagel = original[k*(col+1)];
-    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;
-    i=row;
-    while(i--){
-      *a2++*=i_diagel;
-      *b2++*=i_diagel;
-    }
-
-    /* eliminate the k-th element in each row by adding the weighted normalized row */
-    a2=original+k*row;
-    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);
-        }
-      }
-  }
-  if(err!=0)*err=ok;
-
-  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 ? */
diff --git a/src/mtx_mul.c b/src/mtx_mul.c
index bd6aace..a79b374 100644
--- a/src/mtx_mul.c
+++ b/src/mtx_mul.c
@@ -23,23 +23,6 @@
 
 */
 
-t_matrixfloat*mtx_doMultiply(int rowA, t_matrixfloat*A, int colArowB, t_matrixfloat*B, int colB){
-  t_matrixfloat*result=0;
-  int r, c, n;
-
-  if(!A || !B || !rowA || !colArowB || !colB)return 0;
-  result=(t_matrixfloat*)getbytes(sizeof(t_matrixfloat)*rowA*colB);
-
-  for(r=0; r<rowA; r++){
-    for(c=0; c<colB; c++){
-      t_matrixfloat sum=0.f;
-      for(n=0;n<colArowB; n++)
-        sum+=A[colArowB*r+n]*B[colB*n+c];
-      result[colB*r+c]=sum;
-    }
-  }
-  return result;
-}
 
 /* mtx_mul */
 static t_class *mtx_mul_class, *mtx_mulelement_class, *mtx_mulscalar_class;
diff --git a/src/mtx_transpose.c b/src/mtx_transpose.c
index c2ffc4e..4fcc87c 100644
--- a/src/mtx_transpose.c
+++ b/src/mtx_transpose.c
@@ -16,20 +16,6 @@
 /* mtx_transpose */
 static t_class *mtx_transpose_class;
 
-t_matrixfloat*mtx_doTranspose(t_matrixfloat*transposee, int row, int col){
-  int r,c;
-  t_matrixfloat*transposed=0;
-  if(!transposee||!row||!col)return 0;
-  transposed=(t_matrixfloat*)getbytes(sizeof(t_matrixfloat)*row*col);
-  r=row;
-  while(r--){
-    c=col;
-    while(c--)
-      transposed[c*row+r]=transposee[r*col+c];
-  }
-  return transposed;
-}
-
 static void mtx_transpose_matrix(t_matrix *x, t_symbol *s, int argc, t_atom *argv)
 {
   int row=atom_getfloat(argv++);