Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Minres solver #975

Merged
merged 3 commits into from
Mar 3, 2025
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
18 changes: 11 additions & 7 deletions benchmark/solver/solver_common.hpp
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 Expand Up @@ -200,6 +201,9 @@ std::unique_ptr<gko::LinOpFactory> generate_solver(
gko::solver::Gmres<etype>::build().with_krylov_dim(
FLAGS_gmres_restart),
exec, precond, max_iters);
} else if (description == "minres") {
return add_criteria_precond_finalize<gko::solver::Minres<etype>>(
exec, precond, max_iters);
} else if (description == "lower_trs") {
return gko::solver::LowerTrs<etype>::build()
.with_num_rhs(FLAGS_nrhs)
Expand Down
1 change: 1 addition & 0 deletions common/unified/CMakeLists.txt
Original file line number Diff line number Diff line change
Expand Up @@ -30,6 +30,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)
193 changes: 193 additions & 0 deletions common/unified/solver/minres_kernels.cpp
Original file line number Diff line number Diff line change
@@ -0,0 +1,193 @@
// SPDX-FileCopyrightText: 2017 - 2025 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>
GKO_INLINE GKO_ATTRIBUTES void swap(T& a, T& b)
{
T 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]);

Check warning on line 57 in common/unified/solver/minres_kernels.cpp

View check run for this annotation

Codecov / codecov/patch

common/unified/solver/minres_kernels.cpp#L57

Added line #L57 was not covered by tests
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));
Comment on lines +93 to +94
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Is this to protect against overflow?

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Honestly, I'm not sure about that. This is just the same computation as for the gmres/cb_gmres.

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,
matrix::Dense<ValueType>* 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];

Check warning on line 119 in common/unified/solver/minres_kernels.cpp

View check run for this annotation

Codecov / codecov/patch

common/unified/solver/minres_kernels.cpp#L118-L119

Added lines #L118 - L119 were not covered by tests
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;

Check warning on line 123 in common/unified/solver/minres_kernels.cpp

View check run for this annotation

Codecov / codecov/patch

common/unified/solver/minres_kernels.cpp#L123

Added line #L123 was not covered by tests
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] = sin[col] * sin[col] * tau[col];

Check warning on line 132 in common/unified/solver/minres_kernels.cpp

View check run for this annotation

Codecov / codecov/patch

common/unified/solver/minres_kernels.cpp#L132

Added line #L132 was not covered by tests
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
Original file line number Diff line number Diff line change
Expand Up @@ -119,6 +119,7 @@ target_sources(
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
1 change: 1 addition & 0 deletions core/config/config_helper.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -49,6 +49,7 @@ enum class LinOpFactoryType : int {
Gcr,
Gmres,
CbGmres,
Minres,
Direct,
LowerTrs,
UpperTrs,
Expand Down
3 changes: 2 additions & 1 deletion core/config/registry.cpp
Original file line number Diff line number Diff line change
@@ -1,4 +1,4 @@
// SPDX-FileCopyrightText: 2017 - 2024 The Ginkgo authors
// SPDX-FileCopyrightText: 2017 - 2025 The Ginkgo authors
//
// SPDX-License-Identifier: BSD-3-Clause

Expand Down Expand Up @@ -29,6 +29,7 @@ configuration_map generate_config_map()
{"solver::Gcr", parse<LinOpFactoryType::Gcr>},
{"solver::Gmres", parse<LinOpFactoryType::Gmres>},
{"solver::CbGmres", parse<LinOpFactoryType::CbGmres>},
{"solver::Minres", parse<LinOpFactoryType::Minres>},
{"solver::Direct", parse<LinOpFactoryType::Direct>},
{"solver::LowerTrs", parse<LinOpFactoryType::LowerTrs>},
{"solver::UpperTrs", parse<LinOpFactoryType::UpperTrs>},
Expand Down
4 changes: 3 additions & 1 deletion core/config/solver_config.cpp
Original file line number Diff line number Diff line change
@@ -1,4 +1,4 @@
// SPDX-FileCopyrightText: 2017 - 2024 The Ginkgo authors
// SPDX-FileCopyrightText: 2017 - 2025 The Ginkgo authors
//
// SPDX-License-Identifier: BSD-3-Clause

Expand All @@ -18,6 +18,7 @@
#include <ginkgo/core/solver/gmres.hpp>
#include <ginkgo/core/solver/idr.hpp>
#include <ginkgo/core/solver/ir.hpp>
#include <ginkgo/core/solver/minres.hpp>
#include <ginkgo/core/solver/multigrid.hpp>
#include <ginkgo/core/solver/triangular.hpp>

Expand All @@ -40,6 +41,7 @@ GKO_PARSE_VALUE_TYPE(Idr, gko::solver::Idr);
GKO_PARSE_VALUE_TYPE(Gcr, gko::solver::Gcr);
GKO_PARSE_VALUE_TYPE(Gmres, gko::solver::Gmres);
GKO_PARSE_VALUE_TYPE_BASE(CbGmres, gko::solver::CbGmres);
GKO_PARSE_VALUE_TYPE(Minres, gko::solver::Minres);
GKO_PARSE_VALUE_AND_INDEX_TYPE(Direct, gko::experimental::solver::Direct);
GKO_PARSE_VALUE_AND_INDEX_TYPE(LowerTrs, gko::solver::LowerTrs);
GKO_PARSE_VALUE_AND_INDEX_TYPE(UpperTrs, gko::solver::UpperTrs);
Expand Down
12 changes: 12 additions & 0 deletions core/device_hooks/common_kernels.inc.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -70,6 +70,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 @@ -696,6 +697,17 @@
} // 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);

Check warning on line 705 in core/device_hooks/common_kernels.inc.cpp

View check run for this annotation

Codecov / codecov/patch

core/device_hooks/common_kernels.inc.cpp#L703-L705

Added lines #L703 - L705 were not covered by tests


} // namespace minres


namespace sparsity_csr {


Expand Down
Loading
Loading