BlockKrylovSchur/BlockKrylovSchurEpetraExGenBelos.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 Belos 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 //  Belos as the linear solver within Anasazi.  An Ifpack preconditioner 
00005 //  is constructed to precondition the linear solver.  This operator is 
00006 //  discretized using linear finite elements and constructed as an Epetra 
00007 //  matrix, then passed into the Belos solver to perform 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 #include "Epetra_InvOperator.h"
00027 
00028 // Include header for Belos solver and solver interface for Epetra_Operator
00029 #include "BelosEpetraOperator.h"
00030 #include "BelosEpetraAdapter.hpp"
00031 
00032 // Include header for Ifpack incomplete Cholesky preconditioner
00033 #include "Ifpack_CrsIct.h"
00034 
00035 // Include header for Teuchos serial dense matrix
00036 #include "Teuchos_SerialDenseMatrix.hpp"
00037 
00038 // Include header for the problem definition
00039 #include "ModeLaplace2DQ2.h"
00040 
00041 // Include selected communicator class and map required by Epetra objects
00042 #ifdef EPETRA_MPI
00043 #include "Epetra_MpiComm.h"
00044 #else
00045 #include "Epetra_SerialComm.h"
00046 #endif
00047 #include "Epetra_Map.h"
00048 
00049 int main(int argc, char *argv[]) {
00050   int i, info;
00051 
00052 #ifdef EPETRA_MPI
00053   // Initialize MPI
00054   MPI_Init(&argc,&argv);
00055   Epetra_MpiComm Comm(MPI_COMM_WORLD);
00056 #else
00057   Epetra_SerialComm Comm;
00058 #endif
00059 
00060   int MyPID = Comm.MyPID();
00061 
00062   // Number of dimension of the domain
00063   int space_dim = 2;
00064   
00065   // Size of each of the dimensions of the domain
00066   std::vector<double> brick_dim( space_dim );
00067   brick_dim[0] = 1.0;
00068   brick_dim[1] = 1.0;
00069   
00070   // Number of elements in each of the dimensions of the domain
00071   std::vector<int> elements( space_dim );
00072   elements[0] = 10;
00073   elements[1] = 10;
00074   
00075   // Create problem
00076   Teuchos::RCP<ModalProblem> testCase = Teuchos::rcp( new ModeLaplace2DQ2(Comm, brick_dim[0], elements[0], brick_dim[1], elements[1]) );
00077   
00078   // Get the stiffness and mass matrices
00079   Teuchos::RCP<Epetra_CrsMatrix> K = Teuchos::rcp( const_cast<Epetra_CrsMatrix *>(testCase->getStiffness()), false );
00080   Teuchos::RCP<Epetra_CrsMatrix> M = Teuchos::rcp( const_cast<Epetra_CrsMatrix *>(testCase->getMass()), false );
00081   
00082   //
00083   //*****Select the Preconditioner*****
00084   //
00085   if (MyPID==0) cout << endl << endl;
00086   if (MyPID==0) cout << "Constructing ICT preconditioner" << endl;
00087   int Lfill = 0;
00088   if (argc > 1) Lfill = atoi(argv[1]);
00089   if (MyPID==0) cout << "Using Lfill = " << Lfill << endl;
00090   int Overlap = 0;
00091   if (argc > 2) Overlap = atoi(argv[2]);
00092   if (MyPID==0) cout << "Using Level Overlap = " << Overlap << endl;
00093   double Athresh = 0.0;
00094   if (argc > 3) Athresh = atof(argv[3]);
00095   if (MyPID==0) cout << "Using Absolute Threshold Value of " << Athresh << endl;
00096   double Rthresh = 1.0;
00097   if (argc >4) Rthresh = atof(argv[4]);
00098   if (MyPID==0) cout << "Using Relative Threshold Value of " << Rthresh << endl;
00099   double dropTol = 1.0e-6;
00100   //
00101   Teuchos::RCP<Ifpack_CrsIct> ICT;
00102   //
00103   if (Lfill > -1) {
00104     ICT = Teuchos::rcp( new Ifpack_CrsIct(*K, dropTol, Lfill) );
00105     ICT->SetAbsoluteThreshold(Athresh);
00106     ICT->SetRelativeThreshold(Rthresh);
00107     int initerr = ICT->InitValues(*K);
00108     if (initerr != 0) cout << "InitValues error = " << initerr;
00109     info = ICT->Factor();
00110     assert( info==0 );
00111   } 
00112   //
00113   bool transA = false;
00114   double Cond_Est;
00115   ICT->Condest(transA, Cond_Est);
00116   if (MyPID==0) {
00117     cout << "Condition number estimate for this preconditoner = " << Cond_Est << endl;
00118     cout << endl;
00119   } 
00120   //
00121   //*******************************************************/
00122   // Set up Belos Block GMRES operator for inner iteration
00123   //*******************************************************/
00124   //
00125   int blockSize = 3; // block size used by linear solver and eigensolver [ not required to be the same ]
00126   int maxits = K->NumGlobalRows(); // maximum number of iterations to run
00127   //
00128   // Create the Belos::LinearProblem
00129   //
00130   Teuchos::RCP<Belos::LinearProblem<double,Epetra_MultiVector,Epetra_Operator> > 
00131     My_LP = Teuchos::rcp( new Belos::LinearProblem<double,Epetra_MultiVector,Epetra_Operator>() );
00132   My_LP->setOperator( K );
00133 
00134   // Create the Belos preconditioned operator from the Ifpack preconditioner.
00135   // NOTE:  This is necessary because Belos expects an operator to apply the 
00136   //        preconditioner with Apply() NOT ApplyInverse().
00137   Teuchos::RCP<Epetra_Operator> belosPrec = Teuchos::rcp( new Epetra_InvOperator( &*ICT ) );
00138   My_LP->setLeftPrec( belosPrec );
00139   //
00140   // Create the ParameterList for the Belos Operator
00141   // 
00142   Teuchos::RCP<Teuchos::ParameterList> My_List = Teuchos::rcp( new Teuchos::ParameterList() );
00143   My_List->set( "Solver", "BlockCG" );
00144   My_List->set( "Maximum Iterations", maxits );
00145   My_List->set( "Block Size", 1 );
00146   My_List->set( "Convergence Tolerance", 1e-12 );
00147   //
00148   // Create the Belos::EpetraOperator
00149   //
00150   Teuchos::RCP<Belos::EpetraOperator> BelosOp = 
00151     Teuchos::rcp( new Belos::EpetraOperator( My_LP, My_List ));
00152   //
00153   // ************************************
00154   // Start the block Arnoldi iteration
00155   // ************************************
00156   //
00157   //  Variables used for the Block Arnoldi Method
00158   //
00159   double tol = 1.0e-8;
00160   int nev = 10;
00161   int numBlocks = 3*nev/blockSize;
00162   int maxRestarts = 5;
00163   //int step = 5;
00164   std::string which = "LM";
00165   int verbosity = Anasazi::Errors + Anasazi::Warnings + Anasazi::FinalSummary;
00166   //
00167   // Create parameter list to pass into solver
00168   //
00169   Teuchos::ParameterList MyPL;
00170   MyPL.set( "Verbosity", verbosity );
00171   MyPL.set( "Which", which );
00172   MyPL.set( "Block Size", blockSize );
00173   MyPL.set( "Num Blocks", numBlocks );
00174   MyPL.set( "Maximum Restarts", maxRestarts );
00175   MyPL.set( "Convergence Tolerance", tol );
00176   //MyPL.set( "Step Size", step );
00177   
00178   typedef Epetra_MultiVector MV;
00179   typedef Epetra_Operator OP;
00180   typedef Anasazi::MultiVecTraits<double, MV> MVT;
00181   typedef Anasazi::OperatorTraits<double, MV, OP> OPT;
00182   
00183   // Create an Epetra_MultiVector for an initial vector to start the solver.
00184   // Note:  This needs to have the same number of columns as the blocksize.
00185   Teuchos::RCP<Epetra_MultiVector> ivec = Teuchos::rcp( new Epetra_MultiVector(K->Map(), blockSize) );
00186   MVT::MvRandom( *ivec );
00187   
00188   // Call the ctor that calls the petra ctor for a matrix
00189   Teuchos::RCP<Anasazi::EpetraGenOp> Aop = Teuchos::rcp( new Anasazi::EpetraGenOp(BelosOp, M, false) );
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 }

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