BlockKrylovSchur/BlockKrylovSchurEpetraExGenAmesos.cpp

This is an example of how to use the Anasazi::BlockKrylovSchurSolMgr solver manager to solve a generalized eigenvalue problem, using Epetra data stuctures and the Amesos solver package.

00001 //  This example computes the eigenvalues of smallest magnitude of the 
00002 //  discretized 2D Laplacian operator using the block Krylov-Schur method.  
00003 //  This problem shows the construction of an inner-outer iteration using 
00004 //  Amesos as the linear solver within Anasazi.  This operator is 
00005 //  discretized using linear finite elements and constructed as an Epetra 
00006 //  matrix, then passed into the Amesos solver to perform a factorization.
00007 //  The factorization is then used to apply the shift-invert
00008 //  operation to be used within the Krylov decomposition.  The specifics 
00009 //  of the block Krylov-Schur method can be set by the user.
00010 
00011 // Include autoconfigured header
00012 #include "AnasaziConfigDefs.hpp"
00013 
00014 // Include header for block Krylov-Schur solver
00015 #include "AnasaziBlockKrylovSchurSolMgr.hpp"
00016 
00017 // Include header to define basic eigenproblem Ax = \lambda*Bx
00018 #include "AnasaziBasicEigenproblem.hpp"
00019 
00020 // Include header to provide Anasazi with Epetra adapters
00021 #include "AnasaziEpetraAdapter.hpp"
00022 
00023 // Include header for Epetra compressed-row storage matrix and linear problem
00024 #include "Epetra_CrsMatrix.h"
00025 #include "Epetra_LinearProblem.h"
00026 
00027 // Include header for Amesos solver and solver interface for Epetra_Operator
00028 #include "Amesos.h"
00029 
00030 // Include header for Teuchos serial dense matrix
00031 #include "Teuchos_SerialDenseMatrix.hpp"
00032 
00033 // Include header for the problem definition
00034 #include "ModeLaplace2DQ2.h"
00035 
00036 // Include selected communicator class and map required by Epetra objects
00037 #ifdef EPETRA_MPI
00038 #include "Epetra_MpiComm.h"
00039 #else
00040 #include "Epetra_SerialComm.h"
00041 #endif
00042 #include "Epetra_Map.h"
00043 
00044 // ****************************************************************************
00045 // Class:  AmesosGenOp
00046 // Purpose:  Applies A^{-1}*B*X = Y or A^{-T}*B*X = Y where A is an Amesos solver
00047 //           and B is an Epetra_Operator.  Class supports Epetra_Operator interface.
00048 // Date: Jan 9, 2006
00049 // Author:  Heidi K. Thornquist
00050 // ****************************************************************************
00051 
00052 class AmesosGenOp : public virtual Epetra_Operator
00053 {
00054 public:
00055   // Basic constructor
00056   AmesosGenOp( Epetra_LinearProblem& problem,
00057                const Teuchos::RCP<Amesos_BaseSolver>& solver,
00058                const Teuchos::RCP<Epetra_Operator>& massMtx,
00059                bool useTranspose = false );
00060   // Destructor
00061   ~AmesosGenOp() {};
00062   
00063   // Methods for supporting Epetra_Operator interface
00064   int Apply(const Epetra_MultiVector& X, Epetra_MultiVector& Y ) const;
00065   const char* Label() const { return "Amesos direct solver for applying A^{-1}M"; };
00066   bool UseTranspose() const { return useTranspose_; };
00067   int SetUseTranspose( bool useTranspose );
00068   const Epetra_Comm& Comm() const { return solver_->Comm(); };
00069   const Epetra_Map& OperatorDomainMap() const { return massMtx_->OperatorDomainMap(); };
00070   const Epetra_Map& OperatorRangeMap() const { return massMtx_->OperatorRangeMap(); };
00071     
00072   // Epetra_Operator interface methods that are not supported.
00073   // Note:  ApplyInverse not defined because M not guaranteed to have an inverse.
00074   int ApplyInverse(const Epetra_MultiVector& X, Epetra_MultiVector& Y ) const { return -1; };
00075   bool HasNormInf() const { return false; };
00076   double NormInf() const { return -1.0; };
00077    
00078 private:  
00079   // Default constructor
00080   AmesosGenOp () {};
00081  
00082   // Copy constructor
00083   AmesosGenOp ( const AmesosGenOp& genOp ) {};
00084 
00085   // Epetra_LinearProblem contained in the Amesos_BaseSolver
00086   bool useTranspose_;
00087   Teuchos::RCP<Amesos_BaseSolver> solver_;
00088   Teuchos::RCP<Epetra_Operator> massMtx_;
00089   Epetra_LinearProblem* problem_;
00090   
00091 };
00092 
00093 // ****************************************************************************
00094 // BEGIN MAIN ROUTINE
00095 // ****************************************************************************
00096 
00097 int main(int argc, char *argv[]) {
00098   int i;
00099 
00100 #ifdef EPETRA_MPI
00101   // Initialize MPI
00102   MPI_Init(&argc,&argv);
00103   Epetra_MpiComm Comm(MPI_COMM_WORLD);
00104 #else
00105   Epetra_SerialComm Comm;
00106 #endif
00107 
00108   int MyPID = Comm.MyPID();
00109 
00110   // Number of dimension of the domain
00111   int space_dim = 2;
00112   
00113   // Size of each of the dimensions of the domain
00114   std::vector<double> brick_dim( space_dim );
00115   brick_dim[0] = 1.0;
00116   brick_dim[1] = 1.0;
00117   
00118   // Number of elements in each of the dimensions of the domain
00119   std::vector<int> elements( space_dim );
00120   elements[0] = 10;
00121   elements[1] = 10;
00122   
00123   // Create problem
00124   Teuchos::RCP<ModalProblem> testCase = Teuchos::rcp( new ModeLaplace2DQ2(Comm, brick_dim[0], elements[0], brick_dim[1], elements[1]) );
00125   
00126   // Get the stiffness and mass matrices
00127   Teuchos::RCP<Epetra_CrsMatrix> K = Teuchos::rcp( const_cast<Epetra_CrsMatrix *>(testCase->getStiffness()), false );
00128   Teuchos::RCP<Epetra_CrsMatrix> M = Teuchos::rcp( const_cast<Epetra_CrsMatrix *>(testCase->getMass()), false );
00129 
00130   //
00131   // *******************************************************
00132   // Set up Amesos direct solver for inner iteration
00133   // *******************************************************
00134   //
00135 
00136   // Create Epetra linear problem class to solve "Kx = b"
00137   Epetra_LinearProblem AmesosProblem;
00138   AmesosProblem.SetOperator(K.get());
00139   
00140   // Create Amesos factory and solver for solving "Kx = b" using a direct factorization
00141   Amesos amesosFactory;
00142   Teuchos::RCP<Amesos_BaseSolver> AmesosSolver = 
00143     Teuchos::rcp( amesosFactory.Create( "Klu", AmesosProblem ) );
00144 
00145   // The AmesosGenOp class assumes that the symbolic/numeric factorizations have already
00146   // been performed on the linear problem.
00147   AmesosSolver->SymbolicFactorization();
00148   AmesosSolver->NumericFactorization();
00149 
00150   //
00151   // ************************************
00152   // Start the block Arnoldi iteration
00153   // ************************************
00154   //
00155   //  Variables used for the Block Arnoldi Method
00156   //
00157   int nev = 10;
00158   int blockSize = 3;  
00159   int numBlocks = 3*nev/blockSize;
00160   int maxRestarts = 5;
00161   //int step = 5;
00162   double tol = 1.0e-8;
00163   std::string which = "LM";
00164   int verbosity = Anasazi::Errors + Anasazi::Warnings + Anasazi::FinalSummary;
00165   //
00166   // Create parameter list to pass into solver
00167   //
00168   Teuchos::ParameterList MyPL;
00169   MyPL.set( "Verbosity", verbosity );
00170   MyPL.set( "Which", which );
00171   MyPL.set( "Block Size", blockSize );
00172   MyPL.set( "Num Blocks", numBlocks );
00173   MyPL.set( "Maximum Restarts", maxRestarts );
00174   MyPL.set( "Convergence Tolerance", tol );
00175   //MyPL.set( "Step Size", step );
00176   
00177   typedef Epetra_MultiVector MV;
00178   typedef Epetra_Operator OP;
00179   typedef Anasazi::MultiVecTraits<double, MV> MVT;
00180   typedef Anasazi::OperatorTraits<double, MV, OP> OPT;
00181   
00182   // Create an Epetra_MultiVector for an initial vector to start the solver.
00183   // Note:  This needs to have the same number of columns as the blocksize.
00184   Teuchos::RCP<Epetra_MultiVector> ivec = Teuchos::rcp( new Epetra_MultiVector(K->Map(), blockSize) );
00185   MVT::MvRandom( *ivec );
00186   
00187   // Create the Epetra_Operator for the spectral transformation using the Amesos direct solver.
00188   Teuchos::RCP<AmesosGenOp> Aop 
00189     = Teuchos::rcp( new AmesosGenOp(AmesosProblem, AmesosSolver, M) );
00190   
00191   Teuchos::RCP<Anasazi::BasicEigenproblem<double,MV,OP> > MyProblem = 
00192     Teuchos::rcp( new Anasazi::BasicEigenproblem<double,MV,OP>(Aop, M, ivec) );
00193   
00194   // Inform the eigenproblem that the matrix pencil (K,M) is symmetric
00195   MyProblem->setHermitian(true);
00196   
00197   // Set the number of eigenvalues requested 
00198   MyProblem->setNEV( nev );
00199   
00200   // Inform the eigenproblem that you are finished passing it information
00201   bool boolret = MyProblem->setProblem();
00202   if (boolret != true) {
00203     if (MyPID == 0) {
00204       cout << "Anasazi::BasicEigenproblem::setProblem() returned with error." << endl;
00205     }
00206 #ifdef HAVE_MPI
00207     MPI_Finalize() ;
00208 #endif
00209     return -1;
00210   }
00211 
00212   // Initialize the Block Arnoldi solver
00213   Anasazi::BlockKrylovSchurSolMgr<double, MV, OP> MySolverMgr(MyProblem, MyPL);
00214   
00215   // Solve the problem to the specified tolerances or length
00216   Anasazi::ReturnType returnCode = MySolverMgr.solve();
00217   if (returnCode != Anasazi::Converged && MyPID==0) {
00218     cout << "Anasazi::EigensolverMgr::solve() returned unconverged." << endl;
00219   }
00220   
00221   // Get the eigenvalues and eigenvectors from the eigenproblem
00222   Anasazi::Eigensolution<double,MV> sol = MyProblem->getSolution();
00223   std::vector<Anasazi::Value<double> > evals = sol.Evals;
00224   Teuchos::RCP<MV> evecs = sol.Evecs;
00225   int numev = sol.numVecs;
00226   
00227   if (numev > 0) {
00228     
00229     Teuchos::SerialDenseMatrix<int,double> dmatr(numev,numev);
00230     Epetra_MultiVector tempvec(K->Map(), MVT::GetNumberVecs( *evecs ));
00231     OPT::Apply( *K, *evecs, tempvec );
00232     MVT::MvTransMv( 1.0, tempvec, *evecs, dmatr );
00233     
00234     if (MyPID==0) {
00235       double compeval = 0.0;
00236       cout.setf(std::ios_base::right, std::ios_base::adjustfield);
00237       cout<<"Actual Eigenvalues (obtained by Rayleigh quotient) : "<<endl;
00238       cout<<"------------------------------------------------------"<<endl;
00239       cout<<std::setw(16)<<"Real Part"
00240         <<std::setw(16)<<"Rayleigh Error"<<endl;
00241       cout<<"------------------------------------------------------"<<endl;
00242       for (i=0; i<numev; i++) {
00243         compeval = dmatr(i,i);
00244         cout<<std::setw(16)<<compeval
00245           <<std::setw(16)<<Teuchos::ScalarTraits<double>::magnitude(compeval-1.0/evals[i].realpart)
00246           <<endl;
00247       }
00248       cout<<"------------------------------------------------------"<<endl;
00249     }
00250     
00251   }
00252 
00253 #ifdef EPETRA_MPI
00254   MPI_Finalize();
00255 #endif
00256 
00257   return 0;
00258 }
00259 
00260 // ****************************************************************************
00261 // Class:  AmesosGenOp
00262 // Purpose:  Applies A^{-1}*B*X = Y or A^{-T}*B*X = Y where A is an Amesos solver
00263 //           and B is an Epetra_Operator.  Class supports Epetra_Operator interface.
00264 // Date: Jan 9, 2006
00265 // Author:  Heidi K. Thornquist
00266 // ****************************************************************************
00267 
00268 
00269 AmesosGenOp::AmesosGenOp( Epetra_LinearProblem& problem,
00270                           const Teuchos::RCP<Amesos_BaseSolver>& solver,
00271                           const Teuchos::RCP<Epetra_Operator>& massMtx,
00272                           bool useTranspose )
00273   : useTranspose_(useTranspose),
00274     solver_(solver),
00275     massMtx_(massMtx)
00276 {
00277   problem_ = const_cast<Epetra_LinearProblem*>( solver->GetProblem() );
00278   
00279   if ( solver_->UseTranspose() )
00280     solver_->SetUseTranspose(!useTranspose); 
00281   else
00282     solver_->SetUseTranspose(useTranspose);
00283   
00284   if ( massMtx_->UseTranspose() )
00285     massMtx_->SetUseTranspose(!useTranspose);
00286   else
00287     massMtx_->SetUseTranspose(useTranspose);    
00288 }
00289 
00290 // Methods for supporting Epetra_Operator interface
00291 int AmesosGenOp::SetUseTranspose(bool useTranspose) 
00292 { 
00293   useTranspose_ = useTranspose; 
00294   if ( solver_->UseTranspose() )
00295     solver_->SetUseTranspose(!useTranspose); 
00296   else
00297     solver_->SetUseTranspose(useTranspose);
00298   
00299   if ( massMtx_->UseTranspose() )
00300     massMtx_->SetUseTranspose(!useTranspose);
00301   else
00302     massMtx_->SetUseTranspose(useTranspose);
00303 
00304   return 0;
00305 }
00306 
00307 int AmesosGenOp::Apply(const Epetra_MultiVector& X, Epetra_MultiVector& Y ) const 
00308 {
00309   if (!useTranspose_) {
00310     
00311     // Storage for M*X
00312     Epetra_MultiVector MX(X.Map(),X.NumVectors());
00313     
00314     // Apply M*X
00315     massMtx_->Apply(X, MX);
00316     Y.PutScalar(0.0);
00317     
00318     // Set the LHS and RHS
00319     problem_->SetRHS(&MX);
00320     problem_->SetLHS(&Y);
00321 
00322     // Solve the linear system A*Y = MX
00323     solver_->Solve();
00324   }
00325   else {
00326     // Storage for A^{-T}*X
00327     Epetra_MultiVector ATX(X.Map(),X.NumVectors());
00328     Epetra_MultiVector tmpX = const_cast<Epetra_MultiVector&>(X);
00329     
00330     // Set the LHS and RHS
00331     problem_->SetRHS(&tmpX);
00332     problem_->SetLHS(&ATX);
00333     
00334     // Solve the linear system A^T*Y = X 
00335     solver_->Solve();
00336     
00337     // Apply M*ATX
00338     massMtx_->Apply(ATX, Y);
00339   }
00340   
00341   return 0;
00342 }

Generated on Wed May 12 21:40:22 2010 for Anasazi by  doxygen 1.4.7