test_single_stratimikos_solver.cpp

The bottom of this page shows a simple test program for how the Stratimikos::DefaultLinearSolverBuilder class is used to take a parameter list (e.g. read from a file) and use it to build Thyra::LinearOpWithSolveFactory and Thyra::PreconditionerFactoryBase objects.

The arguments that the testing program driver that calls this example are:

Teuchos::GlobalMPISession::GlobalMPISession(): started serial run

Echoing the command-line:

/home/rabartl/PROJECTS/Trilinos.base/BUILDS/CMAKE/SERIAL_DEBUG/packages/stratimikos/test/Stratimikos_test_single_stratimikos_solver_driver.exe --echo-command-line --help 

Usage: /home/rabartl/PROJECTS/Trilinos.base/BUILDS/CMAKE/SERIAL_DEBUG/packages/stratimikos/test/Stratimikos_test_single_stratimikos_solver_driver.exe [options]
  options:
  --help                         Prints this help message
  --pause-for-debugging          Pauses for user input to allow attaching a debugger
  --echo-command-line            Echo the command-line but continue as normal
  --input-file           string  Input file [Required].
                                 (default: --input-file="")
  --extra-params         string  Extra parameters overriding the parameters read in from --input-file
                                 (default: --extra-params="")

DETAILED DOCUMENTATION:

Testing program for Trilinos (and non-Trilinos) linear solvers access through Thyra.

Sample output looks like:

Teuchos::GlobalMPISession::GlobalMPISession(): started serial run

Echoing the command-line:

/home/rabartl/PROJECTS/Trilinos.base/BUILDS/CMAKE/SERIAL_DEBUG/packages/stratimikos/test/Stratimikos_test_single_stratimikos_solver_driver.exe --echo-command-line --input-file=../test/FourByFour.belos.ifpack.xml 


Reading parameters from XML file "../test/FourByFour.belos.ifpack.xml" ...

Echoing input parameters ...
 Matrix File : string = FourByFour.mtx
 Linear Solver Builder -> 
  Linear Solver Type : string = Belos
  Preconditioner Type : string = Ifpack
 LinearOpWithSolveTester -> 
  [empty list]

Validating top-level input parameters ...

Reading in an epetra matrix A from the file 'FourByFour.mtx' ...

Creating a Stratimikos::DefaultLinearSolverBuilder object ...

