MFEM  v4.4.0
Finite element discretization library
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Pages
ginkgo.cpp
Go to the documentation of this file.
1 // Copyright (c) 2010-2022, Lawrence Livermore National Security, LLC. Produced
2 // at the Lawrence Livermore National Laboratory. All Rights reserved. See files
3 // LICENSE and NOTICE for details. LLNL-CODE-806117.
4 //
5 // This file is part of the MFEM library. For more information and source code
6 // availability visit https://mfem.org.
7 //
8 // MFEM is free software; you can redistribute it and/or modify it under the
9 // terms of the BSD-3 license. We welcome feedback and contributions, see file
10 // CONTRIBUTING.md for details.
11 
12 #include "../config/config.hpp"
13 
14 #ifdef MFEM_USE_GINKGO
15 
16 #include "ginkgo.hpp"
17 #include "sparsemat.hpp"
18 #include "../general/globals.hpp"
19 #include "../general/error.hpp"
20 #include <iostream>
21 #include <iomanip>
22 #include <algorithm>
23 #include <cmath>
24 
25 namespace mfem
26 {
27 
28 namespace Ginkgo
29 {
30 
32 {
33  switch (exec_type)
34  {
36  {
37  executor = gko::ReferenceExecutor::create();
38  break;
39  }
41  {
42  executor = gko::OmpExecutor::create();
43  break;
44  }
46  {
47  if (gko::CudaExecutor::get_num_devices() > 0)
48  {
49 #ifdef MFEM_USE_CUDA
50  int current_device = 0;
51  MFEM_GPU_CHECK(cudaGetDevice(&current_device));
52  executor = gko::CudaExecutor::create(current_device,
53  gko::OmpExecutor::create());
54 #endif
55  }
56  else
57  MFEM_ABORT("gko::CudaExecutor::get_num_devices() did not report "
58  "any valid devices.");
59  break;
60  }
62  {
63  if (gko::HipExecutor::get_num_devices() > 0)
64  {
65 #ifdef MFEM_USE_HIP
66  int current_device = 0;
67  MFEM_GPU_CHECK(hipGetDevice(&current_device));
68  executor = gko::HipExecutor::create(current_device,
69  gko::OmpExecutor::create());
70 #endif
71  }
72  else
73  mfem::err << "gko::HipExecutor::get_num_devices() did not report "
74  << "any valid devices" << std::endl;
75  break;
76  }
77  default:
78  mfem::err << "Invalid ExecType specified" << std::endl;
79  }
80 }
81 
83 {
84 
85  // Pick "best match" Executor based on MFEM device configuration.
86  if (mfem_device.Allows(Backend::CUDA_MASK))
87  {
88  if (gko::CudaExecutor::get_num_devices() > 0)
89  {
90 #ifdef MFEM_USE_CUDA
91  int current_device = 0;
92  MFEM_GPU_CHECK(cudaGetDevice(&current_device));
93  executor = gko::CudaExecutor::create(current_device,
94  gko::OmpExecutor::create());
95 #endif
96  }
97  else
98  MFEM_ABORT("gko::CudaExecutor::get_num_devices() did not report "
99  "any valid devices.");
100  }
101  else if (mfem_device.Allows(Backend::HIP_MASK))
102  {
103  if (gko::HipExecutor::get_num_devices() > 0)
104  {
105 #ifdef MFEM_USE_HIP
106  int current_device = 0;
107  MFEM_GPU_CHECK(hipGetDevice(&current_device));
108  executor = gko::HipExecutor::create(current_device, gko::OmpExecutor::create());
109 #endif
110  }
111  else
112  MFEM_ABORT("gko::HipExecutor::get_num_devices() did not report "
113  "any valid devices.");
114  }
115  else
116  {
117  executor = gko::OmpExecutor::create();
118  }
119 }
120 
122  bool use_implicit_res_norm)
123  : Solver(),
124  use_implicit_res_norm(use_implicit_res_norm)
125 {
126  executor = exec.GetExecutor();
127  print_level = -1;
128 
129  // Build default stopping criterion factory
130  max_iter = 10;
131  rel_tol = 0.0;
132  abs_tol = 0.0;
133  this->update_stop_factory();
134 
135  needs_wrapped_vecs = false;
137 }
138 
140 {
141  using ResidualCriterionFactory = gko::stop::ResidualNorm<>;
142  using ImplicitResidualCriterionFactory = gko::stop::ImplicitResidualNorm<>;
143 
145  {
146  imp_rel_criterion = ImplicitResidualCriterionFactory::build()
147  .with_reduction_factor(sqrt(rel_tol))
148  .with_baseline(gko::stop::mode::initial_resnorm)
149  .on(executor);
150  imp_abs_criterion = ImplicitResidualCriterionFactory::build()
151  .with_reduction_factor(sqrt(abs_tol))
152  .with_baseline(gko::stop::mode::absolute)
153  .on(executor);
155  gko::stop::Combined::build()
156  .with_criteria(imp_rel_criterion,
158  gko::stop::Iteration::build()
159  .with_max_iters(static_cast<unsigned long>(max_iter))
160  .on(executor))
161  .on(executor);
162  }
163  else
164  {
165  rel_criterion = ResidualCriterionFactory::build()
166  .with_reduction_factor(rel_tol)
167  .with_baseline(gko::stop::mode::initial_resnorm)
168  .on(executor);
169  abs_criterion = ResidualCriterionFactory::build()
170  .with_reduction_factor(abs_tol)
171  .with_baseline(gko::stop::mode::absolute)
172  .on(executor);
174  gko::stop::Combined::build()
175  .with_criteria(rel_criterion,
177  gko::stop::Iteration::build()
178  .with_max_iters(static_cast<unsigned long>(max_iter))
179  .on(executor))
180  .on(executor);
181  }
182 }
183 
184 void
185 GinkgoIterativeSolver::initialize_ginkgo_log(gko::matrix::Dense<double>* b)
186 const
187 {
188  // Add the logger object. See the different masks available in Ginkgo's
189  // documentation
190  convergence_logger = gko::log::Convergence<>::create(
191  executor, gko::log::Logger::criterion_check_completed_mask);
192  residual_logger = std::make_shared<ResidualLogger<>>(executor,
193  gko::lend(system_oper),b);
194 
195 }
196 
197 void OperatorWrapper::apply_impl(const gko::LinOp *b, gko::LinOp *x) const
198 {
199 
200  // Cast to VectorWrapper; only accept this type for this impl
201  const VectorWrapper *mfem_b = gko::as<const VectorWrapper>(b);
202  VectorWrapper *mfem_x = gko::as<VectorWrapper>(x);
203 
204  this->wrapped_oper->Mult(mfem_b->get_mfem_vec_const_ref(),
205  mfem_x->get_mfem_vec_ref());
206 }
207 void OperatorWrapper::apply_impl(const gko::LinOp *alpha,
208  const gko::LinOp *b,
209  const gko::LinOp *beta,
210  gko::LinOp *x) const
211 {
212  // x = alpha * op (b) + beta * x
213  // Cast to VectorWrapper; only accept this type for this impl
214  const VectorWrapper *mfem_b = gko::as<const VectorWrapper>(b);
215  VectorWrapper *mfem_x = gko::as<VectorWrapper>(x);
216 
217  // Check that alpha and beta are Dense<double> of size (1,1):
218  if (alpha->get_size()[0] > 1 || alpha->get_size()[1] > 1)
219  {
220  throw gko::BadDimension(
221  __FILE__, __LINE__, __func__, "alpha", alpha->get_size()[0],
222  alpha->get_size()[1],
223  "Expected an object of size [1 x 1] for scaling "
224  " in this operator's apply_impl");
225  }
226  if (beta->get_size()[0] > 1 || beta->get_size()[1] > 1)
227  {
228  throw gko::BadDimension(
229  __FILE__, __LINE__, __func__, "beta", beta->get_size()[0],
230  beta->get_size()[1],
231  "Expected an object of size [1 x 1] for scaling "
232  " in this operator's apply_impl");
233  }
234  double alpha_f;
235  double beta_f;
236 
237  if (alpha->get_executor() == alpha->get_executor()->get_master())
238  {
239  // Access value directly
240  alpha_f = gko::as<gko::matrix::Dense<double>>(alpha)->at(0, 0);
241  }
242  else
243  {
244  // Copy from device to host
245  this->get_executor()->get_master().get()->copy_from(
246  this->get_executor().get(),
247  1, gko::as<gko::matrix::Dense<double>>(alpha)->get_const_values(),
248  &alpha_f);
249  }
250  if (beta->get_executor() == beta->get_executor()->get_master())
251  {
252  // Access value directly
253  beta_f = gko::as<gko::matrix::Dense<double>>(beta)->at(0, 0);
254  }
255  else
256  {
257  // Copy from device to host
258  this->get_executor()->get_master().get()->copy_from(
259  this->get_executor().get(),
260  1, gko::as<gko::matrix::Dense<double>>(beta)->get_const_values(),
261  &beta_f);
262  }
263  // Scale x by beta
264  mfem_x->get_mfem_vec_ref() *= beta_f;
265  // Multiply operator with b and store in tmp
266  mfem::Vector mfem_tmp =
267  mfem::Vector(mfem_x->get_size()[0],
268  mfem_x->get_mfem_vec_ref().GetMemory().GetMemoryType());
269  // Set UseDevice flag to match mfem_x (not automatically done through
270  // MemoryType)
271  mfem_tmp.UseDevice(mfem_x->get_mfem_vec_ref().UseDevice());
272 
273  // Apply the operator
274  this->wrapped_oper->Mult(mfem_b->get_mfem_vec_const_ref(), mfem_tmp);
275  // Scale tmp by alpha and add
276  mfem_x->get_mfem_vec_ref().Add(alpha_f, mfem_tmp);
277 
278  mfem_tmp.Destroy();
279 }
280 
281 void
283 {
284 
285  MFEM_VERIFY(system_oper, "System matrix or operator not initialized");
286  MFEM_VERIFY(executor, "executor is not initialized");
287  MFEM_VERIFY(y.Size() == x.Size(),
288  "Mismatching sizes for rhs and solution");
289 
290  using vec = gko::matrix::Dense<double>;
291  if (!iterative_mode)
292  {
293  y = 0.0;
294  }
295 
296  // Create x and y vectors in Ginkgo's format. Wrap MFEM's data directly,
297  // on CPU or GPU.
298  bool on_device = false;
299  if (executor != executor->get_master())
300  {
301  on_device = true;
302  }
303  std::unique_ptr<vec> gko_x;
304  std::unique_ptr<vec> gko_y;
305 
306  // If we do not have an OperatorWrapper for the system operator or
307  // preconditioner, or have an inner solver using VectorWrappers (as
308  // for IR), then directly create Ginkgo vectors from MFEM's data.
309  if (!needs_wrapped_vecs)
310  {
311  gko_x = vec::create(executor, gko::dim<2>(x.Size(), 1),
312  gko::Array<double>::view(executor,
313  x.Size(), const_cast<double *>(
314  x.Read(on_device))), 1);
315  gko_y = vec::create(executor, gko::dim<2>(y.Size(), 1),
316  gko::Array<double>::view(executor,
317  y.Size(),
318  y.ReadWrite(on_device)), 1);
319  }
320  else // We have at least one wrapped MFEM operator; need wrapped vectors
321  {
322  gko_x = std::unique_ptr<vec>(
323  new VectorWrapper(executor, x.Size(),
324  const_cast<Vector *>(&x), false));
325  gko_y = std::unique_ptr<vec>(
326  new VectorWrapper(executor, y.Size(), &y,
327  false));
328  }
329 
330  // Create the logger object to log some data from the solvers to confirm
331  // convergence.
332  initialize_ginkgo_log(gko::lend(gko_x));
333 
334  MFEM_VERIFY(convergence_logger, "convergence logger not initialized" );
335  if (print_level==1)
336  {
337  MFEM_VERIFY(residual_logger, "residual logger not initialized" );
338  solver->clear_loggers(); // Clear any loggers from previous Mult() calls
339  solver->add_logger(residual_logger);
340  }
341 
342  // Add the convergence logger object to the combined factory to retrieve the
343  // solver and other data
344  combined_factory->clear_loggers();
346 
347  // Finally, apply the solver to x and get the solution in y.
348  solver->apply(gko::lend(gko_x), gko::lend(gko_y));
349 
350  // Get the number of iterations taken to converge to the solution.
351  final_iter = convergence_logger->get_num_iterations();
352 
353  // Some residual norm and convergence print outs.
354  double final_res_norm = 0.0;
355 
356  // The convergence_logger object contains the residual vector after the
357  // solver has returned. use this vector to compute the residual norm of the
358  // solution. Get the residual norm from the logger. As the convergence logger
359  // returns a `linop`, it is necessary to convert it to a Dense matrix.
360  // Additionally, if the logger is logging on the gpu, it is necessary to copy
361  // the data to the host and hence the `residual_norm_d_master`
362  auto residual_norm = convergence_logger->get_residual_norm();
363  auto imp_residual_norm = convergence_logger->get_implicit_sq_resnorm();
364 
366  {
367  auto imp_residual_norm_d =
368  gko::as<gko::matrix::Dense<double>>(imp_residual_norm);
369  auto imp_residual_norm_d_master =
370  gko::matrix::Dense<double>::create(executor->get_master(),
371  gko::dim<2> {1, 1});
372  imp_residual_norm_d_master->copy_from(imp_residual_norm_d);
373 
374  final_res_norm = imp_residual_norm_d_master->at(0,0);
375  }
376  else
377  {
378  auto residual_norm_d =
379  gko::as<gko::matrix::Dense<double>>(residual_norm);
380  auto residual_norm_d_master =
381  gko::matrix::Dense<double>::create(executor->get_master(),
382  gko::dim<2> {1, 1});
383  residual_norm_d_master->copy_from(residual_norm_d);
384 
385  final_res_norm = residual_norm_d_master->at(0,0);
386  }
387 
388  converged = 0;
389  if (convergence_logger->has_converged())
390  {
391  converged = 1;
392  }
393 
394  if (print_level == 1)
395  {
396  residual_logger->write();
397  }
398  if (converged == 0)
399  {
400  mfem::err << "No convergence!" << '\n';
401  }
402  if (print_level >=2 && converged==1 )
403  {
404  mfem::out << "Converged in " << final_iter <<
405  " iterations with final residual norm "
406  << final_res_norm << '\n';
407  }
408 }
409 
411 {
412 
413  if (system_oper)
414  {
415  // If the solver currently needs VectorWrappers, but not due to a
416  // "sub-operator" (preconditioner or inner solver), then it's
417  // because the current system_oper needs them. Reset the property
418  // to false in case the new op is a SparseMatrix.
419  if (needs_wrapped_vecs == true && sub_op_needs_wrapped_vecs == false)
420  {
421  needs_wrapped_vecs = false;
422  }
423  // Reset the pointer
424  system_oper.reset();
425  // Reset the solver generated for the previous operator
426  solver.reset();
427  }
428 
429  // Check for SparseMatrix:
430  SparseMatrix *op_mat = const_cast<SparseMatrix*>(
431  dynamic_cast<const SparseMatrix*>(&op));
432  if (op_mat != NULL)
433  {
434  // Needs to be a square matrix
435  MFEM_VERIFY(op_mat->Height() == op_mat->Width(),
436  "System matrix is not square");
437 
438  bool on_device = false;
439  if (executor != executor->get_master())
440  {
441  on_device = true;
442  }
443 
444  using mtx = gko::matrix::Csr<double, int>;
445  const int nnz = op_mat->GetMemoryData().Capacity();
446  system_oper = mtx::create(
447  executor, gko::dim<2>(op_mat->Height(), op_mat->Width()),
448  gko::Array<double>::view(executor,
449  nnz,
450  op_mat->ReadWriteData(on_device)),
451  gko::Array<int>::view(executor,
452  nnz,
453  op_mat->ReadWriteJ(on_device)),
454  gko::Array<int>::view(executor, op_mat->Height() + 1,
455  op_mat->ReadWriteI(on_device)));
456 
457  }
458  else
459  {
460  needs_wrapped_vecs = true;
461  system_oper = std::shared_ptr<OperatorWrapper>(
462  new OperatorWrapper(executor, op.Height(), &op));
463  }
464 
465  // Generate the solver from the solver using the system matrix or operator.
466  solver = solver_gen->generate(system_oper);
467 }
468 
469 /* ---------------------- CGSolver ------------------------ */
471  : EnableGinkgoSolver(exec, true)
472 {
473  using cg = gko::solver::Cg<double>;
474  this->solver_gen =
475  cg::build().with_criteria(this->combined_factory).on(this->executor);
476 }
477 
479  const GinkgoPreconditioner &preconditioner)
480  : EnableGinkgoSolver(exec, true)
481 {
482  using cg = gko::solver::Cg<double>;
483  // Check for a previously-generated preconditioner (for a specific matrix)
484  if (preconditioner.HasGeneratedPreconditioner())
485  {
486  this->solver_gen = cg::build()
487  .with_criteria(this->combined_factory)
488  .with_generated_preconditioner(
489  preconditioner.GetGeneratedPreconditioner())
490  .on(this->executor);
491  if (dynamic_cast<const OperatorWrapper*>(preconditioner.
492  GetGeneratedPreconditioner().get()))
493  {
494  this->sub_op_needs_wrapped_vecs = true;
495  this->needs_wrapped_vecs = true;
496  }
497  }
498  else // Pass a preconditioner factory (will use same matrix as the solver)
499  {
500  this->solver_gen = cg::build()
501  .with_criteria(this->combined_factory)
502  .with_preconditioner(preconditioner.GetFactory())
503  .on(this->executor);
504  }
505 }
506 
507 
508 /* ---------------------- BICGSTABSolver ------------------------ */
510  : EnableGinkgoSolver(exec, true)
511 {
512  using bicgstab = gko::solver::Bicgstab<double>;
513  this->solver_gen = bicgstab::build()
514  .with_criteria(this->combined_factory)
515  .on(this->executor);
516 }
517 
519  const GinkgoPreconditioner &preconditioner)
520  : EnableGinkgoSolver(exec, true)
521 {
522  using bicgstab = gko::solver::Bicgstab<double>;
523  if (preconditioner.HasGeneratedPreconditioner())
524  {
525  this->solver_gen = bicgstab::build()
526  .with_criteria(this->combined_factory)
527  .with_generated_preconditioner(
528  preconditioner.GetGeneratedPreconditioner())
529  .on(this->executor);
530  if (dynamic_cast<const OperatorWrapper*>(preconditioner.
531  GetGeneratedPreconditioner().get()))
532  {
533  this->sub_op_needs_wrapped_vecs = true;
534  this->needs_wrapped_vecs = true;
535  }
536  }
537  else
538  {
539  this->solver_gen = bicgstab::build()
540  .with_criteria(this->combined_factory)
541  .with_preconditioner(preconditioner.GetFactory())
542  .on(this->executor);
543  }
544 }
545 
546 
547 /* ---------------------- CGSSolver ------------------------ */
549  : EnableGinkgoSolver(exec, true)
550 {
551  using cgs = gko::solver::Cgs<double>;
552  this->solver_gen =
553  cgs::build().with_criteria(this->combined_factory).on(this->executor);
554 }
555 
557  const GinkgoPreconditioner &preconditioner)
558  : EnableGinkgoSolver(exec, true)
559 {
560  using cgs = gko::solver::Cgs<double>;
561  if (preconditioner.HasGeneratedPreconditioner())
562  {
563  this->solver_gen = cgs::build()
564  .with_criteria(this->combined_factory)
565  .with_generated_preconditioner(
566  preconditioner.GetGeneratedPreconditioner())
567  .on(this->executor);
568  if (dynamic_cast<const OperatorWrapper*>(preconditioner.
569  GetGeneratedPreconditioner().get()))
570  {
571  this->sub_op_needs_wrapped_vecs = true;
572  this->needs_wrapped_vecs = true;
573  }
574  }
575  else
576  {
577  this->solver_gen = cgs::build()
578  .with_criteria(this->combined_factory)
579  .with_preconditioner(preconditioner.GetFactory())
580  .on(this->executor);
581  }
582 }
583 
584 
585 /* ---------------------- FCGSolver ------------------------ */
587  : EnableGinkgoSolver(exec, true)
588 {
589  using fcg = gko::solver::Fcg<double>;
590  this->solver_gen =
591  fcg::build().with_criteria(this->combined_factory).on(this->executor);
592 }
593 
595  const GinkgoPreconditioner &preconditioner)
596  : EnableGinkgoSolver(exec, true)
597 {
598  using fcg = gko::solver::Fcg<double>;
599  if (preconditioner.HasGeneratedPreconditioner())
600  {
601  this->solver_gen = fcg::build()
602  .with_criteria(this->combined_factory)
603  .with_generated_preconditioner(
604  preconditioner.GetGeneratedPreconditioner())
605  .on(this->executor);
606  if (dynamic_cast<const OperatorWrapper*>(preconditioner.
607  GetGeneratedPreconditioner().get()))
608  {
609  this->sub_op_needs_wrapped_vecs = true;
610  this->needs_wrapped_vecs = true;
611  }
612  }
613  else
614  {
615  this->solver_gen = fcg::build()
616  .with_criteria(this->combined_factory)
617  .with_preconditioner(preconditioner.GetFactory())
618  .on(this->executor);
619  }
620 }
621 
622 
623 /* ---------------------- GMRESSolver ------------------------ */
625  : EnableGinkgoSolver(exec, false),
626  m{dim}
627 {
628  using gmres = gko::solver::Gmres<double>;
629  if (this->m == 0) // Don't set a dimension, but let Ginkgo use its default
630  {
631  this->solver_gen = gmres::build()
632  .with_criteria(this->combined_factory)
633  .on(this->executor);
634  }
635  else
636  {
637  this->solver_gen = gmres::build()
638  .with_krylov_dim(static_cast<unsigned long>(m))
639  .with_criteria(this->combined_factory)
640  .on(this->executor);
641  }
642 }
643 
645  const GinkgoPreconditioner &preconditioner, int dim)
646  : EnableGinkgoSolver(exec, false),
647  m{dim}
648 {
649  using gmres = gko::solver::Gmres<double>;
650  // Check for a previously-generated preconditioner (for a specific matrix)
651  if (this->m == 0) // Don't set a dimension, but let Ginkgo use its default
652  {
653  if (preconditioner.HasGeneratedPreconditioner())
654  {
655  this->solver_gen = gmres::build()
656  .with_criteria(this->combined_factory)
657  .with_generated_preconditioner(
658  preconditioner.GetGeneratedPreconditioner())
659  .on(this->executor);
660  if (dynamic_cast<const OperatorWrapper*>(preconditioner.
661  GetGeneratedPreconditioner().get()))
662  {
663  this->sub_op_needs_wrapped_vecs = true;
664  this->needs_wrapped_vecs = true;
665  }
666  }
667  else
668  {
669  this->solver_gen = gmres::build()
670  .with_criteria(this->combined_factory)
671  .with_preconditioner(preconditioner.GetFactory())
672  .on(this->executor);
673  }
674  }
675  else
676  {
677  if (preconditioner.HasGeneratedPreconditioner())
678  {
679  this->solver_gen = gmres::build()
680  .with_krylov_dim(static_cast<unsigned long>(m))
681  .with_criteria(this->combined_factory)
682  .with_generated_preconditioner(
683  preconditioner.GetGeneratedPreconditioner())
684  .on(this->executor);
685  if (dynamic_cast<const OperatorWrapper*>(preconditioner.
686  GetGeneratedPreconditioner().get()))
687  {
688  this->sub_op_needs_wrapped_vecs = true;
689  this->needs_wrapped_vecs = true;
690  }
691  }
692  else
693  {
694  this->solver_gen = gmres::build()
695  .with_krylov_dim(static_cast<unsigned long>(m))
696  .with_criteria(this->combined_factory)
697  .with_preconditioner(preconditioner.GetFactory())
698  .on(this->executor);
699  }
700  }
701 }
702 
704 {
705  m = dim;
706  using gmres_type = gko::solver::Gmres<double>;
707  gko::as<gmres_type::Factory>(solver_gen)->get_parameters().krylov_dim = m;
708  if (solver)
709  {
710  gko::as<gmres_type>(solver)->set_krylov_dim(static_cast<unsigned long>(m));
711  }
712 }
713 
714 /* ---------------------- CBGMRESSolver ------------------------ */
716  storage_precision prec)
717  : EnableGinkgoSolver(exec, false),
718  m{dim}
719 {
720  using gmres = gko::solver::CbGmres<double>;
721  if (this->m == 0) // Don't set a dimension, but let Ginkgo use its default
722  {
723  this->solver_gen = gmres::build()
724  .with_criteria(this->combined_factory)
725  .with_storage_precision(prec)
726  .on(this->executor);
727  }
728  else
729  {
730  this->solver_gen = gmres::build()
731  .with_krylov_dim(static_cast<unsigned long>(m))
732  .with_criteria(this->combined_factory)
733  .with_storage_precision(prec)
734  .on(this->executor);
735  }
736 }
737 
739  const GinkgoPreconditioner &preconditioner,
740  int dim, storage_precision prec)
741  : EnableGinkgoSolver(exec, false),
742  m{dim}
743 {
744  using gmres = gko::solver::CbGmres<double>;
745  // Check for a previously-generated preconditioner (for a specific matrix)
746  if (this->m == 0) // Don't set a dimension, but let Ginkgo use its default
747  {
748  if (preconditioner.HasGeneratedPreconditioner())
749  {
750  this->solver_gen = gmres::build()
751  .with_criteria(this->combined_factory)
752  .with_storage_precision(prec)
753  .with_generated_preconditioner(
754  preconditioner.GetGeneratedPreconditioner())
755  .on(this->executor);
756  if (dynamic_cast<const OperatorWrapper*>(preconditioner.
757  GetGeneratedPreconditioner().get()))
758  {
759  this->sub_op_needs_wrapped_vecs = true;
760  this->needs_wrapped_vecs = true;
761  }
762  }
763  else
764  {
765  this->solver_gen = gmres::build()
766  .with_criteria(this->combined_factory)
767  .with_storage_precision(prec)
768  .with_preconditioner(preconditioner.GetFactory())
769  .on(this->executor);
770  }
771  }
772  else
773  {
774  if (preconditioner.HasGeneratedPreconditioner())
775  {
776  this->solver_gen = gmres::build()
777  .with_krylov_dim(static_cast<unsigned long>(m))
778  .with_criteria(this->combined_factory)
779  .with_storage_precision(prec)
780  .with_generated_preconditioner(
781  preconditioner.GetGeneratedPreconditioner())
782  .on(this->executor);
783  if (dynamic_cast<const OperatorWrapper*>(preconditioner.
784  GetGeneratedPreconditioner().get()))
785  {
786  this->sub_op_needs_wrapped_vecs = true;
787  this->needs_wrapped_vecs = true;
788  }
789  }
790  else
791  {
792  this->solver_gen = gmres::build()
793  .with_krylov_dim(static_cast<unsigned long>(m))
794  .with_criteria(this->combined_factory)
795  .with_storage_precision(prec)
796  .with_preconditioner(preconditioner.GetFactory())
797  .on(this->executor);
798  }
799  }
800 }
801 
803 {
804  m = dim;
805  using gmres_type = gko::solver::CbGmres<double>;
806  gko::as<gmres_type::Factory>(solver_gen)->get_parameters().krylov_dim = m;
807  if (solver)
808  {
809  gko::as<gmres_type>(solver)->set_krylov_dim(static_cast<unsigned long>(m));
810  }
811 }
812 
813 /* ---------------------- IRSolver ------------------------ */
815  : EnableGinkgoSolver(exec, false)
816 {
817  using ir = gko::solver::Ir<double>;
818  this->solver_gen =
819  ir::build().with_criteria(this->combined_factory).on(this->executor);
820 }
821 
823  const GinkgoIterativeSolver &inner_solver)
824  : EnableGinkgoSolver(exec, false)
825 {
826  using ir = gko::solver::Ir<double>;
827  this->solver_gen = ir::build()
828  .with_criteria(this->combined_factory)
829  .with_solver(inner_solver.GetFactory())
830  .on(this->executor);
831  if (inner_solver.UsesVectorWrappers())
832  {
833  this->sub_op_needs_wrapped_vecs = true;
834  this->needs_wrapped_vecs = true;
835  }
836 }
837 
838 /* --------------------------------------------------------------- */
839 /* ---------------------- Preconditioners ------------------------ */
841  GinkgoExecutor &exec)
842  : Solver()
843 {
844  has_generated_precond = false;
845  executor = exec.GetExecutor();
846 }
847 
848 void
850 {
851 
852  MFEM_VERIFY(generated_precond, "Preconditioner not initialized");
853  MFEM_VERIFY(executor, "executor is not initialized");
854 
855  using vec = gko::matrix::Dense<double>;
856  if (!iterative_mode)
857  {
858  y = 0.0;
859  }
860 
861  // Create x and y vectors in Ginkgo's format. Wrap MFEM's data directly,
862  // on CPU or GPU.
863  bool on_device = false;
864  if (executor != executor->get_master())
865  {
866  on_device = true;
867  }
868  auto gko_x = vec::create(executor, gko::dim<2>(x.Size(), 1),
869  gko::Array<double>::view(executor,
870  x.Size(), const_cast<double *>(
871  x.Read(on_device))), 1);
872  auto gko_y = vec::create(executor, gko::dim<2>(y.Size(), 1),
873  gko::Array<double>::view(executor,
874  y.Size(),
875  y.ReadWrite(on_device)), 1);
876  generated_precond.get()->apply(gko::lend(gko_x), gko::lend(gko_y));
877 }
878 
880 {
881 
883  {
884  generated_precond.reset();
885  has_generated_precond = false;
886  }
887 
888  // Only accept SparseMatrix for this type.
889  SparseMatrix *op_mat = const_cast<SparseMatrix*>(
890  dynamic_cast<const SparseMatrix*>(&op));
891  MFEM_VERIFY(op_mat != NULL,
892  "GinkgoPreconditioner::SetOperator : not a SparseMatrix!");
893 
894  bool on_device = false;
895  if (executor != executor->get_master())
896  {
897  on_device = true;
898  }
899 
900  using mtx = gko::matrix::Csr<double, int>;
901  const int nnz = op_mat->GetMemoryData().Capacity();
902  auto gko_matrix = mtx::create(
903  executor, gko::dim<2>(op_mat->Height(), op_mat->Width()),
904  gko::Array<double>::view(executor,
905  nnz,
906  op_mat->ReadWriteData(on_device)),
907  gko::Array<int>::view(executor,
908  nnz,
909  op_mat->ReadWriteJ(on_device)),
910  gko::Array<int>::view(executor, op_mat->Height() + 1,
911  op_mat->ReadWriteI(on_device)));
912 
913  generated_precond = precond_gen->generate(gko::give(gko_matrix));
914  has_generated_precond = true;
915 }
916 
917 
918 /* ---------------------- JacobiPreconditioner ------------------------ */
920  GinkgoExecutor &exec,
921  const std::string &storage_opt,
922  const double accuracy,
923  const int max_block_size
924 )
925  : GinkgoPreconditioner(exec)
926 {
927 
928  if (storage_opt == "auto")
929  {
930  precond_gen = gko::preconditioner::Jacobi<double, int>::build()
931  .with_storage_optimization(
932  gko::precision_reduction::autodetect())
933  .with_accuracy(accuracy)
934  .with_max_block_size(static_cast<unsigned int>(max_block_size))
935  .on(executor);
936  }
937  else
938  {
939  precond_gen = gko::preconditioner::Jacobi<double, int>::build()
940  .with_storage_optimization(
941  gko::precision_reduction(0, 0))
942  .with_accuracy(accuracy)
943  .with_max_block_size(static_cast<unsigned int>(max_block_size))
944  .on(executor);
945  }
946 
947 }
948 
949 /* ---------------------- Ilu/IluIsaiPreconditioner ------------------------ */
951  GinkgoExecutor &exec,
952  const std::string &factorization_type,
953  const int sweeps,
954  const bool skip_sort
955 )
956  : GinkgoPreconditioner(exec)
957 {
958  if (factorization_type == "exact")
959  {
960  using ilu_fact_type = gko::factorization::Ilu<double, int>;
961  std::shared_ptr<ilu_fact_type::Factory> fact_factory =
962  ilu_fact_type::build()
963  .with_skip_sorting(skip_sort)
964  .on(executor);
965  precond_gen = gko::preconditioner::Ilu<>::build()
966  .with_factorization_factory(fact_factory)
967  .on(executor);
968  }
969  else
970  {
971  using ilu_fact_type = gko::factorization::ParIlu<double, int>;
972  std::shared_ptr<ilu_fact_type::Factory> fact_factory =
973  ilu_fact_type::build()
974  .with_iterations(static_cast<unsigned long>(sweeps))
975  .with_skip_sorting(skip_sort)
976  .on(executor);
977  precond_gen = gko::preconditioner::Ilu<>::build()
978  .with_factorization_factory(fact_factory)
979  .on(executor);
980  }
981 
982 }
983 
985  GinkgoExecutor &exec,
986  const std::string &factorization_type,
987  const int sweeps,
988  const int sparsity_power,
989  const bool skip_sort
990 )
991  : GinkgoPreconditioner(exec)
992 {
993  using l_solver_type = gko::preconditioner::LowerIsai<>;
994  using u_solver_type = gko::preconditioner::UpperIsai<>;
995 
996  std::shared_ptr<l_solver_type::Factory> l_solver_factory =
997  l_solver_type::build()
998  .with_sparsity_power(sparsity_power)
999  .on(executor);
1000  std::shared_ptr<u_solver_type::Factory> u_solver_factory =
1001  u_solver_type::build()
1002  .with_sparsity_power(sparsity_power)
1003  .on(executor);
1004 
1005 
1006 
1007  if (factorization_type == "exact")
1008  {
1009  using ilu_fact_type = gko::factorization::Ilu<double, int>;
1010  std::shared_ptr<ilu_fact_type::Factory> fact_factory =
1011  ilu_fact_type::build()
1012  .with_skip_sorting(skip_sort)
1013  .on(executor);
1014  precond_gen = gko::preconditioner::Ilu<l_solver_type,
1015  u_solver_type>::build()
1016  .with_factorization_factory(fact_factory)
1017  .with_l_solver_factory(l_solver_factory)
1018  .with_u_solver_factory(u_solver_factory)
1019  .on(executor);
1020 
1021  }
1022  else
1023  {
1024  using ilu_fact_type = gko::factorization::ParIlu<double, int>;
1025  std::shared_ptr<ilu_fact_type::Factory> fact_factory =
1026  ilu_fact_type::build()
1027  .with_iterations(static_cast<unsigned long>(sweeps))
1028  .with_skip_sorting(skip_sort)
1029  .on(executor);
1030  precond_gen = gko::preconditioner::Ilu<l_solver_type,
1031  u_solver_type>::build()
1032  .with_factorization_factory(fact_factory)
1033  .with_l_solver_factory(l_solver_factory)
1034  .with_u_solver_factory(u_solver_factory)
1035  .on(executor);
1036  }
1037 }
1038 
1039 
1040 /* ---------------------- Ic/IcIsaiPreconditioner ------------------------ */
1042  GinkgoExecutor &exec,
1043  const std::string &factorization_type,
1044  const int sweeps,
1045  const bool skip_sort
1046 )
1047  : GinkgoPreconditioner(exec)
1048 {
1049 
1050  if (factorization_type == "exact")
1051  {
1052  using ic_fact_type = gko::factorization::Ic<double, int>;
1053  std::shared_ptr<ic_fact_type::Factory> fact_factory =
1054  ic_fact_type::build()
1055  .with_both_factors(false)
1056  .with_skip_sorting(skip_sort)
1057  .on(executor);
1058  precond_gen = gko::preconditioner::Ic<>::build()
1059  .with_factorization_factory(fact_factory)
1060  .on(executor);
1061  }
1062  else
1063  {
1064  using ic_fact_type = gko::factorization::ParIc<double, int>;
1065  std::shared_ptr<ic_fact_type::Factory> fact_factory =
1066  ic_fact_type::build()
1067  .with_both_factors(false)
1068  .with_iterations(static_cast<unsigned long>(sweeps))
1069  .with_skip_sorting(skip_sort)
1070  .on(executor);
1071  precond_gen = gko::preconditioner::Ic<>::build()
1072  .with_factorization_factory(fact_factory)
1073  .on(executor);
1074  }
1075 }
1076 
1078  GinkgoExecutor &exec,
1079  const std::string &factorization_type,
1080  const int sweeps,
1081  const int sparsity_power,
1082  const bool skip_sort
1083 )
1084  : GinkgoPreconditioner(exec)
1085 {
1086 
1087  using l_solver_type = gko::preconditioner::LowerIsai<>;
1088  std::shared_ptr<l_solver_type::Factory> l_solver_factory =
1089  l_solver_type::build()
1090  .with_sparsity_power(sparsity_power)
1091  .on(executor);
1092  if (factorization_type == "exact")
1093  {
1094  using ic_fact_type = gko::factorization::Ic<double, int>;
1095  std::shared_ptr<ic_fact_type::Factory> fact_factory =
1096  ic_fact_type::build()
1097  .with_both_factors(false)
1098  .with_skip_sorting(skip_sort)
1099  .on(executor);
1100  precond_gen = gko::preconditioner::Ic<l_solver_type>::build()
1101  .with_factorization_factory(fact_factory)
1102  .with_l_solver_factory(l_solver_factory)
1103  .on(executor);
1104  }
1105  else
1106  {
1107  using ic_fact_type = gko::factorization::ParIc<double, int>;
1108  std::shared_ptr<ic_fact_type::Factory> fact_factory =
1109  ic_fact_type::build()
1110  .with_both_factors(false)
1111  .with_iterations(static_cast<unsigned long>(sweeps))
1112  .with_skip_sorting(skip_sort)
1113  .on(executor);
1114  precond_gen = gko::preconditioner::Ic<l_solver_type>::build()
1115  .with_factorization_factory(fact_factory)
1116  .with_l_solver_factory(l_solver_factory)
1117  .on(executor);
1118  }
1119 }
1120 
1121 /* ---------------------- MFEMPreconditioner ------------------------ */
1123  GinkgoExecutor &exec,
1124  const Solver &mfem_precond
1125 )
1126  : GinkgoPreconditioner(exec)
1127 {
1128  generated_precond = std::shared_ptr<OperatorWrapper>(
1130  mfem_precond.Height(), &mfem_precond));
1131  has_generated_precond = true;
1132 }
1133 
1134 } // namespace Ginkgo
1135 
1136 } // namespace mfem
1137 
1138 #endif // MFEM_USE_GINKGO
std::shared_ptr< gko::log::Convergence<> > convergence_logger
Definition: ginkgo.hpp:709
std::shared_ptr< gko::Executor > GetExecutor() const
Definition: ginkgo.hpp:472
std::shared_ptr< ResidualLogger<> > residual_logger
Definition: ginkgo.hpp:715
virtual void SetOperator(const Operator &op)
Definition: ginkgo.cpp:410
std::shared_ptr< gko::stop::ResidualNorm<>::Factory > abs_criterion
Definition: ginkgo.hpp:689
GinkgoIterativeSolver(GinkgoExecutor &exec, bool use_implicit_res_norm)
Definition: ginkgo.cpp:121
CBGMRESSolver(GinkgoExecutor &exec, int dim=0, storage_precision prec=storage_precision::reduce1)
Definition: ginkgo.cpp:715
std::shared_ptr< gko::LinOpFactory > precond_gen
Definition: ginkgo.hpp:557
int Width() const
Get the width (size of input) of the Operator. Synonym with NumCols().
Definition: operator.hpp:72
CGSSolver(GinkgoExecutor &exec)
Definition: ginkgo.cpp:548
Biwise-OR of all HIP backends.
Definition: device.hpp:90
MFEMPreconditioner(GinkgoExecutor &exec, const Solver &mfem_precond)
Definition: ginkgo.cpp:1122
JacobiPreconditioner(GinkgoExecutor &exec, const std::string &storage_opt="none", const double accuracy=1.e-1, const int max_block_size=32)
Definition: ginkgo.cpp:919
int Size() const
Returns the size of the vector.
Definition: vector.hpp:199
std::shared_ptr< gko::Executor > executor
Definition: ginkgo.hpp:577
GinkgoExecutor(ExecType exec_type)
Definition: ginkgo.cpp:31
virtual void Mult(const Vector &x, Vector &y) const =0
Operator application: y=A(x).
void SetKDim(int dim)
Definition: ginkgo.cpp:703
bool iterative_mode
If true, use the second argument of Mult() as an initial guess.
Definition: operator.hpp:655
virtual void UseDevice(bool use_dev) const
Enable execution of Vector operations using the mfem::Device.
Definition: vector.hpp:117
Memory< double > & GetMemoryData()
Definition: sparsemat.hpp:256
int Capacity() const
Return the size of the allocated memory.
BICGSTABSolver(GinkgoExecutor &exec)
Definition: ginkgo.cpp:509
virtual void Mult(const Vector &x, Vector &y) const
Definition: ginkgo.cpp:849
std::shared_ptr< gko::stop::ImplicitResidualNorm<>::Factory > imp_rel_criterion
Definition: ginkgo.hpp:696
bool HasGeneratedPreconditioner() const
Definition: ginkgo.hpp:554
IluPreconditioner(GinkgoExecutor &exec, const std::string &factorization_type="exact", const int sweeps=0, const bool skip_sort=false)
Definition: ginkgo.cpp:950
Data type sparse matrix.
Definition: sparsemat.hpp:46
int Height() const
Get the height (size of output) of the Operator. Synonym with NumRows().
Definition: operator.hpp:66
double b
Definition: lissajous.cpp:42
const std::shared_ptr< gko::LinOpFactory > GetFactory() const
Definition: ginkgo.hpp:604
int * ReadWriteI(bool on_dev=true)
Definition: sparsemat.hpp:230
std::shared_ptr< gko::Executor > executor
Definition: ginkgo.hpp:728
CGSolver(GinkgoExecutor &exec)
Definition: ginkgo.cpp:470
Biwise-OR of all CUDA backends.
Definition: device.hpp:88
GMRESSolver(GinkgoExecutor &exec, int dim=0)
Definition: ginkgo.cpp:624
OutStream err(std::cerr)
Global stream used by the library for standard error output. Initially it uses the same std::streambu...
Definition: globals.hpp:71
const std::shared_ptr< gko::LinOpFactory > GetFactory() const
Definition: ginkgo.hpp:536
OpenMP CPU Executor.
Definition: ginkgo.hpp:446
std::shared_ptr< gko::stop::ResidualNorm<>::Factory > rel_criterion
Definition: ginkgo.hpp:682
static bool Allows(unsigned long b_mask)
Return true if any of the backends in the backend mask, b_mask, are allowed.
Definition: device.hpp:258
IRSolver(GinkgoExecutor &exec)
Definition: ginkgo.cpp:814
const Vector & get_mfem_vec_const_ref() const
Definition: ginkgo.hpp:110
IcPreconditioner(GinkgoExecutor &exec, const std::string &factorization_type="exact", const int sweeps=0, const bool skip_sort=false)
Definition: ginkgo.cpp:1041
std::shared_ptr< gko::stop::Combined::Factory > combined_factory
Definition: ginkgo.hpp:721
virtual double * ReadWrite(bool on_dev=true)
Shortcut for mfem::ReadWrite(vec.GetMemory(), vec.Size(), on_dev).
Definition: vector.hpp:454
std::shared_ptr< gko::LinOpFactory > solver_gen
Definition: ginkgo.hpp:670
FCGSolver(GinkgoExecutor &exec)
Definition: ginkgo.cpp:586
FDualNumber< tbase > sqrt(const FDualNumber< tbase > &f)
sqrt([dual number])
Definition: fdual.hpp:600
int dim
Definition: ex24.cpp:53
IcIsaiPreconditioner(GinkgoExecutor &exec, const std::string &factorization_type="exact", const int sweeps=0, const int sparsity_power=1, const bool skip_sort=false)
Definition: ginkgo.cpp:1077
void Destroy()
Destroy a vector.
Definition: vector.hpp:598
virtual void Mult(const Vector &x, Vector &y) const
Definition: ginkgo.cpp:282
IluIsaiPreconditioner(GinkgoExecutor &exec, const std::string &factorization_type="exact", const int sweeps=0, const int sparsity_power=1, const bool skip_sort=false)
Definition: ginkgo.cpp:984
GinkgoPreconditioner(GinkgoExecutor &exec)
Definition: ginkgo.cpp:840
std::shared_ptr< gko::LinOp > generated_precond
Definition: ginkgo.hpp:570
const double alpha
Definition: ex15.cpp:369
double * ReadWriteData(bool on_dev=true)
Definition: sparsemat.hpp:262
int * ReadWriteJ(bool on_dev=true)
Definition: sparsemat.hpp:246
Vector data type.
Definition: vector.hpp:60
std::shared_ptr< gko::stop::ImplicitResidualNorm<>::Factory > imp_abs_criterion
Definition: ginkgo.hpp:703
Base class for solvers.
Definition: operator.hpp:651
Reference CPU Executor.
Definition: ginkgo.hpp:444
virtual void SetOperator(const Operator &op)
Definition: ginkgo.cpp:879
OutStream out(std::cout)
Global stream used by the library for standard output. Initially it uses the same std::streambuf as s...
Definition: globals.hpp:66
The MFEM Device class abstracts hardware devices such as GPUs, as well as programming models such as ...
Definition: device.hpp:121
Abstract operator.
Definition: operator.hpp:24
virtual const double * Read(bool on_dev=true) const
Shortcut for mfem::Read(vec.GetMemory(), vec.Size(), on_dev).
Definition: vector.hpp:438
const std::shared_ptr< gko::LinOp > GetGeneratedPreconditioner() const
Definition: ginkgo.hpp:545
void apply_impl(const gko::LinOp *b, gko::LinOp *x) const override
Definition: ginkgo.cpp:197
std::shared_ptr< gko::LinOp > solver
Definition: ginkgo.hpp:675