Thyra  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
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
52 #include "Teuchos_CommandLineProcessor.hpp"
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"
58 #include "Teuchos_GlobalMPISession.hpp"
59 
60 
103 template<class Scalar>
104 int 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;
123  using Thyra::LinearOpBase;
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 //
303 int 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_THYRA_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 #ifdef HAVE_THYRA_COMPLEX
357 
358 #if defined(HAVE_THYRA_FLOAT)
359  // Run using std::complex<float>
360  result = exampleImplicitlyComposedLinearOperators<std::complex<float> >(
361  n0, n1, n2, *out, verbLevel, 1e-5, false );
362  if (!result) success = false;
363  // 2007/11/02: rabartl: Above, I skip the test of the adjoint since it
364  // fails but a lot. On my machine, the relative error is:
365  // rel_err((-3.00939,-0.836347),(-0.275689,1.45244)) = 1.14148.
366  // Since this works just fine for the next complex<double> case, I am
367  // going to just skip this test.
368 #endif // defined(HAVE_THYRA_FLOAT)
369 
370 
371  // Run using std::complex<double>
372  result = exampleImplicitlyComposedLinearOperators<std::complex<double> >(
373  n0, n1, n2, *out, verbLevel, 1e-12, true );
374  if (!result) success = false;
375 
376 #endif // HAVE_THYRA_COMPLEX
377 
378 #ifdef HAVE_TEUCHOS_GNU_MP
379 
380  // Run using mpf_class
381  result = exampleImplicitlyComposedLinearOperators<mpf_class>(
382  n0, n1, n2, *out, verbLevel, 1e-20, true );
383  if (!result) success = false;
384 
385 #endif // HAVE_TEUCHOS_GNU_MP
386 
387  }
388  TEUCHOS_STANDARD_CATCH_STATEMENTS(true,*out,success)
389 
390  if(success)
391  *out << "\nCongratulations! All of the tests checked out!\n";
392  else
393  *out << "\nOh no! At least one of the tests failed!\n";
394 
395  return success ? 0 : 1;
396 
397 } // end main()
Testing class for LinearOpBase.
basic_OSTab< char > OSTab
RCP< const LinearOpBase< Scalar > > diagonal(const RCP< VectorBase< Scalar > > &diag, const std::string &label="")
Nonmember constructor function.
Teuchos::RCP< const LinearOpBase< Scalar > > block2x1(const Teuchos::RCP< const LinearOpBase< Scalar > > &A00, const Teuchos::RCP< const LinearOpBase< Scalar > > &A10, const std::string &label="")
Form an implicit block 2x1 linear operator [ A00; A10 ].
Abstract interface for objects that represent a space for vectors.
RCP< const LinearOpBase< Scalar > > scale(const Scalar &scalar, const RCP< const LinearOpBase< Scalar > > &Op, const std::string &label="")
Build an implicit const scaled linear operator.
Teuchos::RCP< const LinearOpBase< Scalar > > block1x2(const Teuchos::RCP< const LinearOpBase< Scalar > > &A00, const Teuchos::RCP< const LinearOpBase< Scalar > > &A01, const std::string &label="")
Form an implicit block 1x2 linear operator [ A00, A01 ].
RCP< const LinearOpBase< Scalar > > add(const RCP< const LinearOpBase< Scalar > > &A, const RCP< const LinearOpBase< Scalar > > &B, const std::string &label="")
Form an implicit addition of two linear operators: M = A + B.
Interface for a collection of column vectors called a multi-vector.
static RCP< FancyOStream > getDefaultOStream()
RCP< VectorBase< Scalar > > createMember(const RCP< const VectorSpaceBase< Scalar > > &vs, const std::string &label="")
Create a vector member from the vector space.
RCP< const LinearOpBase< Scalar > > subtract(const RCP< const LinearOpBase< Scalar > > &A, const RCP< const LinearOpBase< Scalar > > &B, const std::string &label="")
Form an implicit subtraction of two linear operators: M = A - B.
Abstract interface for finite-dimensional dense vectors.
Base class for all linear operators.
RCP< const LinearOpBase< Scalar > > adjoint(const RCP< const LinearOpBase< Scalar > > &Op, const std::string &label="")
Build an implicit const adjoined linear operator.
void randomize(Scalar l, Scalar u, const Ptr< MultiVectorBase< Scalar > > &V)
Generate a random multi-vector with elements uniformly distributed elements.
Teuchos::RCP< const LinearOpBase< Scalar > > block2x2(const Teuchos::RCP< const LinearOpBase< Scalar > > &A00, const Teuchos::RCP< const LinearOpBase< Scalar > > &A01, const Teuchos::RCP< const LinearOpBase< Scalar > > &A10, const Teuchos::RCP< const LinearOpBase< Scalar > > &A11, const std::string &label="")
Form an implicit block 2x2 linear operator [ A00, A01; A10, A11 ].
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.
RCP< const LinearOpBase< Scalar > > identity(const RCP< const VectorSpaceBase< Scalar > > &space, const std::string &label="")
Create an identity linear operator with given a vector space.
TypeTo as(const TypeFrom &t)
#define TEUCHOS_STANDARD_CATCH_STATEMENTS(VERBOSE, ERR_STREAM, SUCCESS_FLAG)
RCP< MultiVectorBase< Scalar > > createMembers(const RCP< const VectorSpaceBase< Scalar > > &vs, int numMembers, const std::string &label="")
Create a set of vector members (a MultiVectorBase) from the vector space.
void set_all_error_tol(const ScalarMag error_tol)
Set all the error tolerances to the same value.
RCP< const LinearOpBase< Scalar > > multiply(const RCP< const LinearOpBase< Scalar > > &A, const RCP< const LinearOpBase< Scalar > > &B, const std::string &M_label="")
Form an implicit multiplication of two linear operators: M = A.