Anasazi: BlockDavidson/BlockDavidsonEpetraExGenPrecIfpack.cpp

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