BelosStatusTestFactory.hpp

Go to the documentation of this file.

//@HEADER
// ************************************************************************
//
//                 Belos: Block Linear Solvers Package
//                  Copyright 2004 Sandia Corporation
//
// Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
// the U.S. Government retains certain rights in this software.
//
// 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 Michael A. Heroux (maherou@sandia.gov)
//
// ************************************************************************
//@HEADER

#include <BelosStatusTestCombo.hpp>
#include <BelosStatusTestGenResNorm.hpp>
#include <BelosStatusTestImpResNorm.hpp>
#include <BelosStatusTestMaxIters.hpp>
#include <Teuchos_ParameterList.hpp>

namespace Belos {

  template<class Scalar, class MV, class OP>
  class StatusTestFactory {
  public:
    typedef typename Teuchos::ScalarTraits<Scalar>::magnitudeType magnitude_type;
    typedef StatusTest<Scalar,MV,OP> base_test;
    typedef StatusTestGenResNorm<Scalar,MV,OP> res_norm_test;
    typedef StatusTestMaxIters<Scalar,MV,OP> max_iter_test;
    typedef StatusTestCombo<Scalar,MV,OP> combo_test;

    static Teuchos::RCP<base_test>
    gmresTest (const magnitude_type convTol,
         const int maxIterCount,
         const bool haveLeftPreconditioner,
         const ScaleType implicitScaleType = Belos::NormOfPrecInitRes,
         const ScaleType explicitScaleType = Belos::NormOfInitRes,
         const int blockSize = 1,
         const int defQuorum = -1,
         const bool showMaxResNormOnly = false);

    static Teuchos::RCP<base_test>
    gmresTest (const bool haveLeftPreconditioner, 
         const Teuchos::RCP<Teuchos::ParameterList>& params);

    static std::pair<bool, bool>
    changeConvTolAndMaxIters (const Teuchos::RCP<base_test>& test, 
            const magnitude_type convTol,
            const int maxIterCount);

    static bool
    changeMaxNumIters (const Teuchos::RCP<base_test>& test, 
           const int maxIterCount);

    static bool
    changeConvTol (const Teuchos::RCP<base_test>& test, 
       const magnitude_type convTol);

    static ScaleType
    stringToScaleType (const std::string& scaleType);

    Teuchos::RCP<const Teuchos::ParameterList> getValidGmresParameters ();

  private:
    Teuchos::RCP<const Teuchos::ParameterList> defaultGmresParams_;
  };


  template<class Scalar, class MV, class OP>
  Teuchos::RCP<typename StatusTestFactory<Scalar, MV, OP>::base_test>
  StatusTestFactory<Scalar, MV, OP>::
  gmresTest (const typename StatusTestFactory<Scalar, MV, OP>::magnitude_type convTol,
       const int maxIterCount,
       const bool haveLeftPreconditioner,
       const ScaleType implicitScaleType = NormOfPrecInitRes,
       const ScaleType explicitScaleType = NormOfInitRes,
       const int blockSize = 1,
       const int defQuorum = -1,
       const bool showMaxResNormOnly = false)
  {
    using Teuchos::null;
    using Teuchos::ParameterList;
    using Teuchos::RCP;
    using Teuchos::rcp;

    TEUCHOS_TEST_FOR_EXCEPTION(blockSize < 1, std::invalid_argument,
           "blockSize (= " << blockSize << ") must be >= 1.");
    TEUCHOS_TEST_FOR_EXCEPTION(defQuorum > blockSize, std::invalid_argument,
           "defQuorum (= " << defQuorum << ") may be no larger "
           "than blockSize (= " << blockSize << ").");

    // The "implicit" residual test checks the "native" residual
    // norm to determine if convergence was achieved.  It is less
    // expensive than the "explicit" residual test.
    RCP<res_norm_test> implicitTest 
      = rcp (new res_norm_test (convTol, defQuorum));
    implicitTest->defineScaleForm (stringToScaleType (implicitScaleType), 
           Belos::TwoNorm);
    implicitTest->setShowMaxResNormOnly (showMaxResNormOnly);

    // If there's a left preconditioner, create a combined status
    // test that check first the "explicit," then the "implicit"
    // residual norm, requiring that both have converged to within
    // the specified tolerance.  Otherwise, we only perform the
    // "implicit" test.
    RCP<res_norm_test> explicitTest;
    if (haveLeftPreconditioner) { // ! problem_->getLeftPrec().is_null()
      explicitTest = rcp (new res_norm_test (convTol, defQuorum));
      explicitTest->defineResForm (res_norm_test::Explicit, Belos::TwoNorm);
      explicitTest->defineScaleForm (stringToScaleType (explicitScaleType),
             Belos::TwoNorm);
      explicitTest->setShowMaxResNormOnly (showMaxResNormOnly);
    }
    else {
      explicitTest = null;
    }

    // Full convergence test:
    //
    // First, the implicit residual norm test,
    // Followed by the explicit residual norm test if applicable,
    // Followed by the user-defined convergence test if supplied.
    RCP<base_test> convTest;
    if (explicitTest.is_null()) {
      convTest = implicitTest;
    } 
    else {
      // The "explicit" residual test is only performed once the
      // native ("implicit") residual is below the convergence
      // tolerance.
      convTest = rcp (new combo_test (combo_test::SEQ, 
              implicitTest, 
              explicitTest));
    }

    // Stopping criterion for maximum number of iterations.
    RCP<max_iter_test> maxIterTest = rcp (new max_iter_test (maxIterCount));

    // The "final" stopping criterion:
    //
    // Either we've run out of iterations, OR we've converged.
    return rcp (new combo_test (combo_test::OR, maxIterTest, convTest));
  }


