cxx_main.cpp

Go to the documentation of this file.
00001 //@HEADER
00002 // ************************************************************************
00003 // 
00004 //
00005 //                 Belos: Block Linear Solvers 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 test uses the MVOPTester.hpp functions to test the Belos adapters
00031 //  to Epetra and Thyra.
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 "BelosConfigDefs.hpp"
00047 #include "BelosMVOPTester.hpp"
00048 #include "BelosEpetraAdapter.hpp"
00049 
00050 #ifdef HAVE_EPETRA_THYRA
00051 #include "BelosThyraAdapter.hpp"
00052 #include "Thyra_EpetraThyraWrappers.hpp"
00053 #include "Thyra_EpetraLinearOp.hpp"
00054 #endif
00055 
00056 
00057 
00058 int main(int argc, char *argv[])
00059 {
00060   
00061   using Teuchos::rcp_implicit_cast;
00062 
00063   int i, ierr, gerr;
00064   gerr = 0;
00065 
00066 #ifdef HAVE_MPI
00067   // Initialize MPI and setup an Epetra communicator
00068   MPI_Init(&argc,&argv);
00069   Teuchos::RCP<Epetra_MpiComm> Comm = Teuchos::rcp( new Epetra_MpiComm(MPI_COMM_WORLD) );
00070 #else
00071   // If we aren't using MPI, then setup a serial communicator.
00072   Teuchos::RCP<Epetra_SerialComm> Comm = Teuchos::rcp( new Epetra_SerialComm() );
00073 #endif
00074 
00075 
00076    // number of global elements
00077   int dim = 100;
00078   int blockSize = 3;
00079 
00080   // PID info
00081   int MyPID = Comm->MyPID();
00082   bool verbose = 0;
00083 
00084   if (argc>1) {
00085     if (argv[1][0]=='-' && argv[1][1]=='v') {
00086       verbose = true;
00087     }
00088   }
00089 
00090   // Construct a Map that puts approximately the same number of 
00091   // equations on each processor.
00092   Teuchos::RCP<Epetra_Map> Map = Teuchos::rcp( new Epetra_Map(dim, 0, *Comm) );
00093   
00094   // Get update list and number of local equations from newly created Map.
00095   int NumMyElements = Map->NumMyElements();
00096   std::vector<int> MyGlobalElements(NumMyElements);
00097   Map->MyGlobalElements(&MyGlobalElements[0]);
00098 
00099   // Create an integer std::vector NumNz that is used to build the Petra Matrix.
00100   // NumNz[i] is the Number of OFF-DIAGONAL term for the ith global equation 
00101   // on this processor
00102   std::vector<int> NumNz(NumMyElements);
00103 
00104   // We are building a tridiagonal matrix where each row has (-1 2 -1)
00105   // So we need 2 off-diagonal terms (except for the first and last equation)
00106   for (i=0; i<NumMyElements; i++) {
00107     if (MyGlobalElements[i]==0 || MyGlobalElements[i] == dim-1) {
00108       NumNz[i] = 2;
00109     }
00110     else {
00111       NumNz[i] = 3;
00112     }
00113   }
00114 
00115   // Create an Epetra_Matrix
00116   Teuchos::RCP<Epetra_CrsMatrix> A = Teuchos::rcp( new Epetra_CrsMatrix(Copy, *Map, &NumNz[0]) );
00117    
00118   // Add  rows one-at-a-time
00119   // Need some vectors to help
00120   // Off diagonal Values will always be -1
00121   std::vector<double> Values(2);
00122   Values[0] = -1.0; Values[1] = -1.0;
00123   std::vector<int> Indices(2);
00124   double two = 2.0;
00125   int NumEntries;
00126   for (i=0; i<NumMyElements; i++) {
00127     if (MyGlobalElements[i]==0) {
00128       Indices[0] = 1;
00129       NumEntries = 1;
00130     }
00131     else if (MyGlobalElements[i] == dim-1) {
00132       Indices[0] = dim-2;
00133       NumEntries = 1;
00134     }
00135     else {
00136       Indices[0] = MyGlobalElements[i]-1;
00137       Indices[1] = MyGlobalElements[i]+1;
00138       NumEntries = 2;
00139     }
00140     ierr = A->InsertGlobalValues(MyGlobalElements[i],NumEntries,&Values[0],&Indices[0]);
00141     assert(ierr==0);
00142     // Put in the diagonal entry
00143     ierr = A->InsertGlobalValues(MyGlobalElements[i],1,&two,&MyGlobalElements[i]);
00144     assert(ierr==0);
00145   }
00146    
00147   // Finish building the epetra matrix A
00148   ierr = A->FillComplete();
00149   assert(ierr==0);
00150 
00151   // Create an Belos::EpetraOp from this Epetra_CrsMatrix
00152   Teuchos::RCP<Belos::EpetraOp> op = Teuchos::rcp(new Belos::EpetraOp(A));
00153 
00154   // Issue several useful typedefs;
00155   typedef Belos::MultiVec<double> EMV;
00156   typedef Belos::Operator<double> EOP;
00157 
00158   // Create an Epetra_MultiVector for an initial std::vector to start the solver.
00159   // Note that this needs to have the same number of columns as the blocksize.
00160   Teuchos::RCP<Belos::EpetraMultiVec> ivec = Teuchos::rcp( new Belos::EpetraMultiVec(*Map, blockSize) );
00161   ivec->Random();
00162 
00163   // Create an output manager to handle the I/O from the solver
00164   Teuchos::RCP<Belos::OutputManager<double> > MyOM = Teuchos::rcp( new Belos::OutputManager<double>( MyPID ) );
00165   if (verbose) {
00166     MyOM->setVerbosity( Belos::Errors + Belos::Warnings );
00167   }
00168 
00169 #ifdef HAVE_EPETRA_THYRA
00170   typedef Thyra::MultiVectorBase<double> TMVB;
00171   typedef Thyra::LinearOpBase<double>    TLOB;
00172   // create thyra objects from the epetra objects
00173 
00174   // first, a Thyra::VectorSpaceBase
00175   Teuchos::RCP<const Thyra::VectorSpaceBase<double> > epetra_vs = 
00176     Thyra::create_VectorSpace(Map);
00177 
00178   // then, a MultiVectorBase (from the Epetra_MultiVector)
00179   Teuchos::RCP<Thyra::MultiVectorBase<double> > thyra_ivec = 
00180     Thyra::create_MultiVector(rcp_implicit_cast<Epetra_MultiVector>(ivec),epetra_vs);
00181 
00182   // then, a LinearOpBase (from the Epetra_CrsMatrix)
00183   Teuchos::RCP<Thyra::LinearOpBase<double> > thyra_op = 
00184     Teuchos::rcp( new Thyra::EpetraLinearOp(A) );
00185 
00186 
00187   // test the Thyra adapter multivector
00188   ierr = Belos::TestMultiVecTraits<double,TMVB>(MyOM,thyra_ivec);
00189   gerr |= ierr;
00190   switch (ierr) {
00191   case Belos::Ok:
00192     if ( verbose && MyPID==0 ) {
00193       std::cout << "*** ThyraAdapter PASSED TestMultiVecTraits()" << std::endl;
00194     }
00195     break;
00196   case Belos::Error:
00197     if ( verbose && MyPID==0 ) {
00198       std::cout << "*** ThyraAdapter FAILED TestMultiVecTraits() ***" 
00199            << std::endl << std::endl;
00200     }
00201     break;
00202   }
00203 
00204   // test the Thyra adapter operator 
00205   ierr = Belos::TestOperatorTraits<double,TMVB,TLOB>(MyOM,thyra_ivec,thyra_op);
00206   gerr |= ierr;
00207   switch (ierr) {
00208   case Belos::Ok:
00209     if ( verbose && MyPID==0 ) {
00210       std::cout << "*** ThyraAdapter PASSED TestOperatorTraits()" << std::endl;
00211     }
00212     break;
00213   case Belos::Error:
00214     if ( verbose && MyPID==0 ) {
00215       std::cout << "*** ThyraAdapter FAILED TestOperatorTraits() ***" 
00216            << std::endl << std::endl;
00217     }
00218     break;
00219   }
00220 #endif
00221 
00222 #ifdef HAVE_MPI
00223   MPI_Finalize();
00224 #endif
00225 
00226   if (gerr) {
00227     if (verbose && MyPID==0)
00228       std::cout << "End Result: TEST FAILED" << std::endl;  
00229     return -1;
00230   }
00231   //
00232   // Default return value
00233   //
00234   if (verbose && MyPID==0)
00235     std::cout << "End Result: TEST PASSED" << std::endl;
00236   return 0;
00237 
00238 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends
Generated on Wed Apr 13 10:19:51 2011 for Stratimikos Package Browser (Single Doxygen Collection) by  doxygen 1.6.3