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