Update of /cvsroot/pure-data/externals/iem/iemmatrix/src In directory sc8-pr-cvs1.sourceforge.net:/tmp/cvs-serv24854
Modified Files: mtx_inverse.c Log Message: for non-square matrices, we now calculate automatically the (correct) pseudoinverse
Index: mtx_inverse.c =================================================================== RCS file: /cvsroot/pure-data/externals/iem/iemmatrix/src/mtx_inverse.c,v retrieving revision 1.3 retrieving revision 1.4 diff -C2 -d -r1.3 -r1.4 *** mtx_inverse.c 11 May 2005 13:05:29 -0000 1.3 --- mtx_inverse.c 11 May 2005 14:07:18 -0000 1.4 *************** *** 18,22 ****
! t_matrixfloat* mtx_doInvert(t_matrixfloat*input, int rowcol){ /* * row==col==rowclo --- 18,22 ----
! t_matrixfloat* mtx_doInvert(t_matrixfloat*input, int rowcol, int*error){ /* * row==col==rowclo *************** *** 31,38 **** int col=rowcol, row=rowcol, row2=row*col; t_matrixfloat *original=input;
! // 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; --- 31,43 ---- 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; *************** *** 42,48 **** while(i--)b1[i*(row+1)]=1;
! // 2. do the Gauss-Jordan for (k=0;k<row;k++) { ! // 2. adjust current row t_matrixfloat diagel = original[k*(col+1)]; t_matrixfloat i_diagel = diagel?1./diagel:0; --- 47,53 ---- 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; *************** *** 57,62 **** *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; --- 62,67 ---- *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; *************** *** 73,78 **** } } - if (ok)post("mtx_inverse: couldn't really invert the matrix !!! %d error%c", ok, (ok-1)?'s':0);
return inverted; --- 78,83 ---- } } if (ok)post("mtx_inverse: couldn't really invert the matrix !!! %d error%c", ok, (ok-1)?'s':0); + if(error!=0)*error=ok;
return inverted; *************** *** 92,123 **** }
! // 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)*row*col); - freebytes(inverted, sizeof(t_matrixfloat)*row*row);
! // 3c output the atombuf; matrix_bang(x); } --- 97,137 ---- }
! /* reserve memory for outputting afterwards */ adjustsize(x, col, row);
! /* 1. extract values of A to float-buf */ ! original=matrix2float(argv);
! if (row==col){ ! /* fine, the matrix is square */ ! inverted=mtx_doInvert(original, row, 0); } else { ! /* we'll have to do the pseudo-inverse: ! * P=A'*inv(A*A') if row<col ! * P=inv(A'*A)*A' if col<row ! */ ! t_matrixfloat*transposed, *invertee; ! int inverteeCol=0; ! transposed=mtx_doTranspose(original, row, col); ! if(row>col){ ! inverteeCol=col; ! invertee =mtx_doMultiply(col, transposed, row, original, col); ! inverted =mtx_doMultiply(col, mtx_doInvert(invertee, col, 0), col, transposed, row); ! } else { ! inverteeCol=row; ! invertee =mtx_doMultiply(row, original, col, transposed, row); ! inverted =mtx_doMultiply(col, transposed, row, mtx_doInvert(invertee, row, 0), row); ! } ! freebytes(transposed, sizeof(t_matrixfloat)*col*row); ! freebytes(invertee , sizeof(t_matrixfloat)*inverteeCol*inverteeCol); }
! /* 3. output the matrix */ ! /* 3a convert the floatbuf to an atombuf; */ float2matrix(x->atombuffer, inverted); ! /* 3b destroy the buffers */ freebytes(original, sizeof(t_matrixfloat)*row*col);
! /* 3c output the atombuf; */ matrix_bang(x); }