  template<class Scalar, class MV, class OP>
  Teuchos::RCP<typename StatusTestFactory<Scalar, MV, OP>::base_test>
  StatusTestFactory<Scalar, MV, OP>::
  gmresTest (const bool haveLeftPreconditioner, 
       const Teuchos::RCP<Teuchos::ParameterList>& params)
  {
    using Teuchos::Exceptions::InvalidParameter;
    using std::string;
    typedef Teuchos::ScalarTraits<magnitude_type> STM;

    // "Convergence Tolerance" is a required parameter and must be
    // nonnegative.
    const magnitude_type convTol = 
      params->get<magnitude_type> ("Convergence Tolerance");
    TEUCHOS_TEST_FOR_EXCEPTION(convTol < STM::zero(), std::invalid_argument,
           "Convergence tolerance " << convTol 
           << " is negative.");
    // "Maximum Iterations" is a required parameter and must be nonnegative.
    const int maxIterCount = params->get<int> ("Maximum Iterations");
    TEUCHOS_TEST_FOR_EXCEPTION(maxIterCount < 0, std::invalid_argument,
           "Maximum number of iterations " << maxIterCount
           << " is negative.");

    // PseudoBlockGmresSolMgr uses as defaults the preconditioned
    // initial residual for the implicit test, and the
    // unpreconditioned initial residual for the explicit test.
    ScaleType implicitScaleType = NormOfPrecInitRes;
    {
      const std::string defaultImplicitScaleType ("Norm of Preconditioned Initial Residual");
      implicitScaleType = stringToScaleType (params->get ("Implicit Residual Scaling", defaultImplicitScaleType));
    }
    ScaleType explicitScaleType = NormOfInitRes;
    {
      const std::string defaultExplicitScaleType ("Norm of Initial Residual");
      explicitScaleType = stringToScaleType (params->get ("Explicit Residual Scaling", defaultExplicitScaleType));
    }
    const int defaultBlockSize = 1;
    int blockSize = params->get ("Block Size", defaultBlockSize);

    const int defaultDefQuorum = -1;
    int defQuorum = params->get ("Deflation Quorum", defaultDefQuorum);

    const bool defaultShowMaxResNormOnly = false;
    bool showMaxResNormOnly = params->get ("Show Maximum Residual Norm Only", defaultShowMaxResNormOnly);

    return gmresTest (convTol, maxIterCount, haveLeftPreconditioner,
          implicitScaleType, explicitScaleType, blockSize,
          defQuorum, showMaxResNormOnly);
  }

