Merge pull request #1661 from jgfouca/jgfouca/par_ilut_test

Adds a better parilut test with gmres

blas/src/KokkosBlas3_trsm.hpp

@@ -42,8 +42,9 @@ namespace KokkosBlas {
 ///                   "L" or "l" indicates matrix A lower part is stored, the
 ///                   other part is not referenced
 /// \param trans [in] "N" or "n" for non-transpose, "T" or "t" for transpose,
-/// "C" or "c" for conjugate transpose. \param diag  [in] "U" or "u" indicates
-/// the diagonal of A is assumed to be unit
+/// "C" or "c" for conjugate transpose.
+/// \param diag  [in] "U" or "u" indicates the diagonal of A is assumed to be
+/// unit
 //                    "N" or "n" indicated the diagonal of A is assumed to be
 //                    non-unit
 /// \param alpha [in] Input coefficient used for multiplication with B

common/src/KokkosKernels_Error.hpp

@@ -53,4 +53,60 @@ inline void hip_internal_safe_call(hipError_t e, const char *name,
 }  // namespace Impl
 }  // namespace KokkosKernels
 
+/*
+ * Asserts and error checking macros/functions.
+ *
+ * KK_KERNEL** are for error checking within kokkos kernels.
+ *
+ * Any check with "assert" in the name is disabled for release builds
+ *
+ * For _MSG checks, the msg argument can contain '<<' if not a kernel check.
+ *
+ * This code is adapted from EKAT/src/ekat/ekat_assert.hpp
+ */
+
+// Internal do not call directly
+#define IMPL_THROW(condition, msg, exception_type)                      \
+  do {                                                                  \
+    if (!(condition)) {                                                 \
+      std::stringstream _ss_;                                           \
+      _ss_ << __FILE__ << ":" << __LINE__ << ": FAIL:\n" << #condition; \
+      _ss_ << "\n" << msg;                                              \
+      throw exception_type(_ss_.str());                                 \
+    }                                                                   \
+  } while (0)
+
+// SYCL cannot printf like the other backends quite yet
+#define IMPL_KERNEL_THROW(condition, msg)                                   \
+  do {                                                                      \
+    if (!(condition)) {                                                     \
+      KOKKOS_IMPL_DO_NOT_USE_PRINTF("KERNEL CHECK FAILED:\n   %s\n   %s\n", \
+                                    #condition, msg);                       \
+      Kokkos::abort("");                                                    \
+    }                                                                       \
+  } while (0)
+
+#ifndef NDEBUG
+#define KK_ASSERT(condition) IMPL_THROW(condition, "", std::logic_error)
+#define KK_ASSERT_MSG(condition, msg) \
+  IMPL_THROW(condition, msg, std::logic_error)
+#define KK_KERNEL_ASSERT(condition) IMPL_KERNEL_THROW(condition, "")
+#define KK_KERNEL_ASSERT_MSG(condition, msg) IMPL_KERNEL_THROW(condition, msg)
+#else
+#define KK_ASSERT(condition) ((void)(0))
+#define KK_ASSERT_MSG(condition, msg) ((void)(0))
+#define KK_KERNEL_ASSERT(condition) ((void)(0))
+#define KK_KERNEL_ASSERT_MSG(condition, msg) ((void)(0))
+#endif
+
+#define KK_REQUIRE(condition) IMPL_THROW(condition, "", std::logic_error)
+#define KK_REQUIRE_MSG(condition, msg) \
+  IMPL_THROW(condition, msg, std::logic_error)
+
+#define KK_KERNEL_REQUIRE(condition) IMPL_KERNEL_THROW(condition, "")
+#define KK_KERNEL_REQUIRE_MSG(condition, msg) IMPL_KERNEL_THROW(condition, msg)
+
+#define KK_ERROR_MSG(msg) KK_REQUIRE_MSG(false, msg)
+#define KK_KERNEL_ERROR_MSG(msg) KK_KERNEL_REQUIRE_MSG(false, msg)
+
 #endif  // KOKKOSKERNELS_ERROR_HPP

sparse/impl/KokkosSparse_par_ilut_numeric_impl.hpp

@@ -331,11 +331,17 @@ struct IlutWrap {
             const auto out_val = lpu_col == col_idx ? lpu_val : r_val / diag;
             // store output entries
             if (row_idx >= col_idx) {
+              KK_KERNEL_ASSERT_MSG(
+                  l_new_nnz < L_new_row_map(row_idx + 1),
+                  "add_candidates: Overflowed L_new, is your A matrix sorted?");
               L_new_entries(l_new_nnz) = col_idx;
               L_new_values(l_new_nnz)  = row_idx == col_idx ? 1. : out_val;
               ++l_new_nnz;
             }
             if (row_idx <= col_idx) {
+              KK_KERNEL_ASSERT_MSG(
+                  u_new_nnz < U_new_row_map(row_idx + 1),
+                  "add_candidates: Overflowed U_new, is your A matrix sorted?");
               U_new_entries(u_new_nnz) = col_idx;
               U_new_values(u_new_nnz)  = out_val;
               ++u_new_nnz;

sparse/src/KokkosSparse_LUPrec.hpp

@@ -0,0 +1,135 @@
+//@HEADER
+// ************************************************************************
+//
+//                        Kokkos v. 4.0
+//       Copyright (2022) National Technology & Engineering
+//               Solutions of Sandia, LLC (NTESS).
+//
+// Under the terms of Contract DE-NA0003525 with NTESS,
+// the U.S. Government retains certain rights in this software.
+//
+// Part of Kokkos, under the Apache License v2.0 with LLVM Exceptions.
+// See https://kokkos.org/LICENSE for license information.
+// SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception
+//
+// ************************************************************************
+//@HEADER
+
+/// @file KokkosSparse_LUPrec.hpp
+
+#ifndef KK_LU_PREC_HPP
+#define KK_LU_PREC_HPP
+
+#include <KokkosSparse_Preconditioner.hpp>
+#include <Kokkos_Core.hpp>
+#include <KokkosSparse_spmv.hpp>
+#include <KokkosSparse_sptrsv.hpp>
+
+namespace KokkosSparse {
+namespace Experimental {
+
+/// \class LUPrec
+/// \brief  This class is for applying LU preconditioning.
+///         It takes L and U and the apply method returns U^inv L^inv x
+/// \tparam CRS the CRS type of L and U
+///
+/// Preconditioner provides the following methods
+///   - initialize() Does nothing; members initialized upon object construction.
+///   - isInitialized() returns true
+///   - compute() Does nothing; members initialized upon object construction.
+///   - isComputed() returns true
+///
+template <class CRS, class KernelHandle>
+class LUPrec : public KokkosSparse::Experimental::Preconditioner<CRS> {
+ public:
+  using ScalarType = typename std::remove_const<typename CRS::value_type>::type;
+  using EXSP       = typename CRS::execution_space;
+  using karith     = typename Kokkos::ArithTraits<ScalarType>;
+  using View1d = typename Kokkos::View<ScalarType *, typename CRS::device_type>;
+
+ private:
+  // trsm takes host views
+  CRS _L, _U;
+  View1d _tmp, _tmp2;
+  mutable KernelHandle _khL;
+  mutable KernelHandle _khU;
+
+ public:
+  //! Constructor:
+  template <class CRSArg>
+  LUPrec(const CRSArg &L, const CRSArg &U)
+      : _L(L),
+        _U(U),
+        _tmp("LUPrec::_tmp", L.numRows()),
+        _tmp2("LUPrec::_tmp", L.numRows()),
+        _khL(),
+        _khU() {
+    KK_REQUIRE_MSG(L.numRows() == U.numRows(),
+                   "LUPrec: L.numRows() != U.numRows()");
+
+    _khL.create_sptrsv_handle(SPTRSVAlgorithm::SEQLVLSCHD_TP1, L.numRows(),
+                              true);
+    _khU.create_sptrsv_handle(SPTRSVAlgorithm::SEQLVLSCHD_TP1, U.numRows(),
+                              false);
+  }
+
+  //! Destructor.
+  virtual ~LUPrec() {
+    _khL.destroy_sptrsv_handle();
+    _khU.destroy_sptrsv_handle();
+  }
+
+  ///// \brief Apply the preconditioner to X, putting the result in Y.
+  /////
+  ///// \tparam XViewType Input vector, as a 1-D Kokkos::View
+  ///// \tparam YViewType Output vector, as a nonconst 1-D Kokkos::View
+  /////
+  ///// \param transM [in] Not used.
+  ///// \param alpha [in] Not used
+  ///// \param beta [in] Not used.
+  /////
+  ///// It takes L and U and the stores U^inv L^inv X in Y
+  //
+  virtual void apply(const Kokkos::View<const ScalarType *, EXSP> &X,
+                     const Kokkos::View<ScalarType *, EXSP> &Y,
+                     const char transM[] = "N",
+                     ScalarType alpha    = karith::one(),
+                     ScalarType beta     = karith::zero()) const {
+    // tmp = trsv(L, x); //Apply L^inv to x
+    // y = trsv(U, tmp); //Apply U^inv to tmp
+
+    KK_REQUIRE_MSG(transM[0] == NoTranspose[0],
+                   "LUPrec::apply only supports 'N' for transM");
+
+    sptrsv_symbolic(&_khL, _L.graph.row_map, _L.graph.entries);
+    sptrsv_solve(&_khL, _L.graph.row_map, _L.graph.entries, _L.values, X, _tmp);
+
+    sptrsv_symbolic(&_khU, _U.graph.row_map, _U.graph.entries);
+    sptrsv_solve(&_khU, _U.graph.row_map, _U.graph.entries, _U.values, _tmp,
+                 _tmp2);
+
+    KokkosBlas::axpby(alpha, _tmp2, beta, Y);
+  }
+  //@}
+
+  //! Set this preconditioner's parameters.
+  void setParameters() {}
+
+  void initialize() {}
+
+  //! True if the preconditioner has been successfully initialized, else false.
+  bool isInitialized() const { return true; }
+
+  void compute() {}
+
+  //! True if the preconditioner has been successfully computed, else false.
+  bool isComputed() const { return true; }
+
+  //! True if the preconditioner implements a transpose operator apply.
+  bool hasTransposeApply() const { return true; }
+};
+
+}  // namespace Experimental
+}  // End namespace KokkosSparse
+
+#endif

sparse/src/KokkosSparse_MatrixPrec.hpp

@@ -34,13 +34,9 @@ namespace Experimental {
 ///         already has a matrix representation of their
 ///         preconditioner M.  The class applies an
 ///         SpMV with M as the preconditioning step.
-/// \tparam ScalarType Type of the matrix's entries
-/// \tparam Layout Kokkos layout of vectors X and Y to which
-///         the preconditioner is applied
-/// \tparam EXSP Execution space for the preconditioner apply
-/// \tparam Ordinal Type of the matrix's indices;
+/// \tparam CRS the type of compressed matrix
 ///
-/// Preconditioner provides the following methods
+/// MatrixPrec provides the following methods
 ///   - initialize() Does nothing; Matrix initialized upon object construction.
 ///   - isInitialized() returns true
 ///   - compute() Does nothing; Matrix initialized upon object construction.
@@ -49,10 +45,7 @@ namespace Experimental {
 template <class CRS>
 class MatrixPrec : public KokkosSparse::Experimental::Preconditioner<CRS> {
  private:
-  CRS A;
-
-  bool isInitialized_ = true;
-  bool isComputed_    = true;
+  CRS _A;
 
  public:
   using ScalarType = typename std::remove_const<typename CRS::value_type>::type;
@@ -61,7 +54,7 @@ class MatrixPrec : public KokkosSparse::Experimental::Preconditioner<CRS> {
 
   //! Constructor:
   template <class CRSArg>
-  MatrixPrec(const CRSArg &mat) : A(mat) {}
+  MatrixPrec(const CRSArg &mat) : _A(mat) {}
 
   //! Destructor.
   virtual ~MatrixPrec() {}
@@ -87,7 +80,7 @@ class MatrixPrec : public KokkosSparse::Experimental::Preconditioner<CRS> {
                      const char transM[] = "N",
                      ScalarType alpha    = karith::one(),
                      ScalarType beta     = karith::zero()) const {
-    KokkosSparse::spmv(transM, alpha, A, X, beta, Y);
+    KokkosSparse::spmv(transM, alpha, _A, X, beta, Y);
   }
   //@}
 
@@ -97,12 +90,12 @@ class MatrixPrec : public KokkosSparse::Experimental::Preconditioner<CRS> {
   void initialize() {}
 
   //! True if the preconditioner has been successfully initialized, else false.
-  bool isInitialized() const { return isInitialized_; }
+  bool isInitialized() const { return true; }
 
   void compute() {}
 
   //! True if the preconditioner has been successfully computed, else false.
-  bool isComputed() const { return isComputed_; }
+  bool isComputed() const { return true; }
 
   //! True if the preconditioner implements a transpose operator apply.
   bool hasTransposeApply() const { return true; }

sparse/src/KokkosSparse_Preconditioner.hpp

@@ -26,14 +26,9 @@
 namespace KokkosSparse {
 namespace Experimental {
 
-/// @file KokkosSparse_Preconditioner.hpp
 /// \class Preconditioner
 /// \brief Interface for KokkosKernels preconditioners
-/// \tparam ScalarType Type of the matrix's entries
-/// \tparam Layout Kokkos layout of vectors X and Y to which
-///         the preconditioner is applied
-/// \tparam EXSP Execution space for the preconditioner apply
-/// \tparam Ordinal Type of the matrix's indices;
+/// \tparam CRS Type of the compressed matrix
 ///
 /// Preconditioner provides the following methods
 ///   - initialize() performs all operations based on the graph of the

sparse/src/KokkosSparse_par_ilut.hpp

@@ -348,6 +348,11 @@ void par_ilut_numeric(KernelHandle* handle, ARowMapType& A_rowmap,
     KokkosKernels::Impl::throw_runtime_exception(os.str());
   }
 
+  KK_REQUIRE_MSG(KokkosSparse::Impl::isCrsGraphSorted(L_rowmap, L_entries),
+                 "L is not sorted");
+  KK_REQUIRE_MSG(KokkosSparse::Impl::isCrsGraphSorted(U_rowmap, U_entries),
+                 "U is not sorted");
+
   using c_size_t   = typename KernelHandle::const_size_type;
   using c_lno_t    = typename KernelHandle::const_nnz_lno_t;
   using c_scalar_t = typename KernelHandle::const_nnz_scalar_t;

sparse/unit_test/Test_Sparse_gmres.hpp

@@ -86,14 +86,13 @@ void run_test_gmres() {
   ViewVectorType Wj("Wj", n);  // For checking residuals at end.
   ViewVectorType B(Kokkos::view_alloc(Kokkos::WithoutInitializing, "B"),
                    n);  // right-hand side vec
+  // Make rhs ones so that results are repeatable:
+  Kokkos::deep_copy(B, 1.0);
 
   gmres_handle->set_verbose(verbose);
 
   // Test CGS2
   {
-    // Make rhs ones so that results are repeatable:
-    Kokkos::deep_copy(B, 1.0);
-
     gmres(&kh, A, B, X);
 
     // Double check residuals at end of solve:
@@ -137,12 +136,12 @@ void run_test_gmres() {
     gmres_handle->set_verbose(verbose);
 
     // Make precond
-    auto myPrec = new KokkosSparse::Experimental::MatrixPrec<sp_matrix_type>(A);
+    KokkosSparse::Experimental::MatrixPrec<sp_matrix_type> myPrec(A);
 
     // reset X for next gmres call
     Kokkos::deep_copy(X, 0.0);
 
-    gmres(&kh, A, B, X, myPrec);
+    gmres(&kh, A, B, X, &myPrec);
 
     // Double check residuals at end of solve:
     float_t nrmB = KokkosBlas::nrm2(B);
@@ -154,8 +153,6 @@ void run_test_gmres() {
 
     EXPECT_LT(endRes, gmres_handle->get_tol());
     EXPECT_EQ(conv_flag, GMRESHandle::Flag::Conv);
-
-    delete myPrec;
   }
 }