Thyra Version of the Day
exampleImplicitlyComposedLinearOperators.cpp
00001 /*
00002 // @HEADER
00003 // ***********************************************************************
00004 // 
00005 //    Thyra: Interfaces and Support for Abstract Numerical Algorithms
00006 //                 Copyright (2004) Sandia Corporation
00007 // 
00008 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
00009 // license for use of this work by or on behalf of the U.S. Government.
00010 // 
00011 // Redistribution and use in source and binary forms, with or without
00012 // modification, are permitted provided that the following conditions are
00013 // met:
00014 //
00015 // 1. Redistributions of source code must retain the above copyright
00016 // notice, this list of conditions and the following disclaimer.
00017 //
00018 // 2. Redistributions in binary form must reproduce the above copyright
00019 // notice, this list of conditions and the following disclaimer in the
00020 // documentation and/or other materials provided with the distribution.
00021 //
00022 // 3. Neither the name of the Corporation nor the names of the
00023 // contributors may be used to endorse or promote products derived from
00024 // this software without specific prior written permission.
00025 //
00026 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
00027 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
00028 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
00029 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
00030 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
00031 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
00032 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
00033 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
00034 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
00035 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
00036 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
00037 //
00038 // Questions? Contact Roscoe A. Bartlett (bartlettra@ornl.gov) 
00039 // 
00040 // ***********************************************************************
00041 // @HEADER
00042 */
00043 
00044 
00045 // Bring in all of the operator/vector ANA client support software
00046 #include "Thyra_OperatorVectorClientSupport.hpp"
00047 
00048 // Just use a default SPMD space for this example
00049 #include "Thyra_DefaultSpmdVectorSpace.hpp"
00050 
00051 // Some utilities from Teuchos
00052 #include "Teuchos_CommandLineProcessor.hpp"
00053 #include "Teuchos_VerboseObject.hpp"
00054 #include "Teuchos_VerbosityLevelCommandLineProcessorHelpers.hpp"
00055 #include "Teuchos_Time.hpp"
00056 #include "Teuchos_StandardCatchMacros.hpp"
00057 #include "Teuchos_as.hpp"
00058 #include "Teuchos_GlobalMPISession.hpp"
00059 
00060 
00103 template<class Scalar>
00104 int exampleImplicitlyComposedLinearOperators(
00105   const int n0,
00106   const int n1,
00107   const int n2,
00108   Teuchos::FancyOStream &out,
00109   const Teuchos::EVerbosityLevel verbLevel,
00110   typename Teuchos::ScalarTraits<Scalar>::magnitudeType errorTol,
00111   const bool testAdjoint
00112   )
00113 {
00114 
00115   // Using and other declarations
00116   typedef Teuchos::ScalarTraits<Scalar> ST;
00117   using Teuchos::as;
00118   using Teuchos::RCP;
00119   using Teuchos::OSTab;
00120   using Thyra::VectorSpaceBase;
00121   using Thyra::VectorBase;
00122   using Thyra::MultiVectorBase;
00123   using Thyra::LinearOpBase;
00124   using Thyra::defaultSpmdVectorSpace;
00125   using Thyra::randomize;
00126   using Thyra::identity;
00127   using Thyra::diagonal;
00128   using Thyra::multiply;
00129   using Thyra::add;
00130   using Thyra::subtract;
00131   using Thyra::scale;
00132   using Thyra::adjoint;
00133   using Thyra::block1x2;
00134   using Thyra::block2x2;
00135   using Thyra::block2x2;
00136 
00137   out << "\n***"
00138       << "\n*** Demonstrating building linear operators for scalar type "
00139       << ST::name()
00140       << "\n***\n";
00141 
00142   OSTab tab(out);
00143 
00144   //
00145   // A) Set up the basic objects and other inputs to build the implicitly
00146   // composed linear operators.
00147   //
00148   
00149   // Create serial vector spaces in this case
00150   const RCP<const VectorSpaceBase<Scalar> >
00151     space0 = defaultSpmdVectorSpace<Scalar>(n0),
00152     space1 = defaultSpmdVectorSpace<Scalar>(n1),
00153     space2 = defaultSpmdVectorSpace<Scalar>(n2);
00154 
00155   // Create the component linear operators first as multi-vectors
00156   const RCP<MultiVectorBase<Scalar> >
00157     mvA = createMembers(space2, n0, "A"),
00158     mvB = createMembers(space0, n2, "B"),
00159     mvC = createMembers(space0, n0, "C"),
00160     mvE = createMembers(space0, n1, "E"),
00161     mvF = createMembers(space0, n1, "F"),
00162     mvJ = createMembers(space2, n1, "J"),
00163     mvK = createMembers(space1, n2, "K"),
00164     mvL = createMembers(space2, n1, "L"),
00165     mvN = createMembers(space0, n1, "N"),
00166     mvP = createMembers(space2, n1, "P"),
00167     mvQ = createMembers(space0, n2, "Q");
00168 
00169   // Create the vector diagonal for D
00170   const RCP<VectorBase<Scalar> > d = createMember(space2);
00171 
00172   // Get the constants
00173   const Scalar
00174     one = 1.0,
00175     beta = 2.0,
00176     gamma = 3.0,
00177     eta = 4.0;
00178 
00179   // Randomize the values in the Multi-Vector
00180   randomize( -one, +one, mvA.ptr() );
00181   randomize( -one, +one, mvB.ptr() );
00182   randomize( -one, +one, mvC.ptr() );
00183   randomize( -one, +one, d.ptr() );
00184   randomize( -one, +one, mvE.ptr() );
00185   randomize( -one, +one, mvF.ptr() );
00186   randomize( -one, +one, mvJ.ptr() );
00187   randomize( -one, +one, mvK.ptr() );
00188   randomize( -one, +one, mvL.ptr() );
00189   randomize( -one, +one, mvN.ptr() );
00190   randomize( -one, +one, mvP.ptr() );
00191   randomize( -one, +one, mvQ.ptr() );
00192 
00193   // Get the linear operator forms of the basic component linear operators
00194   const RCP<const LinearOpBase<Scalar> >
00195     A = mvA,
00196     B = mvB,
00197     C = mvC,
00198     E = mvE,
00199     F = mvF,
00200     J = mvJ,
00201     K = mvK,
00202     L = mvL,
00203     N = mvN,
00204     P = mvP,
00205     Q = mvQ;
00206 
00207   out << describe(*A, verbLevel);
00208   out << describe(*B, verbLevel);
00209   out << describe(*C, verbLevel);
00210   out << describe(*E, verbLevel);
00211   out << describe(*F, verbLevel);
00212   out << describe(*J, verbLevel);
00213   out << describe(*K, verbLevel);
00214   out << describe(*L, verbLevel);
00215   out << describe(*N, verbLevel);
00216   out << describe(*P, verbLevel);
00217   out << describe(*Q, verbLevel);
00218 
00219   //
00220   // B) Create the composed linear operators
00221   //
00222 
00223   // I
00224   const RCP<const LinearOpBase<Scalar> > I = identity(space1, "I");
00225 
00226   // D = diag(d)
00227   const RCP<const LinearOpBase<Scalar> > D = diagonal(d, "D");
00228 
00229   // M00 = [ gama*B*A + C,  E + F ] ^H
00230   //       [ J^H * A,       I     ]
00231   const RCP<const LinearOpBase<Scalar> > M00 =
00232     adjoint(
00233       block2x2(
00234         add( scale(gamma,multiply(B,A)), C ),  add( E, F ),
00235         multiply(adjoint(J),A),                I
00236         ),
00237       "M00"
00238       );
00239 
00240   out << "\nM00 = " << describe(*M00, verbLevel);
00241 
00242   // M01 = beta * [ Q ]
00243   //              [ K ]
00244   const RCP<const LinearOpBase<Scalar> > M01 =
00245     scale(
00246       beta,
00247       block2x1( Q, K ),
00248       "M01"
00249       );
00250 
00251   out << "\nM01 = "  << describe(*M01, verbLevel);
00252             
00253   // M10 = [ L * N^H,  eta*P ]
00254   const RCP<const LinearOpBase<Scalar> > M10 =
00255     block1x2(
00256       multiply(L,adjoint(N)),  scale(eta,P),
00257       "M10"
00258       );
00259 
00260   out << "\nM10 = " << describe(*M10, verbLevel);
00261 
00262   // M11 = D - Q^H*Q
00263   const RCP<const LinearOpBase<Scalar> > M11 =
00264     subtract( D, multiply(adjoint(Q),Q), "M11" );
00265 
00266   out << "\nM11 = "  << describe(*M11, verbLevel);
00267   
00268 
00269   // M = [ M00, M01 ]
00270   //     [ M10, M11 ]
00271   const RCP<const LinearOpBase<Scalar> > M =
00272     block2x2(
00273       M00, M01,
00274       M10, M11,
00275       "M"
00276       );
00277 
00278   out << "\nM = " << describe(*M, verbLevel);
00279 
00280   
00281   //
00282   // C) Test the final composed operator
00283   //
00284 
00285   Thyra::LinearOpTester<Scalar> linearOpTester;
00286   linearOpTester.set_all_error_tol(errorTol);
00287   linearOpTester.check_adjoint(testAdjoint);
00288   if (as<int>(verbLevel) >= as<int>(Teuchos::VERB_HIGH))
00289     linearOpTester.show_all_tests(true);
00290   if (as<int>(verbLevel) >= as<int>(Teuchos::VERB_EXTREME))
00291     linearOpTester.dump_all(true);
00292 
00293   const bool result = linearOpTester.check(*M, Teuchos::inOutArg(out));
00294 
00295   return result;
00296 
00297 }
00298 
00299 
00300 //
00301 // Main driver program
00302 //
00303 int main( int argc, char *argv[] )
00304 {
00305   
00306   using Teuchos::CommandLineProcessor;
00307 
00308   bool success = true;
00309   bool result;
00310 
00311   Teuchos::GlobalMPISession mpiSession(&argc,&argv);
00312   // Above is needed to run in an MPI build with some MPI implementations
00313 
00314   Teuchos::RCP<Teuchos::FancyOStream>
00315     out = Teuchos::VerboseObjectBase::getDefaultOStream();
00316 
00317   try {
00318 
00319     //
00320     // Read in command-line options
00321     //
00322 
00323     CommandLineProcessor  clp;
00324     clp.throwExceptions(false);
00325     clp.addOutputSetupOptions(true);
00326 
00327     int n0 = 2;
00328     clp.setOption( "n0", &n0 );
00329 
00330     int n1 = 3;
00331     clp.setOption( "n1", &n1 );
00332 
00333     int n2 = 4;
00334     clp.setOption( "n2", &n2 );
00335 
00336     Teuchos::EVerbosityLevel verbLevel = Teuchos::VERB_MEDIUM;
00337     setVerbosityLevelOption( "verb-level", &verbLevel,
00338       "Top-level verbosity level.  By default, this gets deincremented as you go deeper into numerical objects.",
00339       &clp );
00340 
00341     CommandLineProcessor::EParseCommandLineReturn parse_return = clp.parse(argc,argv);
00342     if( parse_return != CommandLineProcessor::PARSE_SUCCESSFUL ) return parse_return;
00343 
00344 #if defined(HAVE_THYRA_FLOAT)
00345     // Run using float
00346     result = exampleImplicitlyComposedLinearOperators<float>(
00347       n0, n1, n2, *out, verbLevel, 1e-5, true );
00348     if (!result) success = false;
00349 #endif
00350 
00351     // Run using double
00352     result = exampleImplicitlyComposedLinearOperators<double>(
00353       n0, n1, n2, *out, verbLevel, 1e-12, true );
00354     if (!result) success = false;
00355 
00356 #ifdef HAVE_THYRA_COMPLEX
00357 
00358 #if defined(HAVE_THYRA_FLOAT)
00359     // Run using std::complex<float>
00360     result = exampleImplicitlyComposedLinearOperators<std::complex<float> >(
00361       n0, n1, n2, *out, verbLevel, 1e-5, false );
00362     if (!result) success = false;
00363     // 2007/11/02: rabartl: Above, I skip the test of the adjoint since it
00364     // fails but a lot.  On my machine, the relative error is:
00365     //    rel_err((-3.00939,-0.836347),(-0.275689,1.45244)) = 1.14148.
00366     // Since this works just fine for the next complex<double> case, I am
00367     // going to just skip this test.
00368 #endif // defined(HAVE_THYRA_FLOAT)
00369 
00370 
00371     // Run using std::complex<double>
00372     result = exampleImplicitlyComposedLinearOperators<std::complex<double> >(
00373       n0, n1, n2, *out, verbLevel, 1e-12, true );
00374     if (!result) success = false;
00375 
00376 #endif // HAVE_THYRA_COMPLEX
00377 
00378 #ifdef HAVE_TEUCHOS_GNU_MP
00379 
00380     // Run using mpf_class
00381     result = exampleImplicitlyComposedLinearOperators<mpf_class>(
00382       n0, n1, n2, *out, verbLevel, 1e-20, true );
00383     if (!result) success = false;
00384 
00385 #endif // HAVE_TEUCHOS_GNU_MP
00386 
00387   }
00388   TEUCHOS_STANDARD_CATCH_STATEMENTS(true,*out,success)
00389     
00390   if(success)
00391     *out << "\nCongratulations! All of the tests checked out!\n";
00392   else
00393     *out << "\nOh no! At least one of the tests failed!\n";
00394   
00395   return success ? 0 : 1;
00396 
00397 } // end main()
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines