BelosSolverFactory.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

#ifndef __Belos_SolverFactory_hpp
#define __Belos_SolverFactory_hpp

#include <BelosConfigDefs.hpp>
#include <BelosOutputManager.hpp>
#include <BelosSolverManager.hpp>

#include <BelosBlockCGSolMgr.hpp>
#include <BelosBlockGmresSolMgr.hpp>
#include <BelosGCRODRSolMgr.hpp>
#include <BelosPseudoBlockCGSolMgr.hpp>
#include <BelosPseudoBlockGmresSolMgr.hpp>
#include <BelosPseudoBlockStochasticCGSolMgr.hpp>
#include <BelosLSQRSolMgr.hpp>
#include <BelosMinresSolMgr.hpp>
#include <BelosGmresPolySolMgr.hpp>
#include <BelosPCPGSolMgr.hpp>
#include <BelosRCGSolMgr.hpp>
#include <BelosTFQMRSolMgr.hpp>
#include <BelosPseudoBlockTFQMRSolMgr.hpp>
#include <BelosFixedPointSolMgr.hpp>
#include <BelosBiCGStabSolMgr.hpp>

#include <Teuchos_Array.hpp>
#include <Teuchos_Describable.hpp>
#include <Teuchos_StandardCatchMacros.hpp>
#include <Teuchos_TypeNameTraits.hpp>

#include <algorithm>
#include <locale>
#include <map>
#include <sstream>
#include <stdexcept>
#include <vector>

namespace Belos {

namespace details {

enum EBelosSolverType {
  SOLVER_TYPE_BLOCK_GMRES,
  SOLVER_TYPE_PSEUDO_BLOCK_GMRES,
  SOLVER_TYPE_BLOCK_CG,
  SOLVER_TYPE_PSEUDO_BLOCK_CG,
  SOLVER_TYPE_GCRODR,
  SOLVER_TYPE_RCG,
  SOLVER_TYPE_MINRES,
  SOLVER_TYPE_LSQR,
  SOLVER_TYPE_STOCHASTIC_CG,
  SOLVER_TYPE_TFQMR,
  SOLVER_TYPE_PSEUDO_BLOCK_TFQMR,
  SOLVER_TYPE_GMRES_POLY,
  SOLVER_TYPE_PCPG,
  SOLVER_TYPE_FIXED_POINT,
  SOLVER_TYPE_BICGSTAB
};

} // namespace details

template<class Scalar, class MV, class OP>
class SolverFactory : public Teuchos::Describable {
public:
  typedef SolverManager<Scalar, MV, OP> solver_base_type;

  SolverFactory ();

  Teuchos::RCP<solver_base_type>
  create (const std::string& solverName,
          const Teuchos::RCP<Teuchos::ParameterList>& solverParams);

  int numSupportedSolvers () const;

  Teuchos::Array<std::string> supportedSolverNames () const;

  bool isSupported (const std::string& solverName) const;



  std::string description() const;

  void describe (Teuchos::FancyOStream& out,
                 const Teuchos::EVerbosityLevel verbLevel = Teuchos::Describable::verbLevel_default) const;

private:
  std::map<std::string, std::string> aliasToCanonicalName_;

  std::map<std::string, details::EBelosSolverType> canonicalNameToEnum_;

  void
  reviseParameterListForAlias (const std::string& aliasName,
                               Teuchos::ParameterList& solverParams);

  Teuchos::Array<std::string> canonicalSolverNames () const;

  Teuchos::Array<std::string> solverNameAliases () const;

