Reference documentation for deal.II version GIT relicensing-216-gcda843ca70 2024-03-28 12:20:02+00:00
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
Loading...
Searching...
No Matches
ginkgo_solver.cc
Go to the documentation of this file.
1// ------------------------------------------------------------------------
2//
3// SPDX-License-Identifier: LGPL-2.1-or-later
4// Copyright (C) 2019 - 2023 by the deal.II authors
5//
6// This file is part of the deal.II library.
7//
8// Part of the source code is dual licensed under Apache-2.0 WITH
9// LLVM-exception OR LGPL-2.1-or-later. Detailed license information
10// governing the source code and code contributions can be found in
11// LICENSE.md and CONTRIBUTING.md at the top level directory of deal.II.
12//
13// ------------------------------------------------------------------------
14
15
17
19
20#ifdef DEAL_II_WITH_GINKGO
21
23
24# include <cmath>
25# include <memory>
26
27
29
30namespace GinkgoWrappers
31{
32 template <typename ValueType, typename IndexType>
34 const std::string &exec_type)
35 : solver_control(solver_control)
36 , exec_type(exec_type)
37 {
38 if (exec_type == "reference")
39 {
40 executor = gko::ReferenceExecutor::create();
41 }
42 else if (exec_type == "omp")
43 {
44 executor = gko::OmpExecutor::create();
45 }
46 else if (exec_type == "cuda" && gko::CudaExecutor::get_num_devices() > 0)
47 {
48 executor = gko::CudaExecutor::create(0, gko::OmpExecutor::create());
49 }
50 else
51 {
52 Assert(
53 false,
55 " exec_type needs to be one of the three strings: \"reference\", \"cuda\" or \"omp\" "));
56 }
57 using ResidualCriterionFactory = gko::stop::ResidualNormReduction<>;
58 residual_criterion = ResidualCriterionFactory::build()
59 .with_reduction_factor(solver_control.tolerance())
60 .on(executor);
61
63 gko::stop::Combined::build()
64 .with_criteria(residual_criterion,
65 gko::stop::Iteration::build()
66 .with_max_iters(solver_control.max_steps())
67 .on(executor))
68 .on(executor);
69 }
70
71
72
73 template <typename ValueType, typename IndexType>
74 void
76 {
77 // Add the logger object. See the different masks available in Ginkgo's
78 // documentation
79 convergence_logger = gko::log::Convergence<>::create(
80 executor, gko::log::Logger::criterion_check_completed_mask);
81 }
82
83
84
85 template <typename ValueType, typename IndexType>
86 void
88 const Vector<ValueType> &rhs)
89 {
90 // some shortcuts.
91 using val_array = gko::Array<ValueType>;
92 using vec = gko::matrix::Dense<ValueType>;
93
94 Assert(system_matrix, ExcNotInitialized());
95 Assert(executor, ExcNotInitialized());
96 Assert(rhs.size() == solution.size(),
97 ExcDimensionMismatch(rhs.size(), solution.size()));
98
99 // Generate the solver from the solver using the system matrix.
100 auto solver = solver_gen->generate(system_matrix);
101
102 // Create the rhs vector in Ginkgo's format.
103 std::vector<ValueType> f(rhs.size());
104 std::copy(rhs.begin(), rhs.begin() + rhs.size(), f.begin());
105 auto b =
106 vec::create(executor,
107 gko::dim<2>(rhs.size(), 1),
108 val_array::view(executor->get_master(), rhs.size(), f.data()),
109 1);
110
111 // Create the solution vector in Ginkgo's format.
112 std::vector<ValueType> u(solution.size());
113 std::copy(solution.begin(), solution.begin() + solution.size(), u.begin());
114 auto x = vec::create(executor,
115 gko::dim<2>(solution.size(), 1),
116 val_array::view(executor->get_master(),
117 solution.size(),
118 u.data()),
119 1);
120
121 // Create the logger object to log some data from the solvers to confirm
122 // convergence.
123 initialize_ginkgo_log();
124
125 Assert(convergence_logger, ExcNotInitialized());
126 // Add the convergence logger object to the combined factory to retrieve the
127 // solver and other data
128 combined_factory->add_logger(convergence_logger);
129
130 // Finally, apply the solver to b and get the solution in x.
131 solver->apply(gko::lend(b), gko::lend(x));
132
133 // The convergence_logger object contains the residual vector after the
134 // solver has returned. use this vector to compute the residual norm of the
135 // solution. Get the residual norm from the logger. As the convergence
136 // logger returns a `linop`, it is necessary to convert it to a Dense
137 // matrix. Additionally, if the logger is logging on the gpu, it is
138 // necessary to copy the data to the host and hence the
139 // `residual_norm_d_parent`
140 auto residual_norm = convergence_logger->get_residual_norm();
141 auto residual_norm_d =
142 gko::as<gko::matrix::Dense<ValueType>>(residual_norm);
143 auto residual_norm_d_parent =
144 gko::matrix::Dense<ValueType>::create(executor->get_master(),
145 gko::dim<2>{1, 1});
146 residual_norm_d_parent->copy_from(residual_norm_d);
147
148 // Get the number of iterations taken to converge to the solution.
149 auto num_iteration = convergence_logger->get_num_iterations();
150
151 // Ginkgo works with a relative residual norm through its
152 // ResidualNormReduction criterion. Therefore, to get the normalized
153 // residual, we divide by the norm of the rhs.
154 auto b_norm = gko::matrix::Dense<ValueType>::create(executor->get_master(),
155 gko::dim<2>{1, 1});
156 if (executor != executor->get_master())
157 {
158 auto b_master = vec::create(executor->get_master(),
159 gko::dim<2>(rhs.size(), 1),
160 val_array::view(executor->get_master(),
161 rhs.size(),
162 f.data()),
163 1);
164 b_master->compute_norm2(b_norm.get());
165 }
166 else
167 {
168 b->compute_norm2(b_norm.get());
169 }
170
171 Assert(b_norm.get()->at(0, 0) != 0.0, ExcDivideByZero());
172 // Pass the number of iterations and residual norm to the solver_control
173 // object. As both `residual_norm_d_parent` and `b_norm` are seen as Dense
174 // matrices, we use the `at` function to get the first value here. In case
175 // of multiple right hand sides, this will need to be modified.
176 const SolverControl::State state =
177 solver_control.check(num_iteration,
178 residual_norm_d_parent->at(0, 0) / b_norm->at(0, 0));
179
180 // in case of failure: throw exception
181 if (state != SolverControl::success)
182 AssertThrow(false,
183 SolverControl::NoConvergence(solver_control.last_step(),
184 solver_control.last_value()));
185
186 // Check if the solution is on a CUDA device, if so, copy it over to the
187 // host.
188 if (executor != executor->get_master())
189 {
190 auto x_master = vec::create(executor->get_master(),
191 gko::dim<2>(solution.size(), 1),
192 val_array::view(executor,
193 solution.size(),
194 x->get_values()),
195 1);
196 x.reset(x_master.release());
197 }
198 // Finally copy over the solution vector to deal.II's solution vector.
199 std::copy(x->get_values(),
200 x->get_values() + solution.size(),
201 solution.begin());
202 }
203
204
205
206 template <typename ValueType, typename IndexType>
209 {
210 return solver_control;
211 }
212
213
214
215 template <typename ValueType, typename IndexType>
216 void
218 const SparseMatrix<ValueType> &matrix)
219 {
220 // Needs to be a square matrix
221 Assert(matrix.m() == matrix.n(), ExcNotQuadratic());
222
223 using size_type = ::types::global_dof_index;
224 const size_type N = matrix.m();
225
226 using mtx = gko::matrix::Csr<ValueType, IndexType>;
227 std::shared_ptr<mtx> system_matrix_compute;
228 system_matrix_compute = mtx::create(executor->get_master(),
229 gko::dim<2>(N),
230 matrix.n_nonzero_elements());
231 ValueType *mat_values = system_matrix_compute->get_values();
232 IndexType *mat_row_ptrs = system_matrix_compute->get_row_ptrs();
233 IndexType *mat_col_idxs = system_matrix_compute->get_col_idxs();
234
235 // Copy over the data from the matrix to the data structures Ginkgo needs.
236 //
237 // Final note: if the matrix has entries in the sparsity pattern that are
238 // actually occupied by entries that have a zero numerical value, then we
239 // keep them anyway. people are supposed to provide accurate sparsity
240 // patterns.
241
242 // first fill row lengths array
243 mat_row_ptrs[0] = 0;
244 for (size_type row = 1; row <= N; ++row)
245 mat_row_ptrs[row] =
246 mat_row_ptrs[row - 1] + matrix.get_row_length(row - 1);
247
248 // Copy over matrix elements. note that for sparse matrices,
249 // iterators are sorted so that they traverse each row from start to end
250 // before moving on to the next row. however, this isn't true for block
251 // matrices, so we have to do a bit of book keeping
252 {
253 // Have an array that for each row points to the first entry not yet
254 // written to
255 std::vector<IndexType> row_pointers(N + 1);
256 std::copy(system_matrix_compute->get_row_ptrs(),
257 system_matrix_compute->get_row_ptrs() + N + 1,
258 row_pointers.begin());
259
260 // Loop over the elements of the matrix row by row, as suggested in the
261 // documentation of the sparse matrix iterator class
262 for (size_type row = 0; row < N; ++row)
263 {
265 matrix.begin(row);
266 p != matrix.end(row);
267 ++p)
268 {
269 // Write entry into the first free one for this row
270 mat_col_idxs[row_pointers[row]] = p->column();
271 mat_values[row_pointers[row]] = p->value();
272
273 // Then move pointer ahead
274 ++row_pointers[row];
275 }
276 }
277
278 // At the end, we should have written all rows completely
279 for (size_type i = 0; i < N - 1; ++i)
280 Assert(row_pointers[i] == mat_row_ptrs[i + 1], ExcInternalError());
281 }
282 system_matrix =
283 mtx::create(executor, gko::dim<2>(N), matrix.n_nonzero_elements());
284 system_matrix->copy_from(system_matrix_compute.get());
285 }
286
287
288
289 template <typename ValueType, typename IndexType>
290 void
292 Vector<ValueType> &solution,
293 const Vector<ValueType> &rhs)
294 {
295 initialize(matrix);
296 apply(solution, rhs);
297 }
298
299
300
301 /* ---------------------- SolverCG ------------------------ */
302 template <typename ValueType, typename IndexType>
304 const std::string &exec_type,
305 const AdditionalData &data)
306 : SolverBase<ValueType, IndexType>(solver_control, exec_type)
307 , additional_data(data)
308 {
309 using cg = gko::solver::Cg<ValueType>;
310 this->solver_gen =
311 cg::build().with_criteria(this->combined_factory).on(this->executor);
312 }
313
314
315
316 template <typename ValueType, typename IndexType>
318 SolverControl &solver_control,
319 const std::string &exec_type,
320 const std::shared_ptr<gko::LinOpFactory> &preconditioner,
321 const AdditionalData &data)
322 : SolverBase<ValueType, IndexType>(solver_control, exec_type)
323 , additional_data(data)
324 {
325 using cg = gko::solver::Cg<ValueType>;
326 this->solver_gen = cg::build()
327 .with_criteria(this->combined_factory)
328 .with_preconditioner(preconditioner)
329 .on(this->executor);
330 }
331
332
333
334 /* ---------------------- SolverBicgstab ------------------------ */
335 template <typename ValueType, typename IndexType>
337 SolverControl &solver_control,
338 const std::string &exec_type,
339 const AdditionalData &data)
340 : SolverBase<ValueType, IndexType>(solver_control, exec_type)
341 , additional_data(data)
342 {
343 using bicgstab = gko::solver::Bicgstab<ValueType>;
344 this->solver_gen = bicgstab::build()
345 .with_criteria(this->combined_factory)
346 .on(this->executor);
347 }
348
349
350
351 template <typename ValueType, typename IndexType>
353 SolverControl &solver_control,
354 const std::string &exec_type,
355 const std::shared_ptr<gko::LinOpFactory> &preconditioner,
356 const AdditionalData &data)
357 : SolverBase<ValueType, IndexType>(solver_control, exec_type)
358 , additional_data(data)
359 {
360 using bicgstab = gko::solver::Bicgstab<ValueType>;
361 this->solver_gen = bicgstab::build()
362 .with_criteria(this->combined_factory)
363 .with_preconditioner(preconditioner)
364 .on(this->executor);
365 }
366
367
368
369 /* ---------------------- SolverCGS ------------------------ */
370 template <typename ValueType, typename IndexType>
372 const std::string &exec_type,
373 const AdditionalData &data)
374 : SolverBase<ValueType, IndexType>(solver_control, exec_type)
375 , additional_data(data)
376 {
377 using cgs = gko::solver::Cgs<ValueType>;
378 this->solver_gen =
379 cgs::build().with_criteria(this->combined_factory).on(this->executor);
380 }
381
382
383
384 template <typename ValueType, typename IndexType>
386 SolverControl &solver_control,
387 const std::string &exec_type,
388 const std::shared_ptr<gko::LinOpFactory> &preconditioner,
389 const AdditionalData &data)
390 : SolverBase<ValueType, IndexType>(solver_control, exec_type)
391 , additional_data(data)
392 {
393 using cgs = gko::solver::Cgs<ValueType>;
394 this->solver_gen = cgs::build()
395 .with_criteria(this->combined_factory)
396 .with_preconditioner(preconditioner)
397 .on(this->executor);
398 }
399
400
401
402 /* ---------------------- SolverFCG ------------------------ */
403 template <typename ValueType, typename IndexType>
405 const std::string &exec_type,
406 const AdditionalData &data)
407 : SolverBase<ValueType, IndexType>(solver_control, exec_type)
408 , additional_data(data)
409 {
410 using fcg = gko::solver::Fcg<ValueType>;
411 this->solver_gen =
412 fcg::build().with_criteria(this->combined_factory).on(this->executor);
413 }
414
415
416
417 template <typename ValueType, typename IndexType>
419 SolverControl &solver_control,
420 const std::string &exec_type,
421 const std::shared_ptr<gko::LinOpFactory> &preconditioner,
422 const AdditionalData &data)
423 : SolverBase<ValueType, IndexType>(solver_control, exec_type)
424 , additional_data(data)
425 {
426 using fcg = gko::solver::Fcg<ValueType>;
427 this->solver_gen = fcg::build()
428 .with_criteria(this->combined_factory)
429 .with_preconditioner(preconditioner)
430 .on(this->executor);
431 }
432
433
434
435 /* ---------------------- SolverGMRES ------------------------ */
436 template <typename ValueType, typename IndexType>
438 const unsigned int restart_parameter)
439 : restart_parameter(restart_parameter)
440 {}
441
442
443
444 template <typename ValueType, typename IndexType>
446 const std::string &exec_type,
447 const AdditionalData &data)
448 : SolverBase<ValueType, IndexType>(solver_control, exec_type)
449 , additional_data(data)
450 {
451 using gmres = gko::solver::Gmres<ValueType>;
452 this->solver_gen = gmres::build()
453 .with_krylov_dim(additional_data.restart_parameter)
454 .with_criteria(this->combined_factory)
455 .on(this->executor);
456 }
457
458
459
460 template <typename ValueType, typename IndexType>
462 SolverControl &solver_control,
463 const std::string &exec_type,
464 const std::shared_ptr<gko::LinOpFactory> &preconditioner,
465 const AdditionalData &data)
466 : SolverBase<ValueType, IndexType>(solver_control, exec_type)
467 , additional_data(data)
468 {
469 using gmres = gko::solver::Gmres<ValueType>;
470 this->solver_gen = gmres::build()
471 .with_krylov_dim(additional_data.restart_parameter)
472 .with_criteria(this->combined_factory)
473 .with_preconditioner(preconditioner)
474 .on(this->executor);
475 }
476
477
478
479 /* ---------------------- SolverIR ------------------------ */
480 template <typename ValueType, typename IndexType>
482 const std::string &exec_type,
483 const AdditionalData &data)
484 : SolverBase<ValueType, IndexType>(solver_control, exec_type)
485 , additional_data(data)
486 {
487 using ir = gko::solver::Ir<ValueType>;
488 this->solver_gen =
489 ir::build().with_criteria(this->combined_factory).on(this->executor);
490 }
491
492
493
494 template <typename ValueType, typename IndexType>
496 SolverControl &solver_control,
497 const std::string &exec_type,
498 const std::shared_ptr<gko::LinOpFactory> &inner_solver,
499 const AdditionalData &data)
500 : SolverBase<ValueType, IndexType>(solver_control, exec_type)
501 , additional_data(data)
502 {
503 using ir = gko::solver::Ir<ValueType>;
504 this->solver_gen = ir::build()
505 .with_criteria(this->combined_factory)
506 .with_solver(inner_solver)
507 .on(this->executor);
508 }
509
510
511
512 // Explicit instantiations in GinkgoWrappers
513# define DEALII_INSTANTIATE_FOR_EACH_VALUE_AND_INDEX_TYPE(_macro) \
514 template _macro(float, int32_t); \
515 template _macro(double, int32_t); \
516 template _macro(float, int64_t); \
517 template _macro(double, int64_t);
518
519# define DECLARE_SOLVER_BASE(ValueType, IndexType) \
520 class SolverBase<ValueType, IndexType>
522# undef DECLARE_SOLVER_BASE
523
524# define DECLARE_SOLVER_CG(ValueType, IndexType) \
525 class SolverCG<ValueType, IndexType>
527# undef DECLARE_SOLVER_CG
528
529# define DECLARE_SOLVER_Bicgstab(ValueType, IndexType) \
530 class SolverBicgstab<ValueType, IndexType>
532# undef DECLARE_SOLVER_Bicgstab
533
534# define DECLARE_SOLVER_CGS(ValueType, IndexType) \
535 class SolverCGS<ValueType, IndexType>
537# undef DECLARE_SOLVER_CGS
538
539# define DECLARE_SOLVER_FCG(ValueType, IndexType) \
540 class SolverFCG<ValueType, IndexType>
542# undef DECLARE_SOLVER_FCG
543
544# define DECLARE_SOLVER_GMRES(ValueType, IndexType) \
545 class SolverGMRES<ValueType, IndexType>
547# undef DECLARE_SOLVER_GMRES
548
549# define DECLARE_SOLVER_IR(ValueType, IndexType) \
550 class SolverIR<ValueType, IndexType>
552# undef DECLARE_SOLVER_IR
553
554} // namespace GinkgoWrappers
555
556
558
559#endif // DEAL_II_WITH_GINKGO
SolverControl & control() const
const std::string exec_type
SolverControl & solver_control
void initialize(const SparseMatrix< ValueType > &matrix)
SolverBase(SolverControl &solver_control, const std::string &exec_type)
void apply(Vector< ValueType > &solution, const Vector< ValueType > &rhs)
std::shared_ptr< gko::Executor > executor
std::shared_ptr< gko::LinOpFactory > solver_gen
void solve(const SparseMatrix< ValueType > &matrix, Vector< ValueType > &solution, const Vector< ValueType > &rhs)
std::shared_ptr< gko::stop::ResidualNormReduction<>::Factory > residual_criterion
std::shared_ptr< gko::stop::Combined::Factory > combined_factory
SolverBicgstab(SolverControl &solver_control, const std::string &exec_type, const AdditionalData &data=AdditionalData())
SolverCGS(SolverControl &solver_control, const std::string &exec_type, const AdditionalData &data=AdditionalData())
SolverCG(SolverControl &solver_control, const std::string &exec_type, const AdditionalData &data=AdditionalData())
SolverFCG(SolverControl &solver_control, const std::string &exec_type, const AdditionalData &data=AdditionalData())
const AdditionalData additional_data
SolverGMRES(SolverControl &solver_control, const std::string &exec_type, const AdditionalData &data=AdditionalData())
SolverIR(SolverControl &solver_control, const std::string &exec_type, const AdditionalData &data=AdditionalData())
unsigned int max_steps() const
double tolerance() const
@ success
Stop iteration, goal reached.
virtual size_type size() const override
iterator begin()
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:502
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:503
#define DECLARE_SOLVER_Bicgstab(ValueType, IndexType)
#define DECLARE_SOLVER_FCG(ValueType, IndexType)
#define DECLARE_SOLVER_CG(ValueType, IndexType)
#define DECLARE_SOLVER_BASE(ValueType, IndexType)
#define DECLARE_SOLVER_IR(ValueType, IndexType)
#define DEALII_INSTANTIATE_FOR_EACH_VALUE_AND_INDEX_TYPE(_macro)
#define DECLARE_SOLVER_GMRES(ValueType, IndexType)
#define DECLARE_SOLVER_CGS(ValueType, IndexType)
#define Assert(cond, exc)
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcDivideByZero()
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
static ::ExceptionBase & ExcNotInitialized()
static ::ExceptionBase & ExcNotQuadratic()
static ::ExceptionBase & ExcMessage(std::string arg1)
#define AssertThrow(cond, exc)
unsigned int global_dof_index
Definition types.h:81
AdditionalData(const unsigned int restart_parameter=30)