[Trilinos-Users] Anasazi:: sign of EigenVector

Holger Brandsmeier holger.brandsmeier at sam.math.ethz.ch
Fri May 31 04:01:36 MDT 2013


Sunghwan,

the nature of eigenvlaue problems is that eigenvectors are _not_
unique. Given the EVP
 A u = lambda u
If u is an eigenvector then also alpha * u is an eigenvector to the
same eigenvalue for _any_ alpha from the field of complex numbers.

Aas eigenvectors are not unique every library can give you a different
one. As long as they are all scalar multiples of each other, there is
no error. Usually libraries chose to give you and eigenvectors whose
l_2-norm is 1. For real eigenvectos that makes it unique up toa
multiplication by -1. That is what you described.

Just accept this fact. If you need a certain sign, just multiply the
eigenvector with any scalar. E.g. you can always multiply the vector
such that its first component is positive.

-Holger

On Fri, May 31, 2013 at 6:14 AM, Sunghwan Choi <sunghwanchoi91 at gmail.com> wrote:
> Dear all,
>
> Hi I had a test on Anasazi. I got test matrix from galeri and used
> BlockKrylovSchurSol.
>
> The eigenvalues that I got had no problem but the eigenvectors has problem.
> The sign is opposite. The magnitude in vectors is the same to the exact
> solution but it is positive values;;;;
>
> I don’t know how to solve this problems if you know solutions of it or got
> same troubles before, please let me know.
>
> I attached my test code. If you run this code, you can see first 5
> eigenvalues and first eigenvectors, whose sign is opposite to exact values.
>
> I got eigenvector and eigenvalues from scipy, which is python numeric
> library. Here is what I got from scipy package.
>
> eigenvalues
>
> 0.205758492492
>
> 0.823046639485
>
> 0.823046639485
>
> 1.82082931366
>
> 2.39125149816
>
> First eigenvector
>
> [-0.02106803 -0.04866181 -0.07041662 -0.07855108 -0.07041662 -0.04866181
>
> -0.02106803 -0.04866181 -0.11036498 -0.15856626 -0.17654728 -0.15856626
>
> -0.11036498 -0.04866181 -0.07041662 -0.15856626 -0.22708042 -0.25259643
>
> -0.22708042 -0.15856626 -0.07041662 -0.07855108 -0.17654728 -0.25259643
>
> -0.28090257 -0.25259643 -0.17654728 -0.07855108 -0.07041662 -0.15856626
>
> -0.22708042 -0.25259643 -0.22708042 -0.15856626 -0.07041662 -0.04866181
>
> -0.11036498 -0.15856626 -0.17654728 -0.15856626 -0.11036498 -0.04866181
>
> -0.02106803 -0.04866181 -0.07041662 -0.07855108 -0.07041662 -0.04866181
>
> -0.02106803]
>
>
>
> I think there is no error in scipy because I already tested the solution
> satisfy with \lambda V = M V
>
> Please help me. To me, this is quite important problem
>
>
>
> Sunghwan Choi
>
>
>
>
>
>
>
>
>
> #include "AnasaziBasicEigenproblem.hpp"
>
> #include "AnasaziEpetraAdapter.hpp"
>
> #include "AnasaziBlockKrylovSchurSolMgr.hpp"
>
>
>
> typedef Epetra_MultiVector MV;
>
> typedef Epetra_Operator OP;
>
> typedef Anasazi::MultiVecTraits<double, Epetra_MultiVector> MVT;
>
>
>
> #include "Galeri_Maps.h"
>
> #include "Galeri_CrsMatrices.h"
>
> #include "Galeri_Utils.h"
>
> #ifdef HAVE_MPI
>
> #include "Epetra_MpiComm.h"
>
> #include "mpi.h"
>
> #else
>
> #include "Epetra_SerialComm.h"
>
> #endif
>
> #include "Epetra_Map.h"
>
> #include "Epetra_CrsMatrix.h"
>
> #include "Teuchos_ParameterList.hpp"
>
>
>
> using namespace Galeri;
>
>
>
> // =========== //
>
> // main driver //
>
> // =========== //
>
>
>
> int main(int argc, char* argv[])
>
> {
>
> #ifdef HAVE_MPI
>
>   MPI_Init(&argc, &argv);
>
>   Epetra_MpiComm Comm(MPI_COMM_WORLD);
>
> #else
>
>   Epetra_SerialComm Comm;
>
> #endif
>
>
>
>   // pointer to the map to be created
>
>   Epetra_Map*            Map;
>
>   // pointer to the matrix to be created
>
>   Epetra_CrsMatrix*      Matrix;
>
>   // container for parameters
>
>   Teuchos::ParameterList GaleriList;
>
>   // here we specify the global dimension of the problem
>
>   cout << "NumProc?:" <<Comm.NumProc() <<endl;
>
>
>
>   int nx = 7 * Comm.NumProc();
>
>   int ny = 7 * Comm.NumProc();
>
>
>
>   GaleriList.set("nx", nx);
>
> GaleriList.set("ny", ny);
>
>
>
> try
>
>   {
>
>     // Creates a simple linear map; for more details on the map creation
>
>     // refer to the documentation
>
>     Map = CreateMap("Cartesian2D", Comm, GaleriList);
>
>
>
>     // Creates a diagonal matrix with 1's on the diagonal
>
>     Matrix   = CreateCrsMatrix("Biharmonic2D", Map, GaleriList);
>
>
>
>     // print out the matrix
>
>     cout << *Matrix;
>
>
>
>   }
>
>   catch (Exception& rhs)
>
>   {
>
>     if (Comm.MyPID() == 0)
>
>     {
>
>       cerr << "Caught exception: ";
>
>       rhs.Print();
>
>     }
>
>   }
>
>
>
>   Teuchos::ParameterList* parameters = new Teuchos::ParameterList();
>
>   parameters->sublist("Pure_Diagonalize").set("NumberOfEigenvalues",5);
>
>   parameters->sublist("Pure_Diagonalize").set("Tolerance",0.00001);
>
>   parameters->sublist("Pure_Diagonalize").set("BlockSize",7);
>
>   parameters->sublist("Pure_Diagonalize").set("NumBlocks",5);
>
>
>
>   Teuchos::RCP < Anasazi::BasicEigenproblem<double,MV,OP> > MyProblem
> =Teuchos::rcp( new Anasazi::BasicEigenproblem<double,MV,OP>() );
>
>   MyProblem->setNEV(
> parameters->sublist("Pure_Diagonalize").get<int>("NumberOfEigenvalues"));
>
>
>
> //  MyProblem->setHermitian(true);
>
>   Teuchos::RCP<Epetra_MultiVector> initial_eigenvector = Teuchos::rcp(new
> Epetra_MultiVector(Matrix->RowMap(),parameters->sublist("Pure_Diagonalize").get<int>("NumberOfEigenvalues")
> ) );
>
>   Teuchos::RCP<Epetra_Operator> matrix = Teuchos::rcp(Matrix);
>
>   MyProblem->setOperator(matrix);
>
>   MyProblem->setInitVec(initial_eigenvector);
>
>   MyProblem->setProblem();
>
>   bool verbose = (true);
>
>   int verb = Anasazi::TimingDetails;
>
>   Teuchos::ParameterList MyPL;
>
>   MyPL.set( "Verbosity", verb );
>
>   MyPL.set( "Which", "SR" );
>
>   MyPL.set( "Block Size",
> parameters->sublist("Pure_Diagonalize").get<int>("BlockSize"));
>
>   MyPL.set( "Num Blocks",
> parameters->sublist("Pure_Diagonalize").get<int>("NumBlocks"));
>
>
>
>   MyPL.set( "Convergence Tolerance",
> parameters->sublist("Pure_Diagonalize").get<double>("Tolerance") );
>
>
>
>   Anasazi::BlockKrylovSchurSolMgr<double,MV,OP>
> MyBlockKrylovSchur(MyProblem, MyPL );
>
>
>
> Anasazi::ReturnType solverreturn = MyBlockKrylovSchur.solve();
>
>
>
>   switch (solverreturn) {
>
>       case Anasazi::Unconverged:
>
>           if (verbose){
>
>               cout << "Anasazi::BlockKrylovSchur::solve() did not converge!"
> << endl;
>
>           }
>
>           return -1;
>
>           break;
>
>       case Anasazi::Converged:
>
>           if (verbose){
>
>               cout << "Anasazi::BlockKrylovSchur::solve() converged!" <<
> endl;
>
>           }
>
>           break;
>
>   }
>
>   Anasazi::Eigensolution<double, Epetra_MultiVector> solution;
>
>
>
>   solution = MyProblem->getSolution();
>
>   for (int i
> =0;i<parameters->sublist("Pure_Diagonalize").get<int>("NumberOfEigenvalues");i++){
>
>     cout << i <<"th " << solution.Evals[i].realpart <<"+"<<
> solution.Evals[i].imagpart<<"i" <<endl;
>
>   }
>
>   for (int i=0;i<49 ; i++){
>
>     cout << solution.Evecs->operator[](0)[i] <<endl;
>
>   }
>
>
>
> //    int GID = -1; // set GID to the global node to use to compute
>
>                   // the stencil, or put -1 to automatically select it.
>
> //    PrintStencil2D(Matrix, nx, ny, GID);
>
>
>
>     // To created objects must be free'd using delete
>
>   delete Map;
>
> //  delete Matrix;
>
>
>
> #ifdef HAVE_MPI
>
>   MPI_Finalize();
>
> #endif
>
>
>
>   return 0;
>
> }
>
>
>
>
>
>
> _______________________________________________
> Trilinos-Users mailing list
> Trilinos-Users at software.sandia.gov
> http://software.sandia.gov/mailman/listinfo/trilinos-users
>



--
Holger Brandsmeier, SAM, ETH Zürich
http://www.sam.math.ethz.ch/people/bholger


-- 
Holger Brandsmeier, SAM, ETH Zürich
http://www.sam.math.ethz.ch/people/bholger



More information about the Trilinos-Users mailing list