MVOPTester/MVOPTesterEx.cpp

This is an example of how to use the Anasazi::TestMultiVecTraits() and Anasazi::TestOperatorTraits() methods.

00001 #include "AnasaziConfigDefs.hpp"
00002 #include "AnasaziEpetraAdapter.hpp"
00003 #include "AnasaziMVOPTester.hpp"
00004 #include "AnasaziBasicOutputManager.hpp"
00005 
00006 #include "Epetra_Map.h"
00007 #include "Epetra_CrsMatrix.h"
00008 #include "Epetra_Comm.h"
00009 #include "Teuchos_StandardCatchMacros.hpp"
00010 
00011 #ifdef HAVE_MPI
00012 #include "mpi.h"
00013 #include "Epetra_MpiComm.h"
00014 #else
00015 #include "Epetra_SerialComm.h"
00016 #endif
00017 
00018 int main(int argc, char *argv[])
00019 {
00020 
00021 #ifdef HAVE_MPI
00022   // Initialize MPI and setup an Epetra communicator
00023   MPI_Init(&argc,&argv);
00024   Epetra_MpiComm Comm(MPI_COMM_WORLD);
00025 #else
00026   // If we aren't using MPI, then setup a serial communicator.
00027   Epetra_SerialComm Comm;
00028 #endif
00029 
00030   bool success = true;
00031   try {
00032 
00033     // Create an output manager to handle the I/O from the solver
00034     Teuchos::RCP<Anasazi::OutputManager<double> > MyOM = Teuchos::rcp( new Anasazi::BasicOutputManager<double>() );
00035     MyOM->setVerbosity( Anasazi::Warnings );
00036     MyOM->stream(Anasazi::Errors) << Anasazi::Anasazi_Version() << endl << endl;
00037 
00038     // Number of global elements
00039     int dim = 100;
00040     int blockSize = 5;
00041 
00042     // Construct a Map that puts approximately the same number of 
00043     // equations on each processor.
00044     Epetra_Map Map(dim, 0, Comm);
00045 
00046     // Get update list and number of local equations from newly created Map.
00047     int NumMyElements = Map.NumMyElements();
00048     int * MyGlobalElements = new int[NumMyElements];
00049     Map.MyGlobalElements(MyGlobalElements);
00050 
00051     // Create an integer vector NumNz that is used to build the Petra Matrix.
00052     // NumNz[i] is the Number of OFF-DIAGONAL term for the ith global equation 
00053     // on this processor
00054     int * NumNz = new int[NumMyElements];
00055 
00056     // We are building a tridiagonal matrix where each row has (-1 2 -1)
00057     // So we need 2 off-diagonal terms (except for the first and last equation)
00058     for (int i=0; i<NumMyElements; i++) {
00059       if (MyGlobalElements[i]==0 || MyGlobalElements[i] == dim-1) {
00060         NumNz[i] = 2;
00061       }
00062       else {
00063         NumNz[i] = 3;
00064       }
00065     }
00066 
00067     // Create an Epetra_Matrix
00068     Teuchos::RCP<Epetra_CrsMatrix> A 
00069       = Teuchos::rcp( new Epetra_CrsMatrix(Copy, Map, NumNz) );
00070 
00071     // Add  rows one-at-a-time
00072     // Need some vectors to help
00073     // Off diagonal Values will always be -1
00074     double *Values = new double[2];
00075     Values[0] = -1.0; Values[1] = -1.0;
00076     int *Indices = new int[2];
00077     double two = 2.0;
00078     int NumEntries;
00079     for (int i=0; i<NumMyElements; i++) {
00080       if (MyGlobalElements[i]==0) {
00081         Indices[0] = 1;
00082         NumEntries = 1;
00083       }
00084       else if (MyGlobalElements[i] == dim-1) {
00085         Indices[0] = dim-2;
00086         NumEntries = 1;
00087       }
00088       else {
00089         Indices[0] = MyGlobalElements[i]-1;
00090         Indices[1] = MyGlobalElements[i]+1;
00091         NumEntries = 2;
00092       }
00093       int ierr = A->InsertGlobalValues(MyGlobalElements[i],NumEntries,Values,Indices);
00094       assert(ierr==0);
00095       // Put in the diagonal entry
00096       ierr = A->InsertGlobalValues(MyGlobalElements[i],1,&two,MyGlobalElements+i);
00097       assert(ierr==0);
00098     }
00099 
00100     // Finish building the epetra matrix A
00101     int ierr = A->FillComplete();
00102     assert(ierr==0);
00103 
00104     // Issue several useful typedefs;
00105     // The MultiVecTraits class is for defining ....
00106     typedef Epetra_MultiVector MV;
00107     typedef Epetra_Operator OP;
00108     typedef Anasazi::OperatorTraits<double, MV, OP> OPT;
00109 
00110     // Create an Epetra_MultiVector for an initial vector to start the solver.
00111     // Note that this needs to have the same number of columns as the blocksize.
00112     Teuchos::RCP<MV> ivec 
00113       = Teuchos::rcp( new MV(Map, blockSize) );
00114     ivec->Random();
00115 
00116     // test the multivector class
00117     bool testret;
00118     testret = Anasazi::TestMultiVecTraits<double,MV>(MyOM,ivec);
00119     if (testret) {
00120       MyOM->print(Anasazi::Warnings,"*** PASSED TestMultiVecTraits() ***\n");
00121     }
00122     else {
00123       MyOM->print(Anasazi::Warnings,"*** FAILED TestMultiVecTraits() ***\n");
00124     }
00125 
00126 
00127     // test the operator class
00128     testret = Anasazi::TestOperatorTraits<double,MV,OP>(MyOM,ivec,A);
00129     if (testret) {
00130       MyOM->print(Anasazi::Warnings,"*** PASSED TestOperatorTraits() ***\n");
00131     }
00132     else {
00133       MyOM->print(Anasazi::Warnings,"*** FAILED TestOperatorTraits() ***\n");
00134     }
00135 
00136     // Release all objects
00137     delete [] NumNz;
00138     delete [] Values;
00139     delete [] Indices;
00140     delete [] MyGlobalElements;
00141   }
00142   TEUCHOS_STANDARD_CATCH_STATEMENTS(true,cout,success);
00143 
00144 #ifdef HAVE_MPI
00145   MPI_Finalize();
00146 #endif
00147 
00148   return 0;
00149 }

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