Skip to content

Commit

Permalink
MueLu: Add matrix-free P and R factories and associated operators
Browse files Browse the repository at this point in the history
My magnum opus... for now :)
  • Loading branch information
GrahamBenHarper committed Aug 24, 2022
1 parent ba4d352 commit a2ed36d
Show file tree
Hide file tree
Showing 15 changed files with 1,254 additions and 21 deletions.
2 changes: 2 additions & 0 deletions packages/muelu/src/CMakeLists.txt
Original file line number Diff line number Diff line change
Expand Up @@ -50,6 +50,7 @@ INCLUDE_DIRECTORIES(${DIR}/Transfers/Energy-Minimization/Solvers)
INCLUDE_DIRECTORIES(${DIR}/Transfers/GeneralGeometric)
INCLUDE_DIRECTORIES(${DIR}/Transfers/Generic)
INCLUDE_DIRECTORIES(${DIR}/Transfers/Geometric-Interpolation)
INCLUDE_DIRECTORIES(${DIR}/Transfers/Matrix-Free)
INCLUDE_DIRECTORIES(${DIR}/Transfers/Petrov-Galerkin-SA)
INCLUDE_DIRECTORIES(${DIR}/Transfers/SemiCoarsen)
INCLUDE_DIRECTORIES(${DIR}/Transfers/Smoothed-Aggregation)
Expand Down Expand Up @@ -397,6 +398,7 @@ TRILINOS_CREATE_CLIENT_TEMPLATE_HEADERS(${DIR}/Transfers/Energy-Minimization NOS
TRILINOS_CREATE_CLIENT_TEMPLATE_HEADERS(${DIR}/Transfers/Energy-Minimization/Solvers NOSIERRABJAM)
TRILINOS_CREATE_CLIENT_TEMPLATE_HEADERS(${DIR}/Transfers/GeneralGeometric NOSIERRABJAM)
TRILINOS_CREATE_CLIENT_TEMPLATE_HEADERS(${DIR}/Transfers/Geometric-Interpolation NOSIERRABJAM)
TRILINOS_CREATE_CLIENT_TEMPLATE_HEADERS(${DIR}/Transfers/Matrix-Free NOSIERRABJAM)
TRILINOS_CREATE_CLIENT_TEMPLATE_HEADERS(${DIR}/Transfers/Petrov-Galerkin-SA NOSIERRABJAM)
TRILINOS_CREATE_CLIENT_TEMPLATE_HEADERS(${DIR}/Transfers/Smoothed-Aggregation NOSIERRABJAM)
TRILINOS_CREATE_CLIENT_TEMPLATE_HEADERS(${DIR}/Transfers/User NOSIERRABJAM)
Expand Down
6 changes: 6 additions & 0 deletions packages/muelu/src/Headers/MueLu_UseShortNamesScalar.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -338,6 +338,12 @@ typedef MueLu::TentativePFactory<Scalar,LocalOrdinal,GlobalOrdinal,Node> Tentati
#ifdef MUELU_TENTATIVEPFACTORY_KOKKOS_SHORT
typedef MueLu::TentativePFactory_kokkos<Scalar,LocalOrdinal,GlobalOrdinal,Node> TentativePFactory_kokkos;
#endif
#ifdef MUELU_MATRIXFREETENTATIVEP_KOKKOS_SHORT
typedef MueLu::MatrixFreeTentativeP_kokkos<Scalar,LocalOrdinal,GlobalOrdinal,Node> MatrixFreeTentativeP_kokkos;
#endif
#ifdef MUELU_MATRIXFREETENTATIVEPFACTORY_KOKKOS_SHORT
typedef MueLu::MatrixFreeTentativePFactory_kokkos<Scalar,LocalOrdinal,GlobalOrdinal,Node> MatrixFreeTentativePFactory_kokkos;
#endif
#ifdef MUELU_THRESHOLDAFILTERFACTORY_SHORT
typedef MueLu::ThresholdAFilterFactory<Scalar,LocalOrdinal,GlobalOrdinal,Node> ThresholdAFilterFactory;
#endif
Expand Down
1 change: 1 addition & 0 deletions packages/muelu/src/Interface/CMakeLists.txt
Original file line number Diff line number Diff line change
Expand Up @@ -38,6 +38,7 @@ INCLUDE_DIRECTORIES(${DIR}/../Transfers/Energy-Minimization)
INCLUDE_DIRECTORIES(${DIR}/../Transfers/Energy-Minimization/Solvers)
INCLUDE_DIRECTORIES(${DIR}/../Transfers/GeneralGeometric)
INCLUDE_DIRECTORIES(${DIR}/../Transfers/Geometric-Interpolation)
INCLUDE_DIRECTORIES(${DIR}/../Transfers/Matrix-Free)
INCLUDE_DIRECTORIES(${DIR}/../Transfers/Smoothed-Aggregation)
INCLUDE_DIRECTORIES(${DIR}/../Transfers/Petrov-Galerkin-SA)
INCLUDE_DIRECTORIES(${DIR}/../Transfers/User)
Expand Down
2 changes: 2 additions & 0 deletions packages/muelu/src/Interface/MueLu_FactoryFactory_decl.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -170,6 +170,7 @@
#include "MueLu_SemiCoarsenPFactory_kokkos.hpp"
#include "MueLu_StructuredAggregationFactory_kokkos.hpp"
#include "MueLu_TentativePFactory_kokkos.hpp"
#include "MueLu_MatrixFreeTentativePFactory_kokkos.hpp"
#include "MueLu_UncoupledAggregationFactory_kokkos.hpp"
#endif

Expand Down Expand Up @@ -320,6 +321,7 @@ namespace MueLu {
if (factoryName == "SemiCoarsenPFactory_kokkos") return Build2<SemiCoarsenPFactory_kokkos> (paramList, factoryMapIn, factoryManagersIn);
if (factoryName == "StructuredAggregationFactory_kokkos") return Build2<StructuredAggregationFactory_kokkos> (paramList, factoryMapIn, factoryManagersIn);
if (factoryName == "TentativePFactory_kokkos") return Build2<TentativePFactory_kokkos> (paramList, factoryMapIn, factoryManagersIn);
if (factoryName == "MatrixFreeTentativePFactory_kokkos") return Build2<MatrixFreeTentativePFactory_kokkos> (paramList, factoryMapIn, factoryManagersIn);
if (factoryName == "UncoupledAggregationFactory_kokkos") return Build2<UncoupledAggregationFactory_kokkos> (paramList, factoryMapIn, factoryManagersIn);
#endif

Expand Down
Original file line number Diff line number Diff line change
@@ -0,0 +1,157 @@
// @HEADER
//
// ***********************************************************************
//
// MueLu: A package for multigrid based preconditioning
// 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
// Jonathan Hu ([email protected])
// Andrey Prokopenko ([email protected])
// Ray Tuminaro ([email protected])
//
// ***********************************************************************
//
// @HEADER
#ifndef MUELU_MATRIXFREETENTATIVEPFACTORY_KOKKOS_DECL_HPP
#define MUELU_MATRIXFREETENTATIVEPFACTORY_KOKKOS_DECL_HPP

#include "MueLu_ConfigDefs.hpp"
#ifdef HAVE_MUELU_KOKKOS_REFACTOR

#include "MueLu_MatrixFreeTentativePFactory_kokkos_fwd.hpp"

#include <KokkosCompat_ClassicNodeAPI_Wrapper.hpp>

#include "Teuchos_ScalarTraits.hpp"

#include "MueLu_Aggregates_kokkos_fwd.hpp"
#include "MueLu_AmalgamationFactory_kokkos_fwd.hpp"
#include "MueLu_AmalgamationInfo_kokkos_fwd.hpp"
#include "MueLu_Level_fwd.hpp"
#include "MueLu_PerfUtils_fwd.hpp"
#include "MueLu_PFactory.hpp"
#include "MueLu_Utilities_kokkos_fwd.hpp"

namespace MueLu {

/*!
@class MatrixFreeTentativePFactory class.
@brief Factory for building the matrix-free tentative restrictor.
Factory for creating the matrix-free tentative restrictor. Nullspace vectors are split across aggregates so that they
have local support. The vectors with local support are factored via LAPACK QR. The Q becomes the
tentative prolongator, and the R becomes the coarse nullspace.
Note that the MatrixFreeTentativePFactory also creates the coarse null space vectors, that is, it serves as generating
factory for the Nullspace vectors on the coarse levels. To feed in the near null space vectors on the finest
level one has to use the NullspaceFactory.
@ingroup MueLuTransferClasses
## Input/output of MatrixFreeTentativePFactory ##
### User parameters of MatrixFreeTentativePFactory ###
Parameter | type | default | master.xml | validated | requested | description
----------|------|---------|:----------:|:---------:|:---------:|------------
Aggregates | Factory | null | | * | * | Generating factory of the aggregates (of type "Aggregates" produced, e.g., by the UncoupledAggregationFactory)
Nullspace | Factory | null | | * | * | Generating factory of the fine nullspace vectors (of type "MultiVector"). In the default case the same instance of the MatrixFreeTentativePFactory is also the generating factory for the null space vectors (on the next coarser levels). Therefore, it is recommended to declare the MatrixFreeTentativePFactory to be the generating factory of the "Nullspace" variable globally using the FactoryManager object! For defining the near null space vectors on the finest level one should use the NullspaceFactory.
UnAmalgamationInfo | Factory | null | | * | * | Generating factory of UnAmalgamationInfo. This data (of type "AmalgamationInfo") is produced by the AmalgamationFactory class. The default option should be fine for most standard cases though.
CoarseMap | Factory | null | | * | * | Generating factory of the coarse map. The map generates the coarse domain map of the prolongation operator. The default choice should be fine as long as you do not wish some offset in the domain map (e.g. for block operators in multiphysics problems). The data is generated by the CoarseMapFactory.
The * in the @c master.xml column denotes that the parameter is defined in the @c master.xml file.<br>
The * in the @c validated column means that the parameter is declared in the list of valid input parameters (see MatrixFreeTentativePFactory::GetValidParameters).<br>
The * in the @c requested column states that the data is requested as input with all dependencies (see MatrixFreeTentativePFactory::DeclareInput).
### Variables provided by MatrixFreeTentativePFactory ###
After MatrixFreeTentativePFactory::Build the following data is available (if requested)
Parameter | generated by | description
----------|--------------|------------
| R | MatrixFreeTentativePFactory | Non-smoothed "tentative" prolongation operator (with piece-wise constant transfer operator basis functions)
| Nullspace | MatrixFreeTentativePFactory | Coarse near null space vectors. Please also check the documentation of the NullspaceFactory for the special dependency tree of the "Nullspace" variable throughout all multigrid levels.
*/
template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
class MatrixFreeTentativePFactory_kokkos;

template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class DeviceType>
class MatrixFreeTentativePFactory_kokkos<Scalar, LocalOrdinal, GlobalOrdinal, Kokkos::Compat::KokkosDeviceWrapperNode<DeviceType> > : public PFactory {
public:
typedef LocalOrdinal local_ordinal_type;
typedef GlobalOrdinal global_ordinal_type;
typedef typename DeviceType::execution_space execution_space;
typedef Kokkos::RangePolicy<local_ordinal_type, execution_space> range_type;
typedef Kokkos::Compat::KokkosDeviceWrapperNode<DeviceType> node_type;
typedef typename Teuchos::ScalarTraits<Scalar>::coordinateType real_type;

private:
// For compatibility
typedef node_type Node;
#undef MUELU_MATRIXFREETENTATIVEPFACTORY_KOKKOS_SHORT
#include "MueLu_UseShortNames.hpp"

public:
//! @name Constructors/Destructors.
//@{

//! Constructor
MatrixFreeTentativePFactory_kokkos() { }

//! Destructor.
virtual ~MatrixFreeTentativePFactory_kokkos() { }
//@}

RCP<const ParameterList> GetValidParameterList() const;

//! Input
//@{

void DeclareInput(Level& fineLevel, Level& coarseLevel) const;

//@}

//! @name Build methods.
//@{

void Build (Level& fineLevel, Level& coarseLevel) const;
void BuildP(Level& fineLevel, Level& coarseLevel) const;

//@}
};

} //namespace MueLu

#define MUELU_MATRIXFREETENTATIVEPFACTORY_KOKKOS_SHORT
#endif // HAVE_MUELU_KOKKOS_REFACTOR
#endif // MUELU_MATRIXFREETENTATIVEPFACTORY_KOKKOS_DECL_HPP
Original file line number Diff line number Diff line change
@@ -0,0 +1,141 @@
// @HEADER
//
// ***********************************************************************
//
// MueLu: A package for multigrid based preconditioning
// 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
// Jonathan Hu ([email protected])
// Andrey Prokopenko ([email protected])
// Ray Tuminaro ([email protected])
//
// ***********************************************************************
//
// @HEADER
#ifndef MUELU_MATRIXFREETENTATIVEPFACTORY_KOKKOS_DEF_HPP
#define MUELU_MATRIXFREETENTATIVEPFACTORY_KOKKOS_DEF_HPP

#ifdef HAVE_MUELU_KOKKOS_REFACTOR

#include "Kokkos_UnorderedMap.hpp"

#include "MueLu_MatrixFreeTentativePFactory_kokkos_decl.hpp"

#include "MueLu_Aggregates_kokkos.hpp"
#include "MueLu_AmalgamationFactory_kokkos.hpp"
#include "MueLu_AmalgamationInfo_kokkos.hpp"
#include "MueLu_CoarseMapFactory_kokkos.hpp"
#include "MueLu_MasterList.hpp"
#include "MueLu_NullspaceFactory_kokkos.hpp"
#include "MueLu_PerfUtils.hpp"
#include "MueLu_Monitor.hpp"
#include "MueLu_MatrixFreeTentativeP_kokkos.hpp"
#include "MueLu_Utilities_kokkos.hpp"

namespace MueLu {

template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class DeviceType>
RCP<const ParameterList> MatrixFreeTentativePFactory_kokkos<Scalar,LocalOrdinal,GlobalOrdinal,Kokkos::Compat::KokkosDeviceWrapperNode<DeviceType>>::GetValidParameterList() const {
RCP<ParameterList> validParamList = rcp(new ParameterList());

validParamList->set< RCP<const FactoryBase> >("A", Teuchos::null, "Generating factory of the matrix A");
validParamList->set< RCP<const FactoryBase> >("Aggregates", Teuchos::null, "Generating factory of the aggregates");
validParamList->set< RCP<const FactoryBase> >("Nullspace", Teuchos::null, "Generating factory of the nullspace");
validParamList->set< RCP<const FactoryBase> >("Scaled Nullspace", Teuchos::null, "Generating factory of the scaled nullspace");
validParamList->set< RCP<const FactoryBase> >("UnAmalgamationInfo", Teuchos::null, "Generating factory of UnAmalgamationInfo");
validParamList->set< RCP<const FactoryBase> >("CoarseMap", Teuchos::null, "Generating factory of the coarse map");
validParamList->set< RCP<const FactoryBase> >("Coordinates", Teuchos::null, "Generating factory of the coordinates");

// Make sure we don't recursively validate options for the matrixmatrix kernels
ParameterList norecurse;
norecurse.disableRecursiveValidation();
validParamList->set<ParameterList> ("matrixmatrix: kernel params", norecurse, "MatrixMatrix kernel parameters");

return validParamList;
}

template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class DeviceType>
void MatrixFreeTentativePFactory_kokkos<Scalar,LocalOrdinal,GlobalOrdinal,Kokkos::Compat::KokkosDeviceWrapperNode<DeviceType>>::DeclareInput(Level& fineLevel, Level& /* coarseLevel */) const {

const ParameterList& pL = GetParameterList();
// NOTE: This guy can only either be 'Nullspace' or 'Scaled Nullspace' or else the validator above will cause issues
std::string nspName = "Nullspace";
if(pL.isParameter("Nullspace name")) nspName = pL.get<std::string>("Nullspace name");

Input(fineLevel, "Aggregates");
Input(fineLevel, nspName);
Input(fineLevel, "UnAmalgamationInfo");
Input(fineLevel, "CoarseMap");
}

template <class Scalar,class LocalOrdinal, class GlobalOrdinal, class DeviceType>
void MatrixFreeTentativePFactory_kokkos<Scalar,LocalOrdinal,GlobalOrdinal,Kokkos::Compat::KokkosDeviceWrapperNode<DeviceType>>::Build(Level& fineLevel, Level& coarseLevel) const {
return BuildP(fineLevel, coarseLevel);
}

template <class Scalar,class LocalOrdinal, class GlobalOrdinal, class DeviceType>
void MatrixFreeTentativePFactory_kokkos<Scalar,LocalOrdinal,GlobalOrdinal,Kokkos::Compat::KokkosDeviceWrapperNode<DeviceType>>::BuildP(Level& fineLevel, Level& coarseLevel) const {
FactoryMonitor m(*this, "Build", coarseLevel);

typedef typename Teuchos::ScalarTraits<Scalar>::coordinateType coordinate_type;
typedef Xpetra::MultiVectorFactory<coordinate_type,LO,GO,NO> RealValuedMultiVectorFactory;
const ParameterList& pL = GetParameterList();
std::string nspName = "Nullspace";
if(pL.isParameter("Nullspace name")) nspName = pL.get<std::string>("Nullspace name");

auto aggregates = Get< RCP<Aggregates_kokkos> > (fineLevel, "Aggregates");
auto amalgInfo = Get< RCP<AmalgamationInfo_kokkos> > (fineLevel, "UnAmalgamationInfo");
auto fineNullspace = Get< RCP<MultiVector> > (fineLevel, nspName);
auto coarseMap = Get< RCP<const Map> > (fineLevel, "CoarseMap");
Teuchos::RCP<const Map> fineMap = fineNullspace->getMap();

// Matrix-free should never run with aggregates that cross processors
if (aggregates->AggregatesCrossProcessors())
TEUCHOS_TEST_FOR_EXCEPTION(true,Exceptions::RuntimeError,"MatrixFreeTentativePFactory does not support aggregates that cross processors!");

size_t NSDim = fineNullspace->getNumVectors();
RCP<MultiVector> coarseNullspace = MultiVectorFactory::Build(coarseMap, NSDim);

Teuchos::RCP<Operator> P = Teuchos::rcp(new MatrixFreeTentativeP_kokkos<Scalar,LocalOrdinal,GlobalOrdinal,Kokkos::Compat::KokkosDeviceWrapperNode<DeviceType>>(coarseMap, fineMap, aggregates));
P->apply(*fineNullspace,*coarseNullspace,Teuchos::TRANS,1.0,0.0); // coarse = alpha*R*fine + beta*coarse

Set(coarseLevel, "Nullspace", coarseNullspace);
Set(coarseLevel, "P", P);
}

} //namespace MueLu

#define MUELU_MATRIXFREETENTATIVEPFACTORY_KOKKOS_SHORT
#endif // HAVE_MUELU_KOKKOS_REFACTOR
#endif // MUELU_MATRIXFREETENTATIVEPFACTORY_KOKKOS_DEF_HPP
Loading

0 comments on commit a2ed36d

Please sign in to comment.