Thyra Package Browser (Single Doxygen Collection)  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
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 
45 #include "Thyra_EpetraExtAddTransformer.hpp"
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"
53 #include "Thyra_EpetraThyraWrappers.hpp"
54 #include "Thyra_EpetraLinearOp.hpp"
55 #include "Thyra_get_Epetra_Operator.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;
78 using Thyra::EpetraExtAddTransformer;
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 
93 TEUCHOS_STATIC_SETUP()
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  }
126 
127  return M;
128 }
129 
130 TEUCHOS_UNIT_TEST( EpetraExtAddTransformer, diag_mat_Add )
131 {
132 
133  //
134  // A) Read in problem matrices
135  //
136 
137  out << "\nReading linear system in Epetra format from the file \'"<<matrixFile<<"\' ...\n";
138 
139 #ifdef HAVE_MPI
140  Epetra_MpiComm comm(MPI_COMM_WORLD);
141 #else
142  Epetra_SerialComm comm;
143 #endif
144  RCP<Epetra_CrsMatrix> crsMat;
145  EpetraExt::readEpetraLinearSystem( matrixFile, comm, &crsMat, NULL, NULL, NULL );
146 
147  //
148  // B) Create the Thyra wrapped version
149  //
150  double scaleMat=3.7;
151  //double scaleDiag=-2.9;
152 
153 
154  const RCP<const Thyra::LinearOpBase<double> > A = Thyra::scale<double>(scaleMat,Thyra::epetraLinearOp(crsMat));
155  RCP<VectorBase<double> > d = createMember(A->domain(), "d");
156  V_S( d.ptr(), 3.0 ); // ToDo: Set ton != 1.0 and generalize
157  const RCP<const Thyra::LinearOpBase<double> > B = diagonal(d);
158 
159  out << "\nA = " << *A;
160  out << "\nB = " << *B;
161 
162  for(int scenario=0;scenario<4;scenario++) {
163  //
164  // C) Create implicit A+B operator
165  //
166 
167  const RCP<const Thyra::LinearOpBase<double> > M = buildAddOperator(scenario,A,B);
168 
169  //
170  // D) Do the transformation
171  //
172 
173  const RCP<EpetraExtAddTransformer> ApB_transformer = epetraExtAddTransformer();
174 
175  TEST_ASSERT(ApB_transformer != null);
176 
177  const RCP<LinearOpBase<double> > M_explicit = ApB_transformer->createOutputOp();
178  ApB_transformer->transform( *M, M_explicit.ptr() );
179 
180  out << "\nM_explicit = " << *M_explicit;
181 
182  //
183  // E) Check the explicit operator
184  //
185 
186  LinearOpTester<double> M_explicit_tester;
187  M_explicit_tester.show_all_tests(true);;
188 
189  const bool result = M_explicit_tester.compare( *M, *M_explicit, Teuchos::inOutArg(out) );
190  if (!result) success = false;
191  }
192 
193  for(int scenario=0;scenario<4;scenario++) {
194  //
195  // C) Create implicit A+B operator
196  //
197 
198  const RCP<const Thyra::LinearOpBase<double> > M = buildAddOperator(scenario,B,A);
199 
200  //
201  // D) Do the transformation
202  //
203 
204  const RCP<EpetraExtAddTransformer> ApB_transformer = epetraExtAddTransformer();
205 
206  TEST_ASSERT(ApB_transformer != null);
207 
208  const RCP<LinearOpBase<double> > M_explicit = ApB_transformer->createOutputOp();
209  ApB_transformer->transform( *M, M_explicit.ptr() );
210 
211  out << "\nM_explicit = " << *M_explicit;
212 
213  //
214  // E) Check the explicit operator
215  //
216 
217  LinearOpTester<double> M_explicit_tester;
218  M_explicit_tester.show_all_tests(true);;
219 
220  const bool result = M_explicit_tester.compare( *M, *M_explicit, Teuchos::inOutArg(out) );
221  if (!result) success = false;
222  }
223 }
224 
225 TEUCHOS_UNIT_TEST( EpetraExtAddTransformer, id_mat_Add )
226 {
227 
228  //
229  // A) Read in problem matrices
230  //
231 
232  out << "\nReading linear system in Epetra format from the file \'"<<matrixFile<<"\' ...\n";
233 
234 #ifdef HAVE_MPI
235  Epetra_MpiComm comm(MPI_COMM_WORLD);
236 #else
237  Epetra_SerialComm comm;
238 #endif
239  RCP<Epetra_CrsMatrix> crsMat;
240  EpetraExt::readEpetraLinearSystem( matrixFile, comm, &crsMat, NULL, NULL, NULL );
241 
242  //
243  // B) Create the Thyra wrapped version
244  //
245  double scaleMat=3.7;
246  //double scaleDiag=-2.9;
247 
248 
249  const RCP<const Thyra::LinearOpBase<double> > A = Thyra::scale<double>(scaleMat,Thyra::epetraLinearOp(crsMat));
250  const RCP<const Thyra::LinearOpBase<double> > B = identity(A->domain(),"id");
251 
252  out << "\nA = " << *A;
253  out << "\nB = " << *B;
254 
255  for(int scenario=0;scenario<4;scenario++) {
256  //
257  // C) Create implicit A+B operator
258  //
259 
260  const RCP<const Thyra::LinearOpBase<double> > M = buildAddOperator(scenario,A,B);
261 
262  //
263  // D) Do the transformation
264  //
265 
266  const RCP<EpetraExtAddTransformer> ApB_transformer = epetraExtAddTransformer();
267 
268  TEST_ASSERT(ApB_transformer != null);
269 
270  const RCP<LinearOpBase<double> > M_explicit = ApB_transformer->createOutputOp();
271  ApB_transformer->transform( *M, M_explicit.ptr() );
272 
273  out << "\nM_explicit = " << *M_explicit;
274 
275  //
276  // E) Check the explicit operator
277  //
278 
279  LinearOpTester<double> M_explicit_tester;
280  M_explicit_tester.show_all_tests(true);;
281 
282  const bool result = M_explicit_tester.compare( *M, *M_explicit, Teuchos::inOutArg(out) );
283  if (!result) success = false;
284  }
285 
286  for(int scenario=0;scenario<4;scenario++) {
287  //
288  // C) Create implicit A+B operator
289  //
290 
291  const RCP<const Thyra::LinearOpBase<double> > M = buildAddOperator(scenario,B,A);
292 
293  //
294  // D) Do the transformation
295  //
296 
297  const RCP<EpetraExtAddTransformer> ApB_transformer = epetraExtAddTransformer();
298 
299  TEST_ASSERT(ApB_transformer != null);
300 
301  const RCP<LinearOpBase<double> > M_explicit = ApB_transformer->createOutputOp();
302  ApB_transformer->transform( *M, M_explicit.ptr() );
303 
304  out << "\nM_explicit = " << *M_explicit;
305 
306  //
307  // E) Check the explicit operator
308  //
309 
310  LinearOpTester<double> M_explicit_tester;
311  M_explicit_tester.show_all_tests(true);;
312 
313  const bool result = M_explicit_tester.compare( *M, *M_explicit, Teuchos::inOutArg(out) );
314  if (!result) success = false;
315  }
316 }
317 
318 TEUCHOS_UNIT_TEST( EpetraExtAddTransformer, basic_Add )
319 {
320 
321  //
322  // A) Read in problem matrices
323  //
324 
325  out << "\nReading linear system in Epetra format from the file \'"<<matrixFile<<"\' ...\n";
326  out << "\nReading linear system in Epetra format from the file \'"<<matrixFile2<<"\' ...\n";
327 
328 #ifdef HAVE_MPI
329  Epetra_MpiComm comm(MPI_COMM_WORLD);
330 #else
331  Epetra_SerialComm comm;
332 #endif
333  RCP<Epetra_CrsMatrix> epetra_A;
334  RCP<Epetra_CrsMatrix> epetra_B;
335  EpetraExt::readEpetraLinearSystem( matrixFile, comm, &epetra_A, NULL, NULL, NULL );
336  EpetraExt::readEpetraLinearSystem( matrixFile2, comm, &epetra_B, NULL, NULL, NULL );
337 
338  //
339  // B) Create the Thyra wrapped version
340  //
341  double scaleA=3.7;
342  double scaleB=-2.9;
343 
344  const RCP<const Thyra::LinearOpBase<double> > A = Thyra::scale<double>(scaleA,Thyra::epetraLinearOp(epetra_B));
345  const RCP<const Thyra::LinearOpBase<double> > B = Thyra::scale<double>(scaleB,Thyra::epetraLinearOp(epetra_B));
346 
347  out << "\nA = " << *A;
348  out << "\nB = " << *B;
349 
350  for(int scenario=0;scenario<4;scenario++) {
351  //
352  // C) Create implicit A+B operator
353  //
354 
355  const RCP<const Thyra::LinearOpBase<double> > M = buildAddOperator(scenario,A,B);
356 
357  //
358  // D) Do the transformation
359  //
360 
361  const RCP<EpetraExtAddTransformer> ApB_transformer = epetraExtAddTransformer();
362 
363  TEST_ASSERT(ApB_transformer != null);
364 
365  const RCP<LinearOpBase<double> > M_explicit = ApB_transformer->createOutputOp();
366  ApB_transformer->transform( *M, M_explicit.ptr() );
367 
368  out << "\nM_explicit = " << *M_explicit;
369 
370  //
371  // E) Check the explicit operator
372  //
373 
374  LinearOpTester<double> M_explicit_tester;
375  M_explicit_tester.show_all_tests(true);;
376 
377  const bool result = M_explicit_tester.compare( *M, *M_explicit, Teuchos::inOutArg(out) );
378  if (!result) success = false;
379  }
380 }
381 
382 TEUCHOS_UNIT_TEST( EpetraExtAddTransformer, mod_Add )
383 {
384 
385  //
386  // A) Read in problem matrices
387  //
388 
389  out << "\nReading linear system in Epetra format from the file \'"<<matrixFile<<"\' ...\n";
390  out << "\nReading linear system in Epetra format from the file \'"<<matrixFile2<<"\' ...\n";
391 
392 #ifdef HAVE_MPI
393  Epetra_MpiComm comm(MPI_COMM_WORLD);
394 #else
395  Epetra_SerialComm comm;
396 #endif
397  RCP<Epetra_CrsMatrix> epetra_A;
398  RCP<Epetra_CrsMatrix> epetra_B;
399  EpetraExt::readEpetraLinearSystem( matrixFile, comm, &epetra_A, NULL, NULL, NULL );
400  EpetraExt::readEpetraLinearSystem( matrixFile2, comm, &epetra_B, NULL, NULL, NULL );
401 
402  //
403  // B) Create the Thyra wrapped version
404  //
405  double scaleA=3.7;
406  double scaleB=-2.9;
407 
408  const RCP<const Thyra::LinearOpBase<double> > A = Thyra::scale<double>(scaleA,Thyra::epetraLinearOp(epetra_B));
409  const RCP<const Thyra::LinearOpBase<double> > B = Thyra::scale<double>(scaleB,Thyra::epetraLinearOp(epetra_B));
410 
411  out << "\nA = " << *A;
412  out << "\nB = " << *B;
413 
414  for(int scenario=0;scenario<4;scenario++) {
415  //
416  // C) Create implicit A+B operator
417  //
418 
419  const RCP<const Thyra::LinearOpBase<double> > M = buildAddOperator(scenario,A,B);
420 
421  //
422  // D) Do the transformation
423  //
424 
425  const RCP<EpetraExtAddTransformer> ApB_transformer = epetraExtAddTransformer();
426 
427  TEST_ASSERT(ApB_transformer != null);
428 
429  const RCP<LinearOpBase<double> > M_explicit = ApB_transformer->createOutputOp();
430  const RCP<LinearOpBase<double> > M_explicit_orig = M_explicit;
431  ApB_transformer->transform( *M, M_explicit.ptr() );
432 
433  // do some violence to one of the operators
434  double * view; int numEntries;
435  epetra_B->Scale(3.2);
436  TEUCHOS_ASSERT(epetra_B->ExtractMyRowView(3,numEntries,view)==0);
437  for(int i=0;i<numEntries;i++) view[i] += view[i]*double(i+1);
438 
439  // compute multiply again
440  ApB_transformer->transform( *M, M_explicit.ptr() );
441 
442  out << "\nM_explicit = " << *M_explicit;
443 
444  TEUCHOS_ASSERT(Thyra::get_Epetra_Operator(*M_explicit)
445  ==Thyra::get_Epetra_Operator(*M_explicit_orig));
446 
447  //
448  // E) Check the explicit operator
449  //
450 
451  LinearOpTester<double> M_explicit_tester;
452  M_explicit_tester.show_all_tests(true);;
453 
454  const bool result = M_explicit_tester.compare( *M, *M_explicit, Teuchos::inoutArg(out) );
455  if (!result) success = false;
456  }
457 }
458 
459 } // end namespace
Teuchos_VerboseObject.hpp
Thyra_EpetraExtAddTransformer.hpp
Thyra_EpetraLinearOp.hpp
TEST_ASSERT
#define TEST_ASSERT(v1)
Thyra::EpetraExtAddTransformer
Transformer subclass for adding Epetra/Thyra operators using EpetraExt::MatrixMatrix.
Definition: Thyra_EpetraExtAddTransformer.hpp:56
Teuchos::UnitTestRepository::getCLP
static CommandLineProcessor & getCLP()
Thyra::TEUCHOS_UNIT_TEST
TEUCHOS_UNIT_TEST(EpetraOperatorWrapper, basic)
Definition: EpetraOperatorWrapper_UnitTests.cpp:83
Thyra_get_Epetra_Operator.hpp
Teuchos::CommandLineProcessor::setOption
void setOption(const char option_true[], const char option_false[], bool *option_val, const char documentation[]=NULL)
Thyra_EpetraThyraWrappers.hpp
Teuchos_UnitTestHarness.hpp
Teuchos_StandardCatchMacros.hpp
Teuchos_XMLParameterListHelpers.hpp
Teuchos_GlobalMPISession.hpp
Thyra::get_Epetra_Operator
Teuchos::RCP< Epetra_Operator > get_Epetra_Operator(LinearOpBase< double > &op)
Full specialization for Scalar=double.
Definition: Thyra_get_Epetra_Operator.cpp:50
Teuchos_CommandLineProcessor.hpp
TEUCHOS_ASSERT
#define TEUCHOS_ASSERT(assertion_test)
Generated by   doxygen 1.8.5