Valid parameters for DefaultLinearSolverBuilder ...
 Enable Delayed Solver Construction : bool = 0
 Linear Solver Type : string = Amesos
 Preconditioner Type : string = Ifpack
 Linear Solver Types -> 
  Amesos -> 
   Refactorization Policy : string = RepivotOnRefactorization
   Solver Type : string = Klu
   Throw on Preconditioner Input : bool = 1
   Amesos Settings -> 
    AddToDiag : double = 0
    AddZeroToDiag : bool = 0
    ComputeTrueResidual : bool = 0
    ComputeVectorNorms : bool = 0
    DebugLevel : int = 0
    MatrixProperty : string = general
    MaxProcs : int = -1
    NoDestroy : bool = 0
    OutputLevel : int = 1
    PrintTiming : bool = 0
    RcondThreshold : double = 1e-12
    Redistribute : bool = 0
    Refactorize : bool = 0
    Reindex : bool = 0
    ScaleMethod : int = 0
    TrustMe : bool = 0
    Lapack -> 
     Equilibrate : bool = 1
    Mumps -> 
     ColScaling : double* = 0
     Equilibrate : bool = 1
     RowScaling : double* = 0
    Pardiso -> 
     IPARM(1) : int = 0
     IPARM(10) : int = 0
     IPARM(11) : int = 0
     IPARM(18) : int = 0
     IPARM(19) : int = 0
     IPARM(2) : int = 0
     IPARM(21) : int = 0
     IPARM(3) : int = 0
     IPARM(4) : int = 0
     IPARM(8) : int = 0
     MSGLVL : int = 0
    Scalapack -> 
     2D distribution : bool = 1
     grid_nb : int = 32
    Superludist -> 
     ColPerm : string = NOT SET
     Equil : bool = 0
     Fact : string = SamePattern
     IterRefine : string = NOT SET
     PrintNonzeros : bool = 0
     ReplaceTinyPivot : bool = 1
     ReuseSymbolic : bool = 0
     RowPerm : string = NOT SET
     perm_c : int* = 0
     perm_r : int* = 0
   VerboseObject -> 
    Output File : string = none
    Verbosity Level : string = default
  AztecOO -> 
   Output Every RHS : bool = 0
   Adjoint Solve -> 
    Max Iterations : int = 400
    Tolerance : double = 1e-06
    AztecOO Settings -> 
     Aztec Preconditioner : string = ilu
     Aztec Solver : string = GMRES
     Convergence Test : string = r0
     Drop Tolerance : double = 0
     Fill Factor : double = 1
     Graph Fill : int = 0
     Ill-Conditioning Threshold : double = 1e+11
     Orthogonalization : string = Classical
     Output Frequency : int = 0
     Overlap : int = 0
     Polynomial Order : int = 3
     RCM Reordering : string = Disabled
     Size of Krylov Subspace : int = 300
     Steps : int = 3
   Forward Solve -> 
    Max Iterations : int = 400
    Tolerance : double = 1e-06
    AztecOO Settings -> 
     Aztec Preconditioner : string = ilu
     Aztec Solver : string = GMRES
     Convergence Test : string = r0
     Drop Tolerance : double = 0
     Fill Factor : double = 1
     Graph Fill : int = 0
     Ill-Conditioning Threshold : double = 1e+11
     Orthogonalization : string = Classical
     Output Frequency : int = 0
     Overlap : int = 0
     Polynomial Order : int = 3
     RCM Reordering : string = Disabled
     Size of Krylov Subspace : int = 300
     Steps : int = 3
   VerboseObject -> 
    Output File : string = none
    Verbosity Level : string = default
  Belos -> 
   Solver Type : string = Block GMRES
   Solver Types -> 
    Block CG -> 
     Adaptive Block Size : bool = 1
     Block Size : int = 1
     Convergence Tolerance : double = 1e-08
     Maximum Iterations : int = 1000
     Orthogonalization : string = DGKS
     Orthogonalization Constant : double = -1
     Output Frequency : int = -1
     Output Stream : Teuchos::RCP<std::ostream> = Teuchos::RCP<std::ostream>{ptr=0x5d8a60,node=0x5e2990,count=4}
     Show Maximum Residual Norm Only : bool = 0
     Timer Label : string = Belos
     Verbosity : int = 0
    Block GMRES -> 
     Adaptive Block Size : bool = 1
     Block Size : int = 1
     Convergence Tolerance : double = 1e-08
     Explicit Residual Scaling : string = Norm of Initial Residual
     Explicit Residual Test : bool = 0
     Flexible Gmres : bool = 0
     Implicit Residual Scaling : string = Norm of Preconditioned Initial Residual
     Maximum Iterations : int = 1000
     Maximum Restarts : int = 20
     Num Blocks : int = 300
     Orthogonalization : string = DGKS
     Orthogonalization Constant : double = -1
     Output Frequency : int = -1
     Output Stream : Teuchos::RCP<std::ostream> = Teuchos::RCP<std::ostream>{ptr=0x5d8a60,node=0x5e2730,count=4}
     Show Maximum Residual Norm Only : bool = 0
     Timer Label : string = Belos
     Verbosity : int = 0
    Pseudo Block GMRES -> 
     Adaptive Block Size : bool = 1
     Block Size : int = 1
     Convergence Tolerance : double = 1e-08
     Deflation Quorum : int = 1
     Explicit Residual Scaling : string = Norm of Initial Residual
     Implicit Residual Scaling : string = Norm of Preconditioned Initial Residual
     Maximum Iterations : int = 1000
     Maximum Restarts : int = 20
     Num Blocks : int = 300
     Orthogonalization : string = DGKS
     Orthogonalization Constant : double = -1
     Output Frequency : int = -1
     Output Stream : Teuchos::RCP<std::ostream> = Teuchos::RCP<std::ostream>{ptr=0x5d8a60,node=0x5e2860,count=4}
     Show Maximum Residual Norm Only : bool = 0
     Timer Label : string = Belos
     Verbosity : int = 0
   VerboseObject -> 
    Output File : string = none
    Verbosity Level : string = default
 Preconditioner Types -> 
  Ifpack -> 
   Overlap : int = 0
   Prec Type : string = ILU
   Ifpack Settings -> 
    amesos: solver type : string = Amesos_Klu
    fact: absolute threshold : double = 0
    fact: drop tolerance : double = 0
    fact: ict level-of-fill : double = 1
    fact: ilut level-of-fill : double = 1
    fact: level-of-fill : int = 0
    fact: relative threshold : double = 1
    fact: relax value : double = 0
    fact: sparskit: alph : double = 0
    fact: sparskit: droptol : double = 0
    fact: sparskit: lfil : int = 0
    fact: sparskit: mbloc : int = -1
    fact: sparskit: permtol : double = 0.1
    fact: sparskit: tol : double = 0
    fact: sparskit: type : string = ILUT
    partitioner: local parts : int = 1
    partitioner: overlap : int = 0
    partitioner: print level : int = 0
    partitioner: type : string = greedy
    partitioner: use symmetric graph : bool = 1
    relaxation: damping factor : double = 1
    relaxation: min diagonal value : double = 1
    relaxation: sweeps : int = 1
    relaxation: type : string = Jacobi
    relaxation: zero starting solution : bool = 1
    schwarz: combine mode : string = Zero
    schwarz: compute condest : bool = 1
    schwarz: filter singletons : bool = 0
    schwarz: reordering type : string = none
   VerboseObject -> 
    Output File : string = none
    Verbosity Level : string = default
  ML -> 
   Base Method Defaults : string = DD
   ML Settings -> 
    aggregation: damping factor : double = 1.333
    aggregation: edge prolongator drop threshold : double = 0
    aggregation: local aggregates : int = 1
    aggregation: next-level aggregates per process : int = 128
    aggregation: nodes per aggregate : int = 512
    aggregation: type : string = Uncoupled-MIS
    coarse: max size : int = 128
    coarse: pre or post : string = post
    coarse: sweeps : int = 1
    coarse: type : string = Amesos-KLU
    default values : string = maxwell
    eigen-analysis: iterations : int = 10
    eigen-analysis: type : string = cg
    increasing or decreasing : string = decreasing
    max levels : int = 10
    prec type : string = MGV
    smoother: Aztec as solver : bool = 0
    smoother: Aztec options : Teuchos::RCP<__gnu_debug_def::vector<int, std::allocator<int> > > = Teuchos::RCP<__gnu_debug_def::vector<int, std::allocator<int> > >{ptr=0x5f95e0,node=0x5fd8a0,count=2}
    smoother: Aztec params : Teuchos::RCP<__gnu_debug_def::vector<double, std::allocator<double> > > = Teuchos::RCP<__gnu_debug_def::vector<double, std::allocator<double> > >{ptr=0x5fd320,node=0x5fdfa0,count=2}
    smoother: Hiptmair efficient symmetric : bool = 1
    smoother: damping factor : double = 1
    smoother: pre or post : string = both
    smoother: sweeps : int = 1
    smoother: type : string = Hiptmair
    subsmoother: Chebyshev alpha : double = 20
    subsmoother: edge sweeps : int = 4
    subsmoother: node sweeps : int = 4
    subsmoother: type : string = Chebyshev

