44 #include "Thyra_TestingTools.hpp" 45 #include "Thyra_LinearOpTester.hpp" 46 #include "Thyra_LinearOpWithSolveTester.hpp" 48 #include "Thyra_LinearOpWithSolveFactoryHelpers.hpp" 49 #include "Thyra_DefaultSpmdVectorSpace.hpp" 50 #include "Thyra_MultiVectorStdOps.hpp" 51 #include "Epetra_SerialComm.h" 52 #include "Epetra_LocalMap.h" 53 #include "Epetra_CrsMatrix.h" 54 #include "Epetra_MultiVector.h" 55 #include "Epetra_Vector.h" 56 #include "Teuchos_dyn_cast.hpp" 57 #include "Teuchos_GlobalMPISession.hpp" 58 #include "Teuchos_CommandLineProcessor.hpp" 59 #include "Teuchos_OpaqueWrapper.hpp" 60 #include "Teuchos_TimeMonitor.hpp" 61 #include "Teuchos_StandardCatchMacros.hpp" 63 # include "Epetra_MpiComm.h" 72 void print_performance_stats(
73 const int num_time_samples
74 ,
const double raw_epetra_time
75 ,
const double thyra_wrapped_time
82 <<
"\nAverage times (out of " << num_time_samples <<
" samples):\n" 83 <<
" Raw Epetra = " << (raw_epetra_time/num_time_samples) << std::endl
84 <<
" Thyra Wrapped Epetra = " << (thyra_wrapped_time/num_time_samples) << std::endl
85 <<
"\nRelative performance of Thyra wrapped verses raw Epetra:\n" 86 <<
" ( raw epetra time / thyra wrapped time ) = ( " << raw_epetra_time <<
" / " << thyra_wrapped_time <<
" ) = " 87 << (raw_epetra_time/thyra_wrapped_time) << std::endl;
91 double sum(
const Epetra_MultiVector &ev )
93 std::vector<double> meanValues(ev.NumVectors());
94 ev.MeanValue(&meanValues[0]);
96 for(
int i = 0; i < ev.NumVectors(); ++i ) sum += meanValues[i];
97 return (ev.Map().NumGlobalElements()*sum);
115 int main(
int argc,
char* argv[] )
120 typedef double Scalar;
121 typedef Teuchos::ScalarTraits<Scalar> ST;
125 using Teuchos::dyn_cast;
126 using Teuchos::CommandLineProcessor;
129 using Teuchos::Array;
130 using Teuchos::arcpFromArray;
132 using Teuchos::rcp_static_cast;
133 using Teuchos::rcp_const_cast;
134 using Teuchos::testRelErr;
136 using Thyra::passfail;
137 using Thyra::NOTRANS;
151 Teuchos::GlobalMPISession mpiSession(&argc,&argv);
153 RCP<Teuchos::FancyOStream>
154 out = Teuchos::VerboseObjectBase::getDefaultOStream();
158 Teuchos::Time timer(
"");
164 int local_dim = 1000;
166 double max_rel_err = 1e-13;
167 double max_rel_warn = 1e-15;
169 double max_flop_rate = 1.0;
173 CommandLineProcessor clp;
174 clp.throwExceptions(
false);
175 clp.addOutputSetupOptions(
true);
176 clp.setOption(
"verbose",
"quiet", &verbose,
"Determines if any output is printed or not." );
177 clp.setOption(
"dump-all",
"no-dump", &
dumpAll,
"Determines if quantities are dumped or not." );
178 clp.setOption(
"local-dim", &local_dim,
"Number of vector elements per process." );
179 clp.setOption(
"num-mv-cols", &num_mv_cols,
"Number columns in each multi-vector (>=4)." );
180 clp.setOption(
"max-rel-err-tol", &max_rel_err,
"Maximum relative error tolerance for tests." );
181 clp.setOption(
"max-rel-warn-tol", &max_rel_warn,
"Maximum relative warning tolerance for tests." );
182 clp.setOption(
"scalar", &scalar,
"A scalar used in all computations." );
183 clp.setOption(
"max-flop-rate", &max_flop_rate,
"Approx flop rate used for loop timing." );
185 clp.setOption(
"use-mpi",
"no-use-mpi", &useMPI,
"Actually use MPI or just run independent serial programs." );
187 CommandLineProcessor::EParseCommandLineReturn parse_return = clp.parse(argc,argv);
188 if( parse_return != CommandLineProcessor::PARSE_SUCCESSFUL )
return parse_return;
190 TEUCHOS_TEST_FOR_EXCEPTION(
191 num_mv_cols < 4, std::logic_error
192 ,
"Error, num-mv-cols must be >= 4!" 203 mpiComm = MPI_COMM_WORLD;
204 MPI_Comm_size( mpiComm, &numProc );
205 MPI_Comm_rank( mpiComm, &procRank );
208 mpiComm = MPI_COMM_NULL;
217 <<
"\n*** (A) Creating two vector spaces (an Epetra-based and a non-Epetra-based)" 224 RCP<const Epetra_Comm> epetra_comm;
225 RCP<const Epetra_Map> epetra_map;
226 RCP<const Thyra::VectorSpaceBase<Scalar> > epetra_vs;
227 RCP<const Thyra::VectorSpaceBase<Scalar> > non_epetra_vs;
234 if(verbose) *out <<
"\nCreating vector space using Epetra_MpiComm ...\n";
235 epetra_comm = rcp(
new Epetra_MpiComm(mpiComm));
236 epetra_map = rcp(
new Epetra_Map(-1,local_dim,0,*epetra_comm));
239 if(verbose) *out <<
"\nCreating Thyra::DefaultSpmdVectorSpace ...\n";
241 new Thyra::DefaultSpmdVectorSpace<Scalar>(
253 if(verbose) *out <<
"\nCreating vector space using Epetra_SerialComm ...\n";
254 epetra_comm = rcp(
new Epetra_SerialComm);
255 epetra_map = rcp(
new Epetra_LocalMap(local_dim,0,*epetra_comm));
258 if(verbose) *out <<
"\nCreating Thyra::DefaultSpmdVectorSpace ...\n";
259 non_epetra_vs = Thyra::defaultSpmdVectorSpace<Scalar>(local_dim);
262 #endif // end create vector spacdes [Doxygen looks for this!] 265 const int global_dim = local_dim * numProc;
267 const int global_dim = local_dim;
272 <<
"\nscalar = " << scalar
273 <<
"\nlocal_dim = " << local_dim
274 <<
"\nglobal_dim = " << global_dim
275 <<
"\nnum_mv_cols = " << num_mv_cols
276 <<
"\nepetra_vs.dim() = " << epetra_vs->dim()
277 <<
"\nnon_epetra_vs.dim() = " << non_epetra_vs->dim()
284 RCP<Thyra::VectorBase<Scalar> >
285 ev1 = createMember(epetra_vs),
286 ev2 = createMember(epetra_vs);
287 RCP<Thyra::VectorBase<Scalar> >
288 nev1 = createMember(non_epetra_vs),
289 nev2 = createMember(non_epetra_vs);
291 RCP<Thyra::MultiVectorBase<Scalar> >
292 eV1 = createMembers(epetra_vs,num_mv_cols),
293 eV2 = createMembers(epetra_vs,num_mv_cols);
294 RCP<Thyra::MultiVectorBase<Scalar> >
295 neV1 = createMembers(non_epetra_vs,num_mv_cols),
296 neV2 = createMembers(non_epetra_vs,num_mv_cols);
301 <<
"\n*** (B) Testing Epetra and non-Epetra Thyra wrapped objects" 309 if(verbose) *out <<
"\n*** (B.1) Testing individual vector/multi-vector RTOps\n";
311 Thyra::assign( ev1.ptr(), 0.0 );
312 Thyra::assign( ev2.ptr(), scalar );
313 Thyra::assign( nev1.ptr(), 0.0 );
314 Thyra::assign( nev2.ptr(), scalar );
315 Thyra::assign( eV1.ptr(), 0.0 );
316 Thyra::assign( eV2.ptr(), scalar );
317 Thyra::assign( neV1.ptr(), 0.0 );
318 Thyra::assign( neV2.ptr(), scalar );
321 ev1_nrm = Thyra::norm_1(*ev1),
322 ev2_nrm = Thyra::norm_1(*ev2),
323 eV1_nrm = Thyra::norm_1(*eV1),
324 eV2_nrm = Thyra::norm_1(*eV2),
325 nev1_nrm = Thyra::norm_1(*nev1),
326 nev2_nrm = Thyra::norm_1(*nev2),
327 neV1_nrm = Thyra::norm_1(*neV1),
328 neV2_nrm = Thyra::norm_1(*neV2);
330 const std::string s1_n =
"fabs(scalar)*global_dim";
331 const Scalar s1 = fabs(scalar)*global_dim;
333 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;
334 if(verbose &&
dumpAll) *out <<
"\nev1 =\n" << *ev1;
335 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;
336 if(verbose &&
dumpAll) *out <<
"\nev2 =\n" << *ev2;
337 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;
338 if(verbose &&
dumpAll) *out <<
"\nnev2 =\n" << *ev1;
339 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;
340 if(verbose &&
dumpAll) *out <<
"\nnev2 =\n" << *nev2;
341 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;
342 if(verbose &&
dumpAll) *out <<
"\neV1 =\n" << *eV1;
343 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;
344 if(verbose &&
dumpAll) *out <<
"\neV2 =\n" << *eV2;
345 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;
346 if(verbose &&
dumpAll) *out <<
"\nneV1 =\n" << *neV1;
347 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;
348 if(verbose &&
dumpAll) *out <<
"\nneV2 =\n" << *neV2;
350 if(verbose) *out <<
"\n*** (B.2) Test RTOps with two or more arguments\n";
352 if(verbose) *out <<
"\nPerforming ev1 = ev2 ...\n";
354 Thyra::assign( ev1.ptr(), *ev2 );
356 if(verbose) *out <<
" time = " << timer.totalElapsedTime() <<
" sec\n";
357 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;
358 if(verbose &&
dumpAll) *out <<
"\nev1 =\n" << *ev1;
360 if(verbose) *out <<
"\nPerforming eV1 = eV2 ...\n";
362 Thyra::assign( eV1.ptr(), *eV2 );
364 if(verbose) *out <<
" time = " << timer.totalElapsedTime() <<
" sec\n";
365 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;
366 if(verbose &&
dumpAll) *out <<
"\neV1 =\n" << *eV1;
368 if(verbose) *out <<
"\nPerforming ev1 = nev2 ...\n";
370 Thyra::assign( ev1.ptr(), *nev2 );
372 if(verbose) *out <<
" time = " << timer.totalElapsedTime() <<
" sec\n";
373 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;
374 if(verbose &&
dumpAll) *out <<
"\nev1 =\n" << *ev1;
376 if(verbose) *out <<
"\nPerforming nev1 = ev2 ...\n";
378 Thyra::assign( nev1.ptr(), *ev2 );
380 if(verbose) *out <<
" time = " << timer.totalElapsedTime() <<
" sec\n";
381 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;
382 if(verbose &&
dumpAll) *out <<
"\nnev1 =\n" << *nev1;
384 if(verbose) *out <<
"\nPerforming nev1 = nev2 ...\n";
386 Thyra::assign( nev1.ptr(), *nev2 );
388 if(verbose) *out <<
" time = " << timer.totalElapsedTime() <<
" sec\n";
389 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;
390 if(verbose &&
dumpAll) *out <<
"\nnev1 =\n" << *nev1;
392 if(verbose) *out <<
"\nPerforming eV1 = neV2 ...\n";
394 Thyra::assign( eV1.ptr(), *neV2 );
396 if(verbose) *out <<
" time = " << timer.totalElapsedTime() <<
" sec\n";
397 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;
398 if(verbose &&
dumpAll) *out <<
"\neV1 =\n" << *eV1;
400 if(verbose) *out <<
"\nPerforming neV1 = eV2 ...\n";
402 Thyra::assign( neV1.ptr(), *eV2 );
404 if(verbose) *out <<
" time = " << timer.totalElapsedTime() <<
" sec\n";
405 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;
406 if(verbose &&
dumpAll) *out <<
"\nneV1 =\n" << *neV1;
408 if(verbose) *out <<
"\nPerforming neV1 = neV2 ...\n";
410 Thyra::assign( neV1.ptr(), *neV2 );
412 if(verbose) *out <<
" time = " << timer.totalElapsedTime() <<
" sec\n";
413 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;
414 if(verbose &&
dumpAll) *out <<
"\nneV1 =\n" << *neV1;
416 Thyra::LinearOpTester<Scalar> linearOpTester;
417 linearOpTester.set_all_warning_tol(max_rel_warn);
418 linearOpTester.set_all_error_tol(max_rel_err);
419 linearOpTester.dump_all(
dumpAll);
421 Thyra::LinearOpWithSolveTester<Scalar> linearOpWithSolveTester;
422 linearOpWithSolveTester.set_all_solve_tol(max_rel_err);
423 linearOpWithSolveTester.set_all_slack_error_tol(max_rel_err);
424 linearOpWithSolveTester.set_all_slack_warning_tol(max_rel_warn);
425 linearOpWithSolveTester.dump_all(
dumpAll);
428 if(verbose) *out <<
"\n*** (B.3) Test Vector linear operator interface\n";
430 if(verbose) *out <<
"\nChecking *out linear operator interface of ev1 ...\n";
431 result = linearOpTester.check(*ev1,out.ptr());
432 if(!result) success =
false;
434 if(verbose) *out <<
"\nChecking *out linear operator interface of nev1 ...\n";
435 result = linearOpTester.check(*nev1,out.ptr());
436 if(!result) success =
false;
439 if(verbose) *out <<
"\n*** (B.4) Test MultiVector linear operator interface\n";
441 if(verbose) *out <<
"\nChecking *out linear operator interface of eV1 ...\n";
442 result = linearOpTester.check(*eV1,out.ptr());
443 if(!result) success =
false;
445 if(verbose) *out <<
"\nChecking *out linear operator interface of neV1 ...\n";
446 result = linearOpTester.check(*neV1,out.ptr());
447 if(!result) success =
false;
449 const std::string s2_n =
"scalar^2*global_dim*num_mv_cols";
450 const Scalar s2 = scalar*scalar*global_dim*num_mv_cols;
452 RCP<Thyra::MultiVectorBase<Scalar> >
453 T = createMembers(eV1->domain(),num_mv_cols);
456 if(verbose) *out <<
"\n*** (B.5) Test MultiVector::apply(...)\n";
458 if(verbose) *out <<
"\nPerforming eV1'*eV2 ...\n";
460 apply( *eV1,
TRANS, *eV2, T.ptr() );
462 if(verbose) *out <<
" time = " << timer.totalElapsedTime() <<
" sec\n";
463 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;
464 if(verbose &&
dumpAll) *out <<
"\neV1'*eV2 =\n" << *T;
466 if(verbose) *out <<
"\nPerforming neV1'*eV2 ...\n";
468 apply( *neV1,
TRANS, *eV2, T.ptr() );
470 if(verbose) *out <<
" time = " << timer.totalElapsedTime() <<
" sec\n";
471 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;
472 if(verbose &&
dumpAll) *out <<
"\nneV1'*eV2 =\n" << *T;
474 if(verbose) *out <<
"\nPerforming eV1'*neV2 ...\n";
476 apply( *eV1,
TRANS, *neV2, T.ptr() );
478 if(verbose) *out <<
" time = " << timer.totalElapsedTime() <<
" sec\n";
479 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;
480 if(verbose &&
dumpAll) *out <<
"\neV1'*neV2 =\n" << *T;
482 if(verbose) *out <<
"\nPerforming neV1'*neV2 ...\n";
484 apply( *neV1,
TRANS, *neV2, T.ptr() );
486 if(verbose) *out <<
" time = " << timer.totalElapsedTime() <<
" sec\n";
487 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;
488 if(verbose &&
dumpAll) *out <<
"\nneV1'*neV2 =\n" << *T;
491 if(verbose) *out <<
"\n*** (B.6) Creating a diagonal Epetra_Operator Op\n";
493 RCP<Epetra_Operator> epetra_op;
497 RCP<Epetra_CrsMatrix>
498 epetra_mat = rcp(
new Epetra_CrsMatrix(::Copy,*epetra_map,1));
499 Scalar values[1] = { scalar };
501 const int IB = epetra_map->IndexBase(), offset = procRank*local_dim;
502 for(
int k = 0; k < local_dim; ++k ) {
503 indices[0] = offset + k + IB;
504 epetra_mat->InsertGlobalValues(
511 epetra_mat->FillComplete();
512 epetra_op = epetra_mat;
515 RCP<const Thyra::LinearOpBase<Scalar> >
516 Op = Thyra::epetraLinearOp(epetra_op);
518 if(verbose &&
dumpAll) *out <<
"\nOp=\n" << *Op;
521 if(verbose) *out <<
"\n*** (B.6b) Going through partial then full initialization of EpetraLinearOp ...\n";
525 <<
"\nChecking isFullyUninitialized(*nonconstEpetraLinearOp())==true : ";
526 RCP<Thyra::EpetraLinearOp>
527 thyraEpetraOp = Thyra::nonconstEpetraLinearOp();
528 result = isFullyUninitialized(*thyraEpetraOp);
529 if(!result) success =
false;
530 if(verbose) *out << Thyra::passfail(result) << endl;
536 <<
"\nthyraEpetraOp = partialNonconstEpetraLinearOp(...)\n";
537 RCP<Thyra::EpetraLinearOp> thyraEpetraOp =
538 Thyra::partialNonconstEpetraLinearOp(
539 epetra_vs, epetra_vs, epetra_op
543 <<
"\nChecking isPartiallyInitialized(*thyraEpetraOp)==true : ";
544 result = isPartiallyInitialized(*thyraEpetraOp);
545 if(!result) success =
false;
546 if(verbose) *out << Thyra::passfail(result) << endl;
549 <<
"\nthyraEpetraOp->setFullyInitialized(true)\n";
550 thyraEpetraOp->setFullyInitialized(
true);
553 <<
"\nChecking isFullyInitialized(*thyraEpetraOp)==true : ";
554 result = isFullyInitialized(*thyraEpetraOp);
555 if(!result) success =
false;
556 if(verbose) *out << Thyra::passfail(result) << endl;
561 if(verbose) *out <<
"\n*** (B.7) Test EpetraLinearOp linear operator interface\n";
563 if(verbose) *out <<
"\nChecking *out linear operator interface of Op ...\n";
564 result = linearOpTester.check(*Op,out.ptr());
565 if(!result) success =
false;
567 RCP<Thyra::VectorBase<Scalar> >
568 ey = createMember(epetra_vs);
569 RCP<Thyra::MultiVectorBase<Scalar> >
570 eY = createMembers(epetra_vs,num_mv_cols);
571 RCP<Thyra::VectorBase<Scalar> >
572 ney = createMember(non_epetra_vs);
573 RCP<Thyra::MultiVectorBase<Scalar> >
574 neY = createMembers(non_epetra_vs,num_mv_cols);
577 if(verbose) *out <<
"\n*** (B.8) Mix and match vector and Multi-vectors with Epetra opeator\n";
579 const std::string s3_n =
"2*scalar^2*global_dim";
580 const Scalar s3 = 2*scalar*scalar*global_dim;
582 if(verbose) *out <<
"\nPerforming ey = 2*Op*ev1 ...\n";
584 apply( *Op,
NOTRANS, *ev1, ey.ptr(), 2.0 );
586 if(verbose) *out <<
" time = " << timer.totalElapsedTime() <<
" sec\n";
587 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;
589 if(verbose) *out <<
"\nPerforming eY = 2*Op*eV1 ...\n";
591 apply( *Op,
NOTRANS, *eV1, eY.ptr(), 2.0 );
593 if(verbose) *out <<
" time = " << timer.totalElapsedTime() <<
" sec\n";
594 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;
596 if(verbose) *out <<
"\nPerforming ney = 2*Op*ev1 ...\n";
598 apply( *Op,
NOTRANS, *ev1, ney.ptr(), 2.0 );
600 if(verbose) *out <<
" time = " << timer.totalElapsedTime() <<
" sec\n";
601 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;
603 if(verbose) *out <<
"\nPerforming neY = 2*Op*eV1 ...\n";
605 apply( *Op,
NOTRANS, *eV1, neY.ptr(), 2.0 );
607 if(verbose) *out <<
" time = " << timer.totalElapsedTime() <<
" sec\n";
608 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;
610 if(verbose) *out <<
"\nPerforming ey = 2*Op*nev1 ...\n";
612 apply( *Op,
NOTRANS, *nev1, ey.ptr(), 2.0 );
614 if(verbose) *out <<
" time = " << timer.totalElapsedTime() <<
" sec\n";
615 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;
617 if(verbose) *out <<
"\nPerforming eY = 2*Op*neV1 ...\n";
619 apply( *Op,
NOTRANS, *neV1, eY.ptr(), 2.0 );
621 if(verbose) *out <<
" time = " << timer.totalElapsedTime() <<
" sec\n";
622 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;
624 if(verbose) *out <<
"\nPerforming ney = 2*Op*nev1 ...\n";
626 apply( *Op,
NOTRANS, *nev1, ney.ptr(), 2.0 );
628 if(verbose) *out <<
" time = " << timer.totalElapsedTime() <<
" sec\n";
629 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;
631 if(verbose) *out <<
"\nPerforming ney = 2*Op*nev1 through MultiVector interface ...\n";
633 apply( *Op,
NOTRANS,
static_cast<const Thyra::MultiVectorBase<Scalar>&
>(*nev1), Ptr<Thyra::MultiVectorBase<Scalar> >(ney.ptr()), 2.0 );
635 if(verbose) *out <<
" time = " << timer.totalElapsedTime() <<
" sec\n";
636 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;
638 if(verbose) *out <<
"\nPerforming neY = 2*Op*neV1 ...\n";
640 apply( *Op,
NOTRANS, *neV1, neY.ptr(), 2.0 );
642 if(verbose) *out <<
" time = " << timer.totalElapsedTime() <<
" sec\n";
643 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;
646 if(verbose) *out <<
"\n*** (B.9) Testing Multi-vector views with Epetra operator\n";
648 const Thyra::Range1D col_rng(0,1);
649 const Teuchos::Tuple<int, 2> cols = Teuchos::tuple<int>(2, 3);
651 RCP<const Thyra::MultiVectorBase<Scalar> >
652 eV1_v1 = rcp_static_cast<
const Thyra::MultiVectorBase<Scalar> >(eV1)->subView(col_rng),
653 eV1_v2 = rcp_static_cast<
const Thyra::MultiVectorBase<Scalar> >(eV1)->subView(cols);
654 RCP<const Thyra::MultiVectorBase<Scalar> >
655 neV1_v1 = rcp_static_cast<
const Thyra::MultiVectorBase<Scalar> >(neV1)->subView(col_rng),
656 neV1_v2 = rcp_static_cast<
const Thyra::MultiVectorBase<Scalar> >(neV1)->subView(cols);
658 *out <<
"\neV1_v1=\n" << *eV1_v1;
659 *out <<
"\neV1_v2=\n" << *eV1_v2;
660 *out <<
"\nneV1_v1=\n" << *neV1_v1;
661 *out <<
"\nneV1_v2=\n" << *neV1_v2;
664 if(verbose) *out <<
"\nPerforming eY_v1 = 2*Op*eV1_v1 ...\n";
666 apply( *Op,
NOTRANS, *eV1_v1, eY->subView(col_rng).ptr(), 2.0 );
668 if(verbose) *out <<
" time = " << timer.totalElapsedTime() <<
" sec\n";
669 if(verbose &&
dumpAll) *out <<
"\neV_v1=\n" << *eY->subView(col_rng);
670 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;
672 if(verbose) *out <<
"\nPerforming eY_v2 = 2*Op*eV1_v2 ...\n";
674 apply( *Op,
NOTRANS, *eV1_v2, eY->subView(cols).ptr(), 2.0 );
676 if(verbose) *out <<
" time = " << timer.totalElapsedTime() <<
" sec\n";
677 if(verbose &&
dumpAll) *out <<
"\neV_v2=\n" << *eY->subView(cols);
678 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;
680 if(verbose) *out <<
"\nPerforming neY_v1 = 2*Op*eV1_v1 ...\n";
682 apply( *Op,
NOTRANS, *eV1_v1, neY->subView(col_rng).ptr(), 2.0 );
684 if(verbose) *out <<
" time = " << timer.totalElapsedTime() <<
" sec\n";
685 if(verbose &&
dumpAll) *out <<
"\neV_v1=\n" << *eY->subView(col_rng);
686 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;
688 if(verbose) *out <<
"\nPerforming eY_v1 = 2*Op*neV1_v1 ...\n";
690 apply( *Op,
NOTRANS, *neV1_v1, eY->subView(col_rng).ptr(), 2.0 );
692 if(verbose) *out <<
" time = " << timer.totalElapsedTime() <<
" sec\n";
693 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;
695 if(verbose) *out <<
"\nPerforming neY_v2 = 2*Op*eV1_v2 ...\n";
697 apply( *Op,
NOTRANS, *eV1_v2, neY->subView(cols).ptr(), 2.0 );
699 if(verbose) *out <<
" time = " << timer.totalElapsedTime() <<
" sec\n";
700 if(verbose &&
dumpAll) *out <<
"\neV_v2=\n" << *eY->subView(cols);
701 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;
703 if(verbose) *out <<
"\nPerforming eY_v2 = 2*Op*neV1_v2 ...\n";
705 apply( *Op,
NOTRANS, *neV1_v2, eY->subView(cols).ptr(), 2.0 );
707 if(verbose) *out <<
" time = " << timer.totalElapsedTime() <<
" sec\n";
708 if(verbose &&
dumpAll) *out <<
"\neV_v2=\n" << *eY->subView(cols);
709 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;
712 if(verbose) *out <<
"\n*** (B.10) Testing Vector and MultiVector view creation functions\n";
716 const std::string s_n =
"fabs(scalar)*num_mv_cols";
717 const Scalar s = fabs(scalar)*num_mv_cols;
719 Array<Scalar> t_raw_values( num_mv_cols );
721 arcpFromArray(t_raw_values), 1 );
723 std::fill_n( t_raw_values.begin(), t_raw_values.size(), ST::zero() );
724 Thyra::assign( createMemberView(T->range(),t_raw).ptr(), scalar );
726 Scalar t_nrm = Thyra::norm_1(*t_view);
727 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;
728 if(verbose &&
dumpAll) *out <<
"\nt_view =\n" << *t_view;
741 Array<Scalar> T_raw_values( num_mv_cols * num_mv_cols );
743 arcpFromArray(T_raw_values), num_mv_cols );
745 std::fill_n( T_raw_values.begin(), T_raw_values.size(), ST::zero() );
746 Thyra::assign( createMembersView(T->range(),T_raw).ptr(), scalar );
747 Teuchos::RCP<const Thyra::MultiVectorBase<Scalar> >
749 Scalar T_nrm = Thyra::norm_1(*T_view);
750 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;
751 if(verbose &&
dumpAll) *out <<
"\nT_view =\n" << *T_view;
767 if(verbose) *out <<
"\n*** (B.11) Testing Epetra_Vector and Epetra_MultiVector wrappers\n";
771 Teuchos::RCP<const Thyra::SpmdVectorSpaceBase<Scalar> >
772 mpi_vs = Teuchos::rcp_dynamic_cast<
const Thyra::SpmdVectorSpaceBase<Scalar> >(epetra_vs,
true);
774 if(verbose) *out <<
"\nCreating temporary Epetra_Vector et1 and Epetra_MultiVector eT1 objects ...\n";
775 Teuchos::RCP<Epetra_Vector>
776 et1 = Teuchos::rcp(
new Epetra_Vector(*epetra_map));
777 Teuchos::RCP<Epetra_MultiVector>
778 eT1 = Teuchos::rcp(
new Epetra_MultiVector(*epetra_map,num_mv_cols));
780 if(verbose) *out <<
"\nCreating non-const VectorBase t1 and MultiVectorBase T1 objects from et1 and eT2 ...\n";
781 Teuchos::RCP<Thyra::VectorBase<Scalar> >
783 Teuchos::RCP<Thyra::MultiVectorBase<Scalar> >
786 if(verbose) *out <<
"\nPerforming t1 = ev1 ...\n";
787 assign( t1.ptr(), *ev1 );
789 "sum(t1)",Thyra::sum(*t1),
"sum(ev1)",sum(*ev1)
790 ,
"max_rel_err",max_rel_err,
"max_rel_warn",max_rel_warn,out.ptr())
793 if(verbose) *out <<
"\nPerforming T1 = eV1 ...\n";
794 assign( T1.ptr(), *eV1 );
796 "norm_1(T1)",Thyra::norm_1(*T1),
"norm_1(eV1)",norm_1(*eV1)
797 ,
"max_rel_err",max_rel_err,
"max_rel_warn",max_rel_warn,out.ptr())
800 if(verbose) *out <<
"\nChecking that t1 and T1 yield the same objects as et1 and eT2 ...\n";
802 Teuchos::RCP<Epetra_Vector>
804 result = et1_v.get() == et1.get();
805 if(verbose) *out <<
"\net1_v.get() = " << et1_v.get() <<
" == et1.get() = " << et1.get() <<
" : " << passfail(result) << endl;
806 if(!result) success =
false;
808 Teuchos::RCP<Epetra_MultiVector>
810 result = eT1_v.get() == eT1.get();
811 if(verbose) *out <<
"\neT1_v.get() = " << eT1_v.get() <<
" == eT1.get() = " << eT1.get() <<
" : " << passfail(result) << endl;
812 if(!result) success =
false;
814 if(verbose) *out <<
"\nCreating const VectorBase ct1 and MultiVectorBase cT1 objects from et1 and eT2 ...\n";
815 Teuchos::RCP<const Thyra::VectorBase<Scalar> >
816 ct1 =
create_Vector(Teuchos::rcp_implicit_cast<const Epetra_Vector>(et1),mpi_vs);
817 Teuchos::RCP<const Thyra::MultiVectorBase<Scalar> >
818 cT1 =
create_MultiVector(Teuchos::rcp_implicit_cast<const Epetra_MultiVector>(eT1),mpi_vs);
821 "sum(ct1)",Thyra::sum(*ct1),
"sum(ev1)",sum(*ev1)
822 ,
"max_rel_err",max_rel_err,
"max_rel_warn",max_rel_warn,out.ptr())
826 "norm_1(cT1)",Thyra::norm_1(*cT1),
"norm_1(eV1)",norm_1(*eV1)
827 ,
"max_rel_err",max_rel_err,
"max_rel_warn",max_rel_warn,out.ptr())
830 if(verbose) *out <<
"\nChecking that ct1 and cT1 yield the same objects as et1 and eT2 ...\n";
832 Teuchos::RCP<const Epetra_Vector>
834 result = cet1_v.get() == et1.get();
835 if(verbose) *out <<
"\ncet1_v.get() = " << cet1_v.get() <<
" == et1.get() = " << et1.get() <<
" : " << passfail(result) << endl;
836 if(!result) success =
false;
838 Teuchos::RCP<const Epetra_MultiVector>
840 result = ceT1_v.get() == eT1.get();
841 if(verbose) *out <<
"\nceT1_v.get() = " << ceT1_v.get() <<
" == eT1.get() = " << eT1.get() <<
" : " << passfail(result) << endl;
842 if(!result) success =
false;
844 if(verbose) *out <<
"\nCreating non-const Epetra_Vector ett1 and Epetra_MultiVector etT1 objects from clones of t1 and T2 ...\n";
845 Teuchos::RCP<Epetra_Vector>
847 Teuchos::RCP<Epetra_MultiVector>
850 if(verbose) *out <<
"\nChecking that ett1 and etT1 yield objects with the same value as et1 and eT2 ...\n";
853 "sum(ett1)",sum(*ett1),
"sum(et1)",sum(*et1)
854 ,
"max_rel_err",max_rel_err,
"max_rel_warn",max_rel_warn,out.ptr())
858 "sum(etT1)",sum(*etT1),
"sum(eT1)",sum(*eT1)
859 ,
"max_rel_err",max_rel_err,
"max_rel_warn",max_rel_warn,out.ptr())
862 if(verbose) *out <<
"\nCreating const Epetra_Vector cett1 and Epetra_MultiVector cetT1 objects from clones of t1 and T2 ...\n";
863 Teuchos::RCP<const Epetra_Vector>
864 cett1 =
get_Epetra_Vector(*epetra_map,Teuchos::rcp_implicit_cast<
const Thyra::VectorBase<Scalar> >(t1->clone_v()));
865 Teuchos::RCP<const Epetra_MultiVector>
866 cetT1 =
get_Epetra_MultiVector(*epetra_map,Teuchos::rcp_implicit_cast<
const Thyra::MultiVectorBase<Scalar> >(T1->clone_mv()));
868 if(verbose) *out <<
"\nChecking that cett1 and cetT1 yield objects with the same value as et1 and eT2 ...\n";
871 "sum(cett1)",sum(*cett1),
"sum(et1)",sum(*et1)
872 ,
"max_rel_err",max_rel_err,
"max_rel_warn",max_rel_warn,out.ptr())
876 "sum(cetT1)",sum(*cetT1),
"sum(eT1)",sum(*eT1)
877 ,
"max_rel_err",max_rel_err,
"max_rel_warn",max_rel_warn,out.ptr())
883 if(verbose) *out <<
"\n*** (B.12) Test DiagonalEpetraLinearOpWithSolveFactory \n";
887 if(verbose) *out <<
"\nUsing DiagonalEpetraLinearOpWithSolveFactory to create diagLOWS from Op ...\n";
891 Teuchos::RCP<Thyra::LinearOpWithSolveBase<Scalar> >
892 diagLOWS = Thyra::linearOpWithSolve<double>(diagLOWSFactory,Op);
894 if(verbose) *out <<
"\nTesting LinearOpBase interface of diagLOWS ...\n";
896 result = linearOpTester.check(*diagLOWS, out.ptr());
897 if(!result) success =
false;
899 if(verbose) *out <<
"\nTesting LinearOpWithSolveBase interface of diagLOWS ...\n";
901 result = linearOpWithSolveTester.check(*diagLOWS, &*out);
902 if(!result) success =
false;
910 <<
"\n*** (C) Comparing the speed of Thyra adapted Epetra objects verses raw Epetra objects" 934 double raw_epetra_time, thyra_wrapped_time;
937 if(verbose) *out <<
"\n*** (C.1) Comparing the speed of RTOp verses raw Epetra_Vector operations\n";
939 const double flop_adjust_factor_1 = 3.0;
940 const int num_time_loops_1 = int( max_flop_rate / ( flop_adjust_factor_1 * local_dim * num_mv_cols ) ) + 1;
949 *out <<
"\nPerforming eeV1 = eeV2 (using raw Epetra_MultiVector::operator=(...)) " << num_time_loops_1 <<
" times ...\n";
951 for(
int k = 0; k < num_time_loops_1; ++k ) {
955 raw_epetra_time = timer.totalElapsedTime();
956 if(verbose) *out <<
" total time = " << raw_epetra_time <<
" sec\n";
962 *out <<
"\nPerforming eV1 = eV2 (using Thyra::SpmdMultiVectorBase::applyOp(...)) " << num_time_loops_1 <<
" times ...\n";
964 for(
int k = 0; k < num_time_loops_1; ++k ) {
965 Thyra::assign( eV1.ptr(), *eV2 );
968 thyra_wrapped_time = timer.totalElapsedTime();
969 if(verbose) *out <<
" total time = " << thyra_wrapped_time <<
" sec\n";
971 print_performance_stats( num_time_loops_1, raw_epetra_time, thyra_wrapped_time, verbose, *out );
987 <<
"\n*** (C.2) Comparing Thyra::SpmdMultiVectorBase::apply() verses raw Epetra_MultiVector::Multiply()\n";
989 Teuchos::TimeMonitor::zeroOutTimers();
991 const double flop_adjust_factor_2 = 2.0;
992 const int num_time_loops_2 = int( max_flop_rate / ( flop_adjust_factor_2* local_dim * num_mv_cols * num_mv_cols ) ) + 1;
1000 Epetra_LocalMap eT_map((
int) T->range()->dim(),0,*epetra_comm);
1001 Epetra_MultiVector eT(eT_map,T->domain()->dim());
1004 *out <<
"\nPerforming eeV1'*eeV2 (using raw Epetra_MultiVector::Multiply(...)) " << num_time_loops_2 <<
" times ...\n";
1006 for(
int k = 0; k < num_time_loops_2; ++k ) {
1007 eT.Multiply(
'T',
'N', 1.0, *eeV1, *eeV2, 0.0 );
1010 raw_epetra_time = timer.totalElapsedTime();
1011 if(verbose) *out <<
" total time = " << raw_epetra_time <<
" sec\n";
1016 *out <<
"\nPerforming eV1'*eV2 (using Thyra::SpmdMultiVectorBase::apply(...)) " << num_time_loops_2 <<
" times ...\n";
1018 for(
int k = 0; k < num_time_loops_2; ++k ) {
1019 apply( *eV1,
TRANS, *eV2, T.ptr() );
1022 thyra_wrapped_time = timer.totalElapsedTime();
1023 if(verbose) *out <<
" total time = " << thyra_wrapped_time <<
" sec\n";
1025 print_performance_stats( num_time_loops_2, raw_epetra_time, thyra_wrapped_time, verbose, *out );
1036 if(verbose) *out <<
"\n*** (C.3) Comparing Thyra::EpetraLinearOp::apply() verses raw Epetra_Operator::apply()\n";
1038 Teuchos::TimeMonitor::zeroOutTimers();
1040 const double flop_adjust_factor_3 = 10.0;
1041 const int num_time_loops_3 = int( max_flop_rate / ( flop_adjust_factor_3 * local_dim * num_mv_cols ) ) + 1;
1050 *out <<
"\nPerforming eeY = eOp*eeV1 (using raw Epetra_Operator::apply(...)) " << num_time_loops_3 <<
" times ...\n";
1052 Teuchos::TimeMonitor::zeroOutTimers();
1055 epetra_op->SetUseTranspose(
false);
1056 for(
int k = 0; k < num_time_loops_3; ++k ) {
1057 epetra_op->Apply( *eeV1, *eeY );
1062 raw_epetra_time = timer.totalElapsedTime();
1063 if(verbose) *out <<
" total time = " << raw_epetra_time <<
" sec\n";
1066 Teuchos::TimeMonitor::summarize(*out <<
"\n");
1071 *out <<
"\nPerforming eY = Op*eV1 (using Thyra::EpetraLinearOp::apply(...)) " << num_time_loops_3 <<
" times ...\n";
1073 Teuchos::TimeMonitor::zeroOutTimers();
1076 for(
int k = 0; k < num_time_loops_3; ++k ) {
1077 apply( *Op,
NOTRANS, *eV1, eY.ptr() );
1081 thyra_wrapped_time = timer.totalElapsedTime();
1082 if(verbose) *out <<
" total time = " << thyra_wrapped_time <<
" sec\n";
1085 Teuchos::TimeMonitor::summarize(*out <<
"\n");
1087 print_performance_stats( num_time_loops_3, raw_epetra_time, thyra_wrapped_time, verbose, *out );
1096 TEUCHOS_STANDARD_CATCH_STATEMENTS(
true,*out,success)
1099 if(success) *out <<
"\nCongratulations! All of the tests seem to have run sucessfully!\n";
1100 else *out <<
"\nOh no! at least one of the tests did not check out!\n";
1103 return (success ? 0 : -1);
RCP< Epetra_MultiVector > get_Epetra_MultiVector(const Epetra_Map &map, const RCP< MultiVectorBase< double > > &mv)
Get a non-const Epetra_MultiVector view from a non-const MultiVectorBase object if possible...
RCP< const VectorSpaceBase< double > > create_VectorSpace(const RCP< const Epetra_Map > &epetra_map)
Create an VectorSpaceBase object given an Epetra_Map object.
RCP< MultiVectorBase< double > > create_MultiVector(const RCP< Epetra_MultiVector > &epetra_mv, const RCP< const VectorSpaceBase< double > > &range, const RCP< const VectorSpaceBase< double > > &domain=Teuchos::null)
Create a non-const MultiVectorBase object from a non-const Epetra_MultiVector object.
RCP< VectorBase< double > > create_Vector(const RCP< Epetra_Vector > &epetra_v, const RCP< const VectorSpaceBase< double > > &space)
Create a non-const VectorBase object from a non-const Epetra_Vector object.
RCP< Epetra_Vector > get_Epetra_Vector(const Epetra_Map &map, const RCP< VectorBase< double > > &v)
Get a non-const Epetra_Vector view from a non-const VectorBase object if possible.
int main(int argc, char *argv[])
Create a DefaultDiagonalLinearOpWithSolve out of a diagonal Epetra_RowMatrix object.
RCP< const Teuchos::Comm< Ordinal > > create_Comm(const RCP< const Epetra_Comm > &epetraComm)
Given an Epetra_Comm object, return an equivalent Teuchos::Comm object.