KokkosSparse::gauss_seidel_symbolic

Defined in header KokkosSparse_gauss_seidel.hpp

template <typename ExecutionSpace, typename KernelHandle, typename lno_row_view_t_, typename lno_nnz_view_t_>
void gauss_seidel_symbolic(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, bool is_graph_symmetric = true);

template <typename KernelHandle, typename lno_row_view_t_, typename lno_nnz_view_t_>
void gauss_seidel_symbolic(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, bool is_graph_symmetric = true);

Performs the symbolic phase of the Gauss-Seidel preconditioner setup. The details of the setup depend on the algorithm selected by the handle.

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.

is_graph_symmetric:

control parameters indicating if the graph of the input matrix is symmetric. This information may be used by the implementation to chose faster implementations.

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>

Example

#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;
}