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 // @HEADER
2 // *****************************************************************************
3 // Thyra: Interfaces and Support for Abstract Numerical Algorithms
4 //
5 // Copyright 2004 NTESS and the Thyra contributors.
6 // SPDX-License-Identifier: BSD-3-Clause
7 // *****************************************************************************
8 // @HEADER
9 
10 #include "Thyra_EpetraExtAddTransformer.hpp"
11 #include "Thyra_DefaultScaledAdjointLinearOp.hpp"
12 #include "Thyra_DefaultMultipliedLinearOp.hpp"
13 #include "Thyra_DefaultAddedLinearOp.hpp"
14 #include "Thyra_DefaultDiagonalLinearOp.hpp"
15 #include "Thyra_VectorStdOps.hpp"
16 #include "Thyra_TestingTools.hpp"
17 #include "Thyra_LinearOpTester.hpp"
18 #include "Thyra_EpetraThyraWrappers.hpp"
19 #include "Thyra_EpetraLinearOp.hpp"
20 #include "Thyra_get_Epetra_Operator.hpp"
21 #include "Thyra_DefaultIdentityLinearOp.hpp"
22 #include "EpetraExt_readEpetraLinearSystem.h"
23 #include "Teuchos_GlobalMPISession.hpp"
24 #include "Teuchos_VerboseObject.hpp"
25 #include "Teuchos_XMLParameterListHelpers.hpp"
26 #include "Teuchos_CommandLineProcessor.hpp"
27 #include "Teuchos_StandardCatchMacros.hpp"
28 
29 #ifdef HAVE_MPI
30 # include "Epetra_MpiComm.h"
31 #else
32 # include "Epetra_SerialComm.h"
33 #endif
34 
35 #include "Teuchos_UnitTestHarness.hpp"
36 
37 
38 namespace {
39 
40 
41 using Teuchos::null;
42 using Teuchos::RCP;
43 using Thyra::EpetraExtAddTransformer;
44 using Thyra::epetraExtAddTransformer;
45 using Thyra::VectorBase;
46 using Thyra::LinearOpBase;
47 using Thyra::createMember;
48 using Thyra::LinearOpTester;
49 using Thyra::adjoint;
50 using Thyra::multiply;
51 using Thyra::diagonal;
52 
53 
54 std::string matrixFile = "";
55 std::string matrixFile2 = "";
56 
57 
58 TEUCHOS_STATIC_SETUP()
59 {
60  Teuchos::UnitTestRepository::getCLP().setOption(
61  "matrix-file", &matrixFile,
62  "Defines the Epetra_CrsMatrix to read in." );
63  Teuchos::UnitTestRepository::getCLP().setOption(
64  "matrix-file-2", &matrixFile2,
65  "Defines the Epetra_CrsMatrix to read in." );
66 }
67 
68 const Teuchos::RCP<const Thyra::LinearOpBase<double> >
69 buildAddOperator(int scenario,const Teuchos::RCP<const Thyra::LinearOpBase<double> > & A,
70  const Teuchos::RCP<const Thyra::LinearOpBase<double> > & B)
71 {
72  // build operators for the various addition/adjoint scenarios
73  RCP<const Thyra::LinearOpBase<double> > M;
74 
75  switch(scenario) {
76  case 0:
77  M = Thyra::add(A,B,"A+B");
78  break;
79  case 1:
80  M = Thyra::add(A,B,"A+adj(B)");
81  break;
82  case 2:
83  M = Thyra::add(A,B,"adj(A)+B");
84  break;
85  case 3:
86  M = Thyra::add(A,B,"adb(A)+adb(B)");
87  break;
88  default:
89  TEUCHOS_ASSERT(false);
90  }
91 
92  return M;
93 }
94 
95 TEUCHOS_UNIT_TEST( EpetraExtAddTransformer, diag_mat_Add )
96 {
97 
98  //
99  // A) Read in problem matrices
100  //
101 
102  out << "\nReading linear system in Epetra format from the file \'"<<matrixFile<<"\' ...\n";
103 
104 #ifdef HAVE_MPI
105  Epetra_MpiComm comm(MPI_COMM_WORLD);
106 #else
107  Epetra_SerialComm comm;
108 #endif
109  RCP<Epetra_CrsMatrix> crsMat;
110  EpetraExt::readEpetraLinearSystem( matrixFile, comm, &crsMat, NULL, NULL, NULL );
111 
112  //
113  // B) Create the Thyra wrapped version
114  //
115  double scaleMat=3.7;
116  //double scaleDiag=-2.9;
117 
118 
119  const RCP<const Thyra::LinearOpBase<double> > A = Thyra::scale<double>(scaleMat,Thyra::epetraLinearOp(crsMat));
120  RCP<VectorBase<double> > d = createMember(A->domain(), "d");
121  V_S( d.ptr(), 3.0 ); // ToDo: Set ton != 1.0 and generalize
122  const RCP<const Thyra::LinearOpBase<double> > B = diagonal(d);
123 
124  out << "\nA = " << *A;
125  out << "\nB = " << *B;
126 
127  for(int scenario=0;scenario<4;scenario++) {
128  //
129  // C) Create implicit A+B operator
130  //
131 
132  const RCP<const Thyra::LinearOpBase<double> > M = buildAddOperator(scenario,A,B);
133 
134  //
135  // D) Do the transformation
136  //
137 
138  const RCP<EpetraExtAddTransformer> ApB_transformer = epetraExtAddTransformer();
139 
140  TEST_ASSERT(ApB_transformer != null);
141 
142  const RCP<LinearOpBase<double> > M_explicit = ApB_transformer->createOutputOp();
143  ApB_transformer->transform( *M, M_explicit.ptr() );
144 
145  out << "\nM_explicit = " << *M_explicit;
146 
147  //
148  // E) Check the explicit operator
149  //
150 
151  LinearOpTester<double> M_explicit_tester;
152  M_explicit_tester.show_all_tests(true);;
153 
154  const bool result = M_explicit_tester.compare( *M, *M_explicit, Teuchos::inOutArg(out) );
155  if (!result) success = false;
156  }
157 
158  for(int scenario=0;scenario<4;scenario++) {
159  //
160  // C) Create implicit A+B operator
161  //
162 
163  const RCP<const Thyra::LinearOpBase<double> > M = buildAddOperator(scenario,B,A);
164 
165  //
166  // D) Do the transformation
167  //
168 
169  const RCP<EpetraExtAddTransformer> ApB_transformer = epetraExtAddTransformer();
170 
171  TEST_ASSERT(ApB_transformer != null);
172 
173  const RCP<LinearOpBase<double> > M_explicit = ApB_transformer->createOutputOp();
174  ApB_transformer->transform( *M, M_explicit.ptr() );
175 
176  out << "\nM_explicit = " << *M_explicit;
177 
178  //
179  // E) Check the explicit operator
180  //
181 
182  LinearOpTester<double> M_explicit_tester;
183  M_explicit_tester.show_all_tests(true);;
184 
185  const bool result = M_explicit_tester.compare( *M, *M_explicit, Teuchos::inOutArg(out) );
186  if (!result) success = false;
187  }
188 }
189 
190 TEUCHOS_UNIT_TEST( EpetraExtAddTransformer, id_mat_Add )
191 {
192 
193  //
194  // A) Read in problem matrices
195  //
196 
197  out << "\nReading linear system in Epetra format from the file \'"<<matrixFile<<"\' ...\n";
198 
199 #ifdef HAVE_MPI
200  Epetra_MpiComm comm(MPI_COMM_WORLD);
201 #else
202  Epetra_SerialComm comm;
203 #endif
204  RCP<Epetra_CrsMatrix> crsMat;
205  EpetraExt::readEpetraLinearSystem( matrixFile, comm, &crsMat, NULL, NULL, NULL );
206 
207  //
208  // B) Create the Thyra wrapped version
209  //
210  double scaleMat=3.7;
211  //double scaleDiag=-2.9;
212 
213 
214  const RCP<const Thyra::LinearOpBase<double> > A = Thyra::scale<double>(scaleMat,Thyra::epetraLinearOp(crsMat));
215  const RCP<const Thyra::LinearOpBase<double> > B = identity(A->domain(),"id");
216 
217  out << "\nA = " << *A;
218  out << "\nB = " << *B;
219 
220  for(int scenario=0;scenario<4;scenario++) {
221  //
222  // C) Create implicit A+B operator
223  //
224 
225  const RCP<const Thyra::LinearOpBase<double> > M = buildAddOperator(scenario,A,B);
226 
227  //
228  // D) Do the transformation
229  //
230 
231  const RCP<EpetraExtAddTransformer> ApB_transformer = epetraExtAddTransformer();
232 
233  TEST_ASSERT(ApB_transformer != null);
234 
235  const RCP<LinearOpBase<double> > M_explicit = ApB_transformer->createOutputOp();
236  ApB_transformer->transform( *M, M_explicit.ptr() );
237 
238  out << "\nM_explicit = " << *M_explicit;
239 
240  //
241  // E) Check the explicit operator
242  //
243 
244  LinearOpTester<double> M_explicit_tester;
245  M_explicit_tester.show_all_tests(true);;
246 
247  const bool result = M_explicit_tester.compare( *M, *M_explicit, Teuchos::inOutArg(out) );
248  if (!result) success = false;
249  }
250 
251  for(int scenario=0;scenario<4;scenario++) {
252  //
253  // C) Create implicit A+B operator
254  //
255 
256  const RCP<const Thyra::LinearOpBase<double> > M = buildAddOperator(scenario,B,A);
257 
258  //
259  // D) Do the transformation
260  //
261 
262  const RCP<EpetraExtAddTransformer> ApB_transformer = epetraExtAddTransformer();
263 
264  TEST_ASSERT(ApB_transformer != null);
265 
266  const RCP<LinearOpBase<double> > M_explicit = ApB_transformer->createOutputOp();
267  ApB_transformer->transform( *M, M_explicit.ptr() );
268 
269  out << "\nM_explicit = " << *M_explicit;
270 
271  //
272  // E) Check the explicit operator
273  //
274 
275  LinearOpTester<double> M_explicit_tester;
276  M_explicit_tester.show_all_tests(true);;
277 
278  const bool result = M_explicit_tester.compare( *M, *M_explicit, Teuchos::inOutArg(out) );
279  if (!result) success = false;
280  }
281 }
282 
283 TEUCHOS_UNIT_TEST( EpetraExtAddTransformer, basic_Add )
284 {
285 
286  //
287  // A) Read in problem matrices
288  //
289 
290  out << "\nReading linear system in Epetra format from the file \'"<<matrixFile<<"\' ...\n";
291  out << "\nReading linear system in Epetra format from the file \'"<<matrixFile2<<"\' ...\n";
292 
293 #ifdef HAVE_MPI
294  Epetra_MpiComm comm(MPI_COMM_WORLD);
295 #else
296  Epetra_SerialComm comm;
297 #endif
298  RCP<Epetra_CrsMatrix> epetra_A;
299  RCP<Epetra_CrsMatrix> epetra_B;
300  EpetraExt::readEpetraLinearSystem( matrixFile, comm, &epetra_A, NULL, NULL, NULL );
301  EpetraExt::readEpetraLinearSystem( matrixFile2, comm, &epetra_B, NULL, NULL, NULL );
302 
303  //
304  // B) Create the Thyra wrapped version
305  //
306  double scaleA=3.7;
307  double scaleB=-2.9;
308 
309  const RCP<const Thyra::LinearOpBase<double> > A = Thyra::scale<double>(scaleA,Thyra::epetraLinearOp(epetra_B));
310  const RCP<const Thyra::LinearOpBase<double> > B = Thyra::scale<double>(scaleB,Thyra::epetraLinearOp(epetra_B));
311 
312  out << "\nA = " << *A;
313  out << "\nB = " << *B;
314 
315  for(int scenario=0;scenario<4;scenario++) {
316  //
317  // C) Create implicit A+B operator
318  //
319 
320  const RCP<const Thyra::LinearOpBase<double> > M = buildAddOperator(scenario,A,B);
321 
322  //
323  // D) Do the transformation
324  //
325 
326  const RCP<EpetraExtAddTransformer> ApB_transformer = epetraExtAddTransformer();
327 
328  TEST_ASSERT(ApB_transformer != null);
329 
330  const RCP<LinearOpBase<double> > M_explicit = ApB_transformer->createOutputOp();
331  ApB_transformer->transform( *M, M_explicit.ptr() );
332 
333  out << "\nM_explicit = " << *M_explicit;
334 
335  //
336  // E) Check the explicit operator
337  //
338 
339  LinearOpTester<double> M_explicit_tester;
340  M_explicit_tester.show_all_tests(true);;
341 
342  const bool result = M_explicit_tester.compare( *M, *M_explicit, Teuchos::inOutArg(out) );
343  if (!result) success = false;
344  }
345 }
346 
347 TEUCHOS_UNIT_TEST( EpetraExtAddTransformer, mod_Add )
348 {
349 
350  //
351  // A) Read in problem matrices
352  //
353 
354  out << "\nReading linear system in Epetra format from the file \'"<<matrixFile<<"\' ...\n";
355  out << "\nReading linear system in Epetra format from the file \'"<<matrixFile2<<"\' ...\n";
356 
357 #ifdef HAVE_MPI
358  Epetra_MpiComm comm(MPI_COMM_WORLD);
359 #else
360  Epetra_SerialComm comm;
361 #endif
362  RCP<Epetra_CrsMatrix> epetra_A;
363  RCP<Epetra_CrsMatrix> epetra_B;
364  EpetraExt::readEpetraLinearSystem( matrixFile, comm, &epetra_A, NULL, NULL, NULL );
365  EpetraExt::readEpetraLinearSystem( matrixFile2, comm, &epetra_B, NULL, NULL, NULL );
366 
367  //
368  // B) Create the Thyra wrapped version
369  //
370  double scaleA=3.7;
371  double scaleB=-2.9;
372 
373  const RCP<const Thyra::LinearOpBase<double> > A = Thyra::scale<double>(scaleA,Thyra::epetraLinearOp(epetra_B));
374  const RCP<const Thyra::LinearOpBase<double> > B = Thyra::scale<double>(scaleB,Thyra::epetraLinearOp(epetra_B));
375 
376  out << "\nA = " << *A;
377  out << "\nB = " << *B;
378 
379  for(int scenario=0;scenario<4;scenario++) {
380  //
381  // C) Create implicit A+B operator
382  //
383 
384  const RCP<const Thyra::LinearOpBase<double> > M = buildAddOperator(scenario,A,B);
385 
386  //
387  // D) Do the transformation
388  //
389 
390  const RCP<EpetraExtAddTransformer> ApB_transformer = epetraExtAddTransformer();
391 
392  TEST_ASSERT(ApB_transformer != null);
393 
394  const RCP<LinearOpBase<double> > M_explicit = ApB_transformer->createOutputOp();
395  const RCP<LinearOpBase<double> > M_explicit_orig = M_explicit;
396  ApB_transformer->transform( *M, M_explicit.ptr() );
397 
398  // do some violence to one of the operators
399  double * view; int numEntries;
400  epetra_B->Scale(3.2);
401  TEUCHOS_ASSERT(epetra_B->ExtractMyRowView(3,numEntries,view)==0);
402  for(int i=0;i<numEntries;i++) view[i] += view[i]*double(i+1);
403 
404  // compute multiply again
405  ApB_transformer->transform( *M, M_explicit.ptr() );
406 
407  out << "\nM_explicit = " << *M_explicit;
408 
409  TEUCHOS_ASSERT(Thyra::get_Epetra_Operator(*M_explicit)
410  ==Thyra::get_Epetra_Operator(*M_explicit_orig));
411 
412  //
413  // E) Check the explicit operator
414  //
415 
416  LinearOpTester<double> M_explicit_tester;
417  M_explicit_tester.show_all_tests(true);;
418 
419  const bool result = M_explicit_tester.compare( *M, *M_explicit, Teuchos::inoutArg(out) );
420  if (!result) success = false;
421  }
422 }
423 
424 } // 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:24
Teuchos::UnitTestRepository::getCLP
static CommandLineProcessor & getCLP()
Thyra::TEUCHOS_UNIT_TEST
TEUCHOS_UNIT_TEST(EpetraOperatorWrapper, basic)
Definition: EpetraOperatorWrapper_UnitTests.cpp:49
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:18
Teuchos_CommandLineProcessor.hpp
TEUCHOS_ASSERT
#define TEUCHOS_ASSERT(assertion_test)
Generated on Wed Dec 25 2024 09:21:43 for Thyra Package Browser (Single Doxygen Collection) by   doxygen 1.8.5