test_composite_linear_ops.cpp

Go to the documentation of this file.
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::RefCountPtr<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::RefCountPtr;
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   RefCountPtr<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   RefCountPtr<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   RefCountPtr<Thyra::MultiVectorBase<Scalar> >
00103     mvOrigA = createMembers(space,n/2);
00104   Thyra::seed_randomize<Scalar>(0);
00105   Thyra::randomize( Scalar(Scalar(-1)*ST::one()), Scalar(Scalar(+1)*ST::one()), &*mvOrigA );
00106   RefCountPtr<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   RefCountPtr<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   if(1) {
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   if(1) {
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   RefCountPtr<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   RefCountPtr<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   RefCountPtr<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   RefCountPtr<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   RefCountPtr<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     RefCountPtr<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   RefCountPtr<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   RefCountPtr<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   RefCountPtr<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     RefCountPtr<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   if(out.get()) *out << "\nCreating a blocked 2x2 linear operator A9 = [ A6, A1^H; A1, null ] ...\n";
00332   RefCountPtr<const Thyra::LinearOpBase<Scalar> >
00333     A9 = Thyra::block2x2<Scalar>(
00334       A6,  adjoint(A1)
00335       ,A1, null
00336       );
00337   if(out.get()) *out << "\nA9 =\n" << describe(*A9,verbLevel);
00338   
00339   if(out.get()) *out << "\nTesting A9 ...\n";
00340   Thyra::seed_randomize<Scalar>(0);
00341   result = symLinearOpTester.check(*A9,out.get());
00342   if(!result) success = false;
00343   // Note that testing the symmetry above helps to check the transpose mode
00344   // against the non-transpose mode!
00345 
00346   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";
00347   RefCountPtr<Thyra::PhysicallyBlockedLinearOpBase<Scalar> >
00348     A9_a = rcp(new Thyra::DefaultBlockedLinearOp<Scalar>());
00349   A9_a->beginBlockFill(
00350     rcp_dynamic_cast<const Thyra::BlockedLinearOpBase<Scalar> >(A9,true)->productRange()
00351     ,rcp_dynamic_cast<const Thyra::BlockedLinearOpBase<Scalar> >(A9,true)->productDomain()
00352     );
00353   A9_a->setBlock(0,0,A6);
00354   A9_a->setBlock(0,1,adjoint(A1));
00355   A9_a->setBlock(1,0,A1);
00356   A9_a->endBlockFill();
00357   if(out.get()) *out << "\nA9_a =\n" << describe(*A9_a,verbLevel);
00358   
00359   if(out.get()) *out << "\nTesting A9_a ...\n";
00360   Thyra::seed_randomize<Scalar>(0);
00361   result = symLinearOpTester.check(*A9_a,out.get());
00362   if(!result) success = false;
00363   // Note that testing the symmetry above helps to check the transpose mode
00364   // against the non-transpose mode!
00365   
00366   if(out.get()) *out << "\nComparing A9 == A9_a ...\n";
00367   Thyra::seed_randomize<Scalar>(0);
00368   result = linearOpTester.compare(*A9,*A9_a,out.get());
00369   if(!result) success = false;
00370 
00371   if(out.get()) *out << "\nCreating a blocked 2x2 linear operator A9_b = [ A6, A1^H; A1, null ] using flexible fill ...\n";
00372   RefCountPtr<Thyra::PhysicallyBlockedLinearOpBase<Scalar> >
00373     A9_b = rcp(new Thyra::DefaultBlockedLinearOp<Scalar>());
00374   A9_b->beginBlockFill();
00375   A9_b->setBlock(0,0,A6);
00376   A9_b->setBlock(0,1,adjoint(A1));
00377   A9_b->setBlock(1,0,A1);
00378   A9_b->endBlockFill();
00379   if(out.get()) *out << "\nA9_b =\n" << describe(*A9_b,verbLevel);
00380   
00381   if(out.get()) *out << "\nTesting A9_b ...\n";
00382   Thyra::seed_randomize<Scalar>(0);
00383   result = symLinearOpTester.check(*A9_b,out.get());
00384   if(!result) success = false;
00385   // Note that testing the symmetry above helps to check the transpose mode
00386   // against the non-transpose mode!
00387   
00388   if(out.get()) *out << "\nComparing A9 == A9_b ...\n";
00389   Thyra::seed_randomize<Scalar>(0);
00390   result = linearOpTester.compare(*A9,*A9_b,out.get());
00391   if(!result) success = false;
00392 
00393   if(out.get()) *out << "\nCreating a blocked 2x2 linear operator A9a = [ null, A1^H; A1, null ] ...\n";
00394   RefCountPtr<const Thyra::LinearOpBase<Scalar> >
00395     A9a = Thyra::block2x2<Scalar>(
00396       null,  adjoint(A1)
00397       ,A1,   null
00398       );
00399   if(out.get()) *out << "\nA9a =\n" << describe(*A9a,verbLevel);
00400   
00401   if(out.get()) *out << "\nTesting A9a ...\n";
00402   Thyra::seed_randomize<Scalar>(0);
00403   result = symLinearOpTester.check(*A9a,out.get());
00404   if(!result) success = false;
00405   // Note that testing the symmetry above helps to check the transpose mode
00406   // against the non-transpose mode!
00407   
00408 #ifdef TEUCHOS_DEBUG
00409   if(out.get()) *out << "\nCreating an invalid blocked 2x2 operator A9b = [ A6, A1^H; A1, A1 ] (should throw an exception) ...\n\n";
00410   try {
00411     RefCountPtr<const Thyra::LinearOpBase<Scalar> >
00412       A9b = Thyra::block2x2<Scalar>(
00413         A6,  adjoint(A1)
00414         ,A1, A1
00415         );
00416     result = true;
00417   }
00418   TEUCHOS_STANDARD_CATCH_STATEMENTS(true,out.get()?*out:std::cerr,result)
00419   if(out.get())
00420     *out << "\nCaught expected exception : " << (result?"failed\n":"passed\n");
00421   if(result) success = false;
00422 #endif // TEUCHOS_DEBUG
00423 
00424 #ifdef TEUCHOS_DEBUG
00425   if(out.get()) *out << "\nCreating an invalid blocked 2x2 operator A9c = [ A1, A1 ; null, null ] (should throw an exception) ...\n\n";
00426   try {
00427     RefCountPtr<const Thyra::LinearOpBase<Scalar> >
00428       A9c = Thyra::block2x2<Scalar>(
00429         A1,    A1
00430         ,null, null
00431         );
00432     result = true;
00433   }
00434   TEUCHOS_STANDARD_CATCH_STATEMENTS(true,out.get()?*out:std::cerr,result)
00435   if(out.get())
00436     *out << "\nCaught expected exception : " << (result?"failed\n":"passed\n");
00437   if(result) success = false;
00438 #endif // TEUCHOS_DEBUG
00439 
00440 #ifdef TEUCHOS_DEBUG
00441   if(out.get()) *out << "\nCreating an invalid blocked 2x2 operator A9d = [ A1, null; A1, null ] (should throw an exception) ...\n\n";
00442   try {
00443     RefCountPtr<const Thyra::LinearOpBase<Scalar> >
00444       A9d = Thyra::block2x2<Scalar>(
00445         A1,  null
00446         ,A1, null
00447         );
00448     result = true;
00449   }
00450   TEUCHOS_STANDARD_CATCH_STATEMENTS(true,out.get()?*out:std::cerr,result)
00451   if(out.get())
00452     *out << "\nCaught expected exception : " << (result?"failed\n":"passed\n");
00453   if(result) success = false;
00454 #endif // TEUCHOS_DEBUG
00455 
00456   if(out.get()) *out << "\nCreating a blocked 2x1 linear operator A10 = [ A6; A1 ] ...\n";
00457   RefCountPtr<const Thyra::LinearOpBase<Scalar> >
00458     A10 = Thyra::block2x1<Scalar>(
00459       A6
00460       ,A1
00461       );
00462   if(out.get()) *out << "\nA10 =\n" << describe(*A10,verbLevel);
00463   
00464   if(out.get()) *out << "\nTesting A10 ...\n";
00465   Thyra::seed_randomize<Scalar>(0);
00466   result = linearOpTester.check(*A10,out.get());
00467   if(!result) success = false;
00468 
00469   if(out.get()) *out << "\nCreating a blocked 1x2 linear operator A11 = [ A9, A10 ] ...\n";
00470   RefCountPtr<const Thyra::LinearOpBase<Scalar> >
00471     A11 = Thyra::block1x2<Scalar>( A9, A10 );
00472   if(out.get()) *out << "\nA11 =\n" << describe(*A11,verbLevel);
00473   
00474   if(out.get()) *out << "\nTesting A11 ...\n";
00475   Thyra::seed_randomize<Scalar>(0);
00476   result = linearOpTester.check(*A11,out.get());
00477   if(!result) success = false;
00478 
00479   if(out.get()) *out << "\nCreating a zero linear operator A12 = 0 (range and domain spaces of origA) ...\n";
00480   RefCountPtr<const Thyra::LinearOpBase<Scalar> >
00481     A12 = Thyra::zero(origA->range(),origA->domain());
00482   if(out.get()) *out << "\nA12 =\n" << describe(*A12,verbLevel);
00483 
00484   if(out.get()) *out << "\nTesting A12 ...\n";
00485   Thyra::seed_randomize<Scalar>(0);
00486   result = linearOpTester.check(*A12,out.get());
00487   if(!result) success = false;
00488 
00489   if(out.get()) *out << "\nCreating a blocked 2x2 linear operator A13 = [ zero, A1^H; A1, zero ] ...\n";
00490   RefCountPtr<const Thyra::LinearOpBase<Scalar> >
00491     A13 = Thyra::block2x2<Scalar>(
00492       Thyra::zero(A1->domain(),A1->domain()),  adjoint(A1)
00493       ,A1,                                     Thyra::zero(A1->range(),A1->range())
00494       );
00495   if(out.get()) *out << "\nA13 =\n" << describe(*A13,verbLevel);
00496   
00497   if(out.get()) *out << "\nComparing A9a == A13 ...\n";
00498   Thyra::seed_randomize<Scalar>(0);
00499   result = linearOpTester.compare(*A9a,*A13,out.get());
00500   if(!result) success = false;
00501 
00502   if(out.get()) *out << "\nCreating a zero linear operator A14 = I (range space of origA) ...\n";
00503   RefCountPtr<const Thyra::LinearOpBase<Scalar> >
00504     A14 = Thyra::identity(origA->range());
00505   if(out.get()) *out << "\nA14 =\n" << describe(*A14,verbLevel);
00506   
00507   if(out.get()) *out << "\nTesting A14 ...\n";
00508   Thyra::seed_randomize<Scalar>(0);
00509   result = symLinearOpTester.check(*A14,out.get());
00510   if(!result) success = false;
00511   
00512   if(out.get()) *out << "\n*** Leaving run_composite_linear_ops_tests<"<<ST::name()<<">(...) ...\n";
00513 
00514   return success;
00515 
00516 } // end run_composite_linear_ops_tests() [Doxygen looks for this!]
00517 
00518 int main( int argc, char* argv[] ) {
00519 
00520   using Teuchos::CommandLineProcessor;
00521 
00522   bool success = true;
00523   bool verbose = true;
00524 
00525   Teuchos::GlobalMPISession mpiSession(&argc,&argv);
00526 
00527   const Teuchos::RefCountPtr<const Teuchos::Comm<Thyra::Index> >
00528     comm = Teuchos::DefaultComm<Thyra::Index>::getComm();
00529 
00530   Teuchos::RefCountPtr<Teuchos::FancyOStream>
00531     out = Teuchos::VerboseObjectBase::getDefaultOStream();
00532   
00533   try {
00534 
00535     //
00536     // Read options from command-line
00537     //
00538 
00539     int n         = 4;
00540     bool useSpmd   = true;
00541     bool dumpAll  = false;
00542 
00543     CommandLineProcessor  clp;
00544     clp.throwExceptions(false);
00545     clp.addOutputSetupOptions(true);
00546     clp.setOption( "verbose", "quiet", &verbose, "Set if output is printed or not." );
00547     clp.setOption( "local-dim", &n, "Local number of elements in each constituent vector." );
00548     clp.setOption( "use-spmd", "use-serial", &useSpmd, "Determines if MPI or serial vector space is used." );
00549     clp.setOption( "dump-all", "no-dump-all", &dumpAll, "Determines if vectors are printed or not." );
00550     CommandLineProcessor::EParseCommandLineReturn parse_return = clp.parse(argc,argv);
00551     if( parse_return != CommandLineProcessor::PARSE_SUCCESSFUL ) return parse_return;
00552 
00553     //
00554     // Run the tests
00555     //
00556 
00557     if( !run_composite_linear_ops_tests<float>(comm,n,useSpmd,float(1e-4),dumpAll,verbose?&*out:NULL) ) success = false;
00558     if( !run_composite_linear_ops_tests<double>(comm,n,useSpmd,double(1e-12),dumpAll,verbose?&*out:NULL) ) success = false;
00559 #if defined(HAVE_COMPLEX) && defined(HAVE_TEUCHOS_COMPLEX)
00560     if( !run_composite_linear_ops_tests<std::complex<float> >(comm,n,useSpmd,float(1e-4),dumpAll,verbose?&*out:NULL) ) success = false;
00561     if( !run_composite_linear_ops_tests<std::complex<double> >(comm,n,useSpmd,double(1e-12),dumpAll,verbose?&*out:NULL) ) success = false;
00562 #endif
00563 #if defined(HAVE_TEUCHOS_GNU_MP) && !defined(USE_MPI) // mpf_class can not be used with MPI yet!
00564     if( !run_composite_linear_ops_tests<mpf_class>(comm,n,useSpmd,mpf_class(1e-12),dumpAll,verbose?&*out:NULL) ) success = false;
00565 #endif
00566 
00567   } // end try
00568   TEUCHOS_STANDARD_CATCH_STATEMENTS(true,*out,success)
00569 
00570   if( verbose ) {
00571     if(success) *out << "\nAll of the tests seem to have run successfully!\n";
00572     else        *out << "\nOh no! at least one of the tests failed!\n"; 
00573   }
00574   
00575   return success ? 0 : 1;
00576 
00577 } // end main() [Doxygen looks for this!]

Generated on Thu Sep 18 12:33:02 2008 for Thyra Package Browser (Single Doxygen Collection) by doxygen 1.3.9.1