1515#include < ginkgo/core/base/exception_helpers.hpp>
1616#include < ginkgo/core/base/lin_op.hpp>
1717#include < ginkgo/core/base/precision_dispatch.hpp>
18+ #include < ginkgo/core/base/utils_helper.hpp>
1819#include < ginkgo/core/config/config.hpp>
1920#include < ginkgo/core/config/registry.hpp>
2021#include < ginkgo/core/factorization/par_ic.hpp>
2122#include < ginkgo/core/matrix/dense.hpp>
22- #include < ginkgo/core/preconditioner/isai.hpp>
23- #include < ginkgo/core/preconditioner/utils.hpp>
24- #include < ginkgo/core/solver/gmres.hpp>
25- #include < ginkgo/core/solver/ir.hpp>
2623#include < ginkgo/core/solver/solver_traits.hpp>
2724#include < ginkgo/core/solver/triangular.hpp>
2825#include < ginkgo/core/stop/combined.hpp>
@@ -35,93 +32,14 @@ namespace preconditioner {
3532namespace detail {
3633
3734
38- // helper for handle the transposed type of concrete type and LinOp
39- template <typename Type>
40- struct transposed_type_impl {
41- using type = typename Type::transposed_type;
42- };
43-
44- template <>
45- struct transposed_type_impl <LinOp> {
46- using type = LinOp;
47- };
48-
49-
50- template <typename Type>
51- using transposed_type = typename transposed_type_impl<Type>::type;
52-
53-
54- // helper to get factory type of concrete type or LinOp
55- template <typename Type>
56- struct factory_type_impl {
57- using type = typename Type::Factory;
58- };
59-
60- template <>
61- struct factory_type_impl <LinOp> {
62- using type = LinOpFactory;
63- };
64-
65-
66- template <typename Type>
67- using factory_type = typename factory_type_impl<Type>::type;
68-
69- template <typename Type>
70- constexpr bool is_ginkgo_linop = std::is_convertible_v<Type*, LinOp*>;
71-
72- template <typename Type>
73- struct get_solver_type_impl {
74- using type = std::conditional_t <is_ginkgo_linop<Type>, Type, LinOp>;
75- };
76-
77- template <typename Type>
78- using get_solver_type = typename get_solver_type_impl<Type>::type;
79-
80-
81- // helper to get value_type of concrete type or void for LinOp
82- template <typename Type, typename = void >
83- struct get_value_type_impl {
84- using type = typename Type::value_type;
85- };
86-
87- // We need to use SFINAE not conditional_t because both type needs to be valid
88- // in conditional_t
89- template <typename Type>
90- struct get_value_type_impl <Type, std::enable_if_t <!is_ginkgo_linop<Type>>> {
91- using type = Type;
92- };
93-
94- template <>
95- struct get_value_type_impl <LinOp> {
96- using type = void ;
97- };
98-
99-
100- template <typename Type>
101- using get_value_type = typename get_value_type_impl<Type>::type;
102-
103-
104- // get_first_template is to get the first template argument of class.
105- // It can be easily done by introducing another member type of IC to alias the
106- // first template argument, but it introduces another public interface.
107- template <class >
108- struct get_first_template {};
109-
110- template <template <typename ...> class Base , class First , class ... Rest>
111- struct get_first_template <Base<First, Rest...>> {
112- using type = First;
113- };
114-
115-
116- // reuse the above the condition
11735template <typename SolverTypeOrValueType>
11836constexpr bool support_ic_parse =
119- std::is_same_v<get_solver_type<SolverTypeOrValueType>, LinOp> &&
120- !std::is_same_v<get_value_type<SolverTypeOrValueType>, void >;
37+ std::is_same_v<gko::detail::get_solver_type<SolverTypeOrValueType>, LinOp>;
12138
12239
123- template <typename Ic, std::enable_if_t <!support_ic_parse<
124- typename get_first_template<Ic>::type>>* = nullptr >
40+ template <typename Ic,
41+ std::enable_if_t <!support_ic_parse<
42+ typename gko::detail::get_first_template<Ic>::type>>* = nullptr >
12543typename Ic::parameters_type ic_parse (
12644 const config::pnode& config, const config::registry& context,
12745 const config::type_descriptor& td_for_child)
@@ -130,15 +48,17 @@ typename Ic::parameters_type ic_parse(
13048 " preconditioner::Ic only supports limited type for parse." );
13149}
13250
133- template <typename Ic, std::enable_if_t <support_ic_parse<
134- typename get_first_template<Ic>::type>>* = nullptr >
51+ template <typename Ic,
52+ std::enable_if_t <support_ic_parse<
53+ typename gko::detail::get_first_template<Ic>::type>>* = nullptr >
13554typename Ic::parameters_type ic_parse (
13655 const config::pnode& config, const config::registry& context,
13756 const config::type_descriptor& td_for_child);
13857
13958
14059} // namespace detail
14160
61+
14262/* *
14363 * The Incomplete Cholesky (IC) preconditioner solves the equation $LL^H*x = b$
14464 * for a given lower triangular matrix L and the right hand side b (can contain
@@ -196,14 +116,13 @@ class Ic : public EnableLinOp<Ic<LSolverTypeOrValueType, IndexType>>,
196116 friend class EnablePolymorphicObject <Ic, LinOp>;
197117
198118public:
199- using l_solver_type = detail::get_solver_type<LSolverTypeOrValueType>;
200- static_assert (
201- std::is_same<
202- detail::transposed_type<detail::transposed_type<l_solver_type>>,
203- l_solver_type>::value,
204- " l_solver_type::transposed_type must be symmetric" );
205- using value_type = detail::get_value_type<LSolverTypeOrValueType>;
206- using lh_solver_type = detail::transposed_type<l_solver_type>;
119+ using l_solver_type = gko::detail::get_solver_type<LSolverTypeOrValueType>;
120+ static_assert (std::is_same<gko::detail::transposed_type<
121+ gko::detail::transposed_type<l_solver_type>>,
122+ l_solver_type>::value,
123+ " l_solver_type::transposed_type must be symmetric" );
124+ using value_type = gko::detail::get_value_type<LSolverTypeOrValueType>;
125+ using lh_solver_type = gko::detail::transposed_type<l_solver_type>;
207126 using index_type = IndexType;
208127 using transposed_type = Ic<LSolverTypeOrValueType, IndexType>;
209128
@@ -214,7 +133,7 @@ class Ic : public EnableLinOp<Ic<LSolverTypeOrValueType, IndexType>>,
214133 /* *
215134 * Factory for the L solver
216135 */
217- std::shared_ptr<const detail::factory_type<l_solver_type>>
136+ std::shared_ptr<const gko:: detail::factory_type<l_solver_type>>
218137 l_solver_factory{};
219138
220139 /* *
@@ -225,15 +144,15 @@ class Ic : public EnableLinOp<Ic<LSolverTypeOrValueType, IndexType>>,
225144 GKO_DEPRECATED (" use with_l_solver instead" )
226145 parameters_type& with_l_solver_factory (
227146 deferred_factory_parameter<
228- const detail::factory_type<l_solver_type>>
147+ const gko:: detail::factory_type<l_solver_type>>
229148 solver)
230149 {
231150 return with_l_solver (std::move (solver));
232151 }
233152
234153 parameters_type& with_l_solver (
235154 deferred_factory_parameter<
236- const detail::factory_type<l_solver_type>>
155+ const gko:: detail::factory_type<l_solver_type>>
237156 solver)
238157 {
239158 this ->l_solver_generator = std::move (solver);
@@ -269,7 +188,8 @@ class Ic : public EnableLinOp<Ic<LSolverTypeOrValueType, IndexType>>,
269188 }
270189
271190 private:
272- deferred_factory_parameter<const detail::factory_type<l_solver_type>>
191+ deferred_factory_parameter<
192+ const gko::detail::factory_type<l_solver_type>>
273193 l_solver_generator;
274194
275195 deferred_factory_parameter<const LinOpFactory> factorization_generator;
@@ -331,10 +251,10 @@ class Ic : public EnableLinOp<Ic<LSolverTypeOrValueType, IndexType>>,
331251 new transposed_type{this ->get_executor ()}};
332252 transposed->set_size (gko::transpose (this ->get_size ()));
333253 transposed->l_solver_ =
334- share (as<detail::transposed_type<lh_solver_type>>(
254+ share (as<gko:: detail::transposed_type<lh_solver_type>>(
335255 as<Transposable>(this ->get_lh_solver ())->transpose ()));
336256 transposed->lh_solver_ =
337- share (as<detail::transposed_type<l_solver_type>>(
257+ share (as<gko:: detail::transposed_type<l_solver_type>>(
338258 as<Transposable>(this ->get_l_solver ())->transpose ()));
339259
340260 return std::move (transposed);
@@ -346,10 +266,10 @@ class Ic : public EnableLinOp<Ic<LSolverTypeOrValueType, IndexType>>,
346266 new transposed_type{this ->get_executor ()}};
347267 transposed->set_size (gko::transpose (this ->get_size ()));
348268 transposed->l_solver_ =
349- share (as<detail::transposed_type<lh_solver_type>>(
269+ share (as<gko:: detail::transposed_type<lh_solver_type>>(
350270 as<Transposable>(this ->get_lh_solver ())->conj_transpose ()));
351271 transposed->lh_solver_ =
352- share (as<detail::transposed_type<l_solver_type>>(
272+ share (as<gko:: detail::transposed_type<l_solver_type>>(
353273 as<Transposable>(this ->get_l_solver ())->conj_transpose ()));
354274
355275 return std::move (transposed);
@@ -451,24 +371,17 @@ class Ic : public EnableLinOp<Ic<LSolverTypeOrValueType, IndexType>>,
451371 // build factorization if we weren't passed a composition
452372 if (!comp) {
453373 auto exec = lin_op->get_executor ();
454- if constexpr (std::is_same_v<value_type, void >) {
455- GKO_NOT_IMPLEMENTED;
456- } else {
457- if (!parameters_.factorization_factory ) {
458- parameters_.factorization_factory =
459- factorization::ParIc<value_type, index_type>::build ()
460- .with_both_factors (false )
461- .on (exec);
462- }
463- auto fact = std::shared_ptr<const LinOp>(
464- parameters_.factorization_factory ->generate (lin_op));
465- // ensure that the result is a composition
466- comp = std::dynamic_pointer_cast<const Composition<value_type>>(
467- fact);
468- if (!comp) {
469- GKO_NOT_SUPPORTED (comp);
470- }
374+
375+ if (!parameters_.factorization_factory ) {
376+ parameters_.factorization_factory =
377+ factorization::ParIc<value_type, index_type>::build ()
378+ .with_both_factors (false )
379+ .on (exec);
471380 }
381+ auto fact = std::shared_ptr<const LinOp>(
382+ parameters_.factorization_factory ->generate (lin_op));
383+ // ensure that the result is a composition
384+ comp = gko::as<const Composition<value_type>>(fact);
472385 }
473386 // comp must contain one or two factors
474387 if (comp->get_operators ().size () > 2 || comp->get_operators ().empty ()) {
@@ -519,8 +432,6 @@ class Ic : public EnableLinOp<Ic<LSolverTypeOrValueType, IndexType>>,
519432 cache_.intermediate =
520433 matrix::Dense<value_type>::create (this ->get_executor ());
521434 }
522- // Use b as the initial guess for the first triangular solve
523- cache_.intermediate ->copy_from (b);
524435 }
525436
526437 /* *
0 commit comments