aboutsummaryrefslogtreecommitdiff
diff options
context:
space:
mode:
authorGeorg Holzmann <grholzi@users.sourceforge.net>2007-01-15 20:06:42 +0000
committerGeorg Holzmann <grholzi@users.sourceforge.net>2007-01-15 20:06:42 +0000
commit3a22d644da666c934c4cf218687da5814c7b43bb (patch)
tree4d3c52f8bdabd5c1b98ab25354625814da976c53
parent445cc78fc9572a924af9b04d742fcdb46b281bfb (diff)
small fix
svn path=/trunk/externals/iem/iemmatrix/; revision=7338
-rw-r--r--src/iemmatrix_utility.c104
-rw-r--r--src/mtx_inverse.c66
-rw-r--r--src/mtx_mul.c17
-rw-r--r--src/mtx_transpose.c14
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++);