Creating the LinearOpWithSolveFactoryBase object lowsFactory ...

lowsFactory described as:
 Thyra::BelosLinearOpWithSolveFactory


Running example use cases for not externally preconditioned ...
 
 Running example use cases for a LinearOpWithSolveFactoryBase object ...
  
  Performing a single linear solve ...
  
  Solve status:
   solveStatus = SOLVE_STATUS_CONVERGED
   achievedTol = 1e-08
   message:
    The Belos solver of type "Belos::BlockGmresSolMgr<...,double>{Ortho Type='DGKS', Block Size=1, Num Blocks=300, Max Restarts=20}" returned a solve status of "SOLVE_STATUS_CONVERGED with total CPU time of 0.032548 sec
   extraParameters: NONE
  
  Performing a solve, changing the operator, then performing another solve ...
  
  Performing a solve, changing the operator in a very small way, then performing another solve ...
  
  Performing a solve, changing the operator in a major way, then performing another solve ...

Checking the LOWSB interface ...
 
 *** Entering LinearOpWithSolveTester<double>::check(op,...) ...
 
 describe op: Thyra::BelosLinearOpWithSolve<double>{iterativeSolver='Belos::BlockGmresSolMgr<...,double>{Ortho Type='DGKS', Block Size=1, Num Blocks=300, Max Restarts=20}',fwdOp='Thyra::EpetraLinearOp{op='Epetra_CrsMatrix',rangeDim=4,domainDim=4}',rightPrecOp='Thyra::EpetraLinearOp{op='Ifpack_AdditiveSchwarz<Ifpack_ILU>',rangeDim=4,domainDim=4}'}
 
 this->check_forward_default()==true: Checking the default forward solve ... passed!
 
 this->check_forward_residual()==true: Checking the forward solve with a tolerance on the residual ... passed!
 
 this->check_adjoint_default()==false: Skipping the check of the adjoint solve with a default tolerance!
 
 this->check_adjoint_residual()==false: Skipping the check of the adjoint solve with a tolerance on the residual!
 
 Congratulations, this LinearOpWithSolveBase object seems to check out!
 
 *** Leaving LinearOpWithSolveTester<double>::check(...)

