Skip to content

Commit

Permalink
add minres implementation
Browse files Browse the repository at this point in the history 
based on 'Iterative Methods for Solving Linear Systems' (https://doi.org/10.1137/1.9781611970937) and 'Iterative Methods for Singular Linear Equations and Least-Squares Problems' (PhD thesis, Stanford University)

Co-authored-by: Aditya Kashi <aditya.kashi@kit.edu>
Co-authored-by: Hartwig Anzt <hanzt@icl.utk.edu>
  • Loading branch information
@MarcelKoch @Slaedr @hartwiganzt
3 people committed Dec 15, 2023
1 parent a47745f commit 2af34a1
Show file tree
Hide file tree
Showing 16 changed files with 1,640 additions and 10 deletions.
15 changes: 8 additions & 7 deletions benchmark/solver/solver_common.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -34,13 +34,14 @@ DEFINE_bool(
rel_residual, false,
"Use relative residual instead of residual reduction stopping criterion");

DEFINE_string(solvers, "cg",
"A comma-separated list of solvers to run. "
"Supported values are: bicgstab, bicg, cb_gmres_keep, "
"cb_gmres_reduce1, cb_gmres_reduce2, cb_gmres_integer, "
"cb_gmres_ireduce1, cb_gmres_ireduce2, cg, cgs, fcg, gmres, idr, "
"lower_trs, upper_trs, spd_direct, symm_direct, "
"near_symm_direct, direct, overhead");
DEFINE_string(
solvers, "cg",
"A comma-separated list of solvers to run. "
"Supported values are: bicgstab, bicg, cb_gmres_keep, "
"cb_gmres_reduce1, cb_gmres_reduce2, cb_gmres_integer, "
"cb_gmres_ireduce1, cb_gmres_ireduce2, cg, cgs, direct, fcg, gmres, idr, "
"lower_trs, minres, near_symm_direct, upper_trs, spd_direct, symm_direct, "
"overhead");

DEFINE_uint32(
nrhs, 1,
Expand Down
3 changes: 2 additions & 1 deletion common/unified/CMakeLists.txt
View file
Original file line number Diff line number Diff line change
Expand Up @@ -28,6 +28,7 @@ set(UNIFIED_SOURCES
solver/gcr_kernels.cpp
solver/gmres_kernels.cpp
solver/ir_kernels.cpp
solver/minres_kernels.cpp
)
list(TRANSFORM UNIFIED_SOURCES PREPEND ${CMAKE_CURRENT_SOURCE_DIR}/)
set(GKO_UNIFIED_COMMON_SOURCES ${UNIFIED_SOURCES} PARENT_SCOPE)
set(GKO_UNIFIED_COMMON_SOURCES ${UNIFIED_SOURCES} PARENT_SCOPE)
195 changes: 195 additions & 0 deletions common/unified/solver/minres_kernels.cpp
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,195 @@
// SPDX-FileCopyrightText: 2017-2023 The Ginkgo authors
//
// SPDX-License-Identifier: BSD-3-Clause

#include "core/solver/minres_kernels.hpp"


#include <ginkgo/core/base/executor.hpp>


#include "common/unified/base/kernel_launch_solver.hpp"


namespace gko {
namespace kernels {
namespace GKO_DEVICE_NAMESPACE {
/**
* @brief The Minres solver namespace.
*
* @ingroup minres
*/
namespace minres {
namespace detail {


template <typename T, typename U>
GKO_INLINE GKO_ATTRIBUTES void swap(T& a, U& b)
{
U tmp{b};
b = a;
a = tmp;
}


} // namespace detail


template <typename ValueType>
void initialize(
std::shared_ptr<const DefaultExecutor> exec,
const matrix::Dense<ValueType>* r, matrix::Dense<ValueType>* z,
matrix::Dense<ValueType>* p, matrix::Dense<ValueType>* p_prev,
matrix::Dense<ValueType>* q, matrix::Dense<ValueType>* q_prev,
matrix::Dense<ValueType>* v, matrix::Dense<ValueType>* beta,
matrix::Dense<ValueType>* gamma, matrix::Dense<ValueType>* delta,
matrix::Dense<ValueType>* cos_prev, matrix::Dense<ValueType>* cos,
matrix::Dense<ValueType>* sin_prev, matrix::Dense<ValueType>* sin,
matrix::Dense<ValueType>* eta_next, matrix::Dense<ValueType>* eta,
array<stopping_status>* stop_status)
{
run_kernel(
exec,
[] GKO_KERNEL(auto col, auto beta, auto gamma, auto delta,
auto cos_prev, auto cos, auto sin_prev, auto sin,
auto eta_next, auto eta, auto stop) {
delta[col] = gamma[col] = cos_prev[col] = sin_prev[col] = sin[col] =
zero(*delta);
cos[col] = one(*delta);
eta_next[col] = eta[col] = beta[col] = sqrt(beta[col]);
stop[col].reset();
},
beta->get_num_stored_elements(), row_vector(beta), row_vector(gamma),
row_vector(delta), row_vector(cos_prev), row_vector(cos),
row_vector(sin_prev), row_vector(sin), row_vector(eta_next),
row_vector(eta), *stop_status);

run_kernel_solver(
exec,
[] GKO_KERNEL(auto row, auto col, auto r, auto z, auto p, auto p_prev,
auto q, auto q_prev, auto v, auto beta, auto stop) {
q(row, col) = safe_divide(r(row, col), beta[col]);
z(row, col) = safe_divide(z(row, col), beta[col]);
p(row, col) = p_prev(row, col) = q_prev(row, col) = v(row, col) =
zero(p(row, col));
},
r->get_size(), r->get_stride(), default_stride(r), default_stride(z),
default_stride(p), default_stride(p_prev), default_stride(q),
default_stride(q_prev), default_stride(v), row_vector(beta),
*stop_status);
}

GKO_INSTANTIATE_FOR_EACH_VALUE_TYPE(GKO_DECLARE_MINRES_INITIALIZE_KERNEL);


template <typename ValueType>
GKO_KERNEL void update_givens_rotation(ValueType& alpha, const ValueType& beta,
ValueType& cos, ValueType& sin)
{
if (alpha == zero(alpha)) {
cos = zero(cos);
sin = one(sin);
} else {
const auto scale = abs(alpha) + abs(beta);
const auto hypotenuse =
scale * sqrt(abs(alpha / scale) * abs(alpha / scale) +
abs(beta / scale) * abs(beta / scale));
cos = conj(alpha) / hypotenuse;
sin = conj(beta) / hypotenuse;
}
alpha = cos * alpha + sin * beta;
}


template <typename ValueType>
void step_1(std::shared_ptr<const DefaultExecutor> exec,
matrix::Dense<ValueType>* alpha, matrix::Dense<ValueType>* beta,
matrix::Dense<ValueType>* gamma, matrix::Dense<ValueType>* delta,
matrix::Dense<ValueType>* cos_prev, matrix::Dense<ValueType>* cos,
matrix::Dense<ValueType>* sin_prev, matrix::Dense<ValueType>* sin,
matrix::Dense<ValueType>* eta, matrix::Dense<ValueType>* eta_next,
typename matrix::Dense<ValueType>::absolute_type* tau,
const array<stopping_status>* stop_status)
{
run_kernel(
exec,
[] GKO_KERNEL(auto col, auto alpha, auto beta, auto gamma, auto delta,
auto cos_prev, auto cos, auto sin_prev, auto sin,
auto eta_next, auto eta, auto tau, auto stop) {
if (!stop[col].has_stopped()) {
beta[col] = sqrt(beta[col]);
delta[col] = sin_prev[col] * gamma[col];
const auto tmp_d = gamma[col];
const auto tmp_a = alpha[col];
gamma[col] =
cos_prev[col] * cos[col] * tmp_d + sin[col] * tmp_a;
alpha[col] =
-conj(sin[col]) * cos_prev[col] * tmp_d + cos[col] * tmp_a;

detail::swap(cos[col], cos_prev[col]);
detail::swap(sin[col], sin_prev[col]);
update_givens_rotation(alpha[col], beta[col], cos[col],
sin[col]);

tau[col] = abs(sin[col]) * tau[col];
eta[col] = eta_next[col];
eta_next[col] = -conj(sin[col]) * eta[col];
}
},
alpha->get_num_stored_elements(), row_vector(alpha), row_vector(beta),
row_vector(gamma), row_vector(delta), row_vector(cos_prev),
row_vector(cos), row_vector(sin_prev), row_vector(sin),
row_vector(eta_next), row_vector(eta), row_vector(tau), *stop_status);
}

GKO_INSTANTIATE_FOR_EACH_VALUE_TYPE(GKO_DECLARE_MINRES_STEP_1_KERNEL);


template <typename ValueType>
void step_2(std::shared_ptr<const DefaultExecutor> exec,
matrix::Dense<ValueType>* x, matrix::Dense<ValueType>* p,
const matrix::Dense<ValueType>* p_prev, matrix::Dense<ValueType>* z,
const matrix::Dense<ValueType>* z_tilde,
matrix::Dense<ValueType>* q, matrix::Dense<ValueType>* q_prev,
matrix::Dense<ValueType>* v, const matrix::Dense<ValueType>* alpha,
const matrix::Dense<ValueType>* beta,
const matrix::Dense<ValueType>* gamma,
const matrix::Dense<ValueType>* delta,
const matrix::Dense<ValueType>* cos,
const matrix::Dense<ValueType>* eta,
const array<stopping_status>* stop_status)
{
run_kernel_solver(
exec,
[] GKO_KERNEL(auto row, auto col, auto x, auto p, auto p_prev, auto q,
auto q_prev, auto v, auto z, auto z_tilde, auto alpha,
auto beta, auto gamma, auto delta, auto cos, auto eta,
auto stop) {
if (!stop[col].has_stopped()) {
p(row, col) =
safe_divide(z(row, col) - gamma[col] * p_prev(row, col) -
delta[col] * p(row, col),
alpha[col]);
x(row, col) = x(row, col) + cos[col] * eta[col] * p(row, col);

q_prev(row, col) = v(row, col);
const auto tmp = q(row, col);
z(row, col) = safe_divide(z_tilde(row, col), beta[col]);
q(row, col) = safe_divide(v(row, col), beta[col]);
v(row, col) = tmp * beta[col];
}
},
x->get_size(), p->get_stride(), x, default_stride(p),
default_stride(p_prev), default_stride(q), default_stride(q_prev),
default_stride(v), default_stride(z), default_stride(z_tilde),
row_vector(alpha), row_vector(beta), row_vector(gamma),
row_vector(delta), row_vector(cos), row_vector(eta), *stop_status);
}

GKO_INSTANTIATE_FOR_EACH_VALUE_TYPE(GKO_DECLARE_MINRES_STEP_2_KERNEL);


} // namespace minres
} // namespace GKO_DEVICE_NAMESPACE
} // namespace kernels
} // namespace gko
1 change: 1 addition & 0 deletions core/CMakeLists.txt
View file
Original file line number Diff line number Diff line change
Expand Up @@ -78,6 +78,7 @@ target_sources(ginkgo
solver/gmres.cpp
solver/idr.cpp
solver/ir.cpp
solver/minres.cpp
solver/lower_trs.cpp
solver/multigrid.cpp
solver/upper_trs.cpp
Expand Down
12 changes: 12 additions & 0 deletions core/device_hooks/common_kernels.inc.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -61,6 +61,7 @@
#include "core/solver/idr_kernels.hpp"
#include "core/solver/ir_kernels.hpp"
#include "core/solver/lower_trs_kernels.hpp"
#include "core/solver/minres_kernels.hpp"
#include "core/solver/multigrid_kernels.hpp"
#include "core/solver/upper_trs_kernels.hpp"
#include "core/stop/criterion_kernels.hpp"
Expand Down Expand Up @@ -562,6 +563,17 @@ GKO_STUB_NON_COMPLEX_VALUE_TYPE(GKO_DECLARE_MULTIGRID_KCYCLE_CHECK_STOP_KERNEL);
} // namespace multigrid


namespace minres {


GKO_STUB_VALUE_TYPE(GKO_DECLARE_MINRES_INITIALIZE_KERNEL);
GKO_STUB_VALUE_TYPE(GKO_DECLARE_MINRES_STEP_1_KERNEL);
GKO_STUB_VALUE_TYPE(GKO_DECLARE_MINRES_STEP_2_KERNEL);


} // namespace minres


namespace sparsity_csr {


Expand Down
Loading

0 comments on commit 2af34a1

Please sign in to comment.