Thyra Package Browser (Single Doxygen Collection)  Version of the Day
EpetraExtAddTransformer_UnitTests.cpp
Go to the documentation of this file.
1 /*
2 // @HEADER
3 // ***********************************************************************
4 //
5 // Thyra: Interfaces and Support for Abstract Numerical Algorithms
6 // Copyright (2004) Sandia Corporation
7 //
8 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
9 // license for use of this work by or on behalf of the U.S. Government.
10 //
11 // Redistribution and use in source and binary forms, with or without
12 // modification, are permitted provided that the following conditions are
13 // met:
14 //
15 // 1. Redistributions of source code must retain the above copyright
16 // notice, this list of conditions and the following disclaimer.
17 //
18 // 2. Redistributions in binary form must reproduce the above copyright
19 // notice, this list of conditions and the following disclaimer in the
20 // documentation and/or other materials provided with the distribution.
21 //
22 // 3. Neither the name of the Corporation nor the names of the
23 // contributors may be used to endorse or promote products derived from
24 // this software without specific prior written permission.
25 //
26 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
27 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
28 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
29 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
30 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
31 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
32 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
33 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
34 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
35 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
36 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
37 //
38 // Questions? Contact Roscoe A. Bartlett (bartlettra@ornl.gov)
39 //
40 // ***********************************************************************
41 // @HEADER
42 */
43 
44 
46 #include "Thyra_DefaultScaledAdjointLinearOp.hpp"
47 #include "Thyra_DefaultMultipliedLinearOp.hpp"
48 #include "Thyra_DefaultAddedLinearOp.hpp"
49 #include "Thyra_DefaultDiagonalLinearOp.hpp"
50 #include "Thyra_VectorStdOps.hpp"
51 #include "Thyra_TestingTools.hpp"
52 #include "Thyra_LinearOpTester.hpp"
54 #include "Thyra_EpetraLinearOp.hpp"
56 #include "Thyra_DefaultIdentityLinearOp.hpp"
57 #include "EpetraExt_readEpetraLinearSystem.h"
58 #include "Teuchos_GlobalMPISession.hpp"
59 #include "Teuchos_VerboseObject.hpp"
60 #include "Teuchos_XMLParameterListHelpers.hpp"
61 #include "Teuchos_CommandLineProcessor.hpp"
62 #include "Teuchos_StandardCatchMacros.hpp"
63 
64 #ifdef HAVE_MPI
65 # include "Epetra_MpiComm.h"
66 #else
67 # include "Epetra_SerialComm.h"
68 #endif
69 
70 #include "Teuchos_UnitTestHarness.hpp"
71 
72 
73 namespace {
74 
75 
76 using Teuchos::null;
77 using Teuchos::RCP;
79 using Thyra::epetraExtAddTransformer;
80 using Thyra::VectorBase;
81 using Thyra::LinearOpBase;
82 using Thyra::createMember;
83 using Thyra::LinearOpTester;
84 using Thyra::adjoint;
85 using Thyra::multiply;
86 using Thyra::diagonal;
87 
88 
89 std::string matrixFile = "";
90 std::string matrixFile2 = "";
91 
92 
94 {
95  Teuchos::UnitTestRepository::getCLP().setOption(
96  "matrix-file", &matrixFile,
97  "Defines the Epetra_CrsMatrix to read in." );
98  Teuchos::UnitTestRepository::getCLP().setOption(
99  "matrix-file-2", &matrixFile2,
100  "Defines the Epetra_CrsMatrix to read in." );
101 }
102 
103 const Teuchos::RCP<const Thyra::LinearOpBase<double> >
104 buildAddOperator(int scenario,const Teuchos::RCP<const Thyra::LinearOpBase<double> > & A,
105  const Teuchos::RCP<const Thyra::LinearOpBase<double> > & B)
106 {
107  // build operators for the various addition/adjoint scenarios
108  RCP<const Thyra::LinearOpBase<double> > M;
109 
110  switch(scenario) {
111  case 0:
112  M = Thyra::add(A,B,"A+B");
113  break;
114  case 1:
115  M = Thyra::add(A,B,"A+adj(B)");
116  break;
117  case 2:
118  M = Thyra::add(A,B,"adj(A)+B");
119  break;
120  case 3:
121  M = Thyra::add(A,B,"adb(A)+adb(B)");
122  break;
123  default:
124  TEUCHOS_ASSERT(false);
125  break;
126  }
127 
128  return M;
129 }
130 
131 TEUCHOS_UNIT_TEST( EpetraExtAddTransformer, diag_mat_Add )
132 {
133 
134  //
135  // A) Read in problem matrices
136  //
137 
138  out << "\nReading linear system in Epetra format from the file \'"<<matrixFile<<"\' ...\n";
139 
140 #ifdef HAVE_MPI
141  Epetra_MpiComm comm(MPI_COMM_WORLD);
142 #else
143  Epetra_SerialComm comm;
144 #endif
145  RCP<Epetra_CrsMatrix> crsMat;
146  EpetraExt::readEpetraLinearSystem( matrixFile, comm, &crsMat, NULL, NULL, NULL );
147 
148  //
149  // B) Create the Thyra wrapped version
150  //
151  double scaleMat=3.7;
152  //double scaleDiag=-2.9;
153 
154 
155  const RCP<const Thyra::LinearOpBase<double> > A = Thyra::scale<double>(scaleMat,Thyra::epetraLinearOp(crsMat));
156  RCP<VectorBase<double> > d = createMember(A->domain(), "d");
157  V_S( d.ptr(), 3.0 ); // ToDo: Set ton != 1.0 and generalize
158  const RCP<const Thyra::LinearOpBase<double> > B = diagonal(d);
159 
160  out << "\nA = " << *A;
161  out << "\nB = " << *B;
162 
163  for(int scenario=0;scenario<4;scenario++) {
164  //
165  // C) Create implicit A+B operator
166  //
167 
168  const RCP<const Thyra::LinearOpBase<double> > M = buildAddOperator(scenario,A,B);
169 
170  //
171  // D) Do the transformation
172  //
173 
174  const RCP<EpetraExtAddTransformer> ApB_transformer = epetraExtAddTransformer();
175 
176  TEST_ASSERT(ApB_transformer != null);
177 
178  const RCP<LinearOpBase<double> > M_explicit = ApB_transformer->createOutputOp();
179  ApB_transformer->transform( *M, M_explicit.ptr() );
180 
181  out << "\nM_explicit = " << *M_explicit;
182 
183  //
184  // E) Check the explicit operator
185  //
186 
187  LinearOpTester<double> M_explicit_tester;
188  M_explicit_tester.show_all_tests(true);;
189 
190  const bool result = M_explicit_tester.compare( *M, *M_explicit, Teuchos::inOutArg(out) );
191  if (!result) success = false;
192  }
193 
194  for(int scenario=0;scenario<4;scenario++) {
195  //
196  // C) Create implicit A+B operator
197  //
198 
199  const RCP<const Thyra::LinearOpBase<double> > M = buildAddOperator(scenario,B,A);
200 
201  //
202  // D) Do the transformation
203  //
204 
205  const RCP<EpetraExtAddTransformer> ApB_transformer = epetraExtAddTransformer();
206 
207  TEST_ASSERT(ApB_transformer != null);
208 
209  const RCP<LinearOpBase<double> > M_explicit = ApB_transformer->createOutputOp();
210  ApB_transformer->transform( *M, M_explicit.ptr() );
211 
212  out << "\nM_explicit = " << *M_explicit;
213 
214  //
215  // E) Check the explicit operator
216  //
217 
218  LinearOpTester<double> M_explicit_tester;
219  M_explicit_tester.show_all_tests(true);;
220 
221  const bool result = M_explicit_tester.compare( *M, *M_explicit, Teuchos::inOutArg(out) );
222  if (!result) success = false;
223  }
224 }
225 
226 TEUCHOS_UNIT_TEST( EpetraExtAddTransformer, id_mat_Add )
227 {
228 
229  //
230  // A) Read in problem matrices
231  //
232 
233  out << "\nReading linear system in Epetra format from the file \'"<<matrixFile<<"\' ...\n";
234 
235 #ifdef HAVE_MPI
236  Epetra_MpiComm comm(MPI_COMM_WORLD);
237 #else
238  Epetra_SerialComm comm;
239 #endif
240  RCP<Epetra_CrsMatrix> crsMat;
241  EpetraExt::readEpetraLinearSystem( matrixFile, comm, &crsMat, NULL, NULL, NULL );
242 
243  //
244  // B) Create the Thyra wrapped version
245  //
246  double scaleMat=3.7;
247  //double scaleDiag=-2.9;
248 
249 
250  const RCP<const Thyra::LinearOpBase<double> > A = Thyra::scale<double>(scaleMat,Thyra::epetraLinearOp(crsMat));
251  const RCP<const Thyra::LinearOpBase<double> > B = identity(A->domain(),"id");
252 
253  out << "\nA = " << *A;
254  out << "\nB = " << *B;
255 
256  for(int scenario=0;scenario<4;scenario++) {
257  //
258  // C) Create implicit A+B operator
259  //
260 
261  const RCP<const Thyra::LinearOpBase<double> > M = buildAddOperator(scenario,A,B);
262 
263  //
264  // D) Do the transformation
265  //
266 
267  const RCP<EpetraExtAddTransformer> ApB_transformer = epetraExtAddTransformer();
268 
269  TEST_ASSERT(ApB_transformer != null);
270 
271  const RCP<LinearOpBase<double> > M_explicit = ApB_transformer->createOutputOp();
272  ApB_transformer->transform( *M, M_explicit.ptr() );
273 
274  out << "\nM_explicit = " << *M_explicit;
275 
276  //
277  // E) Check the explicit operator
278  //
279 
280  LinearOpTester<double> M_explicit_tester;
281  M_explicit_tester.show_all_tests(true);;
282 
283  const bool result = M_explicit_tester.compare( *M, *M_explicit, Teuchos::inOutArg(out) );
284  if (!result) success = false;
285  }
286 
287  for(int scenario=0;scenario<4;scenario++) {
288  //
289  // C) Create implicit A+B operator
290  //
291 
292  const RCP<const Thyra::LinearOpBase<double> > M = buildAddOperator(scenario,B,A);
293 
294  //
295  // D) Do the transformation
296  //
297 
298  const RCP<EpetraExtAddTransformer> ApB_transformer = epetraExtAddTransformer();
299 
300  TEST_ASSERT(ApB_transformer != null);
301 
302  const RCP<LinearOpBase<double> > M_explicit = ApB_transformer->createOutputOp();
303  ApB_transformer->transform( *M, M_explicit.ptr() );
304 
305  out << "\nM_explicit = " << *M_explicit;
306 
307  //
308  // E) Check the explicit operator
309  //
310 
311  LinearOpTester<double> M_explicit_tester;
312  M_explicit_tester.show_all_tests(true);;
313 
314  const bool result = M_explicit_tester.compare( *M, *M_explicit, Teuchos::inOutArg(out) );
315  if (!result) success = false;
316  }
317 }
318 
319 TEUCHOS_UNIT_TEST( EpetraExtAddTransformer, basic_Add )
320 {
321 
322  //
323  // A) Read in problem matrices
324  //
325 
326  out << "\nReading linear system in Epetra format from the file \'"<<matrixFile<<"\' ...\n";
327  out << "\nReading linear system in Epetra format from the file \'"<<matrixFile2<<"\' ...\n";
328 
329 #ifdef HAVE_MPI
330  Epetra_MpiComm comm(MPI_COMM_WORLD);
331 #else
332  Epetra_SerialComm comm;
333 #endif
334  RCP<Epetra_CrsMatrix> epetra_A;
335  RCP<Epetra_CrsMatrix> epetra_B;
336  EpetraExt::readEpetraLinearSystem( matrixFile, comm, &epetra_A, NULL, NULL, NULL );
337  EpetraExt::readEpetraLinearSystem( matrixFile2, comm, &epetra_B, NULL, NULL, NULL );
338 
339  //
340  // B) Create the Thyra wrapped version
341  //
342  double scaleA=3.7;
343  double scaleB=-2.9;
344 
345  const RCP<const Thyra::LinearOpBase<double> > A = Thyra::scale<double>(scaleA,Thyra::epetraLinearOp(epetra_B));
346  const RCP<const Thyra::LinearOpBase<double> > B = Thyra::scale<double>(scaleB,Thyra::epetraLinearOp(epetra_B));
347 
348  out << "\nA = " << *A;
349  out << "\nB = " << *B;
350 
351  for(int scenario=0;scenario<4;scenario++) {
352  //
353  // C) Create implicit A+B operator
354  //
355 
356  const RCP<const Thyra::LinearOpBase<double> > M = buildAddOperator(scenario,A,B);
357 
358  //
359  // D) Do the transformation
360  //
361 
362  const RCP<EpetraExtAddTransformer> ApB_transformer = epetraExtAddTransformer();
363 
364  TEST_ASSERT(ApB_transformer != null);
365 
366  const RCP<LinearOpBase<double> > M_explicit = ApB_transformer->createOutputOp();
367  ApB_transformer->transform( *M, M_explicit.ptr() );
368 
369  out << "\nM_explicit = " << *M_explicit;
370 
371  //
372  // E) Check the explicit operator
373  //
374 
375  LinearOpTester<double> M_explicit_tester;
376  M_explicit_tester.show_all_tests(true);;
377 
378  const bool result = M_explicit_tester.compare( *M, *M_explicit, Teuchos::inOutArg(out) );
379  if (!result) success = false;
380  }
381 }
382 
383 TEUCHOS_UNIT_TEST( EpetraExtAddTransformer, mod_Add )
384 {
385 
386  //
387  // A) Read in problem matrices
388  //
389 
390  out << "\nReading linear system in Epetra format from the file \'"<<matrixFile<<"\' ...\n";
391  out << "\nReading linear system in Epetra format from the file \'"<<matrixFile2<<"\' ...\n";
392 
393 #ifdef HAVE_MPI
394  Epetra_MpiComm comm(MPI_COMM_WORLD);
395 #else
396  Epetra_SerialComm comm;
397 #endif
398  RCP<Epetra_CrsMatrix> epetra_A;
399  RCP<Epetra_CrsMatrix> epetra_B;
400  EpetraExt::readEpetraLinearSystem( matrixFile, comm, &epetra_A, NULL, NULL, NULL );
401  EpetraExt::readEpetraLinearSystem( matrixFile2, comm, &epetra_B, NULL, NULL, NULL );
402 
403  //
404  // B) Create the Thyra wrapped version
405  //
406  double scaleA=3.7;
407  double scaleB=-2.9;
408 
409  const RCP<const Thyra::LinearOpBase<double> > A = Thyra::scale<double>(scaleA,Thyra::epetraLinearOp(epetra_B));
410  const RCP<const Thyra::LinearOpBase<double> > B = Thyra::scale<double>(scaleB,Thyra::epetraLinearOp(epetra_B));
411 
412  out << "\nA = " << *A;
413  out << "\nB = " << *B;
414 
415  for(int scenario=0;scenario<4;scenario++) {
416  //
417  // C) Create implicit A+B operator
418  //
419 
420  const RCP<const Thyra::LinearOpBase<double> > M = buildAddOperator(scenario,A,B);
421 
422  //
423  // D) Do the transformation
424  //
425 
426  const RCP<EpetraExtAddTransformer> ApB_transformer = epetraExtAddTransformer();
427 
428  TEST_ASSERT(ApB_transformer != null);
429 
430  const RCP<LinearOpBase<double> > M_explicit = ApB_transformer->createOutputOp();
431  const RCP<LinearOpBase<double> > M_explicit_orig = M_explicit;
432  ApB_transformer->transform( *M, M_explicit.ptr() );
433 
434  // do some violence to one of the operators
435  double * view; int numEntries;
436  epetra_B->Scale(3.2);
437  TEUCHOS_ASSERT(epetra_B->ExtractMyRowView(3,numEntries,view)==0);
438  for(int i=0;i<numEntries;i++) view[i] += view[i]*double(i+1);
439 
440  // compute multiply again
441  ApB_transformer->transform( *M, M_explicit.ptr() );
442 
443  out << "\nM_explicit = " << *M_explicit;
444 
445  TEUCHOS_ASSERT(Thyra::get_Epetra_Operator(*M_explicit)
446  ==Thyra::get_Epetra_Operator(*M_explicit_orig));
447 
448  //
449  // E) Check the explicit operator
450  //
451 
452  LinearOpTester<double> M_explicit_tester;
453  M_explicit_tester.show_all_tests(true);;
454 
455  const bool result = M_explicit_tester.compare( *M, *M_explicit, Teuchos::inoutArg(out) );
456  if (!result) success = false;
457  }
458 }
459 
460 } // end namespace
Transformer subclass for adding Epetra/Thyra operators using EpetraExt::MatrixMatrix.
TEUCHOS_UNIT_TEST(EpetraOperatorWrapper, basic)
Teuchos::RCP< Epetra_Operator > get_Epetra_Operator(LinearOpBase< double > &op)
Full specialization for Scalar=double.