Printing the parameter list (showing what was used) ...
 Matrix File : string = FourByFour.mtx
 Solve Adjoint : bool = 0   [default]
 Linear Solver Builder -> 
  Enable Delayed Solver Construction : bool = 0   [default]
  Linear Solver Type : string = Belos
  Preconditioner Type : string = Ifpack
  Linear Solver Types -> 
   Belos -> 
    Solver Type : string = Block GMRES   [default]
    Solver Types -> 
     Block CG -> 
      Adaptive Block Size : bool = 1   [unused]
      Block Size : int = 1   [unused]
      Convergence Tolerance : double = 1e-08   [unused]
      Maximum Iterations : int = 1000   [unused]
      Orthogonalization : string = DGKS   [unused]
      Orthogonalization Constant : double = -1   [unused]
      Output Frequency : int = -1   [unused]
      Output Stream : Teuchos::RCP<std::ostream> = Teuchos::RCP<std::ostream>{ptr=0x5d8a60,node=0x5e2990,count=6}   [unused]
      Show Maximum Residual Norm Only : bool = 0   [unused]
      Timer Label : string = Belos   [unused]
      Verbosity : int = 0   [unused]
     Block GMRES -> 
      Adaptive Block Size : bool = 1   [unused]
      Block Size : int = 1
      Convergence Tolerance : double = 1e-08
      Explicit Residual Scaling : string = Norm of Initial Residual
      Explicit Residual Test : bool = 0
      Flexible Gmres : bool = 0
      Implicit Residual Scaling : string = Norm of Preconditioned Initial Residual
      Maximum Iterations : int = 1000
      Maximum Restarts : int = 20
      Num Blocks : int = 300
      Orthogonalization : string = DGKS
      Orthogonalization Constant : double = -1
      Output Frequency : int = -1   [unused]
      Output Stream : Teuchos::RCP<std::ostream> = Teuchos::RCP<std::ostream>{ptr=0x5d8a60,node=0x5e2730,count=6}
      Show Maximum Residual Norm Only : bool = 0
      Timer Label : string = Belos
      Verbosity : int = 0
     Pseudo Block GMRES -> 
      Adaptive Block Size : bool = 1   [unused]
      Block Size : int = 1   [unused]
      Convergence Tolerance : double = 1e-08   [unused]
      Deflation Quorum : int = 1   [unused]
      Explicit Residual Scaling : string = Norm of Initial Residual   [unused]
      Implicit Residual Scaling : string = Norm of Preconditioned Initial Residual   [unused]
      Maximum Iterations : int = 1000   [unused]
      Maximum Restarts : int = 20   [unused]
      Num Blocks : int = 300   [unused]
      Orthogonalization : string = DGKS   [unused]
      Orthogonalization Constant : double = -1   [unused]
      Output Frequency : int = -1   [unused]
      Output Stream : Teuchos::RCP<std::ostream> = Teuchos::RCP<std::ostream>{ptr=0x5d8a60,node=0x5e2860,count=6}   [unused]
      Show Maximum Residual Norm Only : bool = 0   [unused]
      Timer Label : string = Belos   [unused]
      Verbosity : int = 0   [unused]
    VerboseObject -> 
     Output File : string = none
     Verbosity Level : string = default
  Preconditioner Types -> 
   Ifpack -> 
    Overlap : int = 0   [default]
    Prec Type : string = ILU   [default]
    Ifpack Settings -> 
     schwarz: combine mode : string = Zero   [default]
     schwarz: compute condest : bool = 1   [default]
     schwarz: filter singletons : bool = 0   [default]
     schwarz: reordering type : string = none   [default]
    VerboseObject -> 
     Output File : string = none   [default]
     Verbosity Level : string = default   [default]
 LinearOpWithSolveTester -> 
  All Slack Error Tol : double = 1e-05   [default]
  All Slack Warning Tol : double = 1e-06   [default]
  All Solve Tol : double = 1e-05   [default]
  Dump All : bool = 0   [default]
  Show All Tests : bool = 0   [default]