  template<class Scalar, class MV, class OP>
  Teuchos::RCP<const Teuchos::ParameterList>
  getValidGmresParameters ()
  {
    using Teuchos::ParameterList;
    using Teuchos::parameterList;
    using Teuchos::RCP;
    typedef Teuchos::ScalarTraits<Scalar> STS;
    typedef Teuchos::ScalarTraits<magnitude_type> STM;

    if (defaultGmresParams_.is_null()) {
      // These two parameters are required.  We supply defaults only
      // as examples.
      const magnitude_type convTol = STS::squareroot (STS::eps());
      const int maxIterCount = 200;

      // PseudoBlockGmresSolMgr uses as defaults the preconditioned
      // initial residual for the implicit test, and the
      // unpreconditioned initial residual for the explicit test.  We
      // take these as the standard for all GMRES variants.
      const std::string defaultImplicitScaleType ("Norm of Preconditioned Initial Residual");
      const std::string defaultExplicitScaleType ("Norm of Initial Residual");
      const int defaultBlockSize = 1;
      const int defaultDefQuorum = -1;
      const bool defaultShowMaxResNormOnly = false;

      RCP<ParameterList> params = parameterList ();
      params->set ("Convergence Tolerance", convTol);
      params->set ("Maximum Iterations", maxIterCount);
      params->set ("Implicit Residual Scaling", defaultImplicitScaleType);
      params->set ("Explicit Residual Scaling", defaultExplicitScaleType);
      params->set ("Block Size", defaultBlockSize);
      params->set ("Deflation Quorum", defaultDefQuorum);
      params->set ("Show Maximum Residual Norm Only", defaultShowMaxResNormOnly);

      defaultGmresParams_ = params;
    }
    return defaultGmresParams_;
  }


  template<class Scalar, class MV, class OP>
  bool
  StatusTestFactory<Scalar, MV, OP>::
  changeMaxNumIters (const Teuchos::RCP<typename StatusTestFactory<Scalar, MV, OP>::base_test>& test,
         const int maxIterCount)
  {
    using Teuchos::rcp_dynamic_cast;
    using Teuchos::nonnull;
    using Teuchos::RCP;
    using Teuchos::rcp;

    // We declare "success" if the test or at least one of its
    // children (if it's a combo_test) is a max_iter_test.
    bool success = false;
    RCP<max_iter_test> maxIterTest = rcp_dynamic_cast<max_iter_test> (test);
    if (nonnull (maxIterTest))
      {
  test->setMaxIters (maxIterCount);
  success = true;
      }
    else
      {
  RCP<combo_test> comboTest = rcp_dynamic_cast<combo_test> (test);
  if (nonnull (comboTest))
    {
      typedef typename combo_test::st_vector st_vector;
      typedef typename st_vector::size_type size_type;
      st_vector tests = test->getStatusTests ();
      // We could use boost lambda to remove this for loop.
      // Standard STL doesn't have a version of mem_fun for
      // member functions that take > 1 argument.
      for (size_type k = 0; result || k < tests.end(); ++k)
        { // Recurse on all children, since it's possible for
    // more than one child to be a max_iter_test.
    const bool result =
      changeMaxNumIters (tests[k], maxIterCount);
    success = result || success;
        }
    }
      }
    return success;
  }

