BlockKrylovSchur/BlockKrylovSchurEpetraExGenAztecOO.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 AztecOO 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 //  AztecOO 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 AztecOO 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 manager
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 AztecOO solver and solver interface for Epetra_Operator
00028 #include "AztecOO.h"
00029 #include "AztecOO_Operator.h"
00030 
00031 // Include header for Ifpack incomplete Cholesky preconditioner
00032 #include "Ifpack_CrsIct.h"
00033 
00034 // Include header for Teuchos serial dense matrix
00035 #include "Teuchos_SerialDenseMatrix.hpp"
00036 
00037 // Include header for the problem definition
00038 #include "ModeLaplace2DQ2.h"
00039 
00040 // Include selected communicator class and map required by Epetra objects
00041 #ifdef EPETRA_MPI
00042 #include "Epetra_MpiComm.h"
00043 #else
00044 #include "Epetra_SerialComm.h"
00045 #endif
00046 #include "Epetra_Map.h"
00047 
00048 int main(int argc, char *argv[]) {
00049   int i, info;
00050 
00051 #ifdef EPETRA_MPI
00052   // Initialize MPI
00053   MPI_Init(&argc,&argv);
00054   Epetra_MpiComm Comm(MPI_COMM_WORLD);
00055 #else
00056   Epetra_SerialComm Comm;
00057 #endif
00058 
00059   int MyPID = Comm.MyPID();
00060 
00061   // Number of dimension of the domain
00062   int space_dim = 2;
00063   
00064   // Size of each of the dimensions of the domain
00065   std::vector<double> brick_dim( space_dim );
00066   brick_dim[0] = 1.0;
00067   brick_dim[1] = 1.0;
00068   
00069   // Number of elements in each of the dimensions of the domain
00070   std::vector<int> elements( space_dim );
00071   elements[0] = 10;
00072   elements[1] = 10;
00073   
00074   // Create problem
00075   Teuchos::RCP<ModalProblem> testCase = Teuchos::rcp( new ModeLaplace2DQ2(Comm, brick_dim[0], elements[0], brick_dim[1], elements[1]) );
00076   
00077   // Get the stiffness and mass matrices
00078   Teuchos::RCP<Epetra_CrsMatrix> K = Teuchos::rcp( const_cast<Epetra_CrsMatrix *>(testCase->getStiffness()), false );
00079   Teuchos::RCP<Epetra_CrsMatrix> M = Teuchos::rcp( const_cast<Epetra_CrsMatrix *>(testCase->getMass()), false );
00080   
00081   //
00082   //*****Select the Preconditioner*****
00083   //
00084   if (MyPID==0) cout << endl << endl;
00085   if (MyPID==0) cout << "Constructing ICT preconditioner" << endl;
00086   int Lfill = 0;
00087   if (argc > 1) Lfill = atoi(argv[1]);
00088   if (MyPID==0) cout << "Using Lfill = " << Lfill << endl;
00089   int Overlap = 0;
00090   if (argc > 2) Overlap = atoi(argv[2]);
00091   if (MyPID==0) cout << "Using Level Overlap = " << Overlap << endl;
00092   double Athresh = 0.0;
00093   if (argc > 3) Athresh = atof(argv[3]);
00094   if (MyPID==0) cout << "Using Absolute Threshold Value of " << Athresh << endl;
00095   double Rthresh = 1.0;
00096   if (argc >4) Rthresh = atof(argv[4]);
00097   if (MyPID==0) cout << "Using Relative Threshold Value of " << Rthresh << endl;
00098   double dropTol = 1.0e-6;
00099   //
00100   Teuchos::RCP<Ifpack_CrsIct> ICT;
00101   //
00102   if (Lfill > -1) {
00103     ICT = Teuchos::rcp( new Ifpack_CrsIct(*K, dropTol, Lfill) );
00104     ICT->SetAbsoluteThreshold(Athresh);
00105     ICT->SetRelativeThreshold(Rthresh);
00106     int initerr = ICT->InitValues(*K);
00107     if (initerr != 0) cout << "InitValues error = " << initerr;
00108     info = ICT->Factor();
00109     assert( info==0 );
00110   } 
00111   //
00112   bool transA = false;
00113   double Cond_Est;
00114   ICT->Condest(transA, Cond_Est);
00115   if (MyPID==0) {
00116     cout << "Condition number estimate for this preconditoner = " << Cond_Est << endl;
00117     cout << endl;
00118   } 
00119   //
00120   // *******************************************************
00121   // Set up AztecOO CG operator for inner iteration
00122   // *******************************************************
00123   //
00124   // Create Epetra linear problem class to solve "Ax = b"
00125   Epetra_LinearProblem precProblem;
00126   precProblem.SetOperator(K.get());
00127   
00128   // Create AztecOO solver for solving "Ax = b" using an incomplete cholesky preconditioner
00129   AztecOO precSolver(precProblem);
00130   precSolver.SetPrecOperator(ICT.get());
00131   precSolver.SetAztecOption(AZ_output, AZ_none);
00132   precSolver.SetAztecOption(AZ_solver, AZ_cg);
00133   
00134   // Use AztecOO solver to create the AztecOO_Operator
00135   Teuchos::RCP<AztecOO_Operator> precOperator =
00136     Teuchos::rcp( new AztecOO_Operator(&precSolver, K->NumGlobalRows(), 1e-12) ); 
00137   //
00138   // ************************************
00139   // Start the block Arnoldi iteration
00140   // ************************************
00141   //
00142   //  Variables used for the Block Arnoldi Method
00143   //
00144   double tol = 1.0e-8;
00145   int nev = 10;
00146   int blockSize = 3;  
00147   int numBlocks = 3*nev/blockSize;
00148   int maxRestarts = 10;
00149   //int step = 5;
00150   std::string which = "LM";
00151   int verbosity = Anasazi::Errors + Anasazi::Warnings + Anasazi::FinalSummary;
00152   //
00153   // Create parameter list to pass into solver
00154   //
00155   Teuchos::ParameterList MyPL;
00156   MyPL.set( "Verbosity", verbosity );
00157   MyPL.set( "Which", which );
00158   MyPL.set( "Block Size", blockSize );
00159   MyPL.set( "Num Blocks", numBlocks );
00160   MyPL.set( "Maximum Restarts", maxRestarts );
00161   MyPL.set( "Convergence Tolerance", tol );
00162   //MyPL.set( "Step Size", step );
00163   
00164   typedef Epetra_MultiVector MV;
00165   typedef Epetra_Operator OP;
00166   typedef Anasazi::MultiVecTraits<double, MV> MVT;
00167   typedef Anasazi::OperatorTraits<double, MV, OP> OPT;
00168   
00169   // Create an Epetra_MultiVector for an initial vector to start the solver.
00170   // Note:  This needs to have the same number of columns as the blocksize.
00171   Teuchos::RCP<Epetra_MultiVector> ivec = Teuchos::rcp( new Epetra_MultiVector(K->Map(), blockSize) );
00172   MVT::MvRandom( *ivec );
00173   
00174   // Call the ctor that calls the petra ctor for a matrix
00175   Teuchos::RCP<Anasazi::EpetraGenOp> Aop = Teuchos::rcp( new Anasazi::EpetraGenOp(precOperator, M) );
00176   
00177   Teuchos::RCP<Anasazi::BasicEigenproblem<double,MV,OP> > MyProblem = 
00178     Teuchos::rcp( new Anasazi::BasicEigenproblem<double,MV,OP>(Aop, M, ivec) );
00179   
00180   // Inform the eigenproblem that the matrix pencil (K,M) is symmetric
00181   MyProblem->setHermitian(true);
00182   
00183   // Set the number of eigenvalues requested 
00184   MyProblem->setNEV( nev );
00185 
00186   // Inform the eigenproblem that you are finished passing it information
00187   bool boolret = MyProblem->setProblem();
00188   if (boolret != true) {
00189     if (MyPID == 0) {
00190       cout << "Anasazi::BasicEigenproblem::setProblem() returned with error." << endl;
00191     }
00192 #ifdef HAVE_MPI
00193     MPI_Finalize() ;
00194 #endif
00195     return -1;
00196   }
00197 
00198   // Initialize the Block Arnoldi solver
00199   Anasazi::BlockKrylovSchurSolMgr<double, MV, OP> MySolverMgr(MyProblem, MyPL);
00200 
00201   // Solve the problem to the specified tolerances or length
00202   Anasazi::ReturnType returnCode = MySolverMgr.solve();
00203   if (returnCode != Anasazi::Converged && MyPID==0) {
00204     cout << "Anasazi::EigensolverMgr::solve() returned unconverged." << endl;
00205   }
00206 
00207   // Get the eigenvalues and eigenvectors from the eigenproblem
00208   Anasazi::Eigensolution<double,MV> sol = MyProblem->getSolution();
00209   std::vector<Anasazi::Value<double> > evals = sol.Evals;
00210   Teuchos::RCP<MV> evecs = sol.Evecs;
00211   int numev = sol.numVecs;
00212 
00213   if (numev > 0) {
00214 
00215     Teuchos::SerialDenseMatrix<int,double> dmatr(numev,numev);
00216     Epetra_MultiVector tempvec(K->Map(), MVT::GetNumberVecs( *evecs ));
00217     OPT::Apply( *K, *evecs, tempvec );
00218     MVT::MvTransMv( 1.0, tempvec, *evecs, dmatr );
00219     
00220     if (MyPID==0) {
00221       double compeval = 0.0;
00222       cout.setf(std::ios_base::right, std::ios_base::adjustfield);
00223       cout<<"Actual Eigenvalues (obtained by Rayleigh quotient) : "<<endl;
00224       cout<<"------------------------------------------------------"<<endl;
00225       cout<<std::setw(16)<<"Real Part"
00226         <<std::setw(16)<<"Rayleigh Error"<<endl;
00227       cout<<"------------------------------------------------------"<<endl;
00228       for (i=0; i<numev; i++) {
00229         compeval = dmatr(i,i);
00230         cout<<std::setw(16)<<compeval
00231           <<std::setw(16)<<Teuchos::ScalarTraits<double>::magnitude(compeval-1.0/evals[i].realpart)
00232           <<endl;
00233       }
00234       cout<<"------------------------------------------------------"<<endl;
00235     }
00236     
00237   }
00238   
00239 #ifdef EPETRA_MPI
00240   MPI_Finalize();
00241 #endif
00242 
00243   return 0;
00244 }

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