Congratulations! All of the tests checked out!

Here is the main driver itself:

// @HEADER
// ***********************************************************************
//
//         Stratimikos: Thyra-based strategies for linear solvers
//                Copyright (2006) Sandia Corporation
//
// Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
// license for use of this work by or on behalf of the U.S. Government.
//
// Redistribution and use in source and binary forms, with or without
// modification, are permitted provided that the following conditions are
// met:
//
// 1. Redistributions of source code must retain the above copyright
// notice, this list of conditions and the following disclaimer.
//
// 2. Redistributions in binary form must reproduce the above copyright
// notice, this list of conditions and the following disclaimer in the
// documentation and/or other materials provided with the distribution.
//
// 3. Neither the name of the Corporation nor the names of the
// contributors may be used to endorse or promote products derived from
// this software without specific prior written permission.
//
// THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
// EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
// PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
// PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
// NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
//
// Questions? Contact Roscoe A. Bartlett (rabartl@sandia.gov)
//
// ***********************************************************************
// @HEADER

#include "test_single_stratimikos_solver.hpp"
#include "Stratimikos_DefaultLinearSolverBuilder.hpp"
#include "Thyra_EpetraLinearOp.hpp"
#include "Thyra_LinearOpTester.hpp"
#include "Thyra_LinearOpWithSolveTester.hpp"
#include "Thyra_LinearOpWithSolveFactoryExamples.hpp"
#include "Thyra_LinearOpWithSolveFactoryHelpers.hpp"
#include "EpetraExt_readEpetraLinearSystem.h"
#include "Teuchos_ParameterList.hpp"

#ifdef HAVE_MPI
#  include "Epetra_MpiComm.h"
#else
#  include "Epetra_SerialComm.h"
#endif