  template<class Scalar, class MV, class OP>
  std::pair<bool, bool>
  StatusTestFactory<Scalar, MV, OP>::
  changeConvTolAndMaxNumIters (const Teuchos::RCP<typename StatusTestFactory<Scalar, MV, OP>::base_test>& test,
             const typename StatusTestFactory<Scalar, MV, OP>::magnitude_type convTol,
             const int maxIterCount)
  {
    using Teuchos::rcp_dynamic_cast;
    using Teuchos::nonnull;
    using Teuchos::RCP;
    using Teuchos::rcp;
    typedef StatusTestResNorm<Scalar,MV,OP> res_norm_base_test;
    RCP<max_iter_test> maxIterTest = rcp_dynamic_cast<max_iter_test> (test);

    // We declare "success" if we encountered a res_norm_base_test
    // _and_ a max_iter_test somewhere along the recursion path.
    bool foundResNormTest = false;
    bool foundMaxIterTest = false;

    RCP<res_norm_base_test> normTest = 
      rcp_dynamic_cast<res_norm_base_test> (test);
    if (nonnull (normTest))
      { 
  // NOTE (mfh 03 Mar 2011) setTolerance() returns an int
  // result.  However, all subclasses' implementations return 0
  // here, and all of them always (re)set the tolerance, so I
  // think it's OK to ignore the result.
  (void) test->setTolerance (convTol);
  foundResNormTest = true;
      }
    else 
      {
  RCP<max_iter_test> maxIterTest = 
    rcp_dynamic_cast<max_iter_test> (test);
  if (nonnull (maxIterTest))
    {
      test->setMaxIters (maxIterCount);
      foundMaxIterTest = true;
    }
      }
    if (! foundResNormTest && ! foundMaxIterTest)
      {
  RCP<combo_test> comboTest = rcp_dynamic_cast<combo_test> (test);
  if (nonnull (comboTest))
    {
      typedef typename combo_test::st_vector st_vector;
      typedef typename st_vector::size_type size_type;
      st_vector tests = test->getStatusTests ();
      // We could use boost lambda to remove this for loop.
      // Standard STL doesn't have a version of mem_fun for
      // member functions that take > 1 argument.
      for (size_type k = 0; result || k < tests.end(); ++k)
        { // Recurse on all children.
    const std::pair<bool, bool> result = 
      changeConvTolAndMaxIters (tests[k], convTol, maxIterCount);
    foundResNormTest = result.first || foundResNormTest;
    foundMaxIterTest = result.second || foundMaxIterTest;     
        }
    }
      }
    return std::make_pair (foundResNormTest, foundMaxIterTest);
  }

  template<class Scalar, class MV, class OP>
  bool
  StatusTestFactory<Scalar, MV, OP>::
  changeConvTol (const Teuchos::RCP<typename StatusTestFactory<Scalar, MV, OP>::base_test>& test,
     const typename StatusTestFactory<Scalar, MV, OP>::magnitude_type convTol)
  {
    using Teuchos::rcp_dynamic_cast;
    using Teuchos::nonnull;
    using Teuchos::RCP;
    using Teuchos::rcp;
    typedef StatusTestResNorm<Scalar,MV,OP> res_norm_base_test;

    // We declare "success" if the test or at least one of its
    // children (if it's a combo_test) is a res_norm_base_test.
    bool success = false;
    RCP<res_norm_base_test> normTest = 
      rcp_dynamic_cast<res_norm_base_test> (test);
    if (nonnull (normTest))
      { 
  // NOTE (mfh 03 Mar 2011) setTolerance() returns an int
  // result.  However, all subclasses' implementations return 0
  // here, and all of them always (re)set the tolerance, so I
  // think it's OK to ignore the result.
  (void) test->setTolerance (convTol);
  success = true;
      }
    else
      {
  RCP<combo_test> comboTest = rcp_dynamic_cast<combo_test> (test);
  if (nonnull (comboTest))
    {
      typedef typename combo_test::st_vector st_vector;
      typedef typename st_vector::size_type size_type;
      st_vector tests = test->getStatusTests ();
      // We could use boost lambda to remove this for loop.
      // Standard STL doesn't have a version of mem_fun for
      // member functions that take > 1 argument.
      for (size_type k = 0; result || k < tests.end(); ++k)
        { // Recurse on all children, since it's possible for
    // more than one child to be a res_norm_base_test.
    const bool result = changeConvTol (tests[k], convTol);
    success = result || success;
        }
    }
      }
    return success;
  }


  template<class Scalar, class MV, class OP>
  static ScaleType
  StatusTestFactory<Scalar, MV, OP>::
  stringToScaleType (const std::string& scaleType) 
  {
    const char* validNames[] = {
      "Norm of Initial Residual", 
      "Norm of Preconditioned Initial Residual",
      "Norm of RHS",
      "Norm of Right-Hand Side",
      "None"
    };
    const int numValidNames = 5;
    const ScaleType correspondingOutputs[] = {
      Belos::NormOfInitRes,
      Belos::NormOfPrecInitRes,
      Belos::NormOfRHS,
      Belos::NormOfRHS,
      Belos::None
    };
    for (int k = 0; k < numValidNames; ++k)
      {
  if (scaleType == validNames[k])
    return correspondingOutputs[k];
      }
    TEUCHOS_TEST_FOR_EXCEPTION (true, std::logic_error,
      "Invalid residual scaling type \"" << scaleType 
      << "\".");
  }

} // namespace Belos