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