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