test_composite_linear_ops.cpp

00001 // @HEADER
00002 // ***********************************************************************
00003 // 
00004 //    Thyra: Interfaces and Support for Abstract Numerical Algorithms
00005 //                 Copyright (2004) Sandia Corporation
00006 // 
00007 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
00008 // license for use of this work by or on behalf of the U.S. Government.
00009 // 
00010 // This library is free software; you can redistribute it and/or modify
00011 // it under the terms of the GNU Lesser General Public License as
00012 // published by the Free Software Foundation; either version 2.1 of the
00013 // License, or (at your option) any later version.
00014 //  
00015 // This library is distributed in the hope that it will be useful, but
00016 // WITHOUT ANY WARRANTY; without even the implied warranty of
00017 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00018 // Lesser General Public License for more details.
00019 //  
00020 // You should have received a copy of the GNU Lesser General Public
00021 // License along with this library; if not, write to the Free Software
00022 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
00023 // USA
00024 // Questions? Contact Michael A. Heroux (maherou@sandia.gov) 
00025 // 
00026 // ***********************************************************************
00027 // @HEADER
00028 
00029 #include "Thyra_DefaultSpmdVectorSpace.hpp"
00030 #include "Thyra_DefaultZeroLinearOp.hpp"
00031 #include "Thyra_DefaultIdentityLinearOp.hpp"
00032 #include "Thyra_DefaultScaledAdjointLinearOp.hpp"
00033 #include "Thyra_DefaultAddedLinearOp.hpp"
00034 #include "Thyra_DefaultMultipliedLinearOp.hpp"
00035 #include "Thyra_DefaultBlockedLinearOp.hpp"
00036 #include "Thyra_VectorStdOps.hpp"
00037 #include "Thyra_MultiVectorStdOps.hpp"
00038 #include "Thyra_DefaultDiagonalLinearOp.hpp"
00039 #include "Thyra_TestingTools.hpp"
00040 #include "Thyra_LinearOpTester.hpp"
00041 #include "Teuchos_CommandLineProcessor.hpp"
00042 #include "Teuchos_GlobalMPISession.hpp"
00043 #include "Teuchos_VerboseObject.hpp"
00044 #include "Teuchos_DefaultComm.hpp"
00045 #include "Teuchos_dyn_cast.hpp"
00046 #include "Teuchos_StandardCatchMacros.hpp"
00047 
00050 template <class Scalar>
00051 bool run_composite_linear_ops_tests(
00052   const Teuchos::RCP<const Teuchos::Comm<Thyra::Index> >   comm
00053   ,const int                                                       n
00054   ,const bool                                                      useSpmd
00055   ,const typename Teuchos::ScalarTraits<Scalar>::magnitudeType     &tol
00056   ,const bool                                                      dumpAll
00057   ,Teuchos::FancyOStream                                           *out_arg
00058   )
00059 {
00060 
00061   using Thyra::relErr;
00062   using Thyra::passfail;
00063   typedef Teuchos::ScalarTraits<Scalar> ST;
00064   typedef typename ST::magnitudeType    ScalarMag;
00065   typedef Teuchos::ScalarTraits<ScalarMag> STM;
00066   using Teuchos::RCP;
00067   using Teuchos::rcp;
00068   using Teuchos::null;
00069   using Teuchos::rcp_const_cast;
00070   using Teuchos::rcp_dynamic_cast;
00071   using Teuchos::dyn_cast;
00072   using Teuchos::OSTab;
00073 
00074   RCP<Teuchos::FancyOStream>
00075     out = rcp(new Teuchos::FancyOStream(rcp(out_arg,false)));
00076 
00077   const Teuchos::EVerbosityLevel
00078     verbLevel = dumpAll?Teuchos::VERB_EXTREME:Teuchos::VERB_HIGH;
00079 
00080   if(out.get()) *out << "\n*** Entering run_composite_linear_ops_tests<"<<ST::name()<<">(...) ...\n";
00081 
00082   bool success = true, result;
00083 
00084   const ScalarMag warning_tol = ScalarMag(1e-2)*tol, error_tol = tol;
00085   Thyra::LinearOpTester<Scalar> linearOpTester;
00086   linearOpTester.linear_properties_warning_tol(warning_tol);
00087   linearOpTester.linear_properties_error_tol(error_tol);
00088   linearOpTester.adjoint_warning_tol(warning_tol);
00089   linearOpTester.adjoint_error_tol(error_tol);
00090   linearOpTester.dump_all(dumpAll);
00091   Thyra::LinearOpTester<Scalar> symLinearOpTester(linearOpTester);
00092   symLinearOpTester.check_for_symmetry(true);
00093   symLinearOpTester.symmetry_warning_tol(STM::squareroot(warning_tol));
00094   symLinearOpTester.symmetry_error_tol(STM::squareroot(error_tol));
00095 
00096   RCP<const Thyra::VectorSpaceBase<Scalar> > space;
00097   if(useSpmd) space = rcp(new Thyra::DefaultSpmdVectorSpace<Scalar>(comm,n,-1));
00098   else       space = rcp(new Thyra::DefaultSpmdVectorSpace<Scalar>(n));
00099   if(out.get()) *out << "\nUsing a basic vector space described as " << describe(*space,verbLevel) << " ...\n";
00100   
00101   if(out.get()) *out << "\nCreating random n x (n/2) multi-vector origA ...\n";
00102   RCP<Thyra::MultiVectorBase<Scalar> >
00103     mvOrigA = createMembers(space,n/2,"origA");
00104   Thyra::seed_randomize<Scalar>(0);
00105   Thyra::randomize( Scalar(Scalar(-1)*ST::one()), Scalar(Scalar(+1)*ST::one()), &*mvOrigA );
00106   RCP<const Thyra::LinearOpBase<Scalar> >
00107     origA = mvOrigA;
00108   if(out.get()) *out << "\norigA =\n" << describe(*origA,verbLevel);
00109 
00110   if(out.get()) *out << "\nTesting origA ...\n";
00111   Thyra::seed_randomize<Scalar>(0);
00112   result = linearOpTester.check(*origA,out.get());
00113   if(!result) success = false;
00114   
00115   if(out.get()) *out << "\nCreating implicit scaled linear operator A1 = scale(0.5,origA) ...\n";
00116   RCP<const Thyra::LinearOpBase<Scalar> >
00117     A1 = scale(Scalar(0.5),origA);
00118   if(out.get()) *out << "\nA1 =\n" << describe(*A1,verbLevel);
00119 
00120   if(out.get()) *out << "\nTesting A1 ...\n";
00121   Thyra::seed_randomize<Scalar>(0);
00122   result = linearOpTester.check(*A1,out.get());
00123   if(!result) success = false;
00124 
00125   if(out.get()) *out << "\nTesting that A1.getOp() == origA ...\n";
00126   Thyra::seed_randomize<Scalar>(0);
00127   result = linearOpTester.compare(*dyn_cast<const Thyra::DefaultScaledAdjointLinearOp<Scalar> >(*A1).getOp(),*origA,out.get());
00128   if(!result) success = false;
00129 
00130   {
00131 
00132     if(out.get()) *out << "\nUnwrapping origA to get non-persisting pointer to origA_1, scalar and transp ...\n";
00133     Scalar  scalar;
00134     Thyra::ETransp transp;
00135     const Thyra::LinearOpBase<Scalar> *origA_1 = NULL;
00136     unwrap( *origA, &scalar, &transp, &origA_1 );
00137     TEST_FOR_EXCEPT( origA_1 == NULL );
00138 
00139     if(out.get()) *out << "\nscalar = " << scalar << " == 1 ? ";
00140     result = (scalar == ST::one());
00141     if(!result) success = false;
00142     if(out.get()) *out << passfail(result) << std::endl;
00143 
00144     if(out.get()) *out << "\ntransp = " << toString(transp) << " == NOTRANS ? ";
00145     result = (transp == Thyra::NOTRANS);
00146     if(!result) success = false;
00147     if(out.get()) *out << passfail(result) << std::endl;
00148     
00149     if(out.get()) *out << "\nTesting that origA_1 == origA ...\n";
00150     Thyra::seed_randomize<Scalar>(0);
00151     result = linearOpTester.compare(*origA_1,*origA,out.get());
00152     if(!result) success = false;
00153     
00154   }
00155 
00156   {
00157 
00158     if(out.get()) *out << "\nUnwrapping A1 to get non-persisting pointer to origA_2 ...\n";
00159     Scalar  scalar;
00160     Thyra::ETransp transp;
00161     const Thyra::LinearOpBase<Scalar> *origA_2 = NULL;
00162     unwrap( *A1, &scalar, &transp, &origA_2 );
00163     TEST_FOR_EXCEPT( origA_2 == NULL );
00164 
00165     if(out.get()) *out << "\nscalar = " << scalar << " == 0.5 ? ";
00166     result = (scalar == Scalar(0.5));
00167     if(!result) success = false;
00168     if(out.get()) *out << passfail(result) << std::endl;
00169 
00170     if(out.get()) *out << "\ntransp = " << toString(transp) << " == NOTRANS ? ";
00171     result = (transp == Thyra::NOTRANS);
00172     if(!result) success = false;
00173     if(out.get()) *out << passfail(result) << std::endl;
00174 
00175     if(out.get()) *out << "\nTesting that origA_2 == origA ...\n";
00176     Thyra::seed_randomize<Scalar>(0);
00177     result = linearOpTester.compare(*origA_2,*origA,out.get());
00178     if(!result) success = false;
00179 
00180   }
00181   
00182   if(out.get()) *out << "\nCreating implicit scaled linear operator A2 = adjoint(A1) ...\n";
00183   RCP<const Thyra::LinearOpBase<Scalar> >
00184     A2 = adjoint(A1);
00185   if(out.get()) *out << "\nA2 =\n" << describe(*A2,verbLevel);
00186 
00187   if(out.get()) *out << "\nTesting A2 ...\n";
00188   Thyra::seed_randomize<Scalar>(0);
00189   result = linearOpTester.check(*A2,out.get());
00190   if(!result) success = false;
00191 
00192   if(out.get()) *out << "\nTesting that A2.getOp() == A1 ...\n";
00193   Thyra::seed_randomize<Scalar>(0);
00194   result = linearOpTester.compare(*dyn_cast<const Thyra::DefaultScaledAdjointLinearOp<Scalar> >(*A2).getOp(),*A1,out.get());
00195   if(!result) success = false;
00196   
00197   if(out.get()) *out << "\nCreating implicit scaled, adjoined linear operator A3 = adjoint(scale(2.0,(A2)) ...\n";
00198   RCP<const Thyra::LinearOpBase<Scalar> >
00199     A3 = adjoint(scale(Scalar(2.0),A2));
00200   if(out.get()) *out << "\nA3 =\n" << describe(*A3,verbLevel);
00201 
00202   if(out.get()) *out << "\nTesting A3 ...\n";
00203   Thyra::seed_randomize<Scalar>(0);
00204   result = linearOpTester.check(*A3,out.get());
00205   if(!result) success = false;
00206 
00207   if(out.get()) *out << "\nTesting that A3 == origA ...\n";
00208   Thyra::seed_randomize<Scalar>(0);
00209   result = linearOpTester.compare(*A3,*origA,out.get());
00210   if(!result) success = false;
00211 
00212   if(out.get()) *out << "\nCalling all of the rest of the functions for non-const just to test them ...\n";
00213   RCP<Thyra::LinearOpBase<Scalar> >
00214     A4 = nonconstScale(
00215       Scalar(0.25)
00216       ,nonconstAdjoint(
00217         nonconstTranspose(
00218           nonconstAdjoint(
00219             nonconstScaleAndAdjoint(
00220               Scalar(4.0)
00221               ,Thyra::TRANS
00222               ,Teuchos::rcp_const_cast<Thyra::LinearOpBase<Scalar> >(origA)
00223               )
00224             )
00225           )
00226         )
00227       );
00228   if(!ST::isComplex) A4 = nonconstTranspose(nonconstAdjoint(A4)); // Should result in CONJ
00229   if(out.get()) *out << "\nA4 =\n" << describe(*A4,verbLevel);
00230 
00231   if(out.get()) *out << "\nTesting A4 ...\n";
00232   Thyra::seed_randomize<Scalar>(0);
00233   result = linearOpTester.check(*A4,out.get());
00234   if(!result) success = false;
00235 
00236   if(out.get()) *out << "\nCalling all of the rest of the functions for const just to test them ...\n";
00237   RCP<const Thyra::LinearOpBase<Scalar> >
00238     A5 = scale(
00239       Scalar(0.25)
00240       ,adjoint(
00241         transpose(
00242           adjoint(
00243             scaleAndAdjoint(
00244               Scalar(4.0)
00245               ,Thyra::TRANS
00246               ,origA
00247               )
00248             )
00249           )
00250         )
00251       );
00252   if(!ST::isComplex) A5 = transpose(adjoint(A5)); // Should result in CONJ
00253   if(out.get()) *out << "\nA5 =\n" << describe(*A5,verbLevel);
00254 
00255   if(out.get()) *out << "\nTesting A5 ...\n";
00256   Thyra::seed_randomize<Scalar>(0);
00257   result = linearOpTester.check(*A5,out.get());
00258   if(!result) success = false;
00259 
00260   if(out.get()) *out << "\nCreating a multiplied operator A6 = origA^H*A1 ...\n";
00261   RCP<const Thyra::LinearOpBase<Scalar> >
00262     A6 = multiply(adjoint(origA),A1);
00263   if(out.get()) *out << "\nA6 =\n" << describe(*A6,verbLevel);
00264 
00265   if(out.get()) *out << "\nTesting A6 ...\n";
00266   Thyra::seed_randomize<Scalar>(0);
00267   result = symLinearOpTester.check(*A6,out.get());
00268   if(!result) success = false;
00269   // Note that testing the symmetry above helps to check the transpose mode
00270   // against the non-transpose mode!
00271 
00272 #ifdef TEUCHOS_DEBUG
00273   if(out.get()) *out << "\nCreating an invalid multiplied operator A6b = origA*origA (should throw an exception) ...\n\n";
00274   try {
00275     RCP<const Thyra::LinearOpBase<Scalar> >
00276       A6b = multiply(origA,origA);
00277     result = true;
00278   }
00279   TEUCHOS_STANDARD_CATCH_STATEMENTS(true,out.get()?*out:std::cerr,result)
00280   if(out.get())
00281     *out << "\nCaught expected exception : " << (result?"failed\n":"passed\n");
00282   if(result) success = false;
00283 #endif // TEUCHOS_DEBUG
00284 
00285   if(out.get()) *out << "\nCreating a non-const multiplied operator A7 = origA^H*A1 ...\n";
00286   RCP<Thyra::LinearOpBase<Scalar> >
00287     A7 = nonconstMultiply(
00288       rcp_const_cast<Thyra::LinearOpBase<Scalar> >(adjoint(origA))
00289       ,rcp_const_cast<Thyra::LinearOpBase<Scalar> >(A1)
00290       );
00291   if(out.get()) *out << "\nA7 =\n" << describe(*A7,verbLevel);
00292 
00293   if(out.get()) *out << "\nTesting A7 ...\n";
00294   Thyra::seed_randomize<Scalar>(0);
00295   result = symLinearOpTester.check(*A7,out.get());
00296   if(!result) success = false;
00297 
00298   if(out.get()) *out << "\nCreating an added operator A8 = origA + A1 ...\n";
00299   RCP<const Thyra::LinearOpBase<Scalar> >
00300     A8 = add(origA,A1);
00301   if(out.get()) *out << "\nA8 =\n" << describe(*A8,verbLevel);
00302 
00303   if(out.get()) *out << "\nTesting A8 ...\n";
00304   Thyra::seed_randomize<Scalar>(0);
00305   result = linearOpTester.check(*A8,out.get());
00306   if(!result) success = false;
00307 
00308   if(out.get()) *out << "\nCreating a symmetric added operator A8b = A6 + adjoint(origA)*origA ...\n";
00309   RCP<const Thyra::LinearOpBase<Scalar> >
00310     A8b = add(A6,multiply(adjoint(origA),origA));
00311   if(out.get()) *out << "\nA8b =\n" << describe(*A8b,verbLevel);
00312 
00313   if(out.get()) *out << "\nTesting A8b ...\n";
00314   Thyra::seed_randomize<Scalar>(0);
00315   result = symLinearOpTester.check(*A8b,out.get());
00316   if(!result) success = false;
00317 
00318 #ifdef TEUCHOS_DEBUG
00319   if(out.get()) *out << "\nCreating an invalid added operator A8c = origA + adjoint(origA) (should throw an exception) ...\n\n";
00320   try {
00321     RCP<const Thyra::LinearOpBase<Scalar> >
00322       A8c = add(origA,adjoint(origA));
00323     result = true;
00324   }
00325   TEUCHOS_STANDARD_CATCH_STATEMENTS(true,out.get()?*out:std::cerr,result)
00326   if(out.get())
00327     *out << "\nCaught expected exception : " << (result?"failed\n":"passed\n");
00328   if(result) success = false;
00329 #endif // TEUCHOS_DEBUG
00330 
00331   RCP<const Thyra::LinearOpBase<Scalar> >
00332     nullOp = null;
00333 
00334   if(out.get()) *out << "\nCreating a blocked 2x2 linear operator A9 = [ A6, A1^H; A1, null ] ...\n";
00335   RCP<const Thyra::LinearOpBase<Scalar> >
00336     A9 = Thyra::block2x2<Scalar>(
00337       A6,  adjoint(A1)
00338       ,A1, nullOp
00339       );
00340   if(out.get()) *out << "\nA9 =\n" << describe(*A9,verbLevel);
00341   
00342   if(out.get()) *out << "\nTesting A9 ...\n";
00343   Thyra::seed_randomize<Scalar>(0);
00344   result = symLinearOpTester.check(*A9,out.get());
00345   if(!result) success = false;
00346   // Note that testing the symmetry above helps to check the transpose mode
00347   // against the non-transpose mode!
00348 
00349   if(out.get()) *out << "\nCreating a blocked 2x2 linear operator A9_a = [ A6, A1^H; A1, null ] using pre-formed range and domain product spaces ...\n";
00350   RCP<Thyra::PhysicallyBlockedLinearOpBase<Scalar> >
00351     A9_a = rcp(new Thyra::DefaultBlockedLinearOp<Scalar>());
00352   A9_a->beginBlockFill(
00353     rcp_dynamic_cast<const Thyra::BlockedLinearOpBase<Scalar> >(A9,true)->productRange()
00354     ,rcp_dynamic_cast<const Thyra::BlockedLinearOpBase<Scalar> >(A9,true)->productDomain()
00355     );
00356   A9_a->setBlock(0,0,A6);
00357   A9_a->setBlock(0,1,adjoint(A1));
00358   A9_a->setBlock(1,0,A1);
00359   A9_a->endBlockFill();
00360   if(out.get()) *out << "\nA9_a =\n" << describe(*A9_a,verbLevel);
00361   
00362   if(out.get()) *out << "\nTesting A9_a ...\n";
00363   Thyra::seed_randomize<Scalar>(0);
00364   result = symLinearOpTester.check(*A9_a,out.get());
00365   if(!result) success = false;
00366   // Note that testing the symmetry above helps to check the transpose mode
00367   // against the non-transpose mode!
00368   
00369   if(out.get()) *out << "\nComparing A9 == A9_a ...\n";
00370   Thyra::seed_randomize<Scalar>(0);
00371   result = linearOpTester.compare(*A9,*A9_a,out.get());
00372   if(!result) success = false;
00373 
00374   if(out.get()) *out << "\nCreating a blocked 2x2 linear operator A9_b = [ A6, A1^H; A1, null ] using flexible fill ...\n";
00375   RCP<Thyra::PhysicallyBlockedLinearOpBase<Scalar> >
00376     A9_b = rcp(new Thyra::DefaultBlockedLinearOp<Scalar>());
00377   A9_b->beginBlockFill();
00378   A9_b->setBlock(0,0,A6);
00379   A9_b->setBlock(0,1,adjoint(A1));
00380   A9_b->setBlock(1,0,A1);
00381   A9_b->endBlockFill();
00382   if(out.get()) *out << "\nA9_b =\n" << describe(*A9_b,verbLevel);
00383   
00384   if(out.get()) *out << "\nTesting A9_b ...\n";
00385   Thyra::seed_randomize<Scalar>(0);
00386   result = symLinearOpTester.check(*A9_b,out.get());
00387   if(!result) success = false;
00388   // Note that testing the symmetry above helps to check the transpose mode
00389   // against the non-transpose mode!
00390   
00391   if(out.get()) *out << "\nComparing A9 == A9_b ...\n";
00392   Thyra::seed_randomize<Scalar>(0);
00393   result = linearOpTester.compare(*A9,*A9_b,out.get());
00394   if(!result) success = false;
00395 
00396   if(out.get()) *out << "\nCreating a blocked 2x2 linear operator A9a = [ null, A1^H; A1, null ] ...\n";
00397   RCP<const Thyra::LinearOpBase<Scalar> >
00398     A9a = Thyra::block2x2<Scalar>(
00399       nullOp,  adjoint(A1),
00400       A1,      nullOp
00401       );
00402   if(out.get()) *out << "\nA9a =\n" << describe(*A9a,verbLevel);
00403   
00404   if(out.get()) *out << "\nTesting A9a ...\n";
00405   Thyra::seed_randomize<Scalar>(0);
00406   result = symLinearOpTester.check(*A9a,out.get());
00407   if(!result) success = false;
00408   // Note that testing the symmetry above helps to check the transpose mode
00409   // against the non-transpose mode!
00410   
00411 #ifdef TEUCHOS_DEBUG
00412   if(out.get()) *out << "\nCreating an invalid blocked 2x2 operator A9b = [ A6, A1^H; A1, A1 ] (should throw an exception) ...\n\n";
00413   try {
00414     RCP<const Thyra::LinearOpBase<Scalar> >
00415       A9b = Thyra::block2x2<Scalar>(
00416         A6,  adjoint(A1),
00417         A1,  A1
00418         );
00419     result = true;
00420   }
00421   TEUCHOS_STANDARD_CATCH_STATEMENTS(true,out.get()?*out:std::cerr,result)
00422   if(out.get())
00423     *out << "\nCaught expected exception : " << (result?"failed\n":"passed\n");
00424   if(result) success = false;
00425 #endif // TEUCHOS_DEBUG
00426 
00427 #ifdef TEUCHOS_DEBUG
00428   if(out.get()) *out << "\nCreating an invalid blocked 2x2 operator A9c = [ A1, A1 ; null, null ] (should throw an exception) ...\n\n";
00429   try {
00430     RCP<const Thyra::LinearOpBase<Scalar> >
00431       A9c = Thyra::block2x2<Scalar>(
00432         A1,       A1,
00433         nullOp,   nullOp
00434         );
00435     result = true;
00436   }
00437   TEUCHOS_STANDARD_CATCH_STATEMENTS(true,out.get()?*out:std::cerr,result)
00438   if(out.get())
00439     *out << "\nCaught expected exception : " << (result?"failed\n":"passed\n");
00440   if(result) success = false;
00441 #endif // TEUCHOS_DEBUG
00442 
00443 #ifdef TEUCHOS_DEBUG
00444   if(out.get()) *out << "\nCreating an invalid blocked 2x2 operator A9d = [ A1, null; A1, null ] (should throw an exception) ...\n\n";
00445   try {
00446     RCP<const Thyra::LinearOpBase<Scalar> >
00447       A9d = Thyra::block2x2<Scalar>(
00448         A1,  nullOp,
00449         A1,  nullOp
00450         );
00451     result = true;
00452   }
00453   TEUCHOS_STANDARD_CATCH_STATEMENTS(true,out.get()?*out:std::cerr,result)
00454   if(out.get())
00455     *out << "\nCaught expected exception : " << (result?"failed\n":"passed\n");
00456   if(result) success = false;
00457 #endif // TEUCHOS_DEBUG
00458 
00459   if(out.get()) *out << "\nCreating a blocked 2x1 linear operator A10 = [ A6; A1 ] ...\n";
00460   RCP<const Thyra::LinearOpBase<Scalar> >
00461     A10 = Thyra::block2x1<Scalar>(
00462       A6,
00463       A1
00464       );
00465   if(out.get()) *out << "\nA10 =\n" << describe(*A10,verbLevel);
00466   
00467   if(out.get()) *out << "\nTesting A10 ...\n";
00468   Thyra::seed_randomize<Scalar>(0);
00469   result = linearOpTester.check(*A10,out.get());
00470   if(!result) success = false;
00471 
00472   if(out.get()) *out << "\nCreating a blocked 1x2 linear operator A11 = [ A9, A10 ] ...\n";
00473   RCP<const Thyra::LinearOpBase<Scalar> >
00474     A11 = Thyra::block1x2<Scalar>( A9, A10 );
00475   if(out.get()) *out << "\nA11 =\n" << describe(*A11,verbLevel);
00476   
00477   if(out.get()) *out << "\nTesting A11 ...\n";
00478   Thyra::seed_randomize<Scalar>(0);
00479   result = linearOpTester.check(*A11,out.get());
00480   if(!result) success = false;
00481 
00482   if(out.get()) *out << "\nCreating a zero linear operator A12 = 0 (range and domain spaces of origA) ...\n";
00483   RCP<const Thyra::LinearOpBase<Scalar> >
00484     A12 = Thyra::zero(origA->range(),origA->domain());
00485   if(out.get()) *out << "\nA12 =\n" << describe(*A12,verbLevel);
00486 
00487   if(out.get()) *out << "\nTesting A12 ...\n";
00488   Thyra::seed_randomize<Scalar>(0);
00489   result = linearOpTester.check(*A12,out.get());
00490   if(!result) success = false;
00491 
00492   if(out.get()) *out << "\nCreating a blocked 2x2 linear operator A13 = [ zero, A1^H; A1, zero ] ...\n";
00493   RCP<const Thyra::LinearOpBase<Scalar> >
00494     A13 = Thyra::block2x2<Scalar>(
00495       Thyra::zero(A1->domain(),A1->domain()),   adjoint(A1),
00496       A1,                                       Thyra::zero(A1->range(),A1->range())
00497       );
00498   if(out.get()) *out << "\nA13 =\n" << describe(*A13,verbLevel);
00499   
00500   if(out.get()) *out << "\nComparing A9a == A13 ...\n";
00501   Thyra::seed_randomize<Scalar>(0);
00502   result = linearOpTester.compare(*A9a,*A13,out.get());
00503   if(!result) success = false;
00504 
00505   if(out.get()) *out << "\nCreating a zero linear operator A14 = I (range space of origA) ...\n";
00506   RCP<const Thyra::LinearOpBase<Scalar> >
00507     A14 = Thyra::identity(origA->range());
00508   if(out.get()) *out << "\nA14 =\n" << describe(*A14,verbLevel);
00509   
00510   if(out.get()) *out << "\nTesting A14 ...\n";
00511   Thyra::seed_randomize<Scalar>(0);
00512   result = symLinearOpTester.check(*A14,out.get());
00513   if(!result) success = false;
00514   
00515   if(out.get()) *out << "\n*** Leaving run_composite_linear_ops_tests<"<<ST::name()<<">(...) ...\n";
00516 
00517   return success;
00518 
00519 } // end run_composite_linear_ops_tests() [Doxygen looks for this!]
00520 
00521 int main( int argc, char* argv[] ) {
00522 
00523   using Teuchos::CommandLineProcessor;
00524 
00525   bool success = true;
00526   bool verbose = true;
00527 
00528   Teuchos::GlobalMPISession mpiSession(&argc,&argv);
00529 
00530   const Teuchos::RCP<const Teuchos::Comm<Thyra::Index> >
00531     comm = Teuchos::DefaultComm<Thyra::Index>::getComm();
00532 
00533   Teuchos::RCP<Teuchos::FancyOStream>
00534     out = Teuchos::VerboseObjectBase::getDefaultOStream();
00535   
00536   try {
00537 
00538     //
00539     // Read options from command-line
00540     //
00541 
00542     int n         = 4;
00543     bool useSpmd   = true;
00544     bool dumpAll  = false;
00545 
00546     CommandLineProcessor  clp;
00547     clp.throwExceptions(false);
00548     clp.addOutputSetupOptions(true);
00549     clp.setOption( "verbose", "quiet", &verbose, "Set if output is printed or not." );
00550     clp.setOption( "local-dim", &n, "Local number of elements in each constituent vector." );
00551     clp.setOption( "use-spmd", "use-serial", &useSpmd, "Determines if MPI or serial vector space is used." );
00552     clp.setOption( "dump-all", "no-dump-all", &dumpAll, "Determines if vectors are printed or not." );
00553     CommandLineProcessor::EParseCommandLineReturn parse_return = clp.parse(argc,argv);
00554     if( parse_return != CommandLineProcessor::PARSE_SUCCESSFUL ) return parse_return;
00555 
00556     //
00557     // Run the tests
00558     //
00559 
00560 #ifdef HAVE_THYRA_TEUCHOS_BLASFLOAT
00561     if( !run_composite_linear_ops_tests<float>(comm,n,useSpmd,float(1e-4),dumpAll,verbose?&*out:NULL) ) success = false;
00562 #endif // HAVE_THYRA_TEUCHOS_BLASFLOAT
00563     if( !run_composite_linear_ops_tests<double>(comm,n,useSpmd,double(1e-12),dumpAll,verbose?&*out:NULL) ) success = false;
00564 #if defined(HAVE_COMPLEX) && defined(HAVE_TEUCHOS_COMPLEX)
00565 #ifdef HAVE_THYRA_TEUCHOS_BLASFLOAT
00566     if( !run_composite_linear_ops_tests<std::complex<float> >(comm,n,useSpmd,float(1e-4),dumpAll,verbose?&*out:NULL) ) success = false;
00567 #endif // HAVE_THYRA_TEUCHOS_BLASFLOAT
00568     if( !run_composite_linear_ops_tests<std::complex<double> >(comm,n,useSpmd,double(1e-12),dumpAll,verbose?&*out:NULL) ) success = false;
00569 #endif
00570 #if defined(HAVE_TEUCHOS_GNU_MP) && !defined(USE_MPI) // mpf_class can not be used with MPI yet!
00571     if( !run_composite_linear_ops_tests<mpf_class>(comm,n,useSpmd,mpf_class(1e-12),dumpAll,verbose?&*out:NULL) ) success = false;
00572 #endif
00573 
00574   } // end try
00575   TEUCHOS_STANDARD_CATCH_STATEMENTS(true,*out,success)
00576 
00577   if( verbose ) {
00578     if(success) *out << "\nAll of the tests seem to have run successfully!\n";
00579     else        *out << "\nOh no! at least one of the tests failed!\n"; 
00580   }
00581   
00582   return success ? 0 : 1;
00583 
00584 } // end main() [Doxygen looks for this!]

Generated on Tue Oct 20 12:46:58 2009 for Thyra Operator/Vector Support by doxygen 1.4.7