exampleImplicitlyComposedLinearOperators.cpp

00001 
00002 // Bring in all of the operator/vector ANA client support software
00003 #include "Thyra_OperatorVectorClientSupport.hpp"
00004 
00005 // Just use a default SPMD space for this example
00006 #include "Thyra_DefaultSpmdVectorSpace.hpp"
00007 
00008 // Some utilities from Teuchos
00009 #include "Teuchos_CommandLineProcessor.hpp"
00010 #include "Teuchos_VerboseObject.hpp"
00011 #include "Teuchos_VerbosityLevelCommandLineProcessorHelpers.hpp"
00012 #include "Teuchos_Time.hpp"
00013 #include "Teuchos_StandardCatchMacros.hpp"
00014 #include "Teuchos_as.hpp"
00015 
00058 template<class Scalar>
00059 int exampleImplicitlyComposedLinearOperators(
00060   const int n0,
00061   const int n1,
00062   const int n2,
00063   Teuchos::FancyOStream &out,
00064   const Teuchos::EVerbosityLevel verbLevel,
00065   typename Teuchos::ScalarTraits<Scalar>::magnitudeType errorTol,
00066   const bool testAdjoint
00067   )
00068 {
00069 
00070   // Using and other declarations
00071   typedef Teuchos::ScalarTraits<Scalar> ST;
00072   using Teuchos::as;
00073   using Teuchos::RCP;
00074   using Teuchos::OSTab;
00075   using Thyra::VectorSpaceBase;
00076   using Thyra::VectorBase;
00077   using Thyra::MultiVectorBase;
00078   using Thyra::LinearOpBase;
00079   using Thyra::defaultSpmdVectorSpace;
00080   using Thyra::randomize;
00081   using Thyra::identity;
00082   using Thyra::diagonal;
00083   using Thyra::multiply;
00084   using Thyra::add;
00085   using Thyra::subtract;
00086   using Thyra::scale;
00087   using Thyra::adjoint;
00088   using Thyra::block1x2;
00089   using Thyra::block2x2;
00090   using Thyra::block2x2;
00091 
00092   out << "\n***"
00093       << "\n*** Demonstrating building linear operators for scalar type "
00094       << ST::name()
00095       << "\n***\n";
00096 
00097   OSTab tab(out);
00098 
00099   //
00100   // A) Set up the basic objects and other inputs to build the implicitly
00101   // composed linear operators.
00102   //
00103   
00104   // Create serial vector spaces in this case
00105   const RCP<const VectorSpaceBase<Scalar> >
00106     space0 = defaultSpmdVectorSpace<Scalar>(n0),
00107     space1 = defaultSpmdVectorSpace<Scalar>(n1),
00108     space2 = defaultSpmdVectorSpace<Scalar>(n2);
00109 
00110   // Create the component linear operators first as multi-vectors
00111   const RCP<MultiVectorBase<Scalar> >
00112     mvA = createMembers(space2,n0,"A"),
00113     mvB = createMembers(space0,n2,"B"),
00114     mvC = createMembers(space0,n0,"C"),
00115     mvE = createMembers(space0,n1,"E"),
00116     mvF = createMembers(space0,n1,"F"),
00117     mvJ = createMembers(space2,n1,"J"),
00118     mvK = createMembers(space1,n2,"K"),
00119     mvL = createMembers(space2,n1,"L"),
00120     mvN = createMembers(space0,n1,"N"),
00121     mvP = createMembers(space2,n1,"P"),
00122     mvQ = createMembers(space0,n2,"Q");
00123 
00124   // Create the vector diagonal for D
00125   const RCP<VectorBase<Scalar> >
00126     d = createMember(space2);
00127 
00128   // Get the constants
00129   const Scalar
00130     one = 1.0,
00131     beta = 2.0,
00132     gamma = 3.0,
00133     eta = 4.0;
00134 
00135   // Randomize the values in the Multi-Vector
00136   randomize( -one, +one, &*mvA );
00137   randomize( -one, +one, &*mvB );
00138   randomize( -one, +one, &*mvC );
00139   randomize( -one, +one, &*d );
00140   randomize( -one, +one, &*mvE );
00141   randomize( -one, +one, &*mvF );
00142   randomize( -one, +one, &*mvJ );
00143   randomize( -one, +one, &*mvK );
00144   randomize( -one, +one, &*mvL );
00145   randomize( -one, +one, &*mvN );
00146   randomize( -one, +one, &*mvP );
00147   randomize( -one, +one, &*mvQ );
00148 
00149   // Get the linear operator forms of the basic component linear operators
00150   const RCP<const LinearOpBase<Scalar> >
00151     A = mvA,
00152     B = mvB,
00153     C = mvC,
00154     E = mvE,
00155     F = mvF,
00156     J = mvJ,
00157     K = mvK,
00158     L = mvL,
00159     N = mvN,
00160     P = mvP,
00161     Q = mvQ;
00162   
00163   //
00164   // B) Create the composed linear operators
00165   //
00166 
00167   // I
00168   const RCP<const LinearOpBase<Scalar> > I = identity(space1,"I");
00169 
00170   // D = diag(d)
00171   const RCP<const LinearOpBase<Scalar> > D = diagonal(d,"D");
00172 
00173   // M00 = [ gama*B*A + C,  E + F ] ^H
00174   //       [ J^H * A,       I     ]
00175   const RCP<const LinearOpBase<Scalar> > M00 =
00176     adjoint(
00177       block2x2(
00178         add( scale(gamma,multiply(B,A)), C ),  add( E, F ),
00179         multiply(adjoint(J),A),                I
00180         ),
00181       "M00"
00182       );
00183 
00184   out << "\nM00 = " << describe(*M00,verbLevel);
00185 
00186   // M01 = beta * [ Q ]
00187   //              [ K ]
00188   const RCP<const LinearOpBase<Scalar> > M01 =
00189     scale(
00190       beta,
00191       block2x1( Q, K ),
00192       "M01"
00193       );
00194 
00195   out << "\nM01 = "  << describe(*M01,verbLevel);
00196             
00197   // M10 = [ L * N^H,  eta*P ]
00198   const RCP<const LinearOpBase<Scalar> > M10 =
00199     block1x2(
00200       multiply(L,adjoint(N)),  scale(eta,P),
00201       "M10"
00202       );
00203 
00204   out << "\nM10 = " << describe(*M10,verbLevel);
00205 
00206   // M11 = D - Q^H*Q
00207   const RCP<const LinearOpBase<Scalar> > M11 =
00208     subtract( D, multiply(adjoint(Q),Q), "M11" );
00209 
00210   out << "\nM11 = "  << describe(*M11,verbLevel);
00211   
00212 
00213   // M = [ M00, M01 ]
00214   //     [ M10, M11 ]
00215   const RCP<const LinearOpBase<Scalar> > M =
00216     block2x2(
00217       M00, M01,
00218       M10, M11,
00219       "M"
00220       );
00221 
00222   out << "\nM = " << describe(*M,verbLevel);
00223 
00224   
00225   //
00226   // C) Test the final composed operator
00227   //
00228 
00229   Thyra::LinearOpTester<Scalar> linearOpTester;
00230   linearOpTester.set_all_error_tol(errorTol);
00231   linearOpTester.check_adjoint(testAdjoint);
00232   if (as<int>(verbLevel) >= as<int>(Teuchos::VERB_HIGH))
00233     linearOpTester.show_all_tests(true);
00234   if (as<int>(verbLevel) >= as<int>(Teuchos::VERB_EXTREME))
00235     linearOpTester.dump_all(true);
00236 
00237   const bool result = linearOpTester.check(*M,&out);
00238 
00239   return result;
00240 
00241 }
00242 
00243 
00244 //
00245 // Main driver program
00246 //
00247 int main( int argc, char *argv[] )
00248 {
00249   
00250   using Teuchos::CommandLineProcessor;
00251 
00252   bool success = true;
00253   bool result;
00254 
00255   Teuchos::RCP<Teuchos::FancyOStream>
00256     out = Teuchos::VerboseObjectBase::getDefaultOStream();
00257 
00258   try {
00259 
00260     //
00261     // Read in command-line options
00262     //
00263 
00264     CommandLineProcessor  clp;
00265     clp.throwExceptions(false);
00266     clp.addOutputSetupOptions(true);
00267 
00268     int n0 = 2;
00269     clp.setOption( "n0", &n0 );
00270 
00271     int n1 = 3;
00272     clp.setOption( "n1", &n1 );
00273 
00274     int n2 = 4;
00275     clp.setOption( "n2", &n2 );
00276 
00277     Teuchos::EVerbosityLevel verbLevel = Teuchos::VERB_MEDIUM;
00278     setVerbosityLevelOption( "verb-level", &verbLevel,
00279       "Top-level verbosity level.  By default, this gets deincremented as you go deeper into numerical objects.",
00280       &clp );
00281 
00282     CommandLineProcessor::EParseCommandLineReturn parse_return = clp.parse(argc,argv);
00283     if( parse_return != CommandLineProcessor::PARSE_SUCCESSFUL ) return parse_return;
00284 
00285     // Run using float
00286     result = exampleImplicitlyComposedLinearOperators<float>(
00287       n0, n1, n2, *out, verbLevel, 1e-5, true );
00288     if (!result) success = false;
00289 
00290     // Run using double
00291     result = exampleImplicitlyComposedLinearOperators<double>(
00292       n0, n1, n2, *out, verbLevel, 1e-12, true );
00293     if (!result) success = false;
00294 
00295 #if defined(HAVE_COMPLEX) && defined(HAVE_TEUCHOS_COMPLEX)
00296 
00297     // Run using std::complex<float>
00298     result = exampleImplicitlyComposedLinearOperators<std::complex<float> >(
00299       n0, n1, n2, *out, verbLevel, 1e-5, false );
00300     if (!result) success = false;
00301     // 2007/11/02: rabartl: Above, I skip the test of the adjoint since it
00302     // fails but a lot.  On my machine, the relative error is:
00303     //    rel_err((-3.00939,-0.836347),(-0.275689,1.45244)) = 1.14148.
00304     // Since this works just fine for the next complex<double> case, I am
00305     // going to just skip this test.
00306 
00307     // Run using std::complex<double>
00308     result = exampleImplicitlyComposedLinearOperators<std::complex<double> >(
00309       n0, n1, n2, *out, verbLevel, 1e-12, true );
00310     if (!result) success = false;
00311 
00312 #endif
00313 
00314 #ifdef HAVE_TEUCHOS_GNU_MP
00315 
00316     // Run using mpf_class
00317     result = exampleImplicitlyComposedLinearOperators<mpf_class>(
00318       n0, n1, n2, *out, verbLevel, 1e-20, true );
00319     if (!result) success = false;
00320 
00321 #endif // HAVE_TEUCHOS_GNU_MP
00322 
00323   }
00324   TEUCHOS_STANDARD_CATCH_STATEMENTS(true,*out,success)
00325     
00326   if(success)
00327     *out << "\nCongratulations! All of the tests checked out!\n";
00328   else
00329     *out << "\nOh no! At least one of the tests failed!\n";
00330   
00331   return success ? 0 : 1;
00332 
00333 } // end main()

Generated on Wed May 12 21:42:27 2010 for Thyra Operator/Vector Support by  doxygen 1.4.7