diff options
Diffstat (limited to 'src/mtx_inverse.c')
-rw-r--r-- | src/mtx_inverse.c | 66 |
1 files changed, 0 insertions, 66 deletions
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 ? */ |