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