Reference documentation for deal.II version 9.1.1
\(\newcommand{\dealcoloneq}{\mathrel{\vcenter{:}}=}\)
ginkgo_solver.cc
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2018 - 2019 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 
17 #include <deal.II/base/logstream.h>
18 #include <deal.II/base/std_cxx14/memory.h>
19 
20 #include <deal.II/lac/ginkgo_solver.h>
21 
22 #ifdef DEAL_II_WITH_GINKGO
23 
24 # include <deal.II/lac/exceptions.h>
25 
26 # include <cmath>
27 
28 
29 DEAL_II_NAMESPACE_OPEN
30 
31 namespace 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,
55  ExcMessage(
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
88  SolverBase<ValueType, IndexType>::apply(Vector<ValueType> & solution,
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_master`
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_master =
145  gko::matrix::Dense<ValueType>::create(executor->get_master(),
146  gko::dim<2>{1, 1});
147  residual_norm_d_master->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_master` 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_master->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>
208  SolverControl &
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>
522  DEALII_INSTANTIATE_FOR_EACH_VALUE_AND_INDEX_TYPE(DECLARE_SOLVER_BASE)
523 # undef DECLARE_SOLVER_BASE
524 
525 # define DECLARE_SOLVER_CG(ValueType, IndexType) \
526  class SolverCG<ValueType, IndexType>
527  DEALII_INSTANTIATE_FOR_EACH_VALUE_AND_INDEX_TYPE(DECLARE_SOLVER_CG)
528 # undef DECLARE_SOLVER_CG
529 
530 # define DECLARE_SOLVER_Bicgstab(ValueType, IndexType) \
531  class SolverBicgstab<ValueType, IndexType>
532  DEALII_INSTANTIATE_FOR_EACH_VALUE_AND_INDEX_TYPE(DECLARE_SOLVER_Bicgstab)
533 # undef DECLARE_SOLVER_Bicgstab
534 
535 # define DECLARE_SOLVER_CGS(ValueType, IndexType) \
536  class SolverCGS<ValueType, IndexType>
537  DEALII_INSTANTIATE_FOR_EACH_VALUE_AND_INDEX_TYPE(DECLARE_SOLVER_CGS)
538 # undef DECLARE_SOLVER_CGS
539 
540 # define DECLARE_SOLVER_FCG(ValueType, IndexType) \
541  class SolverFCG<ValueType, IndexType>
542  DEALII_INSTANTIATE_FOR_EACH_VALUE_AND_INDEX_TYPE(DECLARE_SOLVER_FCG)
543 # undef DECLARE_SOLVER_FCG
544 
545 # define DECLARE_SOLVER_GMRES(ValueType, IndexType) \
546  class SolverGMRES<ValueType, IndexType>
547  DEALII_INSTANTIATE_FOR_EACH_VALUE_AND_INDEX_TYPE(DECLARE_SOLVER_GMRES)
548 # undef DECLARE_SOLVER_GMRES
549 
550 # define DECLARE_SOLVER_IR(ValueType, IndexType) \
551  class SolverIR<ValueType, IndexType>
552  DEALII_INSTANTIATE_FOR_EACH_VALUE_AND_INDEX_TYPE(DECLARE_SOLVER_IR)
553 # undef DECLARE_SOLVER_IR
554 
555 } // namespace GinkgoWrappers
556 
557 
558 DEAL_II_NAMESPACE_CLOSE
559 
560 #endif // DEAL_II_WITH_GINKGO
std::shared_ptr< gko::stop::ResidualNormReduction<>::Factory > residual_criterion
const std::string exec_type
SolverBicgstab(SolverControl &solver_control, const std::string &exec_type, const AdditionalData &data=AdditionalData())
static ::ExceptionBase & ExcNotInitialized()
#define AssertThrow(cond, exc)
Definition: exceptions.h:1519
SolverCG(SolverControl &solver_control, const std::string &exec_type, const AdditionalData &data=AdditionalData())
static ::ExceptionBase & ExcDivideByZero()
SolverGMRES(SolverControl &solver_control, const std::string &exec_type, const AdditionalData &data=AdditionalData())
static ::ExceptionBase & ExcMessage(std::string arg1)
Stop iteration, goal reached.
#define Assert(cond, exc)
Definition: exceptions.h:1407
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
const AdditionalData additional_data
void initialize(const SparseMatrix< ValueType > &matrix)
SolverControl & solver_control
void solve(const SparseMatrix< ValueType > &matrix, Vector< ValueType > &solution, const Vector< ValueType > &rhs)
AdditionalData(const unsigned int restart_parameter=30)
SolverFCG(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())
double tolerance() const
static ::ExceptionBase & ExcNotQuadratic()
unsigned int global_dof_index
Definition: types.h:89
std::shared_ptr< gko::LinOpFactory > solver_gen
void apply(Vector< ValueType > &solution, const Vector< ValueType > &rhs)
unsigned int max_steps() const
std::shared_ptr< gko::Executor > executor
std::shared_ptr< gko::stop::Combined::Factory > combined_factory
SolverControl & control() const
SolverBase(SolverControl &solver_control, const std::string &exec_type)
static ::ExceptionBase & ExcInternalError()
SolverCGS(SolverControl &solver_control, const std::string &exec_type, const AdditionalData &data=AdditionalData())