Skip to content

Commit

Permalink
Stratimikos belos prec (trilinos#11222)
Browse files Browse the repository at this point in the history
* Add Belos as a preconditioner option to Stratimikos. Also clean up RCP usage.
  • Loading branch information
jennloe committed Dec 7, 2022
1 parent 27ef56d commit 3b43ef4
Show file tree
Hide file tree
Showing 15 changed files with 889 additions and 90 deletions.
Original file line number Diff line number Diff line change
Expand Up @@ -76,7 +76,6 @@

namespace Thyra {


// Constructors/initializers/accessors


Expand All @@ -93,17 +92,13 @@ bool Ifpack2PreconditionerFactory<MatrixType>::isCompatible(
const LinearOpSourceBase<scalar_type> &fwdOpSrc
) const
{
const Teuchos::RCP<const LinearOpBase<scalar_type> > fwdOp = fwdOpSrc.getOp();

typedef typename MatrixType::local_ordinal_type local_ordinal_type;
typedef typename MatrixType::global_ordinal_type global_ordinal_type;
typedef typename MatrixType::node_type node_type;

typedef Thyra::TpetraLinearOp<scalar_type, local_ordinal_type, global_ordinal_type, node_type> ThyraTpetraLinOp;
const Teuchos::RCP<const ThyraTpetraLinOp> thyraTpetraFwdOp = Teuchos::rcp_dynamic_cast<const ThyraTpetraLinOp>(fwdOp);

typedef Tpetra::Operator<scalar_type, local_ordinal_type, global_ordinal_type, node_type> TpetraLinOp;
const Teuchos::RCP<const TpetraLinOp> tpetraFwdOp = Teuchos::nonnull(thyraTpetraFwdOp) ? thyraTpetraFwdOp->getConstTpetraOperator() : Teuchos::null;
const Teuchos::RCP<const LinearOpBase<scalar_type> > fwdOp = fwdOpSrc.getOp();
using TpetraExtractHelper = TpetraOperatorVectorExtraction<scalar_type, local_ordinal_type, global_ordinal_type, node_type>;
const auto tpetraFwdOp = TpetraExtractHelper::getConstTpetraOperator(fwdOp);

const Teuchos::RCP<const MatrixType> tpetraFwdMatrix = Teuchos::rcp_dynamic_cast<const MatrixType>(tpetraFwdOp);

Expand All @@ -126,6 +121,9 @@ void Ifpack2PreconditionerFactory<MatrixType>::initializePrec(
const ESupportSolveUse /* supportSolveUse */
) const
{
using Teuchos::rcp;
using Teuchos::RCP;

// Check precondition

TEUCHOS_ASSERT(Teuchos::nonnull(fwdOpSrc));
Expand All @@ -135,7 +133,7 @@ void Ifpack2PreconditionerFactory<MatrixType>::initializePrec(
Teuchos::Time totalTimer(""), timer("");
totalTimer.start(true);

const Teuchos::RCP<Teuchos::FancyOStream> out = this->getOStream();
const RCP<Teuchos::FancyOStream> out = this->getOStream();
const Teuchos::EVerbosityLevel verbLevel = this->getVerbLevel();
Teuchos::OSTab tab(out);
if (Teuchos::nonnull(out) && Teuchos::includesVerbLevel(verbLevel, Teuchos::VERB_MEDIUM)) {
Expand All @@ -144,22 +142,19 @@ void Ifpack2PreconditionerFactory<MatrixType>::initializePrec(

// Retrieve wrapped concrete Tpetra matrix from FwdOp

const Teuchos::RCP<const LinearOpBase<scalar_type> > fwdOp = fwdOpSrc->getOp();
const RCP<const LinearOpBase<scalar_type> > fwdOp = fwdOpSrc->getOp();
TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(fwdOp));

typedef typename MatrixType::local_ordinal_type local_ordinal_type;
typedef typename MatrixType::global_ordinal_type global_ordinal_type;
typedef typename MatrixType::node_type node_type;

typedef Thyra::TpetraLinearOp<scalar_type, local_ordinal_type, global_ordinal_type, node_type> ThyraTpetraLinOp;
const Teuchos::RCP<const ThyraTpetraLinOp> thyraTpetraFwdOp = Teuchos::rcp_dynamic_cast<const ThyraTpetraLinOp>(fwdOp);
TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(thyraTpetraFwdOp));

typedef Tpetra::Operator<scalar_type, local_ordinal_type, global_ordinal_type, node_type> TpetraLinOp;
const Teuchos::RCP<const TpetraLinOp> tpetraFwdOp = thyraTpetraFwdOp->getConstTpetraOperator();
using TpetraExtractHelper = TpetraOperatorVectorExtraction<scalar_type, local_ordinal_type, global_ordinal_type, node_type>;
const auto tpetraFwdOp = TpetraExtractHelper::getConstTpetraOperator(fwdOp);
TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(tpetraFwdOp));

const Teuchos::RCP<const MatrixType> tpetraFwdMatrix = Teuchos::rcp_dynamic_cast<const MatrixType>(tpetraFwdOp);
const RCP<const MatrixType> tpetraFwdMatrix = Teuchos::rcp_dynamic_cast<const MatrixType>(tpetraFwdOp);
TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(tpetraFwdMatrix));

// Retrieve concrete preconditioner object
Expand All @@ -168,59 +163,57 @@ void Ifpack2PreconditionerFactory<MatrixType>::initializePrec(
Teuchos::ptr(dynamic_cast<DefaultPreconditioner<scalar_type> *>(prec));
TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(defaultPrec));

// Process parameter list
// This check needed to address Issue #535.
// Corresponding fix exists in stratimikos/adapters/belos/tpetra/Thyra_BelosTpetraPreconditionerFactory_def.hpp
RCP<Teuchos::ParameterList> innerParamList;
if (paramList_.is_null ()) {
innerParamList = rcp(new Teuchos::ParameterList(*getValidParameters()));
}
else {
innerParamList = paramList_;
}

bool useHalfPrecision = false;
if (paramList_->isParameter("half precision"))
useHalfPrecision = paramList_->get<bool>("half precision");
if (innerParamList->isParameter("half precision"))
useHalfPrecision = Teuchos::getParameter<bool>(*innerParamList, "half precision");

Teuchos::RCP<const Teuchos::ParameterList> constParamList = paramList_;
if (constParamList.is_null ()) {
constParamList = getValidParameters ();
}
const std::string preconditionerType = Teuchos::getParameter<std::string>(*constParamList, "Prec Type");
const Teuchos::RCP<const Teuchos::ParameterList> packageParamList = Teuchos::sublist(constParamList, "Ifpack2 Settings");
const std::string preconditionerType = Teuchos::getParameter<std::string>(*innerParamList, "Prec Type");
const RCP<Teuchos::ParameterList> packageParamList = Teuchos::rcpFromRef(innerParamList->sublist("Ifpack2 Settings"));

// precTypeUpper is the upper-case version of preconditionerType.
std::string precTypeUpper (preconditionerType);
if (precTypeUpper.size () > 0) {
for (size_t k = 0; k < precTypeUpper.size (); ++k) {
precTypeUpper[k] = ::toupper(precTypeUpper[k]);
}
}
std::transform(precTypeUpper.begin(), precTypeUpper.end(),precTypeUpper.begin(), ::toupper);

// mfh 09 Nov 2013: If the Ifpack2 list doesn't already have the
// "schwarz: overlap level" parameter, then override it with the
// value of "Overlap". This avoids use of the newly deprecated
// three-argument version of Ifpack2::Factory::create() that takes
// the overlap as an integer.
if (constParamList->isType<int> ("Overlap") && ! packageParamList.is_null () && ! packageParamList->isType<int> ("schwarz: overlap level") &&
if (innerParamList->isType<int> ("Overlap") && ! packageParamList.is_null () && ! packageParamList->isType<int> ("schwarz: overlap level") &&
precTypeUpper == "SCHWARZ") {
const int overlap = constParamList->get<int> ("Overlap");
Teuchos::RCP<Teuchos::ParameterList> nonconstPackageParamList =
Teuchos::sublist (paramList_, "Ifpack2 Settings");
nonconstPackageParamList->set ("schwarz: overlap level", overlap);
const int overlap = innerParamList->get<int> ("Overlap");
packageParamList->set ("schwarz: overlap level", overlap);
}

// Create the initial preconditioner

if (Teuchos::nonnull(out) && Teuchos::includesVerbLevel(verbLevel, Teuchos::VERB_LOW)) {
if (Teuchos::nonnull(out) && Teuchos::includesVerbLevel(verbLevel, Teuchos::VERB_MEDIUM)) {
*out << "\nCreating a new Ifpack2::Preconditioner object...\n";
}
timer.start(true);

Teuchos::RCP<LinearOpBase<scalar_type> > thyraPrecOp;
RCP<LinearOpBase<scalar_type> > thyraPrecOp;

typedef Ifpack2::Preconditioner<scalar_type, local_ordinal_type, global_ordinal_type, node_type> Ifpack2Prec;
Teuchos::RCP<Ifpack2Prec> concretePrecOp;
RCP<Ifpack2Prec> concretePrecOp;

#ifdef THYRA_IFPACK2_ENABLE_HALF_PRECISION
// CAG: There is nothing special about the combination double-float,
// except that I feel somewhat confident that Trilinos builds
// with both scalar types.
typedef typename Teuchos::ScalarTraits<scalar_type>::halfPrecision half_scalar_type;
typedef Ifpack2::Preconditioner<half_scalar_type, local_ordinal_type, global_ordinal_type, node_type> HalfIfpack2Prec;
Teuchos::RCP<HalfIfpack2Prec> concretePrecOpHalf;
RCP<HalfIfpack2Prec> concretePrecOpHalf;
#endif

if (useHalfPrecision) {
Expand All @@ -236,7 +229,7 @@ void Ifpack2PreconditionerFactory<MatrixType>::initializePrec(
concretePrecOpHalf =
Ifpack2::Factory::create<row_matrix_type> (preconditionerType, tpetraFwdMatrixHalf);
#else
TEUCHOS_TEST_FOR_EXCEPT(true);
TEUCHOS_TEST_FOR_EXCEPTION(true, std::logic_error, "Ifpack2 does not have correct precisions enabled to use half precision.")
#endif
} else {
typedef Tpetra::RowMatrix<scalar_type, local_ordinal_type,
Expand Down Expand Up @@ -268,7 +261,7 @@ void Ifpack2PreconditionerFactory<MatrixType>::initializePrec(
concretePrecOp->compute();

// Wrap concrete preconditioner
thyraPrecOp = Thyra::createLinearOp(Teuchos::RCP<TpetraLinOp>(concretePrecOp));
thyraPrecOp = Thyra::createLinearOp(RCP<TpetraLinOp>(concretePrecOp));
}

defaultPrec->initializeUnspecified(thyraPrecOp);
Expand Down Expand Up @@ -325,7 +318,7 @@ void Ifpack2PreconditionerFactory<MatrixType>::setParameterList(
{
TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(paramList));

const Teuchos::RCP<const Teuchos::ParameterList> validParamList = this->getValidParameters();
const auto validParamList = this->getValidParameters();
paramList->validateParametersAndSetDefaults(*validParamList, 0);
paramList->validateParameters(*validParamList, 1);

Expand Down
16 changes: 12 additions & 4 deletions packages/stratimikos/adapters/belos/src/CMakeLists.txt
Original file line number Diff line number Diff line change
Expand Up @@ -4,10 +4,18 @@ INCLUDE(TrilinosCreateClientTemplateHeaders)
SET(HEADERS "")
SET(SOURCES "")

SET_AND_INC_DIRS(DIR ${CMAKE_CURRENT_SOURCE_DIR})
APPEND_GLOB(HEADERS ${DIR}/*.hpp)
APPEND_GLOB(SOURCES ${DIR}/*.cpp)
TRILINOS_CREATE_CLIENT_TEMPLATE_HEADERS(${DIR})
SET_AND_INC_DIRS(srcDir ${CMAKE_CURRENT_SOURCE_DIR})
APPEND_GLOB(HEADERS ${srcDir}/*.hpp)
APPEND_GLOB(SOURCES ${srcDir}/*.cpp)
TRILINOS_CREATE_CLIENT_TEMPLATE_HEADERS(${srcDir})

IF(${PACKAGE_NAME}_ENABLE_ThyraTpetraAdapters)
SET_AND_INC_DIRS(tpetraSrcDir ${CMAKE_CURRENT_SOURCE_DIR}/tpetra)
ENDIF()

APPEND_GLOB(HEADERS ${tpetraSrcDir}/*.hpp)
APPEND_GLOB(SOURCES ${tpetraSrcDir}/*.cpp)

APPEND_GLOB(HEADERS ${CMAKE_CURRENT_BINARY_DIR}/*.hpp)
INCLUDE_DIRECTORIES(${CMAKE_CURRENT_BINARY_DIR} ${CMAKE_CURRENT_BINARY_DIR}/../../../src)

Expand Down
36 changes: 36 additions & 0 deletions packages/stratimikos/adapters/belos/src/tpetra/README.md
Original file line number Diff line number Diff line change
@@ -0,0 +1,36 @@
# Belos Tpetra Preconditioning in Stratimikos

The adapters in this directory allow a (unpreconditioned) Belos solver using Tpetra linear algebra to be called from a Stratimikos input deck as a preconditioner to another Belos solver.

This can be used, for example, to apply Belos' polynomial preconditioner using Stratimikos.
It can also be used with the fixed-point solver as the outer solver to perform iterative refinement around an inner solver.
(See example in...)



## Usage

To enable the Belos-as-preconditioner functionality, simply include the following in your Stratimikos code:

```
#include <Stratimikos_BelosTpetraPrecHelpers.hpp>
```

Then, after creating your Stratimikos linearSolverBuilder, enable Belos-as-preconditioner as follows:

```
// Create Stratimikos::LinearSolverBuilder
Stratimikos::LinearSolverBuilder<Scalar> linearSolverBuilder;
// Register Belos+Tpetra as preconditioner:
// (Note: One could use Tpetra::Operator instead of Tpetra::CrsMatrix.)
Stratimikos::enableBelosPrecTpetra<Tpetra::CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node>>(linearSolverBuilder);
```

## Examples

One can find the following example parameter list inputs in `trilinos/packages/stratimikos/test/`

+ list1.xml
+ list2.xml
+ list3.xml
Original file line number Diff line number Diff line change
@@ -0,0 +1,74 @@
// @HEADER
// ***********************************************************************
//
// Stratimikos: Thyra-based strategies for linear solvers
// Copyright 2012 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 Jennifer A. Loe (jloe@sandia.gov)
//
// ***********************************************************************
// @HEADER
#ifndef STRATIMIKOS_BELOS_PREC_TPETRA_HELPERS_HPP
#define STRATIMIKOS_BELOS_PREC_TPETRA_HELPERS_HPP

#include "Stratimikos_LinearSolverBuilder.hpp"
#include "Thyra_BelosTpetraPreconditionerFactory_decl.hpp"
#include "Thyra_BelosTpetraPreconditionerFactory_def.hpp"

#include "Teuchos_RCP.hpp"
#include "Teuchos_ParameterList.hpp"
#include "Teuchos_TestForException.hpp"
#include "Teuchos_AbstractFactoryStd.hpp"

#include <string>

namespace Stratimikos {

template <typename MatrixType>
void enableBelosPrecTpetra(LinearSolverBuilder<typename MatrixType::scalar_type>& builder, const std::string& stratName = "BelosPrecTpetra")
{
const Teuchos::RCP<const Teuchos::ParameterList> precValidParams = Teuchos::sublist(builder.getValidParameters(), "Preconditioner Types");

TEUCHOS_TEST_FOR_EXCEPTION(precValidParams->isParameter(stratName), std::logic_error,
"Stratimikos::enableBelosPrecTpetra cannot add \"" + stratName +"\" because it is already included in builder!");

typedef typename MatrixType::scalar_type scalar_type;
typedef Thyra::PreconditionerFactoryBase<scalar_type> Base;
typedef Thyra::BelosTpetraPreconditionerFactory<MatrixType> Impl;

builder.setPreconditioningStrategyFactory(Teuchos::abstractFactoryStd<Base, Impl>(), stratName);
}

} // namespace Stratimikos

#endif
Loading

0 comments on commit 3b43ef4

Please sign in to comment.