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_EpetraLinearOp.hpp"
00035 #include "Thyra_LinearOpTester.hpp"
00036 #include "Thyra_LinearOpWithSolveBase.hpp"
00037 #include "Thyra_LinearOpWithSolveTester.hpp"
00038 #include "EpetraExt_readEpetraLinearSystem.h"
00039 #include "Epetra_SerialComm.h"
00040 #include "Teuchos_ParameterList.hpp"
00041 #ifdef HAVE_AZTECOO_IFPACK
00042 #  include "Ifpack_CrsRiluk.h"
00043 #  include "Ifpack_PreconditionerFactory.hpp"
00044 #endif
00045 
00046 #endif // __sun
00047 
00048 bool Thyra::test_single_aztecoo_thyra_solver(
00049   const std::string                       matrixFile
00050   ,const bool                             testTranspose
00051   ,const int                              numRandomVectors
00052   ,const double                           maxFwdError
00053   ,const int                              maxIterations
00054   ,const double                           maxResid
00055   ,const double                           maxSolutionError
00056   ,const bool                             showAllTests
00057   ,const bool                             dumpAll
00058   ,Teuchos::ParameterList                 *fwdSolveParamList
00059   ,Teuchos::ParameterList                 *adjSolveParamList
00060   ,std::ostream                           *out
00061   )
00062 {
00063   bool result, success = true;
00064 
00065   try {
00066 
00067 #ifndef __sun
00068 
00069     if(out) {
00070       *out << "\n***"
00071            << "\n*** Testing Thyra::AztecOOLinearOpWithSolveFactory (and Thyra::AztecOOLinearOpWithSolve)"
00072            << "\n***\n"
00073            << "\nEchoing input options:"
00074            << "\n  matrixFile             = " << matrixFile
00075            << "\n  testTranspose          = " << testTranspose
00076            << "\n  numRandomVectors       = " << numRandomVectors
00077            << "\n  maxFwdError            = " << maxFwdError
00078            << "\n  maxIterations          = " << maxIterations
00079            << "\n  maxResid               = " << maxResid
00080            << "\n  showAllTests           = " << showAllTests
00081            << "\n  dumpAll                = " << dumpAll
00082            << std::endl;
00083     }
00084   
00085     const bool useAztecPrec = ( 
00086       fwdSolveParamList && fwdSolveParamList->isParameter("AZ_precond") && fwdSolveParamList->get<std::string>("AZ_precond")!="none"
00087       );
00088 
00089     if(out) {
00090       if(useAztecPrec)
00091         *out << "\nUsing aztec preconditioning so we will not test adjoint solves using internal preconditioning ...\n";
00092     }
00093     
00094     if(out) *out << "\nA) Reading in an epetra matrix A from the file \'"<<matrixFile<<"\' ...\n";
00095   
00096     const std::string indentSpacer = "  ";
00097   
00098     Epetra_SerialComm comm;
00099     Teuchos::RefCountPtr<Epetra_CrsMatrix> epetra_A;
00100     EpetraExt::readEpetraLinearSystem( matrixFile, comm, &epetra_A );
00101 
00102     Teuchos::RefCountPtr<LinearOpBase<double> > A = Teuchos::rcp(new EpetraLinearOp(epetra_A));
00103 
00104     if(out && dumpAll) *out << "\ndescribe(A) =\n" << describe(*A,Teuchos::VERB_EXTREME,indentSpacer,indentSpacer);
00105 
00106     if(out) *out << "\nB) Creating a AztecOOLinearOpWithSolveFactory object opFactory ...\n";
00107 
00108     Teuchos::RefCountPtr<const LinearOpWithSolveFactoryBase<double> >
00109       opFactory;
00110     if(1) {
00111       Teuchos::RefCountPtr<AztecOOLinearOpWithSolveFactory>
00112         aztecOpFactory = Teuchos::rcp(new AztecOOLinearOpWithSolveFactory());
00113       aztecOpFactory->fwdDefaultMaxIterations(maxIterations);
00114       aztecOpFactory->fwdDefaultTol(maxResid);
00115       aztecOpFactory->adjDefaultMaxIterations(maxIterations);
00116       aztecOpFactory->adjDefaultTol(maxResid);
00117       if(fwdSolveParamList) aztecOpFactory->setFwdAztecSolveParameters(Teuchos::rcp(fwdSolveParamList,false),true);
00118       if(adjSolveParamList) aztecOpFactory->setAdjAztecSolveParameters(Teuchos::rcp(adjSolveParamList,false),true);
00119       opFactory = aztecOpFactory;
00120     }
00121 
00122     if(out) *out << "\nC) Creating a AztecOOLinearOpWithSolve object nsA from A ...\n";
00123 
00124     Teuchos::RefCountPtr<LinearOpWithSolveBase<double> >
00125       nsA = opFactory->createOp();
00126 
00127     opFactory->initializeOp( A, &*nsA );
00128 
00129     if(out) *out << "\nD) Testing the LinearOpBase interface of nsA ...\n";
00130 
00131     LinearOpTester<double> linearOpTester;
00132     linearOpTester.check_adjoint(testTranspose);
00133     linearOpTester.num_random_vectors(numRandomVectors);
00134     linearOpTester.set_all_error_tol(maxFwdError);
00135     linearOpTester.set_all_warning_tol(1e-2*maxFwdError);
00136     linearOpTester.show_all_tests(showAllTests);
00137     linearOpTester.dump_all(dumpAll);
00138     Thyra::seed_randomize<double>(0);
00139     result = linearOpTester.check(*nsA,out,indentSpacer,indentSpacer);
00140     if(!result) success = false;
00141 
00142     if(out) *out << "\nE) Testing the LinearOpWithSolveBase interface of nsA ...\n";
00143     
00144     LinearOpWithSolveTester<double> linearOpWithSolveTester;
00145     linearOpWithSolveTester.turn_off_all_tests();
00146     linearOpWithSolveTester.check_forward_default(true);
00147     linearOpWithSolveTester.check_forward_residual(true);
00148     if(testTranspose && useAztecPrec) {
00149       linearOpWithSolveTester.check_adjoint_default(true);
00150       linearOpWithSolveTester.check_adjoint_residual(true);
00151     }
00152     else {
00153       linearOpWithSolveTester.check_adjoint_default(false);
00154       linearOpWithSolveTester.check_adjoint_residual(false);
00155     }
00156     linearOpWithSolveTester.set_all_solve_tol(maxResid);
00157     linearOpWithSolveTester.set_all_slack_error_tol(maxResid);
00158     linearOpWithSolveTester.set_all_slack_warning_tol(1e+1*maxResid);
00159     linearOpWithSolveTester.forward_default_residual_error_tol(2*maxResid);
00160     linearOpWithSolveTester.forward_default_solution_error_error_tol(maxSolutionError);
00161     linearOpWithSolveTester.adjoint_default_residual_error_tol(2*maxResid);
00162     linearOpWithSolveTester.adjoint_default_solution_error_error_tol(maxSolutionError);
00163     linearOpWithSolveTester.show_all_tests(showAllTests);
00164     linearOpWithSolveTester.dump_all(dumpAll);
00165     Thyra::seed_randomize<double>(0);
00166     result = linearOpWithSolveTester.check(*nsA,out,indentSpacer,indentSpacer);
00167     if(!result) success = false;
00168 
00169     if(out) *out << "\nF) Uninitialize nsA, create precondtioner for diagonal scaled by 0.99 and then reinitialize nsA reusing the old preconditioner ...\n";
00170 
00171     // Scale the diagonal of the matrix and then create the preconditioner for it
00172     opFactory->uninitializeOp(&*nsA); // Not required but a good idea since we are changing the matrix
00173     if(1){
00174       Epetra_Vector diag(epetra_A->RowMap());
00175       epetra_A->ExtractDiagonalCopy(diag);
00176       diag.Scale(0.5);
00177       epetra_A->ReplaceDiagonalValues(diag);
00178     }
00179     opFactory->initializeOp(A,&*nsA);
00180 
00181     // Scale the matrix back again and then reuse the preconditioner
00182     opFactory->uninitializeOp(&*nsA); // Not required but a good idea since we are changing the matrix
00183     if(1){
00184       Epetra_Vector diag(epetra_A->RowMap());
00185       epetra_A->ExtractDiagonalCopy(diag);
00186       diag.Scale(1.0/0.5);
00187       epetra_A->ReplaceDiagonalValues(diag);
00188     }
00189     opFactory->initializeAndReuseOp(A,&*nsA);
00190 
00191     if(out) *out << "\nG) Testing the LinearOpWithSolveBase interface of nsA ...\n";
00192     
00193     Thyra::seed_randomize<double>(0);
00194     result = linearOpWithSolveTester.check(*nsA,out,indentSpacer,indentSpacer);
00195     if(!result) success = false;
00196 
00197     if(useAztecPrec) {
00198 
00199       if(out) *out << "\nH) Reinitialize (A,A,PRECONDITIONER_INPUT_TYPE_AS_MATRIX) => nsA ...\n";
00200       
00201       opFactory->initializePreconditionedOp(A,A,PRECONDITIONER_INPUT_TYPE_AS_MATRIX,&*nsA);
00202 
00203       if(out) *out << "\nI) Testing the LinearOpWithSolveBase interface of nsA ...\n";
00204       
00205       Thyra::seed_randomize<double>(0);
00206       result = linearOpWithSolveTester.check(*nsA,out,indentSpacer,indentSpacer);
00207       if(!result) success = false;
00208 
00209       if(testTranspose && useAztecPrec) {
00210         linearOpWithSolveTester.check_adjoint_default(true);
00211         linearOpWithSolveTester.check_adjoint_residual(true);
00212       }
00213       else {
00214         linearOpWithSolveTester.check_adjoint_default(false);
00215         linearOpWithSolveTester.check_adjoint_residual(false);
00216       }
00217       
00218     }
00219     else {
00220 
00221       if(out) *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";
00222 
00223     }
00224 
00225 #ifdef HAVE_AZTECOO_IFPACK
00226 
00227     if(useAztecPrec) {
00228 
00229       if(testTranspose) {
00230         linearOpWithSolveTester.check_adjoint_default(true);
00231         linearOpWithSolveTester.check_adjoint_residual(true);
00232       }
00233       
00234       if(out) *out << "\nJ) Create an ifpack preconditioner precA for A ...\n";
00235       
00236       Ifpack::PreconditionerFactory ifpackPrecFactory;
00237       ifpackPrecFactory.calcCondEst(true);
00238       Teuchos::RefCountPtr<Epetra_Operator> epetra_precA;
00239       ifpackPrecFactory.setupPrec(*epetra_A,&epetra_precA,out,indentSpacer,indentSpacer);
00240       Teuchos::RefCountPtr<EpetraLinearOp>
00241         precA = Teuchos::rcp(new EpetraLinearOp(epetra_precA,NOTRANS,EPETRA_OP_APPLY_APPLY_INVERSE));
00242       
00243       if(out && dumpAll) *out << "\ndescribe(precA) =\n" << describe(*precA,Teuchos::VERB_EXTREME,indentSpacer,indentSpacer);
00244       
00245       if(out) *out << "\nK) Reinitialize (A,precA,PRECONDITIONER_INPUT_TYPE_AS_OPERATOR) => nsA ...\n";
00246       
00247       opFactory->initializePreconditionedOp(A,precA,PRECONDITIONER_INPUT_TYPE_AS_OPERATOR,&*nsA);
00248       
00249       if(out) *out << "\nL) Testing the LinearOpWithSolveBase interface of nsA ...\n";
00250       
00251       Thyra::seed_randomize<double>(0);
00252       result = linearOpWithSolveTester.check(*nsA,out,indentSpacer,indentSpacer);
00253       if(!result) success = false;
00254 
00255       if(testTranspose && useAztecPrec) {
00256         linearOpWithSolveTester.check_adjoint_default(true);
00257         linearOpWithSolveTester.check_adjoint_residual(true);
00258       }
00259       else {
00260         linearOpWithSolveTester.check_adjoint_default(false);
00261         linearOpWithSolveTester.check_adjoint_residual(false);
00262       }
00263       
00264     }
00265     else {
00266 
00267       if(out) *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";
00268 
00269     }
00270 
00271 #else // HAVE_AZTECOO_IFPACK
00272 
00273     if(out) *out << "\nSkipping testing steps J, K, and L since they require ifpack and ifpack has not been enabled!\n";
00274 
00275 #endif // HAVE_AZTECOO_IFPACK
00276 
00277 
00278     if(out) *out << "\nM) Scale the epetra_A object by 2.5, and then reinitialize nsA with epetra_A ...\n";
00279 
00280     opFactory->uninitializeOp(&*nsA); // Not required but a good idea since we are changing the matrix
00281     epetra_A->Scale(2.5);
00282     opFactory->initializeOp(A,&*nsA);
00283 
00284     if(out) *out << "\nN) Testing the LinearOpWithSolveBase interface of nsA ...\n";
00285     
00286     Thyra::seed_randomize<double>(0);
00287     result = linearOpWithSolveTester.check(*nsA,out,indentSpacer,indentSpacer);
00288     if(!result) success = false;
00289 
00290     if(out) *out << "\nO) Create a scaled (by 2.5) copy epetra_A2 of epetra_A, and then reinitialize nsA with epetra_A2 ...\n";
00291 
00292     Teuchos::RefCountPtr<Epetra_CrsMatrix>
00293       epetra_A2 = Teuchos::rcp(new Epetra_CrsMatrix(*epetra_A));
00294     epetra_A2->Scale(2.5);
00295     Teuchos::RefCountPtr<LinearOpBase<double> >
00296       A2 = Teuchos::rcp(new EpetraLinearOp(epetra_A2));
00297     opFactory->initializeOp(A2,&*nsA);
00298     // Note that it was okay not to uninitialize nsA first here since A, which
00299     // was used to initialize nsA last, was not changed and therefore the
00300     // state of nsA was fine throughout
00301 
00302     if(out) *out << "\nP) Testing the LinearOpWithSolveBase interface of nsA ...\n";
00303     
00304     Thyra::seed_randomize<double>(0);
00305     result = linearOpWithSolveTester.check(*nsA,out,indentSpacer,indentSpacer);
00306     if(!result) success = false;
00307 
00308     if(!useAztecPrec) {
00309 
00310       if(out) *out << "\nQ) Create an implicitly scaled (by 2.5) and transposed matrix A3 = scale(2.5,transpose(A)) and initialize nsA2 ...\n";
00311     
00312       Teuchos::RefCountPtr<const LinearOpBase<double> >
00313         A3 = scale(2.5,transpose(A));
00314       Teuchos::RefCountPtr<LinearOpWithSolveBase<double> >
00315         nsA2 = createAndInitializeLinearOpWithSolve(*opFactory,A3);
00316     
00317       if(out) *out << "\nR) Testing the LinearOpWithSolveBase interface of nsA2 ...\n";
00318     
00319       Thyra::seed_randomize<double>(0);
00320       result = linearOpWithSolveTester.check(*nsA2,out,indentSpacer,indentSpacer);
00321       if(!result) success = false;
00322     
00323       if(out) *out << "\nS) Testing that LinearOpBase interfaces of transpose(nsA) == nsA2 ...\n";
00324     
00325       result = linearOpTester.compare(
00326         *transpose(Teuchos::rcp_implicit_cast<const LinearOpBase<double> >(nsA)),*nsA2
00327         ,out,indentSpacer,indentSpacer
00328         );
00329       if(!result) success = false;
00330 
00331     }
00332     else {
00333 
00334       if(out) *out << "\nSkipping testing steps Q, R, and S because we are using internal AztecOO preconditioners!\n";
00335 
00336     }
00337 
00338 #else // __sun
00339     
00340     if(out) *out << "\nTest failed since is was not even compiled since __sun was defined!\n";
00341     success = false;
00342 
00343 #endif // __sun
00344 
00345   }
00346   catch( const std::exception &excpt ) {
00347     std::cerr << "*** Caught standard exception : " << excpt.what() << std::endl;
00348     success = false;
00349   }
00350    
00351   return success;
00352     
00353 }

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