Thyra Package Browser (Single Doxygen Collection) Version of the Day
Loading...
Searching...
No Matches
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"
55#include "EpetraExt_readEpetraLinearSystem.h"
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
73
74
75namespace {
76
77
78using Teuchos::null;
79using Teuchos::RCP;
81using Thyra::epetraExtDiagScalingTransformer;
82using Thyra::VectorBase;
83using Thyra::LinearOpBase;
84using Thyra::createMember;
85using Thyra::LinearOpTester;
86using Thyra::adjoint;
87using Thyra::multiply;
88using Thyra::diagonal;
89
90
91std::string matrixFile = "";
92
93
95{
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
103buildADOperator(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 }
140 out << "\nM = " << *M;
141
142 return M;
143}
144
145TEUCHOS_UNIT_TEST( EpetraExtDiagScalingTransformer, isCompatible)
146{
147 out << "\nReading linear system in Epetra format from the file \'"<<matrixFile<<"\'\n";
148
149#ifdef HAVE_MPI
150 Epetra_MpiComm comm(MPI_COMM_WORLD);
151#else
152 Epetra_SerialComm comm;
153#endif
154
155 RCP<Epetra_CrsMatrix> epetra_A;
156 RCP<Epetra_CrsMatrix> epetra_B;
157 EpetraExt::readEpetraLinearSystem( matrixFile, comm, &epetra_A, NULL, NULL, NULL );
158 EpetraExt::readEpetraLinearSystem( matrixFile, comm, &epetra_B, NULL, NULL, NULL );
159 const RCP<const Thyra::LinearOpBase<double> > A = Thyra::epetraLinearOp(epetra_A);
160 const RCP<const Thyra::LinearOpBase<double> > B = Thyra::epetraLinearOp(epetra_B);
161
162 RCP<VectorBase<double> > d = createMember(B->domain(), "d");
163 V_S( d.ptr(), 3.0 ); // ToDo: Set ton != 1.0 and generalize
164
165 const RCP<EpetraExtDiagScalingTransformer> BD_transformer =
166 epetraExtDiagScalingTransformer();
167
168 TEST_ASSERT(BD_transformer != null);
169
170 RCP<const LinearOpBase<double> > M;
171
172 M = multiply(A,B);
173 TEST_ASSERT(not BD_transformer->isCompatible(*M));
174
175 M = multiply(A,scale(3.9,diagonal(d)));
176 TEST_ASSERT(BD_transformer->isCompatible(*M));
177
178 M = multiply(scale(3.0,diagonal(d)),scale(9.0,B));
179 TEST_ASSERT(BD_transformer->isCompatible(*M));
180
181 M = multiply(B,scale(3.0,diagonal(d)),scale(9.0,B));
182 TEST_ASSERT(not BD_transformer->isCompatible(*M));
183}
184
185TEUCHOS_UNIT_TEST( EpetraExtDiagScalingTransformer, reapply_scaling)
186{
187 //
188 // A) Read in problem matrices
189 //
190
191 out << "\nReading linear system in Epetra format from the file \'"<<matrixFile<<"\'\n";
192
193#ifdef HAVE_MPI
194 Epetra_MpiComm comm(MPI_COMM_WORLD);
195#else
196 Epetra_SerialComm comm;
197#endif
198 RCP<Epetra_CrsMatrix> epetra_A;
199 EpetraExt::readEpetraLinearSystem( matrixFile, comm, &epetra_A, NULL, NULL, NULL );
200
201 //
202 // B) Create the Thyra wrapped version
203 //
204
205 const RCP<const Thyra::LinearOpBase<double> > A = Thyra::epetraLinearOp(epetra_A);
206
207 RCP<const Thyra::LinearOpBase<double> > M;
208
209 // test scale on right
210 for(int scenario=1;scenario<=2;scenario++) {
211 out << "**********************************************************" << std::endl;
212 out << "**************** Starting Scenario = " << scenario << std::endl;
213 out << "**********************************************************" << std::endl;
214
215 //
216 // D) Check compatibility
217 //
218 const RCP<EpetraExtDiagScalingTransformer> BD_transformer =
219 epetraExtDiagScalingTransformer();
220
221 TEST_ASSERT(BD_transformer != null);
222
223 //
224 // E) Do the transformation (twice)
225 //
226
227 const RCP<LinearOpBase<double> > M_explicit = BD_transformer->createOutputOp();
228 TEST_ASSERT(M_explicit != null);
229
230 // do trans
231 M = buildADOperator(scenario,A,4.5,out);
232 BD_transformer->transform( *M, M_explicit.ptr() );
233 const RCP<const Epetra_Operator> eOp_explicit0 = Thyra::get_Epetra_Operator(*M_explicit);
234
235 M = buildADOperator(scenario,A,7.5,out);
236 BD_transformer->transform( *M, M_explicit.ptr() );
237 const RCP<const Epetra_Operator> eOp_explicit1 = Thyra::get_Epetra_Operator(*M_explicit);
238
239 // Epetra operator should not change!
240 TEST_EQUALITY(eOp_explicit0,eOp_explicit1);
241
242 out << "\nM_explicit = " << *M_explicit;
243
244 //
245 // F) Check the explicit operator
246 //
247
248 LinearOpTester<double> M_explicit_tester;
249 M_explicit_tester.show_all_tests(true);;
250
251 const bool result = M_explicit_tester.compare( *M, *M_explicit, Teuchos::inoutArg(out) );
252 if (!result) success = false;
253 }
254}
255
256TEUCHOS_UNIT_TEST( EpetraExtDiagScalingTransformer, basic_BDG_GDB)
257{
258
259 //
260 // A) Read in problem matrices
261 //
262
263 out << "\nReading linear system in Epetra format from the file \'"<<matrixFile<<"\'\n";
264
265#ifdef HAVE_MPI
266 Epetra_MpiComm comm(MPI_COMM_WORLD);
267#else
268 Epetra_SerialComm comm;
269#endif
270 RCP<Epetra_CrsMatrix> epetra_A;
271 EpetraExt::readEpetraLinearSystem( matrixFile, comm, &epetra_A, NULL, NULL, NULL );
272
273 //
274 // B) Create the Thyra wrapped version
275 //
276
277 const RCP<const Thyra::LinearOpBase<double> > A = Thyra::epetraLinearOp(epetra_A);
278
279 //
280 // C) Create implicit B*D*B operator
281 //
282
283 // build scenario=1 -> B*D*B, scenario=2-> B*D*B', scenario=3 -> B'*D*B
284 for(int scenario=1;scenario<=4;scenario++) {
285 out << "**********************************************************" << std::endl;
286 out << "**************** Starting Scenario = " << scenario << std::endl;
287 out << "**********************************************************" << std::endl;
288
289 RCP<const Thyra::LinearOpBase<double> > M = buildADOperator(scenario,A,4.5,out);
290
291 //
292 // D) Check compatibility
293 //
294 const RCP<EpetraExtDiagScalingTransformer> BD_transformer =
295 epetraExtDiagScalingTransformer();
296
297 TEST_ASSERT(BD_transformer != null);
298
299 BD_transformer->isCompatible(*M);
300
301 //
302 // E) Do the transformation
303 //
304
305 const RCP<LinearOpBase<double> > M_explicit = BD_transformer->createOutputOp();
306 BD_transformer->transform( *M, M_explicit.ptr() );
307
308 out << "\nM_explicit = " << *M_explicit;
309
310 //
311 // F) Check the explicit operator
312 //
313
314 LinearOpTester<double> M_explicit_tester;
315 M_explicit_tester.show_all_tests(true);
316
317 const bool result = M_explicit_tester.compare( *M, *M_explicit, Teuchos::inoutArg(out) );
318 if (!result) success = false;
319 }
320
321}
322
323} // 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 diagonally scaling a Epetra/Thyra operator.
#define TEUCHOS_ASSERT(assertion_test)
TEST_ASSERT(castedDep1->getValuesAndValidators().size()==2)
TEST_EQUALITY(rcp_dynamic_cast< const EnhancedNumberValidator< double > >(castedDep1->getValuesAndValidators().find("val1") ->second, true) ->getMax(), double1Vali->getMax())
Teuchos::RCP< Epetra_Operator > get_Epetra_Operator(LinearOpBase< double > &op)
Full specialization for Scalar=double.