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

Generated on Sun May 20 12:51:22 2007 for Aztecoo/Thyra Linear Solver Adapter Software by doxygen 1.3.9.1