|
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 // 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()
1.7.4
