# Ublas :

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.

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