test_epetra_adapters.cpp

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

Generated on Thu Sep 18 12:37:58 2008 for Epetra Package Browser (Single Doxygen Collection) by doxygen 1.3.9.1