  static void
  printStringArray (std::ostream& out,
                    const Teuchos::ArrayView<const std::string>& array)
  {
    typedef Teuchos::ArrayView<std::string>::const_iterator iter_type;

    out << "[";
    for (iter_type iter = array.begin(); iter != array.end(); ++iter) {
      out << "\"" << *iter << "\"";
      if (iter + 1 != array.end()) {
        out << ", ";
      }
    }
    out << "]";
  }
};


namespace details {

template<class SolverManagerBaseType, class SolverManagerType>
Teuchos::RCP<SolverManagerBaseType>
makeSolverManagerTmpl (const Teuchos::RCP<Teuchos::ParameterList>& params);

template<class Scalar, class MV, class OP>
Teuchos::RCP<SolverManager<Scalar, MV, OP> >
makeSolverManagerFromEnum (const EBelosSolverType solverType,
                           const Teuchos::RCP<Teuchos::ParameterList>& params)
{
  typedef SolverManager<Scalar, MV, OP> base_type;

  switch (solverType) {
  case SOLVER_TYPE_BLOCK_GMRES: {
    typedef BlockGmresSolMgr<Scalar, MV, OP> impl_type;
    return makeSolverManagerTmpl<base_type, impl_type> (params);
    break;
  }
  case SOLVER_TYPE_PSEUDO_BLOCK_GMRES: {
    typedef PseudoBlockGmresSolMgr<Scalar, MV, OP> impl_type;
    return makeSolverManagerTmpl<base_type, impl_type> (params);
    break;
  }
  case SOLVER_TYPE_BLOCK_CG: {
    typedef BlockCGSolMgr<Scalar, MV, OP> impl_type;
    return makeSolverManagerTmpl<base_type, impl_type> (params);
    break;
  }
  case SOLVER_TYPE_PSEUDO_BLOCK_CG: {
    typedef PseudoBlockCGSolMgr<Scalar, MV, OP> impl_type;
    return makeSolverManagerTmpl<base_type, impl_type> (params);
    break;
  }
  case SOLVER_TYPE_GCRODR: {
    typedef GCRODRSolMgr<Scalar, MV, OP> impl_type;
    return makeSolverManagerTmpl<base_type, impl_type> (params);
    break;
  }
  case SOLVER_TYPE_RCG: {
    typedef RCGSolMgr<Scalar, MV, OP> impl_type;
    return makeSolverManagerTmpl<base_type, impl_type> (params);
    break;
  }
  case SOLVER_TYPE_MINRES: {
    typedef MinresSolMgr<Scalar, MV, OP> impl_type;
    return makeSolverManagerTmpl<base_type, impl_type> (params);
    break;
  }
  case SOLVER_TYPE_LSQR: {
    typedef LSQRSolMgr<Scalar, MV, OP> impl_type;
    return makeSolverManagerTmpl<base_type, impl_type> (params);
    break;
  }
  case SOLVER_TYPE_STOCHASTIC_CG: {
    typedef PseudoBlockStochasticCGSolMgr<Scalar, MV, OP> impl_type;
    return makeSolverManagerTmpl<base_type, impl_type> (params);
  }
  case SOLVER_TYPE_TFQMR: {
    typedef TFQMRSolMgr<Scalar, MV, OP> impl_type;
    return makeSolverManagerTmpl<base_type, impl_type> (params);
  }
  case SOLVER_TYPE_PSEUDO_BLOCK_TFQMR: {
    typedef PseudoBlockTFQMRSolMgr<Scalar, MV, OP> impl_type;
    return makeSolverManagerTmpl<base_type, impl_type> (params);
  }
  case SOLVER_TYPE_GMRES_POLY: {
    typedef GmresPolySolMgr<Scalar, MV, OP> impl_type;
    return makeSolverManagerTmpl<base_type, impl_type> (params);
  }
  case SOLVER_TYPE_PCPG: {
    typedef PCPGSolMgr<Scalar, MV, OP> impl_type;
    return makeSolverManagerTmpl<base_type, impl_type> (params);
  }
  case SOLVER_TYPE_FIXED_POINT: {
    typedef FixedPointSolMgr<Scalar, MV, OP> impl_type;
    return makeSolverManagerTmpl<base_type, impl_type> (params);
  }
  case SOLVER_TYPE_BICGSTAB: {
    typedef BiCGStabSolMgr<Scalar, MV, OP> impl_type;
    return makeSolverManagerTmpl<base_type, impl_type> (params);
  }
  default: // Fall through; let the code below handle it.
    TEUCHOS_TEST_FOR_EXCEPTION(
      true, std::logic_error, "Belos::SolverFactory: Invalid EBelosSolverType "
      "enum value " << solverType << ".  Please report this bug to the Belos "
      "developers.");
  }

  // Compiler guard.  This may result in a warning on some compilers
  // for an unreachable statement, but it will prevent a warning on
  // other compilers for a "missing return statement at end of
  // non-void function."
  return Teuchos::null;
}

template<class SolverManagerBaseType, class SolverManagerType>
Teuchos::RCP<SolverManagerBaseType>
makeSolverManagerTmpl (const Teuchos::RCP<Teuchos::ParameterList>& params)
{
  using Teuchos::ParameterList;
  using Teuchos::parameterList;
  using Teuchos::RCP;

  RCP<SolverManagerType> solver = rcp (new SolverManagerType);

  // Some solvers may not like to get a null ParameterList.  If params
  // is null, replace it with an empty parameter list.  The solver
  // will fill in default parameters for that case.  Use the name of
  // the solver's default parameters to name the new empty list.
  RCP<ParameterList> pl;
  if (params.is_null()) {
    pl = parameterList (solver->getValidParameters ()->name ());
  } else {
    pl = params;
  }
  TEUCHOS_TEST_FOR_EXCEPTION(
    pl.is_null(), std::logic_error,
    "Belos::SolverFactory: ParameterList to pass to solver is null.  This "
    "should never happen.  Please report this bug to the Belos developers.");
  solver->setParameters (pl);
  return solver;
}

} // namespace details


template<class Scalar, class MV, class OP>
SolverFactory<Scalar, MV, OP>::SolverFactory()
{
  aliasToCanonicalName_["GMRES"] = "PSEUDOBLOCK GMRES";
  // NOTE (mfh 29 Nov 2011) Accessing the flexible capability requires
  // setting a parameter in the solver's parameter list.  This affects
  // the SolverFactory's interface, since using the "Flexible GMRES"
  // alias requires modifying the user's parameter list if necessary.
  // This is a good idea because users may not know about the
  // parameter, or may have forgotten.
  //
  // NOTE (mfh 12 Aug 2015) The keys and values need to be all uppercase.
  aliasToCanonicalName_["BLOCK GMRES"] = "BLOCK GMRES";
  aliasToCanonicalName_["FLEXIBLE GMRES"] = "BLOCK GMRES";
  aliasToCanonicalName_["CG"] = "PSEUDOBLOCK CG";
  aliasToCanonicalName_["PSEUDOBLOCKCG"] = "PSEUDOBLOCK CG";
  aliasToCanonicalName_["STOCHASTIC CG"] = "PSEUDOBLOCK STOCHASTIC CG";
  aliasToCanonicalName_["RECYCLING CG"] = "RCG";
  aliasToCanonicalName_["RECYCLING GMRES"] = "GCRODR";
  // For compatibility with Stratimikos' Belos adapter.
  aliasToCanonicalName_["PSEUDO BLOCK GMRES"] = "PSEUDOBLOCK GMRES";
  aliasToCanonicalName_["PSEUDOBLOCKGMRES"] = "PSEUDOBLOCK GMRES";
  aliasToCanonicalName_["PSEUDO BLOCK CG"] = "PSEUDOBLOCK CG";
  aliasToCanonicalName_["PSEUDOBLOCKCG"] = "PSEUDOBLOCK CG";
  aliasToCanonicalName_["TRANSPOSE-FREE QMR"] = "TFQMR";
  aliasToCanonicalName_["PSEUDO BLOCK TFQMR"] = "PSEUDOBLOCK TFQMR";
  aliasToCanonicalName_["PSEUDO BLOCK TRANSPOSE-FREE QMR"] = "PSEUDOBLOCK TFQMR";
  aliasToCanonicalName_["GMRESPOLY"] = "HYBRID BLOCK GMRES";
  aliasToCanonicalName_["SEED GMRES"] = "HYBRID BLOCK GMRES";
  aliasToCanonicalName_["CGPOLY"] = "PCPG";
  aliasToCanonicalName_["SEED CG"] = "PCPG";
  aliasToCanonicalName_["FIXED POINT"] = "FIXED POINT";
  aliasToCanonicalName_["BICGSTAB"] = "BICGSTAB";

  // Mapping from canonical solver name (a string) to its
  // corresponding enum value.  This mapping is one-to-one.
  //
  // NOTE (mfh 12 Aug 2015) The keys need to be all uppercase.
  canonicalNameToEnum_["BLOCK GMRES"] = details::SOLVER_TYPE_BLOCK_GMRES;
  canonicalNameToEnum_["PSEUDOBLOCK GMRES"] = details::SOLVER_TYPE_PSEUDO_BLOCK_GMRES;
  canonicalNameToEnum_["BLOCK CG"] = details::SOLVER_TYPE_BLOCK_CG;
  canonicalNameToEnum_["PSEUDOBLOCK CG"] = details::SOLVER_TYPE_PSEUDO_BLOCK_CG;
  canonicalNameToEnum_["PSEUDOBLOCK STOCHASTIC CG"] = details::SOLVER_TYPE_STOCHASTIC_CG;
  canonicalNameToEnum_["GCRODR"] = details::SOLVER_TYPE_GCRODR;
  canonicalNameToEnum_["RCG"] = details::SOLVER_TYPE_RCG;
  canonicalNameToEnum_["MINRES"] = details::SOLVER_TYPE_MINRES;
  canonicalNameToEnum_["LSQR"] = details::SOLVER_TYPE_LSQR;
  canonicalNameToEnum_["TFQMR"] = details::SOLVER_TYPE_TFQMR;
  canonicalNameToEnum_["PSEUDOBLOCK TFQMR"] = details::SOLVER_TYPE_PSEUDO_BLOCK_TFQMR;
  canonicalNameToEnum_["HYBRID BLOCK GMRES"] = details::SOLVER_TYPE_GMRES_POLY;
  canonicalNameToEnum_["PCPG"] = details::SOLVER_TYPE_PCPG;
  canonicalNameToEnum_["FIXED POINT"] = details::SOLVER_TYPE_FIXED_POINT;
  canonicalNameToEnum_["BICGSTAB"] = details::SOLVER_TYPE_BICGSTAB;
}


template<class Scalar, class MV, class OP>
void
SolverFactory<Scalar, MV, OP>::
reviseParameterListForAlias (const std::string& aliasName,
                             Teuchos::ParameterList& solverParams)
{
  if (aliasName == "FLEXIBLE GMRES") {
    // "Gmres" uses title case in this solver's parameter list.  For
    // our alias, we prefer the all-capitals "GMRES" that the
    // algorithm's authors (Saad and Schultz) used.
    solverParams.set ("Flexible Gmres", true);
  }
}


template<class Scalar, class MV, class OP>
Teuchos::RCP<typename SolverFactory<Scalar, MV, OP>::solver_base_type>
SolverFactory<Scalar, MV, OP>::
create (const std::string& solverName,
        const Teuchos::RCP<Teuchos::ParameterList>& solverParams)
{
  const char prefix[] = "Belos::SolverFactory: ";

  // Upper-case version of the input solver name.
  std::string solverNameUC (solverName);
  {
    typedef std::string::value_type char_t;
    typedef std::ctype<char_t> facet_type;
    const facet_type& facet = std::use_facet<facet_type> (std::locale ());

    const std::string::size_type len = solverName.size ();
    for (std::string::size_type k = 0; k < len; ++k) {
      solverNameUC[k] = facet.toupper (solverName[k]);
    }
  }

  // Check whether the given name is an alias.
  std::map<std::string, std::string>::const_iterator aliasIter =
    aliasToCanonicalName_.find (solverNameUC);
  const bool isAnAlias = (aliasIter != aliasToCanonicalName_.end());
  const std::string candidateCanonicalName =
    isAnAlias ? aliasIter->second : solverNameUC;

  // Get the canonical name.
  std::map<std::string, details::EBelosSolverType>::const_iterator canonicalIter =
    canonicalNameToEnum_.find (candidateCanonicalName);
  const bool validCanonicalName = (canonicalIter != canonicalNameToEnum_.end());

  // Check whether we found a canonical name.  If we didn't and the
  // input name is a valid alias, that's a bug.  Otherwise, the input
  // name is invalid.
  TEUCHOS_TEST_FOR_EXCEPTION
    (! validCanonicalName && isAnAlias, std::logic_error,
     prefix << "Valid alias \"" << solverName << "\" has candidate canonical "
     "name \"" << candidateCanonicalName << "\", which is not a canonical "
     "solver name.  Please report this bug to the Belos developers.");
  TEUCHOS_TEST_FOR_EXCEPTION
    (! validCanonicalName && ! isAnAlias, std::invalid_argument,
     prefix << "Invalid solver name \"" << solverName << "\".");

  // If the input list is null, we create a new list and use that.
  // This is OK because the effect of a null parameter list input is
  // to use default parameter values.  Thus, we can always replace a
  // null list with an empty list.
  Teuchos::RCP<Teuchos::ParameterList> pl =
    solverParams.is_null() ? Teuchos::parameterList() : solverParams;

  // Possibly modify the input parameter list as needed.
  if (isAnAlias) {
    reviseParameterListForAlias (solverNameUC, *pl);
  }

  return details::makeSolverManagerFromEnum<Scalar, MV, OP> (canonicalIter->second, pl);
}


template<class Scalar, class MV, class OP>
std::string
SolverFactory<Scalar, MV, OP>::description() const
{
  using Teuchos::TypeNameTraits;

  std::ostringstream out;
  out << "\"Belos::SolverFactory\": {";
  if (this->getObjectLabel () != "") {
    out << "Label: " << this->getObjectLabel () << ", ";
  }
  out << "Scalar: " << TypeNameTraits<Scalar>::name ()
      << ", MV: " << TypeNameTraits<MV>::name ()
      << ", OP: " << TypeNameTraits<OP>::name ()
      << "}";
  return out.str ();
}


template<class Scalar, class MV, class OP>
void
SolverFactory<Scalar, MV, OP>::
describe (Teuchos::FancyOStream& out,
          const Teuchos::EVerbosityLevel verbLevel) const
{
  using Teuchos::TypeNameTraits;
  using std::endl;

  const Teuchos::EVerbosityLevel vl =
    (verbLevel == Teuchos::VERB_DEFAULT) ? Teuchos::VERB_LOW : verbLevel;

  if (vl == Teuchos::VERB_NONE) {
    return;
  }

  // By convention, describe() always begins with a tab.
  Teuchos::OSTab tab0 (out);
  // The description prints in YAML format.  The class name needs to
  // be protected with quotes, so that YAML doesn't get confused
  // between the colons in the class name and the colon separating
  // (key,value) pairs.
  out << "\"Belos::SolverFactory\":" << endl;
  if (this->getObjectLabel () != "") {
    out << "Label: " << this->getObjectLabel () << endl;
  }
  {
    out << "Template parameters:" << endl;
    Teuchos::OSTab tab1 (out);
    out << "Scalar: " << TypeNameTraits<Scalar>::name () << endl
        << "MV: " << TypeNameTraits<MV>::name () << endl
        << "OP: " << TypeNameTraits<OP>::name () << endl;
  }

  // At higher verbosity levels, print out the list of supported solvers.
  if (vl > Teuchos::VERB_LOW) {
    Teuchos::OSTab tab1 (out);
    out << "Number of solvers: " << numSupportedSolvers ()
        << endl;
    out << "Canonical solver names: ";
    printStringArray (out, canonicalSolverNames ());
    out << endl;

    out << "Aliases to canonical names: ";
    printStringArray (out, solverNameAliases ());
    out << endl;
  }
}

template<class Scalar, class MV, class OP>
int
SolverFactory<Scalar, MV, OP>::numSupportedSolvers () const
{
  return static_cast<int> (canonicalNameToEnum_.size());
}

template<class Scalar, class MV, class OP>
Teuchos::Array<std::string>
SolverFactory<Scalar, MV, OP>::canonicalSolverNames () const
{
  Teuchos::Array<std::string> canonicalNames;
  typedef std::map<std::string, details::EBelosSolverType>::const_iterator iter_type;
  for (iter_type iter = canonicalNameToEnum_.begin();
       iter != canonicalNameToEnum_.end(); ++iter) {
    canonicalNames.push_back (iter->first);
  }
  return canonicalNames;
}

template<class Scalar, class MV, class OP>
Teuchos::Array<std::string>
SolverFactory<Scalar, MV, OP>::solverNameAliases () const
{
  Teuchos::Array<std::string> names;
  {
    typedef std::map<std::string, std::string>::const_iterator iter_type;
    for (iter_type iter = aliasToCanonicalName_.begin();
         iter != aliasToCanonicalName_.end(); ++iter) {
      names.push_back (iter->first);
    }
  }
  return names;
}

template<class Scalar, class MV, class OP>
Teuchos::Array<std::string>
SolverFactory<Scalar, MV, OP>::supportedSolverNames () const
{
  Teuchos::Array<std::string> names;
  {
    typedef std::map<std::string, std::string>::const_iterator iter_type;
    for (iter_type iter = aliasToCanonicalName_.begin();
         iter != aliasToCanonicalName_.end(); ++iter) {
      names.push_back (iter->first);
    }
  }
  {
    typedef std::map<std::string, details::EBelosSolverType>::const_iterator iter_type;
    for (iter_type iter = canonicalNameToEnum_.begin();
         iter != canonicalNameToEnum_.end(); ++iter) {
      names.push_back (iter->first);
    }
  }
  return names;
}

} // namespace Belos

#endif // __Belos_SolverFactory_hpp