SPARSE::gauss_seidel

KokkosSparse Gauss-Seidel (& Successive Over-Relaxation)

Applies a parallel Gauss-Seidel preconditioner to a linear system Ax = b. This may be used as a preconditioner for CG/GMRES, as a multigrid smoother or as an iterative solver.

Header File: KokkosSparse_gauss_seidel.hpp

Basic Complete Usage:

handle.create_gs_handle(KokkosSparse::GS_DEFAULT);
KokkosSparse::Experimental::gauss_seidel_symbolic(
  &handle, numRows, numCols, A.graph.row_map, A.graph.entries, graphIsSymmetric);
KokkosSparse::Experimental::gauss_seidel_numeric(
  &handle, numRows, numCols, A.graph.row_map, A.graph.entries, A.values, graphIsSymmetric);
bool firstIter = true;
double relativeResidual = 1.0;
while(relativeResidual > 1e-6)
{
  KokkosSparse::Experimental::forward_sweep_gauss_seidel_apply(
    &handle, numRows, numCols,
    A.graph.row_map, A.graph.entries, A.values,
    x, b, firstIter, firstIter, 1.0, 1);
  firstIter = false;
  /* recompute relativeResidual */
}
handle.destroy_gs_handle();

This is part of a complete working example at kokkos-kernels/example/wiki/sparse/KokkosSparse_wiki_gauss_seidel.cpp.

Handle Setup

Version for point coloring (default) and two stage GS (inner Jacobi-Richardson):

void KokkosKernelsHandle::create_gs_handle(
  KokkosSparse::GSAlgorithm gs_algorithm = KokkosSparse::GS_DEFAULT);

Allowable enum values for gs_algorithm:

KokkosSparse::GS_PERMUTED
KokkosSparse::GS_TEAM
KokkosSparse::GS_TWOSTAGE

For classical (triangular solve) GS:

This variant does a sparse triangular solve on L for forward sweeps, and on R for backward sweeps, rather than implicitly doing forward or backward substitution on rows of A.

handle.create_gs_handle(KokkosSparse::GS_TWOSTAGE);
handle.set_gs_twostage(false, nrows);

where nrows is the number of rows in the matrix that will be passed to symbolic/numeric/apply.

Version for cluster coloring:

This overload of create_gs_handle() always KokkosSparse::GS_CLUSTER as the algorithm.

void KokkosKernelsHandle::create_gs_handle(
    KokkosSparse::ClusteringAlgorithm clusterAlgo, nnz_lno_t verts_per_cluster);

It is strongly recommended to set clusterAlgo = KokkosSparse::CLUSTER_DEFAULT. verts_per_cluster is the coarsening factor. The recommended range is between 4 and 30.

Handle Cleanup

Do this when Gauss-Seidel apply will not be called again. Frees all memory owned by the Gauss-Seidel handle.

void KokkosKernelsHandle::destroy_gs_handle();

Symbolic setup interface

template <typename KernelHandle, typename lno_row_view_t_, typename lno_nnz_view_t_>
void KokkosSparse::Experimental::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);

Parameters:

handle: Pointer to a KokkosKernelsHandle
num_rows: number of rows in the matrix.
num_cols: number of columns in the matrix.
row_map: the CRS row offsets, as a Kokkos::View containing size_type.
entries: the CRS column indices array, as a Kokkos::View containing nnz_lno_t (ordinal).
is_graph_symmetric: true if the graph is known to be symmetric, and false otherwise.

Requirements:

handle has had create_gs_handle() called on it.
num_cols >= num_rows (this is always true of the local part of a distributed square matrix)
row_map.extent(0) == num_rows + 1
row_map(num_rows) == entries.extent(0)
row_map contains elements of type typename KernelHandle::size_type (OK if const)
entries contains elements of type typename KernelHandle::nnz_lno_t (OK if const)

Numeric setup interface

template <typename KernelHandle, typename lno_row_view_t_, typename lno_nnz_view_t_, typename scalar_nnz_view_t_>
void KokkosSparse::Experimental::gauss_seidel_numeric(
    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,
    bool is_graph_symmetric = true);

Parameters:

handle: Pointer to a KokkosKernelsHandle
num_rows: number of rows in the matrix.
num_cols: number of columns in the matrix.
row_map: the CRS row offsets, as a Kokkos::View containing size_type.
entries: the CRS column indices array, as a Kokkos::View containing nnz_lno_t (ordinal).
values: the CRS matrix values, as a Kokkos::View containing nnz_scalar_t.
is_graph_symmetric: true if the graph is known to be symmetric, and false otherwise.

Requirements:

gauss_seidel_symbolic() has already been called with the same handle, and the matrix dimensions/Views are the same as what was passed into symbolic.
row_map, entries, num_rows, num_cols and is_graph_symmetric must all be the same values that were passed into gauss_seidel_symbolic().
row_map(num_rows) == values.extent(0)
values contains elements of type typename KernelHandle::nnz_scalar_t (OK if const)
All other requirements are the same as for gauss_seidel_symbolic

Apply interface (symmetric, forward and backward apply all have the same interface)

template <typename 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 KokkosSparse::Experimental::symmetric_gauss_seidel_apply(
// OR //
void KokkosSparse::Experimental::forward_sweep_gauss_seidel_apply(
// OR //
void KokkosSparse::Experimental::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);

Parameters:

handle: Pointer to a KokkosKernelsHandle
num_rows: number of rows in the matrix.
num_cols: number of columns in the matrix.
row_map: the CRS row offsets, as a Kokkos::View containing size_type.
entries: the CRS column indices array, as a Kokkos::View containing nnz_lno_t (ordinal).
values: the CRS matrix values, as a Kokkos::View containing nnz_scalar_t.
x_lhs_output_vec: the x (unknown) vector as a Kokkos::View containing nnz_scalar_t.
y_rhs_input_vec: the y (right-hand side) vector as a Kokkos::View containing nnz_scalar_t.
init_zero_x_vector: whether to initialize x_lhs_output_vec to 0 as an initial guess. May be used to avoid initializing x at the beginning of a solve.
update_y_vector: true if the right-hand side has changed since the last apply, or if this is the first call to apply with this right-hand side. Should be false when possible to avoid unnecessary work.
omega: the successive over-relaxation damping factor. 1.0 means standard Gauss-Seidel. The optimal value for convergence may be higher or lower than 1.0.
numIter: how many iterations of Gauss-Seidel to perform in this call to apply. For symmetric apply, 1 iteration means a forward followed by a backward sweep. For forward/backward, 1 iteration means one forward/backward sweep.

Requirements:

Both gauss_seidel_symbolic and gauss_seidel_numeric have been called previously with this handle. Apply may be called any number of times as long as the matrix does not change.
If the matrix's dimensions or sparsity pattern change, the handle may be reused but symbolic and numeric must be called again. One exception is that classical GS requires the creation of a new handle if the number of rows changes.
If the matrix's values change, the handle may be reused but numeric must be called again.
num_rows == y_rhs_input_vec.extent(0)
num_cols == x_lhs_output_vec.extent(0)
x_lhs_output_vec.extent(1) == y_rhs_input_vec.extent(1)
x and y have the same rank: x_scalar_view_t::rank == y_scalar_view_t::rank
x and y can be rank 1 (vector) or rank 2 (multivector)