bool Thyra::test_single_stratimikos_solver(
  Teuchos::ParameterList *paramList_inout
  ,const bool dumpAll
  ,Teuchos::FancyOStream *out
  )
{

  using Teuchos::rcp;
  using Teuchos::RCP;
  using Teuchos::OSTab;
  using Teuchos::ParameterList;
  using Teuchos::getParameter;
  typedef double Scalar;

  bool success = true, result = false;

  try {

    TEUCHOS_TEST_FOR_EXCEPT(!paramList_inout);

    RCP<ParameterList>
      paramList = rcp(paramList_inout,false);

    if(out) {
      *out << "\nEchoing input parameters ...\n";
      paramList->print(*out,1,true,false);
    }

    // Create list of valid parameter sublists
    Teuchos::ParameterList validParamList("test_single_stratimikos_solver");
    validParamList.set("Matrix File","fileName");
    validParamList.set("Solve Adjoint",false);
    validParamList.sublist("Linear Solver Builder").disableRecursiveValidation();
    validParamList.sublist("LinearOpWithSolveTester").disableRecursiveValidation();

    if(out) *out << "\nValidating top-level input parameters ...\n";
    paramList->validateParametersAndSetDefaults(validParamList);

    const std::string
      &matrixFile = getParameter<std::string>(*paramList,"Matrix File");
    const bool
      solveAdjoint = getParameter<bool>(*paramList,"Solve Adjoint");
    RCP<ParameterList>
      solverBuilderSL  = sublist(paramList,"Linear Solver Builder",true),
      lowsTesterSL     = sublist(paramList,"LinearOpWithSolveTester",true);

    if(out) *out << "\nReading in an epetra matrix A from the file \'"<<matrixFile<<"\' ...\n";

#ifdef HAVE_MPI
    Epetra_MpiComm comm(MPI_COMM_WORLD);
#else
    Epetra_SerialComm comm;
#endif
    RCP<Epetra_CrsMatrix> epetra_A;
    EpetraExt::readEpetraLinearSystem( matrixFile, comm, &epetra_A );

    RCP<const LinearOpBase<double> >
      A = Thyra::epetraLinearOp(epetra_A);

    if(out) *out << "\nCreating a Stratimikos::DefaultLinearSolverBuilder object ...\n";

    RCP<Thyra::LinearSolverBuilderBase<double> >
      linearSolverBuilder = rcp(new Stratimikos::DefaultLinearSolverBuilder);

    if(out) {
      *out << "\nValid parameters for DefaultLinearSolverBuilder ...\n";
      linearSolverBuilder->getValidParameters()->print(*out,1,true,false);
    }

    linearSolverBuilder->setParameterList(solverBuilderSL);

    if(out) *out << "\nCreating the LinearOpWithSolveFactoryBase object lowsFactory ...\n";
    RCP<LinearOpWithSolveFactoryBase<double> >
      lowsFactory = createLinearSolveStrategy(*linearSolverBuilder);
    if(out) *out << "\nlowsFactory described as:\n" << describe(*lowsFactory,Teuchos::VERB_MEDIUM) << std::endl;

    if(out) *out << "\nRunning example use cases for not externally preconditioned ...\n";

    TEUCHOS_ASSERT(out != NULL);
    nonExternallyPreconditionedLinearSolveUseCases(
      *A, *lowsFactory, solveAdjoint, *out
      );

    Thyra::LinearOpWithSolveTester<Scalar> linearOpWithSolveTester;
    linearOpWithSolveTester.setParameterList(lowsTesterSL);
    linearOpWithSolveTester.turn_off_all_tests();
    linearOpWithSolveTester.check_forward_default(true);
    linearOpWithSolveTester.check_forward_residual(true);
    if (solveAdjoint) {
      linearOpWithSolveTester.check_adjoint_default(true);
      linearOpWithSolveTester.check_adjoint_residual(true);
    }
    // ToDo: Use parameter lists for the above

    if(out) *out << "\nChecking the LOWSB interface ...\n";
    RCP<Thyra::LinearOpWithSolveBase<Scalar> >
      lowsA = Thyra::linearOpWithSolve<Scalar>(*lowsFactory, A);
    result = linearOpWithSolveTester.check(*lowsA, out);
    if (!result) success = false;

    if(out) {
      *out << "\nPrinting the parameter list (showing what was used) ...\n";
      paramList->print(*out,1,true,true);
    }

  }
  catch( const std::exception &excpt ) {
    std::cerr << "*** Caught standard exception : " << excpt.what() << std::endl;
    success = false;
  }

  return success;

}