Thyra Package Browser (Single Doxygen Collection) Version of the Day
test_epetra_adapters.cpp
Go to the documentation of this file.
00001 // @HEADER
00002 // ***********************************************************************
00003 // 
00004 //    Thyra: Interfaces and Support for Abstract Numerical Algorithms
00005 //                 Copyright (2004) Sandia Corporation
00006 // 
00007 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
00008 // license for use of this work by or on behalf of the U.S. Government.
00009 // 
00010 // Redistribution and use in source and binary forms, with or without
00011 // modification, are permitted provided that the following conditions are
00012 // met:
00013 //
00014 // 1. Redistributions of source code must retain the above copyright
00015 // notice, this list of conditions and the following disclaimer.
00016 //
00017 // 2. Redistributions in binary form must reproduce the above copyright
00018 // notice, this list of conditions and the following disclaimer in the
00019 // documentation and/or other materials provided with the distribution.
00020 //
00021 // 3. Neither the name of the Corporation nor the names of the
00022 // contributors may be used to endorse or promote products derived from
00023 // this software without specific prior written permission.
00024 //
00025 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
00026 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
00027 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
00028 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
00029 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
00030 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
00031 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
00032 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
00033 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
00034 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
00035 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
00036 //
00037 // Questions? Contact Roscoe A. Bartlett (bartlettra@ornl.gov) 
00038 // 
00039 // ***********************************************************************
00040 // @HEADER
00041 
00042 #include "Thyra_EpetraThyraWrappers.hpp"
00043 #include "Thyra_EpetraLinearOp.hpp"
00044 #include "Thyra_TestingTools.hpp"
00045 #include "Thyra_LinearOpTester.hpp"
00046 #include "Thyra_LinearOpWithSolveTester.hpp"
00047 #include "Thyra_DiagonalEpetraLinearOpWithSolveFactory.hpp"
00048 #include "Thyra_LinearOpWithSolveFactoryHelpers.hpp"
00049 #include "Thyra_DefaultSpmdVectorSpace.hpp"
00050 #include "Thyra_MultiVectorStdOps.hpp"
00051 #include "Epetra_SerialComm.h"
00052 #include "Epetra_LocalMap.h"
00053 #include "Epetra_CrsMatrix.h"
00054 #include "Epetra_MultiVector.h"
00055 #include "Epetra_Vector.h"
00056 #include "Teuchos_dyn_cast.hpp"
00057 #include "Teuchos_GlobalMPISession.hpp"
00058 #include "Teuchos_CommandLineProcessor.hpp"
00059 #include "Teuchos_OpaqueWrapper.hpp"
00060 #include "Teuchos_TimeMonitor.hpp"
00061 #include "Teuchos_StandardCatchMacros.hpp"
00062 #ifdef RTOp_USE_MPI
00063 #  include "Epetra_MpiComm.h"
00064 #endif
00065 
00066 //
00067 // Some helper functions
00068 //
00069 
00070 namespace {
00071 
00072 void print_performance_stats(
00073   const int        num_time_samples
00074   ,const double    raw_epetra_time
00075   ,const double    thyra_wrapped_time
00076   ,bool            verbose
00077   ,std::ostream    &out
00078   )
00079 {
00080   if(verbose)
00081     out
00082       << "\nAverage times (out of " << num_time_samples << " samples):\n"
00083       << "  Raw Epetra              = " << (raw_epetra_time/num_time_samples) << std::endl
00084       << "  Thyra Wrapped Epetra    = " << (thyra_wrapped_time/num_time_samples) << std::endl
00085       << "\nRelative performance of Thyra wrapped verses raw Epetra:\n"
00086       << "  ( raw epetra time / thyra wrapped time ) = ( " << raw_epetra_time << " / " << thyra_wrapped_time << " ) = "
00087       << (raw_epetra_time/thyra_wrapped_time) << std::endl;
00088 }
00089 
00090 inline
00091 double sum( const Epetra_MultiVector &ev )
00092 {
00093   std::vector<double> meanValues(ev.NumVectors());
00094   ev.MeanValue(&meanValues[0]);
00095   double sum = 0;
00096   for( int i = 0; i < ev.NumVectors(); ++i ) sum += meanValues[i];
00097   return (ev.Map().NumGlobalElements()*sum);
00098 }
00099 
00100 } // namespace
00101 
00102 /* Testing program for Thyra/Epetra adpaters.
00103  *
00104  * This testing program shows how you can easily mix and match
00105  * different implementations of vectors and multi-vectors for serial
00106  * and SPMD MPI implementations.  This code is worth study to show how
00107  * this is done.
00108  *
00109  * Note that the tests performed do not prove that the Epetra adapters
00110  * (or Epetra itself) perform correctly as only a few post conditions
00111  * are checked.  Because of the simple nature of these computations it
00112  * would be possible to put in more exactly component-wise tests if
00113  * that is needed in the future.
00114  */
00115 int main( int argc, char* argv[] )
00116 {
00117 
00118   using std::endl;
00119 
00120   typedef double Scalar;
00121   typedef Teuchos::ScalarTraits<Scalar> ST;
00122   typedef ST::magnitudeType ScalarMag;
00123   typedef Teuchos::ScalarTraits<ScalarMag> SMT;
00124 
00125   using Teuchos::dyn_cast;
00126   using Teuchos::CommandLineProcessor;
00127   using Teuchos::Ptr;
00128   using Teuchos::RCP;
00129   using Teuchos::Array;
00130   using Teuchos::arcpFromArray;
00131   using Teuchos::rcp;
00132   using Teuchos::rcp_static_cast;
00133   using Teuchos::rcp_const_cast;
00134   using Teuchos::testRelErr;
00135 
00136   using Thyra::passfail;
00137   using Thyra::NOTRANS;
00138   using Thyra::TRANS;
00139   using Thyra::apply;
00140   using Thyra::create_VectorSpace;
00141   using Thyra::create_Vector;
00142   using Thyra::create_MultiVector;
00143   
00144   bool verbose = true;
00145   bool dumpAll = false;
00146   bool success = true;
00147   bool result;
00148 
00149   int procRank = 0;
00150 
00151   Teuchos::GlobalMPISession mpiSession(&argc,&argv);
00152 
00153   RCP<Teuchos::FancyOStream>
00154     out = Teuchos::VerboseObjectBase::getDefaultOStream();
00155   
00156   try {
00157 
00158     Teuchos::Time timer("");
00159 
00160     //
00161     // Read options from the commandline
00162     //
00163 
00164     int     local_dim            = 1000;
00165     int     num_mv_cols          = 4;
00166     double  max_rel_err          = 1e-13;
00167     double  max_rel_warn         = 1e-15;
00168     double  scalar               = 1.5;
00169     double  max_flop_rate        = 1.0; // Turn off by default!
00170 #ifdef RTOp_USE_MPI
00171     bool    useMPI               = true;
00172 #endif
00173     CommandLineProcessor  clp;
00174     clp.throwExceptions(false);
00175     clp.addOutputSetupOptions(true);
00176     clp.setOption( "verbose", "quiet", &verbose, "Determines if any output is printed or not." );
00177     clp.setOption( "dump-all", "no-dump", &dumpAll, "Determines if quantities are dumped or not." );
00178     clp.setOption( "local-dim", &local_dim, "Number of vector elements per process." );
00179     clp.setOption( "num-mv-cols", &num_mv_cols, "Number columns in each multi-vector (>=4)." );
00180     clp.setOption( "max-rel-err-tol", &max_rel_err, "Maximum relative error tolerance for tests." );
00181     clp.setOption( "max-rel-warn-tol", &max_rel_warn, "Maximum relative warning tolerance for tests." );
00182     clp.setOption( "scalar", &scalar, "A scalar used in all computations." );
00183     clp.setOption( "max-flop-rate", &max_flop_rate, "Approx flop rate used for loop timing." );
00184 #ifdef RTOp_USE_MPI
00185     clp.setOption( "use-mpi", "no-use-mpi", &useMPI, "Actually use MPI or just run independent serial programs." );
00186 #endif
00187     CommandLineProcessor::EParseCommandLineReturn parse_return = clp.parse(argc,argv);
00188     if( parse_return != CommandLineProcessor::PARSE_SUCCESSFUL ) return parse_return;
00189 
00190     TEUCHOS_TEST_FOR_EXCEPTION(
00191       num_mv_cols < 4, std::logic_error
00192       ,"Error, num-mv-cols must be >= 4!"
00193       );
00194 
00195     //
00196     // Get basic MPI info
00197     //
00198 
00199 #ifdef RTOp_USE_MPI
00200     MPI_Comm mpiComm;
00201     int numProc;
00202     if(useMPI) {
00203       mpiComm = MPI_COMM_WORLD;
00204       MPI_Comm_size( mpiComm, &numProc );
00205       MPI_Comm_rank( mpiComm, &procRank );
00206     }
00207     else {
00208       mpiComm = MPI_COMM_NULL;
00209       numProc = 1;
00210       procRank = 0;
00211     }
00212 #endif
00213 
00214     if(verbose)
00215       *out
00216         << "\n***"
00217         << "\n*** (A) Creating two vector spaces (an Epetra-based and a non-Epetra-based)"
00218         << "\n***\n";
00219 
00220     //
00221     // Create two different vector spaces (one Epetra and one
00222     // non-Epetra) that should be compatible.
00223     //
00224     RCP<const Epetra_Comm> epetra_comm;
00225     RCP<const Epetra_Map> epetra_map;
00226     RCP<const Thyra::VectorSpaceBase<Scalar> > epetra_vs;
00227     RCP<const Thyra::VectorSpaceBase<Scalar> > non_epetra_vs;
00228 #ifdef RTOp_USE_MPI
00229     if(useMPI) {
00230       //
00231       // Create parallel vector spaces with compatible maps
00232       //
00233       // Epetra vector space
00234       if(verbose) *out << "\nCreating vector space using Epetra_MpiComm ...\n";
00235       epetra_comm = rcp(new Epetra_MpiComm(mpiComm));
00236       epetra_map = rcp(new Epetra_Map(-1,local_dim,0,*epetra_comm));
00237       epetra_vs = Thyra::create_VectorSpace(epetra_map);
00238       // Non-Epetra vector space
00239       if(verbose) *out << "\nCreating Thyra::DefaultSpmdVectorSpace ...\n";
00240       non_epetra_vs = rcp(
00241         new Thyra::DefaultSpmdVectorSpace<Scalar>(
00242           Thyra::create_Comm(epetra_comm)
00243           ,local_dim,-1
00244           )
00245         );
00246     }
00247     else {
00248 #endif
00249       //
00250       // Create serial vector spaces (i.e. VectorSpaceBase::isInCore()==true)
00251       //
00252       // Epetra vector space
00253       if(verbose) *out << "\nCreating vector space using Epetra_SerialComm ...\n";
00254       epetra_comm = rcp(new Epetra_SerialComm);
00255       epetra_map = rcp(new Epetra_LocalMap(local_dim,0,*epetra_comm));
00256       epetra_vs = Thyra::create_VectorSpace(epetra_map);
00257       // Non-Epetra vector space
00258       if(verbose) *out << "\nCreating Thyra::DefaultSpmdVectorSpace ...\n";
00259       non_epetra_vs = Thyra::defaultSpmdVectorSpace<Scalar>(local_dim);
00260 #ifdef RTOp_USE_MPI
00261     }
00262 #endif // end create vector spacdes [Doxygen looks for this!]
00263 
00264 #ifdef RTOp_USE_MPI
00265     const int global_dim = local_dim * numProc;
00266 #else
00267     const int global_dim = local_dim;
00268 #endif
00269 
00270     if(verbose)
00271       *out
00272         << "\nscalar              = " << scalar
00273         << "\nlocal_dim           = " << local_dim
00274         << "\nglobal_dim          = " << global_dim
00275         << "\nnum_mv_cols         = " << num_mv_cols
00276         << "\nepetra_vs.dim()     = " << epetra_vs->dim()
00277         << "\nnon_epetra_vs.dim() = " << non_epetra_vs->dim()
00278         << std::endl;
00279 
00280     //
00281     // Create vectors and multi-vectors from each vector space
00282     //
00283 
00284     RCP<Thyra::VectorBase<Scalar> >
00285       ev1 = createMember(epetra_vs),
00286       ev2 = createMember(epetra_vs);
00287     RCP<Thyra::VectorBase<Scalar> >
00288       nev1 = createMember(non_epetra_vs),
00289       nev2 = createMember(non_epetra_vs);
00290 
00291     RCP<Thyra::MultiVectorBase<Scalar> >
00292       eV1 = createMembers(epetra_vs,num_mv_cols),
00293       eV2 = createMembers(epetra_vs,num_mv_cols);
00294     RCP<Thyra::MultiVectorBase<Scalar> >
00295       neV1 = createMembers(non_epetra_vs,num_mv_cols),
00296       neV2 = createMembers(non_epetra_vs,num_mv_cols);
00297 
00298     if(verbose)
00299       *out
00300         << "\n***"
00301         << "\n*** (B) Testing Epetra and non-Epetra Thyra wrapped objects"
00302         << "\n***\n";
00303 
00304     //
00305     // Check for compatibility of the vector and Multi-vectors
00306     // w.r.t. RTOps
00307     //
00308 
00309     if(verbose) *out << "\n*** (B.1) Testing individual vector/multi-vector RTOps\n";
00310 
00311     Thyra::assign( ev1.ptr(), 0.0 );
00312     Thyra::assign( ev2.ptr(), scalar );
00313     Thyra::assign( nev1.ptr(), 0.0 );
00314     Thyra::assign( nev2.ptr(), scalar );
00315     Thyra::assign( eV1.ptr(), 0.0 );
00316     Thyra::assign( eV2.ptr(), scalar );
00317     Thyra::assign( neV1.ptr(), 0.0 );
00318     Thyra::assign( neV2.ptr(), scalar );
00319 
00320     Scalar
00321       ev1_nrm = Thyra::norm_1(*ev1),
00322       ev2_nrm = Thyra::norm_1(*ev2),
00323       eV1_nrm = Thyra::norm_1(*eV1),
00324       eV2_nrm = Thyra::norm_1(*eV2),
00325       nev1_nrm = Thyra::norm_1(*nev1),
00326       nev2_nrm = Thyra::norm_1(*nev2),
00327       neV1_nrm = Thyra::norm_1(*neV1),
00328       neV2_nrm = Thyra::norm_1(*neV2);
00329 
00330     const std::string s1_n = "fabs(scalar)*global_dim";
00331     const Scalar s1 = fabs(scalar)*global_dim;
00332     
00333     if(!testRelErr("Thyra::norm_1(ev1)",ev1_nrm,"0",Scalar(0),"max_rel_err",max_rel_err,"max_rel_warn",max_rel_warn,out.ptr())) success=false;
00334     if(verbose && dumpAll) *out << "\nev1 =\n" << *ev1;
00335     if(!testRelErr("Thyra::norm_1(ev2)",ev2_nrm,s1_n,s1,"max_rel_err",max_rel_err,"max_rel_warn",max_rel_warn,out.ptr())) success=false;
00336     if(verbose && dumpAll) *out << "\nev2 =\n" << *ev2;
00337     if(!testRelErr("Thyra::norm_1(nev1)",nev1_nrm,"0",Scalar(0),"max_rel_err",max_rel_err,"max_rel_warn",max_rel_warn,out.ptr())) success=false;
00338     if(verbose && dumpAll) *out << "\nnev2 =\n" << *ev1;
00339     if(!testRelErr("Thyra::norm_1(nev2)",nev2_nrm,s1_n,s1,"max_rel_err",max_rel_err,"max_rel_warn",max_rel_warn,out.ptr())) success=false;
00340     if(verbose && dumpAll) *out << "\nnev2 =\n" << *nev2;
00341     if(!testRelErr("Thyra::norm_1(eV1)",eV1_nrm,"0",Scalar(0),"max_rel_err",max_rel_err,"max_rel_warn",max_rel_warn,out.ptr())) success=false;
00342     if(verbose && dumpAll) *out << "\neV1 =\n" << *eV1;
00343     if(!testRelErr("Thyra::norm_1(eV2)",eV2_nrm,s1_n,s1,"max_rel_err",max_rel_err,"max_rel_warn",max_rel_warn,out.ptr())) success=false;
00344     if(verbose && dumpAll) *out << "\neV2 =\n" << *eV2;
00345     if(!testRelErr("Thyra::norm_1(neV1)",neV1_nrm,"0",Scalar(0),"max_rel_err",max_rel_err,"max_rel_warn",max_rel_warn,out.ptr())) success=false;
00346     if(verbose && dumpAll) *out << "\nneV1 =\n" << *neV1;
00347     if(!testRelErr("Thyra::norm_1(neV2)",neV2_nrm,s1_n,s1,"max_rel_err",max_rel_err,"max_rel_warn",max_rel_warn,out.ptr())) success=false;
00348     if(verbose && dumpAll) *out << "\nneV2 =\n" << *neV2;
00349 
00350     if(verbose) *out << "\n*** (B.2) Test RTOps with two or more arguments\n";
00351 
00352     if(verbose) *out << "\nPerforming ev1 = ev2 ...\n";
00353     timer.start(true);
00354      Thyra::assign( ev1.ptr(), *ev2 );
00355     timer.stop();
00356     if(verbose) *out << "  time = " << timer.totalElapsedTime() << " sec\n";
00357     if(!testRelErr("Thyra::norm_1(ev1)",Thyra::norm_1(*ev1),"Thyra::norm_1(ev2)",ev2_nrm,"max_rel_err",max_rel_err,"max_rel_warn",max_rel_warn,out.ptr())) success=false;
00358     if(verbose && dumpAll) *out << "\nev1 =\n" << *ev1;
00359 
00360     if(verbose) *out << "\nPerforming eV1 = eV2 ...\n";
00361     timer.start(true);
00362      Thyra::assign( eV1.ptr(), *eV2 );
00363     timer.stop();
00364     if(verbose) *out << "  time = " << timer.totalElapsedTime() << " sec\n";
00365     if(!testRelErr("Thyra::norm_1(eV1)",Thyra::norm_1(*eV1),"Thyra::norm_1(eV2)",eV2_nrm,"max_rel_err",max_rel_err,"max_rel_warn",max_rel_warn,out.ptr())) success=false;
00366     if(verbose && dumpAll) *out << "\neV1 =\n" << *eV1;
00367 
00368     if(verbose) *out << "\nPerforming ev1 = nev2 ...\n";
00369     timer.start(true);
00370      Thyra::assign( ev1.ptr(), *nev2 );
00371     timer.stop();
00372     if(verbose) *out << "  time = " << timer.totalElapsedTime() << " sec\n";
00373     if(!testRelErr("Thyra::norm_1(ev1)",Thyra::norm_1(*ev1),"Thyra::norm_1(nev2)",nev2_nrm,"max_rel_err",max_rel_err,"max_rel_warn",max_rel_warn,out.ptr())) success=false;
00374     if(verbose && dumpAll) *out << "\nev1 =\n" << *ev1;
00375 
00376     if(verbose) *out << "\nPerforming nev1 = ev2 ...\n";
00377     timer.start(true);
00378      Thyra::assign( nev1.ptr(), *ev2 );
00379     timer.stop();
00380     if(verbose) *out << "  time = " << timer.totalElapsedTime() << " sec\n";
00381     if(!testRelErr("Thyra::norm_1(nev1)",Thyra::norm_1(*nev1),"Thyra::norm_1(ev2)",ev2_nrm,"max_rel_err",max_rel_err,"max_rel_warn",max_rel_warn,out.ptr())) success=false;
00382     if(verbose && dumpAll) *out << "\nnev1 =\n" << *nev1;
00383 
00384     if(verbose) *out << "\nPerforming nev1 = nev2 ...\n";
00385     timer.start(true);
00386      Thyra::assign( nev1.ptr(), *nev2 );
00387     timer.stop();
00388     if(verbose) *out << "  time = " << timer.totalElapsedTime() << " sec\n";
00389     if(!testRelErr("Thyra::norm_1(nev1)",Thyra::norm_1(*nev1),"Thyra::norm_1(nev2)",nev2_nrm,"max_rel_err",max_rel_err,"max_rel_warn",max_rel_warn,out.ptr())) success=false;
00390     if(verbose && dumpAll) *out << "\nnev1 =\n" << *nev1;
00391 
00392     if(verbose) *out << "\nPerforming eV1 = neV2 ...\n";
00393     timer.start(true);
00394      Thyra::assign( eV1.ptr(), *neV2 );
00395     timer.stop();
00396     if(verbose) *out << "  time = " << timer.totalElapsedTime() << " sec\n";
00397     if(!testRelErr("Thyra::norm_1(eV1)",Thyra::norm_1(*eV1),"Thyra::norm_1(neV2)",neV2_nrm,"max_rel_err",max_rel_err,"max_rel_warn",max_rel_warn,out.ptr())) success=false;
00398     if(verbose && dumpAll) *out << "\neV1 =\n" << *eV1;
00399 
00400     if(verbose) *out << "\nPerforming neV1 = eV2 ...\n";
00401     timer.start(true);
00402      Thyra::assign( neV1.ptr(), *eV2 );
00403     timer.stop();
00404     if(verbose) *out << "  time = " << timer.totalElapsedTime() << " sec\n";
00405     if(!testRelErr("Thyra::norm_1(neV1)",Thyra::norm_1(*neV1),"Thyra::norm_1(eV2)",eV2_nrm,"max_rel_err",max_rel_err,"max_rel_warn",max_rel_warn,out.ptr())) success=false;
00406     if(verbose && dumpAll) *out << "\nneV1 =\n" << *neV1;
00407 
00408     if(verbose) *out << "\nPerforming neV1 = neV2 ...\n";
00409     timer.start(true);
00410      Thyra::assign( neV1.ptr(), *neV2 );
00411     timer.stop();
00412     if(verbose) *out << "  time = " << timer.totalElapsedTime() << " sec\n";
00413     if(!testRelErr("Thyra::norm_1(neV1)",Thyra::norm_1(*neV1),"Thyra::norm_1(neV2)",neV2_nrm,"max_rel_err",max_rel_err,"max_rel_warn",max_rel_warn,out.ptr())) success=false;
00414     if(verbose && dumpAll) *out << "\nneV1 =\n" << *neV1;
00415 
00416     Thyra::LinearOpTester<Scalar> linearOpTester;
00417     linearOpTester.set_all_warning_tol(max_rel_warn);
00418     linearOpTester.set_all_error_tol(max_rel_err);
00419     linearOpTester.dump_all(dumpAll);
00420 
00421     Thyra::LinearOpWithSolveTester<Scalar> linearOpWithSolveTester;
00422     linearOpWithSolveTester.set_all_solve_tol(max_rel_err);
00423     linearOpWithSolveTester.set_all_slack_error_tol(max_rel_err);
00424     linearOpWithSolveTester.set_all_slack_warning_tol(max_rel_warn);
00425     linearOpWithSolveTester.dump_all(dumpAll);
00426 
00427 
00428     if(verbose) *out << "\n*** (B.3) Test Vector linear operator interface\n";
00429 
00430     if(verbose) *out << "\nChecking *out linear operator interface of ev1 ...\n";
00431     result = linearOpTester.check(*ev1,out.ptr());
00432     if(!result) success = false;
00433 
00434     if(verbose) *out << "\nChecking *out linear operator interface of nev1 ...\n";
00435     result = linearOpTester.check(*nev1,out.ptr());
00436     if(!result) success = false;
00437 
00438 
00439     if(verbose) *out << "\n*** (B.4) Test MultiVector linear operator interface\n";
00440 
00441     if(verbose) *out << "\nChecking *out linear operator interface of eV1 ...\n";
00442     result = linearOpTester.check(*eV1,out.ptr());
00443     if(!result) success = false;
00444 
00445     if(verbose) *out << "\nChecking *out linear operator interface of neV1 ...\n";
00446     result = linearOpTester.check(*neV1,out.ptr());
00447     if(!result) success = false;
00448 
00449     const std::string s2_n = "scalar^2*global_dim*num_mv_cols";
00450     const Scalar s2 = scalar*scalar*global_dim*num_mv_cols;
00451 
00452     RCP<Thyra::MultiVectorBase<Scalar> >
00453       T = createMembers(eV1->domain(),num_mv_cols);
00454 
00455 
00456     if(verbose) *out << "\n*** (B.5) Test MultiVector::apply(...)\n";
00457 
00458     if(verbose) *out << "\nPerforming eV1'*eV2 ...\n";
00459     timer.start(true);
00460     apply( *eV1, TRANS, *eV2, T.ptr() );
00461     timer.stop();
00462     if(verbose) *out << "  time = " << timer.totalElapsedTime() << " sec\n";
00463     if(!testRelErr("Thyra::norm_1(eV1'*eV2)",Thyra::norm_1(*T),s2_n,s2,"max_rel_err",max_rel_err,"max_rel_warn",max_rel_warn,out.ptr())) success=false;
00464     if(verbose && dumpAll) *out << "\neV1'*eV2 =\n" << *T;
00465 
00466     if(verbose) *out << "\nPerforming neV1'*eV2 ...\n";
00467     timer.start(true);
00468     apply( *neV1, TRANS, *eV2, T.ptr() );
00469     timer.stop();
00470     if(verbose) *out << "  time = " << timer.totalElapsedTime() << " sec\n";
00471     if(!testRelErr("Thyra::norm_1(neV1'*eV2)",Thyra::norm_1(*T),s2_n,s2,"max_rel_err",max_rel_err,"max_rel_warn",max_rel_warn,out.ptr())) success=false;
00472     if(verbose && dumpAll) *out << "\nneV1'*eV2 =\n" << *T;
00473 
00474     if(verbose) *out << "\nPerforming eV1'*neV2 ...\n";
00475     timer.start(true);
00476     apply( *eV1, TRANS, *neV2, T.ptr() );
00477     timer.stop();
00478     if(verbose) *out << "  time = " << timer.totalElapsedTime() << " sec\n";
00479     if(!testRelErr("Thyra::norm_1(eV1'*neV2)",Thyra::norm_1(*T),s2_n,s2,"max_rel_err",max_rel_err,"max_rel_warn",max_rel_warn,out.ptr())) success=false;
00480     if(verbose && dumpAll) *out << "\neV1'*neV2 =\n" << *T;
00481 
00482     if(verbose) *out << "\nPerforming neV1'*neV2 ...\n";
00483     timer.start(true);
00484     apply( *neV1, TRANS, *neV2, T.ptr() );
00485     timer.stop();
00486     if(verbose) *out << "  time = " << timer.totalElapsedTime() << " sec\n";
00487     if(!testRelErr("Thyra::norm_1(neV1'*neV2)",Thyra::norm_1(*T),s2_n,s2,"max_rel_err",max_rel_err,"max_rel_warn",max_rel_warn,out.ptr())) success=false;
00488     if(verbose && dumpAll) *out << "\nneV1'*neV2 =\n" << *T;
00489 
00490 
00491     if(verbose) *out << "\n*** (B.6) Creating a diagonal Epetra_Operator Op\n";
00492 
00493     RCP<Epetra_Operator>  epetra_op;
00494 
00495     {
00496       // Create a diagonal matrix with scalar on the diagonal
00497       RCP<Epetra_CrsMatrix>
00498         epetra_mat = rcp(new Epetra_CrsMatrix(::Copy,*epetra_map,1));
00499       Scalar values[1] = { scalar };
00500       int indices[1];
00501       const int IB = epetra_map->IndexBase(), offset = procRank*local_dim;
00502       for( int k = 0; k < local_dim; ++k ) {
00503         indices[0] = offset + k + IB;  // global column
00504         epetra_mat->InsertGlobalValues(
00505           offset + k + IB     // GlobalRow
00506           ,1                  // NumEntries
00507           ,values             // Values
00508           ,indices            // Indices
00509           );
00510       }
00511       epetra_mat->FillComplete();
00512       epetra_op = epetra_mat;
00513     } // end epetra_op
00514 
00515     RCP<const Thyra::LinearOpBase<Scalar> >
00516       Op = Thyra::epetraLinearOp(epetra_op);
00517 
00518     if(verbose && dumpAll) *out << "\nOp=\n" << *Op;
00519 
00520 
00521     if(verbose) *out << "\n*** (B.6b) Going through partial then full initialization of EpetraLinearOp ...\n";
00522 
00523     {
00524       if(verbose) *out
00525         << "\nChecking isFullyUninitialized(*nonconstEpetraLinearOp())==true : ";
00526       RCP<Thyra::EpetraLinearOp>
00527         thyraEpetraOp = Thyra::nonconstEpetraLinearOp();
00528       result = isFullyUninitialized(*thyraEpetraOp);
00529       if(!result) success = false;
00530       if(verbose) *out << Thyra::passfail(result) << endl;
00531     }
00532 
00533     {
00534 
00535       if(verbose) *out
00536         << "\nthyraEpetraOp = partialNonconstEpetraLinearOp(...)\n";
00537       RCP<Thyra::EpetraLinearOp> thyraEpetraOp =
00538         Thyra::partialNonconstEpetraLinearOp(
00539           epetra_vs, epetra_vs, epetra_op
00540           );
00541 
00542       if(verbose) *out
00543         << "\nChecking isPartiallyInitialized(*thyraEpetraOp)==true : ";
00544       result = isPartiallyInitialized(*thyraEpetraOp);
00545       if(!result) success = false;
00546       if(verbose) *out << Thyra::passfail(result) << endl;
00547 
00548       if(verbose) *out
00549         << "\nthyraEpetraOp->setFullyInitialized(true)\n";
00550       thyraEpetraOp->setFullyInitialized(true);
00551 
00552       if(verbose) *out
00553         << "\nChecking isFullyInitialized(*thyraEpetraOp)==true : ";
00554       result = isFullyInitialized(*thyraEpetraOp);
00555       if(!result) success = false;
00556       if(verbose) *out << Thyra::passfail(result) << endl;
00557 
00558     }
00559 
00560 
00561     if(verbose) *out << "\n*** (B.7) Test EpetraLinearOp linear operator interface\n";
00562 
00563     if(verbose) *out << "\nChecking *out linear operator interface of Op ...\n";
00564     result = linearOpTester.check(*Op,out.ptr());
00565     if(!result) success = false;
00566 
00567     RCP<Thyra::VectorBase<Scalar> >
00568       ey  = createMember(epetra_vs);
00569     RCP<Thyra::MultiVectorBase<Scalar> >
00570       eY  = createMembers(epetra_vs,num_mv_cols);
00571     RCP<Thyra::VectorBase<Scalar> >
00572       ney = createMember(non_epetra_vs);
00573     RCP<Thyra::MultiVectorBase<Scalar> >
00574       neY = createMembers(non_epetra_vs,num_mv_cols);
00575 
00576 
00577     if(verbose) *out << "\n*** (B.8) Mix and match vector and Multi-vectors with Epetra opeator\n";
00578 
00579     const std::string s3_n = "2*scalar^2*global_dim";
00580     const Scalar s3 = 2*scalar*scalar*global_dim;
00581     
00582     if(verbose) *out << "\nPerforming ey = 2*Op*ev1 ...\n";
00583     timer.start(true);
00584     apply( *Op, NOTRANS, *ev1, ey.ptr(), 2.0 );
00585     timer.stop();
00586     if(verbose) *out << "  time = " << timer.totalElapsedTime() << " sec\n";
00587     if(!testRelErr("Thyra::norm_1(ey)",Thyra::norm_1(*ey),s3_n,s3,"max_rel_err",max_rel_err,"max_rel_warn",max_rel_warn,out.ptr())) success=false;
00588 
00589     if(verbose) *out << "\nPerforming eY = 2*Op*eV1 ...\n";
00590     timer.start(true);
00591     apply( *Op, NOTRANS, *eV1, eY.ptr(), 2.0 );
00592     timer.stop();
00593     if(verbose) *out << "  time = " << timer.totalElapsedTime() << " sec\n";
00594     if(!testRelErr("Thyra::norm_1(eY)",Thyra::norm_1(*eY),s3_n,s3,"max_rel_err",max_rel_err,"max_rel_warn",max_rel_warn,out.ptr())) success=false;
00595 
00596     if(verbose) *out << "\nPerforming ney = 2*Op*ev1 ...\n";
00597     timer.start(true);
00598     apply( *Op, NOTRANS, *ev1, ney.ptr(), 2.0 );
00599     timer.stop();
00600     if(verbose) *out << "  time = " << timer.totalElapsedTime() << " sec\n";
00601     if(!testRelErr("Thyra::norm_1(ney)",Thyra::norm_1(*ney),s3_n,s3,"max_rel_err",max_rel_err,"max_rel_warn",max_rel_warn,out.ptr())) success=false;
00602 
00603     if(verbose) *out << "\nPerforming neY = 2*Op*eV1 ...\n";
00604     timer.start(true);
00605     apply( *Op, NOTRANS, *eV1, neY.ptr(), 2.0 );
00606     timer.stop();
00607     if(verbose) *out << "  time = " << timer.totalElapsedTime() << " sec\n";
00608     if(!testRelErr("Thyra::norm_1(neY)",Thyra::norm_1(*neY),s3_n,s3,"max_rel_err",max_rel_err,"max_rel_warn",max_rel_warn,out.ptr())) success=false;
00609 
00610     if(verbose) *out << "\nPerforming ey = 2*Op*nev1 ...\n";
00611     timer.start(true);
00612     apply( *Op, NOTRANS, *nev1, ey.ptr(), 2.0 );
00613     timer.stop();
00614     if(verbose) *out << "  time = " << timer.totalElapsedTime() << " sec\n";
00615     if(!testRelErr("Thyra::norm_1(ey)",Thyra::norm_1(*ey),s3_n,s3,"max_rel_err",max_rel_err,"max_rel_warn",max_rel_warn,out.ptr())) success=false;
00616 
00617     if(verbose) *out << "\nPerforming eY = 2*Op*neV1 ...\n";
00618     timer.start(true);
00619     apply( *Op, NOTRANS, *neV1, eY.ptr(), 2.0 );
00620     timer.stop();
00621     if(verbose) *out << "  time = " << timer.totalElapsedTime() << " sec\n";
00622     if(!testRelErr("Thyra::norm_1(eY)",Thyra::norm_1(*eY),s3_n,s3,"max_rel_err",max_rel_err,"max_rel_warn",max_rel_warn,out.ptr())) success=false;
00623 
00624     if(verbose) *out << "\nPerforming ney = 2*Op*nev1 ...\n";
00625     timer.start(true);
00626     apply( *Op, NOTRANS, *nev1, ney.ptr(), 2.0 );
00627     timer.stop();
00628     if(verbose) *out << "  time = " << timer.totalElapsedTime() << " sec\n";
00629     if(!testRelErr("Thyra::norm_1(ney)",Thyra::norm_1(*ney),s3_n,s3,"max_rel_err",max_rel_err,"max_rel_warn",max_rel_warn,out.ptr())) success=false;
00630 
00631     if(verbose) *out << "\nPerforming ney = 2*Op*nev1 through MultiVector interface ...\n";
00632     timer.start(true);
00633     apply( *Op, NOTRANS, static_cast<const Thyra::MultiVectorBase<Scalar>&>(*nev1), Ptr<Thyra::MultiVectorBase<Scalar> >(ney.ptr()), 2.0 );
00634     timer.stop();
00635     if(verbose) *out << "  time = " << timer.totalElapsedTime() << " sec\n";
00636     if(!testRelErr("Thyra::norm_1(ney)",Thyra::norm_1(*ney),s3_n,s3,"max_rel_err",max_rel_err,"max_rel_warn",max_rel_warn,out.ptr())) success=false;
00637 
00638     if(verbose) *out << "\nPerforming neY = 2*Op*neV1 ...\n";
00639     timer.start(true);
00640     apply( *Op, NOTRANS, *neV1, neY.ptr(), 2.0 );
00641     timer.stop();
00642     if(verbose) *out << "  time = " << timer.totalElapsedTime() << " sec\n";
00643     if(!testRelErr("Thyra::norm_1(neY)",Thyra::norm_1(*neY),s3_n,s3,"max_rel_err",max_rel_err,"max_rel_warn",max_rel_warn,out.ptr())) success=false;
00644 
00645 
00646     if(verbose) *out << "\n*** (B.9) Testing Multi-vector views with Epetra operator\n";
00647 
00648     const Thyra::Range1D col_rng(0,1);
00649     const Teuchos::Tuple<int, 2> cols = Teuchos::tuple<int>(2, 3);
00650 
00651     RCP<const Thyra::MultiVectorBase<Scalar> >
00652       eV1_v1  = rcp_static_cast<const Thyra::MultiVectorBase<Scalar> >(eV1)->subView(col_rng),
00653       eV1_v2  = rcp_static_cast<const Thyra::MultiVectorBase<Scalar> >(eV1)->subView(cols);
00654     RCP<const Thyra::MultiVectorBase<Scalar> >
00655       neV1_v1  = rcp_static_cast<const Thyra::MultiVectorBase<Scalar> >(neV1)->subView(col_rng),
00656       neV1_v2  = rcp_static_cast<const Thyra::MultiVectorBase<Scalar> >(neV1)->subView(cols);
00657     if(verbose && dumpAll) {
00658       *out << "\neV1_v1=\n" << *eV1_v1;
00659       *out << "\neV1_v2=\n" << *eV1_v2;
00660       *out << "\nneV1_v1=\n" << *neV1_v1;
00661       *out << "\nneV1_v2=\n" << *neV1_v2;
00662     }
00663 
00664     if(verbose) *out << "\nPerforming eY_v1 = 2*Op*eV1_v1 ...\n";
00665     timer.start(true);
00666     apply( *Op, NOTRANS, *eV1_v1, eY->subView(col_rng).ptr(), 2.0 );
00667     timer.stop();
00668     if(verbose) *out << "  time = " << timer.totalElapsedTime() << " sec\n";
00669     if(verbose && dumpAll) *out << "\neV_v1=\n" << *eY->subView(col_rng);
00670     if(!testRelErr("Thyra::norm_1(eY_v1)",Thyra::norm_1(*eY->subView(col_rng)),s3_n,s3,"max_rel_err",max_rel_err,"max_rel_warn",max_rel_warn,out.ptr())) success=false;
00671 
00672     if(verbose) *out << "\nPerforming eY_v2 = 2*Op*eV1_v2 ...\n";
00673     timer.start(true);
00674     apply( *Op, NOTRANS, *eV1_v2, eY->subView(cols).ptr(), 2.0 );
00675     timer.stop();
00676     if(verbose) *out << "  time = " << timer.totalElapsedTime() << " sec\n";
00677     if(verbose && dumpAll) *out << "\neV_v2=\n" << *eY->subView(cols);
00678     if(!testRelErr("Thyra::norm_1(eY_v2)",Thyra::norm_1(*eY->subView(cols)),s3_n,s3,"max_rel_err",max_rel_err,"max_rel_warn",max_rel_warn,out.ptr())) success=false;
00679 
00680     if(verbose) *out << "\nPerforming neY_v1 = 2*Op*eV1_v1 ...\n";
00681     timer.start(true);
00682     apply( *Op, NOTRANS, *eV1_v1, neY->subView(col_rng).ptr(), 2.0 );
00683     timer.stop();
00684     if(verbose) *out << "  time = " << timer.totalElapsedTime() << " sec\n";
00685     if(verbose && dumpAll) *out << "\neV_v1=\n" << *eY->subView(col_rng);
00686     if(!testRelErr("Thyra::norm_1(neY_v1)",Thyra::norm_1(*neY->subView(col_rng)),s3_n,s3,"max_rel_err",max_rel_err,"max_rel_warn",max_rel_warn,out.ptr())) success=false;
00687 
00688     if(verbose) *out << "\nPerforming eY_v1 = 2*Op*neV1_v1 ...\n";
00689     timer.start(true);
00690     apply( *Op, NOTRANS, *neV1_v1, eY->subView(col_rng).ptr(), 2.0 );
00691     timer.stop();
00692     if(verbose) *out << "  time = " << timer.totalElapsedTime() << " sec\n";
00693     if(!testRelErr("Thyra::norm_1(eY_v1)",Thyra::norm_1(*eY->subView(col_rng)),s3_n,s3,"max_rel_err",max_rel_err,"max_rel_warn",max_rel_warn,out.ptr())) success=false;
00694 
00695     if(verbose) *out << "\nPerforming neY_v2 = 2*Op*eV1_v2 ...\n";
00696     timer.start(true);
00697     apply( *Op, NOTRANS, *eV1_v2, neY->subView(cols).ptr(), 2.0 );
00698     timer.stop();
00699     if(verbose) *out << "  time = " << timer.totalElapsedTime() << " sec\n";
00700     if(verbose && dumpAll) *out << "\neV_v2=\n" << *eY->subView(cols);
00701     if(!testRelErr("Thyra::norm_1(neY_v2)",Thyra::norm_1(*neY->subView(cols)),s3_n,s3,"max_rel_err",max_rel_err,"max_rel_warn",max_rel_warn,out.ptr())) success=false;
00702 
00703     if(verbose) *out << "\nPerforming eY_v2 = 2*Op*neV1_v2 ...\n";
00704     timer.start(true);
00705     apply( *Op, NOTRANS, *neV1_v2, eY->subView(cols).ptr(), 2.0 );
00706     timer.stop();
00707     if(verbose) *out << "  time = " << timer.totalElapsedTime() << " sec\n";
00708     if(verbose && dumpAll) *out << "\neV_v2=\n" << *eY->subView(cols);
00709     if(!testRelErr("Thyra::norm_1(eY_v2)",Thyra::norm_1(*eY->subView(cols)),s3_n,s3,"max_rel_err",max_rel_err,"max_rel_warn",max_rel_warn,out.ptr())) success=false;
00710 
00711 
00712     if(verbose) *out << "\n*** (B.10) Testing Vector and MultiVector view creation functions\n";
00713 
00714     {
00715 
00716       const std::string s_n = "fabs(scalar)*num_mv_cols";
00717       const Scalar s = fabs(scalar)*num_mv_cols;
00718 
00719       Array<Scalar> t_raw_values( num_mv_cols );
00720       RTOpPack::SubVectorView<Scalar> t_raw( 0, num_mv_cols,
00721         arcpFromArray(t_raw_values), 1 );
00722 
00723       std::fill_n( t_raw_values.begin(), t_raw_values.size(), ST::zero() );
00724       Thyra::assign( createMemberView(T->range(),t_raw).ptr(), scalar );
00725       Teuchos::RCP<const Thyra::VectorBase<Scalar> > t_view = createMemberView(T->range(),static_cast<RTOpPack::ConstSubVectorView<Scalar>&>(t_raw));
00726       Scalar t_nrm = Thyra::norm_1(*t_view);
00727       if(!testRelErr("Thyra::norm_1(t_view)",t_nrm,s_n,s,"max_rel_err",max_rel_err,"max_rel_warn",max_rel_warn,out.ptr())) success=false;
00728       if(verbose && dumpAll) *out << "\nt_view =\n" << *t_view;
00729 
00730 /*
00731 #ifndef SUN_CXX // The sun compiler Forte Developer 5.4 does not destory temporaries properly and this does not work
00732       std::fill_n( t_raw_values.begin(), t_raw_values.size(), ST::zero() );
00733       Thyra::assign( T->range().ptr()->Thyra::VectorSpaceBase<Scalar>::createMemberView(t_raw), scalar );
00734       t_view = T->range()->Thyra::VectorSpaceBase<Scalar>::createMemberView(static_cast<RTOpPack::ConstSubVectorView<Scalar>&>(t_raw));
00735       t_nrm = Thyra::norm_1(*t_view);
00736       if(!testRelErr("Thyra::norm_1(t_view)",t_nrm,s_n,s,"max_rel_err",max_rel_err,"max_rel_warn",max_rel_warn,out.ptr())) success=false;
00737       if(verbose && dumpAll) *out << "\nt_view =\n" << *t_view;
00738 #endif
00739 */
00740 
00741       Array<Scalar> T_raw_values( num_mv_cols * num_mv_cols );
00742       RTOpPack::SubMultiVectorView<Scalar> T_raw( 0, num_mv_cols, 0, num_mv_cols,
00743         arcpFromArray(T_raw_values), num_mv_cols );
00744 
00745       std::fill_n( T_raw_values.begin(), T_raw_values.size(), ST::zero() );
00746       Thyra::assign( createMembersView(T->range(),T_raw).ptr(), scalar );
00747       Teuchos::RCP<const Thyra::MultiVectorBase<Scalar> >
00748         T_view = createMembersView(T->range(),static_cast<RTOpPack::ConstSubMultiVectorView<Scalar>&>(T_raw));
00749       Scalar T_nrm = Thyra::norm_1(*T_view);
00750       if(!testRelErr("Thyra::norm_1(T_view)",T_nrm,s_n,s,"max_rel_err",max_rel_err,"max_rel_warn",max_rel_warn,out.ptr())) success=false;
00751       if(verbose && dumpAll) *out << "\nT_view =\n" << *T_view;
00752 
00753 /*
00754 #ifndef SUN_CXX // The sun compiler Forte Developer 5.4 does not destory temporaries properly and this does not work
00755       std::fill_n( T_raw_values.begin(), T_raw_values.size(), ST::zero() );
00756       Thyra::assign( T->range().ptr()->Thyra::VectorSpaceBase<Scalar>::createMembersView(T_raw), scalar );
00757       T_view = T->range()->Thyra::VectorSpaceBase<Scalar>::createMembersView(static_cast<RTOpPack::ConstSubMultiVectorView<Scalar>&>(T_raw));
00758       T_nrm = Thyra::norm_1(*T_view);
00759       if(!testRelErr("Thyra::norm_1(T_view)",T_nrm,s_n,s,"max_rel_err",max_rel_err,"max_rel_warn",max_rel_warn,out.ptr())) success=false;
00760       if(verbose && dumpAll) *out << "\nT_view =\n" << *T_view;
00761 #endif
00762 */
00763 
00764     }
00765 
00766 
00767     if(verbose) *out << "\n*** (B.11) Testing Epetra_Vector and Epetra_MultiVector wrappers\n";
00768 
00769     {
00770 
00771       Teuchos::RCP<const Thyra::SpmdVectorSpaceBase<Scalar> >
00772         mpi_vs = Teuchos::rcp_dynamic_cast<const Thyra::SpmdVectorSpaceBase<Scalar> >(epetra_vs,true);
00773 
00774       if(verbose) *out << "\nCreating temporary Epetra_Vector et1 and Epetra_MultiVector eT1 objects ...\n";
00775       Teuchos::RCP<Epetra_Vector>
00776         et1 = Teuchos::rcp(new Epetra_Vector(*epetra_map));
00777       Teuchos::RCP<Epetra_MultiVector>
00778         eT1 = Teuchos::rcp(new Epetra_MultiVector(*epetra_map,num_mv_cols));
00779 
00780       if(verbose) *out << "\nCreating non-const VectorBase t1 and MultiVectorBase T1 objects from et1 and eT2 ...\n";
00781       Teuchos::RCP<Thyra::VectorBase<Scalar> >
00782         t1 = create_Vector(et1,mpi_vs);
00783       Teuchos::RCP<Thyra::MultiVectorBase<Scalar> >
00784         T1 = create_MultiVector(eT1,mpi_vs);
00785 
00786       if(verbose) *out << "\nPerforming t1 = ev1 ...\n";
00787       assign( t1.ptr(), *ev1 );
00788       if(!testRelErr(
00789            "sum(t1)",Thyra::sum(*t1),"sum(ev1)",sum(*ev1)
00790            ,"max_rel_err",max_rel_err,"max_rel_warn",max_rel_warn,out.ptr())
00791         ) success=false;
00792 
00793       if(verbose) *out << "\nPerforming T1 = eV1 ...\n";
00794       assign( T1.ptr(), *eV1 );
00795       if(!testRelErr(
00796            "norm_1(T1)",Thyra::norm_1(*T1),"norm_1(eV1)",norm_1(*eV1)
00797            ,"max_rel_err",max_rel_err,"max_rel_warn",max_rel_warn,out.ptr())
00798         ) success=false;
00799 
00800       if(verbose) *out << "\nChecking that t1 and T1 yield the same objects as et1 and eT2 ...\n";
00801   
00802       Teuchos::RCP<Epetra_Vector>
00803         et1_v = get_Epetra_Vector(*epetra_map,t1);
00804       result = et1_v.get() == et1.get();
00805       if(verbose) *out << "\net1_v.get() = " << et1_v.get() << " == et1.get() = " << et1.get() << " : " << passfail(result) << endl;
00806       if(!result) success = false;
00807 
00808       Teuchos::RCP<Epetra_MultiVector>
00809         eT1_v = get_Epetra_MultiVector(*epetra_map,T1);
00810       result = eT1_v.get() == eT1.get();
00811       if(verbose) *out << "\neT1_v.get() = " << eT1_v.get() << " == eT1.get() = " << eT1.get() << " : " << passfail(result) << endl;
00812       if(!result) success = false;
00813 
00814       if(verbose) *out << "\nCreating const VectorBase ct1 and MultiVectorBase cT1 objects from et1 and eT2 ...\n";
00815       Teuchos::RCP<const Thyra::VectorBase<Scalar> >
00816         ct1 = create_Vector(Teuchos::rcp_implicit_cast<const Epetra_Vector>(et1),mpi_vs);
00817       Teuchos::RCP<const Thyra::MultiVectorBase<Scalar> >
00818         cT1 = create_MultiVector(Teuchos::rcp_implicit_cast<const Epetra_MultiVector>(eT1),mpi_vs);
00819 
00820       if(!testRelErr(
00821            "sum(ct1)",Thyra::sum(*ct1),"sum(ev1)",sum(*ev1)
00822            ,"max_rel_err",max_rel_err,"max_rel_warn",max_rel_warn,out.ptr())
00823         ) success=false;
00824 
00825       if(!testRelErr(
00826            "norm_1(cT1)",Thyra::norm_1(*cT1),"norm_1(eV1)",norm_1(*eV1)
00827            ,"max_rel_err",max_rel_err,"max_rel_warn",max_rel_warn,out.ptr())
00828         ) success=false;
00829 
00830       if(verbose) *out << "\nChecking that ct1 and cT1 yield the same objects as et1 and eT2 ...\n";
00831   
00832       Teuchos::RCP<const Epetra_Vector>
00833         cet1_v = get_Epetra_Vector(*epetra_map,ct1);
00834       result = cet1_v.get() == et1.get();
00835       if(verbose) *out << "\ncet1_v.get() = " << cet1_v.get() << " == et1.get() = " << et1.get() << " : " << passfail(result) << endl;
00836       if(!result) success = false;
00837 
00838       Teuchos::RCP<const Epetra_MultiVector>
00839         ceT1_v = get_Epetra_MultiVector(*epetra_map,cT1);
00840       result = ceT1_v.get() == eT1.get();
00841       if(verbose) *out << "\nceT1_v.get() = " << ceT1_v.get() << " == eT1.get() = " << eT1.get() << " : " << passfail(result) << endl;
00842       if(!result) success = false;
00843 
00844       if(verbose) *out << "\nCreating non-const Epetra_Vector ett1 and Epetra_MultiVector etT1 objects from clones of t1 and T2 ...\n";
00845       Teuchos::RCP<Epetra_Vector>
00846         ett1 = get_Epetra_Vector(*epetra_map,t1->clone_v());
00847       Teuchos::RCP<Epetra_MultiVector>
00848         etT1 = get_Epetra_MultiVector(*epetra_map,T1->clone_mv());
00849 
00850       if(verbose) *out << "\nChecking that ett1 and etT1 yield objects with the same value as et1 and eT2 ...\n";
00851 
00852       if(!testRelErr(
00853            "sum(ett1)",sum(*ett1),"sum(et1)",sum(*et1)
00854            ,"max_rel_err",max_rel_err,"max_rel_warn",max_rel_warn,out.ptr())
00855         ) success=false;
00856 
00857       if(!testRelErr(
00858            "sum(etT1)",sum(*etT1),"sum(eT1)",sum(*eT1)
00859            ,"max_rel_err",max_rel_err,"max_rel_warn",max_rel_warn,out.ptr())
00860         ) success=false;
00861 
00862       if(verbose) *out << "\nCreating const Epetra_Vector cett1 and Epetra_MultiVector cetT1 objects from clones of t1 and T2 ...\n";
00863       Teuchos::RCP<const Epetra_Vector>
00864         cett1 = get_Epetra_Vector(*epetra_map,Teuchos::rcp_implicit_cast<const Thyra::VectorBase<Scalar> >(t1->clone_v()));
00865       Teuchos::RCP<const Epetra_MultiVector>
00866         cetT1 = get_Epetra_MultiVector(*epetra_map,Teuchos::rcp_implicit_cast<const Thyra::MultiVectorBase<Scalar> >(T1->clone_mv()));
00867 
00868       if(verbose) *out << "\nChecking that cett1 and cetT1 yield objects with the same value as et1 and eT2 ...\n";
00869 
00870       if(!testRelErr(
00871            "sum(cett1)",sum(*cett1),"sum(et1)",sum(*et1)
00872            ,"max_rel_err",max_rel_err,"max_rel_warn",max_rel_warn,out.ptr())
00873         ) success=false;
00874 
00875       if(!testRelErr(
00876            "sum(cetT1)",sum(*cetT1),"sum(eT1)",sum(*eT1)
00877            ,"max_rel_err",max_rel_err,"max_rel_warn",max_rel_warn,out.ptr())
00878         ) success=false;
00879 
00880     }
00881 
00882 
00883     if(verbose) *out << "\n*** (B.12) Test DiagonalEpetraLinearOpWithSolveFactory \n";
00884 
00885     {
00886 
00887       if(verbose) *out << "\nUsing DiagonalEpetraLinearOpWithSolveFactory to create diagLOWS from Op ...\n";
00888       
00889       Thyra::DiagonalEpetraLinearOpWithSolveFactory diagLOWSFactory;
00890 
00891       Teuchos::RCP<Thyra::LinearOpWithSolveBase<Scalar> >
00892         diagLOWS = Thyra::linearOpWithSolve<double>(diagLOWSFactory,Op);
00893 
00894       if(verbose) *out << "\nTesting LinearOpBase interface of diagLOWS ...\n";
00895 
00896       result = linearOpTester.check(*diagLOWS, out.ptr());
00897       if(!result) success = false;
00898 
00899       if(verbose) *out << "\nTesting LinearOpWithSolveBase interface of diagLOWS ...\n";
00900 
00901       result = linearOpWithSolveTester.check(*diagLOWS, &*out);
00902       if(!result) success = false;
00903     
00904     }
00905 
00906 
00907     if(verbose)
00908       *out
00909         << "\n***"
00910         << "\n*** (C) Comparing the speed of Thyra adapted Epetra objects verses raw Epetra objects"
00911         << "\n***\n";
00912 
00913     //
00914     // Setup the number of timing loops to get good timings
00915     //
00916     // Here we try to shoot for timing ab*out 1 second's worth of
00917     // computations and adjust the number of evaluation loops
00918     // accordingly.  Let X be the approximate number of flops per
00919     // loop (or per evaluation).  We then compute the number of
00920     // loops as:
00921     //
00922     // 1.0 sec |  num CPU flops |   1 loop  |
00923     // --------|----------------|-----------|
00924     //         |       sec      |   X flops |
00925     //
00926     // This just comes *out to be:
00927     //
00928     //   num_time_loops_X =  max_flop_rate / (X flops per loop)
00929     //
00930     // In this computation we ignore extra overhead that will be
00931     // an issue when local_dim is small.
00932     //
00933 
00934     double raw_epetra_time, thyra_wrapped_time;
00935 
00936     
00937     if(verbose) *out << "\n*** (C.1) Comparing the speed of RTOp verses raw Epetra_Vector operations\n";
00938 
00939     const double flop_adjust_factor_1 = 3.0;
00940     const int num_time_loops_1 = int( max_flop_rate / ( flop_adjust_factor_1 * local_dim * num_mv_cols ) ) + 1;
00941 
00942     {
00943         
00944       // Get references to Epetra_MultiVector objects in eV1 and eV2
00945       const RCP<Epetra_MultiVector>       eeV1 = get_Epetra_MultiVector(*epetra_map,eV1);
00946       const RCP<const Epetra_MultiVector> eeV2 = get_Epetra_MultiVector(*epetra_map,eV2);
00947       
00948       if(verbose)
00949         *out << "\nPerforming eeV1 = eeV2 (using raw Epetra_MultiVector::operator=(...)) " << num_time_loops_1 << " times ...\n";
00950       timer.start(true);
00951       for(int k = 0; k < num_time_loops_1; ++k ) {
00952         *eeV1 = *eeV2;
00953       }
00954       timer.stop();
00955       raw_epetra_time = timer.totalElapsedTime();
00956       if(verbose) *out << "  total time = " << raw_epetra_time << " sec\n";
00957 
00958       // When this block ends and eeV1 goes *out of scope then eV1 is guaranteed to be updated!
00959     }
00960     
00961     if(verbose)
00962       *out << "\nPerforming eV1 = eV2 (using Thyra::SpmdMultiVectorBase::applyOp(...)) " << num_time_loops_1 << " times ...\n";
00963     timer.start(true);
00964     for(int k = 0; k < num_time_loops_1; ++k ) {
00965       Thyra::assign( eV1.ptr(), *eV2 );
00966     }
00967     timer.stop();
00968     thyra_wrapped_time = timer.totalElapsedTime();
00969     if(verbose) *out << "  total time = " << thyra_wrapped_time << " sec\n";
00970     
00971     print_performance_stats( num_time_loops_1, raw_epetra_time, thyra_wrapped_time, verbose, *out );
00972 
00973     // RAB: 2004/01/05: Note, the above relative performance is likely
00974     // to be the worst of all of the others since RTOp operators are
00975     // applied seperately column by column but the relative
00976     // performance should go to about 1.0 when local_dim is
00977     // sufficiently large!  However, because
00978     // Epetra_MultiVector::Thyra::Assign(...) is implemented using double
00979     // pointer indexing, the RTOp implementation used with the Thyra
00980     // adapters is actually faster in some cases.  However, the extra
00981     // overhead of RTOp is much worse for very very small (order 10)
00982     // sizes.
00983 
00984 
00985     if(verbose)
00986       *out
00987         << "\n*** (C.2) Comparing Thyra::SpmdMultiVectorBase::apply() verses raw Epetra_MultiVector::Multiply()\n";
00988 
00989     Teuchos::TimeMonitor::zeroOutTimers();
00990 
00991     const double flop_adjust_factor_2 = 2.0;
00992     const int num_time_loops_2 = int( max_flop_rate / ( flop_adjust_factor_2* local_dim * num_mv_cols * num_mv_cols ) ) + 1;
00993 
00994     {
00995       
00996       // Get constant references to Epetra_MultiVector objects in eV1 and eV2
00997       const RCP<const Epetra_MultiVector> eeV1 = get_Epetra_MultiVector(*epetra_map,eV1);
00998       const RCP<const Epetra_MultiVector> eeV2 = get_Epetra_MultiVector(*epetra_map,eV2);
00999       
01000       Epetra_LocalMap eT_map((int) T->range()->dim(),0,*epetra_comm);
01001       Epetra_MultiVector eT(eT_map,T->domain()->dim());
01002       
01003       if(verbose)
01004         *out << "\nPerforming eeV1'*eeV2 (using raw Epetra_MultiVector::Multiply(...)) "  << num_time_loops_2 << " times ...\n";
01005       timer.start(true);
01006       for(int k = 0; k < num_time_loops_2; ++k ) {
01007         eT.Multiply( 'T', 'N', 1.0, *eeV1, *eeV2, 0.0 );
01008       }
01009       timer.stop();
01010       raw_epetra_time = timer.totalElapsedTime();
01011       if(verbose) *out << "  total time = " << raw_epetra_time << " sec\n";
01012       
01013     }
01014     
01015     if(verbose)
01016       *out << "\nPerforming eV1'*eV2 (using Thyra::SpmdMultiVectorBase::apply(...)) " << num_time_loops_2 << " times ...\n";
01017     timer.start(true);
01018     for(int k = 0; k < num_time_loops_2; ++k ) {
01019       apply( *eV1, TRANS, *eV2, T.ptr() );
01020     }
01021     timer.stop();
01022     thyra_wrapped_time = timer.totalElapsedTime();
01023     if(verbose) *out << "  total time = " << thyra_wrapped_time << " sec\n";
01024   
01025     print_performance_stats( num_time_loops_2, raw_epetra_time, thyra_wrapped_time, verbose, *out );
01026     
01027     // RAB: 2004/01/05: Note, even though the Thyra adapter does
01028     // not actually call Epetra_MultiVector::Multiply(...), the
01029     // implementation in Thyra::SpmdMultiVectorBase::apply(...)
01030     // performs almost exactly the same flops and calls dgemm(...)
01031     // as well.  Herefore, except for some small overhead, the raw
01032     // Epetra and the Thyra wrapped computations should give
01033     // almost identical times in almost all cases.
01034 
01035 
01036     if(verbose) *out << "\n*** (C.3) Comparing Thyra::EpetraLinearOp::apply() verses raw Epetra_Operator::apply()\n";
01037 
01038     Teuchos::TimeMonitor::zeroOutTimers();
01039 
01040     const double flop_adjust_factor_3 = 10.0; // lots of indirect addressing
01041     const int num_time_loops_3 = int( max_flop_rate / ( flop_adjust_factor_3 * local_dim * num_mv_cols ) ) + 1;
01042 
01043     {
01044       
01045       // Get constant references to Epetra_MultiVector objects in eV1 and eV2
01046       const RCP<const Epetra_MultiVector> eeV1 = get_Epetra_MultiVector(*epetra_map,eV1);
01047       const RCP<Epetra_MultiVector>       eeY  = get_Epetra_MultiVector(*epetra_map,eY);
01048       
01049       if(verbose)
01050         *out << "\nPerforming eeY = eOp*eeV1 (using raw Epetra_Operator::apply(...)) " << num_time_loops_3 << " times ...\n";
01051 
01052       Teuchos::TimeMonitor::zeroOutTimers();
01053 
01054       timer.start(true);
01055       epetra_op->SetUseTranspose(false);
01056       for(int k = 0; k < num_time_loops_3; ++k ) {
01057         epetra_op->Apply( *eeV1, *eeY );
01058         //eeY->Scale(2.0);
01059       }
01060       timer.stop();
01061 
01062       raw_epetra_time = timer.totalElapsedTime();
01063       if(verbose) *out << "  total time = " << raw_epetra_time << " sec\n";
01064 
01065       if(verbose)
01066         Teuchos::TimeMonitor::summarize(*out << "\n");
01067       
01068     }
01069     
01070     if(verbose)
01071       *out << "\nPerforming eY = Op*eV1 (using Thyra::EpetraLinearOp::apply(...)) " << num_time_loops_3 << " times ...\n";
01072 
01073     Teuchos::TimeMonitor::zeroOutTimers();
01074 
01075     timer.start(true);
01076     for(int k = 0; k < num_time_loops_3; ++k ) {
01077       apply( *Op, NOTRANS, *eV1, eY.ptr() );
01078     }
01079     timer.stop();
01080 
01081     thyra_wrapped_time = timer.totalElapsedTime();
01082     if(verbose) *out << "  total time = " << thyra_wrapped_time << " sec\n";
01083 
01084     if(verbose)
01085       Teuchos::TimeMonitor::summarize(*out << "\n");
01086     
01087     print_performance_stats( num_time_loops_3, raw_epetra_time, thyra_wrapped_time, verbose, *out );
01088 
01089     // RAB: 2004/01/05: Note, the above Epetra adapter is a true
01090     // adapter and simply calls Epetra_Operator::apply(...) so, except
01091     // for some small overhead, the raw Epetra and the Thyra wrapped
01092     // computations should give ab*out exactly the same runtime for
01093     // almost all cases.
01094 
01095   } // end try
01096   TEUCHOS_STANDARD_CATCH_STATEMENTS(true,*out,success)
01097 
01098   if(verbose) {
01099     if(success) *out << "\nCongratulations! All of the tests seem to have run sucessfully!\n";
01100     else        *out << "\nOh no! at least one of the tests did not check out!\n";
01101   }
01102 
01103   return (success ? 0 : -1);
01104 
01105 } // end main()
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines