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