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