LOBPCG/LOBPCGEpetraEx.cpp

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

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

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