[ublas] issues with singular value decomposition

From: jk jk (jk.lists.questions_at_[hidden])
Date: 2008-08-28 14:46:40

I wrote the function pasted below to compute the
pseudo-inverse(http://en.wikipedia.org/wiki/Moore-Penrose_pseudoinverse)
via singular value decomposition of "inMatrix".

I've stared at this thing for hours trying to determine why sigma is
filling up with negative values.  Any help would be appreciated.
Thanks in advance.

---BEGIN PASTE-------
typedef ublas::matrix<double> matrixType;
typedef ublas::vector<double> vectorType;

void partialCalc::mpPseudoInverse(matrixType &inMatrix, matrixType &outMatrix){
  vectorType sigma(std::min(inMatrix.size1(), inMatrix.size2()));

  int ldu  = inMatrix.size1();
  int ucol = std::min(inMatrix.size1(),inMatrix.size2());
  matrixType u(ldu,ucol);

  int ldvt = std::min(inMatrix.size1(),inMatrix.size2());
  matrixType vt(ldvt, inMatrix.size1());

  vectorType w(lapack::gesvd_work('O','A', 'A', inMatrix));

  lapack::gesvd('S', 'S', inMatrix, sigma, u, vt, w);
  std::cout << u; exit(0);
  for ( int index =0; index < sigma.size(); index++)
    if ( sigma(index) != 0)
      sigma(index)=1/sigma(index);

  matrixType sigmam(zeroMatrixType(sigma.size(), sigma.size()));

  for ( unsigned int index =0; index < sigma.size(); index++)
    sigmam(index,index)=sigma(index);

  matrixType metaOutMatrix = prec_prod(herm(vt), sigmam);
  outMatrix = prec_prod(metaOutMatrix, herm(u));

  std::cout << outMatrix; exit(0);
}