[Trilinos-Users] Anasazi:: sign of EigenVector

[Sunghwan Choi]

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;

}

 

 

-------------- next part --------------
An HTML attachment was scrubbed...
URL: https://software.sandia.gov/pipermail/trilinos-users/attachments/20130531/d6d83a3c/attachment.html