LOBPCG/LOBPCGEpetraExGen.cpp

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

00001 #include "AnasaziConfigDefs.hpp"
00002 #include "AnasaziBasicEigenproblem.hpp"
00003 #include "AnasaziLOBPCGSolMgr.hpp"
00004 #include "AnasaziBasicOutputManager.hpp"
00005 #include "AnasaziEpetraAdapter.hpp"
00006 #include "Epetra_CrsMatrix.h"
00007 #include "Teuchos_CommandLineProcessor.hpp"
00008 
00009 #ifdef HAVE_MPI
00010 #include "Epetra_MpiComm.h"
00011 #include <mpi.h>
00012 #else
00013 #include "Epetra_SerialComm.h"
00014 #endif
00015 #include "Epetra_Map.h"
00016 
00017 #include "ModeLaplace2DQ2.h"
00018 
00019 using namespace Anasazi;
00020 
00021 int main(int argc, char *argv[]) {
00022 
00023 #ifdef HAVE_MPI
00024   // Initialize MPI
00025   //
00026   MPI_Init(&argc,&argv);
00027 #endif
00028 
00029   // Create an Epetra communicator
00030   //
00031 #ifdef HAVE_MPI
00032   Epetra_MpiComm Comm(MPI_COMM_WORLD);
00033 #else
00034   Epetra_SerialComm Comm;
00035 #endif
00036 
00037   // Create an Anasazi output manager
00038   //
00039   BasicOutputManager<double> printer;
00040   printer.stream(Errors) << Anasazi_Version() << endl << endl;
00041 
00042   // Get the sorting string from the command line
00043   //
00044   std::string which("LM");
00045   Teuchos::CommandLineProcessor cmdp(false,true);
00046   cmdp.setOption("sort",&which,"Targetted eigenvalues (SM or LM).");
00047   if (cmdp.parse(argc,argv) != Teuchos::CommandLineProcessor::PARSE_SUCCESSFUL) {
00048 #ifdef HAVE_MPI
00049     MPI_Finalize();
00050 #endif
00051     return -1;
00052   }
00053 
00054   
00055   typedef Epetra_MultiVector MV;
00056   typedef Epetra_Operator OP;
00057   typedef MultiVecTraits<double, Epetra_MultiVector> MVT;
00058 
00059   // Number of dimension of the domain
00060   const int space_dim = 2;
00061 
00062   // Size of each of the dimensions of the domain
00063   std::vector<double> brick_dim( space_dim );
00064   brick_dim[0] = 1.0;
00065   brick_dim[1] = 1.0;
00066 
00067   // Number of elements in each of the dimensions of the domain
00068   std::vector<int> elements( space_dim );
00069   elements[0] = 10;
00070   elements[1] = 10;
00071   
00072   // Create problem
00073   Teuchos::RCP<ModalProblem> testCase = 
00074     Teuchos::rcp( new ModeLaplace2DQ2(Comm, brick_dim[0], elements[0], brick_dim[1], elements[1]) );
00075   
00076   // Get the stiffness and mass matrices
00077   Teuchos::RCP<Epetra_CrsMatrix> K = Teuchos::rcp( const_cast<Epetra_CrsMatrix *>(testCase->getStiffness()), false );
00078   Teuchos::RCP<Epetra_CrsMatrix> M = Teuchos::rcp( const_cast<Epetra_CrsMatrix *>(testCase->getMass()), false );
00079 
00080   // Eigensolver parameters
00081   int nev = 10;
00082   int blockSize = 5;
00083   int maxIters = 500;
00084   double tol = 1.0e-8;
00085 
00086   Teuchos::RCP<Epetra_MultiVector> ivec = Teuchos::rcp( new Epetra_MultiVector(K->OperatorDomainMap(), blockSize) );
00087   ivec->Random();
00088 
00089   // Create the eigenproblem.
00090   Teuchos::RCP<BasicEigenproblem<double, MV, OP> > MyProblem =
00091     Teuchos::rcp( new BasicEigenproblem<double, MV, OP>(K, M, ivec) );
00092 
00093   // Inform the eigenproblem that the operator A is symmetric
00094   MyProblem->setHermitian(true);
00095 
00096   // Set the number of eigenvalues requested
00097   MyProblem->setNEV( nev );
00098 
00099   // Inform the eigenproblem that you are finishing passing it information
00100   bool boolret = MyProblem->setProblem();
00101   if (boolret != true) {
00102     printer.print(Errors,"Anasazi::BasicEigenproblem::setProblem() returned an error.\n");
00103 #ifdef HAVE_MPI
00104     MPI_Finalize();
00105 #endif
00106     return -1;
00107   }
00108 
00109   // Create parameter list to pass into the solver manager
00110   //
00111   Teuchos::ParameterList MyPL;
00112   MyPL.set( "Which", which );
00113   MyPL.set( "Block Size", blockSize );
00114   MyPL.set( "Maximum Iterations", maxIters );
00115   MyPL.set( "Convergence Tolerance", tol );
00116   MyPL.set( "Full Ortho", true );
00117   MyPL.set( "Use Locking", true );
00118   //
00119   // Create the solver manager
00120   LOBPCGSolMgr<double, MV, OP> MySolverMan(MyProblem, MyPL);
00121 
00122   // Solve the problem
00123   //
00124   ReturnType returnCode = MySolverMan.solve();
00125 
00126   // Get the eigenvalues and eigenvectors from the eigenproblem
00127   //
00128   Eigensolution<double,MV> sol = MyProblem->getSolution();
00129   std::vector<Value<double> > evals = sol.Evals;
00130   Teuchos::RCP<MV> evecs = sol.Evecs;
00131 
00132   // Compute residuals.
00133   //
00134   std::vector<double> normR(sol.numVecs);
00135   if (sol.numVecs > 0) {
00136     Teuchos::SerialDenseMatrix<int,double> T(sol.numVecs, sol.numVecs);
00137     Epetra_MultiVector Kvec( K->OperatorDomainMap(), evecs->NumVectors() );
00138     Epetra_MultiVector Mvec( M->OperatorDomainMap(), evecs->NumVectors() );
00139     T.putScalar(0.0); 
00140     for (int i=0; i<sol.numVecs; i++) {
00141       T(i,i) = evals[i].realpart;
00142     }
00143     K->Apply( *evecs, Kvec );  
00144     M->Apply( *evecs, Mvec );  
00145     MVT::MvTimesMatAddMv( -1.0, Mvec, T, 1.0, Kvec );
00146     MVT::MvNorm( Kvec, normR );
00147   }
00148 
00149   // Print the results
00150   //
00151   std::ostringstream os;
00152   os.setf(std::ios_base::right, std::ios_base::adjustfield);
00153   os<<"Solver manager returned " << (returnCode == Converged ? "converged." : "unconverged.") << endl;
00154   os<<endl;
00155   os<<"------------------------------------------------------"<<endl;
00156   os<<std::setw(16)<<"Eigenvalue"
00157     <<std::setw(18)<<"Direct Residual"
00158     <<endl;
00159   os<<"------------------------------------------------------"<<endl;
00160   for (int i=0; i<sol.numVecs; i++) {
00161     os<<std::setw(16)<<evals[i].realpart
00162       <<std::setw(18)<<normR[i]/evals[i].realpart
00163       <<endl;
00164   }
00165   os<<"------------------------------------------------------"<<endl;
00166   printer.print(Errors,os.str());
00167 
00168 #ifdef HAVE_MPI
00169   MPI_Finalize();
00170 #endif
00171   return 0;
00172 }

Generated on Tue Jul 13 09:22:46 2010 for Anasazi by  doxygen 1.4.7