Thyra Package Browser (Single Doxygen Collection)  Version of the Day
EpetraExtDiagScalingTransformer_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"
55 #include "EpetraExt_readEpetraLinearSystem.h"
56 #include "Teuchos_GlobalMPISession.hpp"
57 #include "Teuchos_VerboseObject.hpp"
58 #include "Teuchos_XMLParameterListHelpers.hpp"
59 #include "Teuchos_CommandLineProcessor.hpp"
60 #include "Teuchos_StandardCatchMacros.hpp"
61 #ifdef _MSC_VER
62 // This is required for operator not to be defined
63 #include <iso646.h>
64 #endif
65 
66 #ifdef HAVE_MPI
67 # include "Epetra_MpiComm.h"
68 #else
69 # include "Epetra_SerialComm.h"
70 #endif
71 
72 #include "Teuchos_UnitTestHarness.hpp"
73 
74 
75 namespace {
76 
77 
78 using Teuchos::null;
79 using Teuchos::RCP;
81 using Thyra::epetraExtDiagScalingTransformer;
82 using Thyra::VectorBase;
83 using Thyra::LinearOpBase;
84 using Thyra::createMember;
85 using Thyra::LinearOpTester;
86 using Thyra::adjoint;
87 using Thyra::multiply;
88 using Thyra::diagonal;
89 
90 
91 std::string matrixFile = "";
92 
93 
95 {
96  Teuchos::UnitTestRepository::getCLP().setOption(
97  "matrix-file", &matrixFile,
98  "Defines the Epetra_CrsMatrix to read in." );
99 }
100 
101 // helper function to excercise all different versions of B*D*G
102 const Teuchos::RCP<const Thyra::LinearOpBase<double> >
103 buildADOperator(int scenario,const Teuchos::RCP<const Thyra::LinearOpBase<double> > & A,
104  const double vecScale, std::ostream & out)
105 {
106  // Create the implicit operator
107  double scalea=-7.0;
108  double scaled=10.0;
109 
110  RCP<const LinearOpBase<double> > M;
111  RCP<VectorBase<double> > d;
112  if(scenario<=2)
113  d = createMember(A->domain(), "d");
114  else
115  d = createMember(A->range(), "d");
116  V_S( d.ptr(), vecScale ); // ToDo: Set ton != 1.0 and generalize
117 
118  // create an operator based on requested scenario
119  // (these are the same as in EpetraExt)
120  switch(scenario) {
121  case 1:
122  M = multiply( scale(scalea,A), scale(scaled,diagonal(d)), "M" );
123  out << " Testing A*D" << std::endl;
124  break;
125  case 2:
126  M = multiply( scale(scaled,diagonal(d)), scale(scalea,A), "M" );
127  out << " Testing D*A" << std::endl;
128  break;
129  case 3:
130  M = multiply( A, scale(scaled,diagonal(d)), "M" );
131  out << " Testing adj(A)*D" << std::endl;
132  break;
133  case 4:
134  M = multiply( scale(scaled,diagonal(d)), A, "M" );
135  out << " Testing D*adj(A)" << std::endl;
136  break;
137  default:
138  TEUCHOS_ASSERT(false);
139  break;
140  }
141 
142 
143  out << "\nM = " << *M;
144 
145  return M;
146 }
147 
148 TEUCHOS_UNIT_TEST( EpetraExtDiagScalingTransformer, isCompatible)
149 {
150  out << "\nReading linear system in Epetra format from the file \'"<<matrixFile<<"\'\n";
151 
152 #ifdef HAVE_MPI
153  Epetra_MpiComm comm(MPI_COMM_WORLD);
154 #else
155  Epetra_SerialComm comm;
156 #endif
157 
158  RCP<Epetra_CrsMatrix> epetra_A;
159  RCP<Epetra_CrsMatrix> epetra_B;
160  EpetraExt::readEpetraLinearSystem( matrixFile, comm, &epetra_A, NULL, NULL, NULL );
161  EpetraExt::readEpetraLinearSystem( matrixFile, comm, &epetra_B, NULL, NULL, NULL );
162  const RCP<const Thyra::LinearOpBase<double> > A = Thyra::epetraLinearOp(epetra_A);
163  const RCP<const Thyra::LinearOpBase<double> > B = Thyra::epetraLinearOp(epetra_B);
164 
165  RCP<VectorBase<double> > d = createMember(B->domain(), "d");
166  V_S( d.ptr(), 3.0 ); // ToDo: Set ton != 1.0 and generalize
167 
168  const RCP<EpetraExtDiagScalingTransformer> BD_transformer =
170 
171  TEST_ASSERT(BD_transformer != null);
172 
173  RCP<const LinearOpBase<double> > M;
174 
175  M = multiply(A,B);
176  TEST_ASSERT(not BD_transformer->isCompatible(*M));
177 
178  M = multiply(A,scale(3.9,diagonal(d)));
179  TEST_ASSERT(BD_transformer->isCompatible(*M));
180 
181  M = multiply(scale(3.0,diagonal(d)),scale(9.0,B));
182  TEST_ASSERT(BD_transformer->isCompatible(*M));
183 
184  M = multiply(B,scale(3.0,diagonal(d)),scale(9.0,B));
185  TEST_ASSERT(not BD_transformer->isCompatible(*M));
186 }
187 
188 TEUCHOS_UNIT_TEST( EpetraExtDiagScalingTransformer, reapply_scaling)
189 {
190  //
191  // A) Read in problem matrices
192  //
193 
194  out << "\nReading linear system in Epetra format from the file \'"<<matrixFile<<"\'\n";
195 
196 #ifdef HAVE_MPI
197  Epetra_MpiComm comm(MPI_COMM_WORLD);
198 #else
199  Epetra_SerialComm comm;
200 #endif
201  RCP<Epetra_CrsMatrix> epetra_A;
202  EpetraExt::readEpetraLinearSystem( matrixFile, comm, &epetra_A, NULL, NULL, NULL );
203 
204  //
205  // B) Create the Thyra wrapped version
206  //
207 
208  const RCP<const Thyra::LinearOpBase<double> > A = Thyra::epetraLinearOp(epetra_A);
209 
210  RCP<const Thyra::LinearOpBase<double> > M;
211 
212  // test scale on right
213  for(int scenario=1;scenario<=2;scenario++) {
214  out << "**********************************************************" << std::endl;
215  out << "**************** Starting Scenario = " << scenario << std::endl;
216  out << "**********************************************************" << std::endl;
217 
218  //
219  // D) Check compatibility
220  //
221  const RCP<EpetraExtDiagScalingTransformer> BD_transformer =
223 
224  TEST_ASSERT(BD_transformer != null);
225 
226  //
227  // E) Do the transformation (twice)
228  //
229 
230  const RCP<LinearOpBase<double> > M_explicit = BD_transformer->createOutputOp();
231  TEST_ASSERT(M_explicit != null);
232 
233  // do trans
234  M = buildADOperator(scenario,A,4.5,out);
235  BD_transformer->transform( *M, M_explicit.ptr() );
236  const RCP<const Epetra_Operator> eOp_explicit0 = Thyra::get_Epetra_Operator(*M_explicit);
237 
238  M = buildADOperator(scenario,A,7.5,out);
239  BD_transformer->transform( *M, M_explicit.ptr() );
240  const RCP<const Epetra_Operator> eOp_explicit1 = Thyra::get_Epetra_Operator(*M_explicit);
241 
242  // Epetra operator should not change!
243  TEST_EQUALITY(eOp_explicit0,eOp_explicit1);
244 
245  out << "\nM_explicit = " << *M_explicit;
246 
247  //
248  // F) 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 TEUCHOS_UNIT_TEST( EpetraExtDiagScalingTransformer, basic_BDG_GDB)
260 {
261 
262  //
263  // A) Read in problem matrices
264  //
265 
266  out << "\nReading linear system in Epetra format from the file \'"<<matrixFile<<"\'\n";
267 
268 #ifdef HAVE_MPI
269  Epetra_MpiComm comm(MPI_COMM_WORLD);
270 #else
271  Epetra_SerialComm comm;
272 #endif
273  RCP<Epetra_CrsMatrix> epetra_A;
274  EpetraExt::readEpetraLinearSystem( matrixFile, comm, &epetra_A, NULL, NULL, NULL );
275 
276  //
277  // B) Create the Thyra wrapped version
278  //
279 
280  const RCP<const Thyra::LinearOpBase<double> > A = Thyra::epetraLinearOp(epetra_A);
281 
282  //
283  // C) Create implicit B*D*B operator
284  //
285 
286  // build scenario=1 -> B*D*B, scenario=2-> B*D*B', scenario=3 -> B'*D*B
287  for(int scenario=1;scenario<=4;scenario++) {
288  out << "**********************************************************" << std::endl;
289  out << "**************** Starting Scenario = " << scenario << std::endl;
290  out << "**********************************************************" << std::endl;
291 
292  RCP<const Thyra::LinearOpBase<double> > M = buildADOperator(scenario,A,4.5,out);
293 
294  //
295  // D) Check compatibility
296  //
297  const RCP<EpetraExtDiagScalingTransformer> BD_transformer =
299 
300  TEST_ASSERT(BD_transformer != null);
301 
302  BD_transformer->isCompatible(*M);
303 
304  //
305  // E) Do the transformation
306  //
307 
308  const RCP<LinearOpBase<double> > M_explicit = BD_transformer->createOutputOp();
309  BD_transformer->transform( *M, M_explicit.ptr() );
310 
311  out << "\nM_explicit = " << *M_explicit;
312 
313  //
314  // F) Check the explicit operator
315  //
316 
317  LinearOpTester<double> M_explicit_tester;
318  M_explicit_tester.show_all_tests(true);
319 
320  const bool result = M_explicit_tester.compare( *M, *M_explicit, Teuchos::inoutArg(out) );
321  if (!result) success = false;
322  }
323 
324 }
325 
326 } // namespace
TEUCHOS_UNIT_TEST(EpetraOperatorWrapper, basic)
RCP< EpetraExtDiagScalingTransformer > epetraExtDiagScalingTransformer()
Nonmember constructor.
Teuchos::RCP< Epetra_Operator > get_Epetra_Operator(LinearOpBase< double > &op)
Full specialization for Scalar=double.
Transformer subclass for diagonally scaling a Epetra/Thyra operator.