Skip to content

Commit

Permalink
[sor] review updates:
Browse files Browse the repository at this point in the history 
- documentation
- tests
- don't build upper solver if not symmetric

Co-authored-by: Yu-Hsiang M. Tsai <yhmtsai@gmail.com>
  • Loading branch information
@MarcelKoch @yhmtsai
MarcelKoch and yhmtsai committed Aug 9, 2024
1 parent 9cdfcd0 commit 6b7329d
Show file tree
Hide file tree
Showing 6 changed files with 97 additions and 40 deletions.
5 changes: 3 additions & 2 deletions core/preconditioner/sor.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -96,10 +96,11 @@ std::unique_ptr<LinOp> Sor<ValueType, IndexType>::generate_impl(

auto l_trs_factory =
parameters_.l_solver ? parameters_.l_solver : LTrs::build().on(exec);
auto u_trs_factory =
parameters_.u_solver ? parameters_.u_solver : UTrs::build().on(exec);

if (parameters_.symmetric) {
auto u_trs_factory = parameters_.u_solver ? parameters_.u_solver
: UTrs::build().on(exec);

array<index_type> l_row_ptrs{exec, size[0] + 1};
array<index_type> u_row_ptrs{exec, size[0] + 1};
exec->run(make_initialize_row_ptrs_l_u(
Expand Down
65 changes: 49 additions & 16 deletions core/test/config/preconditioner.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -33,6 +33,7 @@ struct PreconditionerConfigTest {
using changed_type = ChangedType;
using default_type = DefaultType;
using preconditioner_config_test = PreconditionerConfigTest;

static pnode::map_type setup_base()
{
return {{"type", pnode{"preconditioner::Ic"}}};
Expand Down Expand Up @@ -325,12 +326,23 @@ struct Sor
config_map["relaxation_factor"] = pnode{0.8};
// float can be cast to double without issues
param.with_relaxation_factor(0.8f);
config_map["l_solver"] = pnode{
{{"type", pnode{"solver::Ir"}}, {"value_type", pnode{"float32"}}}};
param.with_l_solver(DummyIr::build());
config_map["u_solver"] = pnode{
{{"type", pnode{"solver::Ir"}}, {"value_type", pnode{"float32"}}}};
param.with_u_solver(DummyIr::build());
if (from_reg) {
config_map["l_solver"] = pnode{"l_solver"};
param.with_l_solver(
detail::registry_accessor::get_data<gko::LinOpFactory>(
reg, "l_solver"));
config_map["u_solver"] = pnode{"u_solver"};
param.with_u_solver(
detail::registry_accessor::get_data<gko::LinOpFactory>(
reg, "u_solver"));
} else {
config_map["l_solver"] = pnode{{{"type", pnode{"solver::Ir"}},
{"value_type", pnode{"float32"}}}};
param.with_l_solver(Ir::build());
config_map["u_solver"] = pnode{{{"type", pnode{"solver::Ir"}},
{"value_type", pnode{"float32"}}}};
param.with_u_solver(Ir::build());
}
}

template <bool from_reg, typename AnswerType>
Expand All @@ -342,8 +354,13 @@ struct Sor
ASSERT_EQ(res_param.skip_sorting, ans_param.skip_sorting);
ASSERT_EQ(res_param.symmetric, ans_param.symmetric);
ASSERT_EQ(res_param.relaxation_factor, ans_param.relaxation_factor);
ASSERT_EQ(typeid(res_param.l_solver), typeid(ans_param.l_solver));
ASSERT_EQ(typeid(res_param.u_solver), typeid(ans_param.u_solver));
if (from_reg) {
ASSERT_EQ(res_param.l_solver, ans_param.l_solver);
ASSERT_EQ(res_param.u_solver, ans_param.u_solver);
} else {
ASSERT_EQ(typeid(res_param.l_solver), typeid(ans_param.l_solver));
ASSERT_EQ(typeid(res_param.u_solver), typeid(ans_param.u_solver));
}
}
};

Expand Down Expand Up @@ -372,12 +389,23 @@ struct GaussSeidel
param.with_skip_sorting(true);
config_map["symmetric"] = pnode{true};
param.with_symmetric(true);
config_map["l_solver"] = pnode{
{{"type", pnode{"solver::Ir"}}, {"value_type", pnode{"float32"}}}};
param.with_l_solver(DummyIr::build());
config_map["u_solver"] = pnode{
{{"type", pnode{"solver::Ir"}}, {"value_type", pnode{"float32"}}}};
param.with_u_solver(DummyIr::build());
if (from_reg) {
config_map["l_solver"] = pnode{"l_solver"};
param.with_l_solver(
detail::registry_accessor::get_data<gko::LinOpFactory>(
reg, "l_solver"));
config_map["u_solver"] = pnode{"u_solver"};
param.with_u_solver(
detail::registry_accessor::get_data<gko::LinOpFactory>(
reg, "u_solver"));
} else {
config_map["l_solver"] = pnode{{{"type", pnode{"solver::Ir"}},
{"value_type", pnode{"float32"}}}};
param.with_l_solver(Ir::build());
config_map["u_solver"] = pnode{{{"type", pnode{"solver::Ir"}},
{"value_type", pnode{"float32"}}}};
param.with_u_solver(Ir::build());
}
}

template <bool from_reg, typename AnswerType>
Expand All @@ -388,8 +416,13 @@ struct GaussSeidel

ASSERT_EQ(res_param.skip_sorting, ans_param.skip_sorting);
ASSERT_EQ(res_param.symmetric, ans_param.symmetric);
ASSERT_EQ(typeid(res_param.l_solver), typeid(ans_param.l_solver));
ASSERT_EQ(typeid(res_param.u_solver), typeid(ans_param.u_solver));
if (from_reg) {
ASSERT_EQ(res_param.l_solver, ans_param.l_solver);
ASSERT_EQ(res_param.u_solver, ans_param.u_solver);
} else {
ASSERT_EQ(typeid(res_param.l_solver), typeid(ans_param.l_solver));
ASSERT_EQ(typeid(res_param.u_solver), typeid(ans_param.u_solver));
}
}
};

