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