Update of /cvsroot/pure-data/externals/iem/iemmatrix/src In directory sc8-pr-cvs1.sourceforge.net:/tmp/cvs-serv29751

Modified Files: iemmatrix_utility.c mtx_inverse.c mtx_mul.c mtx_transpose.c Log Message: small fix

Index: mtx_transpose.c =================================================================== RCS file: /cvsroot/pure-data/externals/iem/iemmatrix/src/mtx_transpose.c,v retrieving revision 1.5 retrieving revision 1.6 diff -C2 -d -r1.5 -r1.6 *** mtx_transpose.c 29 Mar 2006 12:07:52 -0000 1.5 --- mtx_transpose.c 15 Jan 2007 20:06:42 -0000 1.6 *************** *** 17,34 **** 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) { --- 17,20 ----

Index: mtx_inverse.c =================================================================== RCS file: /cvsroot/pure-data/externals/iem/iemmatrix/src/mtx_inverse.c,v retrieving revision 1.8 retrieving revision 1.9 diff -C2 -d -r1.8 -r1.9 *** mtx_inverse.c 12 Nov 2006 12:56:06 -0000 1.8 --- mtx_inverse.c 15 Jan 2007 20:06:42 -0000 1.9 *************** *** 17,86 **** 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) { --- 17,20 ----

Index: mtx_mul.c =================================================================== RCS file: /cvsroot/pure-data/externals/iem/iemmatrix/src/mtx_mul.c,v retrieving revision 1.5 retrieving revision 1.6 diff -C2 -d -r1.5 -r1.6 *** mtx_mul.c 25 Jul 2006 14:15:22 -0000 1.5 --- mtx_mul.c 15 Jan 2007 20:06:42 -0000 1.6 *************** *** 24,44 **** */

- 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 */ --- 24,27 ----

Index: iemmatrix_utility.c =================================================================== RCS file: /cvsroot/pure-data/externals/iem/iemmatrix/src/iemmatrix_utility.c,v retrieving revision 1.1 retrieving revision 1.2 diff -C2 -d -r1.1 -r1.2 *** iemmatrix_utility.c 11 Jan 2007 19:41:20 -0000 1.1 --- iemmatrix_utility.c 15 Jan 2007 20:06:42 -0000 1.2 *************** *** 417,418 **** --- 417,522 ---- 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; + }