Expand Down
16 changes: 13 additions & 3 deletions include/ginkgo/core/preconditioner/gauss_seidel.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -22,7 +22,8 @@ namespace preconditioner {
/**
* This class generates the Gauss-Seidel preconditioner.
*
* This is the special case of $\omega = 1$ of the (S)SOR preconditioner.
* This is the special case of the relaxation factor $\omega = 1$ of the (S)SOR
* preconditioner.
*
* @see Sor
*/
Expand All @@ -49,12 +50,12 @@ class GaussSeidel
// determines if Gauss-Seidel or symmetric Gauss-Seidel should be used
bool GKO_FACTORY_PARAMETER_SCALAR(symmetric, false);

// factory for the lower triangular factor solver
// factory for the lower triangular factor solver, defaults to LowerTrs
std::shared_ptr<const LinOpFactory> GKO_DEFERRED_FACTORY_PARAMETER(
l_solver);

// factory for the upper triangular factor solver, unused if symmetric
// is false
// is false, defaults to UpperTrs
std::shared_ptr<const LinOpFactory> GKO_DEFERRED_FACTORY_PARAMETER(
u_solver);
};
Expand Down Expand Up @@ -84,6 +85,15 @@ class GaussSeidel
/** Creates a new parameter_type to set up the factory. */
static parameters_type build() { return {}; }

/**
* Creates a new parameter_type based on an input configuration.
*
* @param config The parameter tree for the Gauss-Seidel.
* @param context The context for the parsing.
* @param td_for_child The type_descriptor used for child configurations.
*
* @return a parameter_type object with the settings from config
*/
static parameters_type parse(
const config::pnode& config, const config::registry& context,
const config::type_descriptor& td_for_child =
Expand Down
17 changes: 15 additions & 2 deletions include/ginkgo/core/preconditioner/sor.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -36,6 +36,8 @@ namespace preconditioner {
* M = \frac{1}{\omega (2 - \omega)} (D + \omega L) D^{-1} (D + \omega U) ,
* \quad 0 < \omega < 2.
* $$
* A detailed description can be found in Iterative Methods for Sparse Linear
* Systems (Y. Saad) ch. 4.1.
*
* This class is a factory, which will only generate the preconditioner. The
* resulting LinOp will represent the application of $M^{-1}$.
Expand Down Expand Up @@ -69,12 +71,12 @@ class Sor
remove_complex<value_type> GKO_FACTORY_PARAMETER_SCALAR(
relaxation_factor, remove_complex<value_type>(1.2));

// factory for the lower triangular factor solver
// factory for the lower triangular factor solver, defaults to LowerTrs
std::shared_ptr<const LinOpFactory> GKO_DEFERRED_FACTORY_PARAMETER(
l_solver);

// factory for the upper triangular factor solver, unused if symmetric
// is false
// is false, defaults to UpperTrs
std::shared_ptr<const LinOpFactory> GKO_DEFERRED_FACTORY_PARAMETER(
u_solver);
};
Expand Down Expand Up @@ -104,6 +106,15 @@ class Sor
/** Creates a new parameter_type to set up the factory. */
static parameters_type build() { return {}; }

/**
* Creates a new parameter_type based on an input configuration.
*
* @param config The parameter tree for the SOR.
* @param context The context for the parsing.
* @param td_for_child The type_descriptor used for child configurations.
*
* @return a parameter_type object with the settings from config
*/
static parameters_type parse(
const config::pnode& config, const config::registry& context,
const config::type_descriptor& td_for_child =
Expand All @@ -124,6 +135,8 @@ class Sor
private:
parameters_type parameters_;
};


} // namespace preconditioner
} // namespace gko

