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.

Best regards,
Jonas
------------------------------------------------------------------------------
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
Libmesh-users@lists.sourceforge.net
https://lists.sourceforge.net/lists/listinfo/libmesh-users

Reply via email to