Testing Program for the AztecOO/Thyra Linear Solver Adapters
[Assorted AztecOO to Thyra Linear Solver Adapter Example Code]

Here we show a testing program that takes the AztecOO/%Thyra adatpers through the paces. More...
Below is a complete listing of the testing function with line numbers:

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

Generated on Thu Sep 18 12:29:57 2008 for Aztecoo/Thyra Linear Solver Adapter Software by doxygen 1.3.9.1