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