BlockDavidson/BlockDavidsonEpetraEx.cpp

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

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 int main(int argc, char *argv[]) {
00018 
00019 #ifdef HAVE_MPI
00020   // Initialize MPI
00021   //
00022   MPI_Init(&argc,&argv);
00023 #endif
00024 
00025   // Create an Epetra communicator
00026   //
00027 #ifdef HAVE_MPI
00028   Epetra_MpiComm Comm(MPI_COMM_WORLD);
00029 #else
00030   Epetra_SerialComm Comm;
00031 #endif
00032 
00033   // Create an Anasazi output manager
00034   //
00035   Anasazi::BasicOutputManager<double> printer;
00036   printer.stream(Anasazi::Errors) << Anasazi::Anasazi_Version() << std::endl << std::endl;
00037 
00038   // Get the sorting string from the command line
00039   //
00040   std::string which("LM");
00041   Teuchos::CommandLineProcessor cmdp(false,true);
00042   cmdp.setOption("sort",&which,"Targetted eigenvalues (SM or LM).");
00043   if (cmdp.parse(argc,argv) != Teuchos::CommandLineProcessor::PARSE_SUCCESSFUL) {
00044 #ifdef HAVE_MPI
00045     MPI_Finalize();
00046 #endif
00047     return -1;
00048   }
00049 
00050   // Dimension of the matrix
00051   //
00052   // Discretization points in any one direction.
00053   //
00054   const int nx = 10;                    
00055   // Size of matrix nx*nx
00056   //
00057   const int NumGlobalElements = nx*nx;  
00058 
00059   // Construct a Map that puts approximately the same number of
00060   // equations on each processor.
00061   //
00062   Epetra_Map Map(NumGlobalElements, 0, Comm);
00063 
00064   // Get update list and number of local equations from newly created Map.
00065   //
00066   int NumMyElements = Map.NumMyElements();
00067 
00068   std::vector<int> MyGlobalElements(NumMyElements);
00069       Map.MyGlobalElements(&MyGlobalElements[0]);
00070 
00071   // Create an integer vector NumNz that is used to build the Petra Matrix.
00072   // NumNz[i] is the Number of OFF-DIAGONAL term for the ith global equation
00073   // on this processor
00074   //
00075   std::vector<int> NumNz(NumMyElements);
00076 
00077   /* We are building a matrix of block structure:
00078   
00079       | T -I          |
00080       |-I  T -I       |
00081       |   -I  T       |
00082       |        ...  -I|
00083       |           -I T|
00084 
00085    where each block is dimension nx by nx and the matrix is on the order of
00086    nx*nx.  The block T is a tridiagonal matrix. 
00087   */
00088   for (int i=0; i<NumMyElements; i++) {
00089     if (MyGlobalElements[i] == 0 || MyGlobalElements[i] == NumGlobalElements-1 || 
00090         MyGlobalElements[i] == nx-1 || MyGlobalElements[i] == nx*(nx-1) ) {
00091       NumNz[i] = 3;
00092     }
00093     else if (MyGlobalElements[i] < nx || MyGlobalElements[i] > nx*(nx-1) || 
00094              MyGlobalElements[i]%nx == 0 || (MyGlobalElements[i]+1)%nx == 0) {
00095       NumNz[i] = 4;
00096     }
00097     else {
00098       NumNz[i] = 5;
00099     }
00100   }
00101 
00102   // Create an Epetra_Matrix
00103   //
00104   Teuchos::RCP<Epetra_CrsMatrix> A = Teuchos::rcp( new Epetra_CrsMatrix(Copy, Map, &NumNz[0]) );
00105 
00106   // Compute coefficients for discrete convection-diffution operator
00107   //
00108   const double one = 1.0;
00109   std::vector<double> Values(4);
00110   std::vector<int> Indices(4);
00111   double rho = 0.0;  
00112   double h = one /(nx+1);
00113   double h2 = h*h;
00114   double c = 5.0e-01*rho/ h;
00115   Values[0] = -one/h2 - c; Values[1] = -one/h2 + c; Values[2] = -one/h2; Values[3]= -one/h2;
00116   double diag = 4.0 / h2;
00117   int NumEntries;
00118   
00119   for (int i=0; i<NumMyElements; i++)
00120   {
00121     if (MyGlobalElements[i]==0)
00122     {
00123       Indices[0] = 1;
00124       Indices[1] = nx;
00125       NumEntries = 2;
00126       int info = A->InsertGlobalValues(MyGlobalElements[i], NumEntries, &Values[1], &Indices[0]);
00127       TEST_FOR_EXCEPTION( info != 0, std::runtime_error, "Failure in InsertGlobalValues()" );
00128     }
00129     else if (MyGlobalElements[i] == nx*(nx-1))
00130     {
00131       Indices[0] = nx*(nx-1)+1;
00132       Indices[1] = nx*(nx-2);
00133       NumEntries = 2;
00134       int info = A->InsertGlobalValues(MyGlobalElements[i], NumEntries, &Values[1], &Indices[0]);
00135       TEST_FOR_EXCEPTION( info != 0, std::runtime_error, "Failure in InsertGlobalValues()" );
00136     }
00137     else if (MyGlobalElements[i] == nx-1)
00138     {
00139       Indices[0] = nx-2;
00140       NumEntries = 1;
00141       int info = A->InsertGlobalValues(MyGlobalElements[i], NumEntries, &Values[0], &Indices[0]);
00142       assert( info==0 );
00143       Indices[0] = 2*nx-1;
00144       info = A->InsertGlobalValues(MyGlobalElements[i], NumEntries, &Values[2], &Indices[0]);
00145       TEST_FOR_EXCEPTION( info != 0, std::runtime_error, "Failure in InsertGlobalValues()" );
00146     }
00147     else if (MyGlobalElements[i] == NumGlobalElements-1)
00148     {
00149       Indices[0] = NumGlobalElements-2;
00150       NumEntries = 1;
00151       int info = A->InsertGlobalValues(MyGlobalElements[i], NumEntries, &Values[0], &Indices[0]);
00152       assert( info==0 );
00153       Indices[0] = nx*(nx-1)-1;
00154       info = A->InsertGlobalValues(MyGlobalElements[i], NumEntries, &Values[2], &Indices[0]);
00155       TEST_FOR_EXCEPTION( info != 0, std::runtime_error, "Failure in InsertGlobalValues()" );
00156     }
00157     else if (MyGlobalElements[i] < nx)
00158     {
00159       Indices[0] = MyGlobalElements[i]-1;
00160       Indices[1] = MyGlobalElements[i]+1;
00161       Indices[2] = MyGlobalElements[i]+nx;
00162       NumEntries = 3;
00163       int info = A->InsertGlobalValues(MyGlobalElements[i], NumEntries, &Values[0], &Indices[0]);
00164       TEST_FOR_EXCEPTION( info != 0, std::runtime_error, "Failure in InsertGlobalValues()" );
00165     }
00166     else if (MyGlobalElements[i] > nx*(nx-1))
00167     {
00168       Indices[0] = MyGlobalElements[i]-1;
00169       Indices[1] = MyGlobalElements[i]+1;
00170       Indices[2] = MyGlobalElements[i]-nx;
00171       NumEntries = 3;
00172       int info = A->InsertGlobalValues(MyGlobalElements[i], NumEntries, &Values[0], &Indices[0]);
00173       TEST_FOR_EXCEPTION( info != 0, std::runtime_error, "Failure in InsertGlobalValues()" );
00174     }
00175     else if (MyGlobalElements[i]%nx == 0)
00176     {
00177       Indices[0] = MyGlobalElements[i]+1;
00178       Indices[1] = MyGlobalElements[i]-nx;
00179       Indices[2] = MyGlobalElements[i]+nx;
00180       NumEntries = 3;
00181       int info = A->InsertGlobalValues(MyGlobalElements[i], NumEntries, &Values[1], &Indices[0]);
00182       TEST_FOR_EXCEPTION( info != 0, std::runtime_error, "Failure in InsertGlobalValues()" );
00183     }
00184     else if ((MyGlobalElements[i]+1)%nx == 0)
00185     {
00186       Indices[0] = MyGlobalElements[i]-nx;
00187       Indices[1] = MyGlobalElements[i]+nx;
00188       NumEntries = 2;
00189       int info = A->InsertGlobalValues(MyGlobalElements[i], NumEntries, &Values[2], &Indices[0]);
00190       assert( info==0 );
00191       Indices[0] = MyGlobalElements[i]-1;
00192       NumEntries = 1;
00193       info = A->InsertGlobalValues(MyGlobalElements[i], NumEntries, &Values[0], &Indices[0]);
00194       TEST_FOR_EXCEPTION( info != 0, std::runtime_error, "Failure in InsertGlobalValues()" );
00195     }
00196     else
00197     {
00198       Indices[0] = MyGlobalElements[i]-1;
00199       Indices[1] = MyGlobalElements[i]+1;
00200       Indices[2] = MyGlobalElements[i]-nx;
00201       Indices[3] = MyGlobalElements[i]+nx;
00202       NumEntries = 4;
00203       int info = A->InsertGlobalValues(MyGlobalElements[i], NumEntries, &Values[0], &Indices[0]);
00204       TEST_FOR_EXCEPTION( info != 0, std::runtime_error, "Failure in InsertGlobalValues()" );
00205     }
00206     // Put in the diagonal entry
00207     int info = A->InsertGlobalValues(MyGlobalElements[i], 1, &diag, &MyGlobalElements[i]);
00208     TEST_FOR_EXCEPTION( info != 0, std::runtime_error, "Failure in InsertGlobalValues()" );
00209   }
00210 
00211   // Finish up
00212   int info = A->FillComplete();
00213   assert( info==0 );
00214   A->SetTracebackMode(1); // Shutdown Epetra Warning tracebacks
00215 
00216   // Create a identity matrix for the temporary mass matrix
00217   Teuchos::RCP<Epetra_CrsMatrix> M = Teuchos::rcp( new Epetra_CrsMatrix(Copy, Map, 1) );
00218   for (int i=0; i<NumMyElements; i++)
00219   {
00220     Values[0] = one;
00221     Indices[0] = i;
00222     NumEntries = 1;
00223     info = M->InsertGlobalValues(MyGlobalElements[i], NumEntries, &Values[0], &Indices[0]);
00224     assert( info==0 );
00225   }
00226   // Finish up
00227   info = M->FillComplete();
00228   assert( info==0 );
00229   M->SetTracebackMode(1); // Shutdown Epetra Warning tracebacks
00230   
00231   //************************************
00232   // Call the Block Davidson solver manager
00233   //***********************************
00234   //
00235   //  Variables used for the Block Davidson Method
00236   //
00237   const int    nev         = 4;
00238   const int    blockSize   = 5;
00239   const int    numBlocks   = 8;
00240   const int    maxRestarts = 100;
00241   const double tol         = 1.0e-8;
00242 
00243   typedef Epetra_MultiVector MV;
00244   typedef Epetra_Operator OP;
00245   typedef Anasazi::MultiVecTraits<double, Epetra_MultiVector> MVT;
00246 
00247   // Create an Epetra_MultiVector for an initial vector to start the solver.
00248   // Note:  This needs to have the same number of columns as the blocksize.
00249   //
00250   Teuchos::RCP<Epetra_MultiVector> ivec = Teuchos::rcp( new Epetra_MultiVector(Map, blockSize) );
00251   ivec->Random();
00252 
00253   // Create the eigenproblem.
00254   //
00255   Teuchos::RCP<Anasazi::BasicEigenproblem<double, MV, OP> > MyProblem =
00256     Teuchos::rcp( new Anasazi::BasicEigenproblem<double, MV, OP>(A, ivec) );
00257 
00258   // Inform the eigenproblem that the operator A is symmetric
00259   //
00260   MyProblem->setHermitian(true);
00261 
00262   // Set the number of eigenvalues requested
00263   //
00264   MyProblem->setNEV( nev );
00265 
00266   // Inform the eigenproblem that you are finishing passing it information
00267   //
00268   bool boolret = MyProblem->setProblem();
00269   if (boolret != true) {
00270     printer.print(Anasazi::Errors,"Anasazi::BasicEigenproblem::setProblem() returned an error.\n");
00271 #ifdef HAVE_MPI
00272     MPI_Finalize();
00273 #endif
00274     return -1;
00275   }
00276 
00277   // Create parameter list to pass into the solver manager
00278   //
00279   Teuchos::ParameterList MyPL;
00280   MyPL.set( "Which", which );
00281   MyPL.set( "Block Size", blockSize );
00282   MyPL.set( "Num Blocks", numBlocks );
00283   MyPL.set( "Maximum Restarts", maxRestarts );
00284   MyPL.set( "Convergence Tolerance", tol );
00285   //
00286   // Create the solver manager
00287   Anasazi::BlockDavidsonSolMgr<double, MV, OP> MySolverMan(MyProblem, MyPL);
00288 
00289   // Solve the problem
00290   //
00291   Anasazi::ReturnType returnCode = MySolverMan.solve();
00292 
00293   // Get the eigenvalues and eigenvectors from the eigenproblem
00294   //
00295   Anasazi::Eigensolution<double,MV> sol = MyProblem->getSolution();
00296   std::vector<Anasazi::Value<double> > evals = sol.Evals;
00297   Teuchos::RCP<MV> evecs = sol.Evecs;
00298 
00299   // Compute residuals.
00300   //
00301   std::vector<double> normR(sol.numVecs);
00302   if (sol.numVecs > 0) {
00303     Teuchos::SerialDenseMatrix<int,double> T(sol.numVecs, sol.numVecs);
00304     Epetra_MultiVector tempAevec( Map, sol.numVecs );
00305     T.putScalar(0.0); 
00306     for (int i=0; i<sol.numVecs; i++) {
00307       T(i,i) = evals[i].realpart;
00308     }
00309     A->Apply( *evecs, tempAevec );
00310     MVT::MvTimesMatAddMv( -1.0, *evecs, T, 1.0, tempAevec );
00311     MVT::MvNorm( tempAevec, normR );
00312   }
00313 
00314   // Print the results
00315   //
00316   std::ostringstream os;
00317   os.setf(std::ios_base::right, std::ios_base::adjustfield);
00318   os<<"Solver manager returned " << (returnCode == Anasazi::Converged ? "converged." : "unconverged.") << std::endl;
00319   os<<std::endl;
00320   os<<"------------------------------------------------------"<<std::endl;
00321   os<<std::setw(16)<<"Eigenvalue"
00322     <<std::setw(18)<<"Direct Residual"
00323     <<std::endl;
00324   os<<"------------------------------------------------------"<<std::endl;
00325   for (int i=0; i<sol.numVecs; i++) {
00326     os<<std::setw(16)<<evals[i].realpart
00327       <<std::setw(18)<<normR[i]/evals[i].realpart
00328       <<std::endl;
00329   }
00330   os<<"------------------------------------------------------"<<std::endl;
00331   printer.print(Anasazi::Errors,os.str());
00332 
00333 #ifdef HAVE_MPI
00334   MPI_Finalize();
00335 #endif
00336   return 0;
00337 }

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