Thyra Package Browser (Single Doxygen Collection)  Version of the Day
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 
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"
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;
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 
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  break;
145  }
146 
147 
148  out << "\nM = " << *M;
149 
150  return M;
151 }
152 
153 // helper function to excercise all different versions of B*D*G
154 const Teuchos::RCP<const Thyra::LinearOpBase<double> >
155 buildBGOperator(int scenario,const Teuchos::RCP<const Thyra::LinearOpBase<double> > & B,
156  const Teuchos::RCP<const Thyra::LinearOpBase<double> > & G,
157  std::ostream & out)
158 {
159  // Create the implicit operator
160  double scalea=10.0;
161  double scaleb=-7.0;
162  RCP<const LinearOpBase<double> > M;
163 
164  // create an operator based on requested scenario
165  // (these are the same as in EpetraExt)
166  switch(scenario) {
167  case 1:
168  M = multiply( scale(scalea,B), scale(scaleb,G), "M" );
169  out << " Testing B*G" << std::endl;
170  break;
171  case 2:
172  M = multiply( scale(scalea,B), adjoint(G), "M" );
173  out << " Testing B*adj(G)" << std::endl;
174  break;
175  case 3:
176  M = multiply( adjoint(B), scale(scaleb,G), "M" );
177  out << " Testing adj(B)*G" << std::endl;
178  break;
179  case 4:
180  M = multiply( adjoint(B), adjoint(G), "M" );
181  out << " Testing adj(B)*adj(G)" << std::endl;
182  break;
183  default:
184  TEUCHOS_ASSERT(false);
185  break;
186  }
187 
188 
189  out << "\nM = " << *M;
190 
191  return M;
192 }
193 
194 const Teuchos::RCP<const Thyra::LinearOpBase<double> >
195 buildBDBOperator(int scenario,const Teuchos::RCP<const Thyra::LinearOpBase<double> > & B,
196  const double vecScale, std::ostream & out)
197 {
198  return buildBDGOperator(scenario,B,B,vecScale,out);
199 }
200 
201 
202 
203 TEUCHOS_UNIT_TEST( EpetraExtDiagScaledMatProdTransformer, basic_BDB )
204 {
205 
206  //
207  // A) Read in problem matrices
208  //
209 
210  out << "\nReading linear system in Epetra format from the file \'"<<matrixFile<<"\' ...\n";
211 
212 #ifdef HAVE_MPI
213  Epetra_MpiComm comm(MPI_COMM_WORLD);
214 #else
215  Epetra_SerialComm comm;
216 #endif
217  RCP<Epetra_CrsMatrix> epetra_B;
218  EpetraExt::readEpetraLinearSystem( matrixFile, comm, &epetra_B, NULL, NULL, NULL );
219 
220  //
221  // B) Create the Thyra wrapped version
222  //
223 
224  const RCP<const Thyra::LinearOpBase<double> > B = Thyra::epetraLinearOp(epetra_B);
225 
226  //
227  // C) Create implicit B*D*B operator
228  //
229 
230  // build scenario=1 -> B*D*B, scenario=2 -> B*D*B',
231  // scenario=3 -> B'*D*B, scenario=4 -> B'*D*B'
232  //int scenario = 3;
233  for(int scenario=1;scenario<=5;scenario++) {
234  const RCP<const Thyra::LinearOpBase<double> > M = buildBDBOperator(scenario,B,4.5,out);
235 
236  //
237  // D) Do the transformation
238  //
239 
240  const RCP<EpetraExtDiagScaledMatProdTransformer> BtDB_transformer =
242 
243  TEST_ASSERT(BtDB_transformer != null);
244 
245  const RCP<LinearOpBase<double> > M_explicit = BtDB_transformer->createOutputOp();
246 
247  BtDB_transformer->transform( *M, M_explicit.ptr() );
248 
249  out << "\nM_explicit = " << *M_explicit;
250 
251  //
252  // E) Check the explicit operator
253  //
254 
255  LinearOpTester<double> M_explicit_tester;
256  M_explicit_tester.show_all_tests(true);;
257 
258  const bool result = M_explicit_tester.compare( *M, *M_explicit, Teuchos::inoutArg(out) );
259  if (!result) success = false;
260  }
261 
262 }
263 
264 TEUCHOS_UNIT_TEST( EpetraExtDiagScaledMatProdTransformer, basic_BDG_GDB)
265 {
266 
267  //
268  // A) Read in problem matrices
269  //
270 
271  out << "\nReading linear system in Epetra format from the file \'"<<matrixFile<<"\'";
272  out << " and from the file \'"<<matrixFile2<<"\' ...\n";
273 
274 #ifdef HAVE_MPI
275  Epetra_MpiComm comm(MPI_COMM_WORLD);
276 #else
277  Epetra_SerialComm comm;
278 #endif
279  RCP<Epetra_CrsMatrix> epetra_B;
280  RCP<Epetra_CrsMatrix> epetra_G;
281  EpetraExt::readEpetraLinearSystem( matrixFile, comm, &epetra_B, NULL, NULL, NULL );
282  EpetraExt::readEpetraLinearSystem( matrixFile, comm, &epetra_G, NULL, NULL, NULL );
283 
284  //
285  // B) Create the Thyra wrapped version
286  //
287 
288  const RCP<const Thyra::LinearOpBase<double> > B = Thyra::epetraLinearOp(epetra_B);
289  const RCP<const Thyra::LinearOpBase<double> > G = Thyra::epetraLinearOp(epetra_G);
290 
291  //
292  // C) Create implicit B*D*B operator
293  //
294 
295  // build scenario=1 -> B*D*B, scenario=2-> B*D*B', scenario=3 -> B'*D*B
296  for(int scenario=1;scenario<=3;scenario++) {
297  RCP<const Thyra::LinearOpBase<double> > M;
298  if(scenario==1 || scenario==3)
299  M = buildBDGOperator(scenario,B,G,4.5,out);
300  else
301  M = buildBDGOperator(scenario,G,B,4.5,out);
302 
303  //
304  // D) Do the transformation
305  //
306 
307  const RCP<EpetraExtDiagScaledMatProdTransformer> BtDB_transformer =
309 
310  TEST_ASSERT(BtDB_transformer != null);
311 
312  const RCP<LinearOpBase<double> > M_explicit = BtDB_transformer->createOutputOp();
313 
314  BtDB_transformer->transform( *M, M_explicit.ptr() );
315 
316  out << "\nM_explicit = " << *M_explicit;
317 
318  //
319  // E) Check the explicit operator
320  //
321 
322  LinearOpTester<double> M_explicit_tester;
323  M_explicit_tester.show_all_tests(true);;
324 
325  const bool result = M_explicit_tester.compare( *M, *M_explicit, Teuchos::inoutArg(out) );
326  if (!result) success = false;
327  }
328 
329 }
330 
331 TEUCHOS_UNIT_TEST( EpetraExtDiagScaledMatProdTransformer, basic_GDG )
332 {
333 
334  //
335  // A) Read in problem matrices
336  //
337 
338  out << "\nReading linear system in Epetra format from the file \'"<<matrixFile2<<"\' ...\n";
339 
340 #ifdef HAVE_MPI
341  Epetra_MpiComm comm(MPI_COMM_WORLD);
342 #else
343  Epetra_SerialComm comm;
344 #endif
345  RCP<Epetra_CrsMatrix> epetra_G;
346  EpetraExt::readEpetraLinearSystem( matrixFile, comm, &epetra_G, NULL, NULL, NULL );
347 
348  //
349  // B) Create the Thyra wrapped version
350  //
351 
352  const RCP<const Thyra::LinearOpBase<double> > G = Thyra::epetraLinearOp(epetra_G);
353 
354  //
355  // C) Create implicit B*D*B operator
356  //
357 
358  // build scenario=1 -> B*D*B, scenario=3 -> B'*D*B
359  int scenes[] = {1,4};
360  for(int i=0;i<2;i++) {
361  int scenario = scenes[i];
362  const RCP<const Thyra::LinearOpBase<double> > M = buildBDGOperator(scenario,G,G,4.5,out);
363 
364  //
365  // D) Do the transformation
366  //
367 
368  const RCP<EpetraExtDiagScaledMatProdTransformer> BtDB_transformer =
370 
371  TEST_ASSERT(BtDB_transformer != null);
372 
373  const RCP<LinearOpBase<double> > M_explicit = BtDB_transformer->createOutputOp();
374 
375  BtDB_transformer->transform( *M, M_explicit.ptr() );
376 
377  out << "\nM_explicit = " << *M_explicit;
378 
379  //
380  // E) Check the explicit operator
381  //
382 
383  LinearOpTester<double> M_explicit_tester;
384  M_explicit_tester.show_all_tests(true);;
385 
386  const bool result = M_explicit_tester.compare( *M, *M_explicit, Teuchos::inoutArg(out) );
387  if (!result) success = false;
388  }
389 
390 }
391 
392 TEUCHOS_UNIT_TEST( EpetraExtDiagScaledMatProdTransformer, basic_BG_GB_GG)
393 {
394 
395  //
396  // A) Read in problem matrices
397  //
398 
399  out << "\nReading linear system in Epetra format from the file \'"<<matrixFile<<"\'";
400  out << " and from the file \'"<<matrixFile2<<"\' ...\n";
401 
402 #ifdef HAVE_MPI
403  Epetra_MpiComm comm(MPI_COMM_WORLD);
404 #else
405  Epetra_SerialComm comm;
406 #endif
407  RCP<Epetra_CrsMatrix> epetra_B;
408  RCP<Epetra_CrsMatrix> epetra_G;
409  EpetraExt::readEpetraLinearSystem( matrixFile, comm, &epetra_B, NULL, NULL, NULL );
410  EpetraExt::readEpetraLinearSystem( matrixFile, comm, &epetra_G, NULL, NULL, NULL );
411 
412  //
413  // B) Create the Thyra wrapped version
414  //
415 
416  const RCP<const Thyra::LinearOpBase<double> > B = Thyra::epetraLinearOp(epetra_B);
417  const RCP<const Thyra::LinearOpBase<double> > G = Thyra::epetraLinearOp(epetra_G);
418 
419  //
420  // C) Create implicit B*D*B operator
421  //
422 
423  // build scenario=1 -> B*D*B, scenario=2-> B*D*B', scenario=3 -> B'*D*B
424  for(int scenario=1;scenario<=4;scenario++) {
425  RCP<const Thyra::LinearOpBase<double> > M;
426  if(scenario==1 || scenario==3)
427  M = buildBGOperator(scenario,B,G,out);
428  else if(scenario==2)
429  M = buildBGOperator(scenario,G,B,out);
430  else
431  M = buildBGOperator(scenario,G,G,out);
432 
433  //
434  // D) Do the transformation
435  //
436 
437  const RCP<EpetraExtDiagScaledMatProdTransformer> BtDB_transformer =
439 
440  TEST_ASSERT(BtDB_transformer != null);
441 
442  const RCP<LinearOpBase<double> > M_explicit = BtDB_transformer->createOutputOp();
443 
444  BtDB_transformer->transform( *M, M_explicit.ptr() );
445 
446  out << "\nM_explicit = " << *M_explicit;
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 
461 } // namespace
TEUCHOS_UNIT_TEST(EpetraOperatorWrapper, basic)
Transformer subclass for diagonally scaling and multiplying Epetra/Thyra operators.
RCP< EpetraExtDiagScaledMatProdTransformer > epetraExtDiagScaledMatProdTransformer()
Nonmember constructor.