Thyra Package Browser (Single Doxygen Collection) Version of the Day
Loading...
Searching...
No Matches
EpetraExtAddTransformer_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_DefaultAddedLinearOp.hpp"
49#include "Thyra_DefaultDiagonalLinearOp.hpp"
50#include "Thyra_VectorStdOps.hpp"
51#include "Thyra_TestingTools.hpp"
52#include "Thyra_LinearOpTester.hpp"
56#include "Thyra_DefaultIdentityLinearOp.hpp"
57#include "EpetraExt_readEpetraLinearSystem.h"
63
64#ifdef HAVE_MPI
65# include "Epetra_MpiComm.h"
66#else
67# include "Epetra_SerialComm.h"
68#endif
69
71
72
73namespace {
74
75
76using Teuchos::null;
77using Teuchos::RCP;
79using Thyra::epetraExtAddTransformer;
80using Thyra::VectorBase;
81using Thyra::LinearOpBase;
82using Thyra::createMember;
83using Thyra::LinearOpTester;
84using Thyra::adjoint;
85using Thyra::multiply;
86using Thyra::diagonal;
87
88
89std::string matrixFile = "";
90std::string matrixFile2 = "";
91
92
94{
96 "matrix-file", &matrixFile,
97 "Defines the Epetra_CrsMatrix to read in." );
99 "matrix-file-2", &matrixFile2,
100 "Defines the Epetra_CrsMatrix to read in." );
101}
102
104buildAddOperator(int scenario,const Teuchos::RCP<const Thyra::LinearOpBase<double> > & A,
105 const Teuchos::RCP<const Thyra::LinearOpBase<double> > & B)
106{
107 // build operators for the various addition/adjoint scenarios
108 RCP<const Thyra::LinearOpBase<double> > M;
109
110 switch(scenario) {
111 case 0:
112 M = Thyra::add(A,B,"A+B");
113 break;
114 case 1:
115 M = Thyra::add(A,B,"A+adj(B)");
116 break;
117 case 2:
118 M = Thyra::add(A,B,"adj(A)+B");
119 break;
120 case 3:
121 M = Thyra::add(A,B,"adb(A)+adb(B)");
122 break;
123 default:
124 TEUCHOS_ASSERT(false);
125 }
126
127 return M;
128}
129
130TEUCHOS_UNIT_TEST( EpetraExtAddTransformer, diag_mat_Add )
131{
132
133 //
134 // A) Read in problem matrices
135 //
136
137 out << "\nReading linear system in Epetra format from the file \'"<<matrixFile<<"\' ...\n";
138
139#ifdef HAVE_MPI
140 Epetra_MpiComm comm(MPI_COMM_WORLD);
141#else
142 Epetra_SerialComm comm;
143#endif
144 RCP<Epetra_CrsMatrix> crsMat;
145 EpetraExt::readEpetraLinearSystem( matrixFile, comm, &crsMat, NULL, NULL, NULL );
146
147 //
148 // B) Create the Thyra wrapped version
149 //
150 double scaleMat=3.7;
151 //double scaleDiag=-2.9;
152
153
154 const RCP<const Thyra::LinearOpBase<double> > A = Thyra::scale<double>(scaleMat,Thyra::epetraLinearOp(crsMat));
155 RCP<VectorBase<double> > d = createMember(A->domain(), "d");
156 V_S( d.ptr(), 3.0 ); // ToDo: Set ton != 1.0 and generalize
157 const RCP<const Thyra::LinearOpBase<double> > B = diagonal(d);
158
159 out << "\nA = " << *A;
160 out << "\nB = " << *B;
161
162 for(int scenario=0;scenario<4;scenario++) {
163 //
164 // C) Create implicit A+B operator
165 //
166
167 const RCP<const Thyra::LinearOpBase<double> > M = buildAddOperator(scenario,A,B);
168
169 //
170 // D) Do the transformation
171 //
172
173 const RCP<EpetraExtAddTransformer> ApB_transformer = epetraExtAddTransformer();
174
175 TEST_ASSERT(ApB_transformer != null);
176
177 const RCP<LinearOpBase<double> > M_explicit = ApB_transformer->createOutputOp();
178 ApB_transformer->transform( *M, M_explicit.ptr() );
179
180 out << "\nM_explicit = " << *M_explicit;
181
182 //
183 // E) Check the explicit operator
184 //
185
186 LinearOpTester<double> M_explicit_tester;
187 M_explicit_tester.show_all_tests(true);;
188
189 const bool result = M_explicit_tester.compare( *M, *M_explicit, Teuchos::inOutArg(out) );
190 if (!result) success = false;
191 }
192
193 for(int scenario=0;scenario<4;scenario++) {
194 //
195 // C) Create implicit A+B operator
196 //
197
198 const RCP<const Thyra::LinearOpBase<double> > M = buildAddOperator(scenario,B,A);
199
200 //
201 // D) Do the transformation
202 //
203
204 const RCP<EpetraExtAddTransformer> ApB_transformer = epetraExtAddTransformer();
205
206 TEST_ASSERT(ApB_transformer != null);
207
208 const RCP<LinearOpBase<double> > M_explicit = ApB_transformer->createOutputOp();
209 ApB_transformer->transform( *M, M_explicit.ptr() );
210
211 out << "\nM_explicit = " << *M_explicit;
212
213 //
214 // E) Check the explicit operator
215 //
216
217 LinearOpTester<double> M_explicit_tester;
218 M_explicit_tester.show_all_tests(true);;
219
220 const bool result = M_explicit_tester.compare( *M, *M_explicit, Teuchos::inOutArg(out) );
221 if (!result) success = false;
222 }
223}
224
225TEUCHOS_UNIT_TEST( EpetraExtAddTransformer, id_mat_Add )
226{
227
228 //
229 // A) Read in problem matrices
230 //
231
232 out << "\nReading linear system in Epetra format from the file \'"<<matrixFile<<"\' ...\n";
233
234#ifdef HAVE_MPI
235 Epetra_MpiComm comm(MPI_COMM_WORLD);
236#else
237 Epetra_SerialComm comm;
238#endif
239 RCP<Epetra_CrsMatrix> crsMat;
240 EpetraExt::readEpetraLinearSystem( matrixFile, comm, &crsMat, NULL, NULL, NULL );
241
242 //
243 // B) Create the Thyra wrapped version
244 //
245 double scaleMat=3.7;
246 //double scaleDiag=-2.9;
247
248
249 const RCP<const Thyra::LinearOpBase<double> > A = Thyra::scale<double>(scaleMat,Thyra::epetraLinearOp(crsMat));
250 const RCP<const Thyra::LinearOpBase<double> > B = identity(A->domain(),"id");
251
252 out << "\nA = " << *A;
253 out << "\nB = " << *B;
254
255 for(int scenario=0;scenario<4;scenario++) {
256 //
257 // C) Create implicit A+B operator
258 //
259
260 const RCP<const Thyra::LinearOpBase<double> > M = buildAddOperator(scenario,A,B);
261
262 //
263 // D) Do the transformation
264 //
265
266 const RCP<EpetraExtAddTransformer> ApB_transformer = epetraExtAddTransformer();
267
268 TEST_ASSERT(ApB_transformer != null);
269
270 const RCP<LinearOpBase<double> > M_explicit = ApB_transformer->createOutputOp();
271 ApB_transformer->transform( *M, M_explicit.ptr() );
272
273 out << "\nM_explicit = " << *M_explicit;
274
275 //
276 // E) Check the explicit operator
277 //
278
279 LinearOpTester<double> M_explicit_tester;
280 M_explicit_tester.show_all_tests(true);;
281
282 const bool result = M_explicit_tester.compare( *M, *M_explicit, Teuchos::inOutArg(out) );
283 if (!result) success = false;
284 }
285
286 for(int scenario=0;scenario<4;scenario++) {
287 //
288 // C) Create implicit A+B operator
289 //
290
291 const RCP<const Thyra::LinearOpBase<double> > M = buildAddOperator(scenario,B,A);
292
293 //
294 // D) Do the transformation
295 //
296
297 const RCP<EpetraExtAddTransformer> ApB_transformer = epetraExtAddTransformer();
298
299 TEST_ASSERT(ApB_transformer != null);
300
301 const RCP<LinearOpBase<double> > M_explicit = ApB_transformer->createOutputOp();
302 ApB_transformer->transform( *M, M_explicit.ptr() );
303
304 out << "\nM_explicit = " << *M_explicit;
305
306 //
307 // E) Check the explicit operator
308 //
309
310 LinearOpTester<double> M_explicit_tester;
311 M_explicit_tester.show_all_tests(true);;
312
313 const bool result = M_explicit_tester.compare( *M, *M_explicit, Teuchos::inOutArg(out) );
314 if (!result) success = false;
315 }
316}
317
318TEUCHOS_UNIT_TEST( EpetraExtAddTransformer, basic_Add )
319{
320
321 //
322 // A) Read in problem matrices
323 //
324
325 out << "\nReading linear system in Epetra format from the file \'"<<matrixFile<<"\' ...\n";
326 out << "\nReading linear system in Epetra format from the file \'"<<matrixFile2<<"\' ...\n";
327
328#ifdef HAVE_MPI
329 Epetra_MpiComm comm(MPI_COMM_WORLD);
330#else
331 Epetra_SerialComm comm;
332#endif
333 RCP<Epetra_CrsMatrix> epetra_A;
334 RCP<Epetra_CrsMatrix> epetra_B;
335 EpetraExt::readEpetraLinearSystem( matrixFile, comm, &epetra_A, NULL, NULL, NULL );
336 EpetraExt::readEpetraLinearSystem( matrixFile2, comm, &epetra_B, NULL, NULL, NULL );
337
338 //
339 // B) Create the Thyra wrapped version
340 //
341 double scaleA=3.7;
342 double scaleB=-2.9;
343
344 const RCP<const Thyra::LinearOpBase<double> > A = Thyra::scale<double>(scaleA,Thyra::epetraLinearOp(epetra_B));
345 const RCP<const Thyra::LinearOpBase<double> > B = Thyra::scale<double>(scaleB,Thyra::epetraLinearOp(epetra_B));
346
347 out << "\nA = " << *A;
348 out << "\nB = " << *B;
349
350 for(int scenario=0;scenario<4;scenario++) {
351 //
352 // C) Create implicit A+B operator
353 //
354
355 const RCP<const Thyra::LinearOpBase<double> > M = buildAddOperator(scenario,A,B);
356
357 //
358 // D) Do the transformation
359 //
360
361 const RCP<EpetraExtAddTransformer> ApB_transformer = epetraExtAddTransformer();
362
363 TEST_ASSERT(ApB_transformer != null);
364
365 const RCP<LinearOpBase<double> > M_explicit = ApB_transformer->createOutputOp();
366 ApB_transformer->transform( *M, M_explicit.ptr() );
367
368 out << "\nM_explicit = " << *M_explicit;
369
370 //
371 // E) Check the explicit operator
372 //
373
374 LinearOpTester<double> M_explicit_tester;
375 M_explicit_tester.show_all_tests(true);;
376
377 const bool result = M_explicit_tester.compare( *M, *M_explicit, Teuchos::inOutArg(out) );
378 if (!result) success = false;
379 }
380}
381
382TEUCHOS_UNIT_TEST( EpetraExtAddTransformer, mod_Add )
383{
384
385 //
386 // A) Read in problem matrices
387 //
388
389 out << "\nReading linear system in Epetra format from the file \'"<<matrixFile<<"\' ...\n";
390 out << "\nReading linear system in Epetra format from the file \'"<<matrixFile2<<"\' ...\n";
391
392#ifdef HAVE_MPI
393 Epetra_MpiComm comm(MPI_COMM_WORLD);
394#else
395 Epetra_SerialComm comm;
396#endif
397 RCP<Epetra_CrsMatrix> epetra_A;
398 RCP<Epetra_CrsMatrix> epetra_B;
399 EpetraExt::readEpetraLinearSystem( matrixFile, comm, &epetra_A, NULL, NULL, NULL );
400 EpetraExt::readEpetraLinearSystem( matrixFile2, comm, &epetra_B, NULL, NULL, NULL );
401
402 //
403 // B) Create the Thyra wrapped version
404 //
405 double scaleA=3.7;
406 double scaleB=-2.9;
407
408 const RCP<const Thyra::LinearOpBase<double> > A = Thyra::scale<double>(scaleA,Thyra::epetraLinearOp(epetra_B));
409 const RCP<const Thyra::LinearOpBase<double> > B = Thyra::scale<double>(scaleB,Thyra::epetraLinearOp(epetra_B));
410
411 out << "\nA = " << *A;
412 out << "\nB = " << *B;
413
414 for(int scenario=0;scenario<4;scenario++) {
415 //
416 // C) Create implicit A+B operator
417 //
418
419 const RCP<const Thyra::LinearOpBase<double> > M = buildAddOperator(scenario,A,B);
420
421 //
422 // D) Do the transformation
423 //
424
425 const RCP<EpetraExtAddTransformer> ApB_transformer = epetraExtAddTransformer();
426
427 TEST_ASSERT(ApB_transformer != null);
428
429 const RCP<LinearOpBase<double> > M_explicit = ApB_transformer->createOutputOp();
430 const RCP<LinearOpBase<double> > M_explicit_orig = M_explicit;
431 ApB_transformer->transform( *M, M_explicit.ptr() );
432
433 // do some violence to one of the operators
434 double * view; int numEntries;
435 epetra_B->Scale(3.2);
436 TEUCHOS_ASSERT(epetra_B->ExtractMyRowView(3,numEntries,view)==0);
437 for(int i=0;i<numEntries;i++) view[i] += view[i]*double(i+1);
438
439 // compute multiply again
440 ApB_transformer->transform( *M, M_explicit.ptr() );
441
442 out << "\nM_explicit = " << *M_explicit;
443
445 ==Thyra::get_Epetra_Operator(*M_explicit_orig));
446
447 //
448 // E) Check the explicit operator
449 //
450
451 LinearOpTester<double> M_explicit_tester;
452 M_explicit_tester.show_all_tests(true);;
453
454 const bool result = M_explicit_tester.compare( *M, *M_explicit, Teuchos::inoutArg(out) );
455 if (!result) success = false;
456 }
457}
458
459} // end namespace
TEUCHOS_UNIT_TEST(Comm, Issue1029)
TEUCHOS_STATIC_SETUP()
void setOption(const char option_true[], const char option_false[], bool *option_val, const char documentation[]=NULL)
static CommandLineProcessor & getCLP()
Transformer subclass for adding Epetra/Thyra operators using EpetraExt::MatrixMatrix.
#define TEUCHOS_ASSERT(assertion_test)
TEST_ASSERT(castedDep1->getValuesAndValidators().size()==2)
Teuchos::RCP< Epetra_Operator > get_Epetra_Operator(LinearOpBase< double > &op)
Full specialization for Scalar=double.