MVOPTester/MVOPTesterEx.cpp

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

00001 //@HEADER
00002 // ************************************************************************
00003 // 
00004 //
00005 //                 Anasazi: Block Eigensolvers Package
00006 //                 Copyright (2004) Sandia Corporation
00007 // 
00008 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
00009 // license for use of this work by or on behalf of the U.S. Government.
00010 // 
00011 // This library is free software; you can redistribute it and/or modify
00012 // it under the terms of the GNU Lesser General Public License as
00013 // published by the Free Software Foundation; either version 2.1 of the
00014 // License, or (at your option) any later version.
00015 //  
00016 // This library is distributed in the hope that it will be useful, but
00017 // WITHOUT ANY WARRANTY; without even the implied warranty of
00018 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00019 // Lesser General Public License for more details.
00020 //  
00021 // You should have received a copy of the GNU Lesser General Public
00022 // License along with this library; if not, write to the Free Software
00023 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
00024 // USA
00025 // Questions? Contact Michael A. Heroux (maherou@sandia.gov) 
00026 // 
00027 // ************************************************************************
00028 //@HEADER
00029 //
00030 //  This example uses the MVOPTester.hpp functions to test the Anasazi adapters
00031 //  to Epetra.
00032 //
00033 
00034 #include "Epetra_Map.h"
00035 #include "Epetra_CrsMatrix.h"
00036 #ifdef HAVE_MPI
00037 #include "mpi.h"
00038 #include "Epetra_MpiComm.h"
00039 #endif
00040 #ifndef __cplusplus
00041 #define __cplusplus
00042 #endif
00043 #include "Epetra_Comm.h"
00044 #include "Epetra_SerialComm.h"
00045 
00046 #include "AnasaziConfigDefs.hpp"
00047 #include "AnasaziMVOPTester.hpp"
00048 #include "AnasaziEpetraAdapter.hpp"
00049 #include "AnasaziBasicOutputManager.hpp"
00050 
00051 int main(int argc, char *argv[])
00052 {
00053   int i;
00054   bool ierr, gerr;
00055   gerr = true;
00056 
00057 #ifdef HAVE_MPI
00058   // Initialize MPI and setup an Epetra communicator
00059   MPI_Init(&argc,&argv);
00060   Teuchos::RCP<Epetra_MpiComm> Comm = Teuchos::rcp( new Epetra_MpiComm(MPI_COMM_WORLD) );
00061 #else
00062   // If we aren't using MPI, then setup a serial communicator.
00063   Teuchos::RCP<Epetra_SerialComm> Comm = Teuchos::rcp( new Epetra_SerialComm() );
00064 #endif
00065 
00066    // number of global elements
00067   int dim = 100;
00068   int blockSize = 5;
00069 
00070   bool verbose = false;
00071   if (argc>1) {
00072     if (argv[1][0]=='-' && argv[1][1]=='v') {
00073       verbose = true;
00074     }
00075   }
00076 
00077   // Construct a Map that puts approximately the same number of 
00078   // equations on each processor.
00079   Teuchos::RCP<Epetra_Map> Map = Teuchos::rcp( new Epetra_Map(dim, 0, *Comm) );
00080   
00081   // Get update list and number of local equations from newly created Map.
00082   int NumMyElements = Map->NumMyElements();
00083   std::vector<int> MyGlobalElements(NumMyElements);
00084   Map->MyGlobalElements(&MyGlobalElements[0]);
00085 
00086   // Create an integer vector NumNz that is used to build the Petra Matrix.
00087   // NumNz[i] is the Number of OFF-DIAGONAL term for the ith global equation 
00088   // on this processor
00089   std::vector<int> NumNz(NumMyElements);
00090 
00091   // We are building a tridiagonal matrix where each row has (-1 2 -1)
00092   // So we need 2 off-diagonal terms (except for the first and last equation)
00093   for (i=0; i<NumMyElements; i++) {
00094     if (MyGlobalElements[i]==0 || MyGlobalElements[i] == dim-1) {
00095       NumNz[i] = 2;
00096     }
00097     else {
00098       NumNz[i] = 3;
00099     }
00100   }
00101 
00102   // Create an Epetra_Matrix
00103   Teuchos::RCP<Epetra_CrsMatrix> A = Teuchos::rcp( new Epetra_CrsMatrix(Copy, *Map, &NumNz[0]) );
00104    
00105   // Add  rows one-at-a-time
00106   // Need some vectors to help
00107   // Off diagonal Values will always be -1
00108   std::vector<double> Values(2);
00109   Values[0] = -1.0; Values[1] = -1.0;
00110   std::vector<int> Indices(2);
00111   double two = 2.0;
00112   int NumEntries;
00113   for (i=0; i<NumMyElements; i++) {
00114     if (MyGlobalElements[i]==0) {
00115       Indices[0] = 1;
00116       NumEntries = 1;
00117     }
00118     else if (MyGlobalElements[i] == dim-1) {
00119       Indices[0] = dim-2;
00120       NumEntries = 1;
00121     }
00122     else {
00123       Indices[0] = MyGlobalElements[i]-1;
00124       Indices[1] = MyGlobalElements[i]+1;
00125       NumEntries = 2;
00126     }
00127     ierr = A->InsertGlobalValues(MyGlobalElements[i],NumEntries,&Values[0],&Indices[0]);
00128     assert(ierr==0);
00129     // Put in the diagonal entry
00130     ierr = A->InsertGlobalValues(MyGlobalElements[i],1,&two,&MyGlobalElements[i]);
00131     assert(ierr==0);
00132   }
00133    
00134   // Finish building the epetra matrix A
00135   ierr = A->FillComplete();
00136   assert(ierr==0);
00137 
00138   // Create an Anasazi::EpetraSymOp from this Epetra_CrsMatrix
00139   Teuchos::RCP<Anasazi::EpetraSymOp> op = Teuchos::rcp(new Anasazi::EpetraSymOp(A));
00140 
00141   // Issue several useful typedefs;
00142   typedef Anasazi::MultiVec<double> EMV;
00143   typedef Anasazi::Operator<double> EOP;
00144 
00145   // Create an Epetra_MultiVector for an initial vector to start the solver.
00146   // Note that this needs to have the same number of columns as the blocksize.
00147   Teuchos::RCP<Anasazi::EpetraMultiVec> ivec = Teuchos::rcp( new Anasazi::EpetraMultiVec(*Map, blockSize) );
00148   ivec->Random();
00149 
00150   // Create an output manager to handle the I/O from the solver
00151   Teuchos::RCP<Anasazi::OutputManager<double> > MyOM = Teuchos::rcp( new Anasazi::BasicOutputManager<double>() );
00152   if (verbose) {
00153     MyOM->setVerbosity( Anasazi::Warnings );
00154   }
00155 
00156   // test the Epetra adapter multivector
00157   ierr = Anasazi::TestMultiVecTraits<double,EMV>(MyOM,ivec);
00158   gerr &= ierr;
00159   if (ierr) {
00160     MyOM->print(Anasazi::Warnings,"*** EpetraAdapter PASSED TestMultiVecTraits()\n");
00161   }
00162   else {
00163     MyOM->print(Anasazi::Warnings,"*** EpetraAdapter FAILED TestMultiVecTraits() ***\n\n");
00164   }
00165 
00166   // test the Epetra adapter operator 
00167   ierr = Anasazi::TestOperatorTraits<double,EMV,EOP>(MyOM,ivec,op);
00168   gerr &= ierr;
00169   if (ierr) {
00170     MyOM->print(Anasazi::Warnings,"*** EpetraAdapter PASSED TestOperatorTraits()\n");
00171   }
00172   else {
00173     MyOM->print(Anasazi::Warnings,"*** EpetraAdapter FAILED TestOperatorTraits() ***\n\n");
00174   }
00175 
00176 
00177 #ifdef HAVE_MPI
00178   MPI_Finalize();
00179 #endif
00180 
00181   if (gerr == false) {
00182     MyOM->print(Anasazi::Warnings,"End Result: TEST FAILED\n");
00183     return -1;
00184   }
00185   //
00186   // Default return value
00187   //
00188   MyOM->print(Anasazi::Warnings,"End Result: TEST PASSED\n");
00189   return 0;
00190 
00191 }

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