Anasazi: LOBPCG/LOBPCGEpetraExGenPrecIfpack.cpp

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 "Epetra_InvOperator.h"
00008 #include "Teuchos_CommandLineProcessor.hpp"
00009 
00010 // Include header for Ifpack incomplete Cholesky preconditioner
00011 #include "Ifpack.h"
00012 
00013 #ifdef HAVE_MPI
00014 #include "Epetra_MpiComm.h"
00015 #include <mpi.h>
00016 #else
00017 #include "Epetra_SerialComm.h"
00018 #endif
00019 #include "Epetra_Map.h"
00020 
00021 #include "ModeLaplace2DQ2.h"
00022 
00023 
00024 using namespace Anasazi;
00025 
00026 //*****************************************************************************
00027 // Begin main routine
00028 //*****************************************************************************
00029 int main(int argc, char *argv[]) {
00030 
00031 #ifdef HAVE_MPI
00032   // Initialize MPI
00033   MPI_Init(&argc,&argv);
00034 #endif
00035 
00036   // Create an Epetra communicator
00037 #ifdef HAVE_MPI
00038   Epetra_MpiComm Comm(MPI_COMM_WORLD);
00039 #else
00040   Epetra_SerialComm Comm;
00041 #endif
00042 
00043   //************************************
00044   // Get the parameters from the command line
00045   //************************************
00046   //
00047   int    nev       = 10;
00048   int    blockSize = 5;
00049   int    maxIters  = 500;
00050   double tol       = 1.0e-8;
00051   int numElements = 10;
00052   bool verbose = false;
00053   std::string which("LM");
00054   bool usePrec = true;
00055   double prec_dropTol = 1e-4;
00056   int prec_lofill = 0;
00057   Teuchos::CommandLineProcessor cmdp(false,true);
00058   cmdp.setOption("nev",&nev,"Number of eigenpairs to compted.");
00059   cmdp.setOption("maxIters",&maxIters,"Maximum number of iterations.");
00060   cmdp.setOption("blockSize",&blockSize,"Block size.");
00061   cmdp.setOption("tol",&tol,"Relative convergence tolerance.");
00062   cmdp.setOption("numElements",&numElements,"Number of elements in the discretization.");
00063   cmdp.setOption("verbose","quiet",&verbose,"Print messages and results.");
00064   cmdp.setOption("sort",&which,"Targetted eigenvalues (SM or LM).");
00065   cmdp.setOption("usePrec","noPrec",&usePrec,"Use Ifpack for preconditioning.");
00066   cmdp.setOption("prec_dropTol",&prec_dropTol,"Preconditioner: drop tolerance.");
00067   cmdp.setOption("prec_lofill",&prec_lofill,"Preconditioner: level of fill.");
00068   if (cmdp.parse(argc,argv) != Teuchos::CommandLineProcessor::PARSE_SUCCESSFUL) {
00069 #ifdef HAVE_MPI
00070     MPI_Finalize();
00071 #endif
00072     return -1;
00073   }
00074 
00075   //************************************
00076   // Create an Anasazi output manager
00077   //************************************
00078   //
00079   // Set verbosity level
00080   int verbosity = Anasazi::Errors + Anasazi::Warnings;
00081   if (verbose) {
00082     verbosity += Anasazi::FinalSummary + Anasazi::TimingDetails;
00083   }
00084   BasicOutputManager<double> printer(verbosity);
00085   printer.stream(Errors) << Anasazi_Version() << endl << endl;
00086 
00087   //************************************
00088   // Some useful typedefs
00089   //************************************
00090   //
00091   typedef Epetra_MultiVector MV;
00092   typedef Epetra_Operator OP;
00093   typedef MultiVecTraits<double, Epetra_MultiVector> MVT;
00094 
00095   //************************************
00096   // Create the problem matrices
00097   //************************************
00098   //
00099   printer.stream(Errors) << "Generating problem matrices..." << flush;
00100   // Number of dimension of the domain
00101   const int space_dim = 2;
00102   // Size of each of the dimensions of the domain
00103   std::vector<double> brick_dim( space_dim );
00104   brick_dim[0] = 1.0;
00105   brick_dim[1] = 1.0;
00106   // Number of elements in each of the dimensions of the domain
00107   std::vector<int> elements( space_dim );
00108   elements[0] = numElements;
00109   elements[1] = numElements;
00110   // Create problem
00111   Teuchos::RCP<ModalProblem> testCase = 
00112     Teuchos::rcp( new ModeLaplace2DQ2(Comm, brick_dim[0], elements[0], brick_dim[1], elements[1]) );
00113   // Get the stiffness and mass matrices
00114   Teuchos::RCP<Epetra_CrsMatrix> K = Teuchos::rcp( const_cast<Epetra_CrsMatrix *>(testCase->getStiffness()), false );
00115   Teuchos::RCP<Epetra_CrsMatrix> M = Teuchos::rcp( const_cast<Epetra_CrsMatrix *>(testCase->getMass()), false );
00116   // tell the user that we're done
00117   printer.stream(Errors) << " done." << endl;
00118 
00119 
00120   //************************************
00121   // Select the Preconditioner
00122   //************************************
00123   //
00124   Teuchos::RCP<Ifpack_Preconditioner> prec;
00125   Teuchos::RCP<Epetra_Operator> PrecOp;
00126   if (usePrec) {
00127     printer.stream(Errors) << "Constructing Incomplete Cholesky preconditioner..." << flush;
00128     Ifpack precFactory;
00129     // additive-Schwartz incomplete Cholesky with thresholding; see IFPACK documentation
00130     string precType = "IC stand-alone";
00131     int overlapLevel = 0;
00132     prec = Teuchos::rcp( precFactory.Create(precType,K.get(),overlapLevel) );
00133     // parameters for preconditioner
00134     Teuchos::ParameterList precParams;
00135     precParams.set("fact: drop tolerance",prec_dropTol);
00136     precParams.set("fact: level-of-fill",prec_lofill);
00137     IFPACK_CHK_ERR(prec->SetParameters(precParams));
00138     IFPACK_CHK_ERR(prec->Initialize());
00139     IFPACK_CHK_ERR(prec->Compute());
00140     //
00141     printer.stream(Errors) 
00142       << " done." << endl;
00143     // encapsulate this preconditioner into a IFPACKPrecOp class
00144     PrecOp = Teuchos::rcp( new Epetra_InvOperator(&*prec) );
00145   }
00146 
00147 
00148   //************************************
00149   // Call the LOBPCG solver manager
00150   //***********************************
00151   // Create an Epetra_MultiVector for an initial vector to start the solver.
00152   // Note:  This needs to have the same number of columns as the blocksize.
00153   //
00154   Teuchos::RCP<Epetra_MultiVector> ivec 
00155     = Teuchos::rcp( new Epetra_MultiVector(K->OperatorDomainMap(), blockSize) );
00156   ivec->Random();
00157   // Create the eigenproblem.
00158   //
00159   Teuchos::RCP<BasicEigenproblem<double, MV, OP> > MyProblem =
00160     Teuchos::rcp( new BasicEigenproblem<double, MV, OP>(K, M, ivec) );
00161   // Inform the eigenproblem that the operator K is symmetric
00162   //
00163   MyProblem->setHermitian(true);
00164   // Pass the preconditioner to the eigenproblem
00165   //
00166   if (usePrec) {
00167     MyProblem->setPrec(PrecOp);
00168   }
00169   // Set the number of eigenvalues requested
00170   //
00171   MyProblem->setNEV( nev );
00172   // Inform the eigenproblem that you are finishing passing it information
00173   //
00174   bool boolret = MyProblem->setProblem();
00175   if (boolret != true) {
00176     printer.print(Errors,"Anasazi::BasicEigenproblem::setProblem() returned an error.\n");
00177 #ifdef HAVE_MPI
00178     MPI_Finalize();
00179 #endif
00180     return -1;
00181   }
00182 
00183 
00184   //************************************
00185   // Create parameter list to pass into the solver manager
00186   //************************************
00187   //
00188   Teuchos::ParameterList MyPL;
00189   MyPL.set( "Verbosity", verbosity );
00190   MyPL.set( "Which", which );
00191   MyPL.set( "Block Size", blockSize );
00192   MyPL.set( "Maximum Iterations", maxIters );
00193   MyPL.set( "Convergence Tolerance", tol );
00194   MyPL.set( "Full Ortho", true );
00195   MyPL.set( "Use Locking", true );
00196   //
00197   // Create the solver manager
00198   LOBPCGSolMgr<double, MV, OP> MySolverMan(MyProblem, MyPL);
00199 
00200 
00201   //************************************
00202   // Solve the problem
00203   //************************************
00204   //
00205   printer.stream(Errors) << "Solving eigenvalue problem..." << endl;
00206   ReturnType returnCode = MySolverMan.solve();
00207   // print some precond info
00208   if (usePrec) {
00209     printer.stream(FinalSummary) << *prec << endl;
00210   }
00211 
00212 
00213   //************************************
00214   // Get the eigenvalues and eigenvectors from the eigenproblem
00215   //************************************
00216   //
00217   Eigensolution<double,MV> sol = MyProblem->getSolution();
00218   std::vector<Value<double> > evals = sol.Evals;
00219   Teuchos::RCP<MV> evecs = sol.Evecs;
00220 
00221 
00222   //************************************
00223   // Compute residuals, just for funsies
00224   //************************************
00225   //
00226   std::vector<double> normR(sol.numVecs);
00227   if (sol.numVecs > 0) {
00228     Teuchos::SerialDenseMatrix<int,double> T(sol.numVecs, sol.numVecs);
00229     Epetra_MultiVector Kvec( K->OperatorDomainMap(), evecs->NumVectors() );
00230     Epetra_MultiVector Mvec( M->OperatorDomainMap(), evecs->NumVectors() );
00231     T.putScalar(0.0); 
00232     for (int i=0; i<sol.numVecs; i++) {
00233       T(i,i) = evals[i].realpart;
00234     }
00235     K->Apply( *evecs, Kvec );
00236     M->Apply( *evecs, Mvec );  
00237     MVT::MvTimesMatAddMv( -1.0, Mvec, T, 1.0, Kvec );
00238     MVT::MvNorm( Kvec, normR );
00239   }
00240 
00241 
00242   //************************************
00243   // Print the results
00244   //************************************
00245   //
00246   std::ostringstream os;
00247   os.setf(std::ios_base::right, std::ios_base::adjustfield);
00248   os<<"Solver manager returned " << (returnCode == Converged ? "converged." : "unconverged.") << endl;
00249   os<<endl;
00250   os<<"------------------------------------------------------"<<endl;
00251   os<<std::setw(16)<<"Eigenvalue"
00252     <<std::setw(18)<<"Direct Residual"
00253     <<endl;
00254   os<<"------------------------------------------------------"<<endl;
00255   for (int i=0; i<sol.numVecs; i++) {
00256     os<<std::setw(16)<<evals[i].realpart
00257       <<std::setw(18)<<normR[i]/evals[i].realpart
00258       <<endl;
00259   }
00260   os<<"------------------------------------------------------"<<endl;
00261   printer.print(Errors,os.str());
00262 
00263 #ifdef HAVE_MPI
00264   MPI_Finalize();
00265 #endif
00266   return 0;
00267 }