[Trilinos-Users] [EXTERNAL] Anasazi:: sign of EigenVector

[Day, David]

Dear Sunghwan Choi,

If x is an eigenvector of the matrix A,  A x = x lambda, then so is –x, A (-x) = (-x) lambda.  There is no problem with Anasazi.  You might pick some vector y, and always scale the eigenvector so that <y,x> is positive, or you may address this issue in some other way.

--David Day

From: Sunghwan Choi <sunghwanchoi91 at gmail.com<mailto:sunghwanchoi91 at gmail.com>>
Date: Thursday, May 30, 2013 10:14 PM
To: "trilinos-users at software.sandia.gov<mailto:trilinos-users at software.sandia.gov>" <trilinos-users at software.sandia.gov<mailto:trilinos-users at software.sandia.gov>>
Subject: [EXTERNAL] [Trilinos-Users] Anasazi:: sign of EigenVector

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/81802d21/attachment.html