On Mon, Nov 23, 2015 at 9:37 AM, Jonas Ballani <[email protected]>
wrote:
> Dear all,
>
> I am relatively new to libmesh and would like to report a problem using
> the SVD routine of the DenseMatrix class for rectangular matrices. My
> example program looks as follows:
>
> #include <iostream>
> #include <fstream>
> #include <algorithm>
> #include <math.h>
>
> #include "libmesh/libmesh.h"
> #include "libmesh/dense_vector.h"
> #include "libmesh/dense_matrix.h"
>
> #include <Eigen/SVD>
>
> using namespace libMesh;
>
> int main (int argc, char** argv)
> {
> LibMeshInit init (argc, argv);
>
> {
> std::cout << "## Libmesh DenseMatrix SVD" << std::endl;
>
> DenseMatrix< Number > A( 3, 2 );
> DenseMatrix< Number > U, VT;
> DenseVector< Number > sigma;
>
> A(0,0) = 2.0;
> A(1,0) = 0.0;
> A(2,0) = 0.0;
> A(0,1) = 0.0;
> A(1,1) = 1.0;
> A(2,1) = 0.0;
>
> std::cout << "A = " << std::endl << A << std::endl;
>
> A.svd( sigma, U, VT );
>
> std::cout << "U = " << std::endl << U << std::endl;
> std::cout << "VT = " << std::endl << VT << std::endl;
> std::cout << "sigma = " << std::endl << sigma << std::endl;
> }
>
> {
> std::cout << "## Eigen SVD" << std::endl;
>
> Eigen::MatrixXf A( 3, 2 );
>
> A(0,0) = 2.0;
> A(1,0) = 0.0;
> A(2,0) = 0.0;
> A(0,1) = 0.0;
> A(1,1) = 1.0;
> A(2,1) = 0.0;
>
> std::cout << "A = " << std::endl << A << std::endl << std::endl;
>
> Eigen::JacobiSVD< Eigen::MatrixXf > svd( A, Eigen::ComputeThinU |
> Eigen::ComputeThinV );
>
> std::cout << "U = " << std::endl << svd.matrixU() << std::endl <<
> std::endl;
> std::cout << "VT = " << std::endl << svd.matrixV() << std::endl <<
> std::endl;
> std::cout << "sigma = "<< std::endl << svd.singularValues() <<
> std::endl;
> }
>
> return 0;
> }
>
> The output I get is:
>
> ## Libmesh DenseMatrix SVD
> A =
> 2 0
> 0 1
> 0 0
>
> U =
> 1 0
> -0 1
>
> VT =
> 1 0 1
> 0 -0 0
>
> sigma =
> 2
> 1
>
> ## Eigen SVD
> A =
> 2 0
> 0 1
> 0 0
>
> U =
> 1 0
> 0 1
> 0 0
>
> VT =
> 1 0
> 0 1
>
> sigma =
> 2
> 1
>
> Obviously, the SVD from the DenseMatrix class is wrong, whereas Eigen
> performs well. It would be interesting to hear whether this is a known
> problem.
>
Thanks for the test code, we'll have to take a closer look. It looks like
DenseMatrix::svd() calls Lapack's dgesvd, which I doubt is wrong in itself,
but we might be calling it wrong.
I wonder if we ever tested it on non-square matrices...
--
John
------------------------------------------------------------------------------
Go from Idea to Many App Stores Faster with Intel(R) XDK
Give your users amazing mobile app experiences with Intel(R) XDK.
Use one codebase in this all-in-one HTML5 development environment.
Design, debug & build mobile apps & 2D/3D high-impact games for multiple OSs.
http://pubads.g.doubleclick.net/gampad/clk?id=254741551&iu=/4140
_______________________________________________
Libmesh-users mailing list
[email protected]
https://lists.sourceforge.net/lists/listinfo/libmesh-users
