Ifpack2_ILUT_def.hpp

// @HEADER
// *****************************************************************************
//       Ifpack2: Templated Object-Oriented Algebraic Preconditioner Package
//
// Copyright 2009 NTESS and the Ifpack2 contributors.
// SPDX-License-Identifier: BSD-3-Clause
// *****************************************************************************
// @HEADER
 
#ifndef IFPACK2_ILUT_DEF_HPP
#define IFPACK2_ILUT_DEF_HPP
 
#include <type_traits>
#include "Teuchos_TypeNameTraits.hpp"
#include "Teuchos_StandardParameterEntryValidators.hpp"
#include "Teuchos_Time.hpp"
#include "Tpetra_CrsMatrix.hpp"
#include "KokkosSparse_par_ilut.hpp"
 
#include "Ifpack2_Heap.hpp"
#include "Ifpack2_LocalFilter.hpp"
#include "Ifpack2_LocalSparseTriangularSolver.hpp"
#include "Ifpack2_Parameters.hpp"
#include "Ifpack2_Details_getParamTryingTypes.hpp"
 
namespace Ifpack2 {
 
namespace {
 
struct IlutImplType {
  enum Enum {
    Serial,
    PAR_ILUT  
  };
 
  static void loadPLTypeOption(Teuchos::Array<std::string>& type_strs, Teuchos::Array<Enum>& type_enums) {
    type_strs.resize(2);
    type_strs[0] = "serial";
    type_strs[1] = "par_ilut";
    type_enums.resize(2);
    type_enums[0] = Serial;
    type_enums[1] = PAR_ILUT;
  }
};
 
template <class ScalarType>
inline typename Teuchos::ScalarTraits<ScalarType>::magnitudeType
ilutDefaultDropTolerance() {
  typedef Teuchos::ScalarTraits<ScalarType> STS;
  typedef typename STS::magnitudeType magnitude_type;
  typedef Teuchos::ScalarTraits<magnitude_type> STM;
 
  // 1/2.  Hopefully this can be represented in magnitude_type.
  const magnitude_type oneHalf = STM::one() / (STM::one() + STM::one());
 
  // The min ensures that in case magnitude_type has very low
  // precision, we'll at least get some value strictly less than
  // one.
  return std::min(static_cast<magnitude_type>(1000) * STS::magnitude(STS::eps()), oneHalf);
}
 
// Full specialization for ScalarType = double.
// This specialization preserves ILUT's previous default behavior.
template <>
inline Teuchos::ScalarTraits<double>::magnitudeType
ilutDefaultDropTolerance<double>() {
  return 1e-12;
}
 
}  // namespace
 
template <class MatrixType>
ILUT<MatrixType>::ILUT(const Teuchos::RCP<const row_matrix_type>& A)
  : A_(A)
  , Athresh_(Teuchos::ScalarTraits<magnitude_type>::zero())
  , Rthresh_(Teuchos::ScalarTraits<magnitude_type>::one())
  , RelaxValue_(Teuchos::ScalarTraits<magnitude_type>::zero())
  , LevelOfFill_(1.0)
  , DropTolerance_(ilutDefaultDropTolerance<scalar_type>())
  , par_ilut_options_{1, 0., -1, -1, 0.75, false}
  , InitializeTime_(0.0)
  , ComputeTime_(0.0)
  , ApplyTime_(0.0)
  , NumInitialize_(0)
  , NumCompute_(0)
  , NumApply_(0)
  , IsInitialized_(false)
  , IsComputed_(false)
  , useKokkosKernelsParILUT_(false)
 
{
  allocateSolvers();
}
 
template <class MatrixType>
void ILUT<MatrixType>::allocateSolvers() {
  L_solver_ = Teuchos::rcp(new LocalSparseTriangularSolver<row_matrix_type>());
  L_solver_->setObjectLabel("lower");
  U_solver_ = Teuchos::rcp(new LocalSparseTriangularSolver<row_matrix_type>());
  U_solver_->setObjectLabel("upper");
}
 
template <class MatrixType>
void ILUT<MatrixType>::setParameters(const Teuchos::ParameterList& params) {
  using Ifpack2::Details::getParamTryingTypes;
  const char prefix[] = "Ifpack2::ILUT: ";
 
  // Don't actually change the instance variables until we've checked
  // all parameters.  This ensures that setParameters satisfies the
  // strong exception guarantee (i.e., is transactional).
 
  // Parsing implementation type
  IlutImplType::Enum ilutimplType = IlutImplType::Serial;
  do {
    static const char typeName[] = "fact: type";
 
    if (!params.isType<std::string>(typeName)) break;
 
    // Map std::string <-> IlutImplType::Enum.
    Teuchos::Array<std::string> ilutimplTypeStrs;
    Teuchos::Array<IlutImplType::Enum> ilutimplTypeEnums;
    IlutImplType::loadPLTypeOption(ilutimplTypeStrs, ilutimplTypeEnums);
    Teuchos::StringToIntegralParameterEntryValidator<IlutImplType::Enum>
        s2i(ilutimplTypeStrs(), ilutimplTypeEnums(), typeName, false);
 
    ilutimplType = s2i.getIntegralValue(params.get<std::string>(typeName));
  } while (0);
 
  if (ilutimplType == IlutImplType::PAR_ILUT) {
    this->useKokkosKernelsParILUT_ = true;
  } else {
    this->useKokkosKernelsParILUT_ = false;
  }
 
  // Fill level in ILUT is a double, not a magnitude_type, because it
  // depends on LO and GO, not on Scalar.  Also, you can't cast
  // arbitrary magnitude_type (e.g., Sacado::MP::Vector) to double.
  double fillLevel = LevelOfFill_;
  {
    const std::string paramName("fact: ilut level-of-fill");
    TEUCHOS_TEST_FOR_EXCEPTION(
        (params.isParameter(paramName) && this->useKokkosKernelsParILUT_), std::runtime_error,
        "Ifpack2::ILUT: Parameter " << paramName << " is meaningless for algorithm par_ilut.");
    getParamTryingTypes<double, double, float>(fillLevel, params, paramName, prefix);
    TEUCHOS_TEST_FOR_EXCEPTION(fillLevel < 1.0, std::runtime_error,
                               "Ifpack2::ILUT: The \"" << paramName << "\" parameter must be >= "
                                                                       "1.0, but you set it to "
                                                       << fillLevel << ".  For ILUT, the fill level "
                                                                       "means something different than it does for ILU(k).  ILU(0) produces "
                                                                       "factors with the same sparsity structure as the input matrix A. For "
                                                                       "ILUT, level-of-fill = 1.0 will produce factors with nonzeros matching "
                                                                       "the sparsity structure of A. level-of-fill > 1.0 allows for additional "
                                                                       "fill-in.");
  }
 
  magnitude_type absThresh = Athresh_;
  {
    const std::string paramName("fact: absolute threshold");
    getParamTryingTypes<magnitude_type, magnitude_type, double>(absThresh, params, paramName, prefix);
  }
 
  magnitude_type relThresh = Rthresh_;
  {
    const std::string paramName("fact: relative threshold");
    getParamTryingTypes<magnitude_type, magnitude_type, double>(relThresh, params, paramName, prefix);
  }
 
  magnitude_type relaxValue = RelaxValue_;
  {
    const std::string paramName("fact: relax value");
    getParamTryingTypes<magnitude_type, magnitude_type, double>(relaxValue, params, paramName, prefix);
  }
 
  magnitude_type dropTol = DropTolerance_;
  {
    const std::string paramName("fact: drop tolerance");
    getParamTryingTypes<magnitude_type, magnitude_type, double>(dropTol, params, paramName, prefix);
  }
 
  int par_ilut_max_iter                            = 20;
  magnitude_type par_ilut_residual_norm_delta_stop = 1e-2;
  int par_ilut_team_size                           = 0;
  int par_ilut_vector_size                         = 0;
  float par_ilut_fill_in_limit                     = 0.75;
  bool par_ilut_verbose                            = false;
  if (this->useKokkosKernelsParILUT_) {
    par_ilut_max_iter                 = par_ilut_options_.max_iter;
    par_ilut_residual_norm_delta_stop = par_ilut_options_.residual_norm_delta_stop;
    par_ilut_team_size                = par_ilut_options_.team_size;
    par_ilut_vector_size              = par_ilut_options_.vector_size;
    par_ilut_fill_in_limit            = par_ilut_options_.fill_in_limit;
    par_ilut_verbose                  = par_ilut_options_.verbose;
 
    std::string par_ilut_plist_name("parallel ILUT options");
    if (params.isSublist(par_ilut_plist_name)) {
      Teuchos::ParameterList const& par_ilut_plist = params.sublist(par_ilut_plist_name);
 
      std::string paramName("maximum iterations");
      getParamTryingTypes<int, int>(par_ilut_max_iter, par_ilut_plist, paramName, prefix);
 
      paramName = "residual norm delta stop";
      getParamTryingTypes<magnitude_type, magnitude_type, double>(par_ilut_residual_norm_delta_stop, par_ilut_plist, paramName, prefix);
 
      paramName = "team size";
      getParamTryingTypes<int, int>(par_ilut_team_size, par_ilut_plist, paramName, prefix);
 
      paramName = "vector size";
      getParamTryingTypes<int, int>(par_ilut_vector_size, par_ilut_plist, paramName, prefix);
 
      paramName = "fill in limit";
      getParamTryingTypes<float, float, double>(par_ilut_fill_in_limit, par_ilut_plist, paramName, prefix);
 
      paramName = "verbose";
      getParamTryingTypes<bool, bool>(par_ilut_verbose, par_ilut_plist, paramName, prefix);
 
    }  // if (params.isSublist(par_ilut_plist_name))
 
    par_ilut_options_.max_iter                 = par_ilut_max_iter;
    par_ilut_options_.residual_norm_delta_stop = par_ilut_residual_norm_delta_stop;
    par_ilut_options_.team_size                = par_ilut_team_size;
    par_ilut_options_.vector_size              = par_ilut_vector_size;
    par_ilut_options_.fill_in_limit            = par_ilut_fill_in_limit;
    par_ilut_options_.verbose                  = par_ilut_verbose;
 
  }  // if (this->useKokkosKernelsParILUT_)
 
  // Forward to trisolvers.
  L_solver_->setParameters(params);
  U_solver_->setParameters(params);
 
  LevelOfFill_   = fillLevel;
  Athresh_       = absThresh;
  Rthresh_       = relThresh;
  RelaxValue_    = relaxValue;
  DropTolerance_ = dropTol;
}
 
template <class MatrixType>
Teuchos::RCP<const Teuchos::Comm<int> >
ILUT<MatrixType>::getComm() const {
  TEUCHOS_TEST_FOR_EXCEPTION(
      A_.is_null(), std::runtime_error,
      "Ifpack2::ILUT::getComm: "
      "The matrix is null.  Please call setMatrix() with a nonnull input "
      "before calling this method.");
  return A_->getComm();
}
 
template <class MatrixType>
Teuchos::RCP<const typename ILUT<MatrixType>::row_matrix_type>
ILUT<MatrixType>::getMatrix() const {
  return A_;
}
 
template <class MatrixType>
Teuchos::RCP<const typename ILUT<MatrixType>::map_type>
ILUT<MatrixType>::getDomainMap() const {
  TEUCHOS_TEST_FOR_EXCEPTION(
      A_.is_null(), std::runtime_error,
      "Ifpack2::ILUT::getDomainMap: "
      "The matrix is null.  Please call setMatrix() with a nonnull input "
      "before calling this method.");
  return A_->getDomainMap();
}
 
template <class MatrixType>
Teuchos::RCP<const typename ILUT<MatrixType>::map_type>
ILUT<MatrixType>::getRangeMap() const {
  TEUCHOS_TEST_FOR_EXCEPTION(
      A_.is_null(), std::runtime_error,
      "Ifpack2::ILUT::getRangeMap: "
      "The matrix is null.  Please call setMatrix() with a nonnull input "
      "before calling this method.");
  return A_->getRangeMap();
}
 
template <class MatrixType>
bool ILUT<MatrixType>::hasTransposeApply() const {
  return true;
}
 
template <class MatrixType>
int ILUT<MatrixType>::getNumInitialize() const {
  return NumInitialize_;
}
 
template <class MatrixType>
int ILUT<MatrixType>::getNumCompute() const {
  return NumCompute_;
}
 
template <class MatrixType>
int ILUT<MatrixType>::getNumApply() const {
  return NumApply_;
}
 
template <class MatrixType>
double ILUT<MatrixType>::getInitializeTime() const {
  return InitializeTime_;
}
 
template <class MatrixType>
double ILUT<MatrixType>::getComputeTime() const {
  return ComputeTime_;
}
 
template <class MatrixType>
double ILUT<MatrixType>::getApplyTime() const {
  return ApplyTime_;
}
 
template <class MatrixType>
size_t ILUT<MatrixType>::getNodeSmootherComplexity() const {
  TEUCHOS_TEST_FOR_EXCEPTION(
      A_.is_null(), std::runtime_error,
      "Ifpack2::ILUT::getNodeSmootherComplexity: "
      "The input matrix A is null.  Please call setMatrix() with a nonnull "
      "input matrix, then call compute(), before calling this method.");
  // ILUT methods cost roughly one apply + the nnz in the upper+lower triangles
  return A_->getLocalNumEntries() + getLocalNumEntries();
}
 
template <class MatrixType>
global_size_t ILUT<MatrixType>::getGlobalNumEntries() const {
  return L_->getGlobalNumEntries() + U_->getGlobalNumEntries();
}
 
template <class MatrixType>
size_t ILUT<MatrixType>::getLocalNumEntries() const {
  return L_->getLocalNumEntries() + U_->getLocalNumEntries();
}
 
template <class MatrixType>
void ILUT<MatrixType>::setMatrix(const Teuchos::RCP<const row_matrix_type>& A) {
  if (A.getRawPtr() != A_.getRawPtr()) {
    // Check in serial or one-process mode if the matrix is square.
    TEUCHOS_TEST_FOR_EXCEPTION(
        !A.is_null() && A->getComm()->getSize() == 1 &&
            A->getLocalNumRows() != A->getLocalNumCols(),
        std::runtime_error,
        "Ifpack2::ILUT::setMatrix: If A's communicator only "
        "contains one process, then A must be square.  Instead, you provided a "
        "matrix A with "
            << A->getLocalNumRows() << " rows and "
            << A->getLocalNumCols() << " columns.");
 
    // It's legal for A to be null; in that case, you may not call
    // initialize() until calling setMatrix() with a nonnull input.
    // Regardless, setting the matrix invalidates any previous
    // factorization.
    IsInitialized_ = false;
    IsComputed_    = false;
    A_local_       = Teuchos::null;
 
    // The sparse triangular solvers get a triangular factor as their
    // input matrix.  The triangular factors L_ and U_ are getting
    // reset, so we reset the solvers' matrices to null.  Do that
    // before setting L_ and U_ to null, so that latter step actually
    // frees the factors.
    if (!L_solver_.is_null()) {
      L_solver_->setMatrix(Teuchos::null);
    }
    if (!U_solver_.is_null()) {
      U_solver_->setMatrix(Teuchos::null);
    }
 
    L_ = Teuchos::null;
    U_ = Teuchos::null;
    A_ = A;
  }
}
 
template <class MatrixType>
Teuchos::RCP<const typename ILUT<MatrixType>::row_matrix_type>
ILUT<MatrixType>::makeLocalFilter(const Teuchos::RCP<const row_matrix_type>& A) {
  using Teuchos::RCP;
  using Teuchos::rcp;
  using Teuchos::rcp_dynamic_cast;
  using Teuchos::rcp_implicit_cast;
 
  // If A_'s communicator only has one process, or if its column and
  // row Maps are the same, then it is already local, so use it
  // directly.
  if (A->getRowMap()->getComm()->getSize() == 1 ||
      A->getRowMap()->isSameAs(*(A->getColMap()))) {
    return A;
  }
 
  // If A_ is already a LocalFilter, then use it directly.  This
  // should be the case if RILUT is being used through
  // AdditiveSchwarz, for example.
  RCP<const LocalFilter<row_matrix_type> > A_lf_r =
      rcp_dynamic_cast<const LocalFilter<row_matrix_type> >(A);
  if (!A_lf_r.is_null()) {
    return rcp_implicit_cast<const row_matrix_type>(A_lf_r);
  } else {
    // A_'s communicator has more than one process, its row Map and
    // its column Map differ, and A_ is not a LocalFilter.  Thus, we
    // have to wrap it in a LocalFilter.
    return rcp(new LocalFilter<row_matrix_type>(A));
  }
}
 
template <class MatrixType>
void ILUT<MatrixType>::initialize() {
  using Teuchos::Array;
  using Teuchos::RCP;
  using Teuchos::rcp_const_cast;
  Teuchos::Time timer("ILUT::initialize");
  double startTime = timer.wallTime();
  {
    Teuchos::TimeMonitor timeMon(timer);
 
    // Check that the matrix is nonnull.
    TEUCHOS_TEST_FOR_EXCEPTION(
        A_.is_null(), std::runtime_error,
        "Ifpack2::ILUT::initialize: "
        "The matrix to precondition is null.  Please call setMatrix() with a "
        "nonnull input before calling this method.");
 
    // Clear any previous computations.
    IsInitialized_ = false;
    IsComputed_    = false;
    A_local_       = Teuchos::null;
    L_             = Teuchos::null;
    U_             = Teuchos::null;
 
    A_local_ = makeLocalFilter(A_);  // Compute the local filter.
    TEUCHOS_TEST_FOR_EXCEPTION(
        A_local_.is_null(), std::logic_error,
        "Ifpack2::RILUT::initialize: "
        "makeLocalFilter returned null; it failed to compute A_local.  "
        "Please report this bug to the Ifpack2 developers.");
 
    if (this->useKokkosKernelsParILUT_) {
      this->KernelHandle_ = Teuchos::rcp(new kk_handle_type());
      KernelHandle_->create_par_ilut_handle();
      auto par_ilut_handle = KernelHandle_->get_par_ilut_handle();
      par_ilut_handle->set_residual_norm_delta_stop(par_ilut_options_.residual_norm_delta_stop);
      par_ilut_handle->set_team_size(par_ilut_options_.team_size);
      par_ilut_handle->set_vector_size(par_ilut_options_.vector_size);
      par_ilut_handle->set_max_iter(par_ilut_options_.max_iter);
      par_ilut_handle->set_fill_in_limit(par_ilut_options_.fill_in_limit);
      par_ilut_handle->set_verbose(par_ilut_options_.verbose);
      par_ilut_handle->set_async_update(false);
 
      RCP<const crs_matrix_type> A_local_crs = Teuchos::rcp_dynamic_cast<const crs_matrix_type>(A_local_);
      if (A_local_crs.is_null()) {
        // the result is a host-based matrix, which is the same as what happens in RILUK
        local_ordinal_type numRows = A_local_->getLocalNumRows();
        Array<size_t> entriesPerRow(numRows);
        for (local_ordinal_type i = 0; i < numRows; i++) {
          entriesPerRow[i] = A_local_->getNumEntriesInLocalRow(i);
        }
        RCP<crs_matrix_type> A_local_crs_nc =
            rcp(new crs_matrix_type(A_local_->getRowMap(),
                                    A_local_->getColMap(),
                                    entriesPerRow()));
        // copy entries into A_local_crs
        nonconst_local_inds_host_view_type indices("indices", A_local_->getLocalMaxNumRowEntries());
        nonconst_values_host_view_type values("values", A_local_->getLocalMaxNumRowEntries());
        for (local_ordinal_type i = 0; i < numRows; i++) {
          size_t numEntries = 0;
          A_local_->getLocalRowCopy(i, indices, values, numEntries);
          A_local_crs_nc->insertLocalValues(i, numEntries, reinterpret_cast<scalar_type*>(values.data()), indices.data());
        }
        A_local_crs_nc->fillComplete(A_local_->getDomainMap(), A_local_->getRangeMap());
        A_local_crs = rcp_const_cast<const crs_matrix_type>(A_local_crs_nc);
      }
      auto A_local_crs_device = A_local_crs->getLocalMatrixDevice();
 
      // KokkosKernels requires unsigned
      typedef typename Kokkos::View<usize_type*, array_layout, device_type> ulno_row_view_t;
      const int NumMyRows = A_local_crs->getRowMap()->getLocalNumElements();
      L_rowmap_           = ulno_row_view_t("L_row_map", NumMyRows + 1);
      U_rowmap_           = ulno_row_view_t("U_row_map", NumMyRows + 1);
      L_rowmap_orig_      = ulno_row_view_t("L_row_map_orig", NumMyRows + 1);
      U_rowmap_orig_      = ulno_row_view_t("U_row_map_orig", NumMyRows + 1);
 
      KokkosSparse::Experimental::par_ilut_symbolic(KernelHandle_.getRawPtr(),
                                                    A_local_crs_device.graph.row_map, A_local_crs_device.graph.entries,
                                                    L_rowmap_,
                                                    U_rowmap_);
 
      Kokkos::deep_copy(L_rowmap_orig_, L_rowmap_);
      Kokkos::deep_copy(U_rowmap_orig_, U_rowmap_);
    }
 
    IsInitialized_ = true;
    ++NumInitialize_;
  }  // timer scope
  InitializeTime_ += (timer.wallTime() - startTime);
}
 
template <typename ScalarType>
typename Teuchos::ScalarTraits<ScalarType>::magnitudeType
scalar_mag(const ScalarType& s) {
  return Teuchos::ScalarTraits<ScalarType>::magnitude(s);
}
 
template <class MatrixType>
void ILUT<MatrixType>::compute() {
  using Teuchos::Array;
  using Teuchos::ArrayRCP;
  using Teuchos::ArrayView;
  using Teuchos::as;
  using Teuchos::rcp;
  using Teuchos::RCP;
  using Teuchos::rcp_const_cast;
  using Teuchos::reduceAll;
 
  // Don't count initialization in the compute() time.
  if (!isInitialized()) {
    initialize();
  }
 
  Teuchos::Time timer("ILUT::compute");
  double startTime = timer.wallTime();
  {  // Timer scope for timing compute()
    Teuchos::TimeMonitor timeMon(timer, true);
 
    if (!this->useKokkosKernelsParILUT_) {
      //--------------------------------------------------------------------------
      // Ifpack2::ILUT's serial version is a translation of the Aztec ILUT
      // implementation. The Aztec ILUT implementation was written by Ray Tuminaro.
      //
      // This isn't an exact translation of the Aztec ILUT algorithm, for the
      // following reasons:
      // 1. Minor differences result from the fact that Aztec factors a MSR format
      // matrix in place, while the code below factors an input CrsMatrix which
      // remains untouched and stores the resulting factors in separate L and U
      // CrsMatrix objects.
      // Also, the Aztec code begins by shifting the matrix pointers back
      // by one, and the pointer contents back by one, and then using 1-based
      // Fortran-style indexing in the algorithm. This Ifpack2 code uses C-style
      // 0-based indexing throughout.
      // 2. Aztec stores the inverse of the diagonal of U. This Ifpack2 code
      // stores the non-inverted diagonal in U.
      // The triangular solves (in Ifpack2::ILUT::apply()) are performed by
      // calling the Tpetra::CrsMatrix::solve method on the L and U objects, and
      // this requires U to contain the non-inverted diagonal.
      //
      // ABW.
      //--------------------------------------------------------------------------
 
      const scalar_type zero = STS::zero();
      const scalar_type one  = STS::one();
 
      const local_ordinal_type myNumRows = A_local_->getLocalNumRows();
 
      // If this macro is defined, files containing the L and U factors
      // will be written. DON'T CHECK IN THE CODE WITH THIS MACRO ENABLED!!!
      // #define IFPACK2_WRITE_ILUT_FACTORS
#ifdef IFPACK2_WRITE_ILUT_FACTORS
      std::ofstream ofsL("L.ifpack2_ilut.mtx", std::ios::out);
      std::ofstream ofsU("U.ifpack2_ilut.mtx", std::ios::out);
#endif
 
      // Calculate how much fill will be allowed in addition to the
      // space that corresponds to the input matrix entries.
      double local_nnz = static_cast<double>(A_local_->getLocalNumEntries());
      double fill      = ((getLevelOfFill() - 1.0) * local_nnz) / (2 * myNumRows);
 
      // std::ceil gives the smallest integer larger than the argument.
      // this may give a slightly different result than Aztec's fill value in
      // some cases.
      double fill_ceil = std::ceil(fill);
 
      // Similarly to Aztec, we will allow the same amount of fill for each
      // row, half in L and half in U.
      size_type fillL = static_cast<size_type>(fill_ceil);
      size_type fillU = static_cast<size_type>(fill_ceil);
 
      Array<scalar_type> InvDiagU(myNumRows, zero);
 
      Array<Array<local_ordinal_type> > L_tmp_idx(myNumRows);
      Array<Array<scalar_type> > L_tmpv(myNumRows);
      Array<Array<local_ordinal_type> > U_tmp_idx(myNumRows);
      Array<Array<scalar_type> > U_tmpv(myNumRows);
 
      enum { UNUSED,
             ORIG,
             FILL };
      local_ordinal_type max_col = myNumRows;
 
      Array<int> pattern(max_col, UNUSED);
      Array<scalar_type> cur_row(max_col, zero);
      Array<magnitude_type> unorm(max_col);
      magnitude_type rownorm;
      Array<local_ordinal_type> L_cols_heap;
      Array<local_ordinal_type> U_cols;
      Array<local_ordinal_type> L_vals_heap;
      Array<local_ordinal_type> U_vals_heap;
 
      // A comparison object which will be used to create 'heaps' of indices
      // that are ordered according to the corresponding values in the
      // 'cur_row' array.
      greater_indirect<scalar_type, local_ordinal_type> vals_comp(cur_row);
 
      // =================== //
      // start factorization //
      // =================== //
      nonconst_local_inds_host_view_type ColIndicesARCP;
      nonconst_values_host_view_type ColValuesARCP;
      if (!A_local_->supportsRowViews()) {
        const size_t maxnz = A_local_->getLocalMaxNumRowEntries();
        Kokkos::resize(ColIndicesARCP, maxnz);
        Kokkos::resize(ColValuesARCP, maxnz);
      }
 
      for (local_ordinal_type row_i = 0; row_i < myNumRows; ++row_i) {
        local_inds_host_view_type ColIndicesA;
        values_host_view_type ColValuesA;
        size_t RowNnz;
 
        if (A_local_->supportsRowViews()) {
          A_local_->getLocalRowView(row_i, ColIndicesA, ColValuesA);
          RowNnz = ColIndicesA.size();
        } else {
          A_local_->getLocalRowCopy(row_i, ColIndicesARCP, ColValuesARCP, RowNnz);
          ColIndicesA = Kokkos::subview(ColIndicesARCP, std::make_pair((size_t)0, RowNnz));
          ColValuesA  = Kokkos::subview(ColValuesARCP, std::make_pair((size_t)0, RowNnz));
        }
 
        // Always include the diagonal in the U factor. The value should get
        // set in the next loop below.
        U_cols.push_back(row_i);
        cur_row[row_i] = zero;
        pattern[row_i] = ORIG;
 
        size_type L_cols_heaplen = 0;
        rownorm                  = STM::zero();
        for (size_t i = 0; i < RowNnz; ++i) {
          if (ColIndicesA[i] < myNumRows) {
            if (ColIndicesA[i] < row_i) {
              add_to_heap(ColIndicesA[i], L_cols_heap, L_cols_heaplen);
            } else if (ColIndicesA[i] > row_i) {
              U_cols.push_back(ColIndicesA[i]);
            }
 
            cur_row[ColIndicesA[i]] = ColValuesA[i];
            pattern[ColIndicesA[i]] = ORIG;
            rownorm += scalar_mag(ColValuesA[i]);
          }
        }
 
        // Alter the diagonal according to the absolute-threshold and
        // relative-threshold values. If not set, those values default
        // to zero and one respectively.
        const magnitude_type rthresh = getRelativeThreshold();
        const scalar_type& v         = cur_row[row_i];
        cur_row[row_i]               = as<scalar_type>(getAbsoluteThreshold() * IFPACK2_SGN(v)) + rthresh * v;
 
        size_type orig_U_len = U_cols.size();
        RowNnz               = L_cols_heap.size() + orig_U_len;
        rownorm              = getDropTolerance() * rownorm / RowNnz;
 
        // The following while loop corresponds to the 'L30' goto's in Aztec.
        size_type L_vals_heaplen = 0;
        while (L_cols_heaplen > 0) {
          local_ordinal_type row_k = L_cols_heap.front();
 
          scalar_type multiplier  = cur_row[row_k] * InvDiagU[row_k];
          cur_row[row_k]          = multiplier;
          magnitude_type mag_mult = scalar_mag(multiplier);
          if (mag_mult * unorm[row_k] < rownorm) {
            pattern[row_k] = UNUSED;
            rm_heap_root(L_cols_heap, L_cols_heaplen);
            continue;
          }
          if (pattern[row_k] != ORIG) {
            if (L_vals_heaplen < fillL) {
              add_to_heap(row_k, L_vals_heap, L_vals_heaplen, vals_comp);
            } else if (L_vals_heaplen == 0 ||
                       mag_mult < scalar_mag(cur_row[L_vals_heap.front()])) {
              pattern[row_k] = UNUSED;
              rm_heap_root(L_cols_heap, L_cols_heaplen);
              continue;
            } else {
              pattern[L_vals_heap.front()] = UNUSED;
              rm_heap_root(L_vals_heap, L_vals_heaplen, vals_comp);
              add_to_heap(row_k, L_vals_heap, L_vals_heaplen, vals_comp);
            }
          }
 
          /* Reduce current row */
 
          ArrayView<local_ordinal_type> ColIndicesU = U_tmp_idx[row_k]();
          ArrayView<scalar_type> ColValuesU         = U_tmpv[row_k]();
          size_type ColNnzU                         = ColIndicesU.size();
 
          for (size_type j = 0; j < ColNnzU; ++j) {
            if (ColIndicesU[j] > row_k) {
              scalar_type tmp          = multiplier * ColValuesU[j];
              local_ordinal_type col_j = ColIndicesU[j];
              if (pattern[col_j] != UNUSED) {
                cur_row[col_j] -= tmp;
              } else if (scalar_mag(tmp) > rownorm) {
                cur_row[col_j] = -tmp;
                pattern[col_j] = FILL;
                if (col_j > row_i) {
                  U_cols.push_back(col_j);
                } else {
                  add_to_heap(col_j, L_cols_heap, L_cols_heaplen);
                }
              }
            }
          }
 
          rm_heap_root(L_cols_heap, L_cols_heaplen);
        }  // end of while(L_cols_heaplen) loop
 
        // Put indices and values for L into arrays and then into the L_ matrix.
 
        //   first, the original entries from the L section of A:
        for (size_type i = 0; i < (size_type)ColIndicesA.size(); ++i) {
          if (ColIndicesA[i] < row_i) {
            L_tmp_idx[row_i].push_back(ColIndicesA[i]);
            L_tmpv[row_i].push_back(cur_row[ColIndicesA[i]]);
            pattern[ColIndicesA[i]] = UNUSED;
          }
        }
 
        //   next, the L entries resulting from fill:
        for (size_type j = 0; j < L_vals_heaplen; ++j) {
          L_tmp_idx[row_i].push_back(L_vals_heap[j]);
          L_tmpv[row_i].push_back(cur_row[L_vals_heap[j]]);
          pattern[L_vals_heap[j]] = UNUSED;
        }
 
        // L has a one on the diagonal, but we don't explicitly store
        // it.  If we don't store it, then the kernel which performs the
        // triangular solve can assume a unit diagonal, take a short-cut
        // and perform faster.
 
#ifdef IFPACK2_WRITE_ILUT_FACTORS
        for (size_type ii = 0; ii < L_tmp_idx[row_i].size(); ++ii) {
          ofsL << row_i << " " << L_tmp_idx[row_i][ii] << " "
               << L_tmpv[row_i][ii] << std::endl;
        }
#endif
 
        // Pick out the diagonal element, store its reciprocal.
        if (cur_row[row_i] == zero) {
          std::cerr << "Ifpack2::ILUT::Compute: zero pivot encountered! "
                    << "Replacing with rownorm and continuing..."
                    << "(You may need to set the parameter "
                    << "'fact: absolute threshold'.)" << std::endl;
          cur_row[row_i] = rownorm;
        }
        InvDiagU[row_i] = one / cur_row[row_i];
 
        // Non-inverted diagonal is stored for U:
        U_tmp_idx[row_i].push_back(row_i);
        U_tmpv[row_i].push_back(cur_row[row_i]);
        unorm[row_i]   = scalar_mag(cur_row[row_i]);
        pattern[row_i] = UNUSED;
 
        // Now put indices and values for U into arrays and then into the U_ matrix.
        // The first entry in U_cols is the diagonal, which we just handled, so we'll
        // start our loop at j=1.
 
        size_type U_vals_heaplen = 0;
        for (size_type j = 1; j < U_cols.size(); ++j) {
          local_ordinal_type col = U_cols[j];
          if (pattern[col] != ORIG) {
            if (U_vals_heaplen < fillU) {
              add_to_heap(col, U_vals_heap, U_vals_heaplen, vals_comp);
            } else if (U_vals_heaplen != 0 && scalar_mag(cur_row[col]) >
                                                  scalar_mag(cur_row[U_vals_heap.front()])) {
              rm_heap_root(U_vals_heap, U_vals_heaplen, vals_comp);
              add_to_heap(col, U_vals_heap, U_vals_heaplen, vals_comp);
            }
          } else {
            U_tmp_idx[row_i].push_back(col);
            U_tmpv[row_i].push_back(cur_row[col]);
            unorm[row_i] += scalar_mag(cur_row[col]);
          }
          pattern[col] = UNUSED;
        }
 
        for (size_type j = 0; j < U_vals_heaplen; ++j) {
          U_tmp_idx[row_i].push_back(U_vals_heap[j]);
          U_tmpv[row_i].push_back(cur_row[U_vals_heap[j]]);
          unorm[row_i] += scalar_mag(cur_row[U_vals_heap[j]]);
        }
 
        unorm[row_i] /= (orig_U_len + U_vals_heaplen);
 
#ifdef IFPACK2_WRITE_ILUT_FACTORS
        for (int ii = 0; ii < U_tmp_idx[row_i].size(); ++ii) {
          ofsU << row_i << " " << U_tmp_idx[row_i][ii] << " "
               << U_tmpv[row_i][ii] << std::endl;
        }
#endif
 
        L_cols_heap.clear();
        U_cols.clear();
        L_vals_heap.clear();
        U_vals_heap.clear();
      }  // end of for(row_i) loop
 
      // Now allocate and fill the matrices
      Array<size_t> nnzPerRow(myNumRows);
 
      // Make sure to release the old memory for L & U prior to recomputing to
      // avoid bloating the high-water mark.
      L_ = Teuchos::null;
      U_ = Teuchos::null;
      L_solver_->setMatrix(Teuchos::null);
      U_solver_->setMatrix(Teuchos::null);
 
      for (local_ordinal_type row_i = 0; row_i < myNumRows; ++row_i) {
        nnzPerRow[row_i] = L_tmp_idx[row_i].size();
      }
 
      L_ = rcp(new crs_matrix_type(A_local_->getRowMap(), A_local_->getColMap(),
                                   nnzPerRow()));
 
      for (local_ordinal_type row_i = 0; row_i < myNumRows; ++row_i) {
        L_->insertLocalValues(row_i, L_tmp_idx[row_i](), L_tmpv[row_i]());
      }
 
      L_->fillComplete();
 
      for (local_ordinal_type row_i = 0; row_i < myNumRows; ++row_i) {
        nnzPerRow[row_i] = U_tmp_idx[row_i].size();
      }
 
      U_ = rcp(new crs_matrix_type(A_local_->getRowMap(), A_local_->getColMap(),
                                   nnzPerRow()));
 
      for (local_ordinal_type row_i = 0; row_i < myNumRows; ++row_i) {
        U_->insertLocalValues(row_i, U_tmp_idx[row_i](), U_tmpv[row_i]());
      }
 
      U_->fillComplete();
 
      L_solver_->setMatrix(L_);
      L_solver_->initialize();
      L_solver_->compute();
 
      U_solver_->setMatrix(U_);
      U_solver_->initialize();
      U_solver_->compute();
 
    }  // if (!this->useKokkosKernelsParILUT_)
    else {
      // Set L, U rowmaps back to original state. Par_ilut can change them, which invalidates them
      // if compute is called again.
      if (this->isComputed()) {
        Kokkos::resize(L_rowmap_, L_rowmap_orig_.size());
        Kokkos::resize(U_rowmap_, U_rowmap_orig_.size());
        Kokkos::deep_copy(L_rowmap_, L_rowmap_orig_);
        Kokkos::deep_copy(U_rowmap_, U_rowmap_orig_);
      }
 
      RCP<const crs_matrix_type> A_local_crs = Teuchos::rcp_dynamic_cast<const crs_matrix_type>(A_local_);
      {  // Make sure values in A is picked up even in case of pattern reuse
        if (A_local_crs.is_null()) {
          local_ordinal_type numRows = A_local_->getLocalNumRows();
          Array<size_t> entriesPerRow(numRows);
          for (local_ordinal_type i = 0; i < numRows; i++) {
            entriesPerRow[i] = A_local_->getNumEntriesInLocalRow(i);
          }
          RCP<crs_matrix_type> A_local_crs_nc =
              rcp(new crs_matrix_type(A_local_->getRowMap(),
                                      A_local_->getColMap(),
                                      entriesPerRow()));
          // copy entries into A_local_crs
          nonconst_local_inds_host_view_type indices("indices", A_local_->getLocalMaxNumRowEntries());
          nonconst_values_host_view_type values("values", A_local_->getLocalMaxNumRowEntries());
          for (local_ordinal_type i = 0; i < numRows; i++) {
            size_t numEntries = 0;
            A_local_->getLocalRowCopy(i, indices, values, numEntries);
            A_local_crs_nc->insertLocalValues(i, numEntries, reinterpret_cast<scalar_type*>(values.data()), indices.data());
          }
          A_local_crs_nc->fillComplete(A_local_->getDomainMap(), A_local_->getRangeMap());
          A_local_crs = rcp_const_cast<const crs_matrix_type>(A_local_crs_nc);
        }
        auto lclMtx      = A_local_crs->getLocalMatrixDevice();
        A_local_rowmap_  = lclMtx.graph.row_map;
        A_local_entries_ = lclMtx.graph.entries;
        A_local_values_  = lclMtx.values;
      }
 
      // JHU TODO Should allocation of L & U's column (aka entry) and value arrays occur here or in init()?
      auto par_ilut_handle              = KernelHandle_->get_par_ilut_handle();
      auto nnzL                         = par_ilut_handle->get_nnzL();
      static_graph_entries_t L_entries_ = static_graph_entries_t("L_entries", nnzL);
      local_matrix_values_t L_values_   = local_matrix_values_t("L_values", nnzL);
 
      auto nnzU                         = par_ilut_handle->get_nnzU();
      static_graph_entries_t U_entries_ = static_graph_entries_t("U_entries", nnzU);
      local_matrix_values_t U_values_   = local_matrix_values_t("U_values", nnzU);
 
      KokkosSparse::Experimental::par_ilut_numeric(KernelHandle_.getRawPtr(),
                                                   A_local_rowmap_, A_local_entries_, A_local_values_,
                                                   L_rowmap_, L_entries_, L_values_, U_rowmap_, U_entries_, U_values_);
 
      auto L_kokkosCrsGraph = local_graph_device_type(L_entries_, L_rowmap_);
      auto U_kokkosCrsGraph = local_graph_device_type(U_entries_, U_rowmap_);
 
      local_matrix_device_type L_localCrsMatrix_device;
      L_localCrsMatrix_device = local_matrix_device_type("L_Factor_localmatrix",
                                                         A_local_->getLocalNumRows(),
                                                         L_values_,
                                                         L_kokkosCrsGraph);
 
      L_ = rcp(new crs_matrix_type(L_localCrsMatrix_device,
                                   A_local_crs->getRowMap(),
                                   A_local_crs->getColMap(),
                                   A_local_crs->getDomainMap(),
                                   A_local_crs->getRangeMap(),
                                   A_local_crs->getGraph()->getImporter(),
                                   A_local_crs->getGraph()->getExporter()));
 
      local_matrix_device_type U_localCrsMatrix_device;
      U_localCrsMatrix_device = local_matrix_device_type("U_Factor_localmatrix",
                                                         A_local_->getLocalNumRows(),
                                                         U_values_,
                                                         U_kokkosCrsGraph);
 
      U_ = rcp(new crs_matrix_type(U_localCrsMatrix_device,
                                   A_local_crs->getRowMap(),
                                   A_local_crs->getColMap(),
                                   A_local_crs->getDomainMap(),
                                   A_local_crs->getRangeMap(),
                                   A_local_crs->getGraph()->getImporter(),
                                   A_local_crs->getGraph()->getExporter()));
 
      L_solver_->setMatrix(L_);
      L_solver_->compute();  // NOTE: Only do compute if the pointer changed. Otherwise, do nothing
      U_solver_->setMatrix(U_);
      U_solver_->compute();  // NOTE: Only do compute if the pointer changed. Otherwise, do nothing
    }                        // if (!this->useKokkosKernelsParILUT_) ... else ...
 
  }  // Timer scope for timing compute()
  ComputeTime_ += (timer.wallTime() - startTime);
  IsComputed_ = true;
  ++NumCompute_;
}  // compute()
 
template <class MatrixType>
void ILUT<MatrixType>::
    apply(const Tpetra::MultiVector<scalar_type, local_ordinal_type, global_ordinal_type, node_type>& X,
          Tpetra::MultiVector<scalar_type, local_ordinal_type, global_ordinal_type, node_type>& Y,
          Teuchos::ETransp mode,
          scalar_type alpha,
          scalar_type beta) const {
  using Teuchos::RCP;
  using Teuchos::rcp;
  using Teuchos::rcpFromRef;
 
  TEUCHOS_TEST_FOR_EXCEPTION(
      !isComputed(), std::runtime_error,
      "Ifpack2::ILUT::apply: You must call compute() to compute the incomplete "
      "factorization, before calling apply().");
 
  TEUCHOS_TEST_FOR_EXCEPTION(
      X.getNumVectors() != Y.getNumVectors(), std::runtime_error,
      "Ifpack2::ILUT::apply: X and Y must have the same number of columns.  "
      "X has "
          << X.getNumVectors() << " columns, but Y has "
          << Y.getNumVectors() << " columns.");
 
  const scalar_type one  = STS::one();
  const scalar_type zero = STS::zero();
 
  Teuchos::Time timer("ILUT::apply");
  double startTime = timer.wallTime();
  {  // Start timing
    Teuchos::TimeMonitor timeMon(timer, true);
 
    if (alpha == one && beta == zero) {
      if (mode == Teuchos::NO_TRANS) {  // Solve L (U Y) = X for Y.
        // Start by solving L Y = X for Y.
        L_solver_->apply(X, Y, mode);
 
        // Solve U Y = Y.
        U_solver_->apply(Y, Y, mode);
      } else {  // Solve U^P (L^P Y)) = X for Y (where P is * or T).
 
        // Start by solving U^P Y = X for Y.
        U_solver_->apply(X, Y, mode);
 
        // Solve L^P Y = Y.
        L_solver_->apply(Y, Y, mode);
      }
    } else {  // alpha != 1 or beta != 0
      if (alpha == zero) {
        // The special case for beta == 0 ensures that if Y contains Inf
        // or NaN values, we replace them with 0 (following BLAS
        // convention), rather than multiplying them by 0 to get NaN.
        if (beta == zero) {
          Y.putScalar(zero);
        } else {
          Y.scale(beta);
        }
      } else {  // alpha != zero
        MV Y_tmp(Y.getMap(), Y.getNumVectors());
        apply(X, Y_tmp, mode);
        Y.update(alpha, Y_tmp, beta);
      }
    }
  }  // end timing
 
  ++NumApply_;
  ApplyTime_ += (timer.wallTime() - startTime);
}  // apply()
 
template <class MatrixType>
std::string ILUT<MatrixType>::description() const {
  std::ostringstream os;
 
  // Output is a valid YAML dictionary in flow style.  If you don't
  // like everything on a single line, you should call describe()
  // instead.
  os << "\"Ifpack2::ILUT\": {";
  os << "Initialized: " << (isInitialized() ? "true" : "false") << ", "
     << "Computed: " << (isComputed() ? "true" : "false") << ", ";
 
  os << "Level-of-fill: " << getLevelOfFill() << ", "
     << "absolute threshold: " << getAbsoluteThreshold() << ", "
     << "relative threshold: " << getRelativeThreshold() << ", "
     << "relaxation value: " << getRelaxValue() << ", ";
 
  if (A_.is_null()) {
    os << "Matrix: null";
  } else {
    os << "Global matrix dimensions: ["
       << A_->getGlobalNumRows() << ", " << A_->getGlobalNumCols() << "]"
       << ", Global nnz: " << A_->getGlobalNumEntries();
  }
 
  os << "}";
  return os.str();
}
 
template <class MatrixType>
void ILUT<MatrixType>::
    describe(Teuchos::FancyOStream& out,
             const Teuchos::EVerbosityLevel verbLevel) const {
  using std::endl;
  using Teuchos::Comm;
  using Teuchos::OSTab;
  using Teuchos::RCP;
  using Teuchos::TypeNameTraits;
  using Teuchos::VERB_DEFAULT;
  using Teuchos::VERB_EXTREME;
  using Teuchos::VERB_HIGH;
  using Teuchos::VERB_LOW;
  using Teuchos::VERB_MEDIUM;
  using Teuchos::VERB_NONE;
 
  const Teuchos::EVerbosityLevel vl =
      (verbLevel == VERB_DEFAULT) ? VERB_LOW : verbLevel;
  OSTab tab0(out);
 
  if (vl > VERB_NONE) {
    out << "\"Ifpack2::ILUT\":" << endl;
    OSTab tab1(out);
    out << "MatrixType: " << TypeNameTraits<MatrixType>::name() << endl;
    if (this->getObjectLabel() != "") {
      out << "Label: \"" << this->getObjectLabel() << "\"" << endl;
    }
    out << "Initialized: " << (isInitialized() ? "true" : "false")
        << endl
        << "Computed: " << (isComputed() ? "true" : "false")
        << endl
        << "Level of fill: " << getLevelOfFill() << endl
        << "Absolute threshold: " << getAbsoluteThreshold() << endl
        << "Relative threshold: " << getRelativeThreshold() << endl
        << "Relax value: " << getRelaxValue() << endl;
 
    if (isComputed() && vl >= VERB_HIGH) {
      const double fillFraction =
          (double)getGlobalNumEntries() / (double)A_->getGlobalNumEntries();
      const double nnzToRows =
          (double)getGlobalNumEntries() / (double)U_->getGlobalNumRows();
 
      out << "Dimensions of L: [" << L_->getGlobalNumRows() << ", "
          << L_->getGlobalNumRows() << "]" << endl
          << "Dimensions of U: [" << U_->getGlobalNumRows() << ", "
          << U_->getGlobalNumRows() << "]" << endl
          << "Number of nonzeros in factors: " << getGlobalNumEntries() << endl
          << "Fill fraction of factors over A: " << fillFraction << endl
          << "Ratio of nonzeros to rows: " << nnzToRows << endl;
    }
 
    out << "Number of initialize calls: " << getNumInitialize() << endl
        << "Number of compute calls: " << getNumCompute() << endl
        << "Number of apply calls: " << getNumApply() << endl
        << "Total time in seconds for initialize: " << getInitializeTime() << endl
        << "Total time in seconds for compute: " << getComputeTime() << endl
        << "Total time in seconds for apply: " << getApplyTime() << endl;
 
    out << "Local matrix:" << endl;
    A_local_->describe(out, vl);
  }
}
 
}  // namespace Ifpack2
 
// FIXME (mfh 16 Sep 2014) We should really only use RowMatrix here!
// There's no need to instantiate for CrsMatrix too.  All Ifpack2
// preconditioners can and should do dynamic casts if they need a type
// more specific than RowMatrix.
 
#define IFPACK2_ILUT_INSTANT(S, LO, GO, N) \
  template class Ifpack2::ILUT<Tpetra::RowMatrix<S, LO, GO, N> >;
 
#endif /* IFPACK2_ILUT_DEF_HPP */