Stjepan Rajko wrote:
Philippe,
The following code snippets might be helpful - I used singular value decomposition using lapack bindings to get eigenstuff. I don't know if it's the same thing as eig() in your case.
ublas::matrix<double > A(rows, cols);
ublas::matrix<double > ATA; ATA = prod(trans(A),A);
ublas::matrix<double, boost::numeric::ublas::column_major > U(cols, cols), Vt(cols,cols); ublas::vector<double> S(cols); boost::numeric::bindings::lapack::gesvd(ATA, S, U, Vt);
cout << U << endl << endl; cout << Vt << endl << endl; cout << S << endl << endl;
It would be easier and more accurate to calculate the singular value decomposition of A itself, why bother with prod(trans(A),A) ? If you really do want to diagonalize prod(trans(A),A), a symmetric (and non-negative) matrix, the SVD is equivalent to an ordinary eigensolver and you would be better off using syev or syevr. Cheers, Ian