KokkosSparse::backward_sweep_gauss_seidel_apply

Defined in header KokkosSparse_gauss_seidel.hpp

template <class ExecutionSpace, KokkosSparse::SparseMatrixFormat format = KokkosSparse::SparseMatrixFormat::CRS,
          class KernelHandle, typename lno_row_view_t_, typename lno_nnz_view_t_, typename scalar_nnz_view_t_,
          typename x_scalar_view_t, typename y_scalar_view_t>
void backward_sweep_gauss_seidel_apply(const ExecutionSpace &space, KernelHandle *handle,
                                       typename KernelHandle::const_nnz_lno_t num_rows,
                                       typename KernelHandle::const_nnz_lno_t num_cols, lno_row_view_t_ row_map,
                                       lno_nnz_view_t_ entries, scalar_nnz_view_t_ values,
                                       x_scalar_view_t x_lhs_output_vec, y_scalar_view_t y_rhs_input_vec,
                                       bool init_zero_x_vector, bool update_y_vector,
                                       typename KernelHandle::nnz_scalar_t omega, int numIter);

template <KokkosSparse::SparseMatrixFormat format = KokkosSparse::SparseMatrixFormat::CRS, class KernelHandle,
          typename lno_row_view_t_, typename lno_nnz_view_t_, typename scalar_nnz_view_t_, typename x_scalar_view_t,
          typename y_scalar_view_t>
void backward_sweep_gauss_seidel_apply(KernelHandle *handle, typename KernelHandle::const_nnz_lno_t num_rows,
                                       typename KernelHandle::const_nnz_lno_t num_cols, lno_row_view_t_ row_map,
                                       lno_nnz_view_t_ entries, scalar_nnz_view_t_ values,
                                       x_scalar_view_t x_lhs_output_vec, y_scalar_view_t y_rhs_input_vec,
                                       bool init_zero_x_vector, bool update_y_vector,
                                       typename KernelHandle::nnz_scalar_t omega, int numIter)

Apply the Gauss-Seidel preconditioner to a linear system of equations \(Ax=b\).

  1. Perform symbolic operations to setup the Gauss-Seidel/SOR preconditioner.

  2. Same as 1. but uses the resources of the default execution space instance associated with KernelHandle::HandleExecSpace.

Parameters

space:

execution space instance describing the resources available to run the kernels.

handle:

an instance of KokkosKernels::KokkosKernelsHandle from which a gs_handle will be used to extract control parameters.

num_rows, num_cols:

the number of rows and columns of the input matrix.

rowmap:

row map of the input matrix.

entries:

column indices of the input matrix.

values:

values of the input matrix.

x_lhs_output_vec, y_rhs_input_vec:

left hand side and right hand side vectors of the linear system.

init_zero_x_vector:

whether x_lhs_output_vec should be zero initialized.

update_y_vector:

whether y_rhs_input_vec has changed since the last call to one of the Gauss-Seidel apply functions using handle.

omega:

damping parameter for successive over-relaxations.

numIter:

the number of preconditioner iteration to perform on the linear system.

Type Requirements

  • The types of the input parameters must be consistent with the types defined by the KernelHandle:

    • std::is_same_v<typename KernelHandle::const_size_type, typename lno_row_view_t_::const_value_type>

    • std::is_same_v<typename KernelHandle::const_nnz_lno_t, typename lno_nnz_view_t_::const_value_type>

    • std::is_same_v<typename KernelHandle::const_nnz_scalar_t, typename scalar_nnz_view_t_::const_value_type>

    • std::is_same_v<typename KernelHandle::const_nnz_scalar_t, typename y_scalar_view_t::const_value_type>

    • std::is_same_v<typename KernelHandle::nnz_scalar_t, typename x_scalar_view_t::value_type>

  • The views describing the matrix, rowmap, entries and values, should not have a Kokkos::LayoutStride

    • !std::is_same_v<typename lno_row_view_t_::array_layout, Kokkos::LayoutStride>

    • !std::is_same_v<typename lno_nnz_view_t_::array_layout, Kokkos::LayoutStride>

    • !std::is_same_v<typename scalar_nnz_view_t_::array_layout, Kokkos::LayoutStride>

Note

The non-strided requirement in this function should probably be checked on earlier during the symbolic and numeric phases? Also are we really happy with something not layout stride? Should we check on is_contiguous instead?

Example

Note that while the example below uses forward sweeps, calling the backward sweep API is almost identical.

#include "Kokkos_Core.hpp"
#include "KokkosKernels_default_types.hpp"
#include "KokkosKernels_Handle.hpp"
#include "KokkosKernels_IOUtils.hpp"
#include "KokkosSparse_IOUtils.hpp"
#include "KokkosSparse_spmv.hpp"
#include "KokkosSparse_CrsMatrix.hpp"
#include "KokkosSparse_gauss_seidel.hpp"
#include "KokkosBlas1_nrm2.hpp"