Expand Down
2 changes: 1 addition & 1 deletion reference/CMakeLists.txt
View file
Original file line number Diff line number Diff line change
Expand Up @@ -43,7 +43,7 @@ target_sources(ginkgo_reference
matrix/sparsity_csr_kernels.cpp
multigrid/pgm_kernels.cpp
preconditioner/batch_jacobi_kernels.cpp
preconditioner/sor_kernels.cpp
preconditioner/sor_kernels.cpp
preconditioner/isai_kernels.cpp
preconditioner/jacobi_kernels.cpp
reorder/rcm_kernels.cpp
Expand Down
32 changes: 16 additions & 16 deletions reference/test/preconditioner/sor_kernels.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -80,15 +80,15 @@ TYPED_TEST(Sor, CanInitializeLFactorWithWeight)
result->scale(
gko::initialize<gko::matrix::Dense<value_type>>({0.0}, this->exec));
std::shared_ptr<csr_type> expected_l =
gko::initialize<csr_type>({{1, 0, 0, 0, 0},
{-2, 1, 0, 0, 0},
{0, -5, 1, 0, 0},
{-3, 0, 0, 1, 0},
{-4, 0, -6, -7, 1}},
gko::initialize<csr_type>({{2 * this->diag_value, 0, 0, 0, 0},
{-2, 2 * this->diag_value, 0, 0, 0},
{0, -5, 2 * this->diag_value, 0, 0},
{-3, 0, 0, 2 * this->diag_value, 0},
{-4, 0, -6, -7, 2 * this->diag_value}},
this->exec);

gko::kernels::reference::sor::initialize_weighted_l(
this->exec, this->mtx.get(), this->diag_value, result.get());
this->exec, this->mtx.get(), 0.5f, result.get());

GKO_ASSERT_MTX_NEAR(result, expected_l, r<value_type>::value);
}
Expand Down Expand Up @@ -122,15 +122,16 @@ TYPED_TEST(Sor, CanInitializeLAndUFactorWithWeight)
gko::initialize<gko::matrix::Dense<value_type>>({0.0}, this->exec));
result_u->scale(
gko::initialize<gko::matrix::Dense<value_type>>({0.0}, this->exec));
auto diag_weight = static_cast<gko::remove_complex<value_type>>(
1.0 / (2 - this->diag_value));
auto off_diag_weight = this->diag_value * diag_weight;
auto factor = static_cast<gko::remove_complex<value_type>>(0.5);
auto diag_weight =
static_cast<gko::remove_complex<value_type>>(1.0 / (2 - factor));
auto off_diag_weight = factor * diag_weight;
std::shared_ptr<csr_type> expected_l =
gko::initialize<csr_type>({{1, 0, 0, 0, 0},
{-2, 1, 0, 0, 0},
{0, -5, 1, 0, 0},
{-3, 0, 0, 1, 0},
{-4, 0, -6, -7, 1}},
gko::initialize<csr_type>({{2 * this->diag_value, 0, 0, 0, 0},
{-2, 2 * this->diag_value, 0, 0, 0},
{0, -5, 2 * this->diag_value, 0, 0},
{-3, 0, 0, 2 * this->diag_value, 0},
{-4, 0, -6, -7, 2 * this->diag_value}},
this->exec);
std::shared_ptr<csr_type> expected_u = gko::initialize<csr_type>(
{{this->diag_value * diag_weight, 2 * off_diag_weight, 0,
Expand All @@ -142,8 +143,7 @@ TYPED_TEST(Sor, CanInitializeLAndUFactorWithWeight)
this->exec);

gko::kernels::reference::sor::initialize_weighted_l_u(
this->exec, this->mtx.get(), this->diag_value, result_l.get(),
result_u.get());
this->exec, this->mtx.get(), factor, result_l.get(), result_u.get());

GKO_ASSERT_MTX_NEAR(result_l, expected_l, r<value_type>::value);
GKO_ASSERT_MTX_NEAR(result_u, expected_u, r<value_type>::value);
Expand Down

0 comments on commit 6b7329d

Please sign in to comment.