Thyra Package Browser (Single Doxygen Collection)  Version of the Day
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 
43 #include "Thyra_EpetraLinearOp.hpp"
44 #include "Thyra_TestingTools.hpp"
45 #include "Thyra_LinearOpTester.hpp"
46 #include "Thyra_LinearOpWithSolveTester.hpp"
48 #include "Thyra_LinearOpWithSolveFactoryHelpers.hpp"
49 #include "Thyra_DefaultSpmdVectorSpace.hpp"
50 #include "Thyra_MultiVectorStdOps.hpp"
51 #include "Epetra_SerialComm.h"
52 #include "Epetra_LocalMap.h"
53 #include "Epetra_CrsMatrix.h"
54 #include "Epetra_MultiVector.h"
55 #include "Epetra_Vector.h"
56 #include "Teuchos_dyn_cast.hpp"
57 #include "Teuchos_GlobalMPISession.hpp"
58 #include "Teuchos_CommandLineProcessor.hpp"
59 #include "Teuchos_OpaqueWrapper.hpp"
60 #include "Teuchos_TimeMonitor.hpp"
61 #include "Teuchos_StandardCatchMacros.hpp"
62 #ifdef RTOp_USE_MPI
63 # include "Epetra_MpiComm.h"
64 #endif
65 
66 //
67 // Some helper functions
68 //
69 
70 namespace {
71 
72 void 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 
90 inline
91 double 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  */
115 int main( int argc, char* argv[] )
116 {
117 
118  using std::endl;
119 
120  typedef double Scalar;
121  typedef Teuchos::ScalarTraits<Scalar> ST;
122  // typedef ST::magnitudeType ScalarMag; // unused
123  // typedef Teuchos::ScalarTraits<ScalarMag> SMT; // unused
124 
125  using Teuchos::dyn_cast;
126  using Teuchos::CommandLineProcessor;
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;
134  using Teuchos::testRelErr;
135 
136  using Thyra::passfail;
137  using Thyra::NOTRANS;
138  using Thyra::TRANS;
139  using Thyra::apply;
141  using Thyra::create_Vector;
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>
154  out = Teuchos::VerboseObjectBase::getDefaultOStream();
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 
190  TEUCHOS_TEST_FOR_EXCEPTION(
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(
241  new Thyra::DefaultSpmdVectorSpace<Scalar>(
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 for compatibility of the vector and Multi-vectors
306  // w.r.t. RTOps
307  //
308 
309  if(verbose) *out << "\n*** (B.1) Testing individual vector/multi-vector RTOps\n";
310 
311  Thyra::assign( ev1.ptr(), 0.0 );
312  Thyra::assign( ev2.ptr(), scalar );
313  Thyra::assign( nev1.ptr(), 0.0 );
314  Thyra::assign( nev2.ptr(), scalar );
315  Thyra::assign( eV1.ptr(), 0.0 );
316  Thyra::assign( eV2.ptr(), scalar );
317  Thyra::assign( neV1.ptr(), 0.0 );
318  Thyra::assign( neV2.ptr(), scalar );
319 
320  Scalar
321  ev1_nrm = Thyra::norm_1(*ev1),
322  ev2_nrm = Thyra::norm_1(*ev2),
323  eV1_nrm = Thyra::norm_1(*eV1),
324  eV2_nrm = Thyra::norm_1(*eV2),
325  nev1_nrm = Thyra::norm_1(*nev1),
326  nev2_nrm = Thyra::norm_1(*nev2),
327  neV1_nrm = Thyra::norm_1(*neV1),
328  neV2_nrm = Thyra::norm_1(*neV2);
329 
330  const std::string s1_n = "fabs(scalar)*global_dim";
331  const Scalar s1 = fabs(scalar)*global_dim;
332 
333  if(!testRelErr("Thyra::norm_1(ev1)",ev1_nrm,"0",Scalar(0),"max_rel_err",max_rel_err,"max_rel_warn",max_rel_warn,out.ptr())) success=false;
334  if(verbose && dumpAll) *out << "\nev1 =\n" << *ev1;
335  if(!testRelErr("Thyra::norm_1(ev2)",ev2_nrm,s1_n,s1,"max_rel_err",max_rel_err,"max_rel_warn",max_rel_warn,out.ptr())) success=false;
336  if(verbose && dumpAll) *out << "\nev2 =\n" << *ev2;
337  if(!testRelErr("Thyra::norm_1(nev1)",nev1_nrm,"0",Scalar(0),"max_rel_err",max_rel_err,"max_rel_warn",max_rel_warn,out.ptr())) success=false;
338  if(verbose && dumpAll) *out << "\nnev2 =\n" << *ev1;
339  if(!testRelErr("Thyra::norm_1(nev2)",nev2_nrm,s1_n,s1,"max_rel_err",max_rel_err,"max_rel_warn",max_rel_warn,out.ptr())) success=false;
340  if(verbose && dumpAll) *out << "\nnev2 =\n" << *nev2;
341  if(!testRelErr("Thyra::norm_1(eV1)",eV1_nrm,"0",Scalar(0),"max_rel_err",max_rel_err,"max_rel_warn",max_rel_warn,out.ptr())) success=false;
342  if(verbose && dumpAll) *out << "\neV1 =\n" << *eV1;
343  if(!testRelErr("Thyra::norm_1(eV2)",eV2_nrm,s1_n,s1,"max_rel_err",max_rel_err,"max_rel_warn",max_rel_warn,out.ptr())) success=false;
344  if(verbose && dumpAll) *out << "\neV2 =\n" << *eV2;
345  if(!testRelErr("Thyra::norm_1(neV1)",neV1_nrm,"0",Scalar(0),"max_rel_err",max_rel_err,"max_rel_warn",max_rel_warn,out.ptr())) success=false;
346  if(verbose && dumpAll) *out << "\nneV1 =\n" << *neV1;
347  if(!testRelErr("Thyra::norm_1(neV2)",neV2_nrm,s1_n,s1,"max_rel_err",max_rel_err,"max_rel_warn",max_rel_warn,out.ptr())) success=false;
348  if(verbose && dumpAll) *out << "\nneV2 =\n" << *neV2;
349 
350  if(verbose) *out << "\n*** (B.2) Test RTOps with two or more arguments\n";
351 
352  if(verbose) *out << "\nPerforming ev1 = ev2 ...\n";
353  timer.start(true);
354  Thyra::assign( ev1.ptr(), *ev2 );
355  timer.stop();
356  if(verbose) *out << " time = " << timer.totalElapsedTime() << " sec\n";
357  if(!testRelErr("Thyra::norm_1(ev1)",Thyra::norm_1(*ev1),"Thyra::norm_1(ev2)",ev2_nrm,"max_rel_err",max_rel_err,"max_rel_warn",max_rel_warn,out.ptr())) success=false;
358  if(verbose && dumpAll) *out << "\nev1 =\n" << *ev1;
359 
360  if(verbose) *out << "\nPerforming eV1 = eV2 ...\n";
361  timer.start(true);
362  Thyra::assign( eV1.ptr(), *eV2 );
363  timer.stop();
364  if(verbose) *out << " time = " << timer.totalElapsedTime() << " sec\n";
365  if(!testRelErr("Thyra::norm_1(eV1)",Thyra::norm_1(*eV1),"Thyra::norm_1(eV2)",eV2_nrm,"max_rel_err",max_rel_err,"max_rel_warn",max_rel_warn,out.ptr())) success=false;
366  if(verbose && dumpAll) *out << "\neV1 =\n" << *eV1;
367 
368  if(verbose) *out << "\nPerforming ev1 = nev2 ...\n";
369  timer.start(true);
370  Thyra::assign( ev1.ptr(), *nev2 );
371  timer.stop();
372  if(verbose) *out << " time = " << timer.totalElapsedTime() << " sec\n";
373  if(!testRelErr("Thyra::norm_1(ev1)",Thyra::norm_1(*ev1),"Thyra::norm_1(nev2)",nev2_nrm,"max_rel_err",max_rel_err,"max_rel_warn",max_rel_warn,out.ptr())) success=false;
374  if(verbose && dumpAll) *out << "\nev1 =\n" << *ev1;
375 
376  if(verbose) *out << "\nPerforming nev1 = ev2 ...\n";
377  timer.start(true);
378  Thyra::assign( nev1.ptr(), *ev2 );
379  timer.stop();
380  if(verbose) *out << " time = " << timer.totalElapsedTime() << " sec\n";
381  if(!testRelErr("Thyra::norm_1(nev1)",Thyra::norm_1(*nev1),"Thyra::norm_1(ev2)",ev2_nrm,"max_rel_err",max_rel_err,"max_rel_warn",max_rel_warn,out.ptr())) success=false;
382  if(verbose && dumpAll) *out << "\nnev1 =\n" << *nev1;
383 
384  if(verbose) *out << "\nPerforming nev1 = nev2 ...\n";
385  timer.start(true);
386  Thyra::assign( nev1.ptr(), *nev2 );
387  timer.stop();
388  if(verbose) *out << " time = " << timer.totalElapsedTime() << " sec\n";
389  if(!testRelErr("Thyra::norm_1(nev1)",Thyra::norm_1(*nev1),"Thyra::norm_1(nev2)",nev2_nrm,"max_rel_err",max_rel_err,"max_rel_warn",max_rel_warn,out.ptr())) success=false;
390  if(verbose && dumpAll) *out << "\nnev1 =\n" << *nev1;
391 
392  if(verbose) *out << "\nPerforming eV1 = neV2 ...\n";
393  timer.start(true);
394  Thyra::assign( eV1.ptr(), *neV2 );
395  timer.stop();
396  if(verbose) *out << " time = " << timer.totalElapsedTime() << " sec\n";
397  if(!testRelErr("Thyra::norm_1(eV1)",Thyra::norm_1(*eV1),"Thyra::norm_1(neV2)",neV2_nrm,"max_rel_err",max_rel_err,"max_rel_warn",max_rel_warn,out.ptr())) success=false;
398  if(verbose && dumpAll) *out << "\neV1 =\n" << *eV1;
399 
400  if(verbose) *out << "\nPerforming neV1 = eV2 ...\n";
401  timer.start(true);
402  Thyra::assign( neV1.ptr(), *eV2 );
403  timer.stop();
404  if(verbose) *out << " time = " << timer.totalElapsedTime() << " sec\n";
405  if(!testRelErr("Thyra::norm_1(neV1)",Thyra::norm_1(*neV1),"Thyra::norm_1(eV2)",eV2_nrm,"max_rel_err",max_rel_err,"max_rel_warn",max_rel_warn,out.ptr())) success=false;
406  if(verbose && dumpAll) *out << "\nneV1 =\n" << *neV1;
407 
408  if(verbose) *out << "\nPerforming neV1 = neV2 ...\n";
409  timer.start(true);
410  Thyra::assign( neV1.ptr(), *neV2 );
411  timer.stop();
412  if(verbose) *out << " time = " << timer.totalElapsedTime() << " sec\n";
413  if(!testRelErr("Thyra::norm_1(neV1)",Thyra::norm_1(*neV1),"Thyra::norm_1(neV2)",neV2_nrm,"max_rel_err",max_rel_err,"max_rel_warn",max_rel_warn,out.ptr())) success=false;
414  if(verbose && dumpAll) *out << "\nneV1 =\n" << *neV1;
415 
416  Thyra::LinearOpTester<Scalar> linearOpTester;
417  linearOpTester.set_all_warning_tol(max_rel_warn);
418  linearOpTester.set_all_error_tol(max_rel_err);
419  linearOpTester.dump_all(dumpAll);
420 
421  Thyra::LinearOpWithSolveTester<Scalar> linearOpWithSolveTester;
422  linearOpWithSolveTester.set_all_solve_tol(max_rel_err);
423  linearOpWithSolveTester.set_all_slack_error_tol(max_rel_err);
424  linearOpWithSolveTester.set_all_slack_warning_tol(max_rel_warn);
425  linearOpWithSolveTester.dump_all(dumpAll);
426 
427 
428  if(verbose) *out << "\n*** (B.3) Test Vector linear operator interface\n";
429 
430  if(verbose) *out << "\nChecking *out linear operator interface of ev1 ...\n";
431  result = linearOpTester.check(*ev1,out.ptr());
432  if(!result) success = false;
433 
434  if(verbose) *out << "\nChecking *out linear operator interface of nev1 ...\n";
435  result = linearOpTester.check(*nev1,out.ptr());
436  if(!result) success = false;
437 
438 
439  if(verbose) *out << "\n*** (B.4) Test MultiVector linear operator interface\n";
440 
441  if(verbose) *out << "\nChecking *out linear operator interface of eV1 ...\n";
442  result = linearOpTester.check(*eV1,out.ptr());
443  if(!result) success = false;
444 
445  if(verbose) *out << "\nChecking *out linear operator interface of neV1 ...\n";
446  result = linearOpTester.check(*neV1,out.ptr());
447  if(!result) success = false;
448 
449  const std::string s2_n = "scalar^2*global_dim*num_mv_cols";
450  const Scalar s2 = scalar*scalar*global_dim*num_mv_cols;
451 
452  RCP<Thyra::MultiVectorBase<Scalar> >
453  T = createMembers(eV1->domain(),num_mv_cols);
454 
455 
456  if(verbose) *out << "\n*** (B.5) Test MultiVector::apply(...)\n";
457 
458  if(verbose) *out << "\nPerforming eV1'*eV2 ...\n";
459  timer.start(true);
460  apply( *eV1, TRANS, *eV2, T.ptr() );
461  timer.stop();
462  if(verbose) *out << " time = " << timer.totalElapsedTime() << " sec\n";
463  if(!testRelErr("Thyra::norm_1(eV1'*eV2)",Thyra::norm_1(*T),s2_n,s2,"max_rel_err",max_rel_err,"max_rel_warn",max_rel_warn,out.ptr())) success=false;
464  if(verbose && dumpAll) *out << "\neV1'*eV2 =\n" << *T;
465 
466  if(verbose) *out << "\nPerforming neV1'*eV2 ...\n";
467  timer.start(true);
468  apply( *neV1, TRANS, *eV2, T.ptr() );
469  timer.stop();
470  if(verbose) *out << " time = " << timer.totalElapsedTime() << " sec\n";
471  if(!testRelErr("Thyra::norm_1(neV1'*eV2)",Thyra::norm_1(*T),s2_n,s2,"max_rel_err",max_rel_err,"max_rel_warn",max_rel_warn,out.ptr())) success=false;
472  if(verbose && dumpAll) *out << "\nneV1'*eV2 =\n" << *T;
473 
474  if(verbose) *out << "\nPerforming eV1'*neV2 ...\n";
475  timer.start(true);
476  apply( *eV1, TRANS, *neV2, T.ptr() );
477  timer.stop();
478  if(verbose) *out << " time = " << timer.totalElapsedTime() << " sec\n";
479  if(!testRelErr("Thyra::norm_1(eV1'*neV2)",Thyra::norm_1(*T),s2_n,s2,"max_rel_err",max_rel_err,"max_rel_warn",max_rel_warn,out.ptr())) success=false;
480  if(verbose && dumpAll) *out << "\neV1'*neV2 =\n" << *T;
481 
482  if(verbose) *out << "\nPerforming neV1'*neV2 ...\n";
483  timer.start(true);
484  apply( *neV1, TRANS, *neV2, T.ptr() );
485  timer.stop();
486  if(verbose) *out << " time = " << timer.totalElapsedTime() << " sec\n";
487  if(!testRelErr("Thyra::norm_1(neV1'*neV2)",Thyra::norm_1(*T),s2_n,s2,"max_rel_err",max_rel_err,"max_rel_warn",max_rel_warn,out.ptr())) success=false;
488  if(verbose && dumpAll) *out << "\nneV1'*neV2 =\n" << *T;
489 
490 
491  if(verbose) *out << "\n*** (B.6) Creating a diagonal Epetra_Operator Op\n";
492 
493  RCP<Epetra_Operator> epetra_op;
494 
495  {
496  // Create a diagonal matrix with scalar on the diagonal
497  RCP<Epetra_CrsMatrix>
498  epetra_mat = rcp(new Epetra_CrsMatrix(::Copy,*epetra_map,1));
499  Scalar values[1] = { scalar };
500  int indices[1];
501  const int IB = epetra_map->IndexBase(), offset = procRank*local_dim;
502  for( int k = 0; k < local_dim; ++k ) {
503  indices[0] = offset + k + IB; // global column
504  epetra_mat->InsertGlobalValues(
505  offset + k + IB // GlobalRow
506  ,1 // NumEntries
507  ,values // Values
508  ,indices // Indices
509  );
510  }
511  epetra_mat->FillComplete();
512  epetra_op = epetra_mat;
513  } // end epetra_op
514 
515  RCP<const Thyra::LinearOpBase<Scalar> >
516  Op = Thyra::epetraLinearOp(epetra_op);
517 
518  if(verbose && dumpAll) *out << "\nOp=\n" << *Op;
519 
520 
521  if(verbose) *out << "\n*** (B.6b) Going through partial then full initialization of EpetraLinearOp ...\n";
522 
523  {
524  if(verbose) *out
525  << "\nChecking isFullyUninitialized(*nonconstEpetraLinearOp())==true : ";
526  RCP<Thyra::EpetraLinearOp>
527  thyraEpetraOp = Thyra::nonconstEpetraLinearOp();
528  result = isFullyUninitialized(*thyraEpetraOp);
529  if(!result) success = false;
530  if(verbose) *out << Thyra::passfail(result) << endl;
531  }
532 
533  {
534 
535  if(verbose) *out
536  << "\nthyraEpetraOp = partialNonconstEpetraLinearOp(...)\n";
537  RCP<Thyra::EpetraLinearOp> thyraEpetraOp =
538  Thyra::partialNonconstEpetraLinearOp(
539  epetra_vs, epetra_vs, epetra_op
540  );
541 
542  if(verbose) *out
543  << "\nChecking isPartiallyInitialized(*thyraEpetraOp)==true : ";
544  result = isPartiallyInitialized(*thyraEpetraOp);
545  if(!result) success = false;
546  if(verbose) *out << Thyra::passfail(result) << endl;
547 
548  if(verbose) *out
549  << "\nthyraEpetraOp->setFullyInitialized(true)\n";
550  thyraEpetraOp->setFullyInitialized(true);
551 
552  if(verbose) *out
553  << "\nChecking isFullyInitialized(*thyraEpetraOp)==true : ";
554  result = isFullyInitialized(*thyraEpetraOp);
555  if(!result) success = false;
556  if(verbose) *out << Thyra::passfail(result) << endl;
557 
558  }
559 
560 
561  if(verbose) *out << "\n*** (B.7) Test EpetraLinearOp linear operator interface\n";
562 
563  if(verbose) *out << "\nChecking *out linear operator interface of Op ...\n";
564  result = linearOpTester.check(*Op,out.ptr());
565  if(!result) success = false;
566 
567  RCP<Thyra::VectorBase<Scalar> >
568  ey = createMember(epetra_vs);
569  RCP<Thyra::MultiVectorBase<Scalar> >
570  eY = createMembers(epetra_vs,num_mv_cols);
571  RCP<Thyra::VectorBase<Scalar> >
572  ney = createMember(non_epetra_vs);
573  RCP<Thyra::MultiVectorBase<Scalar> >
574  neY = createMembers(non_epetra_vs,num_mv_cols);
575 
576 
577  if(verbose) *out << "\n*** (B.8) Mix and match vector and Multi-vectors with Epetra opeator\n";
578 
579  const std::string s3_n = "2*scalar^2*global_dim";
580  const Scalar s3 = 2*scalar*scalar*global_dim;
581 
582  if(verbose) *out << "\nPerforming ey = 2*Op*ev1 ...\n";
583  timer.start(true);
584  apply( *Op, NOTRANS, *ev1, ey.ptr(), 2.0 );
585  timer.stop();
586  if(verbose) *out << " time = " << timer.totalElapsedTime() << " sec\n";
587  if(!testRelErr("Thyra::norm_1(ey)",Thyra::norm_1(*ey),s3_n,s3,"max_rel_err",max_rel_err,"max_rel_warn",max_rel_warn,out.ptr())) success=false;
588 
589  if(verbose) *out << "\nPerforming eY = 2*Op*eV1 ...\n";
590  timer.start(true);
591  apply( *Op, NOTRANS, *eV1, eY.ptr(), 2.0 );
592  timer.stop();
593  if(verbose) *out << " time = " << timer.totalElapsedTime() << " sec\n";
594  if(!testRelErr("Thyra::norm_1(eY)",Thyra::norm_1(*eY),s3_n,s3,"max_rel_err",max_rel_err,"max_rel_warn",max_rel_warn,out.ptr())) success=false;
595 
596  if(verbose) *out << "\nPerforming ney = 2*Op*ev1 ...\n";
597  timer.start(true);
598  apply( *Op, NOTRANS, *ev1, ney.ptr(), 2.0 );
599  timer.stop();
600  if(verbose) *out << " time = " << timer.totalElapsedTime() << " sec\n";
601  if(!testRelErr("Thyra::norm_1(ney)",Thyra::norm_1(*ney),s3_n,s3,"max_rel_err",max_rel_err,"max_rel_warn",max_rel_warn,out.ptr())) success=false;
602 
603  if(verbose) *out << "\nPerforming neY = 2*Op*eV1 ...\n";
604  timer.start(true);
605  apply( *Op, NOTRANS, *eV1, neY.ptr(), 2.0 );
606  timer.stop();
607  if(verbose) *out << " time = " << timer.totalElapsedTime() << " sec\n";
608  if(!testRelErr("Thyra::norm_1(neY)",Thyra::norm_1(*neY),s3_n,s3,"max_rel_err",max_rel_err,"max_rel_warn",max_rel_warn,out.ptr())) success=false;
609 
610  if(verbose) *out << "\nPerforming ey = 2*Op*nev1 ...\n";
611  timer.start(true);
612  apply( *Op, NOTRANS, *nev1, ey.ptr(), 2.0 );
613  timer.stop();
614  if(verbose) *out << " time = " << timer.totalElapsedTime() << " sec\n";
615  if(!testRelErr("Thyra::norm_1(ey)",Thyra::norm_1(*ey),s3_n,s3,"max_rel_err",max_rel_err,"max_rel_warn",max_rel_warn,out.ptr())) success=false;
616 
617  if(verbose) *out << "\nPerforming eY = 2*Op*neV1 ...\n";
618  timer.start(true);
619  apply( *Op, NOTRANS, *neV1, eY.ptr(), 2.0 );
620  timer.stop();
621  if(verbose) *out << " time = " << timer.totalElapsedTime() << " sec\n";
622  if(!testRelErr("Thyra::norm_1(eY)",Thyra::norm_1(*eY),s3_n,s3,"max_rel_err",max_rel_err,"max_rel_warn",max_rel_warn,out.ptr())) success=false;
623 
624  if(verbose) *out << "\nPerforming ney = 2*Op*nev1 ...\n";
625  timer.start(true);
626  apply( *Op, NOTRANS, *nev1, ney.ptr(), 2.0 );
627  timer.stop();
628  if(verbose) *out << " time = " << timer.totalElapsedTime() << " sec\n";
629  if(!testRelErr("Thyra::norm_1(ney)",Thyra::norm_1(*ney),s3_n,s3,"max_rel_err",max_rel_err,"max_rel_warn",max_rel_warn,out.ptr())) success=false;
630 
631  if(verbose) *out << "\nPerforming ney = 2*Op*nev1 through MultiVector interface ...\n";
632  timer.start(true);
633  apply( *Op, NOTRANS, static_cast<const Thyra::MultiVectorBase<Scalar>&>(*nev1), Ptr<Thyra::MultiVectorBase<Scalar> >(ney.ptr()), 2.0 );
634  timer.stop();
635  if(verbose) *out << " time = " << timer.totalElapsedTime() << " sec\n";
636  if(!testRelErr("Thyra::norm_1(ney)",Thyra::norm_1(*ney),s3_n,s3,"max_rel_err",max_rel_err,"max_rel_warn",max_rel_warn,out.ptr())) success=false;
637 
638  if(verbose) *out << "\nPerforming neY = 2*Op*neV1 ...\n";
639  timer.start(true);
640  apply( *Op, NOTRANS, *neV1, neY.ptr(), 2.0 );
641  timer.stop();
642  if(verbose) *out << " time = " << timer.totalElapsedTime() << " sec\n";
643  if(!testRelErr("Thyra::norm_1(neY)",Thyra::norm_1(*neY),s3_n,s3,"max_rel_err",max_rel_err,"max_rel_warn",max_rel_warn,out.ptr())) success=false;
644 
645 
646  if(verbose) *out << "\n*** (B.9) Testing Multi-vector views with Epetra operator\n";
647 
648  const Thyra::Range1D col_rng(0,1);
649  const Teuchos::Tuple<int, 2> cols = Teuchos::tuple<int>(2, 3);
650 
651  RCP<const Thyra::MultiVectorBase<Scalar> >
652  eV1_v1 = rcp_static_cast<const Thyra::MultiVectorBase<Scalar> >(eV1)->subView(col_rng),
653  eV1_v2 = rcp_static_cast<const Thyra::MultiVectorBase<Scalar> >(eV1)->subView(cols);
654  RCP<const Thyra::MultiVectorBase<Scalar> >
655  neV1_v1 = rcp_static_cast<const Thyra::MultiVectorBase<Scalar> >(neV1)->subView(col_rng),
656  neV1_v2 = rcp_static_cast<const Thyra::MultiVectorBase<Scalar> >(neV1)->subView(cols);
657  if(verbose && dumpAll) {
658  *out << "\neV1_v1=\n" << *eV1_v1;
659  *out << "\neV1_v2=\n" << *eV1_v2;
660  *out << "\nneV1_v1=\n" << *neV1_v1;
661  *out << "\nneV1_v2=\n" << *neV1_v2;
662  }
663 
664  if(verbose) *out << "\nPerforming eY_v1 = 2*Op*eV1_v1 ...\n";
665  timer.start(true);
666  apply( *Op, NOTRANS, *eV1_v1, eY->subView(col_rng).ptr(), 2.0 );
667  timer.stop();
668  if(verbose) *out << " time = " << timer.totalElapsedTime() << " sec\n";
669  if(verbose && dumpAll) *out << "\neV_v1=\n" << *eY->subView(col_rng);
670  if(!testRelErr("Thyra::norm_1(eY_v1)",Thyra::norm_1(*eY->subView(col_rng)),s3_n,s3,"max_rel_err",max_rel_err,"max_rel_warn",max_rel_warn,out.ptr())) success=false;
671 
672  if(verbose) *out << "\nPerforming eY_v2 = 2*Op*eV1_v2 ...\n";
673  timer.start(true);
674  apply( *Op, NOTRANS, *eV1_v2, eY->subView(cols).ptr(), 2.0 );
675  timer.stop();
676  if(verbose) *out << " time = " << timer.totalElapsedTime() << " sec\n";
677  if(verbose && dumpAll) *out << "\neV_v2=\n" << *eY->subView(cols);
678  if(!testRelErr("Thyra::norm_1(eY_v2)",Thyra::norm_1(*eY->subView(cols)),s3_n,s3,"max_rel_err",max_rel_err,"max_rel_warn",max_rel_warn,out.ptr())) success=false;
679 
680  if(verbose) *out << "\nPerforming neY_v1 = 2*Op*eV1_v1 ...\n";
681  timer.start(true);
682  apply( *Op, NOTRANS, *eV1_v1, neY->subView(col_rng).ptr(), 2.0 );
683  timer.stop();
684  if(verbose) *out << " time = " << timer.totalElapsedTime() << " sec\n";
685  if(verbose && dumpAll) *out << "\neV_v1=\n" << *eY->subView(col_rng);
686  if(!testRelErr("Thyra::norm_1(neY_v1)",Thyra::norm_1(*neY->subView(col_rng)),s3_n,s3,"max_rel_err",max_rel_err,"max_rel_warn",max_rel_warn,out.ptr())) success=false;
687 
688  if(verbose) *out << "\nPerforming eY_v1 = 2*Op*neV1_v1 ...\n";
689  timer.start(true);
690  apply( *Op, NOTRANS, *neV1_v1, eY->subView(col_rng).ptr(), 2.0 );
691  timer.stop();
692  if(verbose) *out << " time = " << timer.totalElapsedTime() << " sec\n";
693  if(!testRelErr("Thyra::norm_1(eY_v1)",Thyra::norm_1(*eY->subView(col_rng)),s3_n,s3,"max_rel_err",max_rel_err,"max_rel_warn",max_rel_warn,out.ptr())) success=false;
694 
695  if(verbose) *out << "\nPerforming neY_v2 = 2*Op*eV1_v2 ...\n";
696  timer.start(true);
697  apply( *Op, NOTRANS, *eV1_v2, neY->subView(cols).ptr(), 2.0 );
698  timer.stop();
699  if(verbose) *out << " time = " << timer.totalElapsedTime() << " sec\n";
700  if(verbose && dumpAll) *out << "\neV_v2=\n" << *eY->subView(cols);
701  if(!testRelErr("Thyra::norm_1(neY_v2)",Thyra::norm_1(*neY->subView(cols)),s3_n,s3,"max_rel_err",max_rel_err,"max_rel_warn",max_rel_warn,out.ptr())) success=false;
702 
703  if(verbose) *out << "\nPerforming eY_v2 = 2*Op*neV1_v2 ...\n";
704  timer.start(true);
705  apply( *Op, NOTRANS, *neV1_v2, eY->subView(cols).ptr(), 2.0 );
706  timer.stop();
707  if(verbose) *out << " time = " << timer.totalElapsedTime() << " sec\n";
708  if(verbose && dumpAll) *out << "\neV_v2=\n" << *eY->subView(cols);
709  if(!testRelErr("Thyra::norm_1(eY_v2)",Thyra::norm_1(*eY->subView(cols)),s3_n,s3,"max_rel_err",max_rel_err,"max_rel_warn",max_rel_warn,out.ptr())) success=false;
710 
711 
712  if(verbose) *out << "\n*** (B.10) Testing Vector and MultiVector view creation functions\n";
713 
714  {
715 
716  const std::string s_n = "fabs(scalar)*num_mv_cols";
717  const Scalar s = fabs(scalar)*num_mv_cols;
718 
719  Array<Scalar> t_raw_values( num_mv_cols );
720  RTOpPack::SubVectorView<Scalar> t_raw( 0, num_mv_cols,
721  arcpFromArray(t_raw_values), 1 );
722 
723  std::fill_n( t_raw_values.begin(), t_raw_values.size(), ST::zero() );
724  Thyra::assign( createMemberView(T->range(),t_raw).ptr(), scalar );
725  Teuchos::RCP<const Thyra::VectorBase<Scalar> > t_view = createMemberView(T->range(),static_cast<RTOpPack::ConstSubVectorView<Scalar>&>(t_raw));
726  Scalar t_nrm = Thyra::norm_1(*t_view);
727  if(!testRelErr("Thyra::norm_1(t_view)",t_nrm,s_n,s,"max_rel_err",max_rel_err,"max_rel_warn",max_rel_warn,out.ptr())) success=false;
728  if(verbose && dumpAll) *out << "\nt_view =\n" << *t_view;
729 
730 /*
731 #ifndef SUN_CXX // The sun compiler Forte Developer 5.4 does not destory temporaries properly and this does not work
732  std::fill_n( t_raw_values.begin(), t_raw_values.size(), ST::zero() );
733  Thyra::assign( T->range().ptr()->Thyra::VectorSpaceBase<Scalar>::createMemberView(t_raw), scalar );
734  t_view = T->range()->Thyra::VectorSpaceBase<Scalar>::createMemberView(static_cast<RTOpPack::ConstSubVectorView<Scalar>&>(t_raw));
735  t_nrm = Thyra::norm_1(*t_view);
736  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;
737  if(verbose && dumpAll) *out << "\nt_view =\n" << *t_view;
738 #endif
739 */
740 
741  Array<Scalar> T_raw_values( num_mv_cols * num_mv_cols );
742  RTOpPack::SubMultiVectorView<Scalar> T_raw( 0, num_mv_cols, 0, num_mv_cols,
743  arcpFromArray(T_raw_values), num_mv_cols );
744 
745  std::fill_n( T_raw_values.begin(), T_raw_values.size(), ST::zero() );
746  Thyra::assign( createMembersView(T->range(),T_raw).ptr(), scalar );
747  Teuchos::RCP<const Thyra::MultiVectorBase<Scalar> >
748  T_view = createMembersView(T->range(),static_cast<RTOpPack::ConstSubMultiVectorView<Scalar>&>(T_raw));
749  Scalar T_nrm = Thyra::norm_1(*T_view);
750  if(!testRelErr("Thyra::norm_1(T_view)",T_nrm,s_n,s,"max_rel_err",max_rel_err,"max_rel_warn",max_rel_warn,out.ptr())) success=false;
751  if(verbose && dumpAll) *out << "\nT_view =\n" << *T_view;
752 
753 /*
754 #ifndef SUN_CXX // The sun compiler Forte Developer 5.4 does not destory temporaries properly and this does not work
755  std::fill_n( T_raw_values.begin(), T_raw_values.size(), ST::zero() );
756  Thyra::assign( T->range().ptr()->Thyra::VectorSpaceBase<Scalar>::createMembersView(T_raw), scalar );
757  T_view = T->range()->Thyra::VectorSpaceBase<Scalar>::createMembersView(static_cast<RTOpPack::ConstSubMultiVectorView<Scalar>&>(T_raw));
758  T_nrm = Thyra::norm_1(*T_view);
759  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;
760  if(verbose && dumpAll) *out << "\nT_view =\n" << *T_view;
761 #endif
762 */
763 
764  }
765 
766 
767  if(verbose) *out << "\n*** (B.11) Testing Epetra_Vector and Epetra_MultiVector wrappers\n";
768 
769  {
770 
771  Teuchos::RCP<const Thyra::SpmdVectorSpaceBase<Scalar> >
772  mpi_vs = Teuchos::rcp_dynamic_cast<const Thyra::SpmdVectorSpaceBase<Scalar> >(epetra_vs,true);
773 
774  if(verbose) *out << "\nCreating temporary Epetra_Vector et1 and Epetra_MultiVector eT1 objects ...\n";
775  Teuchos::RCP<Epetra_Vector>
776  et1 = Teuchos::rcp(new Epetra_Vector(*epetra_map));
777  Teuchos::RCP<Epetra_MultiVector>
778  eT1 = Teuchos::rcp(new Epetra_MultiVector(*epetra_map,num_mv_cols));
779 
780  if(verbose) *out << "\nCreating non-const VectorBase t1 and MultiVectorBase T1 objects from et1 and eT2 ...\n";
781  Teuchos::RCP<Thyra::VectorBase<Scalar> >
782  t1 = create_Vector(et1,mpi_vs);
783  Teuchos::RCP<Thyra::MultiVectorBase<Scalar> >
784  T1 = create_MultiVector(eT1,mpi_vs);
785 
786  if(verbose) *out << "\nPerforming t1 = ev1 ...\n";
787  assign( t1.ptr(), *ev1 );
788  if(!testRelErr(
789  "sum(t1)",Thyra::sum(*t1),"sum(ev1)",sum(*ev1)
790  ,"max_rel_err",max_rel_err,"max_rel_warn",max_rel_warn,out.ptr())
791  ) success=false;
792 
793  if(verbose) *out << "\nPerforming T1 = eV1 ...\n";
794  assign( T1.ptr(), *eV1 );
795  if(!testRelErr(
796  "norm_1(T1)",Thyra::norm_1(*T1),"norm_1(eV1)",norm_1(*eV1)
797  ,"max_rel_err",max_rel_err,"max_rel_warn",max_rel_warn,out.ptr())
798  ) success=false;
799 
800  if(verbose) *out << "\nChecking that t1 and T1 yield the same objects as et1 and eT2 ...\n";
801 
802  Teuchos::RCP<Epetra_Vector>
803  et1_v = get_Epetra_Vector(*epetra_map,t1);
804  result = et1_v.get() == et1.get();
805  if(verbose) *out << "\net1_v.get() = " << et1_v.get() << " == et1.get() = " << et1.get() << " : " << passfail(result) << endl;
806  if(!result) success = false;
807 
808  Teuchos::RCP<Epetra_MultiVector>
809  eT1_v = get_Epetra_MultiVector(*epetra_map,T1);
810  result = eT1_v.get() == eT1.get();
811  if(verbose) *out << "\neT1_v.get() = " << eT1_v.get() << " == eT1.get() = " << eT1.get() << " : " << passfail(result) << endl;
812  if(!result) success = false;
813 
814  if(verbose) *out << "\nCreating const VectorBase ct1 and MultiVectorBase cT1 objects from et1 and eT2 ...\n";
815  Teuchos::RCP<const Thyra::VectorBase<Scalar> >
816  ct1 = create_Vector(Teuchos::rcp_implicit_cast<const Epetra_Vector>(et1),mpi_vs);
817  Teuchos::RCP<const Thyra::MultiVectorBase<Scalar> >
818  cT1 = create_MultiVector(Teuchos::rcp_implicit_cast<const Epetra_MultiVector>(eT1),mpi_vs);
819 
820  if(!testRelErr(
821  "sum(ct1)",Thyra::sum(*ct1),"sum(ev1)",sum(*ev1)
822  ,"max_rel_err",max_rel_err,"max_rel_warn",max_rel_warn,out.ptr())
823  ) success=false;
824 
825  if(!testRelErr(
826  "norm_1(cT1)",Thyra::norm_1(*cT1),"norm_1(eV1)",norm_1(*eV1)
827  ,"max_rel_err",max_rel_err,"max_rel_warn",max_rel_warn,out.ptr())
828  ) success=false;
829 
830  if(verbose) *out << "\nChecking that ct1 and cT1 yield the same objects as et1 and eT2 ...\n";
831 
832  Teuchos::RCP<const Epetra_Vector>
833  cet1_v = get_Epetra_Vector(*epetra_map,ct1);
834  result = cet1_v.get() == et1.get();
835  if(verbose) *out << "\ncet1_v.get() = " << cet1_v.get() << " == et1.get() = " << et1.get() << " : " << passfail(result) << endl;
836  if(!result) success = false;
837 
838  Teuchos::RCP<const Epetra_MultiVector>
839  ceT1_v = get_Epetra_MultiVector(*epetra_map,cT1);
840  result = ceT1_v.get() == eT1.get();
841  if(verbose) *out << "\nceT1_v.get() = " << ceT1_v.get() << " == eT1.get() = " << eT1.get() << " : " << passfail(result) << endl;
842  if(!result) success = false;
843 
844  if(verbose) *out << "\nCreating non-const Epetra_Vector ett1 and Epetra_MultiVector etT1 objects from clones of t1 and T2 ...\n";
845  Teuchos::RCP<Epetra_Vector>
846  ett1 = get_Epetra_Vector(*epetra_map,t1->clone_v());
847  Teuchos::RCP<Epetra_MultiVector>
848  etT1 = get_Epetra_MultiVector(*epetra_map,T1->clone_mv());
849 
850  if(verbose) *out << "\nChecking that ett1 and etT1 yield objects with the same value as et1 and eT2 ...\n";
851 
852  if(!testRelErr(
853  "sum(ett1)",sum(*ett1),"sum(et1)",sum(*et1)
854  ,"max_rel_err",max_rel_err,"max_rel_warn",max_rel_warn,out.ptr())
855  ) success=false;
856 
857  if(!testRelErr(
858  "sum(etT1)",sum(*etT1),"sum(eT1)",sum(*eT1)
859  ,"max_rel_err",max_rel_err,"max_rel_warn",max_rel_warn,out.ptr())
860  ) success=false;
861 
862  if(verbose) *out << "\nCreating const Epetra_Vector cett1 and Epetra_MultiVector cetT1 objects from clones of t1 and T2 ...\n";
863  Teuchos::RCP<const Epetra_Vector>
864  cett1 = get_Epetra_Vector(*epetra_map,Teuchos::rcp_implicit_cast<const Thyra::VectorBase<Scalar> >(t1->clone_v()));
865  Teuchos::RCP<const Epetra_MultiVector>
866  cetT1 = get_Epetra_MultiVector(*epetra_map,Teuchos::rcp_implicit_cast<const Thyra::MultiVectorBase<Scalar> >(T1->clone_mv()));
867 
868  if(verbose) *out << "\nChecking that cett1 and cetT1 yield objects with the same value as et1 and eT2 ...\n";
869 
870  if(!testRelErr(
871  "sum(cett1)",sum(*cett1),"sum(et1)",sum(*et1)
872  ,"max_rel_err",max_rel_err,"max_rel_warn",max_rel_warn,out.ptr())
873  ) success=false;
874 
875  if(!testRelErr(
876  "sum(cetT1)",sum(*cetT1),"sum(eT1)",sum(*eT1)
877  ,"max_rel_err",max_rel_err,"max_rel_warn",max_rel_warn,out.ptr())
878  ) success=false;
879 
880  }
881 
882 
883  if(verbose) *out << "\n*** (B.12) Test DiagonalEpetraLinearOpWithSolveFactory \n";
884 
885  {
886 
887  if(verbose) *out << "\nUsing DiagonalEpetraLinearOpWithSolveFactory to create diagLOWS from Op ...\n";
888 
890 
891  Teuchos::RCP<Thyra::LinearOpWithSolveBase<Scalar> >
892  diagLOWS = Thyra::linearOpWithSolve<double>(diagLOWSFactory,Op);
893 
894  if(verbose) *out << "\nTesting LinearOpBase interface of diagLOWS ...\n";
895 
896  result = linearOpTester.check(*diagLOWS, out.ptr());
897  if(!result) success = false;
898 
899  if(verbose) *out << "\nTesting LinearOpWithSolveBase interface of diagLOWS ...\n";
900 
901  result = linearOpWithSolveTester.check(*diagLOWS, &*out);
902  if(!result) success = false;
903 
904  }
905 
906 
907  if(verbose)
908  *out
909  << "\n***"
910  << "\n*** (C) Comparing the speed of Thyra adapted Epetra objects verses raw Epetra objects"
911  << "\n***\n";
912 
913  //
914  // Setup the number of timing loops to get good timings
915  //
916  // Here we try to shoot for timing ab*out 1 second's worth of
917  // computations and adjust the number of evaluation loops
918  // accordingly. Let X be the approximate number of flops per
919  // loop (or per evaluation). We then compute the number of
920  // loops as:
921  //
922  // 1.0 sec | num CPU flops | 1 loop |
923  // --------|----------------|-----------|
924  // | sec | X flops |
925  //
926  // This just comes *out to be:
927  //
928  // num_time_loops_X = max_flop_rate / (X flops per loop)
929  //
930  // In this computation we ignore extra overhead that will be
931  // an issue when local_dim is small.
932  //
933 
934  double raw_epetra_time, thyra_wrapped_time;
935 
936 
937  if(verbose) *out << "\n*** (C.1) Comparing the speed of RTOp verses raw Epetra_Vector operations\n";
938 
939  const double flop_adjust_factor_1 = 3.0;
940  const int num_time_loops_1 = int( max_flop_rate / ( flop_adjust_factor_1 * local_dim * num_mv_cols ) ) + 1;
941 
942  {
943 
944  // Get references to Epetra_MultiVector objects in eV1 and eV2
945  const RCP<Epetra_MultiVector> eeV1 = get_Epetra_MultiVector(*epetra_map,eV1);
946  const RCP<const Epetra_MultiVector> eeV2 = get_Epetra_MultiVector(*epetra_map,eV2);
947 
948  if(verbose)
949  *out << "\nPerforming eeV1 = eeV2 (using raw Epetra_MultiVector::operator=(...)) " << num_time_loops_1 << " times ...\n";
950  timer.start(true);
951  for(int k = 0; k < num_time_loops_1; ++k ) {
952  *eeV1 = *eeV2;
953  }
954  timer.stop();
955  raw_epetra_time = timer.totalElapsedTime();
956  if(verbose) *out << " total time = " << raw_epetra_time << " sec\n";
957 
958  // When this block ends and eeV1 goes *out of scope then eV1 is guaranteed to be updated!
959  }
960 
961  if(verbose)
962  *out << "\nPerforming eV1 = eV2 (using Thyra::SpmdMultiVectorBase::applyOp(...)) " << num_time_loops_1 << " times ...\n";
963  timer.start(true);
964  for(int k = 0; k < num_time_loops_1; ++k ) {
965  Thyra::assign( eV1.ptr(), *eV2 );
966  }
967  timer.stop();
968  thyra_wrapped_time = timer.totalElapsedTime();
969  if(verbose) *out << " total time = " << thyra_wrapped_time << " sec\n";
970 
971  print_performance_stats( num_time_loops_1, raw_epetra_time, thyra_wrapped_time, verbose, *out );
972 
973  // RAB: 2004/01/05: Note, the above relative performance is likely
974  // to be the worst of all of the others since RTOp operators are
975  // applied seperately column by column but the relative
976  // performance should go to about 1.0 when local_dim is
977  // sufficiently large! However, because
978  // Epetra_MultiVector::Thyra::Assign(...) is implemented using double
979  // pointer indexing, the RTOp implementation used with the Thyra
980  // adapters is actually faster in some cases. However, the extra
981  // overhead of RTOp is much worse for very very small (order 10)
982  // sizes.
983 
984 
985  if(verbose)
986  *out
987  << "\n*** (C.2) Comparing Thyra::SpmdMultiVectorBase::apply() verses raw Epetra_MultiVector::Multiply()\n";
988 
989  Teuchos::TimeMonitor::zeroOutTimers();
990 
991  const double flop_adjust_factor_2 = 2.0;
992  const int num_time_loops_2 = int( max_flop_rate / ( flop_adjust_factor_2* local_dim * num_mv_cols * num_mv_cols ) ) + 1;
993 
994  {
995 
996  // Get constant references to Epetra_MultiVector objects in eV1 and eV2
997  const RCP<const Epetra_MultiVector> eeV1 = get_Epetra_MultiVector(*epetra_map,eV1);
998  const RCP<const Epetra_MultiVector> eeV2 = get_Epetra_MultiVector(*epetra_map,eV2);
999 
1000  Epetra_LocalMap eT_map((int) T->range()->dim(),0,*epetra_comm);
1001  Epetra_MultiVector eT(eT_map,T->domain()->dim());
1002 
1003  if(verbose)
1004  *out << "\nPerforming eeV1'*eeV2 (using raw Epetra_MultiVector::Multiply(...)) " << num_time_loops_2 << " times ...\n";
1005  timer.start(true);
1006  for(int k = 0; k < num_time_loops_2; ++k ) {
1007  eT.Multiply( 'T', 'N', 1.0, *eeV1, *eeV2, 0.0 );
1008  }
1009  timer.stop();
1010  raw_epetra_time = timer.totalElapsedTime();
1011  if(verbose) *out << " total time = " << raw_epetra_time << " sec\n";
1012 
1013  }
1014 
1015  if(verbose)
1016  *out << "\nPerforming eV1'*eV2 (using Thyra::SpmdMultiVectorBase::apply(...)) " << num_time_loops_2 << " times ...\n";
1017  timer.start(true);
1018  for(int k = 0; k < num_time_loops_2; ++k ) {
1019  apply( *eV1, TRANS, *eV2, T.ptr() );
1020  }
1021  timer.stop();
1022  thyra_wrapped_time = timer.totalElapsedTime();
1023  if(verbose) *out << " total time = " << thyra_wrapped_time << " sec\n";
1024 
1025  print_performance_stats( num_time_loops_2, raw_epetra_time, thyra_wrapped_time, verbose, *out );
1026 
1027  // RAB: 2004/01/05: Note, even though the Thyra adapter does
1028  // not actually call Epetra_MultiVector::Multiply(...), the
1029  // implementation in Thyra::SpmdMultiVectorBase::apply(...)
1030  // performs almost exactly the same flops and calls dgemm(...)
1031  // as well. Herefore, except for some small overhead, the raw
1032  // Epetra and the Thyra wrapped computations should give
1033  // almost identical times in almost all cases.
1034 
1035 
1036  if(verbose) *out << "\n*** (C.3) Comparing Thyra::EpetraLinearOp::apply() verses raw Epetra_Operator::apply()\n";
1037 
1038  Teuchos::TimeMonitor::zeroOutTimers();
1039 
1040  const double flop_adjust_factor_3 = 10.0; // lots of indirect addressing
1041  const int num_time_loops_3 = int( max_flop_rate / ( flop_adjust_factor_3 * local_dim * num_mv_cols ) ) + 1;
1042 
1043  {
1044 
1045  // Get constant references to Epetra_MultiVector objects in eV1 and eV2
1046  const RCP<const Epetra_MultiVector> eeV1 = get_Epetra_MultiVector(*epetra_map,eV1);
1047  const RCP<Epetra_MultiVector> eeY = get_Epetra_MultiVector(*epetra_map,eY);
1048 
1049  if(verbose)
1050  *out << "\nPerforming eeY = eOp*eeV1 (using raw Epetra_Operator::apply(...)) " << num_time_loops_3 << " times ...\n";
1051 
1052  Teuchos::TimeMonitor::zeroOutTimers();
1053 
1054  timer.start(true);
1055  epetra_op->SetUseTranspose(false);
1056  for(int k = 0; k < num_time_loops_3; ++k ) {
1057  epetra_op->Apply( *eeV1, *eeY );
1058  //eeY->Scale(2.0);
1059  }
1060  timer.stop();
1061 
1062  raw_epetra_time = timer.totalElapsedTime();
1063  if(verbose) *out << " total time = " << raw_epetra_time << " sec\n";
1064 
1065  if(verbose)
1066  Teuchos::TimeMonitor::summarize(*out << "\n");
1067 
1068  }
1069 
1070  if(verbose)
1071  *out << "\nPerforming eY = Op*eV1 (using Thyra::EpetraLinearOp::apply(...)) " << num_time_loops_3 << " times ...\n";
1072 
1073  Teuchos::TimeMonitor::zeroOutTimers();
1074 
1075  timer.start(true);
1076  for(int k = 0; k < num_time_loops_3; ++k ) {
1077  apply( *Op, NOTRANS, *eV1, eY.ptr() );
1078  }
1079  timer.stop();
1080 
1081  thyra_wrapped_time = timer.totalElapsedTime();
1082  if(verbose) *out << " total time = " << thyra_wrapped_time << " sec\n";
1083 
1084  if(verbose)
1085  Teuchos::TimeMonitor::summarize(*out << "\n");
1086 
1087  print_performance_stats( num_time_loops_3, raw_epetra_time, thyra_wrapped_time, verbose, *out );
1088 
1089  // RAB: 2004/01/05: Note, the above Epetra adapter is a true
1090  // adapter and simply calls Epetra_Operator::apply(...) so, except
1091  // for some small overhead, the raw Epetra and the Thyra wrapped
1092  // computations should give ab*out exactly the same runtime for
1093  // almost all cases.
1094 
1095  } // end try
1096  TEUCHOS_STANDARD_CATCH_STATEMENTS(true,*out,success)
1097 
1098  if(verbose) {
1099  if(success) *out << "\nCongratulations! All of the tests seem to have run sucessfully!\n";
1100  else *out << "\nOh no! at least one of the tests did not check out!\n";
1101  }
1102 
1103  return (success ? 0 : -1);
1104 
1105 } // end main()
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...
TRANS
RCP< const VectorSpaceBase< double > > create_VectorSpace(const RCP< const Epetra_Map > &epetra_map)
Create an VectorSpaceBase object given an Epetra_Map object.
RCP< MultiVectorBase< double > > create_MultiVector(const RCP< Epetra_MultiVector > &epetra_mv, const RCP< const VectorSpaceBase< double > > &range, const RCP< const VectorSpaceBase< double > > &domain=Teuchos::null)
Create a non-const MultiVectorBase object from a non-const Epetra_MultiVector object.
RCP< VectorBase< double > > create_Vector(const RCP< Epetra_Vector > &epetra_v, const RCP< const VectorSpaceBase< double > > &space)
Create a non-const VectorBase object from a non-const Epetra_Vector object.
RCP< Epetra_Vector > get_Epetra_Vector(const Epetra_Map &map, const RCP< VectorBase< double > > &v)
Get a non-const Epetra_Vector view from a non-const VectorBase object if possible.
int main(int argc, char *argv[])
Create a DefaultDiagonalLinearOpWithSolve out of a diagonal Epetra_RowMatrix object.
RCP< const Teuchos::Comm< Ordinal > > create_Comm(const RCP< const Epetra_Comm > &epetraComm)
Given an Epetra_Comm object, return an equivalent Teuchos::Comm object.
NOTRANS