Thyra Version of the Day
Loading...
Searching...
No Matches
exampleImplicitlyComposedLinearOperators.cpp
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// Bring in all of the operator/vector ANA client support software
46#include "Thyra_OperatorVectorClientSupport.hpp"
47
48// Just use a default SPMD space for this example
49#include "Thyra_DefaultSpmdVectorSpace.hpp"
50
51// Some utilities from Teuchos
53#include "Teuchos_VerboseObject.hpp"
54#include "Teuchos_VerbosityLevelCommandLineProcessorHelpers.hpp"
55#include "Teuchos_Time.hpp"
56#include "Teuchos_StandardCatchMacros.hpp"
57#include "Teuchos_as.hpp"
59
60
103template<class Scalar>
104int exampleImplicitlyComposedLinearOperators(
105 const int n0,
106 const int n1,
107 const int n2,
109 const Teuchos::EVerbosityLevel verbLevel,
111 const bool testAdjoint
112 )
113{
114
115 // Using and other declarations
117 using Teuchos::as;
118 using Teuchos::RCP;
119 using Teuchos::OSTab;
121 using Thyra::VectorBase;
124 using Thyra::defaultSpmdVectorSpace;
125 using Thyra::randomize;
126 using Thyra::identity;
127 using Thyra::diagonal;
128 using Thyra::multiply;
129 using Thyra::add;
130 using Thyra::subtract;
131 using Thyra::scale;
132 using Thyra::adjoint;
133 using Thyra::block1x2;
134 using Thyra::block2x2;
135 using Thyra::block2x2;
136
137 out << "\n***"
138 << "\n*** Demonstrating building linear operators for scalar type "
139 << ST::name()
140 << "\n***\n";
141
142 OSTab tab(out);
143
144 //
145 // A) Set up the basic objects and other inputs to build the implicitly
146 // composed linear operators.
147 //
148
149 // Create serial vector spaces in this case
150 const RCP<const VectorSpaceBase<Scalar> >
151 space0 = defaultSpmdVectorSpace<Scalar>(n0),
152 space1 = defaultSpmdVectorSpace<Scalar>(n1),
153 space2 = defaultSpmdVectorSpace<Scalar>(n2);
154
155 // Create the component linear operators first as multi-vectors
156 const RCP<MultiVectorBase<Scalar> >
157 mvA = createMembers(space2, n0, "A"),
158 mvB = createMembers(space0, n2, "B"),
159 mvC = createMembers(space0, n0, "C"),
160 mvE = createMembers(space0, n1, "E"),
161 mvF = createMembers(space0, n1, "F"),
162 mvJ = createMembers(space2, n1, "J"),
163 mvK = createMembers(space1, n2, "K"),
164 mvL = createMembers(space2, n1, "L"),
165 mvN = createMembers(space0, n1, "N"),
166 mvP = createMembers(space2, n1, "P"),
167 mvQ = createMembers(space0, n2, "Q");
168
169 // Create the vector diagonal for D
170 const RCP<VectorBase<Scalar> > d = createMember(space2);
171
172 // Get the constants
173 const Scalar
174 one = 1.0,
175 beta = 2.0,
176 gamma = 3.0,
177 eta = 4.0;
178
179 // Randomize the values in the Multi-Vector
180 randomize( -one, +one, mvA.ptr() );
181 randomize( -one, +one, mvB.ptr() );
182 randomize( -one, +one, mvC.ptr() );
183 randomize( -one, +one, d.ptr() );
184 randomize( -one, +one, mvE.ptr() );
185 randomize( -one, +one, mvF.ptr() );
186 randomize( -one, +one, mvJ.ptr() );
187 randomize( -one, +one, mvK.ptr() );
188 randomize( -one, +one, mvL.ptr() );
189 randomize( -one, +one, mvN.ptr() );
190 randomize( -one, +one, mvP.ptr() );
191 randomize( -one, +one, mvQ.ptr() );
192
193 // Get the linear operator forms of the basic component linear operators
194 const RCP<const LinearOpBase<Scalar> >
195 A = mvA,
196 B = mvB,
197 C = mvC,
198 E = mvE,
199 F = mvF,
200 J = mvJ,
201 K = mvK,
202 L = mvL,
203 N = mvN,
204 P = mvP,
205 Q = mvQ;
206
207 out << describe(*A, verbLevel);
208 out << describe(*B, verbLevel);
209 out << describe(*C, verbLevel);
210 out << describe(*E, verbLevel);
211 out << describe(*F, verbLevel);
212 out << describe(*J, verbLevel);
213 out << describe(*K, verbLevel);
214 out << describe(*L, verbLevel);
215 out << describe(*N, verbLevel);
216 out << describe(*P, verbLevel);
217 out << describe(*Q, verbLevel);
218
219 //
220 // B) Create the composed linear operators
221 //
222
223 // I
224 const RCP<const LinearOpBase<Scalar> > I = identity(space1, "I");
225
226 // D = diag(d)
227 const RCP<const LinearOpBase<Scalar> > D = diagonal(d, "D");
228
229 // M00 = [ gama*B*A + C, E + F ] ^H
230 // [ J^H * A, I ]
231 const RCP<const LinearOpBase<Scalar> > M00 =
232 adjoint(
233 block2x2(
234 add( scale(gamma,multiply(B,A)), C ), add( E, F ),
235 multiply(adjoint(J),A), I
236 ),
237 "M00"
238 );
239
240 out << "\nM00 = " << describe(*M00, verbLevel);
241
242 // M01 = beta * [ Q ]
243 // [ K ]
244 const RCP<const LinearOpBase<Scalar> > M01 =
245 scale(
246 beta,
247 block2x1( Q, K ),
248 "M01"
249 );
250
251 out << "\nM01 = " << describe(*M01, verbLevel);
252
253 // M10 = [ L * N^H, eta*P ]
254 const RCP<const LinearOpBase<Scalar> > M10 =
255 block1x2(
256 multiply(L,adjoint(N)), scale(eta,P),
257 "M10"
258 );
259
260 out << "\nM10 = " << describe(*M10, verbLevel);
261
262 // M11 = D - Q^H*Q
263 const RCP<const LinearOpBase<Scalar> > M11 =
264 subtract( D, multiply(adjoint(Q),Q), "M11" );
265
266 out << "\nM11 = " << describe(*M11, verbLevel);
267
268
269 // M = [ M00, M01 ]
270 // [ M10, M11 ]
271 const RCP<const LinearOpBase<Scalar> > M =
272 block2x2(
273 M00, M01,
274 M10, M11,
275 "M"
276 );
277
278 out << "\nM = " << describe(*M, verbLevel);
279
280
281 //
282 // C) Test the final composed operator
283 //
284
285 Thyra::LinearOpTester<Scalar> linearOpTester;
286 linearOpTester.set_all_error_tol(errorTol);
287 linearOpTester.check_adjoint(testAdjoint);
288 if (as<int>(verbLevel) >= as<int>(Teuchos::VERB_HIGH))
289 linearOpTester.show_all_tests(true);
290 if (as<int>(verbLevel) >= as<int>(Teuchos::VERB_EXTREME))
291 linearOpTester.dump_all(true);
292
293 const bool result = linearOpTester.check(*M, Teuchos::inOutArg(out));
294
295 return result;
296
297}
298
299
300//
301// Main driver program
302//
303int main( int argc, char *argv[] )
304{
305
307
308 bool success = true;
309 bool result;
310
311 Teuchos::GlobalMPISession mpiSession(&argc,&argv);
312 // Above is needed to run in an MPI build with some MPI implementations
313
316
317 try {
318
319 //
320 // Read in command-line options
321 //
322
323 CommandLineProcessor clp;
324 clp.throwExceptions(false);
325 clp.addOutputSetupOptions(true);
326
327 int n0 = 2;
328 clp.setOption( "n0", &n0 );
329
330 int n1 = 3;
331 clp.setOption( "n1", &n1 );
332
333 int n2 = 4;
334 clp.setOption( "n2", &n2 );
335
337 setVerbosityLevelOption( "verb-level", &verbLevel,
338 "Top-level verbosity level. By default, this gets deincremented as you go deeper into numerical objects.",
339 &clp );
340
341 CommandLineProcessor::EParseCommandLineReturn parse_return = clp.parse(argc,argv);
342 if( parse_return != CommandLineProcessor::PARSE_SUCCESSFUL ) return parse_return;
343
344#if defined(HAVE_TEUCHOS_INST_FLOAT)
345 // Run using float
346 result = exampleImplicitlyComposedLinearOperators<float>(
347 n0, n1, n2, *out, verbLevel, 1e-5, true );
348 if (!result) success = false;
349#endif
350
351 // Run using double
352 result = exampleImplicitlyComposedLinearOperators<double>(
353 n0, n1, n2, *out, verbLevel, 1e-12, true );
354 if (!result) success = false;
355
356#if defined(HAVE_TEUCHOS_INST_COMPLEX_FLOAT) && defined(HAVE_TEUCHOS_INST_FLOAT)
357 // Run using std::complex<float>
358 result = exampleImplicitlyComposedLinearOperators<std::complex<float> >(
359 n0, n1, n2, *out, verbLevel, 1e-5, false );
360 if (!result) success = false;
361 // 2007/11/02: rabartl: Above, I skip the test of the adjoint since it
362 // fails but a lot. On my machine, the relative error is:
363 // rel_err((-3.00939,-0.836347),(-0.275689,1.45244)) = 1.14148.
364 // Since this works just fine for the next complex<double> case, I am
365 // going to just skip this test.
366#endif
367
368#if defined(HAVE_TEUCHOS_INST_COMPLEX_DOUBLE)
369 // Run using std::complex<double>
370 result = exampleImplicitlyComposedLinearOperators<std::complex<double> >(
371 n0, n1, n2, *out, verbLevel, 1e-12, true );
372 if (!result) success = false;
373#endif
374
375#ifdef HAVE_TEUCHOS_GNU_MP
376
377 // Run using mpf_class
378 result = exampleImplicitlyComposedLinearOperators<mpf_class>(
379 n0, n1, n2, *out, verbLevel, 1e-20, true );
380 if (!result) success = false;
381
382#endif // HAVE_TEUCHOS_GNU_MP
383
384 }
385 TEUCHOS_STANDARD_CATCH_STATEMENTS(true,*out,success)
386
387 if(success)
388 *out << "\nCongratulations! All of the tests checked out!\n";
389 else
390 *out << "\nOh no! At least one of the tests failed!\n";
391
392 return success ? 0 : 1;
393
394} // end main()
static RCP< FancyOStream > getDefaultOStream()
Base class for all linear operators.
Testing class for LinearOpBase.
void set_all_error_tol(const ScalarMag error_tol)
Set all the error tolerances to the same value.
bool check(const LinearOpBase< Scalar > &op, const Ptr< MultiVectorRandomizerBase< Scalar > > &rangeRandomizer, const Ptr< MultiVectorRandomizerBase< Scalar > > &domainRandomizer, const Ptr< FancyOStream > &out) const
Check a linear operator.
Interface for a collection of column vectors called a multi-vector.
Abstract interface for finite-dimensional dense vectors.
Abstract interface for objects that represent a space for vectors.
TypeTo as(const TypeFrom &t)
#define TEUCHOS_STANDARD_CATCH_STATEMENTS(VERBOSE, ERR_STREAM, SUCCESS_FLAG)
basic_OSTab< char > OSTab