BlockDavidson/BlockDavidsonEpetraExGen.cpp

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

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

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