Thyra Package Browser (Single Doxygen Collection) Version of the Day
Loading...
Searching...
No Matches
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
70namespace {
71
72
73using Teuchos::null;
74using Teuchos::RCP;
75using Thyra::EpetraExtDiagScaledMatProdTransformer;
76using Thyra::epetraExtDiagScaledMatProdTransformer;
77using Thyra::VectorBase;
78using Thyra::LinearOpBase;
79using Thyra::createMember;
80using Thyra::LinearOpTester;
81using Thyra::adjoint;
82using Thyra::multiply;
83using Thyra::diagonal;
84
85
86std::string matrixFile = "";
87std::string matrixFile2 = "";
88
89
90TEUCHOS_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
101const Teuchos::RCP<const Thyra::LinearOpBase<double> >
102buildBDGOperator(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
152const Teuchos::RCP<const Thyra::LinearOpBase<double> >
153buildBGOperator(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
190const Teuchos::RCP<const Thyra::LinearOpBase<double> >
191buildBDBOperator(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
199TEUCHOS_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
260TEUCHOS_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
327TEUCHOS_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
388TEUCHOS_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_UNIT_TEST
TEUCHOS_UNIT_TEST(Comm, Issue1029)
TEUCHOS_STATIC_SETUP
TEUCHOS_STATIC_SETUP()
Teuchos_CommandLineProcessor.hpp
Teuchos_GlobalMPISession.hpp
TEST_ASSERT
#define TEST_ASSERT(v1)
Teuchos_StandardCatchMacros.hpp
Teuchos_UnitTestHarness.hpp
Teuchos_VerboseObject.hpp
Teuchos_XMLParameterListHelpers.hpp
Thyra_EpetraExtDiagScaledMatProdTransformer.hpp
Thyra_EpetraLinearOp.hpp
Thyra_EpetraThyraWrappers.hpp
RCP::ptr
Ptr< T > ptr() const
Teuchos::CommandLineProcessor::setOption
void setOption(const char option_true[], const char option_false[], bool *option_val, const char documentation[]=NULL)
Teuchos::UnitTestRepository::getCLP
static CommandLineProcessor & getCLP()
Thyra::EpetraExtDiagScaledMatProdTransformer
Transformer subclass for diagonally scaling and multiplying Epetra/Thyra operators.
Definition Thyra_EpetraExtDiagScaledMatProdTransformer.hpp:58
Thyra::EpetraExtDiagScaledMatProdTransformer::epetraExtDiagScaledMatProdTransformer
RCP< EpetraExtDiagScaledMatProdTransformer > epetraExtDiagScaledMatProdTransformer()
Nonmember constructor.
Definition Thyra_EpetraExtDiagScaledMatProdTransformer.hpp:89
Thyra::EpetraLinearOp::epetraLinearOp
RCP< const EpetraLinearOp > epetraLinearOp(const RCP< const Epetra_Operator > &op, EOpTransp opTrans=NOTRANS, EApplyEpetraOpAs applyAs=EPETRA_OP_APPLY_APPLY, EAdjointEpetraOp adjointSupport=EPETRA_OP_ADJOINT_SUPPORTED, const RCP< const VectorSpaceBase< double > > &range=Teuchos::null, const RCP< const VectorSpaceBase< double > > &domain=Teuchos::null)
Dynamically allocate a nonconst EpetraLinearOp to wrap a const Epetra_Operator object.
TEUCHOS_ASSERT
#define TEUCHOS_ASSERT(assertion_test)
Generated by doxygen 1.13.2