46#include "Thyra_OperatorVectorClientSupport.hpp"
49#include "Thyra_DefaultSpmdVectorSpace.hpp"
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"
103template<
class Scalar>
104int exampleImplicitlyComposedLinearOperators(
111 const bool testAdjoint
124 using Thyra::defaultSpmdVectorSpace;
125 using Thyra::randomize;
126 using Thyra::identity;
127 using Thyra::diagonal;
128 using Thyra::multiply;
130 using Thyra::subtract;
132 using Thyra::adjoint;
133 using Thyra::block1x2;
134 using Thyra::block2x2;
135 using Thyra::block2x2;
138 <<
"\n*** Demonstrating building linear operators for scalar type "
150 const RCP<const VectorSpaceBase<Scalar> >
151 space0 = defaultSpmdVectorSpace<Scalar>(n0),
152 space1 = defaultSpmdVectorSpace<Scalar>(n1),
153 space2 = defaultSpmdVectorSpace<Scalar>(n2);
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");
170 const RCP<VectorBase<Scalar> > d = createMember(space2);
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() );
194 const RCP<const LinearOpBase<Scalar> >
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);
224 const RCP<const LinearOpBase<Scalar> > I = identity(space1,
"I");
227 const RCP<const LinearOpBase<Scalar> > D = diagonal(d,
"D");
231 const RCP<const LinearOpBase<Scalar> > M00 =
234 add( scale(gamma,multiply(B,A)), C ), add( E, F ),
235 multiply(adjoint(J),A), I
240 out <<
"\nM00 = " << describe(*M00, verbLevel);
244 const RCP<const LinearOpBase<Scalar> > M01 =
251 out <<
"\nM01 = " << describe(*M01, verbLevel);
254 const RCP<const LinearOpBase<Scalar> > M10 =
256 multiply(L,adjoint(N)), scale(eta,P),
260 out <<
"\nM10 = " << describe(*M10, verbLevel);
263 const RCP<const LinearOpBase<Scalar> > M11 =
264 subtract( D, multiply(adjoint(Q),Q),
"M11" );
266 out <<
"\nM11 = " << describe(*M11, verbLevel);
271 const RCP<const LinearOpBase<Scalar> > M =
278 out <<
"\nM = " << describe(*M, verbLevel);
287 linearOpTester.check_adjoint(testAdjoint);
289 linearOpTester.show_all_tests(
true);
291 linearOpTester.dump_all(
true);
293 const bool result = linearOpTester.
check(*M, Teuchos::inOutArg(out));
303int main(
int argc,
char *argv[] )
323 CommandLineProcessor clp;
324 clp.throwExceptions(
false);
325 clp.addOutputSetupOptions(
true);
328 clp.setOption(
"n0", &n0 );
331 clp.setOption(
"n1", &n1 );
334 clp.setOption(
"n2", &n2 );
337 setVerbosityLevelOption(
"verb-level", &verbLevel,
338 "Top-level verbosity level. By default, this gets deincremented as you go deeper into numerical objects.",
341 CommandLineProcessor::EParseCommandLineReturn parse_return = clp.parse(argc,argv);
342 if( parse_return != CommandLineProcessor::PARSE_SUCCESSFUL )
return parse_return;
344#if defined(HAVE_TEUCHOS_INST_FLOAT)
346 result = exampleImplicitlyComposedLinearOperators<float>(
347 n0, n1, n2, *out, verbLevel, 1e-5,
true );
348 if (!result) success =
false;
352 result = exampleImplicitlyComposedLinearOperators<double>(
353 n0, n1, n2, *out, verbLevel, 1e-12,
true );
354 if (!result) success =
false;
356#if defined(HAVE_TEUCHOS_INST_COMPLEX_FLOAT) && defined(HAVE_TEUCHOS_INST_FLOAT)
358 result = exampleImplicitlyComposedLinearOperators<std::complex<float> >(
359 n0, n1, n2, *out, verbLevel, 1e-5, false );
360 if (!result) success =
false;
368#if defined(HAVE_TEUCHOS_INST_COMPLEX_DOUBLE)
370 result = exampleImplicitlyComposedLinearOperators<std::complex<double> >(
371 n0, n1, n2, *out, verbLevel, 1e-12, true );
372 if (!result) success =
false;
375#ifdef HAVE_TEUCHOS_GNU_MP
378 result = exampleImplicitlyComposedLinearOperators<mpf_class>(
379 n0, n1, n2, *out, verbLevel, 1e-20,
true );
380 if (!result) success =
false;
388 *out <<
"\nCongratulations! All of the tests checked out!\n";
390 *out <<
"\nOh no! At least one of the tests failed!\n";
392 return success ? 0 : 1;
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