Ginkgo  Generated from pipelines/2017069469 branch based on develop. Ginkgo version 1.11.0
A numerical linear algebra library targeting many-core architectures
ilu.hpp
1 // SPDX-FileCopyrightText: 2017 - 2025 The Ginkgo authors
2 //
3 // SPDX-License-Identifier: BSD-3-Clause
4 
5 #ifndef GKO_PUBLIC_CORE_PRECONDITIONER_ILU_HPP_
6 #define GKO_PUBLIC_CORE_PRECONDITIONER_ILU_HPP_
7 
8 
9 #include <memory>
10 #include <type_traits>
11 
12 #include <ginkgo/core/base/abstract_factory.hpp>
13 #include <ginkgo/core/base/composition.hpp>
14 #include <ginkgo/core/base/exception.hpp>
15 #include <ginkgo/core/base/exception_helpers.hpp>
16 #include <ginkgo/core/base/lin_op.hpp>
17 #include <ginkgo/core/base/precision_dispatch.hpp>
18 #include <ginkgo/core/base/type_traits.hpp>
19 #include <ginkgo/core/config/config.hpp>
20 #include <ginkgo/core/config/registry.hpp>
21 #include <ginkgo/core/factorization/par_ilu.hpp>
22 #include <ginkgo/core/matrix/dense.hpp>
23 #include <ginkgo/core/solver/solver_traits.hpp>
24 #include <ginkgo/core/solver/triangular.hpp>
25 #include <ginkgo/core/stop/combined.hpp>
26 #include <ginkgo/core/stop/iteration.hpp>
27 #include <ginkgo/core/stop/residual_norm.hpp>
28 
29 
30 namespace gko {
31 namespace preconditioner {
32 namespace detail {
33 
34 
35 template <typename Type>
36 constexpr bool support_ilu_parse =
37  std::is_same_v<typename Type::l_solver_type, LinOp>&&
38  std::is_same_v<typename Type::u_solver_type, LinOp>;
39 
40 
41 template <typename Ilu, std::enable_if_t<!support_ilu_parse<Ilu>>* = nullptr>
42 typename Ilu::parameters_type ilu_parse(
43  const config::pnode& config, const config::registry& context,
44  const config::type_descriptor& td_for_child)
45 {
46  GKO_INVALID_STATE(
47  "preconditioner::Ilu only supports limited type for parse.");
48 }
49 
50 template <typename Ilu, std::enable_if_t<support_ilu_parse<Ilu>>* = nullptr>
51 typename Ilu::parameters_type ilu_parse(
52  const config::pnode& config, const config::registry& context,
53  const config::type_descriptor& td_for_child);
54 
55 } // namespace detail
56 
57 
110 template <typename LSolverTypeOrValueType = solver::LowerTrs<>,
111  typename USolverTypeOrValueType =
112  gko::detail::transposed_type<LSolverTypeOrValueType>,
113  bool ReverseApply = false, typename IndexType = int32>
114 class Ilu
115  : public EnableLinOp<Ilu<LSolverTypeOrValueType, USolverTypeOrValueType,
116  ReverseApply, IndexType>>,
117  public Transposable {
118  friend class EnableLinOp<Ilu>;
119  friend class EnablePolymorphicObject<Ilu, LinOp>;
120 
121 public:
122  static_assert(
123  std::is_same_v<gko::detail::get_value_type<LSolverTypeOrValueType>,
124  gko::detail::get_value_type<USolverTypeOrValueType>>,
125  "Both the L- and the U-solver must use the same `value_type`!");
126  using value_type = gko::detail::get_value_type<LSolverTypeOrValueType>;
127  using l_solver_type =
128  std::conditional_t<gko::detail::is_ginkgo_linop<LSolverTypeOrValueType>,
129  LSolverTypeOrValueType, LinOp>;
130  using u_solver_type =
131  std::conditional_t<gko::detail::is_ginkgo_linop<USolverTypeOrValueType>,
132  USolverTypeOrValueType, LinOp>;
133  static constexpr bool performs_reverse_apply = ReverseApply;
134  using index_type = IndexType;
135  using transposed_type =
137  gko::detail::transposed_type<LSolverTypeOrValueType>, ReverseApply,
138  IndexType>;
139 
140  class Factory;
141 
143  : public enable_parameters_type<parameters_type, Factory> {
147  std::shared_ptr<const gko::detail::factory_type<l_solver_type>>
149 
153  std::shared_ptr<const gko::detail::factory_type<u_solver_type>>
155 
159  std::shared_ptr<const LinOpFactory> factorization_factory{};
160 
161  GKO_DEPRECATED("use with_l_solver instead")
162  parameters_type& with_l_solver_factory(
164  const gko::detail::factory_type<l_solver_type>>
165  solver)
166  {
167  return with_l_solver(std::move(solver));
168  }
169 
177  const gko::detail::factory_type<l_solver_type>>
178  solver)
179  {
180  this->l_solver_generator = std::move(solver);
181  this->deferred_factories["l_solver"] = [](const auto& exec,
182  auto& params) {
183  if (!params.l_solver_generator.is_empty()) {
184  params.l_solver_factory =
185  params.l_solver_generator.on(exec);
186  }
187  };
188  return *this;
189  }
190 
191  GKO_DEPRECATED("use with_u_solver instead")
192  parameters_type& with_u_solver_factory(
194  const gko::detail::factory_type<u_solver_type>>
195  solver)
196  {
197  return with_u_solver(std::move(solver));
198  }
199 
207  const gko::detail::factory_type<u_solver_type>>
208  solver)
209  {
210  this->u_solver_generator = std::move(solver);
211  this->deferred_factories["u_solver"] = [](const auto& exec,
212  auto& params) {
213  if (!params.u_solver_generator.is_empty()) {
214  params.u_solver_factory =
215  params.u_solver_generator.on(exec);
216  }
217  };
218  return *this;
219  }
220 
221  GKO_DEPRECATED("use with_factorization instead")
222  parameters_type& with_factorization_factory(
223  deferred_factory_parameter<const LinOpFactory> factorization)
224  {
225  return with_factorization(std::move(factorization));
226  }
227 
228  parameters_type& with_factorization(
230  {
231  this->factorization_generator = std::move(factorization);
232  this->deferred_factories["factorization"] = [](const auto& exec,
233  auto& params) {
234  if (!params.factorization_generator.is_empty()) {
235  params.factorization_factory =
236  params.factorization_generator.on(exec);
237  }
238  };
239  return *this;
240  }
241 
242  private:
243  deferred_factory_parameter<
244  const gko::detail::factory_type<l_solver_type>>
245  l_solver_generator;
246 
247  deferred_factory_parameter<
248  const gko::detail::factory_type<u_solver_type>>
249  u_solver_generator;
250 
251  deferred_factory_parameter<const LinOpFactory> factorization_generator;
252  };
253 
256 
275  const config::pnode& config, const config::registry& context,
276  const config::type_descriptor& td_for_child =
277  config::make_type_descriptor<value_type, index_type>())
278  {
279  // parse is not templated, so we can only use SFINAE later
280  return detail::ilu_parse<Ilu>(config, context, td_for_child);
281  }
282 
288  std::shared_ptr<const l_solver_type> get_l_solver() const
289  {
290  return l_solver_;
291  }
292 
298  std::shared_ptr<const u_solver_type> get_u_solver() const
299  {
300  return u_solver_;
301  }
302 
303  std::unique_ptr<LinOp> transpose() const override
304  {
305  std::unique_ptr<transposed_type> transposed{
306  new transposed_type{this->get_executor()}};
307  transposed->set_size(gko::transpose(this->get_size()));
308  transposed->l_solver_ =
309  share(as<gko::detail::transposed_type<u_solver_type>>(
310  as<Transposable>(this->get_u_solver())->transpose()));
311  transposed->u_solver_ =
312  share(as<gko::detail::transposed_type<l_solver_type>>(
313  as<Transposable>(this->get_l_solver())->transpose()));
314 
315  return std::move(transposed);
316  }
317 
318  std::unique_ptr<LinOp> conj_transpose() const override
319  {
320  std::unique_ptr<transposed_type> transposed{
321  new transposed_type{this->get_executor()}};
322  transposed->set_size(gko::transpose(this->get_size()));
323  transposed->l_solver_ =
324  share(as<gko::detail::transposed_type<u_solver_type>>(
325  as<Transposable>(this->get_u_solver())->conj_transpose()));
326  transposed->u_solver_ =
327  share(as<gko::detail::transposed_type<l_solver_type>>(
328  as<Transposable>(this->get_l_solver())->conj_transpose()));
329 
330  return std::move(transposed);
331  }
332 
338  Ilu& operator=(const Ilu& other)
339  {
340  if (&other != this) {
342  auto exec = this->get_executor();
343  l_solver_ = other.l_solver_;
344  u_solver_ = other.u_solver_;
345  parameters_ = other.parameters_;
346  if (other.get_executor() != exec) {
347  l_solver_ = gko::clone(exec, l_solver_);
348  u_solver_ = gko::clone(exec, u_solver_);
349  }
350  }
351  return *this;
352  }
353 
360  Ilu& operator=(Ilu&& other)
361  {
362  if (&other != this) {
364  auto exec = this->get_executor();
365  l_solver_ = std::move(other.l_solver_);
366  u_solver_ = std::move(other.u_solver_);
367  parameters_ = std::exchange(other.parameters_, parameters_type{});
368  if (other.get_executor() != exec) {
369  l_solver_ = gko::clone(exec, l_solver_);
370  u_solver_ = gko::clone(exec, u_solver_);
371  }
372  }
373  return *this;
374  }
375 
380  Ilu(const Ilu& other) : Ilu{other.get_executor()} { *this = other; }
381 
387  Ilu(Ilu&& other) : Ilu{other.get_executor()} { *this = std::move(other); }
388 
389 protected:
390  void apply_impl(const LinOp* b, LinOp* x) const override
391  {
392  // take care of real-to-complex apply
393  precision_dispatch_real_complex<value_type>(
394  [&](auto dense_b, auto dense_x) {
395  this->set_cache_to(dense_b);
396  if (!ReverseApply) {
397  l_solver_->apply(dense_b, cache_.intermediate);
398  if (u_solver_->apply_uses_initial_guess()) {
399  dense_x->copy_from(cache_.intermediate);
400  }
401  u_solver_->apply(cache_.intermediate, dense_x);
402  } else {
403  u_solver_->apply(dense_b, cache_.intermediate);
404  if (l_solver_->apply_uses_initial_guess()) {
405  dense_x->copy_from(cache_.intermediate);
406  }
407  l_solver_->apply(cache_.intermediate, dense_x);
408  }
409  },
410  b, x);
411  }
412 
413  void apply_impl(const LinOp* alpha, const LinOp* b, const LinOp* beta,
414  LinOp* x) const override
415  {
416  precision_dispatch_real_complex<value_type>(
417  [&](auto dense_alpha, auto dense_b, auto dense_beta, auto dense_x) {
418  this->set_cache_to(dense_b);
419  if (!ReverseApply) {
420  l_solver_->apply(dense_b, cache_.intermediate);
421  u_solver_->apply(dense_alpha, cache_.intermediate,
422  dense_beta, dense_x);
423  } else {
424  u_solver_->apply(dense_b, cache_.intermediate);
425  l_solver_->apply(dense_alpha, cache_.intermediate,
426  dense_beta, dense_x);
427  }
428  },
429  alpha, b, beta, x);
430  }
431 
432  explicit Ilu(std::shared_ptr<const Executor> exec)
433  : EnableLinOp<Ilu>(std::move(exec))
434  {}
435 
436  explicit Ilu(const Factory* factory, std::shared_ptr<const LinOp> lin_op)
437  : EnableLinOp<Ilu>(factory->get_executor(), lin_op->get_size()),
438  parameters_{factory->get_parameters()}
439  {
440  auto comp =
441  std::dynamic_pointer_cast<const Composition<value_type>>(lin_op);
442  std::shared_ptr<const LinOp> l_factor;
443  std::shared_ptr<const LinOp> u_factor;
444 
445  // build factorization if we weren't passed a composition
446  if (!comp) {
447  auto exec = lin_op->get_executor();
448  if (!parameters_.factorization_factory) {
449  parameters_.factorization_factory =
450  factorization::ParIlu<value_type, index_type>::build().on(
451  exec);
452  }
453  auto fact = std::shared_ptr<const LinOp>(
454  parameters_.factorization_factory->generate(lin_op));
455  // ensure that the result is a composition
456  comp = as<const Composition<value_type>>(fact);
457  }
458  if (comp->get_operators().size() == 2) {
459  l_factor = comp->get_operators()[0];
460  u_factor = comp->get_operators()[1];
461  } else {
462  GKO_NOT_SUPPORTED(comp);
463  }
464  GKO_ASSERT_EQUAL_DIMENSIONS(l_factor, u_factor);
465 
466  auto exec = this->get_executor();
467 
468  // If no factories are provided, generate default ones
469  if (!parameters_.l_solver_factory) {
470  // when l_solver_type is LinOp, use LowerTrs as the default one
471  l_solver_ = generate_default_solver<std::conditional_t<
472  std::is_same_v<l_solver_type, LinOp>,
473  solver::LowerTrs<value_type, index_type>, l_solver_type>>(
474  exec, l_factor);
475  } else {
476  l_solver_ = parameters_.l_solver_factory->generate(l_factor);
477  }
478  if (!parameters_.u_solver_factory) {
479  // when u_solver_type is LinOp, use UpperTrs as the default one
480  u_solver_ = generate_default_solver<std::conditional_t<
481  std::is_same_v<u_solver_type, LinOp>,
482  solver::UpperTrs<value_type, index_type>, u_solver_type>>(
483  exec, u_factor);
484  } else {
485  u_solver_ = parameters_.u_solver_factory->generate(u_factor);
486  }
487  }
488 
496  void set_cache_to(const LinOp* b) const
497  {
498  if (cache_.intermediate == nullptr) {
499  cache_.intermediate =
501  }
502  // Use b as the initial guess for the first triangular solve
503  cache_.intermediate->copy_from(b);
504  }
505 
506 
514  template <typename SolverType>
515  static std::enable_if_t<solver::has_with_criteria<SolverType>::value,
516  std::unique_ptr<SolverType>>
517  generate_default_solver(const std::shared_ptr<const Executor>& exec,
518  const std::shared_ptr<const LinOp>& mtx)
519  {
520  // half can not use constexpr constructor
521  const gko::remove_complex<value_type> default_reduce_residual{1e-4};
522  const unsigned int default_max_iters{
523  static_cast<unsigned int>(mtx->get_size()[0])};
524 
525  return SolverType::build()
526  .with_criteria(
527  gko::stop::Iteration::build().with_max_iters(default_max_iters),
529  .with_reduction_factor(default_reduce_residual))
530  .on(exec)
531  ->generate(mtx);
532  }
533 
537  template <typename SolverType>
538  static std::enable_if_t<!solver::has_with_criteria<SolverType>::value,
539  std::unique_ptr<SolverType>>
540  generate_default_solver(const std::shared_ptr<const Executor>& exec,
541  const std::shared_ptr<const LinOp>& mtx)
542  {
543  return SolverType::build().on(exec)->generate(mtx);
544  }
545 
546 private:
547  std::shared_ptr<const l_solver_type> l_solver_{};
548  std::shared_ptr<const u_solver_type> u_solver_{};
559  mutable struct cache_struct {
560  cache_struct() = default;
561  ~cache_struct() = default;
562  cache_struct(const cache_struct&) {}
563  cache_struct(cache_struct&&) {}
564  cache_struct& operator=(const cache_struct&) { return *this; }
565  cache_struct& operator=(cache_struct&&) { return *this; }
566  std::unique_ptr<LinOp> intermediate{};
567  } cache_;
568 };
569 
570 
571 } // namespace preconditioner
572 } // namespace gko
573 
574 
575 #endif // GKO_PUBLIC_CORE_PRECONDITIONER_ILU_HPP_
gko::preconditioner::Ilu::Factory
Definition: ilu.hpp:254
gko::config::pnode
pnode describes a tree of properties.
Definition: property_tree.hpp:28
gko::LinOp
Definition: lin_op.hpp:117
gko::preconditioner::Ilu::operator=
Ilu & operator=(Ilu &&other)
Move-assigns an ILU preconditioner.
Definition: ilu.hpp:360
gko::preconditioner::Ilu::transpose
std::unique_ptr< LinOp > transpose() const override
Returns a LinOp representing the transpose of the Transposable object.
Definition: ilu.hpp:303
gko::matrix::Dense::create
static std::unique_ptr< Dense > create(std::shared_ptr< const Executor > exec, const dim< 2 > &size={}, size_type stride=0)
Creates an uninitialized Dense matrix of the specified size.
gko::log::profile_event_category::factory
LinOpFactory events.
gko::Transposable
Linear operators which support transposition should implement the Transposable interface.
Definition: lin_op.hpp:433
gko::preconditioner::Ilu::operator=
Ilu & operator=(const Ilu &other)
Copy-assigns an ILU preconditioner.
Definition: ilu.hpp:338
gko::preconditioner::Ilu::conj_transpose
std::unique_ptr< LinOp > conj_transpose() const override
Returns a LinOp representing the conjugate transpose of the Transposable object.
Definition: ilu.hpp:318
gko::preconditioner::Ilu::parameters_type::u_solver_factory
std::shared_ptr< const gko::detail::factory_type< u_solver_type > > u_solver_factory
Factory for the U solver.
Definition: ilu.hpp:154
gko::config::type_descriptor
This class describes the value and index types to be used when building a Ginkgo type from a configur...
Definition: type_descriptor.hpp:39
gko::clone
detail::cloned_type< Pointer > clone(const Pointer &p)
Creates a unique clone of the object pointed to by p.
Definition: utils_helper.hpp:173
gko::preconditioner::Ilu::Ilu
Ilu(const Ilu &other)
Copy-constructs an ILU preconditioner.
Definition: ilu.hpp:380
gko
The Ginkgo namespace.
Definition: abstract_factory.hpp:20
gko::stop::ResidualNorm
The ResidualNorm class is a stopping criterion which stops the iteration process when the actual resi...
Definition: residual_norm.hpp:113
gko::as
std::decay_t< T > * as(U *obj)
Performs polymorphic type conversion.
Definition: utils_helper.hpp:307
gko::preconditioner::Ilu::parameters_type::with_u_solver
parameters_type & with_u_solver(deferred_factory_parameter< const gko::detail::factory_type< u_solver_type >> solver)
When USolverTypeOrValueType is a concrete solver type, this only accepts the factory from the same co...
Definition: ilu.hpp:205
gko::preconditioner::Ilu::parameters_type::with_l_solver
parameters_type & with_l_solver(deferred_factory_parameter< const gko::detail::factory_type< l_solver_type >> solver)
When LSolverTypeOrValueType is a concrete solver type, this only accepts the factory from the same co...
Definition: ilu.hpp:175
gko::preconditioner::Ilu::get_l_solver
std::shared_ptr< const l_solver_type > get_l_solver() const
Returns the solver which is used for the provided L matrix.
Definition: ilu.hpp:288
GKO_ENABLE_LIN_OP_FACTORY
#define GKO_ENABLE_LIN_OP_FACTORY(_lin_op, _parameters_name, _factory_name)
This macro will generate a default implementation of a LinOpFactory for the LinOp subclass it is defi...
Definition: lin_op.hpp:1017
gko::share
detail::shared_type< OwningPointer > share(OwningPointer &&p)
Marks the object pointed to by p as shared.
Definition: utils_helper.hpp:224
gko::transpose
batch_dim< 2, DimensionType > transpose(const batch_dim< 2, DimensionType > &input)
Returns a batch_dim object with its dimensions swapped for batched operators.
Definition: batch_dim.hpp:119
gko::preconditioner::Ilu::parameters_type
Definition: ilu.hpp:142
gko::preconditioner::Ilu::parameters_type::l_solver_factory
std::shared_ptr< const gko::detail::factory_type< l_solver_type > > l_solver_factory
Factory for the L solver.
Definition: ilu.hpp:148
gko::config::registry
This class stores additional context for creating Ginkgo objects from configuration files.
Definition: registry.hpp:167
GKO_ENABLE_BUILD_METHOD
#define GKO_ENABLE_BUILD_METHOD(_factory_name)
Defines a build method for the factory, simplifying its construction by removing the repetitive typin...
Definition: abstract_factory.hpp:394
gko::preconditioner::Ilu::parse
static parameters_type parse(const config::pnode &config, const config::registry &context, const config::type_descriptor &td_for_child=config::make_type_descriptor< value_type, index_type >())
Create the parameters from the property_tree.
Definition: ilu.hpp:274
gko::LinOpFactory
A LinOpFactory represents a higher order mapping which transforms one linear operator into another.
Definition: lin_op.hpp:384
gko::PolymorphicObject::get_executor
std::shared_ptr< const Executor > get_executor() const noexcept
Returns the Executor of the object.
Definition: polymorphic_object.hpp:243
gko::preconditioner::Ilu::get_u_solver
std::shared_ptr< const u_solver_type > get_u_solver() const
Returns the solver which is used for the provided U matrix.
Definition: ilu.hpp:298
gko::enable_parameters_type
The enable_parameters_type mixin is used to create a base implementation of the factory parameters st...
Definition: abstract_factory.hpp:211
gko::LinOp::get_size
const dim< 2 > & get_size() const noexcept
Returns the size of the operator.
Definition: lin_op.hpp:210
gko::remove_complex
typename detail::remove_complex_s< T >::type remove_complex
Obtain the type which removed the complex of complex/scalar type or the template parameter of class b...
Definition: math.hpp:264
gko::preconditioner::Ilu::parameters_type::factorization_factory
std::shared_ptr< const LinOpFactory > factorization_factory
Factory for the factorization.
Definition: ilu.hpp:159
gko::deferred_factory_parameter
Represents a factory parameter of factory type that can either initialized by a pre-existing factory ...
Definition: abstract_factory.hpp:309
gko::EnableLinOp
The EnableLinOp mixin can be used to provide sensible default implementations of the majority of the ...
Definition: lin_op.hpp:877
gko::LinOp::LinOp
LinOp(const LinOp &)=default
Copy-constructs a LinOp.
gko::preconditioner::Ilu::Ilu
Ilu(Ilu &&other)
Move-constructs an ILU preconditioner.
Definition: ilu.hpp:387
gko::preconditioner::Ilu
The Incomplete LU (ILU) preconditioner solves the equation for a given lower triangular matrix L,...
Definition: ilu.hpp:114
gko::EnablePolymorphicObject
This mixin inherits from (a subclass of) PolymorphicObject and provides a base implementation of a ne...
Definition: polymorphic_object.hpp:667