MFEM  v4.5.2
Finite element discretization library
ginkgo.cpp
Go to the documentation of this file.
1 // Copyright (c) 2010-2023, 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 #if MFEM_GINKGO_VERSION < 10500
191  convergence_logger = gko::log::Convergence<>::create(
192  executor, gko::log::Logger::criterion_check_completed_mask);
193 #else
194  convergence_logger = gko::log::Convergence<>::create(
195  gko::log::Logger::criterion_check_completed_mask);
196 #endif
197  residual_logger = std::make_shared<ResidualLogger<>>(executor,
198  gko::lend(system_oper),b);
199 
200 }
201 
202 void OperatorWrapper::apply_impl(const gko::LinOp *b, gko::LinOp *x) const
203 {
204 
205  // Cast to VectorWrapper; only accept this type for this impl
206  const VectorWrapper *mfem_b = gko::as<const VectorWrapper>(b);
207  VectorWrapper *mfem_x = gko::as<VectorWrapper>(x);
208 
209  this->wrapped_oper->Mult(mfem_b->get_mfem_vec_const_ref(),
210  mfem_x->get_mfem_vec_ref());
211 }
212 void OperatorWrapper::apply_impl(const gko::LinOp *alpha,
213  const gko::LinOp *b,
214  const gko::LinOp *beta,
215  gko::LinOp *x) const
216 {
217  // x = alpha * op (b) + beta * x
218  // Cast to VectorWrapper; only accept this type for this impl
219  const VectorWrapper *mfem_b = gko::as<const VectorWrapper>(b);
220  VectorWrapper *mfem_x = gko::as<VectorWrapper>(x);
221 
222  // Check that alpha and beta are Dense<double> of size (1,1):
223  if (alpha->get_size()[0] > 1 || alpha->get_size()[1] > 1)
224  {
225  throw gko::BadDimension(
226  __FILE__, __LINE__, __func__, "alpha", alpha->get_size()[0],
227  alpha->get_size()[1],
228  "Expected an object of size [1 x 1] for scaling "
229  " in this operator's apply_impl");
230  }
231  if (beta->get_size()[0] > 1 || beta->get_size()[1] > 1)
232  {
233  throw gko::BadDimension(
234  __FILE__, __LINE__, __func__, "beta", beta->get_size()[0],
235  beta->get_size()[1],
236  "Expected an object of size [1 x 1] for scaling "
237  " in this operator's apply_impl");
238  }
239  double alpha_f;
240  double beta_f;
241 
242  if (alpha->get_executor() == alpha->get_executor()->get_master())
243  {
244  // Access value directly
245  alpha_f = gko::as<gko::matrix::Dense<double>>(alpha)->at(0, 0);
246  }
247  else
248  {
249  // Copy from device to host
250  this->get_executor()->get_master().get()->copy_from(
251  this->get_executor().get(),
252  1, gko::as<gko::matrix::Dense<double>>(alpha)->get_const_values(),
253  &alpha_f);
254  }
255  if (beta->get_executor() == beta->get_executor()->get_master())
256  {
257  // Access value directly
258  beta_f = gko::as<gko::matrix::Dense<double>>(beta)->at(0, 0);
259  }
260  else
261  {
262  // Copy from device to host
263  this->get_executor()->get_master().get()->copy_from(
264  this->get_executor().get(),
265  1, gko::as<gko::matrix::Dense<double>>(beta)->get_const_values(),
266  &beta_f);
267  }
268  // Scale x by beta
269  mfem_x->get_mfem_vec_ref() *= beta_f;
270  // Multiply operator with b and store in tmp
271  mfem::Vector mfem_tmp =
272  mfem::Vector(mfem_x->get_size()[0],
273  mfem_x->get_mfem_vec_ref().GetMemory().GetMemoryType());
274  // Set UseDevice flag to match mfem_x (not automatically done through
275  // MemoryType)
276  mfem_tmp.UseDevice(mfem_x->get_mfem_vec_ref().UseDevice());
277 
278  // Apply the operator
279  this->wrapped_oper->Mult(mfem_b->get_mfem_vec_const_ref(), mfem_tmp);
280  // Scale tmp by alpha and add
281  mfem_x->get_mfem_vec_ref().Add(alpha_f, mfem_tmp);
282 
283  mfem_tmp.Destroy();
284 }
285 
286 void
288 {
289 
290  MFEM_VERIFY(system_oper, "System matrix or operator not initialized");
291  MFEM_VERIFY(executor, "executor is not initialized");
292  MFEM_VERIFY(y.Size() == x.Size(),
293  "Mismatching sizes for rhs and solution");
294 
295  using vec = gko::matrix::Dense<double>;
296  if (!iterative_mode)
297  {
298  y = 0.0;
299  }
300 
301  // Create x and y vectors in Ginkgo's format. Wrap MFEM's data directly,
302  // on CPU or GPU.
303  bool on_device = false;
304  if (executor != executor->get_master())
305  {
306  on_device = true;
307  }
308  std::unique_ptr<vec> gko_x;
309  std::unique_ptr<vec> gko_y;
310 
311  // If we do not have an OperatorWrapper for the system operator or
312  // preconditioner, or have an inner solver using VectorWrappers (as
313  // for IR), then directly create Ginkgo vectors from MFEM's data.
314  if (!needs_wrapped_vecs)
315  {
316  gko_x = vec::create(executor, gko::dim<2>(x.Size(), 1),
318  x.Size(), const_cast<double *>(
319  x.Read(on_device))), 1);
320  gko_y = vec::create(executor, gko::dim<2>(y.Size(), 1),
322  y.Size(),
323  y.ReadWrite(on_device)), 1);
324  }
325  else // We have at least one wrapped MFEM operator; need wrapped vectors
326  {
327  gko_x = std::unique_ptr<vec>(
328  new VectorWrapper(executor, x.Size(),
329  const_cast<Vector *>(&x), false));
330  gko_y = std::unique_ptr<vec>(
331  new VectorWrapper(executor, y.Size(), &y,
332  false));
333  }
334 
335  // Create the logger object to log some data from the solvers to confirm
336  // convergence.
337  initialize_ginkgo_log(gko::lend(gko_x));
338 
339  MFEM_VERIFY(convergence_logger, "convergence logger not initialized" );
340  if (print_level==1)
341  {
342  MFEM_VERIFY(residual_logger, "residual logger not initialized" );
343  solver->clear_loggers(); // Clear any loggers from previous Mult() calls
344  solver->add_logger(residual_logger);
345  }
346 
347  // Add the convergence logger object to the combined factory to retrieve the
348  // solver and other data
349  combined_factory->clear_loggers();
351 
352  // Finally, apply the solver to x and get the solution in y.
353  solver->apply(gko::lend(gko_x), gko::lend(gko_y));
354 
355  // Get the number of iterations taken to converge to the solution.
356  final_iter = convergence_logger->get_num_iterations();
357 
358  // Some residual norm and convergence print outs.
359  double final_res_norm = 0.0;
360 
361  // The convergence_logger object contains the residual vector after the
362  // solver has returned. use this vector to compute the residual norm of the
363  // solution. Get the residual norm from the logger. As the convergence logger
364  // returns a `linop`, it is necessary to convert it to a Dense matrix.
365  // Additionally, if the logger is logging on the gpu, it is necessary to copy
366  // the data to the host and hence the `residual_norm_d_master`
367  auto residual_norm = convergence_logger->get_residual_norm();
368  auto imp_residual_norm = convergence_logger->get_implicit_sq_resnorm();
369 
371  {
372  auto imp_residual_norm_d =
373  gko::as<gko::matrix::Dense<double>>(imp_residual_norm);
374  auto imp_residual_norm_d_master =
375  gko::matrix::Dense<double>::create(executor->get_master(),
376  gko::dim<2> {1, 1});
377  imp_residual_norm_d_master->copy_from(imp_residual_norm_d);
378 
379  final_res_norm = imp_residual_norm_d_master->at(0,0);
380  }
381  else
382  {
383  auto residual_norm_d =
384  gko::as<gko::matrix::Dense<double>>(residual_norm);
385  auto residual_norm_d_master =
386  gko::matrix::Dense<double>::create(executor->get_master(),
387  gko::dim<2> {1, 1});
388  residual_norm_d_master->copy_from(residual_norm_d);
389 
390  final_res_norm = residual_norm_d_master->at(0,0);
391  }
392 
393  converged = 0;
394  if (convergence_logger->has_converged())
395  {
396  converged = 1;
397  }
398 
399  if (print_level == 1)
400  {
401  residual_logger->write();
402  }
403  if (converged == 0)
404  {
405  mfem::err << "No convergence!" << '\n';
406  }
407  if (print_level >=2 && converged==1 )
408  {
409  mfem::out << "Converged in " << final_iter <<
410  " iterations with final residual norm "
411  << final_res_norm << '\n';
412  }
413 }
414 
416 {
417 
418  if (system_oper)
419  {
420  // If the solver currently needs VectorWrappers, but not due to a
421  // "sub-operator" (preconditioner or inner solver), then it's
422  // because the current system_oper needs them. Reset the property
423  // to false in case the new op is a SparseMatrix.
424  if (needs_wrapped_vecs == true && sub_op_needs_wrapped_vecs == false)
425  {
426  needs_wrapped_vecs = false;
427  }
428  // Reset the pointer
429  system_oper.reset();
430  // Reset the solver generated for the previous operator
431  solver.reset();
432  }
433 
434  // Check for SparseMatrix:
435  SparseMatrix *op_mat = const_cast<SparseMatrix*>(
436  dynamic_cast<const SparseMatrix*>(&op));
437  if (op_mat != NULL)
438  {
439  // Needs to be a square matrix
440  MFEM_VERIFY(op_mat->Height() == op_mat->Width(),
441  "System matrix is not square");
442 
443  bool on_device = false;
444  if (executor != executor->get_master())
445  {
446  on_device = true;
447  }
448 
449  using mtx = gko::matrix::Csr<double, int>;
450  const int nnz = op_mat->GetMemoryData().Capacity();
451  system_oper = mtx::create(
452  executor, gko::dim<2>(op_mat->Height(), op_mat->Width()),
454  nnz,
455  op_mat->ReadWriteData(on_device)),
457  nnz,
458  op_mat->ReadWriteJ(on_device)),
459  gko_array<int>::view(executor, op_mat->Height() + 1,
460  op_mat->ReadWriteI(on_device)));
461 
462  }
463  else
464  {
465  needs_wrapped_vecs = true;
466  system_oper = std::shared_ptr<OperatorWrapper>(
467  new OperatorWrapper(executor, op.Height(), &op));
468  }
469 
470  // Generate the solver from the solver using the system matrix or operator.
471  solver = solver_gen->generate(system_oper);
472 }
473 
474 /* ---------------------- CGSolver ------------------------ */
476  : EnableGinkgoSolver(exec, true)
477 {
478  using cg = gko::solver::Cg<double>;
479  this->solver_gen =
480  cg::build().with_criteria(this->combined_factory).on(this->executor);
481 }
482 
484  const GinkgoPreconditioner &preconditioner)
485  : EnableGinkgoSolver(exec, true)
486 {
487  using cg = gko::solver::Cg<double>;
488  // Check for a previously-generated preconditioner (for a specific matrix)
489  if (preconditioner.HasGeneratedPreconditioner())
490  {
491  this->solver_gen = cg::build()
492  .with_criteria(this->combined_factory)
493  .with_generated_preconditioner(
494  preconditioner.GetGeneratedPreconditioner())
495  .on(this->executor);
496  if (dynamic_cast<const OperatorWrapper*>(preconditioner.
497  GetGeneratedPreconditioner().get()))
498  {
499  this->sub_op_needs_wrapped_vecs = true;
500  this->needs_wrapped_vecs = true;
501  }
502  }
503  else // Pass a preconditioner factory (will use same matrix as the solver)
504  {
505  this->solver_gen = cg::build()
506  .with_criteria(this->combined_factory)
507  .with_preconditioner(preconditioner.GetFactory())
508  .on(this->executor);
509  }
510 }
511 
512 
513 /* ---------------------- BICGSTABSolver ------------------------ */
515  : EnableGinkgoSolver(exec, true)
516 {
517  using bicgstab = gko::solver::Bicgstab<double>;
518  this->solver_gen = bicgstab::build()
519  .with_criteria(this->combined_factory)
520  .on(this->executor);
521 }
522 
524  const GinkgoPreconditioner &preconditioner)
525  : EnableGinkgoSolver(exec, true)
526 {
527  using bicgstab = gko::solver::Bicgstab<double>;
528  if (preconditioner.HasGeneratedPreconditioner())
529  {
530  this->solver_gen = bicgstab::build()
531  .with_criteria(this->combined_factory)
532  .with_generated_preconditioner(
533  preconditioner.GetGeneratedPreconditioner())
534  .on(this->executor);
535  if (dynamic_cast<const OperatorWrapper*>(preconditioner.
536  GetGeneratedPreconditioner().get()))
537  {
538  this->sub_op_needs_wrapped_vecs = true;
539  this->needs_wrapped_vecs = true;
540  }
541  }
542  else
543  {
544  this->solver_gen = bicgstab::build()
545  .with_criteria(this->combined_factory)
546  .with_preconditioner(preconditioner.GetFactory())
547  .on(this->executor);
548  }
549 }
550 
551 
552 /* ---------------------- CGSSolver ------------------------ */
554  : EnableGinkgoSolver(exec, true)
555 {
556  using cgs = gko::solver::Cgs<double>;
557  this->solver_gen =
558  cgs::build().with_criteria(this->combined_factory).on(this->executor);
559 }
560 
562  const GinkgoPreconditioner &preconditioner)
563  : EnableGinkgoSolver(exec, true)
564 {
565  using cgs = gko::solver::Cgs<double>;
566  if (preconditioner.HasGeneratedPreconditioner())
567  {
568  this->solver_gen = cgs::build()
569  .with_criteria(this->combined_factory)
570  .with_generated_preconditioner(
571  preconditioner.GetGeneratedPreconditioner())
572  .on(this->executor);
573  if (dynamic_cast<const OperatorWrapper*>(preconditioner.
574  GetGeneratedPreconditioner().get()))
575  {
576  this->sub_op_needs_wrapped_vecs = true;
577  this->needs_wrapped_vecs = true;
578  }
579  }
580  else
581  {
582  this->solver_gen = cgs::build()
583  .with_criteria(this->combined_factory)
584  .with_preconditioner(preconditioner.GetFactory())
585  .on(this->executor);
586  }
587 }
588 
589 
590 /* ---------------------- FCGSolver ------------------------ */
592  : EnableGinkgoSolver(exec, true)
593 {
594  using fcg = gko::solver::Fcg<double>;
595  this->solver_gen =
596  fcg::build().with_criteria(this->combined_factory).on(this->executor);
597 }
598 
600  const GinkgoPreconditioner &preconditioner)
601  : EnableGinkgoSolver(exec, true)
602 {
603  using fcg = gko::solver::Fcg<double>;
604  if (preconditioner.HasGeneratedPreconditioner())
605  {
606  this->solver_gen = fcg::build()
607  .with_criteria(this->combined_factory)
608  .with_generated_preconditioner(
609  preconditioner.GetGeneratedPreconditioner())
610  .on(this->executor);
611  if (dynamic_cast<const OperatorWrapper*>(preconditioner.
612  GetGeneratedPreconditioner().get()))
613  {
614  this->sub_op_needs_wrapped_vecs = true;
615  this->needs_wrapped_vecs = true;
616  }
617  }
618  else
619  {
620  this->solver_gen = fcg::build()
621  .with_criteria(this->combined_factory)
622  .with_preconditioner(preconditioner.GetFactory())
623  .on(this->executor);
624  }
625 }
626 
627 
628 /* ---------------------- GMRESSolver ------------------------ */
630  : EnableGinkgoSolver(exec, false),
631  m{dim}
632 {
633  using gmres = gko::solver::Gmres<double>;
634  if (this->m == 0) // Don't set a dimension, but let Ginkgo use its default
635  {
636  this->solver_gen = gmres::build()
637  .with_criteria(this->combined_factory)
638  .on(this->executor);
639  }
640  else
641  {
642  this->solver_gen = gmres::build()
643  .with_krylov_dim(static_cast<unsigned long>(m))
644  .with_criteria(this->combined_factory)
645  .on(this->executor);
646  }
647 }
648 
650  const GinkgoPreconditioner &preconditioner, int dim)
651  : EnableGinkgoSolver(exec, false),
652  m{dim}
653 {
654  using gmres = gko::solver::Gmres<double>;
655  // Check for a previously-generated preconditioner (for a specific matrix)
656  if (this->m == 0) // Don't set a dimension, but let Ginkgo use its default
657  {
658  if (preconditioner.HasGeneratedPreconditioner())
659  {
660  this->solver_gen = gmres::build()
661  .with_criteria(this->combined_factory)
662  .with_generated_preconditioner(
663  preconditioner.GetGeneratedPreconditioner())
664  .on(this->executor);
665  if (dynamic_cast<const OperatorWrapper*>(preconditioner.
666  GetGeneratedPreconditioner().get()))
667  {
668  this->sub_op_needs_wrapped_vecs = true;
669  this->needs_wrapped_vecs = true;
670  }
671  }
672  else
673  {
674  this->solver_gen = gmres::build()
675  .with_criteria(this->combined_factory)
676  .with_preconditioner(preconditioner.GetFactory())
677  .on(this->executor);
678  }
679  }
680  else
681  {
682  if (preconditioner.HasGeneratedPreconditioner())
683  {
684  this->solver_gen = gmres::build()
685  .with_krylov_dim(static_cast<unsigned long>(m))
686  .with_criteria(this->combined_factory)
687  .with_generated_preconditioner(
688  preconditioner.GetGeneratedPreconditioner())
689  .on(this->executor);
690  if (dynamic_cast<const OperatorWrapper*>(preconditioner.
691  GetGeneratedPreconditioner().get()))
692  {
693  this->sub_op_needs_wrapped_vecs = true;
694  this->needs_wrapped_vecs = true;
695  }
696  }
697  else
698  {
699  this->solver_gen = gmres::build()
700  .with_krylov_dim(static_cast<unsigned long>(m))
701  .with_criteria(this->combined_factory)
702  .with_preconditioner(preconditioner.GetFactory())
703  .on(this->executor);
704  }
705  }
706 }
707 
709 {
710  m = dim;
711  using gmres_type = gko::solver::Gmres<double>;
712  gko::as<gmres_type::Factory>(solver_gen)->get_parameters().krylov_dim = m;
713  if (solver)
714  {
715  gko::as<gmres_type>(solver)->set_krylov_dim(static_cast<unsigned long>(m));
716  }
717 }
718 
719 /* ---------------------- CBGMRESSolver ------------------------ */
721  storage_precision prec)
722  : EnableGinkgoSolver(exec, false),
723  m{dim}
724 {
725  using gmres = gko::solver::CbGmres<double>;
726  if (this->m == 0) // Don't set a dimension, but let Ginkgo use its default
727  {
728  this->solver_gen = gmres::build()
729  .with_criteria(this->combined_factory)
730  .with_storage_precision(prec)
731  .on(this->executor);
732  }
733  else
734  {
735  this->solver_gen = gmres::build()
736  .with_krylov_dim(static_cast<unsigned long>(m))
737  .with_criteria(this->combined_factory)
738  .with_storage_precision(prec)
739  .on(this->executor);
740  }
741 }
742 
744  const GinkgoPreconditioner &preconditioner,
745  int dim, storage_precision prec)
746  : EnableGinkgoSolver(exec, false),
747  m{dim}
748 {
749  using gmres = gko::solver::CbGmres<double>;
750  // Check for a previously-generated preconditioner (for a specific matrix)
751  if (this->m == 0) // Don't set a dimension, but let Ginkgo use its default
752  {
753  if (preconditioner.HasGeneratedPreconditioner())
754  {
755  this->solver_gen = gmres::build()
756  .with_criteria(this->combined_factory)
757  .with_storage_precision(prec)
758  .with_generated_preconditioner(
759  preconditioner.GetGeneratedPreconditioner())
760  .on(this->executor);
761  if (dynamic_cast<const OperatorWrapper*>(preconditioner.
762  GetGeneratedPreconditioner().get()))
763  {
764  this->sub_op_needs_wrapped_vecs = true;
765  this->needs_wrapped_vecs = true;
766  }
767  }
768  else
769  {
770  this->solver_gen = gmres::build()
771  .with_criteria(this->combined_factory)
772  .with_storage_precision(prec)
773  .with_preconditioner(preconditioner.GetFactory())
774  .on(this->executor);
775  }
776  }
777  else
778  {
779  if (preconditioner.HasGeneratedPreconditioner())
780  {
781  this->solver_gen = gmres::build()
782  .with_krylov_dim(static_cast<unsigned long>(m))
783  .with_criteria(this->combined_factory)
784  .with_storage_precision(prec)
785  .with_generated_preconditioner(
786  preconditioner.GetGeneratedPreconditioner())
787  .on(this->executor);
788  if (dynamic_cast<const OperatorWrapper*>(preconditioner.
789  GetGeneratedPreconditioner().get()))
790  {
791  this->sub_op_needs_wrapped_vecs = true;
792  this->needs_wrapped_vecs = true;
793  }
794  }
795  else
796  {
797  this->solver_gen = gmres::build()
798  .with_krylov_dim(static_cast<unsigned long>(m))
799  .with_criteria(this->combined_factory)
800  .with_storage_precision(prec)
801  .with_preconditioner(preconditioner.GetFactory())
802  .on(this->executor);
803  }
804  }
805 }
806 
808 {
809  m = dim;
810  using gmres_type = gko::solver::CbGmres<double>;
811  gko::as<gmres_type::Factory>(solver_gen)->get_parameters().krylov_dim = m;
812  if (solver)
813  {
814  gko::as<gmres_type>(solver)->set_krylov_dim(static_cast<unsigned long>(m));
815  }
816 }
817 
818 /* ---------------------- IRSolver ------------------------ */
820  : EnableGinkgoSolver(exec, false)
821 {
822  using ir = gko::solver::Ir<double>;
823  this->solver_gen =
824  ir::build().with_criteria(this->combined_factory).on(this->executor);
825 }
826 
828  const GinkgoIterativeSolver &inner_solver)
829  : EnableGinkgoSolver(exec, false)
830 {
831  using ir = gko::solver::Ir<double>;
832  this->solver_gen = ir::build()
833  .with_criteria(this->combined_factory)
834  .with_solver(inner_solver.GetFactory())
835  .on(this->executor);
836  if (inner_solver.UsesVectorWrappers())
837  {
838  this->sub_op_needs_wrapped_vecs = true;
839  this->needs_wrapped_vecs = true;
840  }
841 }
842 
843 /* --------------------------------------------------------------- */
844 /* ---------------------- Preconditioners ------------------------ */
846  GinkgoExecutor &exec)
847  : Solver()
848 {
849  has_generated_precond = false;
850  executor = exec.GetExecutor();
851 }
852 
853 void
855 {
856 
857  MFEM_VERIFY(generated_precond, "Preconditioner not initialized");
858  MFEM_VERIFY(executor, "executor is not initialized");
859 
860  using vec = gko::matrix::Dense<double>;
861  if (!iterative_mode)
862  {
863  y = 0.0;
864  }
865 
866  // Create x and y vectors in Ginkgo's format. Wrap MFEM's data directly,
867  // on CPU or GPU.
868  bool on_device = false;
869  if (executor != executor->get_master())
870  {
871  on_device = true;
872  }
873  auto gko_x = vec::create(executor, gko::dim<2>(x.Size(), 1),
875  x.Size(), const_cast<double *>(
876  x.Read(on_device))), 1);
877  auto gko_y = vec::create(executor, gko::dim<2>(y.Size(), 1),
879  y.Size(),
880  y.ReadWrite(on_device)), 1);
881  generated_precond.get()->apply(gko::lend(gko_x), gko::lend(gko_y));
882 }
883 
885 {
886 
888  {
889  generated_precond.reset();
890  has_generated_precond = false;
891  }
892 
893  // Only accept SparseMatrix for this type.
894  SparseMatrix *op_mat = const_cast<SparseMatrix*>(
895  dynamic_cast<const SparseMatrix*>(&op));
896  MFEM_VERIFY(op_mat != NULL,
897  "GinkgoPreconditioner::SetOperator : not a SparseMatrix!");
898 
899  bool on_device = false;
900  if (executor != executor->get_master())
901  {
902  on_device = true;
903  }
904 
905  using mtx = gko::matrix::Csr<double, int>;
906  const int nnz = op_mat->GetMemoryData().Capacity();
907  auto gko_matrix = mtx::create(
908  executor, gko::dim<2>(op_mat->Height(), op_mat->Width()),
910  nnz,
911  op_mat->ReadWriteData(on_device)),
913  nnz,
914  op_mat->ReadWriteJ(on_device)),
915  gko_array<int>::view(executor, op_mat->Height() + 1,
916  op_mat->ReadWriteI(on_device)));
917 
918  generated_precond = precond_gen->generate(gko::give(gko_matrix));
919  has_generated_precond = true;
920 }
921 
922 
923 /* ---------------------- JacobiPreconditioner ------------------------ */
925  GinkgoExecutor &exec,
926  const std::string &storage_opt,
927  const double accuracy,
928  const int max_block_size
929 )
930  : GinkgoPreconditioner(exec)
931 {
932 
933  if (storage_opt == "auto")
934  {
935  precond_gen = gko::preconditioner::Jacobi<double, int>::build()
936  .with_storage_optimization(
937  gko::precision_reduction::autodetect())
938  .with_accuracy(accuracy)
939  .with_max_block_size(static_cast<unsigned int>(max_block_size))
940  .on(executor);
941  }
942  else
943  {
944  precond_gen = gko::preconditioner::Jacobi<double, int>::build()
945  .with_storage_optimization(
946  gko::precision_reduction(0, 0))
947  .with_accuracy(accuracy)
948  .with_max_block_size(static_cast<unsigned int>(max_block_size))
949  .on(executor);
950  }
951 
952 }
953 
954 /* ---------------------- Ilu/IluIsaiPreconditioner ------------------------ */
956  GinkgoExecutor &exec,
957  const std::string &factorization_type,
958  const int sweeps,
959  const bool skip_sort
960 )
961  : GinkgoPreconditioner(exec)
962 {
963  if (factorization_type == "exact")
964  {
965  using ilu_fact_type = gko::factorization::Ilu<double, int>;
966  std::shared_ptr<ilu_fact_type::Factory> fact_factory =
967  ilu_fact_type::build()
968  .with_skip_sorting(skip_sort)
969  .on(executor);
970  precond_gen = gko::preconditioner::Ilu<>::build()
971  .with_factorization_factory(fact_factory)
972  .on(executor);
973  }
974  else
975  {
976  using ilu_fact_type = gko::factorization::ParIlu<double, int>;
977  std::shared_ptr<ilu_fact_type::Factory> fact_factory =
978  ilu_fact_type::build()
979  .with_iterations(static_cast<unsigned long>(sweeps))
980  .with_skip_sorting(skip_sort)
981  .on(executor);
982  precond_gen = gko::preconditioner::Ilu<>::build()
983  .with_factorization_factory(fact_factory)
984  .on(executor);
985  }
986 
987 }
988 
990  GinkgoExecutor &exec,
991  const std::string &factorization_type,
992  const int sweeps,
993  const int sparsity_power,
994  const bool skip_sort
995 )
996  : GinkgoPreconditioner(exec)
997 {
998  using l_solver_type = gko::preconditioner::LowerIsai<>;
999  using u_solver_type = gko::preconditioner::UpperIsai<>;
1000 
1001  std::shared_ptr<l_solver_type::Factory> l_solver_factory =
1002  l_solver_type::build()
1003  .with_sparsity_power(sparsity_power)
1004  .on(executor);
1005  std::shared_ptr<u_solver_type::Factory> u_solver_factory =
1006  u_solver_type::build()
1007  .with_sparsity_power(sparsity_power)
1008  .on(executor);
1009 
1010 
1011 
1012  if (factorization_type == "exact")
1013  {
1014  using ilu_fact_type = gko::factorization::Ilu<double, int>;
1015  std::shared_ptr<ilu_fact_type::Factory> fact_factory =
1016  ilu_fact_type::build()
1017  .with_skip_sorting(skip_sort)
1018  .on(executor);
1019  precond_gen = gko::preconditioner::Ilu<l_solver_type,
1020  u_solver_type>::build()
1021  .with_factorization_factory(fact_factory)
1022  .with_l_solver_factory(l_solver_factory)
1023  .with_u_solver_factory(u_solver_factory)
1024  .on(executor);
1025 
1026  }
1027  else
1028  {
1029  using ilu_fact_type = gko::factorization::ParIlu<double, int>;
1030  std::shared_ptr<ilu_fact_type::Factory> fact_factory =
1031  ilu_fact_type::build()
1032  .with_iterations(static_cast<unsigned long>(sweeps))
1033  .with_skip_sorting(skip_sort)
1034  .on(executor);
1035  precond_gen = gko::preconditioner::Ilu<l_solver_type,
1036  u_solver_type>::build()
1037  .with_factorization_factory(fact_factory)
1038  .with_l_solver_factory(l_solver_factory)
1039  .with_u_solver_factory(u_solver_factory)
1040  .on(executor);
1041  }
1042 }
1043 
1044 
1045 /* ---------------------- Ic/IcIsaiPreconditioner ------------------------ */
1047  GinkgoExecutor &exec,
1048  const std::string &factorization_type,
1049  const int sweeps,
1050  const bool skip_sort
1051 )
1052  : GinkgoPreconditioner(exec)
1053 {
1054 
1055  if (factorization_type == "exact")
1056  {
1057  using ic_fact_type = gko::factorization::Ic<double, int>;
1058  std::shared_ptr<ic_fact_type::Factory> fact_factory =
1059  ic_fact_type::build()
1060  .with_both_factors(false)
1061  .with_skip_sorting(skip_sort)
1062  .on(executor);
1063  precond_gen = gko::preconditioner::Ic<>::build()
1064  .with_factorization_factory(fact_factory)
1065  .on(executor);
1066  }
1067  else
1068  {
1069  using ic_fact_type = gko::factorization::ParIc<double, int>;
1070  std::shared_ptr<ic_fact_type::Factory> fact_factory =
1071  ic_fact_type::build()
1072  .with_both_factors(false)
1073  .with_iterations(static_cast<unsigned long>(sweeps))
1074  .with_skip_sorting(skip_sort)
1075  .on(executor);
1076  precond_gen = gko::preconditioner::Ic<>::build()
1077  .with_factorization_factory(fact_factory)
1078  .on(executor);
1079  }
1080 }
1081 
1083  GinkgoExecutor &exec,
1084  const std::string &factorization_type,
1085  const int sweeps,
1086  const int sparsity_power,
1087  const bool skip_sort
1088 )
1089  : GinkgoPreconditioner(exec)
1090 {
1091 
1092  using l_solver_type = gko::preconditioner::LowerIsai<>;
1093  std::shared_ptr<l_solver_type::Factory> l_solver_factory =
1094  l_solver_type::build()
1095  .with_sparsity_power(sparsity_power)
1096  .on(executor);
1097  if (factorization_type == "exact")
1098  {
1099  using ic_fact_type = gko::factorization::Ic<double, int>;
1100  std::shared_ptr<ic_fact_type::Factory> fact_factory =
1101  ic_fact_type::build()
1102  .with_both_factors(false)
1103  .with_skip_sorting(skip_sort)
1104  .on(executor);
1105  precond_gen = gko::preconditioner::Ic<l_solver_type>::build()
1106  .with_factorization_factory(fact_factory)
1107  .with_l_solver_factory(l_solver_factory)
1108  .on(executor);
1109  }
1110  else
1111  {
1112  using ic_fact_type = gko::factorization::ParIc<double, int>;
1113  std::shared_ptr<ic_fact_type::Factory> fact_factory =
1114  ic_fact_type::build()
1115  .with_both_factors(false)
1116  .with_iterations(static_cast<unsigned long>(sweeps))
1117  .with_skip_sorting(skip_sort)
1118  .on(executor);
1119  precond_gen = gko::preconditioner::Ic<l_solver_type>::build()
1120  .with_factorization_factory(fact_factory)
1121  .with_l_solver_factory(l_solver_factory)
1122  .on(executor);
1123  }
1124 }
1125 
1126 /* ---------------------- MFEMPreconditioner ------------------------ */
1128  GinkgoExecutor &exec,
1129  const Solver &mfem_precond
1130 )
1131  : GinkgoPreconditioner(exec)
1132 {
1133  generated_precond = std::shared_ptr<OperatorWrapper>(
1135  mfem_precond.Height(), &mfem_precond));
1136  has_generated_precond = true;
1137 }
1138 
1139 } // namespace Ginkgo
1140 
1141 } // namespace mfem
1142 
1143 #endif // MFEM_USE_GINKGO
std::shared_ptr< gko::log::Convergence<> > convergence_logger
Definition: ginkgo.hpp:726
std::shared_ptr< ResidualLogger<> > residual_logger
Definition: ginkgo.hpp:732
virtual void SetOperator(const Operator &op)
Definition: ginkgo.cpp:415
std::shared_ptr< gko::stop::ResidualNorm<>::Factory > abs_criterion
Definition: ginkgo.hpp:706
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
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:720
const std::shared_ptr< gko::LinOp > GetGeneratedPreconditioner() const
Definition: ginkgo.hpp:562
virtual void UseDevice(bool use_dev) const
Enable execution of Vector operations using the mfem::Device.
Definition: vector.hpp:117
std::shared_ptr< gko::LinOpFactory > precond_gen
Definition: ginkgo.hpp:574
CGSSolver(GinkgoExecutor &exec)
Definition: ginkgo.cpp:553
int Width() const
Get the width (size of input) of the Operator. Synonym with NumCols().
Definition: operator.hpp:72
Biwise-OR of all HIP backends.
Definition: device.hpp:90
MFEMPreconditioner(GinkgoExecutor &exec, const Solver &mfem_precond)
Definition: ginkgo.cpp:1127
JacobiPreconditioner(GinkgoExecutor &exec, const std::string &storage_opt="none", const double accuracy=1.e-1, const int max_block_size=32)
Definition: ginkgo.cpp:924
int Size() const
Returns the size of the vector.
Definition: vector.hpp:199
std::shared_ptr< gko::Executor > executor
Definition: ginkgo.hpp:594
virtual const double * Read(bool on_dev=true) const
Shortcut for mfem::Read(vec.GetMemory(), vec.Size(), on_dev).
Definition: vector.hpp:448
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:708
bool iterative_mode
If true, use the second argument of Mult() as an initial guess.
Definition: operator.hpp:686
Memory< double > & GetMemoryData()
Definition: sparsemat.hpp:269
virtual void Mult(const Vector &x, Vector &y) const
Definition: ginkgo.cpp:287
BICGSTABSolver(GinkgoExecutor &exec)
Definition: ginkgo.cpp:514
bool HasGeneratedPreconditioner() const
Definition: ginkgo.hpp:571
std::shared_ptr< gko::stop::ImplicitResidualNorm<>::Factory > imp_rel_criterion
Definition: ginkgo.hpp:713
IluPreconditioner(GinkgoExecutor &exec, const std::string &factorization_type="exact", const int sweeps=0, const bool skip_sort=false)
Definition: ginkgo.cpp:955
Data type sparse matrix.
Definition: sparsemat.hpp:50
double b
Definition: lissajous.cpp:42
int * ReadWriteI(bool on_dev=true)
Definition: sparsemat.hpp:243
const std::shared_ptr< gko::LinOpFactory > GetFactory() const
Definition: ginkgo.hpp:621
std::shared_ptr< gko::Executor > executor
Definition: ginkgo.hpp:745
int Capacity() const
Return the size of the allocated memory.
CGSolver(GinkgoExecutor &exec)
Definition: ginkgo.cpp:475
virtual void Mult(const Vector &x, Vector &y) const
Definition: ginkgo.cpp:854
Biwise-OR of all CUDA backends.
Definition: device.hpp:88
GMRESSolver(GinkgoExecutor &exec, int dim=0)
Definition: ginkgo.cpp:629
OpenMP CPU Executor.
Definition: ginkgo.hpp:463
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
std::shared_ptr< gko::Executor > GetExecutor() const
Definition: ginkgo.hpp:489
std::shared_ptr< gko::stop::ResidualNorm<>::Factory > rel_criterion
Definition: ginkgo.hpp:699
int Height() const
Get the height (size of output) of the Operator. Synonym with NumRows().
Definition: operator.hpp:66
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:819
const Vector & get_mfem_vec_const_ref() const
Definition: ginkgo.hpp:122
IcPreconditioner(GinkgoExecutor &exec, const std::string &factorization_type="exact", const int sweeps=0, const bool skip_sort=false)
Definition: ginkgo.cpp:1046
std::shared_ptr< gko::stop::Combined::Factory > combined_factory
Definition: ginkgo.hpp:738
virtual double * ReadWrite(bool on_dev=true)
Shortcut for mfem::ReadWrite(vec.GetMemory(), vec.Size(), on_dev).
Definition: vector.hpp:464
std::shared_ptr< gko::LinOpFactory > solver_gen
Definition: ginkgo.hpp:687
const std::shared_ptr< gko::LinOpFactory > GetFactory() const
Definition: ginkgo.hpp:553
FCGSolver(GinkgoExecutor &exec)
Definition: ginkgo.cpp:591
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:1082
void Destroy()
Destroy a vector.
Definition: vector.hpp:589
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:989
GinkgoPreconditioner(GinkgoExecutor &exec)
Definition: ginkgo.cpp:845
std::shared_ptr< gko::LinOp > generated_precond
Definition: ginkgo.hpp:587
const double alpha
Definition: ex15.cpp:369
double * ReadWriteData(bool on_dev=true)
Definition: sparsemat.hpp:275
int * ReadWriteJ(bool on_dev=true)
Definition: sparsemat.hpp:259
Vector data type.
Definition: vector.hpp:60
std::shared_ptr< gko::stop::ImplicitResidualNorm<>::Factory > imp_abs_criterion
Definition: ginkgo.hpp:720
Base class for solvers.
Definition: operator.hpp:682
Reference CPU Executor.
Definition: ginkgo.hpp:461
virtual void SetOperator(const Operator &op)
Definition: ginkgo.cpp:884
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
gko::Array< T > gko_array
Definition: ginkgo.hpp:43
void apply_impl(const gko::LinOp *b, gko::LinOp *x) const override
Definition: ginkgo.cpp:202
std::shared_ptr< gko::LinOp > solver
Definition: ginkgo.hpp:692