Thyra Package Browser (Single Doxygen Collection) Version of the Day
Loading...
Searching...
No Matches
test_epetra_adapters.cpp
Go to the documentation of this file.
1// @HEADER
2// ***********************************************************************
3//
4// Thyra: Interfaces and Support for Abstract Numerical Algorithms
5// Copyright (2004) Sandia Corporation
6//
7// Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
8// license for use of this work by or on behalf of the U.S. Government.
9//
10// Redistribution and use in source and binary forms, with or without
11// modification, are permitted provided that the following conditions are
12// met:
13//
14// 1. Redistributions of source code must retain the above copyright
15// notice, this list of conditions and the following disclaimer.
16//
17// 2. Redistributions in binary form must reproduce the above copyright
18// notice, this list of conditions and the following disclaimer in the
19// documentation and/or other materials provided with the distribution.
20//
21// 3. Neither the name of the Corporation nor the names of the
22// contributors may be used to endorse or promote products derived from
23// this software without specific prior written permission.
24//
25// THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
26// EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
27// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
28// PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
29// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
30// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
31// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
32// PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
33// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
34// NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
35// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
36//
37// Questions? Contact Roscoe A. Bartlett (bartlettra@ornl.gov)
38//
39// ***********************************************************************
40// @HEADER
41
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"
62#ifdef RTOp_USE_MPI
63# include "Epetra_MpiComm.h"
64#endif
65
66//
67// Some helper functions
68//
69
70namespace {
71
72void print_performance_stats(
73 const int num_time_samples
74 ,const double raw_epetra_time
75 ,const double thyra_wrapped_time
76 ,bool verbose
77 ,std::ostream &out
78 )
79{
80 if(verbose)
81 out
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;
88}
89
90inline
91double sum( const Epetra_MultiVector &ev )
92{
93 std::vector<double> meanValues(ev.NumVectors());
94 ev.MeanValue(&meanValues[0]);
95 double sum = 0;
96 for( int i = 0; i < ev.NumVectors(); ++i ) sum += meanValues[i];
97 return (ev.Map().NumGlobalElements()*sum);
98}
99
100} // namespace
101
102/* Testing program for Thyra/Epetra adpaters.
103 *
104 * This testing program shows how you can easily mix and match
105 * different implementations of vectors and multi-vectors for serial
106 * and SPMD MPI implementations. This code is worth study to show how
107 * this is done.
108 *
109 * Note that the tests performed do not prove that the Epetra adapters
110 * (or Epetra itself) perform correctly as only a few post conditions
111 * are checked. Because of the simple nature of these computations it
112 * would be possible to put in more exactly component-wise tests if
113 * that is needed in the future.
114 */
115int main( int argc, char* argv[] )
116{
117
118 using std::endl;
119
120 typedef double Scalar;
122 // typedef ST::magnitudeType ScalarMag; // unused
123 // typedef Teuchos::ScalarTraits<ScalarMag> SMT; // unused
124
125 using Teuchos::dyn_cast;
127 using Teuchos::Ptr;
128 using Teuchos::RCP;
129 using Teuchos::Array;
130 using Teuchos::arcpFromArray;
131 using Teuchos::rcp;
132 using Teuchos::rcp_static_cast;
133 using Teuchos::rcp_const_cast;
135
136 using Thyra::passfail;
137 using Thyra::NOTRANS;
138 using Thyra::TRANS;
139 using Thyra::apply;
143
144 bool verbose = true;
145 bool dumpAll = false;
146 bool success = true;
147 bool result;
148
149 int procRank = 0;
150
151 Teuchos::GlobalMPISession mpiSession(&argc,&argv);
152
153 RCP<Teuchos::FancyOStream>
155
156 try {
157
158 Teuchos::Time timer("");
159
160 //
161 // Read options from the commandline
162 //
163
164 int local_dim = 1000;
165 int num_mv_cols = 4;
166 double max_rel_err = 1e-13;
167 double max_rel_warn = 1e-15;
168 double scalar = 1.5;
169 double max_flop_rate = 1.0; // Turn off by default!
170#ifdef RTOp_USE_MPI
171 bool useMPI = true;
172#endif
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." );
184#ifdef RTOp_USE_MPI
185 clp.setOption( "use-mpi", "no-use-mpi", &useMPI, "Actually use MPI or just run independent serial programs." );
186#endif
187 CommandLineProcessor::EParseCommandLineReturn parse_return = clp.parse(argc,argv);
188 if( parse_return != CommandLineProcessor::PARSE_SUCCESSFUL ) return parse_return;
189
191 num_mv_cols < 4, std::logic_error
192 ,"Error, num-mv-cols must be >= 4!"
193 );
194
195 //
196 // Get basic MPI info
197 //
198
199#ifdef RTOp_USE_MPI
200 MPI_Comm mpiComm;
201 int numProc;
202 if(useMPI) {
203 mpiComm = MPI_COMM_WORLD;
204 MPI_Comm_size( mpiComm, &numProc );
205 MPI_Comm_rank( mpiComm, &procRank );
206 }
207 else {
208 mpiComm = MPI_COMM_NULL;
209 numProc = 1;
210 procRank = 0;
211 }
212#endif
213
214 if(verbose)
215 *out
216 << "\n***"
217 << "\n*** (A) Creating two vector spaces (an Epetra-based and a non-Epetra-based)"
218 << "\n***\n";
219
220 //
221 // Create two different vector spaces (one Epetra and one
222 // non-Epetra) that should be compatible.
223 //
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;
228#ifdef RTOp_USE_MPI
229 if(useMPI) {
230 //
231 // Create parallel vector spaces with compatible maps
232 //
233 // Epetra vector space
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));
237 epetra_vs = Thyra::create_VectorSpace(epetra_map);
238 // Non-Epetra vector space
239 if(verbose) *out << "\nCreating Thyra::DefaultSpmdVectorSpace ...\n";
240 non_epetra_vs = rcp(
242 Thyra::create_Comm(epetra_comm)
243 ,local_dim,-1
244 )
245 );
246 }
247 else {
248#endif
249 //
250 // Create serial vector spaces (i.e. VectorSpaceBase::isInCore()==true)
251 //
252 // Epetra vector space
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));
256 epetra_vs = Thyra::create_VectorSpace(epetra_map);
257 // Non-Epetra vector space
258 if(verbose) *out << "\nCreating Thyra::DefaultSpmdVectorSpace ...\n";
259 non_epetra_vs = Thyra::defaultSpmdVectorSpace<Scalar>(local_dim);
260#ifdef RTOp_USE_MPI
261 }
262#endif // end create vector spacdes [Doxygen looks for this!]
263
264#ifdef RTOp_USE_MPI
265 const int global_dim = local_dim * numProc;
266#else
267 const int global_dim = local_dim;
268#endif
269
270 if(verbose)
271 *out
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()
278 << std::endl;
279
280 //
281 // Create vectors and multi-vectors from each vector space
282 //
283
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);
290
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);
297
298 if(verbose)
299 *out
300 << "\n***"
301 << "\n*** (B) Testing Epetra and non-Epetra Thyra wrapped objects"
302 << "\n***\n";
303
304 //
305 // Check that the original epetra map can be extracted from the vector space RCP
306 // Note: get_Epetra_Map accepts also an RCP<const Epetra_Comm>, which is used to
307 // create an epetra map if the RCP<const Epetra_Map> is not found in the
308 // extra data of the inptu RCP<const VectorSpaceBase<double>>.
309 // However, in this case we KNOW that that extra data SHOULD be there, so
310 // we don't pass the comm, and let the converter throw if the data is not found.
311 //
312
313 if(verbose) *out << "\n*** (B.0) Testing idempotency of transformation: RCP<const Epetra_Map> => RCP<const VectorSpaceBase<double> => RCP<const Epetra_Map> \n";
314 RCP<const Epetra_Map> epetra_map_from_vs = get_Epetra_Map(epetra_vs);
315
316 result = epetra_map->SameAs(*epetra_map_from_vs);
317 if(!result) success = false;
318
319 //
320 // Check for compatibility of the vector and Multi-vectors
321 // w.r.t. RTOps
322 //
323
324 if(verbose) *out << "\n*** (B.1) Testing individual vector/multi-vector RTOps\n";
325
326 Thyra::assign( ev1.ptr(), 0.0 );
327 Thyra::assign( ev2.ptr(), scalar );
328 Thyra::assign( nev1.ptr(), 0.0 );
329 Thyra::assign( nev2.ptr(), scalar );
330 Thyra::assign( eV1.ptr(), 0.0 );
331 Thyra::assign( eV2.ptr(), scalar );
332 Thyra::assign( neV1.ptr(), 0.0 );
333 Thyra::assign( neV2.ptr(), scalar );
334
335 Scalar
336 ev1_nrm = Thyra::norm_1(*ev1),
337 ev2_nrm = Thyra::norm_1(*ev2),
338 eV1_nrm = Thyra::norm_1(*eV1),
339 eV2_nrm = Thyra::norm_1(*eV2),
340 nev1_nrm = Thyra::norm_1(*nev1),
341 nev2_nrm = Thyra::norm_1(*nev2),
342 neV1_nrm = Thyra::norm_1(*neV1),
343 neV2_nrm = Thyra::norm_1(*neV2);
344
345 const std::string s1_n = "fabs(scalar)*global_dim";
346 const Scalar s1 = fabs(scalar)*global_dim;
347
348 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;
349 if(verbose && dumpAll) *out << "\nev1 =\n" << *ev1;
350 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;
351 if(verbose && dumpAll) *out << "\nev2 =\n" << *ev2;
352 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;
353 if(verbose && dumpAll) *out << "\nnev2 =\n" << *ev1;
354 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;
355 if(verbose && dumpAll) *out << "\nnev2 =\n" << *nev2;
356 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;
357 if(verbose && dumpAll) *out << "\neV1 =\n" << *eV1;
358 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;
359 if(verbose && dumpAll) *out << "\neV2 =\n" << *eV2;
360 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;
361 if(verbose && dumpAll) *out << "\nneV1 =\n" << *neV1;
362 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;
363 if(verbose && dumpAll) *out << "\nneV2 =\n" << *neV2;
364
365 if(verbose) *out << "\n*** (B.2) Test RTOps with two or more arguments\n";
366
367 if(verbose) *out << "\nPerforming ev1 = ev2 ...\n";
368 timer.start(true);
369 Thyra::assign( ev1.ptr(), *ev2 );
370 timer.stop();
371 if(verbose) *out << " time = " << timer.totalElapsedTime() << " sec\n";
372 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;
373 if(verbose && dumpAll) *out << "\nev1 =\n" << *ev1;
374
375 if(verbose) *out << "\nPerforming eV1 = eV2 ...\n";
376 timer.start(true);
377 Thyra::assign( eV1.ptr(), *eV2 );
378 timer.stop();
379 if(verbose) *out << " time = " << timer.totalElapsedTime() << " sec\n";
380 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;
381 if(verbose && dumpAll) *out << "\neV1 =\n" << *eV1;
382
383 if(verbose) *out << "\nPerforming ev1 = nev2 ...\n";
384 timer.start(true);
385 Thyra::assign( ev1.ptr(), *nev2 );
386 timer.stop();
387 if(verbose) *out << " time = " << timer.totalElapsedTime() << " sec\n";
388 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;
389 if(verbose && dumpAll) *out << "\nev1 =\n" << *ev1;
390
391 if(verbose) *out << "\nPerforming nev1 = ev2 ...\n";
392 timer.start(true);
393 Thyra::assign( nev1.ptr(), *ev2 );
394 timer.stop();
395 if(verbose) *out << " time = " << timer.totalElapsedTime() << " sec\n";
396 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;
397 if(verbose && dumpAll) *out << "\nnev1 =\n" << *nev1;
398
399 if(verbose) *out << "\nPerforming nev1 = nev2 ...\n";
400 timer.start(true);
401 Thyra::assign( nev1.ptr(), *nev2 );
402 timer.stop();
403 if(verbose) *out << " time = " << timer.totalElapsedTime() << " sec\n";
404 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;
405 if(verbose && dumpAll) *out << "\nnev1 =\n" << *nev1;
406
407 if(verbose) *out << "\nPerforming eV1 = neV2 ...\n";
408 timer.start(true);
409 Thyra::assign( eV1.ptr(), *neV2 );
410 timer.stop();
411 if(verbose) *out << " time = " << timer.totalElapsedTime() << " sec\n";
412 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;
413 if(verbose && dumpAll) *out << "\neV1 =\n" << *eV1;
414
415 if(verbose) *out << "\nPerforming neV1 = eV2 ...\n";
416 timer.start(true);
417 Thyra::assign( neV1.ptr(), *eV2 );
418 timer.stop();
419 if(verbose) *out << " time = " << timer.totalElapsedTime() << " sec\n";
420 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;
421 if(verbose && dumpAll) *out << "\nneV1 =\n" << *neV1;
422
423 if(verbose) *out << "\nPerforming neV1 = neV2 ...\n";
424 timer.start(true);
425 Thyra::assign( neV1.ptr(), *neV2 );
426 timer.stop();
427 if(verbose) *out << " time = " << timer.totalElapsedTime() << " sec\n";
428 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;
429 if(verbose && dumpAll) *out << "\nneV1 =\n" << *neV1;
430
431 Thyra::LinearOpTester<Scalar> linearOpTester;
432 linearOpTester.set_all_warning_tol(max_rel_warn);
433 linearOpTester.set_all_error_tol(max_rel_err);
434 linearOpTester.dump_all(dumpAll);
435
436 Thyra::LinearOpWithSolveTester<Scalar> linearOpWithSolveTester;
437 linearOpWithSolveTester.set_all_solve_tol(max_rel_err);
438 linearOpWithSolveTester.set_all_slack_error_tol(max_rel_err);
439 linearOpWithSolveTester.set_all_slack_warning_tol(max_rel_warn);
440 linearOpWithSolveTester.dump_all(dumpAll);
441
442
443 if(verbose) *out << "\n*** (B.3) Test Vector linear operator interface\n";
444
445 if(verbose) *out << "\nChecking *out linear operator interface of ev1 ...\n";
446 result = linearOpTester.check(*ev1,out.ptr());
447 if(!result) success = false;
448
449 if(verbose) *out << "\nChecking *out linear operator interface of nev1 ...\n";
450 result = linearOpTester.check(*nev1,out.ptr());
451 if(!result) success = false;
452
453
454 if(verbose) *out << "\n*** (B.4) Test MultiVector linear operator interface\n";
455
456 if(verbose) *out << "\nChecking *out linear operator interface of eV1 ...\n";
457 result = linearOpTester.check(*eV1,out.ptr());
458 if(!result) success = false;
459
460 if(verbose) *out << "\nChecking *out linear operator interface of neV1 ...\n";
461 result = linearOpTester.check(*neV1,out.ptr());
462 if(!result) success = false;
463
464 const std::string s2_n = "scalar^2*global_dim*num_mv_cols";
465 const Scalar s2 = scalar*scalar*global_dim*num_mv_cols;
466
467 RCP<Thyra::MultiVectorBase<Scalar> >
468 T = createMembers(eV1->domain(),num_mv_cols);
469
470
471 if(verbose) *out << "\n*** (B.5) Test MultiVector::apply(...)\n";
472
473 if(verbose) *out << "\nPerforming eV1'*eV2 ...\n";
474 timer.start(true);
475 apply( *eV1, TRANS, *eV2, T.ptr() );
476 timer.stop();
477 if(verbose) *out << " time = " << timer.totalElapsedTime() << " sec\n";
478 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;
479 if(verbose && dumpAll) *out << "\neV1'*eV2 =\n" << *T;
480
481 if(verbose) *out << "\nPerforming neV1'*eV2 ...\n";
482 timer.start(true);
483 apply( *neV1, TRANS, *eV2, T.ptr() );
484 timer.stop();
485 if(verbose) *out << " time = " << timer.totalElapsedTime() << " sec\n";
486 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;
487 if(verbose && dumpAll) *out << "\nneV1'*eV2 =\n" << *T;
488
489 if(verbose) *out << "\nPerforming eV1'*neV2 ...\n";
490 timer.start(true);
491 apply( *eV1, TRANS, *neV2, T.ptr() );
492 timer.stop();
493 if(verbose) *out << " time = " << timer.totalElapsedTime() << " sec\n";
494 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;
495 if(verbose && dumpAll) *out << "\neV1'*neV2 =\n" << *T;
496
497 if(verbose) *out << "\nPerforming neV1'*neV2 ...\n";
498 timer.start(true);
499 apply( *neV1, TRANS, *neV2, T.ptr() );
500 timer.stop();
501 if(verbose) *out << " time = " << timer.totalElapsedTime() << " sec\n";
502 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;
503 if(verbose && dumpAll) *out << "\nneV1'*neV2 =\n" << *T;
504
505
506 if(verbose) *out << "\n*** (B.6) Creating a diagonal Epetra_Operator Op\n";
507
508 RCP<Epetra_Operator> epetra_op;
509
510 {
511 // Create a diagonal matrix with scalar on the diagonal
512 RCP<Epetra_CrsMatrix>
513 epetra_mat = rcp(new Epetra_CrsMatrix(::Copy,*epetra_map,1));
514 Scalar values[1] = { scalar };
515 int indices[1];
516 const int IB = epetra_map->IndexBase(), offset = procRank*local_dim;
517 for( int k = 0; k < local_dim; ++k ) {
518 indices[0] = offset + k + IB; // global column
519 epetra_mat->InsertGlobalValues(
520 offset + k + IB // GlobalRow
521 ,1 // NumEntries
522 ,values // Values
523 ,indices // Indices
524 );
525 }
526 epetra_mat->FillComplete();
527 epetra_op = epetra_mat;
528 } // end epetra_op
529
530 RCP<const Thyra::LinearOpBase<Scalar> >
531 Op = Thyra::epetraLinearOp(epetra_op);
532
533 if(verbose && dumpAll) *out << "\nOp=\n" << *Op;
534
535
536 if(verbose) *out << "\n*** (B.6b) Going through partial then full initialization of EpetraLinearOp ...\n";
537
538 {
539 if(verbose) *out
540 << "\nChecking isFullyUninitialized(*nonconstEpetraLinearOp())==true : ";
541 RCP<Thyra::EpetraLinearOp>
542 thyraEpetraOp = Thyra::nonconstEpetraLinearOp();
543 result = isFullyUninitialized(*thyraEpetraOp);
544 if(!result) success = false;
545 if(verbose) *out << Thyra::passfail(result) << endl;
546 }
547
548 {
549
550 if(verbose) *out
551 << "\nthyraEpetraOp = partialNonconstEpetraLinearOp(...)\n";
552 RCP<Thyra::EpetraLinearOp> thyraEpetraOp =
553 Thyra::partialNonconstEpetraLinearOp(
554 epetra_vs, epetra_vs, epetra_op
555 );
556
557 if(verbose) *out
558 << "\nChecking isPartiallyInitialized(*thyraEpetraOp)==true : ";
559 result = isPartiallyInitialized(*thyraEpetraOp);
560 if(!result) success = false;
561 if(verbose) *out << Thyra::passfail(result) << endl;
562
563 if(verbose) *out
564 << "\nthyraEpetraOp->setFullyInitialized(true)\n";
565 thyraEpetraOp->setFullyInitialized(true);
566
567 if(verbose) *out
568 << "\nChecking isFullyInitialized(*thyraEpetraOp)==true : ";
569 result = isFullyInitialized(*thyraEpetraOp);
570 if(!result) success = false;
571 if(verbose) *out << Thyra::passfail(result) << endl;
572
573 }
574
575
576 if(verbose) *out << "\n*** (B.7) Test EpetraLinearOp linear operator interface\n";
577
578 if(verbose) *out << "\nChecking *out linear operator interface of Op ...\n";
579 result = linearOpTester.check(*Op,out.ptr());
580 if(!result) success = false;
581
582 RCP<Thyra::VectorBase<Scalar> >
583 ey = createMember(epetra_vs);
584 RCP<Thyra::MultiVectorBase<Scalar> >
585 eY = createMembers(epetra_vs,num_mv_cols);
586 RCP<Thyra::VectorBase<Scalar> >
587 ney = createMember(non_epetra_vs);
588 RCP<Thyra::MultiVectorBase<Scalar> >
589 neY = createMembers(non_epetra_vs,num_mv_cols);
590
591
592 if(verbose) *out << "\n*** (B.8) Mix and match vector and Multi-vectors with Epetra opeator\n";
593
594 const std::string s3_n = "2*scalar^2*global_dim";
595 const Scalar s3 = 2*scalar*scalar*global_dim;
596
597 if(verbose) *out << "\nPerforming ey = 2*Op*ev1 ...\n";
598 timer.start(true);
599 apply( *Op, NOTRANS, *ev1, ey.ptr(), 2.0 );
600 timer.stop();
601 if(verbose) *out << " time = " << timer.totalElapsedTime() << " sec\n";
602 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;
603
604 if(verbose) *out << "\nPerforming eY = 2*Op*eV1 ...\n";
605 timer.start(true);
606 apply( *Op, NOTRANS, *eV1, eY.ptr(), 2.0 );
607 timer.stop();
608 if(verbose) *out << " time = " << timer.totalElapsedTime() << " sec\n";
609 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;
610
611 if(verbose) *out << "\nPerforming ney = 2*Op*ev1 ...\n";
612 timer.start(true);
613 apply( *Op, NOTRANS, *ev1, ney.ptr(), 2.0 );
614 timer.stop();
615 if(verbose) *out << " time = " << timer.totalElapsedTime() << " sec\n";
616 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;
617
618 if(verbose) *out << "\nPerforming neY = 2*Op*eV1 ...\n";
619 timer.start(true);
620 apply( *Op, NOTRANS, *eV1, neY.ptr(), 2.0 );
621 timer.stop();
622 if(verbose) *out << " time = " << timer.totalElapsedTime() << " sec\n";
623 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;
624
625 if(verbose) *out << "\nPerforming ey = 2*Op*nev1 ...\n";
626 timer.start(true);
627 apply( *Op, NOTRANS, *nev1, ey.ptr(), 2.0 );
628 timer.stop();
629 if(verbose) *out << " time = " << timer.totalElapsedTime() << " sec\n";
630 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;
631
632 if(verbose) *out << "\nPerforming eY = 2*Op*neV1 ...\n";
633 timer.start(true);
634 apply( *Op, NOTRANS, *neV1, eY.ptr(), 2.0 );
635 timer.stop();
636 if(verbose) *out << " time = " << timer.totalElapsedTime() << " sec\n";
637 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;
638
639 if(verbose) *out << "\nPerforming ney = 2*Op*nev1 ...\n";
640 timer.start(true);
641 apply( *Op, NOTRANS, *nev1, ney.ptr(), 2.0 );
642 timer.stop();
643 if(verbose) *out << " time = " << timer.totalElapsedTime() << " sec\n";
644 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;
645
646 if(verbose) *out << "\nPerforming ney = 2*Op*nev1 through MultiVector interface ...\n";
647 timer.start(true);
648 apply( *Op, NOTRANS, static_cast<const Thyra::MultiVectorBase<Scalar>&>(*nev1), Ptr<Thyra::MultiVectorBase<Scalar> >(ney.ptr()), 2.0 );
649 timer.stop();
650 if(verbose) *out << " time = " << timer.totalElapsedTime() << " sec\n";
651 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;
652
653 if(verbose) *out << "\nPerforming neY = 2*Op*neV1 ...\n";
654 timer.start(true);
655 apply( *Op, NOTRANS, *neV1, neY.ptr(), 2.0 );
656 timer.stop();
657 if(verbose) *out << " time = " << timer.totalElapsedTime() << " sec\n";
658 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;
659
660
661 if(verbose) *out << "\n*** (B.9) Testing Multi-vector views with Epetra operator\n";
662
663 const Thyra::Range1D col_rng(0,1);
665
666 RCP<const Thyra::MultiVectorBase<Scalar> >
667 eV1_v1 = rcp_static_cast<const Thyra::MultiVectorBase<Scalar> >(eV1)->subView(col_rng),
668 eV1_v2 = rcp_static_cast<const Thyra::MultiVectorBase<Scalar> >(eV1)->subView(cols);
669 RCP<const Thyra::MultiVectorBase<Scalar> >
670 neV1_v1 = rcp_static_cast<const Thyra::MultiVectorBase<Scalar> >(neV1)->subView(col_rng),
671 neV1_v2 = rcp_static_cast<const Thyra::MultiVectorBase<Scalar> >(neV1)->subView(cols);
672 if(verbose && dumpAll) {
673 *out << "\neV1_v1=\n" << *eV1_v1;
674 *out << "\neV1_v2=\n" << *eV1_v2;
675 *out << "\nneV1_v1=\n" << *neV1_v1;
676 *out << "\nneV1_v2=\n" << *neV1_v2;
677 }
678
679 if(verbose) *out << "\nPerforming eY_v1 = 2*Op*eV1_v1 ...\n";
680 timer.start(true);
681 apply( *Op, NOTRANS, *eV1_v1, eY->subView(col_rng).ptr(), 2.0 );
682 timer.stop();
683 if(verbose) *out << " time = " << timer.totalElapsedTime() << " sec\n";
684 if(verbose && dumpAll) *out << "\neV_v1=\n" << *eY->subView(col_rng);
685 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;
686
687 if(verbose) *out << "\nPerforming eY_v2 = 2*Op*eV1_v2 ...\n";
688 timer.start(true);
689 apply( *Op, NOTRANS, *eV1_v2, eY->subView(cols).ptr(), 2.0 );
690 timer.stop();
691 if(verbose) *out << " time = " << timer.totalElapsedTime() << " sec\n";
692 if(verbose && dumpAll) *out << "\neV_v2=\n" << *eY->subView(cols);
693 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;
694
695 if(verbose) *out << "\nPerforming neY_v1 = 2*Op*eV1_v1 ...\n";
696 timer.start(true);
697 apply( *Op, NOTRANS, *eV1_v1, neY->subView(col_rng).ptr(), 2.0 );
698 timer.stop();
699 if(verbose) *out << " time = " << timer.totalElapsedTime() << " sec\n";
700 if(verbose && dumpAll) *out << "\neV_v1=\n" << *eY->subView(col_rng);
701 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;
702
703 if(verbose) *out << "\nPerforming eY_v1 = 2*Op*neV1_v1 ...\n";
704 timer.start(true);
705 apply( *Op, NOTRANS, *neV1_v1, eY->subView(col_rng).ptr(), 2.0 );
706 timer.stop();
707 if(verbose) *out << " time = " << timer.totalElapsedTime() << " sec\n";
708 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;
709
710 if(verbose) *out << "\nPerforming neY_v2 = 2*Op*eV1_v2 ...\n";
711 timer.start(true);
712 apply( *Op, NOTRANS, *eV1_v2, neY->subView(cols).ptr(), 2.0 );
713 timer.stop();
714 if(verbose) *out << " time = " << timer.totalElapsedTime() << " sec\n";
715 if(verbose && dumpAll) *out << "\neV_v2=\n" << *eY->subView(cols);
716 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;
717
718 if(verbose) *out << "\nPerforming eY_v2 = 2*Op*neV1_v2 ...\n";
719 timer.start(true);
720 apply( *Op, NOTRANS, *neV1_v2, eY->subView(cols).ptr(), 2.0 );
721 timer.stop();
722 if(verbose) *out << " time = " << timer.totalElapsedTime() << " sec\n";
723 if(verbose && dumpAll) *out << "\neV_v2=\n" << *eY->subView(cols);
724 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;
725
726
727 if(verbose) *out << "\n*** (B.10) Testing Vector and MultiVector view creation functions\n";
728
729 {
730
731 const std::string s_n = "fabs(scalar)*num_mv_cols";
732 const Scalar s = fabs(scalar)*num_mv_cols;
733
734 Array<Scalar> t_raw_values( num_mv_cols );
735 RTOpPack::SubVectorView<Scalar> t_raw( 0, num_mv_cols,
736 arcpFromArray(t_raw_values), 1 );
737
738 std::fill_n( t_raw_values.begin(), t_raw_values.size(), ST::zero() );
739 Thyra::assign( createMemberView(T->range(),t_raw).ptr(), scalar );
740 Teuchos::RCP<const Thyra::VectorBase<Scalar> > t_view = createMemberView(T->range(),static_cast<RTOpPack::ConstSubVectorView<Scalar>&>(t_raw));
741 Scalar t_nrm = Thyra::norm_1(*t_view);
742 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;
743 if(verbose && dumpAll) *out << "\nt_view =\n" << *t_view;
744
745#if 0
746 std::fill_n( t_raw_values.begin(), t_raw_values.size(), ST::zero() );
747 Thyra::assign( T->range().ptr()->Thyra::VectorSpaceBase<Scalar>::createMemberView(t_raw), scalar );
748 t_view = T->range()->Thyra::VectorSpaceBase<Scalar>::createMemberView(static_cast<RTOpPack::ConstSubVectorView<Scalar>&>(t_raw));
749 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;
752#endif
753
754 Array<Scalar> T_raw_values( num_mv_cols * num_mv_cols );
755 RTOpPack::SubMultiVectorView<Scalar> T_raw( 0, num_mv_cols, 0, num_mv_cols,
756 arcpFromArray(T_raw_values), num_mv_cols );
757
758 std::fill_n( T_raw_values.begin(), T_raw_values.size(), ST::zero() );
759 Thyra::assign( createMembersView(T->range(),T_raw).ptr(), scalar );
761 T_view = createMembersView(T->range(),static_cast<RTOpPack::ConstSubMultiVectorView<Scalar>&>(T_raw));
762 Scalar T_nrm = Thyra::norm_1(*T_view);
763 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;
764 if(verbose && dumpAll) *out << "\nT_view =\n" << *T_view;
765
766#if 0
767 std::fill_n( T_raw_values.begin(), T_raw_values.size(), ST::zero() );
768 Thyra::assign( T->range().ptr()->Thyra::VectorSpaceBase<Scalar>::createMembersView(T_raw), scalar );
769 T_view = T->range()->Thyra::VectorSpaceBase<Scalar>::createMembersView(static_cast<RTOpPack::ConstSubMultiVectorView<Scalar>&>(T_raw));
770 T_nrm = Thyra::norm_1(*T_view);
771 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;
772 if(verbose && dumpAll) *out << "\nT_view =\n" << *T_view;
773#endif
774
775 }
776
777
778 if(verbose) *out << "\n*** (B.11) Testing Epetra_Vector and Epetra_MultiVector wrappers\n";
779
780 {
781
784
785 if(verbose) *out << "\nCreating temporary Epetra_Vector et1 and Epetra_MultiVector eT1 objects ...\n";
787 et1 = Teuchos::rcp(new Epetra_Vector(*epetra_map));
789 eT1 = Teuchos::rcp(new Epetra_MultiVector(*epetra_map,num_mv_cols));
790
791 if(verbose) *out << "\nCreating non-const VectorBase t1 and MultiVectorBase T1 objects from et1 and eT2 ...\n";
793 t1 = create_Vector(et1,mpi_vs);
795 T1 = create_MultiVector(eT1,mpi_vs);
796
797 if(verbose) *out << "\nPerforming t1 = ev1 ...\n";
798 assign( t1.ptr(), *ev1 );
799 if(!testRelErr(
800 "sum(t1)",Thyra::sum(*t1),"sum(ev1)",sum(*ev1)
801 ,"max_rel_err",max_rel_err,"max_rel_warn",max_rel_warn,out.ptr())
802 ) success=false;
803
804 if(verbose) *out << "\nPerforming T1 = eV1 ...\n";
805 assign( T1.ptr(), *eV1 );
806 if(!testRelErr(
807 "norm_1(T1)",Thyra::norm_1(*T1),"norm_1(eV1)",norm_1(*eV1)
808 ,"max_rel_err",max_rel_err,"max_rel_warn",max_rel_warn,out.ptr())
809 ) success=false;
810
811 if(verbose) *out << "\nChecking that t1 and T1 yield the same objects as et1 and eT2 ...\n";
812
814 et1_v = get_Epetra_Vector(*epetra_map,t1);
815 result = et1_v.get() == et1.get();
816 if(verbose) *out << "\net1_v.get() = " << et1_v.get() << " == et1.get() = " << et1.get() << " : " << passfail(result) << endl;
817 if(!result) success = false;
818
819 // Like the above, but looks for RCP<Epetra_Vector> in t1's extra data
820 et1_v = get_Epetra_Vector(t1);
821 result = et1_v.get() == et1.get();
822 if(verbose) *out << "\net1_v.get() = " << et1_v.get() << " == et1.get() = " << et1.get() << " : " << passfail(result) << endl;
823 if(!result) success = false;
824
826 eT1_v = get_Epetra_MultiVector(*epetra_map,T1);
827 result = eT1_v.get() == eT1.get();
828 if(verbose) *out << "\neT1_v.get() = " << eT1_v.get() << " == eT1.get() = " << eT1.get() << " : " << passfail(result) << endl;
829 if(!result) success = false;
830
831 // Like the above, but looks for RCP<Epetra_MultiVector> in T1's extra data
832 eT1_v = get_Epetra_MultiVector(T1);
833 result = eT1_v.get() == eT1.get();
834 if(verbose) *out << "\neT1_v.get() = " << eT1_v.get() << " == eT1.get() = " << eT1.get() << " : " << passfail(result) << endl;
835 if(!result) success = false;
836
837 if(verbose) *out << "\nCreating const VectorBase ct1 and MultiVectorBase cT1 objects from et1 and eT2 ...\n";
839 ct1 = create_Vector(Teuchos::rcp_implicit_cast<const Epetra_Vector>(et1),mpi_vs);
841 cT1 = create_MultiVector(Teuchos::rcp_implicit_cast<const Epetra_MultiVector>(eT1),mpi_vs);
842
843 if(!testRelErr(
844 "sum(ct1)",Thyra::sum(*ct1),"sum(ev1)",sum(*ev1)
845 ,"max_rel_err",max_rel_err,"max_rel_warn",max_rel_warn,out.ptr())
846 ) success=false;
847
848 if(!testRelErr(
849 "norm_1(cT1)",Thyra::norm_1(*cT1),"norm_1(eV1)",norm_1(*eV1)
850 ,"max_rel_err",max_rel_err,"max_rel_warn",max_rel_warn,out.ptr())
851 ) success=false;
852
853 if(verbose) *out << "\nChecking that ct1 and cT1 yield the same objects as et1 and eT2 ...\n";
854
856 cet1_v = get_Epetra_Vector(*epetra_map,ct1);
857 result = cet1_v.get() == et1.get();
858 if(verbose) *out << "\ncet1_v.get() = " << cet1_v.get() << " == et1.get() = " << et1.get() << " : " << passfail(result) << endl;
859 if(!result) success = false;
860
861 // Like the above, but looks for RCP<const Epetra_Vector> in ct1's extra data
862 cet1_v = get_Epetra_Vector(ct1);
863 result = cet1_v.get() == et1.get();
864 if(verbose) *out << "\ncet1_v.get() = " << cet1_v.get() << " == et1.get() = " << et1.get() << " : " << passfail(result) << endl;
865 if(!result) success = false;
866
868 ceT1_v = get_Epetra_MultiVector(*epetra_map,cT1);
869 result = ceT1_v.get() == eT1.get();
870 if(verbose) *out << "\nceT1_v.get() = " << ceT1_v.get() << " == eT1.get() = " << eT1.get() << " : " << passfail(result) << endl;
871 if(!result) success = false;
872
873 // Like the above, but looks for RCP<const Epetra_MultiVector> in ct1's extra data
874 ceT1_v = get_Epetra_MultiVector(cT1);
875 result = ceT1_v.get() == eT1.get();
876 if(verbose) *out << "\nceT1_v.get() = " << ceT1_v.get() << " == eT1.get() = " << eT1.get() << " : " << passfail(result) << endl;
877 if(!result) success = false;
878
879 if(verbose) *out << "\nCreating non-const Epetra_Vector ett1 and Epetra_MultiVector etT1 objects from clones of t1 and T2 ...\n";
881 ett1 = get_Epetra_Vector(*epetra_map,t1->clone_v());
883 etT1 = get_Epetra_MultiVector(*epetra_map,T1->clone_mv());
884
885 if(verbose) *out << "\nChecking that ett1 and etT1 yield objects with the same value as et1 and eT2 ...\n";
886
887 if(!testRelErr(
888 "sum(ett1)",sum(*ett1),"sum(et1)",sum(*et1)
889 ,"max_rel_err",max_rel_err,"max_rel_warn",max_rel_warn,out.ptr())
890 ) success=false;
891
892 if(!testRelErr(
893 "sum(etT1)",sum(*etT1),"sum(eT1)",sum(*eT1)
894 ,"max_rel_err",max_rel_err,"max_rel_warn",max_rel_warn,out.ptr())
895 ) success=false;
896
897 if(verbose) *out << "\nCreating const Epetra_Vector cett1 and Epetra_MultiVector cetT1 objects from clones of t1 and T2 ...\n";
899 cett1 = get_Epetra_Vector(*epetra_map,Teuchos::rcp_implicit_cast<const Thyra::VectorBase<Scalar> >(t1->clone_v()));
901 cetT1 = get_Epetra_MultiVector(*epetra_map,Teuchos::rcp_implicit_cast<const Thyra::MultiVectorBase<Scalar> >(T1->clone_mv()));
902
903 if(verbose) *out << "\nChecking that cett1 and cetT1 yield objects with the same value as et1 and eT2 ...\n";
904
905 if(!testRelErr(
906 "sum(cett1)",sum(*cett1),"sum(et1)",sum(*et1)
907 ,"max_rel_err",max_rel_err,"max_rel_warn",max_rel_warn,out.ptr())
908 ) success=false;
909
910 if(!testRelErr(
911 "sum(cetT1)",sum(*cetT1),"sum(eT1)",sum(*eT1)
912 ,"max_rel_err",max_rel_err,"max_rel_warn",max_rel_warn,out.ptr())
913 ) success=false;
914
915 }
916
917
918 if(verbose) *out << "\n*** (B.12) Test DiagonalEpetraLinearOpWithSolveFactory \n";
919
920 {
921
922 if(verbose) *out << "\nUsing DiagonalEpetraLinearOpWithSolveFactory to create diagLOWS from Op ...\n";
923
925
927 diagLOWS = Thyra::linearOpWithSolve<double>(diagLOWSFactory,Op);
928
929 if(verbose) *out << "\nTesting LinearOpBase interface of diagLOWS ...\n";
930
931 result = linearOpTester.check(*diagLOWS, out.ptr());
932 if(!result) success = false;
933
934 if(verbose) *out << "\nTesting LinearOpWithSolveBase interface of diagLOWS ...\n";
935
936 result = linearOpWithSolveTester.check(*diagLOWS, &*out);
937 if(!result) success = false;
938
939 }
940
941
942 if(verbose)
943 *out
944 << "\n***"
945 << "\n*** (C) Comparing the speed of Thyra adapted Epetra objects verses raw Epetra objects"
946 << "\n***\n";
947
948 //
949 // Setup the number of timing loops to get good timings
950 //
951 // Here we try to shoot for timing ab*out 1 second's worth of
952 // computations and adjust the number of evaluation loops
953 // accordingly. Let X be the approximate number of flops per
954 // loop (or per evaluation). We then compute the number of
955 // loops as:
956 //
957 // 1.0 sec | num CPU flops | 1 loop |
958 // --------|----------------|-----------|
959 // | sec | X flops |
960 //
961 // This just comes *out to be:
962 //
963 // num_time_loops_X = max_flop_rate / (X flops per loop)
964 //
965 // In this computation we ignore extra overhead that will be
966 // an issue when local_dim is small.
967 //
968
969 double raw_epetra_time, thyra_wrapped_time;
970
971
972 if(verbose) *out << "\n*** (C.1) Comparing the speed of RTOp verses raw Epetra_Vector operations\n";
973
974 const double flop_adjust_factor_1 = 3.0;
975 const int num_time_loops_1 = int( max_flop_rate / ( flop_adjust_factor_1 * local_dim * num_mv_cols ) ) + 1;
976
977 {
978
979 // Get references to Epetra_MultiVector objects in eV1 and eV2
980 const RCP<Epetra_MultiVector> eeV1 = get_Epetra_MultiVector(*epetra_map,eV1);
981 const RCP<const Epetra_MultiVector> eeV2 = get_Epetra_MultiVector(*epetra_map,eV2);
982
983 if(verbose)
984 *out << "\nPerforming eeV1 = eeV2 (using raw Epetra_MultiVector::operator=(...)) " << num_time_loops_1 << " times ...\n";
985 timer.start(true);
986 for(int k = 0; k < num_time_loops_1; ++k ) {
987 *eeV1 = *eeV2;
988 }
989 timer.stop();
990 raw_epetra_time = timer.totalElapsedTime();
991 if(verbose) *out << " total time = " << raw_epetra_time << " sec\n";
992
993 // When this block ends and eeV1 goes *out of scope then eV1 is guaranteed to be updated!
994 }
995
996 if(verbose)
997 *out << "\nPerforming eV1 = eV2 (using Thyra::SpmdMultiVectorBase::applyOp(...)) " << num_time_loops_1 << " times ...\n";
998 timer.start(true);
999 for(int k = 0; k < num_time_loops_1; ++k ) {
1000 Thyra::assign( eV1.ptr(), *eV2 );
1001 }
1002 timer.stop();
1003 thyra_wrapped_time = timer.totalElapsedTime();
1004 if(verbose) *out << " total time = " << thyra_wrapped_time << " sec\n";
1005
1006 print_performance_stats( num_time_loops_1, raw_epetra_time, thyra_wrapped_time, verbose, *out );
1007
1008 // RAB: 2004/01/05: Note, the above relative performance is likely
1009 // to be the worst of all of the others since RTOp operators are
1010 // applied seperately column by column but the relative
1011 // performance should go to about 1.0 when local_dim is
1012 // sufficiently large! However, because
1013 // Epetra_MultiVector::Thyra::Assign(...) is implemented using double
1014 // pointer indexing, the RTOp implementation used with the Thyra
1015 // adapters is actually faster in some cases. However, the extra
1016 // overhead of RTOp is much worse for very very small (order 10)
1017 // sizes.
1018
1019
1020 if(verbose)
1021 *out
1022 << "\n*** (C.2) Comparing Thyra::SpmdMultiVectorBase::apply() verses raw Epetra_MultiVector::Multiply()\n";
1023
1025
1026 const double flop_adjust_factor_2 = 2.0;
1027 const int num_time_loops_2 = int( max_flop_rate / ( flop_adjust_factor_2* local_dim * num_mv_cols * num_mv_cols ) ) + 1;
1028
1029 {
1030
1031 // Get constant references to Epetra_MultiVector objects in eV1 and eV2
1032 const RCP<const Epetra_MultiVector> eeV1 = get_Epetra_MultiVector(*epetra_map,eV1);
1033 const RCP<const Epetra_MultiVector> eeV2 = get_Epetra_MultiVector(*epetra_map,eV2);
1034
1035 Epetra_LocalMap eT_map((int) T->range()->dim(),0,*epetra_comm);
1036 Epetra_MultiVector eT(eT_map,T->domain()->dim());
1037
1038 if(verbose)
1039 *out << "\nPerforming eeV1'*eeV2 (using raw Epetra_MultiVector::Multiply(...)) " << num_time_loops_2 << " times ...\n";
1040 timer.start(true);
1041 for(int k = 0; k < num_time_loops_2; ++k ) {
1042 eT.Multiply( 'T', 'N', 1.0, *eeV1, *eeV2, 0.0 );
1043 }
1044 timer.stop();
1045 raw_epetra_time = timer.totalElapsedTime();
1046 if(verbose) *out << " total time = " << raw_epetra_time << " sec\n";
1047
1048 }
1049
1050 if(verbose)
1051 *out << "\nPerforming eV1'*eV2 (using Thyra::SpmdMultiVectorBase::apply(...)) " << num_time_loops_2 << " times ...\n";
1052 timer.start(true);
1053 for(int k = 0; k < num_time_loops_2; ++k ) {
1054 apply( *eV1, TRANS, *eV2, T.ptr() );
1055 }
1056 timer.stop();
1057 thyra_wrapped_time = timer.totalElapsedTime();
1058 if(verbose) *out << " total time = " << thyra_wrapped_time << " sec\n";
1059
1060 print_performance_stats( num_time_loops_2, raw_epetra_time, thyra_wrapped_time, verbose, *out );
1061
1062 // RAB: 2004/01/05: Note, even though the Thyra adapter does
1063 // not actually call Epetra_MultiVector::Multiply(...), the
1064 // implementation in Thyra::SpmdMultiVectorBase::apply(...)
1065 // performs almost exactly the same flops and calls dgemm(...)
1066 // as well. Herefore, except for some small overhead, the raw
1067 // Epetra and the Thyra wrapped computations should give
1068 // almost identical times in almost all cases.
1069
1070
1071 if(verbose) *out << "\n*** (C.3) Comparing Thyra::EpetraLinearOp::apply() verses raw Epetra_Operator::apply()\n";
1072
1074
1075 const double flop_adjust_factor_3 = 10.0; // lots of indirect addressing
1076 const int num_time_loops_3 = int( max_flop_rate / ( flop_adjust_factor_3 * local_dim * num_mv_cols ) ) + 1;
1077
1078 {
1079
1080 // Get constant references to Epetra_MultiVector objects in eV1 and eV2
1081 const RCP<const Epetra_MultiVector> eeV1 = get_Epetra_MultiVector(*epetra_map,eV1);
1082 const RCP<Epetra_MultiVector> eeY = get_Epetra_MultiVector(*epetra_map,eY);
1083
1084 if(verbose)
1085 *out << "\nPerforming eeY = eOp*eeV1 (using raw Epetra_Operator::apply(...)) " << num_time_loops_3 << " times ...\n";
1086
1088
1089 timer.start(true);
1090 epetra_op->SetUseTranspose(false);
1091 for(int k = 0; k < num_time_loops_3; ++k ) {
1092 epetra_op->Apply( *eeV1, *eeY );
1093 //eeY->Scale(2.0);
1094 }
1095 timer.stop();
1096
1097 raw_epetra_time = timer.totalElapsedTime();
1098 if(verbose) *out << " total time = " << raw_epetra_time << " sec\n";
1099
1100 if(verbose)
1101 Teuchos::TimeMonitor::summarize(*out << "\n");
1102
1103 }
1104
1105 if(verbose)
1106 *out << "\nPerforming eY = Op*eV1 (using Thyra::EpetraLinearOp::apply(...)) " << num_time_loops_3 << " times ...\n";
1107
1109
1110 timer.start(true);
1111 for(int k = 0; k < num_time_loops_3; ++k ) {
1112 apply( *Op, NOTRANS, *eV1, eY.ptr() );
1113 }
1114 timer.stop();
1115
1116 thyra_wrapped_time = timer.totalElapsedTime();
1117 if(verbose) *out << " total time = " << thyra_wrapped_time << " sec\n";
1118
1119 if(verbose)
1120 Teuchos::TimeMonitor::summarize(*out << "\n");
1121
1122 print_performance_stats( num_time_loops_3, raw_epetra_time, thyra_wrapped_time, verbose, *out );
1123
1124 // RAB: 2004/01/05: Note, the above Epetra adapter is a true
1125 // adapter and simply calls Epetra_Operator::apply(...) so, except
1126 // for some small overhead, the raw Epetra and the Thyra wrapped
1127 // computations should give ab*out exactly the same runtime for
1128 // almost all cases.
1129
1130 } // end try
1131 TEUCHOS_STANDARD_CATCH_STATEMENTS(true,*out,success)
1132
1133 if(verbose) {
1134 if(success) *out << "\nCongratulations! All of the tests seem to have run sucessfully!\n";
1135 else *out << "\nOh no! at least one of the tests did not check out!\n";
1136 }
1137
1138 return (success ? 0 : -1);
1139
1140} // end main()
int main(int argc, char *argv[])
#define TEUCHOS_STANDARD_CATCH_STATEMENTS(VERBOSE, ERR_STREAM, SUCCESS_FLAG)
Ptr< T > ptr() const
T * get() const
static void zeroOutTimers()
static void summarize(Ptr< const Comm< int > > comm, std::ostream &out=std::cout, const bool alwaysWriteLocal=false, const bool writeGlobalStats=true, const bool writeZeroTimers=true, const ECounterSetOp setOp=Intersection, const std::string &filter="", const bool ignoreZeroTimers=false)
double totalElapsedTime(bool readCurrentTime=false) const
void start(bool reset=false)
double stop()
static RCP< FancyOStream > getDefaultOStream()
Create a DefaultDiagonalLinearOpWithSolve out of a diagonal Epetra_RowMatrix object.
RCP< VectorBase< double > > create_Vector(const RCP< Epetra_Vector > &epetra_v, const RCP< const VectorSpaceBase< double > > &space=Teuchos::null)
Create a non-const VectorBase object from a non-const Epetra_Vector object.
RCP< MultiVectorBase< double > > create_MultiVector(const RCP< Epetra_MultiVector > &epetra_mv, const RCP< const VectorSpaceBase< double > > &range=Teuchos::null, const RCP< const VectorSpaceBase< double > > &domain=Teuchos::null)
Create a non-const MultiVectorBase object from a non-const Epetra_MultiVector 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.
RCP< const VectorSpaceBase< double > > create_VectorSpace(const RCP< const Epetra_Map > &epetra_map)
Create an VectorSpaceBase object given an Epetra_Map object.
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
bool testRelErr(const std::string &v1_name, const T1 &v1, const std::string &v2_name, const T2 &v2, const std::string &maxRelErr_error_name, const typename TestRelErr< T1, T2 >::magnitudeType &maxRelErr_error, const std::string &maxRelErr_warning_name, const typename TestRelErr< T1, T2 >::magnitudeType &maxRelErr_warning, const Ptr< std::ostream > &out)
T_To & dyn_cast(T_From &from)