// Parallel Gauss-Seidel Preconditioner/Smoother
//  -Uses graph coloring to find independent row sets,
//   and applies GS to each set in parallel
//  -Here, use to solve a diagonally dominant linear system directly.

// Helper to print out colors in the shape of the grid
int main() {
  using Scalar    = KokkosKernels::default_scalar;
  using Mag       = Kokkos::ArithTraits<Scalar>::mag_type;
  using Ordinal   = KokkosKernels::default_lno_t;
  using Offset    = KokkosKernels::default_size_type;
  using ExecSpace = Kokkos::DefaultExecutionSpace;
  using MemSpace  = typename ExecSpace::memory_space;
  using Device    = Kokkos::Device<ExecSpace, MemSpace>;
  using Handle    = KokkosKernels::Experimental::KokkosKernelsHandle<Offset, Ordinal, KokkosKernels::default_scalar,
                                                                  ExecSpace, MemSpace, MemSpace>;
  using Matrix    = KokkosSparse::CrsMatrix<Scalar, Ordinal, Device, void, Offset>;
  using Vector    = typename Matrix::values_type;
  constexpr Ordinal numRows = 10000;
  const Scalar one          = Kokkos::ArithTraits<Scalar>::one();
  const Mag magOne          = Kokkos::ArithTraits<Mag>::one();
  // Solve tolerance
  const Mag tolerance = 1e-6 * magOne;
  Kokkos::initialize();
  {
    // Generate a square, strictly diagonally dominant, but nonsymmetric matrix
    // on which Gauss-Seidel should converge. Get approx. 20 entries per row
    // Diagonals are 2x the absolute sum of all other entries.
    Offset nnz = numRows * 20;
    Matrix A = KokkosSparse::Impl::kk_generate_diagonally_dominant_sparse_matrix<Matrix>(numRows, numRows, nnz, 2, 100,
                                                                                         1.05 * one);
    std::cout << "Generated a matrix with " << numRows << " rows/cols, and " << nnz << " entries.\n";
    // Create a kernel handle, then a Gauss-Seidel handle with the default
    // algorithm
    Handle handle;
    handle.create_gs_handle(KokkosSparse::GS_DEFAULT);
    // Set up Gauss-Seidel for the graph (matrix sparsity pattern)
    KokkosSparse::gauss_seidel_symbolic(&handle, numRows, numRows, A.graph.row_map, A.graph.entries, false);
    // Set up Gauss-Seidel for the matrix values (numeric)
    // Another matrix with the same sparsity pattern could re-use the handle and
    // symbolic phase, and only call numeric.
    KokkosSparse::gauss_seidel_numeric(&handle, numRows, numRows, A.graph.row_map, A.graph.entries, A.values, false);
    // Now, preconditioner is ready to use. Set up an unknown vector
    // (uninitialized) and randomized right-hand-side vector.
    Vector x(Kokkos::view_alloc(Kokkos::WithoutInitializing, "x"), numRows);
    Vector b(Kokkos::view_alloc(Kokkos::WithoutInitializing, "b"), numRows);
    Vector res(Kokkos::view_alloc(Kokkos::WithoutInitializing, "res"), numRows);
    auto bHost = Kokkos::create_mirror_view(b);
    for (Ordinal i = 0; i < numRows; i++) bHost(i) = 3 * ((one * rand()) / RAND_MAX);
    Kokkos::deep_copy(b, bHost);
    // Measure initial residual norm ||Ax - b||, where x is 0
    Mag initialRes    = KokkosBlas::nrm2(b);
    Mag scaledResNorm = magOne;
    bool firstIter    = true;
    // Iterate until reaching the tolerance
    int numIters = 0;
    while (scaledResNorm > tolerance) {
      // Run one sweep of forward Gauss-Seidel (SOR with omega = 1.0)
      // If this is the first iteration, tell apply:
      //  * to zero out x (it was uninitialized)
      //  * that b has changed since the previous apply (since there was no
      //  previous apply)
      KokkosSparse::forward_sweep_gauss_seidel_apply(&handle, numRows, numRows, A.graph.row_map, A.graph.entries,
                                                     A.values, x, b, firstIter, firstIter, one, 1);
      firstIter = false;
      // Now, compute the new residual norm using SPMV
      Kokkos::deep_copy(res, b);
      // Compute res := Ax - res (since res is now equal to b, this is Ax - b)
      KokkosSparse::spmv("N", one, A, x, -one, res);
      // Recompute the scaled norm
      scaledResNorm = KokkosBlas::nrm2(res) / initialRes;
      numIters++;
      std::cout << "Iteration " << numIters << " scaled residual norm: " << scaledResNorm << '\n';
    }
    std::cout << "SUCCESS: converged in " << numIters << " iterations.\n";
  }
  Kokkos::finalize();
  return 0;
}