Stratimikos Package Browser (Single Doxygen Collection) Version of the Day
test_single_aztecoo_thyra_solver.cpp
Go to the documentation of this file.
00001 /*@HEADER
00002 // ***********************************************************************
00003 // 
00004 //        AztecOO: An Object-Oriented Aztec Linear Solver Package 
00005 //                 Copyright (2002) 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 "test_single_aztecoo_thyra_solver.hpp"
00030 
00031 #ifndef SUN_CXX
00032 
00033 #include "Thyra_AztecOOLinearOpWithSolveFactory.hpp"
00034 #include "Thyra_LinearOpWithSolveFactoryHelpers.hpp"
00035 #include "Thyra_PreconditionerFactoryHelpers.hpp"
00036 #include "Thyra_LinearOpWithSolveFactoryExamples.hpp"
00037 #include "Thyra_DefaultScaledAdjointLinearOp.hpp"
00038 #include "Thyra_EpetraLinearOp.hpp"
00039 #include "Thyra_LinearOpTester.hpp"
00040 #include "Thyra_LinearOpWithSolveBase.hpp"
00041 #include "Thyra_LinearOpWithSolveTester.hpp"
00042 #include "EpetraExt_readEpetraLinearSystem.h"
00043 #include "Epetra_SerialComm.h"
00044 #include "Teuchos_ParameterList.hpp"
00045 #ifdef HAVE_AZTECOO_IFPACK
00046 #  include "Thyra_IfpackPreconditionerFactory.hpp"
00047 #endif
00048 #ifdef HAVE_MPI
00049 #  include "Epetra_MpiComm.h"
00050 #endif
00051 
00052 #endif // SUN_CXX
00053 
00054 bool Thyra::test_single_aztecoo_thyra_solver(
00055   const std::string                       matrixFile
00056   ,const bool                             testTranspose
00057   ,const int                              numRandomVectors
00058   ,const double                           maxFwdError
00059   ,const double                           maxResid
00060   ,const double                           maxSolutionError
00061   ,const bool                             showAllTests
00062   ,const bool                             dumpAll
00063   ,Teuchos::ParameterList                 *aztecooLOWSFPL
00064   ,Teuchos::FancyOStream                  *out_arg
00065   )
00066 {
00067   using Teuchos::rcp;
00068   using Teuchos::RCP;
00069   using Teuchos::OSTab;
00070   typedef Teuchos::ParameterList::PrintOptions PLPrintOptions;
00071   bool result, success = true;
00072 
00073   RCP<Teuchos::FancyOStream>
00074     out = Teuchos::rcp(out_arg,false);
00075 
00076   try {
00077 
00078 #ifndef SUN_CXX
00079 
00080     if(out.get()) {
00081       *out
00082         << "\n***"
00083         << "\n*** Testing Thyra::AztecOOLinearOpWithSolveFactory (and Thyra::AztecOOLinearOpWithSolve)"
00084         << "\n***\n"
00085         << "\nEchoing input options:"
00086         << "\n  matrixFile             = " << matrixFile
00087         << "\n  testTranspose          = " << testTranspose
00088         << "\n  numRandomVectors       = " << numRandomVectors
00089         << "\n  maxFwdError            = " << maxFwdError
00090         << "\n  maxResid               = " << maxResid
00091         << "\n  showAllTests           = " << showAllTests
00092         << "\n  dumpAll                = " << dumpAll
00093         << std::endl;
00094     }
00095     
00096     const bool useAztecPrec = ( 
00097       aztecooLOWSFPL
00098       &&
00099       aztecooLOWSFPL->sublist("Forward Solve")
00100       .sublist("AztecOO Settings")
00101       .get("Aztec Preconditioner","none")!="none"
00102       );
00103 
00104     if(out.get()) {
00105       if(useAztecPrec)
00106         *out << "\nUsing aztec preconditioning so we will not test adjoint solves using internal preconditioning ...\n";
00107     }
00108     
00109     if(out.get()) *out << "\nA) Reading in an epetra matrix A from the file \'"<<matrixFile<<"\' ...\n";
00110   
00111 #ifdef HAVE_MPI
00112     Epetra_MpiComm comm(MPI_COMM_WORLD);
00113 #else
00114     Epetra_SerialComm comm;
00115 #endif
00116     RCP<Epetra_CrsMatrix> epetra_A;
00117     EpetraExt::readEpetraLinearSystem( matrixFile, comm, &epetra_A );
00118 
00119     RCP<const LinearOpBase<double> > A = Thyra::epetraLinearOp(epetra_A);
00120 
00121     if(out.get() && dumpAll) *out << "\ndescribe(A) =\n" << describe(*A,Teuchos::VERB_EXTREME);
00122 
00123     if(out.get()) *out << "\nB) Creating a AztecOOLinearOpWithSolveFactory object opFactory ...\n";
00124 
00125     RCP<LinearOpWithSolveFactoryBase<double> >
00126       lowsFactory = Teuchos::rcp(new AztecOOLinearOpWithSolveFactory());
00127     if(out.get()) {
00128       *out << "\nlowsFactory.getValidParameters() initially:\n";
00129       lowsFactory->getValidParameters()->print(OSTab(out).o(),PLPrintOptions().showTypes(true).showDoc(true));
00130     }
00131     aztecooLOWSFPL->sublist("Forward Solve").set("Tolerance",maxResid);
00132     aztecooLOWSFPL->sublist("Adjoint Solve").set("Tolerance",maxResid);
00133     if(showAllTests) {
00134       aztecooLOWSFPL->set("Output Every RHS",bool(true));
00135     }
00136     if(out.get()) {
00137       *out << "\naztecooLOWSFPL before setting parameters:\n";
00138       aztecooLOWSFPL->print(OSTab(out).o(),0,true);
00139     }
00140     if(aztecooLOWSFPL) lowsFactory->setParameterList(Teuchos::rcp(aztecooLOWSFPL,false));
00141     
00142     if(out.get()) *out << "\nC) Creating a AztecOOLinearOpWithSolve object nsA from A ...\n";
00143 
00144     RCP<LinearOpWithSolveBase<double> >
00145       nsA = lowsFactory->createOp();
00146 
00147     Thyra::initializeOp<double>(*lowsFactory, A, nsA.ptr());
00148 
00149     if(out.get()) *out << "\nD) Testing the LinearOpBase interface of nsA ...\n";
00150 
00151     LinearOpTester<double> linearOpTester;
00152     linearOpTester.check_adjoint(testTranspose);
00153     linearOpTester.num_random_vectors(numRandomVectors);
00154     linearOpTester.set_all_error_tol(maxFwdError);
00155     linearOpTester.set_all_warning_tol(1e-2*maxFwdError);
00156     linearOpTester.show_all_tests(showAllTests);
00157     linearOpTester.dump_all(dumpAll);
00158     Thyra::seed_randomize<double>(0);
00159     result = linearOpTester.check(*nsA,out());
00160     if(!result) success = false;
00161 
00162     if(out.get()) *out << "\nE) Testing the LinearOpWithSolveBase interface of nsA ...\n";
00163     
00164     LinearOpWithSolveTester<double> linearOpWithSolveTester;
00165     linearOpWithSolveTester.turn_off_all_tests();
00166     linearOpWithSolveTester.check_forward_default(true);
00167     linearOpWithSolveTester.check_forward_residual(true);
00168     if(testTranspose && useAztecPrec) {
00169       linearOpWithSolveTester.check_adjoint_default(true);
00170       linearOpWithSolveTester.check_adjoint_residual(true);
00171     }
00172     else {
00173       linearOpWithSolveTester.check_adjoint_default(false);
00174       linearOpWithSolveTester.check_adjoint_residual(false);
00175     }
00176     linearOpWithSolveTester.set_all_solve_tol(maxResid);
00177     linearOpWithSolveTester.set_all_slack_error_tol(maxResid);
00178     linearOpWithSolveTester.set_all_slack_warning_tol(10.0*maxResid);
00179     linearOpWithSolveTester.forward_default_residual_error_tol(2.5*maxResid);
00180     linearOpWithSolveTester.forward_default_solution_error_error_tol(maxSolutionError);
00181     linearOpWithSolveTester.adjoint_default_residual_error_tol(2.5*maxResid);
00182     linearOpWithSolveTester.adjoint_default_solution_error_error_tol(maxSolutionError);
00183     linearOpWithSolveTester.show_all_tests(showAllTests);
00184     linearOpWithSolveTester.dump_all(dumpAll);
00185     Thyra::seed_randomize<double>(0);
00186     result = linearOpWithSolveTester.check(*nsA,out.get());
00187     if(!result) success = false;
00188 
00189     if(out.get()) *out << "\nF) Uninitialize nsA, create preconditioner for diagonal scaled by 0.99 and then reinitialize nsA reusing the old preconditioner ...\n";
00190 
00191     // Scale the diagonal of the matrix and then create the preconditioner for it
00192     Thyra::uninitializeOp<double>(*lowsFactory, nsA.ptr());
00193     // Above is not required but a good idea since we are changing the matrix
00194     {
00195       Epetra_Vector diag(epetra_A->RowMap());
00196       epetra_A->ExtractDiagonalCopy(diag);
00197       diag.Scale(0.5);
00198       epetra_A->ReplaceDiagonalValues(diag);
00199     }
00200     Thyra::initializeOp<double>(*lowsFactory, A, nsA.ptr());
00201 
00202     // Scale the matrix back again and then reuse the preconditioner
00203     Thyra::uninitializeOp<double>(*lowsFactory, nsA.ptr());
00204     // Above is not required but a good idea since we are changing the matrix
00205     {
00206       Epetra_Vector diag(epetra_A->RowMap());
00207       epetra_A->ExtractDiagonalCopy(diag);
00208       diag.Scale(1.0/0.5);
00209       epetra_A->ReplaceDiagonalValues(diag);
00210     }
00211     initializeAndReuseOp<double>(*lowsFactory, A, nsA.ptr());
00212 
00213     if(out.get()) *out << "\nG) Testing the LinearOpWithSolveBase interface of nsA ...\n";
00214     
00215     Thyra::seed_randomize<double>(0);
00216     result = linearOpWithSolveTester.check(*nsA,out.get());
00217     if(!result) success = false;
00218 
00219     if(useAztecPrec) {
00220 
00221       if(out.get()) *out << "\nH) Reinitialize (A,A,PRECONDITIONER_INPUT_TYPE_AS_MATRIX) => nsA ...\n";
00222       
00223       initializeApproxPreconditionedOp<double>(*lowsFactory, A, A, nsA.ptr());
00224 
00225       if(out.get()) *out << "\nI) Testing the LinearOpWithSolveBase interface of nsA ...\n";
00226       
00227       Thyra::seed_randomize<double>(0);
00228       result = linearOpWithSolveTester.check(*nsA,out.get());
00229       if(!result) success = false;
00230 
00231       if(testTranspose && useAztecPrec) {
00232         linearOpWithSolveTester.check_adjoint_default(true);
00233         linearOpWithSolveTester.check_adjoint_residual(true);
00234       }
00235       else {
00236         linearOpWithSolveTester.check_adjoint_default(false);
00237         linearOpWithSolveTester.check_adjoint_residual(false);
00238       }
00239       
00240     }
00241     else {
00242 
00243       if(out.get()) *out << "\nSkipping testing steps H and I since we are not using aztec preconditioning and therefore will not test with an external preconditioner matrix!\n";
00244 
00245     }
00246 
00247     
00248     RCP<PreconditionerFactoryBase<double> >
00249       precFactory;
00250 
00251 #ifdef HAVE_AZTECOO_IFPACK
00252 
00253     if(useAztecPrec) {
00254 
00255       if(testTranspose) {
00256         linearOpWithSolveTester.check_adjoint_default(true);
00257         linearOpWithSolveTester.check_adjoint_residual(true);
00258       }
00259       
00260       if(out.get()) *out << "\nJ) Create an ifpack preconditioner precA for A ...\n";
00261 
00262       precFactory = Teuchos::rcp(new IfpackPreconditionerFactory());
00263 
00264       if(out.get()) {
00265         *out << "\nprecFactory.description() = " << precFactory->description() << std::endl;
00266         *out << "\nprecFactory.getValidParameters() =\n";
00267         precFactory->getValidParameters()->print(OSTab(out).o(),0,true,false);
00268       }
00269 
00270       RCP<Teuchos::ParameterList>
00271         ifpackPFPL = Teuchos::rcp(new Teuchos::ParameterList("IfpackPreconditionerFactory"));
00272       ifpackPFPL->set("Prec Type","ILUT");
00273       ifpackPFPL->set("Overlap",int(1));
00274       if(out.get()) {
00275         *out << "\nifpackPFPL before setting parameters =\n";
00276         ifpackPFPL->print(OSTab(out).o(),0,true);
00277       }
00278 
00279       precFactory->setParameterList(ifpackPFPL);
00280 
00281       RCP<PreconditionerBase<double> >
00282         precA = precFactory->createPrec();
00283       Thyra::initializePrec<double>(*precFactory,A,&*precA);
00284       
00285       if(out.get()) {
00286         *out << "\nifpackPFPL after setting parameters =\n";
00287         ifpackPFPL->print(OSTab(out).o(),0,true);
00288         *out << "\nprecFactory.description() = " << precFactory->description() << std::endl;
00289       }
00290       
00291       if(out.get()) *out << "\nprecA.description() = " << precA->description() << std::endl;
00292       if(out.get() && dumpAll) *out << "\ndescribe(precA) =\n" << describe(*precA,Teuchos::VERB_EXTREME);
00293       
00294       if(out.get()) *out << "\nK) Reinitialize (A,precA->getUnspecifiedPrecOp(),PRECONDITIONER_INPUT_TYPE_AS_OPERATOR) => nsA ...\n";
00295       
00296       Thyra::initializePreconditionedOp<double>(*lowsFactory,A,precA,&*nsA);
00297       
00298       if(out.get()) *out << "\nL) Testing the LinearOpWithSolveBase interface of nsA ...\n";
00299       
00300       Thyra::seed_randomize<double>(0);
00301       result = linearOpWithSolveTester.check(*nsA,out.get());
00302       if(!result) success = false;
00303 
00304       if(testTranspose && useAztecPrec) {
00305         linearOpWithSolveTester.check_adjoint_default(true);
00306         linearOpWithSolveTester.check_adjoint_residual(true);
00307       }
00308       else {
00309         linearOpWithSolveTester.check_adjoint_default(false);
00310         linearOpWithSolveTester.check_adjoint_residual(false);
00311       }
00312       
00313     }
00314     else {
00315 
00316       if(out.get()) *out << "\nSkipping testing steps J, K, and L since we are not using aztec preconditioning and therefore will not test with an ifpack preconditioner!\n";
00317 
00318     }
00319 
00320 #else // HAVE_AZTECOO_IFPACK
00321 
00322     if(out.get()) *out << "\nSkipping testing steps J, K, and L since they require ifpack and ifpack has not been enabled!\n";
00323 
00324 #endif // HAVE_AZTECOO_IFPACK
00325 
00326 
00327     if(out.get()) *out << "\nM) Scale the epetra_A object by 2.5, and then reinitialize nsA with epetra_A ...\n";
00328 
00329     Thyra::uninitializeOp<double>(*lowsFactory, nsA.ptr());
00330     // Not required but a good idea since we are changing the matrix
00331     epetra_A->Scale(2.5);
00332     initializeOp<double>(*lowsFactory, A, nsA.ptr());
00333 
00334     if(out.get()) *out << "\nN) Testing the LinearOpWithSolveBase interface of nsA ...\n";
00335     
00336     Thyra::seed_randomize<double>(0);
00337     result = linearOpWithSolveTester.check(*nsA,out.get());
00338     if(!result) success = false;
00339 
00340     if(out.get()) *out << "\nO) Create a scaled (by 2.5) copy epetra_A2 of epetra_A, and then reinitialize nsA with epetra_A2 ...\n";
00341 
00342     RCP<Epetra_CrsMatrix>
00343       epetra_A2 = Teuchos::rcp(new Epetra_CrsMatrix(*epetra_A));
00344     epetra_A2->Scale(2.5);
00345     RCP<const LinearOpBase<double> >
00346       A2 = Thyra::epetraLinearOp(epetra_A2);
00347     initializeOp<double>(*lowsFactory, A2, nsA.ptr());
00348     // Note that it was okay not to uninitialize nsA first here since A, which
00349     // was used to initialize nsA last, was not changed and therefore the
00350     // state of nsA was fine throughout
00351 
00352     if(out.get()) *out << "\nP) Testing the LinearOpWithSolveBase interface of nsA ...\n";
00353     
00354     Thyra::seed_randomize<double>(0);
00355     result = linearOpWithSolveTester.check(*nsA,out.get());
00356     if(!result) success = false;
00357 
00358     if(!useAztecPrec) {
00359 
00360       if(out.get()) *out << "\nQ) Create an implicitly scaled (by 2.5) and transposed matrix A3 = scale(2.5,transpose(A)) and initialize nsA2 ...\n";
00361     
00362       RCP<const LinearOpBase<double> >
00363         A3 = scale<double>(2.5,transpose<double>(A));
00364       RCP<LinearOpWithSolveBase<double> >
00365         nsA2 = linearOpWithSolve(*lowsFactory,A3);
00366     
00367       if(out.get()) *out << "\nR) Testing the LinearOpWithSolveBase interface of nsA2 ...\n";
00368     
00369       Thyra::seed_randomize<double>(0);
00370       result = linearOpWithSolveTester.check(*nsA2,out.get());
00371       if(!result) success = false;
00372     
00373       if(out.get()) *out << "\nS) Testing that LinearOpBase interfaces of transpose(nsA) == nsA2 ...\n";
00374     
00375       result = linearOpTester.compare(
00376         *transpose(Teuchos::rcp_implicit_cast<const LinearOpBase<double> >(nsA)),*nsA2
00377         ,out()
00378         );
00379       if(!result) success = false;
00380       
00381     }
00382     else {
00383 
00384       if(out.get()) *out << "\nSkipping testing steps Q, R, and S because we are using internal AztecOO preconditioners!\n";
00385 
00386     }
00387 
00388   if(out.get()) *out << "\nT) Running example use cases ...\n";
00389 
00390   nonExternallyPreconditionedLinearSolveUseCases(
00391     *A,*lowsFactory,testTranspose,*out
00392     );
00393 
00394   if(precFactory.get()) {
00395     externallyPreconditionedLinearSolveUseCases(
00396       *A,*lowsFactory,*precFactory,false,true,*out
00397       );
00398   }
00399 
00400 #else // SUN_CXX
00401     
00402     if(out.get()) *out << "\nTest failed since is was not even compiled since SUN_CXX was defined!\n";
00403     success = false;
00404 
00405 #endif // SUN_CXX
00406 
00407   }
00408   catch( const std::exception &excpt ) {
00409     std::cerr << "\n*** Caught standard exception : " << excpt.what() << std::endl;
00410     success = false;
00411   }
00412    
00413   return success;
00414     
00415 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines