Skip to content

Commit

Permalink
Merge Add FGMRES
Browse files Browse the repository at this point in the history 
This PR adds the Flexible GMRES (FGMRES).
It can be used from GMRES with `with_flexible(true)`.

Related PR: #1244
  • Loading branch information
@yhmtsai
yhmtsai authored Feb 8, 2023
2 parents 89b8cec + a3b9699 commit bcf212e
Show file tree
Hide file tree
Showing 3 changed files with 98 additions and 23 deletions.
79 changes: 58 additions & 21 deletions core/solver/gmres.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -76,6 +76,7 @@ std::unique_ptr<LinOp> Gmres<ValueType>::transpose() const
share(as<Transposable>(this->get_preconditioner())->transpose()))
.with_criteria(this->get_stop_criterion_factory())
.with_krylov_dim(this->get_krylov_dim())
.with_flexible(this->get_parameters().flexible)
.on(this->get_executor())
->generate(
share(as<Transposable>(this->get_system_matrix())->transpose()));
Expand All @@ -90,6 +91,7 @@ std::unique_ptr<LinOp> Gmres<ValueType>::conj_transpose() const
as<Transposable>(this->get_preconditioner())->conj_transpose()))
.with_criteria(this->get_stop_criterion_factory())
.with_krylov_dim(this->get_krylov_dim())
.with_flexible(this->get_parameters().flexible)
.on(this->get_executor())
->generate(share(
as<Transposable>(this->get_system_matrix())->conj_transpose()));
Expand Down Expand Up @@ -196,17 +198,25 @@ void Gmres<ValueType>::apply_dense_impl(const VectorType* dense_b,

auto exec = this->get_executor();
this->setup_workspace();

const auto is_flexible = this->get_parameters().flexible;
const auto num_rows = this->get_size()[0];
const auto local_num_rows =
::gko::detail::get_local(dense_b)->get_size()[0];
const auto num_rhs = dense_b->get_size()[1];
const auto krylov_dim = this->get_krylov_dim();
GKO_SOLVER_VECTOR(residual, dense_b);
GKO_SOLVER_VECTOR(preconditioned_vector, dense_b);
// TODO: check whether it can be empty when it is flexible GMRES
auto krylov_bases = this->create_workspace_op_with_type_of(
ws::krylov_bases, dense_b, dim<2>{num_rows * (krylov_dim + 1), num_rhs},
dim<2>{local_num_rows * (krylov_dim + 1), num_rhs});
VectorType* preconditioned_krylov_bases = nullptr;
if (is_flexible) {
preconditioned_krylov_bases = this->create_workspace_op_with_type_of(
ws::preconditioned_krylov_bases, dense_b,
dim<2>{num_rows * (krylov_dim + 1), num_rhs},
dim<2>{local_num_rows * (krylov_dim + 1), num_rhs});
}
// rows: rows of Hessenberg matrix, columns: block for each entry
auto hessenberg = this->template create_workspace_op<LocalVector>(
ws::hessenberg, dim<2>{krylov_dim + 1, krylov_dim * num_rhs});
Expand Down Expand Up @@ -341,9 +351,19 @@ void Gmres<ValueType>::apply_dense_impl(const VectorType* dense_b,
span{local_num_rows * (restart_iter + 1),
local_num_rows * (restart_iter + 2)},
span{0, num_rhs});
// preconditioned_vector = get_preconditioner() * this_krylov
std::unique_ptr<VectorType> preconditioned_krylov;
auto preconditioned_krylov_vector = preconditioned_vector;
if (is_flexible) {
preconditioned_krylov = create_submatrix_helper(
preconditioned_krylov_bases, dim<2>{num_rows, num_rhs},
span{local_num_rows * restart_iter,
local_num_rows * (restart_iter + 1)},
span{0, num_rhs});
preconditioned_krylov_vector = preconditioned_krylov.get();
}
// preconditioned_krylov_vector = get_preconditioner() * this_krylov
this->get_preconditioner()->apply(this_krylov.get(),
preconditioned_vector);
preconditioned_krylov_vector);

// Create view of current column in the hessenberg matrix:
// hessenberg_iter = hessenberg(:, restart_iter);
Expand All @@ -352,8 +372,8 @@ void Gmres<ValueType>::apply_dense_impl(const VectorType* dense_b,
span{num_rhs * restart_iter, num_rhs * (restart_iter + 1)});

// Start of Arnoldi
// next_krylov = A * preconditioned_vector
this->get_system_matrix()->apply(preconditioned_vector,
// next_krylov = A * preconditioned_krylov_vector
this->get_system_matrix()->apply(preconditioned_krylov_vector,
next_krylov.get());

for (size_type i = 0; i <= restart_iter; i++) {
Expand Down Expand Up @@ -411,9 +431,6 @@ void Gmres<ValueType>::apply_dense_impl(const VectorType* dense_b,
restart_iter++;
}

auto krylov_bases_small = create_submatrix_helper(
krylov_bases, dim<2>{num_rows, num_rhs},
span{0, local_num_rows * (restart_iter + 1)}, span{0, num_rhs});
auto hessenberg_small = hessenberg->create_submatrix(
span{0, restart_iter}, span{0, num_rhs * (restart_iter)});

Expand All @@ -422,15 +439,30 @@ void Gmres<ValueType>::apply_dense_impl(const VectorType* dense_b,
exec->run(gmres::make_solve_krylov(
residual_norm_collection, hessenberg_small.get(), y,
final_iter_nums.get_const_data(), stop_status.get_const_data()));
// before_preconditioner = krylov_bases * y
exec->run(gmres::make_multi_axpy(
gko::detail::get_local(krylov_bases_small.get()), y,
gko::detail::get_local(before_preconditioner),
final_iter_nums.get_const_data(), stop_status.get_data()));

// x = x + get_preconditioner() * before_preconditioner
this->get_preconditioner()->apply(before_preconditioner,
after_preconditioner);
if (is_flexible) {
auto preconditioned_krylov_bases_small = create_submatrix_helper(
preconditioned_krylov_bases, dim<2>{num_rows, num_rhs},
span{0, local_num_rows * (restart_iter + 1)}, span{0, num_rhs});
// after_preconditioner = preconditioned_krylov_bases * y
exec->run(gmres::make_multi_axpy(
gko::detail::get_local(preconditioned_krylov_bases_small.get()), y,
gko::detail::get_local(after_preconditioner),
final_iter_nums.get_const_data(), stop_status.get_data()));
} else {
auto krylov_bases_small = create_submatrix_helper(
krylov_bases, dim<2>{num_rows, num_rhs},
span{0, local_num_rows * (restart_iter + 1)}, span{0, num_rhs});
// before_preconditioner = krylov_bases * y
exec->run(gmres::make_multi_axpy(
gko::detail::get_local(krylov_bases_small.get()), y,
gko::detail::get_local(before_preconditioner),
final_iter_nums.get_const_data(), stop_status.get_data()));

// after_preconditioner = get_preconditioner() * before_preconditioner
this->get_preconditioner()->apply(before_preconditioner,
after_preconditioner);
}
// x = x + after_preconditioner
dense_x->add_scaled(one_op, after_preconditioner);
}

Expand Down Expand Up @@ -463,7 +495,7 @@ int workspace_traits<Gmres<ValueType>>::num_arrays(const Solver&)
template <typename ValueType>
int workspace_traits<Gmres<ValueType>>::num_vectors(const Solver&)
{
return 14;
return 15;
}


Expand All @@ -484,7 +516,8 @@ std::vector<std::string> workspace_traits<Gmres<ValueType>>::op_names(
"after_preconditioner",
"one",
"minus_one",
"next_krylov_norm_tmp"};
"next_krylov_norm_tmp",
"preconditioned_krylov_bases"};
}


Expand All @@ -509,8 +542,12 @@ std::vector<int> workspace_traits<Gmres<ValueType>>::scalars(const Solver&)
template <typename ValueType>
std::vector<int> workspace_traits<Gmres<ValueType>>::vectors(const Solver&)
{
return {residual, preconditioned_vector, krylov_bases,
before_preconditioner, after_preconditioner};
return {residual,
preconditioned_vector,
krylov_bases,
before_preconditioner,
after_preconditioner,
preconditioned_krylov_bases};
}


Expand Down
7 changes: 7 additions & 0 deletions include/ginkgo/core/solver/gmres.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -132,6 +132,11 @@ class Gmres
* Krylov dimension factory.
*/
size_type GKO_FACTORY_PARAMETER_SCALAR(krylov_dim, 0u);

/**
* Flexible GMRES
*/
bool GKO_FACTORY_PARAMETER_SCALAR(flexible, false);
};
GKO_ENABLE_LIN_OP_FACTORY(Gmres, parameters, Factory);
GKO_ENABLE_BUILD_METHOD(Factory);
Expand Down Expand Up @@ -208,6 +213,8 @@ struct workspace_traits<Gmres<ValueType>> {
constexpr static int minus_one = 12;
// temporary norm vector of next_krylov to copy into hessenberg matrix
constexpr static int next_krylov_norm_tmp = 13;
// preconditioned krylov basis multivector
constexpr static int preconditioned_krylov_bases = 14;

// stopping status array
constexpr static int stop = 0;
Expand Down
35 changes: 33 additions & 2 deletions test/solver/solver.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -274,6 +274,36 @@ struct Gmres : SimpleSolverTest<gko::solver::Gmres<solver_value_type>> {
};


template <unsigned dimension>
struct FGmres : SimpleSolverTest<gko::solver::Gmres<solver_value_type>> {
static typename solver_type::parameters_type build(
std::shared_ptr<const gko::Executor> exec,
gko::size_type iteration_count)
{
return solver_type::build()
.with_criteria(gko::stop::Iteration::build()
.with_max_iters(iteration_count)
.on(exec))
.with_krylov_dim(dimension)
.with_flexible(true);
}

static typename solver_type::parameters_type build_preconditioned(
std::shared_ptr<const gko::Executor> exec,
gko::size_type iteration_count)
{
return solver_type::build()
.with_criteria(gko::stop::Iteration::build()
.with_max_iters(iteration_count)
.on(exec))
.with_preconditioner(
precond_type::build().with_max_block_size(1u).on(exec))
.with_krylov_dim(dimension)
.with_flexible(true);
}
};


struct LowerTrs : SimpleSolverTest<gko::solver::LowerTrs<solver_value_type>> {
static constexpr bool will_not_allocate() { return false; }

Expand Down Expand Up @@ -855,8 +885,9 @@ using SolverTypes =
::testing::Types<Cg, Cgs, Fcg, Bicg, Bicgstab,
/* "IDR uses different initialization approaches even when
deterministic", Idr<1>, Idr<4>,*/
Ir, CbGmres<2>, CbGmres<10>, Gmres<2>, Gmres<10>, LowerTrs,
UpperTrs, LowerTrsUnitdiag, UpperTrsUnitdiag
Ir, CbGmres<2>, CbGmres<10>, Gmres<2>, Gmres<10>,
FGmres<2>, FGmres<10>, LowerTrs, UpperTrs,
LowerTrsUnitdiag, UpperTrsUnitdiag
#ifdef GKO_COMPILING_CUDA
,
LowerTrsSyncfree, UpperTrsSyncfree,
Expand Down

0 comments on commit bcf212e

Please sign in to comment.