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 );
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.get());
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 precondtioner 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); // Not required but a good idea since we are changing the matrix
00193     {
00194       Epetra_Vector diag(epetra_A->RowMap());
00195       epetra_A->ExtractDiagonalCopy(diag);
00196       diag.Scale(0.5);
00197       epetra_A->ReplaceDiagonalValues(diag);
00198     }
00199     Thyra::initializeOp<double>(*lowsFactory,A,&*nsA);
00200 
00201     // Scale the matrix back again and then reuse the preconditioner
00202     Thyra::uninitializeOp<double>(*lowsFactory,&*nsA); // Not required but a good idea since we are changing the matrix
00203     {
00204       Epetra_Vector diag(epetra_A->RowMap());
00205       epetra_A->ExtractDiagonalCopy(diag);
00206       diag.Scale(1.0/0.5);
00207       epetra_A->ReplaceDiagonalValues(diag);
00208     }
00209     initializeAndReuseOp<double>(*lowsFactory,A,&*nsA);
00210 
00211     if(out.get()) *out << "\nG) Testing the LinearOpWithSolveBase interface of nsA ...\n";
00212     
00213     Thyra::seed_randomize<double>(0);
00214     result = linearOpWithSolveTester.check(*nsA,out.get());
00215     if(!result) success = false;
00216 
00217     if(useAztecPrec) {
00218 
00219       if(out.get()) *out << "\nH) Reinitialize (A,A,PRECONDITIONER_INPUT_TYPE_AS_MATRIX) => nsA ...\n";
00220       
00221       initializeApproxPreconditionedOp<double>(*lowsFactory,A,A,&*nsA);
00222 
00223       if(out.get()) *out << "\nI) Testing the LinearOpWithSolveBase interface of nsA ...\n";
00224       
00225       Thyra::seed_randomize<double>(0);
00226       result = linearOpWithSolveTester.check(*nsA,out.get());
00227       if(!result) success = false;
00228 
00229       if(testTranspose && useAztecPrec) {
00230         linearOpWithSolveTester.check_adjoint_default(true);
00231         linearOpWithSolveTester.check_adjoint_residual(true);
00232       }
00233       else {
00234         linearOpWithSolveTester.check_adjoint_default(false);
00235         linearOpWithSolveTester.check_adjoint_residual(false);
00236       }
00237       
00238     }
00239     else {
00240 
00241       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";
00242 
00243     }
00244 
00245     
00246     RCP<PreconditionerFactoryBase<double> >
00247       precFactory;
00248 
00249 #ifdef HAVE_AZTECOO_IFPACK
00250 
00251     if(useAztecPrec) {
00252 
00253       if(testTranspose) {
00254         linearOpWithSolveTester.check_adjoint_default(true);
00255         linearOpWithSolveTester.check_adjoint_residual(true);
00256       }
00257       
00258       if(out.get()) *out << "\nJ) Create an ifpack preconditioner precA for A ...\n";
00259 
00260       precFactory = Teuchos::rcp(new IfpackPreconditionerFactory());
00261 
00262       if(out.get()) {
00263         *out << "\nprecFactory.description() = " << precFactory->description() << std::endl;
00264         *out << "\nprecFactory.getValidParameters() =\n";
00265         precFactory->getValidParameters()->print(OSTab(out).o(),0,true,false);
00266       }
00267 
00268       RCP<Teuchos::ParameterList>
00269         ifpackPFPL = Teuchos::rcp(new Teuchos::ParameterList("IfpackPreconditionerFactory"));
00270       ifpackPFPL->set("Prec Type","ILUT");
00271       ifpackPFPL->set("Overlap",int(1));
00272       if(out.get()) {
00273         *out << "\nifpackPFPL before setting parameters =\n";
00274         ifpackPFPL->print(OSTab(out).o(),0,true);
00275       }
00276 
00277       precFactory->setParameterList(ifpackPFPL);
00278 
00279       RCP<PreconditionerBase<double> >
00280         precA = precFactory->createPrec();
00281       Thyra::initializePrec<double>(*precFactory,A,&*precA);
00282       
00283       if(out.get()) {
00284         *out << "\nifpackPFPL after setting parameters =\n";
00285         ifpackPFPL->print(OSTab(out).o(),0,true);
00286         *out << "\nprecFactory.description() = " << precFactory->description() << std::endl;
00287       }
00288       
00289       if(out.get()) *out << "\nprecA.description() = " << precA->description() << std::endl;
00290       if(out.get() && dumpAll) *out << "\ndescribe(precA) =\n" << describe(*precA,Teuchos::VERB_EXTREME);
00291       
00292       if(out.get()) *out << "\nK) Reinitialize (A,precA->getUnspecifiedPrecOp(),PRECONDITIONER_INPUT_TYPE_AS_OPERATOR) => nsA ...\n";
00293       
00294       Thyra::initializePreconditionedOp<double>(*lowsFactory,A,precA,&*nsA);
00295       
00296       if(out.get()) *out << "\nL) Testing the LinearOpWithSolveBase interface of nsA ...\n";
00297       
00298       Thyra::seed_randomize<double>(0);
00299       result = linearOpWithSolveTester.check(*nsA,out.get());
00300       if(!result) success = false;
00301 
00302       if(testTranspose && useAztecPrec) {
00303         linearOpWithSolveTester.check_adjoint_default(true);
00304         linearOpWithSolveTester.check_adjoint_residual(true);
00305       }
00306       else {
00307         linearOpWithSolveTester.check_adjoint_default(false);
00308         linearOpWithSolveTester.check_adjoint_residual(false);
00309       }
00310       
00311     }
00312     else {
00313 
00314       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";
00315 
00316     }
00317 
00318 #else // HAVE_AZTECOO_IFPACK
00319 
00320     if(out.get()) *out << "\nSkipping testing steps J, K, and L since they require ifpack and ifpack has not been enabled!\n";
00321 
00322 #endif // HAVE_AZTECOO_IFPACK
00323 
00324 
00325     if(out.get()) *out << "\nM) Scale the epetra_A object by 2.5, and then reinitialize nsA with epetra_A ...\n";
00326 
00327     Thyra::uninitializeOp<double>(*lowsFactory,&*nsA); // Not required but a good idea since we are changing the matrix
00328     epetra_A->Scale(2.5);
00329     initializeOp<double>(*lowsFactory,A,&*nsA);
00330 
00331     if(out.get()) *out << "\nN) Testing the LinearOpWithSolveBase interface of nsA ...\n";
00332     
00333     Thyra::seed_randomize<double>(0);
00334     result = linearOpWithSolveTester.check(*nsA,out.get());
00335     if(!result) success = false;
00336 
00337     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";
00338 
00339     RCP<Epetra_CrsMatrix>
00340       epetra_A2 = Teuchos::rcp(new Epetra_CrsMatrix(*epetra_A));
00341     epetra_A2->Scale(2.5);
00342     RCP<const LinearOpBase<double> >
00343       A2 = Thyra::epetraLinearOp(epetra_A2);
00344     initializeOp<double>(*lowsFactory,A2,&*nsA);
00345     // Note that it was okay not to uninitialize nsA first here since A, which
00346     // was used to initialize nsA last, was not changed and therefore the
00347     // state of nsA was fine throughout
00348 
00349     if(out.get()) *out << "\nP) Testing the LinearOpWithSolveBase interface of nsA ...\n";
00350     
00351     Thyra::seed_randomize<double>(0);
00352     result = linearOpWithSolveTester.check(*nsA,out.get());
00353     if(!result) success = false;
00354 
00355     if(!useAztecPrec) {
00356 
00357       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";
00358     
00359       RCP<const LinearOpBase<double> >
00360         A3 = scale<double>(2.5,transpose<double>(A));
00361       RCP<LinearOpWithSolveBase<double> >
00362         nsA2 = linearOpWithSolve(*lowsFactory,A3);
00363     
00364       if(out.get()) *out << "\nR) Testing the LinearOpWithSolveBase interface of nsA2 ...\n";
00365     
00366       Thyra::seed_randomize<double>(0);
00367       result = linearOpWithSolveTester.check(*nsA2,out.get());
00368       if(!result) success = false;
00369     
00370       if(out.get()) *out << "\nS) Testing that LinearOpBase interfaces of transpose(nsA) == nsA2 ...\n";
00371     
00372       result = linearOpTester.compare(
00373         *transpose(Teuchos::rcp_implicit_cast<const LinearOpBase<double> >(nsA)),*nsA2
00374         ,out.get()
00375         );
00376       if(!result) success = false;
00377       
00378     }
00379     else {
00380 
00381       if(out.get()) *out << "\nSkipping testing steps Q, R, and S because we are using internal AztecOO preconditioners!\n";
00382 
00383     }
00384 
00385   if(out.get()) *out << "\nT) Running example use cases ...\n";
00386 
00387   nonExternallyPreconditionedLinearSolveUseCases(
00388     *A,*lowsFactory,testTranspose,*out
00389     );
00390 
00391   if(precFactory.get()) {
00392     externallyPreconditionedLinearSolveUseCases(
00393       *A,*lowsFactory,*precFactory,false,true,*out
00394       );
00395 
00396   }
00397 
00398 #else // SUN_CXX
00399     
00400     if(out.get()) *out << "\nTest failed since is was not even compiled since SUN_CXX was defined!\n";
00401     success = false;
00402 
00403 #endif // SUN_CXX
00404 
00405   }
00406   catch( const std::exception &excpt ) {
00407     std::cerr << "\n*** Caught standard exception : " << excpt.what() << std::endl;
00408     success = false;
00409   }
00410    
00411   return success;
00412     
00413 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends
Generated on Wed Apr 13 10:19:51 2011 for Stratimikos Package Browser (Single Doxygen Collection) by  doxygen 1.6.3