LOBPCG/LOBPCGEpetraExSimple.cpp

This is an example of how to use the Anasazi::SimpleLOBPCGSolMgr solver manager.

00001 #include "AnasaziConfigDefs.hpp"
00002 #include "AnasaziBasicEigenproblem.hpp"
00003 #include "AnasaziSimpleLOBPCGSolMgr.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     int 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   //
00286   // Create the solver manager
00287   SimpleLOBPCGSolMgr<double, MV, OP> MySolverMan(MyProblem, MyPL);
00288 
00289   // Solve the problem
00290   //
00291   ReturnType returnCode = MySolverMan.solve();
00292 
00293   // Get the eigenvalues and eigenvectors from the eigenproblem
00294   //
00295   Eigensolution<double,MV> sol = MyProblem->getSolution();
00296   std::vector<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 == Converged ? "converged." : "unconverged.") << endl;
00319   os<<endl;
00320   os<<"------------------------------------------------------"<<endl;
00321   os<<std::setw(16)<<"Eigenvalue"
00322     <<std::setw(18)<<"Direct Residual"
00323     <<endl;
00324   os<<"------------------------------------------------------"<<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       <<endl;
00329   }
00330   os<<"------------------------------------------------------"<<endl;
00331   printer.print(Errors,os.str());
00332 
00333 #ifdef HAVE_MPI
00334   MPI_Finalize();
00335 #endif
00336   return 0;
00337 }

Generated on Tue Oct 20 12:45:18 2009 for Anasazi by doxygen 1.4.7