From: jk jk (jk.lists.questions_at_[hidden])
Date: 2008-08-28 15:06:09
It turns out I'm a fool. Sorry to bother you guys. It was the
column_major vs row_major thing again.
On Thu, Aug 28, 2008 at 2:46 PM, jk jk <jk.lists.questions_at_[hidden]> wrote:
> 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);
> }
>