[PD-cvs] externals/iem/iemmatrix/src mtx_inverse.c,1.3,1.4

[IOhannes m zmölnig]

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);
  }