MFEM  v4.5.1
Finite element discretization library
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Pages
hypre.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_MPI
15 
16 #include "linalg.hpp"
17 #include "../fem/fem.hpp"
18 #include "../general/forall.hpp"
19 
20 #include <fstream>
21 #include <iomanip>
22 #include <cmath>
23 #include <cstdlib>
24 
25 using namespace std;
26 
27 namespace mfem
28 {
29 
30 Hypre::Hypre()
31 {
32 #if MFEM_HYPRE_VERSION >= 21900
33  // Initializing hypre
34  HYPRE_Init();
35 #endif
36 
37  // Global hypre options that we set by default
38  SetDefaultOptions();
39 }
40 
41 void Hypre::Finalize()
42 {
43  Hypre &hypre = Instance();
44  if (!hypre.finalized)
45  {
46 #if MFEM_HYPRE_VERSION >= 21900
47  HYPRE_Finalize();
48 #endif
49  hypre.finalized = true;
50  }
51 }
52 
53 void Hypre::SetDefaultOptions()
54 {
55  // Global hypre options, see
56  // https://hypre.readthedocs.io/en/latest/solvers-boomeramg.html#gpu-supported-options
57 
58 #if MFEM_HYPRE_VERSION >= 22100
59 #ifdef HYPRE_USING_CUDA
60  // Use hypre's SpGEMM instead of cuSPARSE for performance reasons
61  HYPRE_SetSpGemmUseCusparse(0);
62 #elif defined(HYPRE_USING_HIP)
63  // Use rocSPARSE instead of hypre's SpGEMM for performance reasons (default)
64  // HYPRE_SetSpGemmUseCusparse(1);
65 #endif
66 #endif
67 
68  // The following options are hypre's defaults as of hypre-2.24
69 
70  // Allocate hypre objects in GPU memory (default)
71  // HYPRE_SetMemoryLocation(HYPRE_MEMORY_DEVICE);
72 
73  // Where to execute when using UVM (default)
74  // HYPRE_SetExecutionPolicy(HYPRE_EXEC_DEVICE);
75 
76  // Use GPU-based random number generator (default)
77  // HYPRE_SetUseGpuRand(1);
78 
79  // The following options are to be used with UMPIRE memory pools
80 
81  // Set Umpire names for device and UVM memory pools. If names are set by
82  // calling these functions, hypre doesn't own the pool and just uses it.If
83  // these functions are not called, hypre will allocate and own the pool
84  // (provided it is configured with --with-umpire).
85  // HYPRE_SetUmpireDevicePoolName("HYPRE_DEVICE_POOL");
86  // HYPRE_SetUmpireUMPoolName("HYPRE_UVM_POOL");
87 }
88 
89 
90 template<typename TargetT, typename SourceT>
91 static TargetT *DuplicateAs(const SourceT *array, int size,
92  bool cplusplus = true)
93 {
94  TargetT *target_array = cplusplus ? (TargetT*) Memory<TargetT>(size)
95  /* */ : mfem_hypre_TAlloc_host(TargetT, size);
96  for (int i = 0; i < size; i++)
97  {
98  target_array[i] = array[i];
99  }
100  return target_array;
101 }
102 
103 
104 /// Return true if the @a src Memory can be used with the MemoryClass @a mc.
105 /** If this function returns true then src.{Read,Write,ReadWrite} can be called
106  safely with the MemoryClass @a mc. */
107 template <typename T>
109 {
110  MemoryType src_h_mt = src.GetHostMemoryType();
111  MemoryType src_d_mt = src.GetDeviceMemoryType();
112  if (src_d_mt == MemoryType::DEFAULT)
113  {
114  src_d_mt = MemoryManager::GetDualMemoryType(src_h_mt);
115  }
116  return (MemoryClassContainsType(mc, src_h_mt) ||
117  MemoryClassContainsType(mc, src_d_mt));
118 }
119 
120 
121 inline void HypreParVector::_SetDataAndSize_()
122 {
123  hypre_Vector *x_loc = hypre_ParVectorLocalVector(x);
124 #if !defined(HYPRE_USING_GPU)
125  SetDataAndSize(hypre_VectorData(x_loc),
126  internal::to_int(hypre_VectorSize(x_loc)));
127 #else
128  size = internal::to_int(hypre_VectorSize(x_loc));
129  MemoryType mt = (hypre_VectorMemoryLocation(x_loc) == HYPRE_MEMORY_HOST
130  ? MemoryType::HOST : GetHypreMemoryType());
131  if (hypre_VectorData(x_loc) != NULL)
132  {
133  data.Wrap(hypre_VectorData(x_loc), size, mt, false);
134  }
135  else
136  {
137  data.Reset();
138  }
139 #endif
140 }
141 
142 HypreParVector::HypreParVector(MPI_Comm comm, HYPRE_BigInt glob_size,
143  HYPRE_BigInt *col) : Vector()
144 {
145  x = hypre_ParVectorCreate(comm,glob_size,col);
146  hypre_ParVectorInitialize(x);
147 #if MFEM_HYPRE_VERSION <= 22200
148  hypre_ParVectorSetPartitioningOwner(x,0);
149 #endif
150  // The data will be destroyed by hypre (this is the default)
151  hypre_ParVectorSetDataOwner(x,1);
152  hypre_SeqVectorSetDataOwner(hypre_ParVectorLocalVector(x),1);
153  _SetDataAndSize_();
154  own_ParVector = 1;
155 }
156 
158  double *data_, HYPRE_BigInt *col,
159  bool is_device_ptr)
160  : Vector()
161 {
162  x = hypre_ParVectorCreate(comm,glob_size,col);
163  hypre_ParVectorSetDataOwner(x,1); // owns the seq vector
164  hypre_Vector *x_loc = hypre_ParVectorLocalVector(x);
165  hypre_SeqVectorSetDataOwner(x_loc,0);
166 #if MFEM_HYPRE_VERSION <= 22200
167  hypre_ParVectorSetPartitioningOwner(x,0);
168 #endif
169  double tmp = 0.0;
170  hypre_VectorData(x_loc) = &tmp;
171 #ifdef HYPRE_USING_GPU
172  hypre_VectorMemoryLocation(x_loc) =
173  is_device_ptr ? HYPRE_MEMORY_DEVICE : HYPRE_MEMORY_HOST;
174 #else
175  (void)is_device_ptr;
176 #endif
177  // If hypre_ParVectorLocalVector(x) and &tmp are non-NULL,
178  // hypre_ParVectorInitialize(x) does not allocate memory!
179  hypre_ParVectorInitialize(x);
180  // Set the internal data array to the one passed in
181  hypre_VectorData(x_loc) = data_;
182  _SetDataAndSize_();
183  own_ParVector = 1;
184 }
185 
186 // Call the move constructor on the "compatible" temp vector
188  y.CreateCompatibleVector())
189 {
190  // Deep copy the local data
191  hypre_SeqVectorCopy(hypre_ParVectorLocalVector(y.x),
192  hypre_ParVectorLocalVector(x));
193 }
194 
196 {
197  own_ParVector = 0;
198  *this = std::move(y);
199 }
200 
202  int transpose) : Vector()
203 {
204  if (!transpose)
205  {
206  x = hypre_ParVectorInDomainOf(const_cast<HypreParMatrix&>(A));
207  }
208  else
209  {
210  x = hypre_ParVectorInRangeOf(const_cast<HypreParMatrix&>(A));
211  }
212  _SetDataAndSize_();
213  own_ParVector = 1;
214 }
215 
217 {
218  x = (hypre_ParVector *) y;
219  _SetDataAndSize_();
220  own_ParVector = 0;
221 }
222 
224 {
225  x = hypre_ParVectorCreate(pfes->GetComm(), pfes->GlobalTrueVSize(),
226  pfes->GetTrueDofOffsets());
227  hypre_ParVectorInitialize(x);
228 #if MFEM_HYPRE_VERSION <= 22200
229  hypre_ParVectorSetPartitioningOwner(x,0);
230 #endif
231  // The data will be destroyed by hypre (this is the default)
232  hypre_ParVectorSetDataOwner(x,1);
233  hypre_SeqVectorSetDataOwner(hypre_ParVectorLocalVector(x),1);
234  _SetDataAndSize_();
235  own_ParVector = 1;
236 }
237 
239 {
240  HypreParVector result;
241  result.x = hypre_ParVectorCreate(x -> comm, x -> global_size,
242  x -> partitioning);
243  hypre_ParVectorInitialize(result.x);
244 #if MFEM_HYPRE_VERSION <= 22200
245  hypre_ParVectorSetPartitioningOwner(result.x,0);
246 #endif
247  hypre_ParVectorSetDataOwner(result.x,1);
248  hypre_SeqVectorSetDataOwner(hypre_ParVectorLocalVector(result.x),1);
249  result._SetDataAndSize_();
250  result.own_ParVector = 1;
251 
252  return result;
253 }
254 
255 void HypreParVector::WrapHypreParVector(hypre_ParVector *y, bool owner)
256 {
257  if (own_ParVector) { hypre_ParVectorDestroy(x); }
258  Destroy();
259  x = y;
260  _SetDataAndSize_();
261  own_ParVector = owner;
262 }
263 
265 {
266  hypre_Vector *hv = hypre_ParVectorToVectorAll(*this);
267  Vector *v = new Vector(hv->data, internal::to_int(hv->size));
268  v->MakeDataOwner();
269  hypre_SeqVectorSetDataOwner(hv,0);
270  hypre_SeqVectorDestroy(hv);
271  return v;
272 }
273 
275 {
277  return *this;
278 }
279 
281 {
282 #ifdef MFEM_DEBUG
283  if (size != y.Size())
284  {
285  mfem_error("HypreParVector::operator=");
286  }
287 #endif
288 
290  return *this;
291 }
292 
294 {
295  // If the argument vector owns its data, then the calling vector will as well
296  WrapHypreParVector(static_cast<hypre_ParVector*>(y), y.own_ParVector);
297  // Either way the argument vector will no longer own its data
298  y.own_ParVector = 0;
299  y.x = nullptr;
300  y.data.Reset();
301  y.size = 0;
302  return *this;
303 }
304 
305 void HypreParVector::SetData(double *data_)
306 {
307  hypre_VectorData(hypre_ParVectorLocalVector(x)) = data_;
308  Vector::SetData(data_);
309 }
310 
312 {
313  hypre_Vector *x_loc = hypre_ParVectorLocalVector(x);
314  hypre_VectorData(x_loc) =
315  const_cast<double*>(data.Read(GetHypreMemoryClass(), size));
316 #ifdef HYPRE_USING_GPU
317  hypre_VectorMemoryLocation(x_loc) = HYPRE_MEMORY_DEVICE;
318 #endif
319 }
320 
322 {
323  hypre_Vector *x_loc = hypre_ParVectorLocalVector(x);
324  hypre_VectorData(x_loc) = data.ReadWrite(GetHypreMemoryClass(), size);
325 #ifdef HYPRE_USING_GPU
326  hypre_VectorMemoryLocation(x_loc) = HYPRE_MEMORY_DEVICE;
327 #endif
328 }
329 
331 {
332  hypre_Vector *x_loc = hypre_ParVectorLocalVector(x);
333  hypre_VectorData(x_loc) = data.Write(GetHypreMemoryClass(), size);
334 #ifdef HYPRE_USING_GPU
335  hypre_VectorMemoryLocation(x_loc) = HYPRE_MEMORY_DEVICE;
336 #endif
337 }
338 
340 {
341  MFEM_ASSERT(CanShallowCopy(mem, GetHypreMemoryClass()), "");
342  MFEM_ASSERT(mem.Capacity() >= size, "");
343 
344  data.Delete();
345  hypre_Vector *x_loc = hypre_ParVectorLocalVector(x);
346  hypre_VectorData(x_loc) =
347  const_cast<double*>(mem.Read(GetHypreMemoryClass(), size));
348 #ifdef HYPRE_USING_GPU
349  hypre_VectorMemoryLocation(x_loc) = HYPRE_MEMORY_DEVICE;
350 #endif
351  data.MakeAlias(mem, 0, size);
352 }
353 
355 {
356  MFEM_ASSERT(CanShallowCopy(mem, GetHypreMemoryClass()), "");
357  MFEM_ASSERT(mem.Capacity() >= size, "");
358 
359  data.Delete();
360  hypre_Vector *x_loc = hypre_ParVectorLocalVector(x);
361  hypre_VectorData(x_loc) = mem.ReadWrite(GetHypreMemoryClass(), size);
362 #ifdef HYPRE_USING_GPU
363  hypre_VectorMemoryLocation(x_loc) = HYPRE_MEMORY_DEVICE;
364 #endif
365  data.MakeAlias(mem, 0, size);
366 }
367 
369 {
370  MFEM_ASSERT(CanShallowCopy(mem, GetHypreMemoryClass()), "");
371  MFEM_ASSERT(mem.Capacity() >= size, "");
372 
373  data.Delete();
374  hypre_Vector *x_loc = hypre_ParVectorLocalVector(x);
375  hypre_VectorData(x_loc) = mem.Write(GetHypreMemoryClass(), size);
376 #ifdef HYPRE_USING_GPU
377  hypre_VectorMemoryLocation(x_loc) = HYPRE_MEMORY_DEVICE;
378 #endif
379  data.MakeAlias(mem, 0, size);
380 }
381 
382 HYPRE_Int HypreParVector::Randomize(HYPRE_Int seed)
383 {
384  return hypre_ParVectorSetRandomValues(x,seed);
385 }
386 
387 void HypreParVector::Print(const char *fname) const
388 {
389  hypre_ParVectorPrint(x,fname);
390 }
391 
392 void HypreParVector::Read(MPI_Comm comm, const char *fname)
393 {
394  if (own_ParVector)
395  {
396  hypre_ParVectorDestroy(x);
397  }
398  data.Delete();
399  x = hypre_ParVectorRead(comm, fname);
400  own_ParVector = true;
401  _SetDataAndSize_();
402 }
403 
405 {
406  if (own_ParVector)
407  {
408  hypre_ParVectorDestroy(x);
409  }
410 }
411 
412 
414 {
415  return hypre_ParVectorInnerProd(*x, *y);
416 }
417 
419 {
420  return hypre_ParVectorInnerProd(x, y);
421 }
422 
423 
424 double ParNormlp(const Vector &vec, double p, MPI_Comm comm)
425 {
426  double norm = 0.0;
427  if (p == 1.0)
428  {
429  double loc_norm = vec.Norml1();
430  MPI_Allreduce(&loc_norm, &norm, 1, MPI_DOUBLE, MPI_SUM, comm);
431  }
432  if (p == 2.0)
433  {
434  double loc_norm = vec*vec;
435  MPI_Allreduce(&loc_norm, &norm, 1, MPI_DOUBLE, MPI_SUM, comm);
436  norm = sqrt(norm);
437  }
438  if (p < infinity())
439  {
440  double sum = 0.0;
441  for (int i = 0; i < vec.Size(); i++)
442  {
443  sum += pow(fabs(vec(i)), p);
444  }
445  MPI_Allreduce(&sum, &norm, 1, MPI_DOUBLE, MPI_SUM, comm);
446  norm = pow(norm, 1.0/p);
447  }
448  else
449  {
450  double loc_norm = vec.Normlinf();
451  MPI_Allreduce(&loc_norm, &norm, 1, MPI_DOUBLE, MPI_MAX, comm);
452  }
453  return norm;
454 }
455 
456 /** @brief Shallow or deep copy @a src to @a dst with the goal to make the
457  array @a src accessible through @a dst with the MemoryClass @a dst_mc. If
458  one of the host/device MemoryType%s of @a src is contained in @a dst_mc,
459  then a shallow copy will be used and @a dst will simply be an alias of
460  @a src. Otherwise, @a dst will be properly allocated and @a src will be deep
461  copied to @a dst. */
462 /** If @a dst_owner is set to true and shallow copy is being used, then @a dst
463  will not be an alias of @a src; instead, @a src is copied to @a dst and all
464  ownership flags of @a src are reset.
465 
466  In both cases (deep or shallow copy), when @a dst is no longer needed,
467  dst.Delete() must be called to ensure all associated memory allocations are
468  freed.
469 
470  The input contents of @a dst, if any, is not used and it is overwritten by
471  this function. In particular, @a dst should be empty or deleted before
472  calling this function. */
473 template <typename T>
474 void CopyMemory(Memory<T> &src, Memory<T> &dst, MemoryClass dst_mc,
475  bool dst_owner)
476 {
477  if (CanShallowCopy(src, dst_mc))
478  {
479  // shallow copy
480  if (!dst_owner)
481  {
482  src.Read(dst_mc, src.Capacity()); // Registers src if on host only
483  dst.MakeAlias(src, 0, src.Capacity());
484  }
485  else
486  {
487  dst = src;
488  src.ClearOwnerFlags();
489  }
490  }
491  else
492  {
493  // deep copy
494  dst.New(src.Capacity(), GetMemoryType(dst_mc));
495  dst.CopyFrom(src, src.Capacity());
496  }
497 }
498 
499 /** @brief Deep copy and convert @a src to @a dst with the goal to make the
500  array @a src accessible through @a dst with the MemoryClass @a dst_mc and
501  convert it from type SrcT to type DstT. */
502 /** When @a dst is no longer needed, dst.Delete() must be called to ensure all
503  associated memory allocations are freed.
504 
505  The input contents of @a dst, if any, is not used and it is overwritten by
506  this function. In particular, @a dst should be empty or deleted before
507  calling this function. */
508 template <typename SrcT, typename DstT>
510 {
511  auto capacity = src.Capacity();
512  dst.New(capacity, GetMemoryType(dst_mc));
513  // Perform the copy using the configured mfem Device
514  auto src_p = mfem::Read(src, capacity);
515  auto dst_p = mfem::Write(dst, capacity);
516  MFEM_FORALL(i, capacity, dst_p[i] = src_p[i];);
517 }
518 
519 
520 void HypreParMatrix::Init()
521 {
522  A = NULL;
523  X = Y = NULL;
524  auxX.Reset(); auxY.Reset();
525  diagOwner = offdOwner = colMapOwner = -1;
526  ParCSROwner = 1;
527  mem_diag.I.Reset();
528  mem_diag.J.Reset();
529  mem_diag.data.Reset();
530  mem_offd.I.Reset();
531  mem_offd.J.Reset();
532  mem_offd.data.Reset();
533 }
534 
535 void HypreParMatrix::Read(MemoryClass mc) const
536 {
537  hypre_CSRMatrix *diag = hypre_ParCSRMatrixDiag(A);
538  hypre_CSRMatrix *offd = hypre_ParCSRMatrixOffd(A);
539  const int num_rows = NumRows();
540  const int diag_nnz = internal::to_int(diag->num_nonzeros);
541  const int offd_nnz = internal::to_int(offd->num_nonzeros);
542  diag->i = const_cast<HYPRE_Int*>(mem_diag.I.Read(mc, num_rows+1));
543  diag->j = const_cast<HYPRE_Int*>(mem_diag.J.Read(mc, diag_nnz));
544  diag->data = const_cast<double*>(mem_diag.data.Read(mc, diag_nnz));
545  offd->i = const_cast<HYPRE_Int*>(mem_offd.I.Read(mc, num_rows+1));
546  offd->j = const_cast<HYPRE_Int*>(mem_offd.J.Read(mc, offd_nnz));
547  offd->data = const_cast<double*>(mem_offd.data.Read(mc, offd_nnz));
548 #if MFEM_HYPRE_VERSION >= 21800
549  decltype(diag->memory_location) ml =
550  (mc != GetHypreMemoryClass() ? HYPRE_MEMORY_HOST : HYPRE_MEMORY_DEVICE);
551  diag->memory_location = ml;
552  offd->memory_location = ml;
553 #endif
554 }
555 
556 void HypreParMatrix::ReadWrite(MemoryClass mc)
557 {
558  hypre_CSRMatrix *diag = hypre_ParCSRMatrixDiag(A);
559  hypre_CSRMatrix *offd = hypre_ParCSRMatrixOffd(A);
560  const int num_rows = NumRows();
561  const int diag_nnz = internal::to_int(diag->num_nonzeros);
562  const int offd_nnz = internal::to_int(offd->num_nonzeros);
563  diag->i = mem_diag.I.ReadWrite(mc, num_rows+1);
564  diag->j = mem_diag.J.ReadWrite(mc, diag_nnz);
565  diag->data = mem_diag.data.ReadWrite(mc, diag_nnz);
566  offd->i = mem_offd.I.ReadWrite(mc, num_rows+1);
567  offd->j = mem_offd.J.ReadWrite(mc, offd_nnz);
568  offd->data = mem_offd.data.ReadWrite(mc, offd_nnz);
569 #if MFEM_HYPRE_VERSION >= 21800
570  decltype(diag->memory_location) ml =
571  (mc != GetHypreMemoryClass() ? HYPRE_MEMORY_HOST : HYPRE_MEMORY_DEVICE);
572  diag->memory_location = ml;
573  offd->memory_location = ml;
574 #endif
575 }
576 
577 void HypreParMatrix::Write(MemoryClass mc, bool set_diag, bool set_offd)
578 {
579  hypre_CSRMatrix *diag = hypre_ParCSRMatrixDiag(A);
580  hypre_CSRMatrix *offd = hypre_ParCSRMatrixOffd(A);
581  if (set_diag)
582  {
583  diag->i = mem_diag.I.Write(mc, mem_diag.I.Capacity());
584  diag->j = mem_diag.J.Write(mc, mem_diag.J.Capacity());
585  diag->data = mem_diag.data.Write(mc, mem_diag.data.Capacity());
586  }
587  if (set_offd)
588  {
589  offd->i = mem_offd.I.Write(mc, mem_offd.I.Capacity());
590  offd->j = mem_offd.J.Write(mc, mem_offd.J.Capacity());
591  offd->data = mem_offd.data.Write(mc, mem_offd.data.Capacity());
592  }
593 #if MFEM_HYPRE_VERSION >= 21800
594  decltype(diag->memory_location) ml =
595  (mc != GetHypreMemoryClass() ? HYPRE_MEMORY_HOST : HYPRE_MEMORY_DEVICE);
596  if (set_diag) { diag->memory_location = ml; }
597  if (set_offd) { offd->memory_location = ml; }
598 #endif
599 }
600 
602 {
603  Init();
604  height = width = 0;
605 }
606 
607 void HypreParMatrix::WrapHypreParCSRMatrix(hypre_ParCSRMatrix *a, bool owner)
608 {
609  Destroy();
610  Init();
611  A = a;
612  ParCSROwner = owner;
613  height = GetNumRows();
614  width = GetNumCols();
615 #if MFEM_HYPRE_VERSION >= 21800
616  MemoryType diag_mt = (A->diag->memory_location == HYPRE_MEMORY_HOST
618  MemoryType offd_mt = (A->offd->memory_location == HYPRE_MEMORY_HOST
620 #else
621  const MemoryType diag_mt = MemoryType::HOST;
622  const MemoryType offd_mt = MemoryType::HOST;
623 #endif
624  diagOwner = HypreCsrToMem(A->diag, diag_mt, false, mem_diag);
625  offdOwner = HypreCsrToMem(A->offd, offd_mt, false, mem_offd);
626  HypreRead();
627 }
628 
629 signed char HypreParMatrix::CopyCSR(SparseMatrix *csr,
630  MemoryIJData &mem_csr,
631  hypre_CSRMatrix *hypre_csr,
632  bool mem_owner)
633 {
634  const MemoryClass hypre_mc = GetHypreMemoryClass();
635 #ifndef HYPRE_BIGINT
636  // code for the case HYPRE_Int == int
637  CopyMemory(csr->GetMemoryI(), mem_csr.I, hypre_mc, mem_owner);
638  CopyMemory(csr->GetMemoryJ(), mem_csr.J, hypre_mc, mem_owner);
639 #else
640  // code for the case HYPRE_Int == long long int
641  CopyConvertMemory(csr->GetMemoryI(), hypre_mc, mem_csr.I);
642  CopyConvertMemory(csr->GetMemoryJ(), hypre_mc, mem_csr.J);
643 #endif
644  CopyMemory(csr->GetMemoryData(), mem_csr.data, hypre_mc, mem_owner);
645 
646  const int num_rows = csr->Height();
647  const int nnz = csr->NumNonZeroElems();
648  hypre_csr->i = const_cast<HYPRE_Int*>(mem_csr.I.Read(hypre_mc, num_rows+1));
649  hypre_csr->j = const_cast<HYPRE_Int*>(mem_csr.J.Read(hypre_mc, nnz));
650  hypre_csr->data = const_cast<double*>(mem_csr.data.Read(hypre_mc, nnz));
651 
652  MFEM_ASSERT(mem_csr.I.OwnsHostPtr() == mem_csr.J.OwnsHostPtr(),
653  "invalid state: host ownership for I and J differ!");
654  return (mem_csr.I.OwnsHostPtr() ? 1 : 0) +
655  (mem_csr.data.OwnsHostPtr() ? 2 : 0);
656 }
657 
658 signed char HypreParMatrix::CopyBoolCSR(Table *bool_csr,
659  MemoryIJData &mem_csr,
660  hypre_CSRMatrix *hypre_csr)
661 {
662  const MemoryClass hypre_mc = GetHypreMemoryClass();
663 #ifndef HYPRE_BIGINT
664  // code for the case HYPRE_Int == int
665  CopyMemory(bool_csr->GetIMemory(), mem_csr.I, hypre_mc, false);
666  CopyMemory(bool_csr->GetJMemory(), mem_csr.J, hypre_mc, false);
667 #else
668  // code for the case HYPRE_Int == long long int
669  CopyConvertMemory(bool_csr->GetIMemory(), hypre_mc, mem_csr.I);
670  CopyConvertMemory(bool_csr->GetJMemory(), hypre_mc, mem_csr.J);
671 #endif
672  const int num_rows = bool_csr->Size();
673  const int nnz = bool_csr->Size_of_connections();
674  mem_csr.data.New(nnz, GetHypreMemoryType());
675  double *data = mfem::HostWrite(mem_csr.data, nnz);
676  for (int i = 0; i < nnz; i++)
677  {
678  data[i] = 1.0;
679  }
680  hypre_csr->i = const_cast<HYPRE_Int*>(mem_csr.I.Read(hypre_mc, num_rows+1));
681  hypre_csr->j = const_cast<HYPRE_Int*>(mem_csr.J.Read(hypre_mc, nnz));
682  hypre_csr->data = const_cast<double*>(mem_csr.data.Read(hypre_mc, nnz));
683 
684  MFEM_ASSERT(mem_csr.I.OwnsHostPtr() == mem_csr.J.OwnsHostPtr(),
685  "invalid state: host ownership for I and J differ!");
686  return (mem_csr.I.OwnsHostPtr() ? 1 : 0) +
687  (mem_csr.data.OwnsHostPtr() ? 2 : 0);
688 }
689 
690 // Copy the j array of a MemoryIJData object to the given dst_J array,
691 // converting the indices from HYPRE_Int to int.
692 #ifdef HYPRE_BIGINT
693 static void CopyCSR_J(const int nnz, const MemoryIJData &mem_csr,
694  Memory<int> &dst_J)
695 {
696  // Perform the copy using the configured mfem Device
697  auto src_p = mfem::Read(mem_csr.J, nnz);
698  auto dst_p = mfem::Write(dst_J, nnz);
699  MFEM_FORALL(i, nnz, dst_p[i] = src_p[i];);
700 }
701 #endif
702 
703 // static method
704 signed char HypreParMatrix::HypreCsrToMem(hypre_CSRMatrix *h_mat,
705  MemoryType h_mat_mt,
706  bool own_ija,
707  MemoryIJData &mem)
708 {
709  const int nr1 = internal::to_int(h_mat->num_rows) + 1;
710  const int nnz = internal::to_int(h_mat->num_nonzeros);
711  mem.I.Wrap(h_mat->i, nr1, h_mat_mt, own_ija);
712  mem.J.Wrap(h_mat->j, nnz, h_mat_mt, own_ija);
713  mem.data.Wrap(h_mat->data, nnz, h_mat_mt, own_ija);
714  const MemoryClass hypre_mc = GetHypreMemoryClass();
715  if (!CanShallowCopy(mem.I, hypre_mc))
716  {
717  const MemoryType hypre_mt = GetHypreMemoryType();
718  MemoryIJData h_mem;
719  h_mem.I.New(nr1, hypre_mt);
720  h_mem.I.CopyFrom(mem.I, nr1);
721  mem.I.Delete();
722  h_mem.J.New(nnz, hypre_mt);
723  h_mem.J.CopyFrom(mem.J, nnz);
724  mem.J.Delete();
725  h_mem.data.New(nnz, hypre_mt);
726  h_mem.data.CopyFrom(mem.data, nnz);
727  mem.data.Delete();
728  mem = h_mem;
729  if (!own_ija)
730  {
731  // FIXME: Even if own_ija == false, it does not necessarily mean we
732  // need to delete h_mat->{i,j,data} even if h_mat->owns_data == true.
733 
734  // h_mat owns i; owns j,data if h_mat->owns_data
735 #if MFEM_HYPRE_VERSION < 21400
736  hypre_TFree(h_mat->i);
737 #elif MFEM_HYPRE_VERSION < 21800
738  hypre_TFree(h_mat->i, HYPRE_MEMORY_SHARED);
739 #else
740  hypre_TFree(h_mat->i, h_mat->memory_location);
741 #endif
742  if (h_mat->owns_data)
743  {
744 #if MFEM_HYPRE_VERSION < 21400
745  hypre_TFree(h_mat->j);
746  hypre_TFree(h_mat->data);
747 #elif MFEM_HYPRE_VERSION < 21800
748  hypre_TFree(h_mat->j, HYPRE_MEMORY_SHARED);
749  hypre_TFree(h_mat->data, HYPRE_MEMORY_SHARED);
750 #else
751  hypre_TFree(h_mat->j, h_mat->memory_location);
752  hypre_TFree(h_mat->data, h_mat->memory_location);
753 #endif
754  }
755  }
756  h_mat->i = mem.I.ReadWrite(hypre_mc, nr1);
757  h_mat->j = mem.J.ReadWrite(hypre_mc, nnz);
758  h_mat->data = mem.data.ReadWrite(hypre_mc, nnz);
759  h_mat->owns_data = 0;
760 #if MFEM_HYPRE_VERSION >= 21800
761  h_mat->memory_location = HYPRE_MEMORY_DEVICE;
762 #endif
763  return 3;
764  }
765  return own_ija ? 3 : (h_mat_mt == GetHypreMemoryType() ? -2 : -1);
766 }
767 
768 // Square block-diagonal constructor (4 arguments, v1)
770  HYPRE_BigInt *row_starts, SparseMatrix *diag)
771  : Operator(diag->Height(), diag->Width())
772 {
773  Init();
774  A = hypre_ParCSRMatrixCreate(comm, glob_size, glob_size, row_starts,
775  row_starts, 0, diag->NumNonZeroElems(), 0);
776  hypre_ParCSRMatrixSetDataOwner(A,1);
777 #if MFEM_HYPRE_VERSION <= 22200
778  hypre_ParCSRMatrixSetRowStartsOwner(A,0);
779  hypre_ParCSRMatrixSetColStartsOwner(A,0);
780 #endif
781 
782  hypre_CSRMatrixSetDataOwner(A->diag,0);
783  diagOwner = CopyCSR(diag, mem_diag, A->diag, false);
784  hypre_CSRMatrixSetRownnz(A->diag);
785 
786  hypre_CSRMatrixSetDataOwner(A->offd,1);
787  hypre_CSRMatrixI(A->offd) = mfem_hypre_CTAlloc(HYPRE_Int, diag->Height()+1);
788  offdOwner = HypreCsrToMem(A->offd, GetHypreMemoryType(), false, mem_offd);
789 
790  /* Don't need to call these, since they allocate memory only
791  if it was not already allocated */
792  // hypre_CSRMatrixInitialize(A->diag);
793  // hypre_ParCSRMatrixInitialize(A);
794 
795  hypre_ParCSRMatrixSetNumNonzeros(A);
796 
797  /* Make sure that the first entry in each row is the diagonal one. */
798  HypreReadWrite();
799  hypre_CSRMatrixReorder(hypre_ParCSRMatrixDiag(A));
800 #ifdef HYPRE_BIGINT
802  {
803  CopyCSR_J(A->diag->num_nonzeros, mem_diag, diag->GetMemoryJ());
804  }
805 #endif
806 
807  hypre_MatvecCommPkgCreate(A);
808 }
809 
810 // Rectangular block-diagonal constructor (6 arguments, v1)
812  HYPRE_BigInt global_num_rows,
813  HYPRE_BigInt global_num_cols,
814  HYPRE_BigInt *row_starts,
815  HYPRE_BigInt *col_starts,
816  SparseMatrix *diag)
817  : Operator(diag->Height(), diag->Width())
818 {
819  Init();
820  A = hypre_ParCSRMatrixCreate(comm, global_num_rows, global_num_cols,
821  row_starts, col_starts,
822  0, diag->NumNonZeroElems(), 0);
823  hypre_ParCSRMatrixSetDataOwner(A,1);
824 #if MFEM_HYPRE_VERSION <= 22200
825  hypre_ParCSRMatrixSetRowStartsOwner(A,0);
826  hypre_ParCSRMatrixSetColStartsOwner(A,0);
827 #endif
828 
829  hypre_CSRMatrixSetDataOwner(A->diag,0);
830  diagOwner = CopyCSR(diag, mem_diag, A->diag, false);
831  hypre_CSRMatrixSetRownnz(A->diag);
832 
833  hypre_CSRMatrixSetDataOwner(A->offd,1);
834  hypre_CSRMatrixI(A->offd) = mfem_hypre_CTAlloc(HYPRE_Int, diag->Height()+1);
835  offdOwner = HypreCsrToMem(A->offd, GetHypreMemoryType(), false, mem_offd);
836 
837  hypre_ParCSRMatrixSetNumNonzeros(A);
838 
839  /* Make sure that the first entry in each row is the diagonal one. */
840  if (row_starts == col_starts)
841  {
842  HypreReadWrite();
843  hypre_CSRMatrixReorder(hypre_ParCSRMatrixDiag(A));
844 #ifdef HYPRE_BIGINT
846  {
847  CopyCSR_J(A->diag->num_nonzeros, mem_diag, diag->GetMemoryJ());
848  }
849 #endif
850  }
851 
852  hypre_MatvecCommPkgCreate(A);
853 }
854 
855 // General rectangular constructor with diagonal and off-diagonal (8+1
856 // arguments)
858  HYPRE_BigInt global_num_rows,
859  HYPRE_BigInt global_num_cols,
860  HYPRE_BigInt *row_starts,
861  HYPRE_BigInt *col_starts,
862  SparseMatrix *diag, SparseMatrix *offd,
863  HYPRE_BigInt *cmap,
864  bool own_diag_offd)
865  : Operator(diag->Height(), diag->Width())
866 {
867  Init();
868  A = hypre_ParCSRMatrixCreate(comm, global_num_rows, global_num_cols,
869  row_starts, col_starts,
870  offd->Width(), diag->NumNonZeroElems(),
871  offd->NumNonZeroElems());
872  hypre_ParCSRMatrixSetDataOwner(A,1);
873 #if MFEM_HYPRE_VERSION <= 22200
874  hypre_ParCSRMatrixSetRowStartsOwner(A,0);
875  hypre_ParCSRMatrixSetColStartsOwner(A,0);
876 #endif
877 
878  hypre_CSRMatrixSetDataOwner(A->diag,0);
879  diagOwner = CopyCSR(diag, mem_diag, A->diag, own_diag_offd);
880  if (own_diag_offd) { delete diag; }
881  hypre_CSRMatrixSetRownnz(A->diag);
882 
883  hypre_CSRMatrixSetDataOwner(A->offd,0);
884  offdOwner = CopyCSR(offd, mem_offd, A->offd, own_diag_offd);
885  if (own_diag_offd) { delete offd; }
886  hypre_CSRMatrixSetRownnz(A->offd);
887 
888  hypre_ParCSRMatrixColMapOffd(A) = cmap;
889  // Prevent hypre from destroying A->col_map_offd
890  colMapOwner = 0;
891 
892  hypre_ParCSRMatrixSetNumNonzeros(A);
893 
894  /* Make sure that the first entry in each row is the diagonal one. */
895  if (row_starts == col_starts)
896  {
897  HypreReadWrite();
898  hypre_CSRMatrixReorder(hypre_ParCSRMatrixDiag(A));
899 #ifdef HYPRE_BIGINT
901  {
902  CopyCSR_J(A->diag->num_nonzeros, mem_diag, diag->GetMemoryJ());
903  }
904 #endif
905  }
906 
907  hypre_MatvecCommPkgCreate(A);
908 }
909 
910 // General rectangular constructor with diagonal and off-diagonal (13+1
911 // arguments)
913  MPI_Comm comm,
914  HYPRE_BigInt global_num_rows, HYPRE_BigInt global_num_cols,
915  HYPRE_BigInt *row_starts, HYPRE_BigInt *col_starts,
916  HYPRE_Int *diag_i, HYPRE_Int *diag_j, double *diag_data,
917  HYPRE_Int *offd_i, HYPRE_Int *offd_j, double *offd_data,
918  HYPRE_Int offd_num_cols, HYPRE_BigInt *offd_col_map,
919  bool hypre_arrays)
920 {
921  Init();
922  A = hypre_ParCSRMatrixCreate(comm, global_num_rows, global_num_cols,
923  row_starts, col_starts, offd_num_cols, 0, 0);
924  hypre_ParCSRMatrixSetDataOwner(A,1);
925 #if MFEM_HYPRE_VERSION <= 22200
926  hypre_ParCSRMatrixSetRowStartsOwner(A,0);
927  hypre_ParCSRMatrixSetColStartsOwner(A,0);
928 #endif
929 
930  HYPRE_Int local_num_rows = hypre_CSRMatrixNumRows(A->diag);
931 
932  hypre_CSRMatrixSetDataOwner(A->diag, hypre_arrays);
933  hypre_CSRMatrixI(A->diag) = diag_i;
934  hypre_CSRMatrixJ(A->diag) = diag_j;
935  hypre_CSRMatrixData(A->diag) = diag_data;
936  hypre_CSRMatrixNumNonzeros(A->diag) = diag_i[local_num_rows];
937 #ifdef HYPRE_USING_GPU
938  hypre_CSRMatrixMemoryLocation(A->diag) = HYPRE_MEMORY_HOST;
939 #endif
940 
941  hypre_CSRMatrixSetDataOwner(A->offd, hypre_arrays);
942  hypre_CSRMatrixI(A->offd) = offd_i;
943  hypre_CSRMatrixJ(A->offd) = offd_j;
944  hypre_CSRMatrixData(A->offd) = offd_data;
945  hypre_CSRMatrixNumNonzeros(A->offd) = offd_i[local_num_rows];
946 #ifdef HYPRE_USING_GPU
947  hypre_CSRMatrixMemoryLocation(A->offd) = HYPRE_MEMORY_HOST;
948 #endif
949 
950  hypre_ParCSRMatrixColMapOffd(A) = offd_col_map;
951  // Prevent hypre from destroying A->col_map_offd, own A->col_map_offd
952  colMapOwner = hypre_arrays ? -1 : 1;
953 
954  hypre_ParCSRMatrixSetNumNonzeros(A);
955 
956  /* Make sure that the first entry in each row is the diagonal one. */
957  if (row_starts == col_starts)
958  {
959  hypre_CSRMatrixReorder(hypre_ParCSRMatrixDiag(A));
960  }
961 
962  hypre_MatvecCommPkgCreate(A);
963 
964  height = GetNumRows();
965  width = GetNumCols();
966 
967  if (!hypre_arrays)
968  {
969  const MemoryType host_mt = Device::GetHostMemoryType();
970  diagOwner = HypreCsrToMem(A->diag, host_mt, true, mem_diag);
971  offdOwner = HypreCsrToMem(A->offd, host_mt, true, mem_offd);
972  }
973  else
974  {
975  const MemoryType host_mt = MemoryType::HOST;
976  diagOwner = HypreCsrToMem(A->diag, host_mt, false, mem_diag);
977  offdOwner = HypreCsrToMem(A->offd, host_mt, false, mem_offd);
978  }
979  HypreRead();
980 
981  hypre_CSRMatrixSetRownnz(A->diag);
982  hypre_CSRMatrixSetRownnz(A->offd);
983 }
984 
985 // Constructor from a CSR matrix on rank 0 (4 arguments, v2)
987  HYPRE_BigInt *row_starts,
988  HYPRE_BigInt *col_starts,
989  SparseMatrix *sm_a)
990 {
991  MFEM_ASSERT(sm_a != NULL, "invalid input");
992  MFEM_VERIFY(!HYPRE_AssumedPartitionCheck(),
993  "this method can not be used with assumed partition");
994 
995  Init();
996 
997  hypre_CSRMatrix *csr_a;
998  csr_a = hypre_CSRMatrixCreate(sm_a -> Height(), sm_a -> Width(),
999  sm_a -> NumNonZeroElems());
1000 
1001  hypre_CSRMatrixSetDataOwner(csr_a,0);
1002  MemoryIJData mem_a;
1003  CopyCSR(sm_a, mem_a, csr_a, false);
1004  hypre_CSRMatrixSetRownnz(csr_a);
1005 
1006  // NOTE: this call creates a matrix on host even when device support is
1007  // enabled in hypre.
1008  hypre_ParCSRMatrix *new_A =
1009  hypre_CSRMatrixToParCSRMatrix(comm, csr_a, row_starts, col_starts);
1010 
1011  mem_a.I.Delete();
1012  mem_a.J.Delete();
1013  mem_a.data.Delete();
1014 
1015  hypre_CSRMatrixI(csr_a) = NULL;
1016  hypre_CSRMatrixDestroy(csr_a);
1017 
1018  /* Make sure that the first entry in each row is the diagonal one. */
1019  if (row_starts == col_starts)
1020  {
1021  hypre_CSRMatrixReorder(hypre_ParCSRMatrixDiag(new_A));
1022  }
1023 
1024  hypre_MatvecCommPkgCreate(A);
1025 
1026  WrapHypreParCSRMatrix(new_A);
1027 }
1028 
1029 // Boolean, rectangular, block-diagonal constructor (6 arguments, v2)
1031  HYPRE_BigInt global_num_rows,
1032  HYPRE_BigInt global_num_cols,
1033  HYPRE_BigInt *row_starts,
1034  HYPRE_BigInt *col_starts,
1035  Table *diag)
1036 {
1037  Init();
1038  int nnz = diag->Size_of_connections();
1039  A = hypre_ParCSRMatrixCreate(comm, global_num_rows, global_num_cols,
1040  row_starts, col_starts, 0, nnz, 0);
1041  hypre_ParCSRMatrixSetDataOwner(A,1);
1042 #if MFEM_HYPRE_VERSION <= 22200
1043  hypre_ParCSRMatrixSetRowStartsOwner(A,0);
1044  hypre_ParCSRMatrixSetColStartsOwner(A,0);
1045 #endif
1046 
1047  hypre_CSRMatrixSetDataOwner(A->diag,0);
1048  diagOwner = CopyBoolCSR(diag, mem_diag, A->diag);
1049  hypre_CSRMatrixSetRownnz(A->diag);
1050 
1051  hypre_CSRMatrixSetDataOwner(A->offd,1);
1052  hypre_CSRMatrixI(A->offd) = mfem_hypre_CTAlloc(HYPRE_Int, diag->Size()+1);
1053  offdOwner = HypreCsrToMem(A->offd, GetHypreMemoryType(), false, mem_offd);
1054 
1055  hypre_ParCSRMatrixSetNumNonzeros(A);
1056 
1057  /* Make sure that the first entry in each row is the diagonal one. */
1058  if (row_starts == col_starts)
1059  {
1060  HypreReadWrite();
1061  hypre_CSRMatrixReorder(hypre_ParCSRMatrixDiag(A));
1062 #ifdef HYPRE_BIGINT
1063  // No need to sync the J array back to the Table diag.
1064  // CopyCSR_J(A->diag->num_nonzeros, mem_diag, diag->GetJMemory());
1065 #endif
1066  }
1067 
1068  hypre_MatvecCommPkgCreate(A);
1069 
1070  height = GetNumRows();
1071  width = GetNumCols();
1072 }
1073 
1074 // Boolean, general rectangular constructor with diagonal and off-diagonal
1075 // (11 arguments)
1076 HypreParMatrix::HypreParMatrix(MPI_Comm comm, int id, int np,
1077  HYPRE_BigInt *row, HYPRE_BigInt *col,
1078  HYPRE_Int *i_diag, HYPRE_Int *j_diag,
1079  HYPRE_Int *i_offd, HYPRE_Int *j_offd,
1080  HYPRE_BigInt *cmap, HYPRE_Int cmap_size)
1081 {
1082  HYPRE_Int diag_nnz, offd_nnz;
1083 
1084  Init();
1085  if (HYPRE_AssumedPartitionCheck())
1086  {
1087  diag_nnz = i_diag[row[1]-row[0]];
1088  offd_nnz = i_offd[row[1]-row[0]];
1089 
1090  A = hypre_ParCSRMatrixCreate(comm, row[2], col[2], row, col,
1091  cmap_size, diag_nnz, offd_nnz);
1092  }
1093  else
1094  {
1095  diag_nnz = i_diag[row[id+1]-row[id]];
1096  offd_nnz = i_offd[row[id+1]-row[id]];
1097 
1098  A = hypre_ParCSRMatrixCreate(comm, row[np], col[np], row, col,
1099  cmap_size, diag_nnz, offd_nnz);
1100  }
1101 
1102  hypre_ParCSRMatrixSetDataOwner(A,1);
1103 #if MFEM_HYPRE_VERSION <= 22200
1104  hypre_ParCSRMatrixSetRowStartsOwner(A,0);
1105  hypre_ParCSRMatrixSetColStartsOwner(A,0);
1106 #endif
1107 
1108  mem_diag.data.New(diag_nnz);
1109  for (HYPRE_Int i = 0; i < diag_nnz; i++)
1110  {
1111  mem_diag.data[i] = 1.0;
1112  }
1113 
1114  mem_offd.data.New(offd_nnz);
1115  for (HYPRE_Int i = 0; i < offd_nnz; i++)
1116  {
1117  mem_offd.data[i] = 1.0;
1118  }
1119 
1120  hypre_CSRMatrixSetDataOwner(A->diag,0);
1121  hypre_CSRMatrixI(A->diag) = i_diag;
1122  hypre_CSRMatrixJ(A->diag) = j_diag;
1123  hypre_CSRMatrixData(A->diag) = mem_diag.data;
1124 #ifdef HYPRE_USING_GPU
1125  hypre_CSRMatrixMemoryLocation(A->diag) = HYPRE_MEMORY_HOST;
1126 #endif
1127 
1128  hypre_CSRMatrixSetDataOwner(A->offd,0);
1129  hypre_CSRMatrixI(A->offd) = i_offd;
1130  hypre_CSRMatrixJ(A->offd) = j_offd;
1131  hypre_CSRMatrixData(A->offd) = mem_offd.data;
1132 #ifdef HYPRE_USING_GPU
1133  hypre_CSRMatrixMemoryLocation(A->offd) = HYPRE_MEMORY_HOST;
1134 #endif
1135 
1136  hypre_ParCSRMatrixColMapOffd(A) = cmap;
1137  // Prevent hypre from destroying A->col_map_offd, own A->col_map_offd
1138  colMapOwner = 1;
1139 
1140  hypre_ParCSRMatrixSetNumNonzeros(A);
1141 
1142  /* Make sure that the first entry in each row is the diagonal one. */
1143  if (row == col)
1144  {
1145  hypre_CSRMatrixReorder(hypre_ParCSRMatrixDiag(A));
1146  }
1147 
1148  hypre_MatvecCommPkgCreate(A);
1149 
1150  height = GetNumRows();
1151  width = GetNumCols();
1152 
1153  const MemoryType host_mt = Device::GetHostMemoryType();
1154  diagOwner = HypreCsrToMem(A->diag, host_mt, true, mem_diag);
1155  offdOwner = HypreCsrToMem(A->offd, host_mt, true, mem_offd);
1156  HypreRead();
1157 
1158  hypre_CSRMatrixSetRownnz(A->diag);
1159  hypre_CSRMatrixSetRownnz(A->offd);
1160 }
1161 
1162 // General rectangular constructor with diagonal and off-diagonal constructed
1163 // from a CSR matrix that contains both diagonal and off-diagonal blocks
1164 // (9 arguments)
1165 HypreParMatrix::HypreParMatrix(MPI_Comm comm, int nrows,
1166  HYPRE_BigInt glob_nrows,
1167  HYPRE_BigInt glob_ncols,
1168  int *I, HYPRE_BigInt *J,
1169  double *data,
1170  HYPRE_BigInt *rows,
1171  HYPRE_BigInt *cols)
1172 {
1173  Init();
1174 
1175  // Determine partitioning size, and my column start and end
1176  int part_size;
1177  HYPRE_BigInt my_col_start, my_col_end; // my range: [my_col_start, my_col_end)
1178  if (HYPRE_AssumedPartitionCheck())
1179  {
1180  part_size = 2;
1181  my_col_start = cols[0];
1182  my_col_end = cols[1];
1183  }
1184  else
1185  {
1186  int myid;
1187  MPI_Comm_rank(comm, &myid);
1188  MPI_Comm_size(comm, &part_size);
1189  part_size++;
1190  my_col_start = cols[myid];
1191  my_col_end = cols[myid+1];
1192  }
1193 
1194  // Copy in the row and column partitionings
1195  HYPRE_BigInt *row_starts, *col_starts;
1196  if (rows == cols)
1197  {
1198  row_starts = col_starts = mfem_hypre_TAlloc_host(HYPRE_BigInt, part_size);
1199  for (int i = 0; i < part_size; i++)
1200  {
1201  row_starts[i] = rows[i];
1202  }
1203  }
1204  else
1205  {
1206  row_starts = mfem_hypre_TAlloc_host(HYPRE_BigInt, part_size);
1207  col_starts = mfem_hypre_TAlloc_host(HYPRE_BigInt, part_size);
1208  for (int i = 0; i < part_size; i++)
1209  {
1210  row_starts[i] = rows[i];
1211  col_starts[i] = cols[i];
1212  }
1213  }
1214 
1215  // Create a map for the off-diagonal indices - global to local. Count the
1216  // number of diagonal and off-diagonal entries.
1217  HYPRE_Int diag_nnz = 0, offd_nnz = 0, offd_num_cols = 0;
1218  map<HYPRE_BigInt, HYPRE_Int> offd_map;
1219  for (HYPRE_Int j = 0, loc_nnz = I[nrows]; j < loc_nnz; j++)
1220  {
1221  HYPRE_BigInt glob_col = J[j];
1222  if (my_col_start <= glob_col && glob_col < my_col_end)
1223  {
1224  diag_nnz++;
1225  }
1226  else
1227  {
1228  offd_map.insert(pair<const HYPRE_BigInt, HYPRE_Int>(glob_col, -1));
1229  offd_nnz++;
1230  }
1231  }
1232  // count the number of columns in the off-diagonal and set the local indices
1233  for (auto it = offd_map.begin(); it != offd_map.end(); ++it)
1234  {
1235  it->second = offd_num_cols++;
1236  }
1237 
1238  // construct the global ParCSR matrix
1239  A = hypre_ParCSRMatrixCreate(comm, glob_nrows, glob_ncols,
1240  row_starts, col_starts, offd_num_cols,
1241  diag_nnz, offd_nnz);
1242  hypre_ParCSRMatrixInitialize(A);
1243 
1244  diagOwner = HypreCsrToMem(A->diag, GetHypreMemoryType(), false, mem_diag);
1245  offdOwner = HypreCsrToMem(A->offd, GetHypreMemoryType(), false, mem_offd);
1246  HostWrite();
1247 
1248  HYPRE_Int *diag_i, *diag_j, *offd_i, *offd_j;
1249  HYPRE_BigInt *offd_col_map;
1250  double *diag_data, *offd_data;
1251  diag_i = A->diag->i;
1252  diag_j = A->diag->j;
1253  diag_data = A->diag->data;
1254  offd_i = A->offd->i;
1255  offd_j = A->offd->j;
1256  offd_data = A->offd->data;
1257  offd_col_map = A->col_map_offd;
1258 
1259  diag_nnz = offd_nnz = 0;
1260  for (HYPRE_Int i = 0, j = 0; i < nrows; i++)
1261  {
1262  diag_i[i] = diag_nnz;
1263  offd_i[i] = offd_nnz;
1264  for (HYPRE_Int j_end = I[i+1]; j < j_end; j++)
1265  {
1266  HYPRE_BigInt glob_col = J[j];
1267  if (my_col_start <= glob_col && glob_col < my_col_end)
1268  {
1269  diag_j[diag_nnz] = glob_col - my_col_start;
1270  diag_data[diag_nnz] = data[j];
1271  diag_nnz++;
1272  }
1273  else
1274  {
1275  offd_j[offd_nnz] = offd_map[glob_col];
1276  offd_data[offd_nnz] = data[j];
1277  offd_nnz++;
1278  }
1279  }
1280  }
1281  diag_i[nrows] = diag_nnz;
1282  offd_i[nrows] = offd_nnz;
1283  for (auto it = offd_map.begin(); it != offd_map.end(); ++it)
1284  {
1285  offd_col_map[it->second] = it->first;
1286  }
1287 
1288  hypre_ParCSRMatrixSetNumNonzeros(A);
1289  /* Make sure that the first entry in each row is the diagonal one. */
1290  if (row_starts == col_starts)
1291  {
1292  hypre_CSRMatrixReorder(hypre_ParCSRMatrixDiag(A));
1293  }
1294 #if MFEM_HYPRE_VERSION > 22200
1295  mfem_hypre_TFree_host(row_starts);
1296  if (rows != cols)
1297  {
1298  mfem_hypre_TFree_host(col_starts);
1299  }
1300 #endif
1301  hypre_MatvecCommPkgCreate(A);
1302 
1303  height = GetNumRows();
1304  width = GetNumCols();
1305 
1306  HypreRead();
1307 }
1308 
1310 {
1311  hypre_ParCSRMatrix *Ph = static_cast<hypre_ParCSRMatrix *>(P);
1312 
1313  Init();
1314 
1315  // Clone the structure
1316  A = hypre_ParCSRMatrixCompleteClone(Ph);
1317  // Make a deep copy of the data from the source
1318  hypre_ParCSRMatrixCopy(Ph, A, 1);
1319 
1320  height = GetNumRows();
1321  width = GetNumCols();
1322 
1323  CopyRowStarts();
1324  CopyColStarts();
1325 
1326  hypre_ParCSRMatrixSetNumNonzeros(A);
1327 
1328  hypre_MatvecCommPkgCreate(A);
1329 
1330  diagOwner = HypreCsrToMem(A->diag, GetHypreMemoryType(), false, mem_diag);
1331  offdOwner = HypreCsrToMem(A->offd, GetHypreMemoryType(), false, mem_offd);
1332 }
1333 
1335 {
1336  Destroy();
1337  Init();
1338  A = master.A;
1339  ParCSROwner = 0;
1340  height = master.GetNumRows();
1341  width = master.GetNumCols();
1342  mem_diag.I.MakeAlias(master.mem_diag.I, 0, master.mem_diag.I.Capacity());
1343  mem_diag.J.MakeAlias(master.mem_diag.J, 0, master.mem_diag.J.Capacity());
1344  mem_diag.data.MakeAlias(master.mem_diag.data, 0,
1345  master.mem_diag.data.Capacity());
1346  mem_offd.I.MakeAlias(master.mem_offd.I, 0, master.mem_offd.I.Capacity());
1347  mem_offd.J.MakeAlias(master.mem_offd.J, 0, master.mem_offd.J.Capacity());
1348  mem_offd.data.MakeAlias(master.mem_offd.data, 0,
1349  master.mem_offd.data.Capacity());
1350 }
1351 
1352 hypre_ParCSRMatrix* HypreParMatrix::StealData()
1353 {
1354  // Only safe when (diagOwner < 0 && offdOwner < 0 && colMapOwner == -1)
1355  // Otherwise, there may be memory leaks or hypre may destroy arrays allocated
1356  // with operator new.
1357  MFEM_ASSERT(diagOwner < 0 && offdOwner < 0 && colMapOwner == -1, "");
1358  MFEM_ASSERT(diagOwner == offdOwner, "");
1359  MFEM_ASSERT(ParCSROwner, "");
1360  hypre_ParCSRMatrix *R = A;
1361 #ifdef HYPRE_USING_GPU
1362  if (diagOwner == -1) { HostReadWrite(); }
1363  else { HypreReadWrite(); }
1364 #endif
1365  ParCSROwner = false;
1366  Destroy();
1367  Init();
1368  return R;
1369 }
1370 
1371 void HypreParMatrix::SetOwnerFlags(signed char diag, signed char offd,
1372  signed char colmap)
1373 {
1374  diagOwner = diag;
1375  mem_diag.I.SetHostPtrOwner((diag >= 0) && (diag & 1));
1376  mem_diag.I.SetDevicePtrOwner((diag >= 0) && (diag & 1));
1377 
1378  mem_diag.J.SetHostPtrOwner((diag >= 0) && (diag & 1));
1379  mem_diag.J.SetDevicePtrOwner((diag >= 0) && (diag & 1));
1380 
1381  mem_diag.data.SetHostPtrOwner((diag >= 0) && (diag & 2));
1382  mem_diag.data.SetDevicePtrOwner((diag >= 0) && (diag & 2));
1383 
1384  offdOwner = offd;
1385  mem_offd.I.SetHostPtrOwner((offd >= 0) && (offd & 1));
1386  mem_offd.J.SetHostPtrOwner((offd >= 0) && (offd & 1));
1387 
1388  mem_offd.I.SetDevicePtrOwner((offd >= 0) && (offd & 1));
1389  mem_offd.J.SetDevicePtrOwner((offd >= 0) && (offd & 1));
1390 
1391  mem_offd.data.SetHostPtrOwner((offd >= 0) && (offd & 2));
1392  mem_offd.data.SetDevicePtrOwner((offd >= 0) && (offd & 2));
1393  colMapOwner = colmap;
1394 }
1395 
1397 {
1398 #if MFEM_HYPRE_VERSION <= 22200
1399  if (!A || hypre_ParCSRMatrixOwnsRowStarts(A) ||
1400  (hypre_ParCSRMatrixRowStarts(A) == hypre_ParCSRMatrixColStarts(A) &&
1401  hypre_ParCSRMatrixOwnsColStarts(A)))
1402  {
1403  return;
1404  }
1405 
1406  int row_starts_size;
1407  if (HYPRE_AssumedPartitionCheck())
1408  {
1409  row_starts_size = 2;
1410  }
1411  else
1412  {
1413  MPI_Comm_size(hypre_ParCSRMatrixComm(A), &row_starts_size);
1414  row_starts_size++; // num_proc + 1
1415  }
1416 
1417  HYPRE_BigInt *old_row_starts = hypre_ParCSRMatrixRowStarts(A);
1418  HYPRE_BigInt *new_row_starts = mfem_hypre_CTAlloc_host(HYPRE_BigInt,
1419  row_starts_size);
1420  for (int i = 0; i < row_starts_size; i++)
1421  {
1422  new_row_starts[i] = old_row_starts[i];
1423  }
1424 
1425  hypre_ParCSRMatrixRowStarts(A) = new_row_starts;
1426  hypre_ParCSRMatrixOwnsRowStarts(A) = 1;
1427 
1428  if (hypre_ParCSRMatrixColStarts(A) == old_row_starts)
1429  {
1430  hypre_ParCSRMatrixColStarts(A) = new_row_starts;
1431  hypre_ParCSRMatrixOwnsColStarts(A) = 0;
1432  }
1433 #endif
1434 }
1435 
1437 {
1438 #if MFEM_HYPRE_VERSION <= 22200
1439  if (!A || hypre_ParCSRMatrixOwnsColStarts(A) ||
1440  (hypre_ParCSRMatrixRowStarts(A) == hypre_ParCSRMatrixColStarts(A) &&
1441  hypre_ParCSRMatrixOwnsRowStarts(A)))
1442  {
1443  return;
1444  }
1445 
1446  int col_starts_size;
1447  if (HYPRE_AssumedPartitionCheck())
1448  {
1449  col_starts_size = 2;
1450  }
1451  else
1452  {
1453  MPI_Comm_size(hypre_ParCSRMatrixComm(A), &col_starts_size);
1454  col_starts_size++; // num_proc + 1
1455  }
1456 
1457  HYPRE_BigInt *old_col_starts = hypre_ParCSRMatrixColStarts(A);
1458  HYPRE_BigInt *new_col_starts = mfem_hypre_CTAlloc_host(HYPRE_BigInt,
1459  col_starts_size);
1460  for (int i = 0; i < col_starts_size; i++)
1461  {
1462  new_col_starts[i] = old_col_starts[i];
1463  }
1464 
1465  hypre_ParCSRMatrixColStarts(A) = new_col_starts;
1466 
1467  if (hypre_ParCSRMatrixRowStarts(A) == old_col_starts)
1468  {
1469  hypre_ParCSRMatrixRowStarts(A) = new_col_starts;
1470  hypre_ParCSRMatrixOwnsRowStarts(A) = 1;
1471  hypre_ParCSRMatrixOwnsColStarts(A) = 0;
1472  }
1473  else
1474  {
1475  hypre_ParCSRMatrixOwnsColStarts(A) = 1;
1476  }
1477 #endif
1478 }
1479 
1481 {
1482  const int size = Height();
1483  diag.SetSize(size);
1484 #ifdef HYPRE_USING_GPU
1486  {
1487  MFEM_ASSERT(A->diag->memory_location == HYPRE_MEMORY_DEVICE, "");
1488  double *d_diag = diag.Write();
1489  const HYPRE_Int *A_diag_i = A->diag->i;
1490  const double *A_diag_d = A->diag->data;
1491  MFEM_FORALL(i, size, d_diag[i] = A_diag_d[A_diag_i[i]];);
1492  }
1493  else
1494 #endif
1495  {
1496  diag.HostWrite();
1497  HostRead();
1498  for (int j = 0; j < size; j++)
1499  {
1500  diag(j) = A->diag->data[A->diag->i[j]];
1501  MFEM_ASSERT(A->diag->j[A->diag->i[j]] == j,
1502  "the first entry in each row must be the diagonal one");
1503  }
1504  HypreRead();
1505  }
1506 }
1507 
1508 static void MakeSparseMatrixWrapper(int nrows, int ncols,
1509  HYPRE_Int *I, HYPRE_Int *J, double *data,
1510  SparseMatrix &wrapper)
1511 {
1512 #ifndef HYPRE_BIGINT
1513  SparseMatrix tmp(I, J, data, nrows, ncols, false, false, false);
1514 #else
1515  int *mI = Memory<int>(nrows + 1);
1516  for (int i = 0; i <= nrows; i++)
1517  {
1518  mI[i] = internal::to_int(I[i]); // checks for overflow in debug mode
1519  }
1520  const int nnz = mI[nrows];
1521  int *mJ = Memory<int>(nnz);
1522  for (int j = 0; j < nnz; j++)
1523  {
1524  mJ[j] = internal::to_int(J[j]); // checks for overflow in debug mode
1525  }
1526  SparseMatrix tmp(mI, mJ, data, nrows, ncols, true, false, false);
1527 #endif
1528  wrapper.Swap(tmp);
1529 }
1530 
1531 static void MakeWrapper(const hypre_CSRMatrix *mat,
1532  const MemoryIJData &mem,
1533  SparseMatrix &wrapper)
1534 {
1535  const int nrows = internal::to_int(hypre_CSRMatrixNumRows(mat));
1536  const int ncols = internal::to_int(hypre_CSRMatrixNumCols(mat));
1537  const int nnz = internal::to_int(mat->num_nonzeros);
1538  const HYPRE_Int *I = mfem::HostRead(mem.I, nrows + 1);
1539  const HYPRE_Int *J = mfem::HostRead(mem.J, nnz);
1540  const double *data = mfem::HostRead(mem.data, nnz);
1541  MakeSparseMatrixWrapper(nrows, ncols,
1542  const_cast<HYPRE_Int*>(I),
1543  const_cast<HYPRE_Int*>(J),
1544  const_cast<double*>(data),
1545  wrapper);
1546 }
1547 
1549 {
1550  MakeWrapper(A->diag, mem_diag, diag);
1551 }
1552 
1554 {
1555  MakeWrapper(A->offd, mem_offd, offd);
1556  cmap = A->col_map_offd;
1557 }
1558 
1560 {
1561  HostRead();
1562  hypre_CSRMatrix *hypre_merged = hypre_MergeDiagAndOffd(A);
1563  HypreRead();
1564  // Wrap 'hypre_merged' as a SparseMatrix 'merged_tmp'
1565  SparseMatrix merged_tmp;
1566 #if MFEM_HYPRE_VERSION >= 21600
1567  hypre_CSRMatrixBigJtoJ(hypre_merged);
1568 #endif
1569  MakeSparseMatrixWrapper(
1570  internal::to_int(hypre_merged->num_rows),
1571  internal::to_int(hypre_merged->num_cols),
1572  hypre_merged->i,
1573  hypre_merged->j,
1574  hypre_merged->data,
1575  merged_tmp);
1576  // Deep copy 'merged_tmp' to 'merged' so that 'merged' does not need
1577  // 'hypre_merged'
1578  merged = merged_tmp;
1579  merged_tmp.Clear();
1580  hypre_CSRMatrixDestroy(hypre_merged);
1581 }
1582 
1584  bool interleaved_rows,
1585  bool interleaved_cols) const
1586 {
1587  int nr = blocks.NumRows();
1588  int nc = blocks.NumCols();
1589 
1590  hypre_ParCSRMatrix **hypre_blocks = new hypre_ParCSRMatrix*[nr * nc];
1591  HostRead();
1592  internal::hypre_ParCSRMatrixSplit(A, nr, nc, hypre_blocks,
1593  interleaved_rows, interleaved_cols);
1594  HypreRead();
1595 
1596  for (int i = 0; i < nr; i++)
1597  {
1598  for (int j = 0; j < nc; j++)
1599  {
1600  blocks[i][j] = new HypreParMatrix(hypre_blocks[i*nc + j]);
1601  }
1602  }
1603 
1604  delete [] hypre_blocks;
1605 }
1606 
1608 {
1609  hypre_ParCSRMatrix * At;
1610  hypre_ParCSRMatrixTranspose(A, &At, 1);
1611  hypre_ParCSRMatrixSetNumNonzeros(At);
1612 
1613  if (!hypre_ParCSRMatrixCommPkg(At)) { hypre_MatvecCommPkgCreate(At); }
1614 
1615  if ( M() == N() )
1616  {
1617  /* If the matrix is square, make sure that the first entry in each
1618  row is the diagonal one. */
1619  hypre_CSRMatrixReorder(hypre_ParCSRMatrixDiag(At));
1620  }
1621 
1622  return new HypreParMatrix(At);
1623 }
1624 
1625 #if MFEM_HYPRE_VERSION >= 21800
1627  double threshold) const
1628 {
1629  // hypre_ParCSRMatrixExtractSubmatrixFC works on host only, so we move this
1630  // matrix to host, temporarily:
1631  HostRead();
1632 
1633  if (!(A->comm))
1634  {
1635  hypre_MatvecCommPkgCreate(A);
1636  }
1637 
1638  hypre_ParCSRMatrix *submat;
1639 
1640  // Get number of rows stored on this processor
1641  int local_num_vars = hypre_CSRMatrixNumRows(hypre_ParCSRMatrixDiag(A));
1642 
1643  // Form hypre CF-splitting array designating submatrix as F-points (-1)
1644 #ifdef hypre_IntArrayData
1645  // hypre_BoomerAMGCoarseParms needs CF_marker to be hypre_IntArray *
1646  hypre_IntArray *CF_marker;
1647 
1648  CF_marker = hypre_IntArrayCreate(local_num_vars);
1649  hypre_IntArrayInitialize_v2(CF_marker, HYPRE_MEMORY_HOST);
1650  hypre_IntArraySetConstantValues(CF_marker, 1);
1651 #else
1652  Array<HYPRE_Int> CF_marker(local_num_vars);
1653  CF_marker = 1;
1654 #endif
1655  for (int j=0; j<indices.Size(); j++)
1656  {
1657  if (indices[j] > local_num_vars)
1658  {
1659  MFEM_WARNING("WARNING : " << indices[j] << " > " << local_num_vars);
1660  }
1661 #ifdef hypre_IntArrayData
1662  hypre_IntArrayData(CF_marker)[indices[j]] = -1;
1663 #else
1664  CF_marker[indices[j]] = -1;
1665 #endif
1666  }
1667 
1668  // Construct cpts_global array on hypre matrix structure
1669 #if (MFEM_HYPRE_VERSION > 22300) || (MFEM_HYPRE_VERSION == 22300 && HYPRE_DEVELOP_NUMBER >=8)
1670  HYPRE_BigInt cpts_global[2];
1671 
1672  hypre_BoomerAMGCoarseParms(MPI_COMM_WORLD, local_num_vars, 1, NULL,
1673  CF_marker, NULL, cpts_global);
1674 #else
1675  HYPRE_BigInt *cpts_global;
1676  hypre_BoomerAMGCoarseParms(MPI_COMM_WORLD, local_num_vars, 1, NULL,
1677  CF_marker, NULL, &cpts_global);
1678 #endif
1679 
1680  // Extract submatrix into *submat
1681 #ifdef hypre_IntArrayData
1682  hypre_ParCSRMatrixExtractSubmatrixFC(A, hypre_IntArrayData(CF_marker),
1683  cpts_global, "FF", &submat,
1684  threshold);
1685 #else
1686  hypre_ParCSRMatrixExtractSubmatrixFC(A, CF_marker, cpts_global,
1687  "FF", &submat, threshold);
1688 #endif
1689 
1690 #if (MFEM_HYPRE_VERSION <= 22300) && !(MFEM_HYPRE_VERSION == 22300 && HYPRE_DEVELOP_NUMBER >=8)
1691  mfem_hypre_TFree(cpts_global);
1692 #endif
1693 #ifdef hypre_IntArrayData
1694  hypre_IntArrayDestroy(CF_marker);
1695 #endif
1696 
1697  HypreRead(); // restore the matrix location to the default hypre location
1698 
1699  return new HypreParMatrix(submat);
1700 }
1701 #endif
1702 
1704 {
1705 #if (MFEM_HYPRE_VERSION == 22500 && HYPRE_DEVELOP_NUMBER >= 1) || \
1706  (MFEM_HYPRE_VERSION > 22500)
1707 #ifdef HYPRE_USING_GPU
1708  hypre_ParCSRMatrixLocalTranspose(A);
1709 #endif
1710 #endif
1711 }
1712 
1714 {
1715 #if (MFEM_HYPRE_VERSION == 22500 && HYPRE_DEVELOP_NUMBER >= 1) || \
1716  (MFEM_HYPRE_VERSION > 22500)
1717 #ifdef HYPRE_USING_GPU
1718  if (A->diagT)
1719  {
1720  hypre_CSRMatrixDestroy(A->diagT);
1721  A->diagT = NULL;
1722  }
1723  if (A->offdT)
1724  {
1725  hypre_CSRMatrixDestroy(A->offdT);
1726  A->offdT = NULL;
1727  }
1728 #endif
1729 #endif
1730 }
1731 
1733  double a, double b) const
1734 {
1735  x.HypreRead();
1736  (b == 0.0) ? y.HypreWrite() : y.HypreReadWrite();
1737  return hypre_ParCSRMatrixMatvec(a, A, x, b, y);
1738 }
1739 
1740 void HypreParMatrix::Mult(double a, const Vector &x, double b, Vector &y) const
1741 {
1742  MFEM_ASSERT(x.Size() == Width(), "invalid x.Size() = " << x.Size()
1743  << ", expected size = " << Width());
1744  MFEM_ASSERT(y.Size() == Height(), "invalid y.Size() = " << y.Size()
1745  << ", expected size = " << Height());
1746 
1747  if (X == NULL)
1748  {
1749  X = new HypreParVector(A->comm,
1750  GetGlobalNumCols(),
1751  nullptr,
1752  GetColStarts());
1753  Y = new HypreParVector(A->comm,
1754  GetGlobalNumRows(),
1755  nullptr,
1756  GetRowStarts());
1757  }
1758 
1759  const bool xshallow = CanShallowCopy(x.GetMemory(), GetHypreMemoryClass());
1760  const bool yshallow = CanShallowCopy(y.GetMemory(), GetHypreMemoryClass());
1761 
1762  if (xshallow)
1763  {
1764  X->WrapMemoryRead(x.GetMemory());
1765  }
1766  else
1767  {
1768  if (auxX.Empty()) { auxX.New(NumCols(), GetHypreMemoryType()); }
1769  auxX.CopyFrom(x.GetMemory(), auxX.Capacity()); // Deep copy
1770  X->WrapMemoryRead(auxX);
1771  }
1772 
1773  if (yshallow)
1774  {
1775  if (b != 0.0) { Y->WrapMemoryReadWrite(y.GetMemory()); }
1776  else { Y->WrapMemoryWrite(y.GetMemory()); }
1777  }
1778  else
1779  {
1780  if (auxY.Empty()) { auxY.New(NumRows(), GetHypreMemoryType()); }
1781  if (b != 0.0)
1782  {
1783  auxY.CopyFrom(y.GetMemory(), auxY.Capacity()); // Deep copy
1784  Y->WrapMemoryReadWrite(auxY);
1785  }
1786  else
1787  {
1788  Y->WrapMemoryWrite(auxY);
1789  }
1790  }
1791 
1792  hypre_ParCSRMatrixMatvec(a, A, *X, b, *Y);
1793 
1794  if (!yshallow) { y = *Y; } // Deep copy
1795 }
1796 
1798  double b, Vector &y) const
1799 {
1800  MFEM_ASSERT(x.Size() == Height(), "invalid x.Size() = " << x.Size()
1801  << ", expected size = " << Height());
1802  MFEM_ASSERT(y.Size() == Width(), "invalid y.Size() = " << y.Size()
1803  << ", expected size = " << Width());
1804 
1805  // Note: x has the dimensions of Y (height), and
1806  // y has the dimensions of X (width)
1807  if (X == NULL)
1808  {
1809  X = new HypreParVector(A->comm,
1810  GetGlobalNumCols(),
1811  nullptr,
1812  GetColStarts());
1813  Y = new HypreParVector(A->comm,
1814  GetGlobalNumRows(),
1815  nullptr,
1816  GetRowStarts());
1817  }
1818 
1819  const bool xshallow = CanShallowCopy(x.GetMemory(), GetHypreMemoryClass());
1820  const bool yshallow = CanShallowCopy(y.GetMemory(), GetHypreMemoryClass());
1821 
1822  // x <--> Y
1823  if (xshallow)
1824  {
1825  Y->WrapMemoryRead(x.GetMemory());
1826  }
1827  else
1828  {
1829  if (auxY.Empty()) { auxY.New(NumRows(), GetHypreMemoryType()); }
1830  auxY.CopyFrom(x.GetMemory(), auxY.Capacity()); // Deep copy
1831  Y->WrapMemoryRead(auxY);
1832  }
1833 
1834  // y <--> X
1835  if (yshallow)
1836  {
1837  if (b != 0.0) { X->WrapMemoryReadWrite(y.GetMemory()); }
1838  else { X->WrapMemoryWrite(y.GetMemory()); }
1839  }
1840  else
1841  {
1842  if (auxX.Empty()) { auxX.New(NumCols(), GetHypreMemoryType()); }
1843  if (b != 0.0)
1844  {
1845  auxX.CopyFrom(y.GetMemory(), auxX.Capacity()); // Deep copy
1846  X->WrapMemoryReadWrite(auxX);
1847  }
1848  else
1849  {
1850  X->WrapMemoryWrite(auxX);
1851  }
1852  }
1853 
1855 
1856  hypre_ParCSRMatrixMatvecT(a, A, *Y, b, *X);
1857 
1858  if (!yshallow) { y = *X; } // Deep copy
1859 }
1860 
1861 HYPRE_Int HypreParMatrix::Mult(HYPRE_ParVector x, HYPRE_ParVector y,
1862  double a, double b) const
1863 {
1864  return hypre_ParCSRMatrixMatvec(a, A, (hypre_ParVector *) x, b,
1865  (hypre_ParVector *) y);
1866 }
1867 
1869  double a, double b) const
1870 {
1872  x.HypreRead();
1873  (b == 0.0) ? y.HypreWrite() : y.HypreReadWrite();
1874  return hypre_ParCSRMatrixMatvecT(a, A, x, b, y);
1875 }
1876 
1877 void HypreParMatrix::AbsMult(double a, const Vector &x,
1878  double b, Vector &y) const
1879 {
1880  MFEM_ASSERT(x.Size() == Width(), "invalid x.Size() = " << x.Size()
1881  << ", expected size = " << Width());
1882  MFEM_ASSERT(y.Size() == Height(), "invalid y.Size() = " << y.Size()
1883  << ", expected size = " << Height());
1884 
1885  auto x_data = x.HostRead();
1886  auto y_data = (b == 0.0) ? y.HostWrite() : y.HostReadWrite();
1887 
1888  HostRead();
1889  internal::hypre_ParCSRMatrixAbsMatvec(A, a, const_cast<double*>(x_data),
1890  b, y_data);
1891  HypreRead();
1892 }
1893 
1895  double b, Vector &y) const
1896 {
1897  MFEM_ASSERT(x.Size() == Height(), "invalid x.Size() = " << x.Size()
1898  << ", expected size = " << Height());
1899  MFEM_ASSERT(y.Size() == Width(), "invalid y.Size() = " << y.Size()
1900  << ", expected size = " << Width());
1901 
1902  auto x_data = x.HostRead();
1903  auto y_data = (b == 0.0) ? y.HostWrite() : y.HostReadWrite();
1904 
1905  HostRead();
1906  internal::hypre_ParCSRMatrixAbsMatvecT(A, a, const_cast<double*>(x_data),
1907  b, y_data);
1908  HypreRead();
1909 }
1910 
1912  HYPRE_BigInt* row_starts) const
1913 {
1914  const bool assumed_partition = HYPRE_AssumedPartitionCheck();
1915  const bool row_starts_given = (row_starts != NULL);
1916  if (!row_starts_given)
1917  {
1918  row_starts = hypre_ParCSRMatrixRowStarts(A);
1919  MFEM_VERIFY(D.Height() == hypre_CSRMatrixNumRows(A->diag),
1920  "the matrix D is NOT compatible with the row starts of"
1921  " this HypreParMatrix, row_starts must be given.");
1922  }
1923  else
1924  {
1925  int offset;
1926  if (assumed_partition)
1927  {
1928  offset = 0;
1929  }
1930  else
1931  {
1932  MPI_Comm_rank(GetComm(), &offset);
1933  }
1934  int local_num_rows = row_starts[offset+1]-row_starts[offset];
1935  MFEM_VERIFY(local_num_rows == D.Height(), "the number of rows in D is "
1936  " not compatible with the given row_starts");
1937  }
1938  // D.Width() will be checked for compatibility by the SparseMatrix
1939  // multiplication function, mfem::Mult(), called below.
1940 
1941  int part_size;
1942  HYPRE_BigInt global_num_rows;
1943  if (assumed_partition)
1944  {
1945  part_size = 2;
1946  if (row_starts_given)
1947  {
1948  global_num_rows = row_starts[2];
1949  // Here, we use row_starts[2], so row_starts must come from the
1950  // methods GetDofOffsets/GetTrueDofOffsets of ParFiniteElementSpace
1951  // (HYPRE's partitions have only 2 entries).
1952  }
1953  else
1954  {
1955  global_num_rows = hypre_ParCSRMatrixGlobalNumRows(A);
1956  }
1957  }
1958  else
1959  {
1960  MPI_Comm_size(GetComm(), &part_size);
1961  global_num_rows = row_starts[part_size];
1962  part_size++;
1963  }
1964 
1965  HYPRE_BigInt *col_starts = hypre_ParCSRMatrixColStarts(A);
1966  HYPRE_BigInt *col_map_offd;
1967 
1968  // get the diag and offd blocks as SparseMatrix wrappers
1969  SparseMatrix A_diag, A_offd;
1970  GetDiag(A_diag);
1971  GetOffd(A_offd, col_map_offd);
1972 
1973  // Multiply the diag and offd blocks with D -- these products will be the
1974  // diag and offd blocks of the output HypreParMatrix, DA.
1975  SparseMatrix* DA_diag = mfem::Mult(D, A_diag);
1976  SparseMatrix* DA_offd = mfem::Mult(D, A_offd);
1977 
1978  // Copy row_starts, col_starts, and col_map_offd; ownership of these arrays
1979  // will be given to the newly constructed output HypreParMatrix, DA.
1980  HYPRE_BigInt *new_row_starts =
1981  DuplicateAs<HYPRE_BigInt>(row_starts, part_size, false);
1982  HYPRE_BigInt *new_col_starts =
1983  (row_starts == col_starts ? new_row_starts :
1984  DuplicateAs<HYPRE_BigInt>(col_starts, part_size, false));
1985  HYPRE_BigInt *new_col_map_offd =
1986  DuplicateAs<HYPRE_BigInt>(col_map_offd, A_offd.Width());
1987 
1988  // Ownership of DA_diag and DA_offd is transferred to the HypreParMatrix
1989  // constructor.
1990  const bool own_diag_offd = true;
1991 
1992  // Create the output HypreParMatrix, DA, from DA_diag and DA_offd
1993  HypreParMatrix* DA =
1994  new HypreParMatrix(GetComm(),
1995  global_num_rows, hypre_ParCSRMatrixGlobalNumCols(A),
1996  new_row_starts, new_col_starts,
1997  DA_diag, DA_offd, new_col_map_offd,
1998  own_diag_offd);
1999 
2000 #if MFEM_HYPRE_VERSION <= 22200
2001  // Give ownership of row_starts, col_starts, and col_map_offd to DA
2002  hypre_ParCSRMatrixSetRowStartsOwner(DA->A, 1);
2003  hypre_ParCSRMatrixSetColStartsOwner(DA->A, 1);
2004 #else
2005  mfem_hypre_TFree_host(new_row_starts);
2006  mfem_hypre_TFree_host(new_col_starts);
2007 #endif
2008  DA->colMapOwner = 1;
2009 
2010  return DA;
2011 }
2012 
2014 {
2015  if (hypre_CSRMatrixNumRows(A->diag) != hypre_CSRMatrixNumRows(A->offd))
2016  {
2017  mfem_error("Row does not match");
2018  }
2019 
2020  if (hypre_CSRMatrixNumRows(A->diag) != diag.Size())
2021  {
2022  mfem_error("Note the Vector diag is not of compatible dimensions with A\n");
2023  }
2024 
2025  HostReadWrite();
2026  diag.HostRead();
2027 
2028  int size = Height();
2029  double *Adiag_data = hypre_CSRMatrixData(A->diag);
2030  HYPRE_Int *Adiag_i = hypre_CSRMatrixI(A->diag);
2031 
2032  double *Aoffd_data = hypre_CSRMatrixData(A->offd);
2033  HYPRE_Int *Aoffd_i = hypre_CSRMatrixI(A->offd);
2034  double val;
2035  HYPRE_Int jj;
2036  for (int i(0); i < size; ++i)
2037  {
2038  val = diag[i];
2039  for (jj = Adiag_i[i]; jj < Adiag_i[i+1]; ++jj)
2040  {
2041  Adiag_data[jj] *= val;
2042  }
2043  for (jj = Aoffd_i[i]; jj < Aoffd_i[i+1]; ++jj)
2044  {
2045  Aoffd_data[jj] *= val;
2046  }
2047  }
2048 
2049  HypreRead();
2050 }
2051 
2053 {
2054  if (hypre_CSRMatrixNumRows(A->diag) != hypre_CSRMatrixNumRows(A->offd))
2055  {
2056  mfem_error("Row does not match");
2057  }
2058 
2059  if (hypre_CSRMatrixNumRows(A->diag) != diag.Size())
2060  {
2061  mfem_error("Note the Vector diag is not of compatible dimensions with A\n");
2062  }
2063 
2064  HostReadWrite();
2065  diag.HostRead();
2066 
2067  int size = Height();
2068  double *Adiag_data = hypre_CSRMatrixData(A->diag);
2069  HYPRE_Int *Adiag_i = hypre_CSRMatrixI(A->diag);
2070 
2071 
2072  double *Aoffd_data = hypre_CSRMatrixData(A->offd);
2073  HYPRE_Int *Aoffd_i = hypre_CSRMatrixI(A->offd);
2074  double val;
2075  HYPRE_Int jj;
2076  for (int i(0); i < size; ++i)
2077  {
2078 #ifdef MFEM_DEBUG
2079  if (0.0 == diag(i))
2080  {
2081  mfem_error("HypreParMatrix::InvDiagScale : Division by 0");
2082  }
2083 #endif
2084  val = 1./diag(i);
2085  for (jj = Adiag_i[i]; jj < Adiag_i[i+1]; ++jj)
2086  {
2087  Adiag_data[jj] *= val;
2088  }
2089  for (jj = Aoffd_i[i]; jj < Aoffd_i[i+1]; ++jj)
2090  {
2091  Aoffd_data[jj] *= val;
2092  }
2093  }
2094 
2095  HypreRead();
2096 }
2097 
2099 {
2100  if (hypre_CSRMatrixNumRows(A->diag) != hypre_CSRMatrixNumRows(A->offd))
2101  {
2102  mfem_error("Row does not match");
2103  }
2104 
2105  HostReadWrite();
2106 
2107  HYPRE_Int size=hypre_CSRMatrixNumRows(A->diag);
2108  HYPRE_Int jj;
2109 
2110  double *Adiag_data = hypre_CSRMatrixData(A->diag);
2111  HYPRE_Int *Adiag_i = hypre_CSRMatrixI(A->diag);
2112  for (jj = 0; jj < Adiag_i[size]; ++jj)
2113  {
2114  Adiag_data[jj] *= s;
2115  }
2116 
2117  double *Aoffd_data = hypre_CSRMatrixData(A->offd);
2118  HYPRE_Int *Aoffd_i = hypre_CSRMatrixI(A->offd);
2119  for (jj = 0; jj < Aoffd_i[size]; ++jj)
2120  {
2121  Aoffd_data[jj] *= s;
2122  }
2123 
2124  HypreRead();
2125 }
2126 
2127 static void get_sorted_rows_cols(const Array<int> &rows_cols,
2128  Array<HYPRE_Int> &hypre_sorted)
2129 {
2130  rows_cols.HostRead();
2131  hypre_sorted.SetSize(rows_cols.Size());
2132  bool sorted = true;
2133  for (int i = 0; i < rows_cols.Size(); i++)
2134  {
2135  hypre_sorted[i] = rows_cols[i];
2136  if (i && rows_cols[i-1] > rows_cols[i]) { sorted = false; }
2137  }
2138  if (!sorted) { hypre_sorted.Sort(); }
2139 }
2140 
2141 void HypreParMatrix::Threshold(double threshold)
2142 {
2143  int ierr = 0;
2144 
2145  MPI_Comm comm;
2146  hypre_CSRMatrix * csr_A;
2147  hypre_CSRMatrix * csr_A_wo_z;
2148  hypre_ParCSRMatrix * parcsr_A_ptr;
2149  HYPRE_BigInt * row_starts = NULL; HYPRE_BigInt * col_starts = NULL;
2150  HYPRE_BigInt row_start = -1; HYPRE_BigInt row_end = -1;
2151  HYPRE_BigInt col_start = -1; HYPRE_BigInt col_end = -1;
2152 
2153  comm = hypre_ParCSRMatrixComm(A);
2154 
2155  ierr += hypre_ParCSRMatrixGetLocalRange(A,
2156  &row_start,&row_end,
2157  &col_start,&col_end );
2158 
2159  row_starts = hypre_ParCSRMatrixRowStarts(A);
2160  col_starts = hypre_ParCSRMatrixColStarts(A);
2161 
2162 #if MFEM_HYPRE_VERSION <= 22200
2163  bool old_owns_row = hypre_ParCSRMatrixOwnsRowStarts(A);
2164  bool old_owns_col = hypre_ParCSRMatrixOwnsColStarts(A);
2165 #endif
2166  HYPRE_BigInt global_num_rows = hypre_ParCSRMatrixGlobalNumRows(A);
2167  HYPRE_BigInt global_num_cols = hypre_ParCSRMatrixGlobalNumCols(A);
2168  parcsr_A_ptr = hypre_ParCSRMatrixCreate(comm, global_num_rows,
2169  global_num_cols,
2170  row_starts, col_starts,
2171  0, 0, 0);
2172 #if MFEM_HYPRE_VERSION <= 22200
2173  hypre_ParCSRMatrixOwnsRowStarts(parcsr_A_ptr) = old_owns_row;
2174  hypre_ParCSRMatrixOwnsColStarts(parcsr_A_ptr) = old_owns_col;
2175  hypre_ParCSRMatrixOwnsRowStarts(A) = 0;
2176  hypre_ParCSRMatrixOwnsColStarts(A) = 0;
2177 #endif
2178 
2179  csr_A = hypre_MergeDiagAndOffd(A);
2180 
2181  // Free A, if owned
2182  Destroy();
2183  Init();
2184 
2185  csr_A_wo_z = hypre_CSRMatrixDeleteZeros(csr_A,threshold);
2186 
2187  /* hypre_CSRMatrixDeleteZeros will return a NULL pointer rather than a usable
2188  CSR matrix if it finds no non-zeros */
2189  if (csr_A_wo_z == NULL)
2190  {
2191  csr_A_wo_z = csr_A;
2192  }
2193  else
2194  {
2195  ierr += hypre_CSRMatrixDestroy(csr_A);
2196  }
2197 
2198  /* TODO: GenerateDiagAndOffd() uses an int array of size equal to the number
2199  of columns in csr_A_wo_z which is the global number of columns in A. This
2200  does not scale well. */
2201  ierr += GenerateDiagAndOffd(csr_A_wo_z,parcsr_A_ptr,
2202  col_start,col_end);
2203 
2204  ierr += hypre_CSRMatrixDestroy(csr_A_wo_z);
2205 
2206  MFEM_VERIFY(ierr == 0, "");
2207 
2208  A = parcsr_A_ptr;
2209 
2210  hypre_ParCSRMatrixSetNumNonzeros(A);
2211  /* Make sure that the first entry in each row is the diagonal one. */
2212 #if MFEM_HYPRE_VERSION <= 22200
2213  if (row_starts == col_starts)
2214 #else
2215  if ((row_starts[0] == col_starts[0]) &&
2216  (row_starts[1] == col_starts[1]))
2217 #endif
2218  {
2219  hypre_CSRMatrixReorder(hypre_ParCSRMatrixDiag(A));
2220  }
2221  if (!hypre_ParCSRMatrixCommPkg(A)) { hypre_MatvecCommPkgCreate(A); }
2222  height = GetNumRows();
2223  width = GetNumCols();
2224 }
2225 
2227 {
2228  HYPRE_Int old_err = hypre_error_flag;
2229  hypre_error_flag = 0;
2230 
2231 #if MFEM_HYPRE_VERSION < 21400
2232 
2233  double threshold = 0.0;
2234  if (tol > 0.0)
2235  {
2236  HYPRE_Int *diag_I = A->diag->i, *offd_I = A->offd->i;
2237  double *diag_d = A->diag->data, *offd_d = A->offd->data;
2238  HYPRE_Int local_num_rows = A->diag->num_rows;
2239  double max_l2_row_norm = 0.0;
2240  Vector row;
2241  for (HYPRE_Int r = 0; r < local_num_rows; r++)
2242  {
2243  row.SetDataAndSize(diag_d + diag_I[r], diag_I[r+1]-diag_I[r]);
2244  double l2_row_norm = row.Norml2();
2245  row.SetDataAndSize(offd_d + offd_I[r], offd_I[r+1]-offd_I[r]);
2246  l2_row_norm = std::hypot(l2_row_norm, row.Norml2());
2247  max_l2_row_norm = std::max(max_l2_row_norm, l2_row_norm);
2248  }
2249  double loc_max_l2_row_norm = max_l2_row_norm;
2250  MPI_Allreduce(&loc_max_l2_row_norm, &max_l2_row_norm, 1, MPI_DOUBLE,
2251  MPI_MAX, A->comm);
2252  threshold = tol * max_l2_row_norm;
2253  }
2254 
2255  Threshold(threshold);
2256 
2257 #elif MFEM_HYPRE_VERSION < 21800
2258 
2259  HYPRE_Int err_flag = hypre_ParCSRMatrixDropSmallEntries(A, tol);
2260  MFEM_VERIFY(!err_flag, "error encountered: error code = " << err_flag);
2261 
2262 #else
2263 
2264  HYPRE_Int err_flag = hypre_ParCSRMatrixDropSmallEntries(A, tol, 2);
2265  MFEM_VERIFY(!err_flag, "error encountered: error code = " << err_flag);
2266 
2267 #endif
2268 
2269  hypre_error_flag = old_err;
2270 }
2271 
2273  const HypreParVector &x,
2274  HypreParVector &b)
2275 {
2276  Array<HYPRE_Int> rc_sorted;
2277  get_sorted_rows_cols(rows_cols, rc_sorted);
2278 
2279  internal::hypre_ParCSRMatrixEliminateAXB(
2280  A, rc_sorted.Size(), rc_sorted.GetData(), x, b);
2281 }
2282 
2284 {
2285  Array<HYPRE_Int> rc_sorted;
2286  get_sorted_rows_cols(rows_cols, rc_sorted);
2287 
2288  hypre_ParCSRMatrix* Ae;
2289  HostReadWrite();
2290  internal::hypre_ParCSRMatrixEliminateAAe(
2291  A, &Ae, rc_sorted.Size(), rc_sorted.GetData());
2292  HypreRead();
2293 
2294  return new HypreParMatrix(Ae, true);
2295 }
2296 
2298 {
2299  Array<HYPRE_Int> rc_sorted;
2300  get_sorted_rows_cols(cols, rc_sorted);
2301 
2302  hypre_ParCSRMatrix* Ae;
2303  HostReadWrite();
2304  internal::hypre_ParCSRMatrixEliminateAAe(
2305  A, &Ae, rc_sorted.Size(), rc_sorted.GetData(), 1);
2306  HypreRead();
2307 
2308  return new HypreParMatrix(Ae, true);
2309 }
2310 
2312 {
2313  if (rows.Size() > 0)
2314  {
2315  Array<HYPRE_Int> r_sorted;
2316  get_sorted_rows_cols(rows, r_sorted);
2317  HostReadWrite();
2318  internal::hypre_ParCSRMatrixEliminateRows(A, r_sorted.Size(),
2319  r_sorted.GetData());
2320  HypreRead();
2321  }
2322 }
2323 
2325  const Array<int> &ess_dof_list,
2326  const Vector &x, Vector &b) const
2327 {
2328  // b -= Ae*x
2329  Ae.Mult(-1.0, x, 1.0, b);
2330 
2331  // All operations below are local, so we can skip them if ess_dof_list is
2332  // empty on this processor to avoid potential host <--> device transfers.
2333  if (ess_dof_list.Size() == 0) { return; }
2334 
2335  HostRead();
2336  hypre_CSRMatrix *A_diag = hypre_ParCSRMatrixDiag(A);
2337  double *data = hypre_CSRMatrixData(A_diag);
2338  HYPRE_Int *I = hypre_CSRMatrixI(A_diag);
2339 #ifdef MFEM_DEBUG
2340  HYPRE_Int *J = hypre_CSRMatrixJ(A_diag);
2341  hypre_CSRMatrix *A_offd = hypre_ParCSRMatrixOffd(A);
2342  HYPRE_Int *I_offd = hypre_CSRMatrixI(A_offd);
2343  double *data_offd = hypre_CSRMatrixData(A_offd);
2344 #endif
2345 
2346  ess_dof_list.HostRead();
2347  x.HostRead();
2348  b.HostReadWrite();
2349 
2350  for (int i = 0; i < ess_dof_list.Size(); i++)
2351  {
2352  int r = ess_dof_list[i];
2353  b(r) = data[I[r]] * x(r);
2354 #ifdef MFEM_DEBUG
2355  MFEM_ASSERT(I[r] < I[r+1], "empty row found!");
2356  // Check that in the rows specified by the ess_dof_list, the matrix A has
2357  // only one entry -- the diagonal.
2358  // if (I[r+1] != I[r]+1 || J[I[r]] != r || I_offd[r] != I_offd[r+1])
2359  if (J[I[r]] != r)
2360  {
2361  MFEM_ABORT("the diagonal entry must be the first entry in the row!");
2362  }
2363  for (int j = I[r]+1; j < I[r+1]; j++)
2364  {
2365  if (data[j] != 0.0)
2366  {
2367  MFEM_ABORT("all off-diagonal entries must be zero!");
2368  }
2369  }
2370  for (int j = I_offd[r]; j < I_offd[r+1]; j++)
2371  {
2372  if (data_offd[j] != 0.0)
2373  {
2374  MFEM_ABORT("all off-diagonal entries must be zero!");
2375  }
2376  }
2377 #endif
2378  }
2379  HypreRead();
2380 }
2381 
2383  DiagonalPolicy diag_policy)
2384 {
2385  hypre_ParCSRMatrix *A_hypre = *this;
2386  HypreReadWrite();
2387 
2388  hypre_CSRMatrix *diag = hypre_ParCSRMatrixDiag(A_hypre);
2389  hypre_CSRMatrix *offd = hypre_ParCSRMatrixOffd(A_hypre);
2390 
2391  HYPRE_Int diag_nrows = hypre_CSRMatrixNumRows(diag);
2392  HYPRE_Int offd_ncols = hypre_CSRMatrixNumCols(offd);
2393 
2394  const int n_ess_dofs = ess_dofs.Size();
2395  const auto ess_dofs_d = ess_dofs.GetMemory().Read(
2396  GetHypreMemoryClass(), n_ess_dofs);
2397 
2398  // Start communication to figure out which columns need to be eliminated in
2399  // the off-diagonal block
2400  hypre_ParCSRCommHandle *comm_handle;
2401  HYPRE_Int *int_buf_data, *eliminate_row, *eliminate_col;
2402  {
2403  eliminate_row = mfem_hypre_CTAlloc(HYPRE_Int, diag_nrows);
2404  eliminate_col = mfem_hypre_CTAlloc(HYPRE_Int, offd_ncols);
2405 
2406  // Make sure A has a communication package
2407  hypre_ParCSRCommPkg *comm_pkg = hypre_ParCSRMatrixCommPkg(A_hypre);
2408  if (!comm_pkg)
2409  {
2410  hypre_MatvecCommPkgCreate(A_hypre);
2411  comm_pkg = hypre_ParCSRMatrixCommPkg(A_hypre);
2412  }
2413 
2414  // Which of the local rows are to be eliminated?
2415  MFEM_HYPRE_FORALL(i, diag_nrows, eliminate_row[i] = 0; );
2416  MFEM_HYPRE_FORALL(i, n_ess_dofs, eliminate_row[ess_dofs_d[i]] = 1; );
2417 
2418  // Use a matvec communication pattern to find (in eliminate_col) which of
2419  // the local offd columns are to be eliminated
2420 
2421  HYPRE_Int num_sends = hypre_ParCSRCommPkgNumSends(comm_pkg);
2422  HYPRE_Int int_buf_sz = hypre_ParCSRCommPkgSendMapStart(comm_pkg, num_sends);
2423  int_buf_data = mfem_hypre_CTAlloc(HYPRE_Int, int_buf_sz);
2424 
2425  HYPRE_Int *send_map_elmts;
2426 #if defined(HYPRE_USING_GPU)
2427  hypre_ParCSRCommPkgCopySendMapElmtsToDevice(comm_pkg);
2428  send_map_elmts = hypre_ParCSRCommPkgDeviceSendMapElmts(comm_pkg);
2429 #else
2430  send_map_elmts = hypre_ParCSRCommPkgSendMapElmts(comm_pkg);
2431 #endif
2432  MFEM_HYPRE_FORALL(i, int_buf_sz,
2433  {
2434  int k = send_map_elmts[i];
2435  int_buf_data[i] = eliminate_row[k];
2436  });
2437 
2438 #if defined(HYPRE_USING_GPU)
2439  // Try to use device-aware MPI for the communication if available
2440  comm_handle = hypre_ParCSRCommHandleCreate_v2(
2441  11, comm_pkg, HYPRE_MEMORY_DEVICE, int_buf_data,
2442  HYPRE_MEMORY_DEVICE, eliminate_col);
2443 #else
2444  comm_handle = hypre_ParCSRCommHandleCreate(
2445  11, comm_pkg, int_buf_data, eliminate_col );
2446 #endif
2447  }
2448 
2449  // Eliminate rows and columns in the diagonal block
2450  {
2451  const auto I = diag->i;
2452  const auto J = diag->j;
2453  auto data = diag->data;
2454 
2455  MFEM_HYPRE_FORALL(i, n_ess_dofs,
2456  {
2457  const int idof = ess_dofs_d[i];
2458  for (int j=I[idof]; j<I[idof+1]; ++j)
2459  {
2460  const int jdof = J[j];
2461  if (jdof == idof)
2462  {
2463  if (diag_policy == DiagonalPolicy::DIAG_ONE)
2464  {
2465  data[j] = 1.0;
2466  }
2467  else if (diag_policy == DiagonalPolicy::DIAG_ZERO)
2468  {
2469  data[j] = 0.0;
2470  }
2471  // else (diag_policy == DiagonalPolicy::DIAG_KEEP)
2472  }
2473  else
2474  {
2475  data[j] = 0.0;
2476  for (int k=I[jdof]; k<I[jdof+1]; ++k)
2477  {
2478  if (J[k] == idof)
2479  {
2480  data[k] = 0.0;
2481  break;
2482  }
2483  }
2484  }
2485  }
2486  });
2487  }
2488 
2489  // Eliminate rows in the off-diagonal block
2490  {
2491  const auto I = offd->i;
2492  auto data = offd->data;
2493  MFEM_HYPRE_FORALL(i, n_ess_dofs,
2494  {
2495  const int idof = ess_dofs_d[i];
2496  for (int j=I[idof]; j<I[idof+1]; ++j)
2497  {
2498  data[j] = 0.0;
2499  }
2500  });
2501  }
2502 
2503  // Wait for MPI communication to finish
2504  hypre_ParCSRCommHandleDestroy(comm_handle);
2505  mfem_hypre_TFree(int_buf_data);
2506  mfem_hypre_TFree(eliminate_row);
2507 
2508  // Eliminate columns in the off-diagonal block
2509  {
2510  const int nrows_offd = hypre_CSRMatrixNumRows(offd);
2511  const auto I = offd->i;
2512  const auto J = offd->j;
2513  auto data = offd->data;
2514  MFEM_HYPRE_FORALL(i, nrows_offd,
2515  {
2516  for (int j=I[i]; j<I[i+1]; ++j)
2517  {
2518  data[j] *= 1 - eliminate_col[J[j]];
2519  }
2520  });
2521  }
2522 
2523  mfem_hypre_TFree(eliminate_col);
2524 }
2525 
2526 void HypreParMatrix::Print(const char *fname, HYPRE_Int offi,
2527  HYPRE_Int offj) const
2528 {
2529  HostRead();
2530  hypre_ParCSRMatrixPrintIJ(A,offi,offj,fname);
2531  HypreRead();
2532 }
2533 
2534 void HypreParMatrix::Read(MPI_Comm comm, const char *fname)
2535 {
2536  Destroy();
2537  Init();
2538 
2539  HYPRE_Int base_i, base_j;
2540  hypre_ParCSRMatrixReadIJ(comm, fname, &base_i, &base_j, &A);
2541  hypre_ParCSRMatrixSetNumNonzeros(A);
2542 
2543  if (!hypre_ParCSRMatrixCommPkg(A)) { hypre_MatvecCommPkgCreate(A); }
2544 
2545  height = GetNumRows();
2546  width = GetNumCols();
2547 }
2548 
2549 void HypreParMatrix::Read_IJMatrix(MPI_Comm comm, const char *fname)
2550 {
2551  Destroy();
2552  Init();
2553 
2554  HYPRE_IJMatrix A_ij;
2555  HYPRE_IJMatrixRead(fname, comm, 5555, &A_ij); // HYPRE_PARCSR = 5555
2556 
2557  HYPRE_ParCSRMatrix A_parcsr;
2558  HYPRE_IJMatrixGetObject(A_ij, (void**) &A_parcsr);
2559 
2560  A = (hypre_ParCSRMatrix*)A_parcsr;
2561 
2562  hypre_ParCSRMatrixSetNumNonzeros(A);
2563 
2564  if (!hypre_ParCSRMatrixCommPkg(A)) { hypre_MatvecCommPkgCreate(A); }
2565 
2566  height = GetNumRows();
2567  width = GetNumCols();
2568 }
2569 
2570 void HypreParMatrix::PrintCommPkg(std::ostream &os) const
2571 {
2572  hypre_ParCSRCommPkg *comm_pkg = A->comm_pkg;
2573  MPI_Comm comm = A->comm;
2574  char c = '\0';
2575  const int tag = 46801;
2576  int myid, nproc;
2577  MPI_Comm_rank(comm, &myid);
2578  MPI_Comm_size(comm, &nproc);
2579 
2580  if (myid != 0)
2581  {
2582  MPI_Recv(&c, 1, MPI_CHAR, myid-1, tag, comm, MPI_STATUS_IGNORE);
2583  }
2584  else
2585  {
2586  os << "\nHypreParMatrix: hypre_ParCSRCommPkg:\n";
2587  }
2588  os << "Rank " << myid << ":\n"
2589  " number of sends = " << comm_pkg->num_sends <<
2590  " (" << sizeof(double)*comm_pkg->send_map_starts[comm_pkg->num_sends] <<
2591  " bytes)\n"
2592  " number of recvs = " << comm_pkg->num_recvs <<
2593  " (" << sizeof(double)*comm_pkg->recv_vec_starts[comm_pkg->num_recvs] <<
2594  " bytes)\n";
2595  if (myid != nproc-1)
2596  {
2597  os << std::flush;
2598  MPI_Send(&c, 1, MPI_CHAR, myid+1, tag, comm);
2599  }
2600  else
2601  {
2602  os << std::endl;
2603  }
2604  MPI_Barrier(comm);
2605 }
2606 
2607 void HypreParMatrix::PrintHash(std::ostream &os) const
2608 {
2609  HashFunction hf;
2610 
2611  os << "global number of rows : " << A->global_num_rows << '\n'
2612  << "global number of columns : " << A->global_num_cols << '\n'
2613  << "first row index : " << A->first_row_index << '\n'
2614  << " last row index : " << A->last_row_index << '\n'
2615  << "first col diag : " << A->first_col_diag << '\n'
2616  << " last col diag : " << A->last_col_diag << '\n'
2617  << "number of nonzeros : " << A->num_nonzeros << '\n';
2618  // diagonal, off-diagonal
2619  hypre_CSRMatrix *csr = A->diag;
2620  const char *csr_name = "diag";
2621  for (int m = 0; m < 2; m++)
2622  {
2623  auto csr_nnz = csr->i[csr->num_rows];
2624  os << csr_name << " num rows : " << csr->num_rows << '\n'
2625  << csr_name << " num cols : " << csr->num_cols << '\n'
2626  << csr_name << " num nnz : " << csr->num_nonzeros << '\n'
2627  << csr_name << " i last : " << csr_nnz
2628  << (csr_nnz == csr->num_nonzeros ?
2629  " [good]" : " [** BAD **]") << '\n';
2630  hf.AppendInts(csr->i, csr->num_rows + 1);
2631  os << csr_name << " i hash : " << hf.GetHash() << '\n';
2632  os << csr_name << " j hash : ";
2633  if (csr->j == nullptr)
2634  {
2635  os << "(null)\n";
2636  }
2637  else
2638  {
2639  hf.AppendInts(csr->j, csr_nnz);
2640  os << hf.GetHash() << '\n';
2641  }
2642 #if MFEM_HYPRE_VERSION >= 21600
2643  os << csr_name << " big j hash : ";
2644  if (csr->big_j == nullptr)
2645  {
2646  os << "(null)\n";
2647  }
2648  else
2649  {
2650  hf.AppendInts(csr->big_j, csr_nnz);
2651  os << hf.GetHash() << '\n';
2652  }
2653 #endif
2654  os << csr_name << " data hash : ";
2655  if (csr->data == nullptr)
2656  {
2657  os << "(null)\n";
2658  }
2659  else
2660  {
2661  hf.AppendDoubles(csr->data, csr_nnz);
2662  os << hf.GetHash() << '\n';
2663  }
2664 
2665  csr = A->offd;
2666  csr_name = "offd";
2667  }
2668 
2669  hf.AppendInts(A->col_map_offd, A->offd->num_cols);
2670  os << "col map offd hash : " << hf.GetHash() << '\n';
2671 }
2672 
2673 inline void delete_hypre_ParCSRMatrixColMapOffd(hypre_ParCSRMatrix *A)
2674 {
2675  HYPRE_BigInt *A_col_map_offd = hypre_ParCSRMatrixColMapOffd(A);
2676  int size = hypre_CSRMatrixNumCols(hypre_ParCSRMatrixOffd(A));
2677  Memory<HYPRE_BigInt>(A_col_map_offd, size, true).Delete();
2678 }
2679 
2680 void HypreParMatrix::Destroy()
2681 {
2682  if ( X != NULL ) { delete X; }
2683  if ( Y != NULL ) { delete Y; }
2684  auxX.Delete();
2685  auxY.Delete();
2686 
2687  if (A == NULL) { return; }
2688 
2689 #ifdef HYPRE_USING_GPU
2690  if (ParCSROwner && (diagOwner < 0 || offdOwner < 0))
2691  {
2692  // Put the "host" or "hypre" pointers in {i,j,data} of A->{diag,offd}, so
2693  // that they can be destroyed by hypre when hypre_ParCSRMatrixDestroy(A)
2694  // is called below.
2695 
2696  // Check that if both diagOwner and offdOwner are negative then they have
2697  // the same value.
2698  MFEM_VERIFY(!(diagOwner < 0 && offdOwner < 0) || diagOwner == offdOwner,
2699  "invalid state");
2700 
2701  MemoryClass mc = (diagOwner == -1 || offdOwner == -1) ?
2703  Write(mc, diagOwner < 0, offdOwner <0);
2704  }
2705 #endif
2706 
2707  mem_diag.I.Delete();
2708  mem_diag.J.Delete();
2709  mem_diag.data.Delete();
2710  if (diagOwner >= 0)
2711  {
2712  hypre_CSRMatrixI(A->diag) = NULL;
2713  hypre_CSRMatrixJ(A->diag) = NULL;
2714  hypre_CSRMatrixData(A->diag) = NULL;
2715  }
2716  mem_offd.I.Delete();
2717  mem_offd.J.Delete();
2718  mem_offd.data.Delete();
2719  if (offdOwner >= 0)
2720  {
2721  hypre_CSRMatrixI(A->offd) = NULL;
2722  hypre_CSRMatrixJ(A->offd) = NULL;
2723  hypre_CSRMatrixData(A->offd) = NULL;
2724  }
2725  if (colMapOwner >= 0)
2726  {
2727  if (colMapOwner & 1)
2728  {
2730  }
2731  hypre_ParCSRMatrixColMapOffd(A) = NULL;
2732  }
2733 
2734  if (ParCSROwner)
2735  {
2736  hypre_ParCSRMatrixDestroy(A);
2737  }
2738 }
2739 
2741 {
2742 #ifndef HYPRE_BIGINT
2743  bool own_i = A_hyp.GetDiagMemoryI().OwnsHostPtr();
2744  bool own_j = A_hyp.GetDiagMemoryJ().OwnsHostPtr();
2745  MFEM_CONTRACT_VAR(own_j);
2746  MFEM_ASSERT(own_i == own_j, "Inconsistent ownership");
2747  if (!own_i)
2748  {
2749  std::swap(A_diag.GetMemoryI(), A_hyp.GetDiagMemoryI());
2750  std::swap(A_diag.GetMemoryJ(), A_hyp.GetDiagMemoryJ());
2751  }
2752 #endif
2753  if (!A_hyp.GetDiagMemoryData().OwnsHostPtr())
2754  {
2755  std::swap(A_diag.GetMemoryData(), A_hyp.GetDiagMemoryData());
2756  }
2757  A_hyp.SetOwnerFlags(3, A_hyp.OwnsOffd(), A_hyp.OwnsColMap());
2758 }
2759 
2760 #if MFEM_HYPRE_VERSION >= 21800
2761 
2763  const Vector *b, HypreParVector *d,
2764  int blocksize, BlockInverseScaleJob job)
2765 {
2766  if (job == BlockInverseScaleJob::MATRIX_ONLY ||
2768  {
2769  hypre_ParCSRMatrix *C_hypre;
2770  hypre_ParcsrBdiagInvScal(*A, blocksize, &C_hypre);
2771  hypre_ParCSRMatrixDropSmallEntries(C_hypre, 1e-15, 1);
2772  C->WrapHypreParCSRMatrix(C_hypre);
2773  }
2774 
2775  if (job == BlockInverseScaleJob::RHS_ONLY ||
2777  {
2778  HypreParVector b_Hypre(A->GetComm(),
2779  A->GetGlobalNumRows(),
2780  b->GetData(), A->GetRowStarts());
2781  hypre_ParVector *d_hypre;
2782  hypre_ParvecBdiagInvScal(b_Hypre, blocksize, &d_hypre, *A);
2783 
2784  d->WrapHypreParVector(d_hypre, true);
2785  }
2786 }
2787 
2788 #endif
2789 
2790 #if MFEM_HYPRE_VERSION < 21400
2791 
2793  double beta, const HypreParMatrix &B)
2794 {
2795  hypre_ParCSRMatrix *C_hypre =
2796  internal::hypre_ParCSRMatrixAdd(const_cast<HypreParMatrix &>(A),
2797  const_cast<HypreParMatrix &>(B));
2798  MFEM_VERIFY(C_hypre, "error in hypre_ParCSRMatrixAdd");
2799 
2800  if (!hypre_ParCSRMatrixCommPkg(C_hypre)) { hypre_MatvecCommPkgCreate(C_hypre); }
2801  HypreParMatrix *C = new HypreParMatrix(C_hypre);
2802  *C = 0.0;
2803  C->Add(alpha, A);
2804  C->Add(beta, B);
2805 
2806  return C;
2807 }
2808 
2810 {
2811  hypre_ParCSRMatrix * C = internal::hypre_ParCSRMatrixAdd(*A,*B);
2812 
2813  if (!hypre_ParCSRMatrixCommPkg(C)) { hypre_MatvecCommPkgCreate(C); }
2814 
2815  return new HypreParMatrix(C);
2816 }
2817 
2818 #else
2819 
2820 HypreParMatrix *Add(double alpha, const HypreParMatrix &A,
2821  double beta, const HypreParMatrix &B)
2822 {
2823  hypre_ParCSRMatrix *C;
2824 #if MFEM_HYPRE_VERSION <= 22000
2825  hypre_ParcsrAdd(alpha, A, beta, B, &C);
2826 #else
2827  hypre_ParCSRMatrixAdd(alpha, A, beta, B, &C);
2828 #endif
2829  if (!hypre_ParCSRMatrixCommPkg(C)) { hypre_MatvecCommPkgCreate(C); }
2830 
2831  return new HypreParMatrix(C);
2832 }
2833 
2834 HypreParMatrix * ParAdd(const HypreParMatrix *A, const HypreParMatrix *B)
2835 {
2836  hypre_ParCSRMatrix *C;
2837 #if MFEM_HYPRE_VERSION <= 22000
2838  hypre_ParcsrAdd(1.0, *A, 1.0, *B, &C);
2839 #else
2840  hypre_ParCSRMatrixAdd(1.0, *A, 1.0, *B, &C);
2841 #endif
2842  if (!hypre_ParCSRMatrixCommPkg(C)) { hypre_MatvecCommPkgCreate(C); }
2843 
2844  return new HypreParMatrix(C);
2845 }
2846 
2847 #endif
2848 
2850  bool own_matrix)
2851 {
2852  hypre_ParCSRMatrix * ab;
2853 #ifdef HYPRE_USING_GPU
2854  ab = hypre_ParCSRMatMat(*A, *B);
2855 #else
2856  ab = hypre_ParMatmul(*A,*B);
2857 #endif
2858  hypre_ParCSRMatrixSetNumNonzeros(ab);
2859 
2860  if (!hypre_ParCSRMatrixCommPkg(ab)) { hypre_MatvecCommPkgCreate(ab); }
2861  HypreParMatrix *C = new HypreParMatrix(ab);
2862  if (own_matrix)
2863  {
2864  C->CopyRowStarts();
2865  C->CopyColStarts();
2866  }
2867  return C;
2868 }
2869 
2871 {
2872  hypre_ParCSRMatrix * rap;
2873 
2874 #ifdef HYPRE_USING_GPU
2875  // FIXME: this way of computing Pt A P can completely eliminate zero rows
2876  // from the sparsity pattern of the product which prevents
2877  // EliminateZeroRows() from working correctly. This issue is observed
2878  // in ex28p.
2879  // Quick fix: add a diagonal matrix with 0 diagonal.
2880  // Maybe use hypre_CSRMatrixCheckDiagFirst to see if we need the fix.
2881  {
2882  hypre_ParCSRMatrix *Q = hypre_ParCSRMatMat(*A,*P);
2883  const bool keepTranspose = false;
2884  rap = hypre_ParCSRTMatMatKT(*P,Q,keepTranspose);
2885  hypre_ParCSRMatrixDestroy(Q);
2886 
2887  // alternative:
2888  // hypre_ParCSRMatrixRAPKT
2889  }
2890 #else
2891 #if MFEM_HYPRE_VERSION <= 22200
2892  HYPRE_Int P_owns_its_col_starts =
2893  hypre_ParCSRMatrixOwnsColStarts((hypre_ParCSRMatrix*)(*P));
2894 #endif
2895 
2896  hypre_BoomerAMGBuildCoarseOperator(*P,*A,*P,&rap);
2897 
2898 #if MFEM_HYPRE_VERSION <= 22200
2899  /* Warning: hypre_BoomerAMGBuildCoarseOperator steals the col_starts
2900  from P (even if it does not own them)! */
2901  hypre_ParCSRMatrixSetRowStartsOwner(rap,0);
2902  hypre_ParCSRMatrixSetColStartsOwner(rap,0);
2903  if (P_owns_its_col_starts)
2904  {
2905  hypre_ParCSRMatrixSetColStartsOwner(*P, 1);
2906  }
2907 #endif
2908 #endif
2909 
2910  hypre_ParCSRMatrixSetNumNonzeros(rap);
2911  // hypre_MatvecCommPkgCreate(rap);
2912 
2913  return new HypreParMatrix(rap);
2914 }
2915 
2917  const HypreParMatrix *P)
2918 {
2919  hypre_ParCSRMatrix * rap;
2920 
2921 #ifdef HYPRE_USING_GPU
2922  {
2923  hypre_ParCSRMatrix *Q = hypre_ParCSRMatMat(*A,*P);
2924  rap = hypre_ParCSRTMatMat(*Rt,Q);
2925  hypre_ParCSRMatrixDestroy(Q);
2926  }
2927 #else
2928 #if MFEM_HYPRE_VERSION <= 22200
2929  HYPRE_Int P_owns_its_col_starts =
2930  hypre_ParCSRMatrixOwnsColStarts((hypre_ParCSRMatrix*)(*P));
2931  HYPRE_Int Rt_owns_its_col_starts =
2932  hypre_ParCSRMatrixOwnsColStarts((hypre_ParCSRMatrix*)(*Rt));
2933 #endif
2934 
2935  hypre_BoomerAMGBuildCoarseOperator(*Rt,*A,*P,&rap);
2936 
2937 #if MFEM_HYPRE_VERSION <= 22200
2938  /* Warning: hypre_BoomerAMGBuildCoarseOperator steals the col_starts
2939  from Rt and P (even if they do not own them)! */
2940  hypre_ParCSRMatrixSetRowStartsOwner(rap,0);
2941  hypre_ParCSRMatrixSetColStartsOwner(rap,0);
2942  if (P_owns_its_col_starts)
2943  {
2944  hypre_ParCSRMatrixSetColStartsOwner(*P, 1);
2945  }
2946  if (Rt_owns_its_col_starts)
2947  {
2948  hypre_ParCSRMatrixSetColStartsOwner(*Rt, 1);
2949  }
2950 #endif
2951 #endif
2952 
2953  hypre_ParCSRMatrixSetNumNonzeros(rap);
2954  // hypre_MatvecCommPkgCreate(rap);
2955 
2956  return new HypreParMatrix(rap);
2957 }
2958 
2959 // Helper function for HypreParMatrixFromBlocks. Note that scalability to
2960 // extremely large processor counts is limited by the use of MPI_Allgather.
2961 void GatherBlockOffsetData(MPI_Comm comm, const int rank, const int nprocs,
2962  const int num_loc, const Array<int> &offsets,
2963  std::vector<int> &all_num_loc, const int numBlocks,
2964  std::vector<std::vector<HYPRE_BigInt>> &blockProcOffsets,
2965  std::vector<HYPRE_BigInt> &procOffsets,
2966  std::vector<std::vector<int>> &procBlockOffsets,
2967  HYPRE_BigInt &firstLocal, HYPRE_BigInt &globalNum)
2968 {
2969  std::vector<std::vector<int>> all_block_num_loc(numBlocks);
2970 
2971  MPI_Allgather(&num_loc, 1, MPI_INT, all_num_loc.data(), 1, MPI_INT, comm);
2972 
2973  for (int j = 0; j < numBlocks; ++j)
2974  {
2975  all_block_num_loc[j].resize(nprocs);
2976  blockProcOffsets[j].resize(nprocs);
2977 
2978  const int blockNumRows = offsets[j + 1] - offsets[j];
2979  MPI_Allgather(&blockNumRows, 1, MPI_INT, all_block_num_loc[j].data(), 1,
2980  MPI_INT, comm);
2981  blockProcOffsets[j][0] = 0;
2982  for (int i = 0; i < nprocs - 1; ++i)
2983  {
2984  blockProcOffsets[j][i + 1] = blockProcOffsets[j][i]
2985  + all_block_num_loc[j][i];
2986  }
2987  }
2988 
2989  firstLocal = 0;
2990  globalNum = 0;
2991  procOffsets[0] = 0;
2992  for (int i = 0; i < nprocs; ++i)
2993  {
2994  globalNum += all_num_loc[i];
2995  if (rank == 0)
2996  {
2997  MFEM_VERIFY(globalNum >= 0, "overflow in global size");
2998  }
2999  if (i < rank)
3000  {
3001  firstLocal += all_num_loc[i];
3002  }
3003 
3004  if (i < nprocs - 1)
3005  {
3006  procOffsets[i + 1] = procOffsets[i] + all_num_loc[i];
3007  }
3008 
3009  procBlockOffsets[i].resize(numBlocks);
3010  procBlockOffsets[i][0] = 0;
3011  for (int j = 1; j < numBlocks; ++j)
3012  {
3013  procBlockOffsets[i][j] = procBlockOffsets[i][j - 1]
3014  + all_block_num_loc[j - 1][i];
3015  }
3016  }
3017 }
3018 
3020  Array2D<double> *blockCoeff)
3021 {
3022  const int numBlockRows = blocks.NumRows();
3023  const int numBlockCols = blocks.NumCols();
3024 
3025  MFEM_VERIFY(numBlockRows > 0 &&
3026  numBlockCols > 0, "Invalid input to HypreParMatrixFromBlocks");
3027 
3028  if (blockCoeff != NULL)
3029  {
3030  MFEM_VERIFY(numBlockRows == blockCoeff->NumRows() &&
3031  numBlockCols == blockCoeff->NumCols(),
3032  "Invalid input to HypreParMatrixFromBlocks");
3033  }
3034 
3035  Array<int> rowOffsets(numBlockRows+1);
3036  Array<int> colOffsets(numBlockCols+1);
3037 
3038  int nonNullBlockRow0 = -1;
3039  for (int j=0; j<numBlockCols; ++j)
3040  {
3041  if (blocks(0,j) != NULL)
3042  {
3043  nonNullBlockRow0 = j;
3044  break;
3045  }
3046  }
3047 
3048  MFEM_VERIFY(nonNullBlockRow0 >= 0, "Null row of blocks");
3049  MPI_Comm comm = blocks(0,nonNullBlockRow0)->GetComm();
3050 
3051  // Set offsets based on the number of rows or columns in each block.
3052  rowOffsets = 0;
3053  colOffsets = 0;
3054  for (int i=0; i<numBlockRows; ++i)
3055  {
3056  for (int j=0; j<numBlockCols; ++j)
3057  {
3058  if (blocks(i,j) != NULL)
3059  {
3060  const int nrows = blocks(i,j)->NumRows();
3061  const int ncols = blocks(i,j)->NumCols();
3062 
3063  MFEM_VERIFY(nrows > 0 &&
3064  ncols > 0, "Invalid block in HypreParMatrixFromBlocks");
3065 
3066  if (rowOffsets[i+1] == 0)
3067  {
3068  rowOffsets[i+1] = nrows;
3069  }
3070  else
3071  {
3072  MFEM_VERIFY(rowOffsets[i+1] == nrows,
3073  "Inconsistent blocks in HypreParMatrixFromBlocks");
3074  }
3075 
3076  if (colOffsets[j+1] == 0)
3077  {
3078  colOffsets[j+1] = ncols;
3079  }
3080  else
3081  {
3082  MFEM_VERIFY(colOffsets[j+1] == ncols,
3083  "Inconsistent blocks in HypreParMatrixFromBlocks");
3084  }
3085  }
3086  }
3087 
3088  MFEM_VERIFY(rowOffsets[i+1] > 0, "Invalid input blocks");
3089  rowOffsets[i+1] += rowOffsets[i];
3090  }
3091 
3092  for (int j=0; j<numBlockCols; ++j)
3093  {
3094  MFEM_VERIFY(colOffsets[j+1] > 0, "Invalid input blocks");
3095  colOffsets[j+1] += colOffsets[j];
3096  }
3097 
3098  const int num_loc_rows = rowOffsets[numBlockRows];
3099  const int num_loc_cols = colOffsets[numBlockCols];
3100 
3101  int nprocs, rank;
3102  MPI_Comm_rank(comm, &rank);
3103  MPI_Comm_size(comm, &nprocs);
3104 
3105  std::vector<int> all_num_loc_rows(nprocs);
3106  std::vector<int> all_num_loc_cols(nprocs);
3107  std::vector<HYPRE_BigInt> procRowOffsets(nprocs);
3108  std::vector<HYPRE_BigInt> procColOffsets(nprocs);
3109  std::vector<std::vector<HYPRE_BigInt>> blockRowProcOffsets(numBlockRows);
3110  std::vector<std::vector<HYPRE_BigInt>> blockColProcOffsets(numBlockCols);
3111  std::vector<std::vector<int>> procBlockRowOffsets(nprocs);
3112  std::vector<std::vector<int>> procBlockColOffsets(nprocs);
3113 
3114  HYPRE_BigInt first_loc_row, glob_nrows, first_loc_col, glob_ncols;
3115  GatherBlockOffsetData(comm, rank, nprocs, num_loc_rows, rowOffsets,
3116  all_num_loc_rows, numBlockRows, blockRowProcOffsets,
3117  procRowOffsets, procBlockRowOffsets, first_loc_row,
3118  glob_nrows);
3119 
3120  GatherBlockOffsetData(comm, rank, nprocs, num_loc_cols, colOffsets,
3121  all_num_loc_cols, numBlockCols, blockColProcOffsets,
3122  procColOffsets, procBlockColOffsets, first_loc_col,
3123  glob_ncols);
3124 
3125  std::vector<int> opI(num_loc_rows + 1);
3126  std::vector<int> cnt(num_loc_rows);
3127 
3128  for (int i = 0; i < num_loc_rows; ++i)
3129  {
3130  opI[i] = 0;
3131  cnt[i] = 0;
3132  }
3133 
3134  opI[num_loc_rows] = 0;
3135 
3136  Array2D<hypre_CSRMatrix *> csr_blocks(numBlockRows, numBlockCols);
3137 
3138  // Loop over all blocks, to determine nnz for each row.
3139  for (int i = 0; i < numBlockRows; ++i)
3140  {
3141  for (int j = 0; j < numBlockCols; ++j)
3142  {
3143  if (blocks(i, j) == NULL)
3144  {
3145  csr_blocks(i, j) = NULL;
3146  }
3147  else
3148  {
3149  blocks(i, j)->HostRead();
3150  csr_blocks(i, j) = hypre_MergeDiagAndOffd(*blocks(i, j));
3151  blocks(i, j)->HypreRead();
3152 
3153  for (int k = 0; k < csr_blocks(i, j)->num_rows; ++k)
3154  {
3155  opI[rowOffsets[i] + k + 1] +=
3156  csr_blocks(i, j)->i[k + 1] - csr_blocks(i, j)->i[k];
3157  }
3158  }
3159  }
3160  }
3161 
3162  // Now opI[i] is nnz for row i-1. Do a partial sum to get offsets.
3163  for (int i = 0; i < num_loc_rows; ++i)
3164  {
3165  opI[i + 1] += opI[i];
3166  }
3167 
3168  const int nnz = opI[num_loc_rows];
3169 
3170  std::vector<HYPRE_BigInt> opJ(nnz);
3171  std::vector<double> data(nnz);
3172 
3173  // Loop over all blocks, to set matrix data.
3174  for (int i = 0; i < numBlockRows; ++i)
3175  {
3176  for (int j = 0; j < numBlockCols; ++j)
3177  {
3178  if (csr_blocks(i, j) != NULL)
3179  {
3180  const int nrows = csr_blocks(i, j)->num_rows;
3181  const double cij = blockCoeff ? (*blockCoeff)(i, j) : 1.0;
3182 #if MFEM_HYPRE_VERSION >= 21600
3183  const bool usingBigJ = (csr_blocks(i, j)->big_j != NULL);
3184 #endif
3185 
3186  for (int k = 0; k < nrows; ++k)
3187  {
3188  const int rowg = rowOffsets[i] + k; // process-local row
3189  const int nnz_k = csr_blocks(i,j)->i[k+1]-csr_blocks(i,j)->i[k];
3190  const int osk = csr_blocks(i, j)->i[k];
3191 
3192  for (int l = 0; l < nnz_k; ++l)
3193  {
3194  // Find the column process offset for the block.
3195 #if MFEM_HYPRE_VERSION >= 21600
3196  const HYPRE_Int bcol = usingBigJ ?
3197  csr_blocks(i, j)->big_j[osk + l] :
3198  csr_blocks(i, j)->j[osk + l];
3199 #else
3200  const HYPRE_Int bcol = csr_blocks(i, j)->j[osk + l];
3201 #endif
3202 
3203  // find the processor 'bcolproc' that holds column 'bcol':
3204  const auto &offs = blockColProcOffsets[j];
3205  const int bcolproc =
3206  std::upper_bound(offs.begin() + 1, offs.end(), bcol)
3207  - offs.begin() - 1;
3208 
3209  opJ[opI[rowg] + cnt[rowg]] = procColOffsets[bcolproc] +
3210  procBlockColOffsets[bcolproc][j]
3211  + bcol
3212  - blockColProcOffsets[j][bcolproc];
3213  data[opI[rowg] + cnt[rowg]] = cij * csr_blocks(i, j)->data[osk + l];
3214  cnt[rowg]++;
3215  }
3216  }
3217  }
3218  }
3219  }
3220 
3221  for (int i = 0; i < numBlockRows; ++i)
3222  {
3223  for (int j = 0; j < numBlockCols; ++j)
3224  {
3225  if (csr_blocks(i, j) != NULL)
3226  {
3227  hypre_CSRMatrixDestroy(csr_blocks(i, j));
3228  }
3229  }
3230  }
3231 
3232  MFEM_VERIFY(HYPRE_AssumedPartitionCheck(),
3233  "only 'assumed partition' mode is supported");
3234 
3235  std::vector<HYPRE_BigInt> rowStarts2(2);
3236  rowStarts2[0] = first_loc_row;
3237  rowStarts2[1] = first_loc_row + all_num_loc_rows[rank];
3238 
3239  int square = std::equal(all_num_loc_rows.begin(), all_num_loc_rows.end(),
3240  all_num_loc_cols.begin());
3241  if (square)
3242  {
3243  return new HypreParMatrix(comm, num_loc_rows, glob_nrows, glob_ncols,
3244  opI.data(), opJ.data(),
3245  data.data(),
3246  rowStarts2.data(),
3247  rowStarts2.data());
3248  }
3249  else
3250  {
3251  std::vector<HYPRE_BigInt> colStarts2(2);
3252  colStarts2[0] = first_loc_col;
3253  colStarts2[1] = first_loc_col + all_num_loc_cols[rank];
3254 
3255  return new HypreParMatrix(comm, num_loc_rows, glob_nrows, glob_ncols,
3256  opI.data(), opJ.data(),
3257  data.data(),
3258  rowStarts2.data(),
3259  colStarts2.data());
3260  }
3261 }
3262 
3263 void EliminateBC(const HypreParMatrix &A, const HypreParMatrix &Ae,
3264  const Array<int> &ess_dof_list,
3265  const Vector &X, Vector &B)
3266 {
3267  A.EliminateBC(Ae, ess_dof_list, X, B);
3268 }
3269 
3270 // Taubin or "lambda-mu" scheme, which alternates between positive and
3271 // negative step sizes to approximate low-pass filter effect.
3272 
3273 int ParCSRRelax_Taubin(hypre_ParCSRMatrix *A, // matrix to relax with
3274  hypre_ParVector *f, // right-hand side
3275  double lambda,
3276  double mu,
3277  int N,
3278  double max_eig,
3279  hypre_ParVector *u, // initial/updated approximation
3280  hypre_ParVector *r // another temp vector
3281  )
3282 {
3283  hypre_CSRMatrix *A_diag = hypre_ParCSRMatrixDiag(A);
3284  HYPRE_Int num_rows = hypre_CSRMatrixNumRows(A_diag);
3285 
3286  double *u_data = hypre_VectorData(hypre_ParVectorLocalVector(u));
3287  double *r_data = hypre_VectorData(hypre_ParVectorLocalVector(r));
3288 
3289  for (int i = 0; i < N; i++)
3290  {
3291  // get residual: r = f - A*u
3292  hypre_ParVectorCopy(f, r);
3293  hypre_ParCSRMatrixMatvec(-1.0, A, u, 1.0, r);
3294 
3295  double coef;
3296  (0 == (i % 2)) ? coef = lambda : coef = mu;
3297 
3298  for (HYPRE_Int j = 0; j < num_rows; j++)
3299  {
3300  u_data[j] += coef*r_data[j] / max_eig;
3301  }
3302  }
3303 
3304  return 0;
3305 }
3306 
3307 // FIR scheme, which uses Chebyshev polynomials and a window function
3308 // to approximate a low-pass step filter.
3309 
3310 int ParCSRRelax_FIR(hypre_ParCSRMatrix *A, // matrix to relax with
3311  hypre_ParVector *f, // right-hand side
3312  double max_eig,
3313  int poly_order,
3314  double* fir_coeffs,
3315  hypre_ParVector *u, // initial/updated approximation
3316  hypre_ParVector *x0, // temporaries
3317  hypre_ParVector *x1,
3318  hypre_ParVector *x2,
3319  hypre_ParVector *x3)
3320 
3321 {
3322  hypre_CSRMatrix *A_diag = hypre_ParCSRMatrixDiag(A);
3323  HYPRE_Int num_rows = hypre_CSRMatrixNumRows(A_diag);
3324 
3325  double *u_data = hypre_VectorData(hypre_ParVectorLocalVector(u));
3326 
3327  double *x0_data = hypre_VectorData(hypre_ParVectorLocalVector(x0));
3328  double *x1_data = hypre_VectorData(hypre_ParVectorLocalVector(x1));
3329  double *x2_data = hypre_VectorData(hypre_ParVectorLocalVector(x2));
3330  double *x3_data = hypre_VectorData(hypre_ParVectorLocalVector(x3));
3331 
3332  hypre_ParVectorCopy(u, x0);
3333 
3334  // x1 = f -A*x0/max_eig
3335  hypre_ParVectorCopy(f, x1);
3336  hypre_ParCSRMatrixMatvec(-1.0, A, x0, 1.0, x1);
3337 
3338  for (HYPRE_Int i = 0; i < num_rows; i++)
3339  {
3340  x1_data[i] /= -max_eig;
3341  }
3342 
3343  // x1 = x0 -x1
3344  for (HYPRE_Int i = 0; i < num_rows; i++)
3345  {
3346  x1_data[i] = x0_data[i] -x1_data[i];
3347  }
3348 
3349  // x3 = f0*x0 +f1*x1
3350  for (HYPRE_Int i = 0; i < num_rows; i++)
3351  {
3352  x3_data[i] = fir_coeffs[0]*x0_data[i] +fir_coeffs[1]*x1_data[i];
3353  }
3354 
3355  for (int n = 2; n <= poly_order; n++)
3356  {
3357  // x2 = f - A*x1/max_eig
3358  hypre_ParVectorCopy(f, x2);
3359  hypre_ParCSRMatrixMatvec(-1.0, A, x1, 1.0, x2);
3360 
3361  for (HYPRE_Int i = 0; i < num_rows; i++)
3362  {
3363  x2_data[i] /= -max_eig;
3364  }
3365 
3366  // x2 = (x1-x0) +(x1-2*x2)
3367  // x3 = x3 +f[n]*x2
3368  // x0 = x1
3369  // x1 = x2
3370 
3371  for (HYPRE_Int i = 0; i < num_rows; i++)
3372  {
3373  x2_data[i] = (x1_data[i]-x0_data[i]) +(x1_data[i]-2*x2_data[i]);
3374  x3_data[i] += fir_coeffs[n]*x2_data[i];
3375  x0_data[i] = x1_data[i];
3376  x1_data[i] = x2_data[i];
3377  }
3378  }
3379 
3380  for (HYPRE_Int i = 0; i < num_rows; i++)
3381  {
3382  u_data[i] = x3_data[i];
3383  }
3384 
3385  return 0;
3386 }
3387 
3389 {
3390  type = default_type;
3391  relax_times = 1;
3392  relax_weight = 1.0;
3393  omega = 1.0;
3394  poly_order = 2;
3395  poly_fraction = .3;
3396  lambda = 0.5;
3397  mu = -0.5;
3398  taubin_iter = 40;
3399 
3400  l1_norms = NULL;
3401  pos_l1_norms = false;
3402  eig_est_cg_iter = 10;
3403  B = X = V = Z = NULL;
3404  auxB.Reset(); auxX.Reset();
3405  X0 = X1 = NULL;
3406  fir_coeffs = NULL;
3407  A_is_symmetric = false;
3408 }
3409 
3411  int relax_times_, double relax_weight_,
3412  double omega_, int poly_order_,
3413  double poly_fraction_, int eig_est_cg_iter_)
3414 {
3415  type = type_;
3416  relax_times = relax_times_;
3417  relax_weight = relax_weight_;
3418  omega = omega_;
3419  poly_order = poly_order_;
3420  poly_fraction = poly_fraction_;
3421  eig_est_cg_iter = eig_est_cg_iter_;
3422 
3423  l1_norms = NULL;
3424  pos_l1_norms = false;
3425  B = X = V = Z = NULL;
3426  auxB.Reset(); auxX.Reset();
3427  X0 = X1 = NULL;
3428  fir_coeffs = NULL;
3429  A_is_symmetric = false;
3430 
3431  SetOperator(A_);
3432 }
3433 
3434 void HypreSmoother::SetType(HypreSmoother::Type type_, int relax_times_)
3435 {
3436  type = static_cast<int>(type_);
3437  relax_times = relax_times_;
3438 }
3439 
3440 void HypreSmoother::SetSOROptions(double relax_weight_, double omega_)
3441 {
3442  relax_weight = relax_weight_;
3443  omega = omega_;
3444 }
3445 
3446 void HypreSmoother::SetPolyOptions(int poly_order_, double poly_fraction_,
3447  int eig_est_cg_iter_)
3448 {
3449  poly_order = poly_order_;
3450  poly_fraction = poly_fraction_;
3451  eig_est_cg_iter = eig_est_cg_iter_;
3452 }
3453 
3454 void HypreSmoother::SetTaubinOptions(double lambda_, double mu_,
3455  int taubin_iter_)
3456 {
3457  lambda = lambda_;
3458  mu = mu_;
3459  taubin_iter = taubin_iter_;
3460 }
3461 
3462 void HypreSmoother::SetWindowByName(const char* name)
3463 {
3464  double a = -1, b, c;
3465  if (!strcmp(name,"Rectangular")) { a = 1.0, b = 0.0, c = 0.0; }
3466  if (!strcmp(name,"Hanning")) { a = 0.5, b = 0.5, c = 0.0; }
3467  if (!strcmp(name,"Hamming")) { a = 0.54, b = 0.46, c = 0.0; }
3468  if (!strcmp(name,"Blackman")) { a = 0.42, b = 0.50, c = 0.08; }
3469  if (a < 0)
3470  {
3471  mfem_error("HypreSmoother::SetWindowByName : name not recognized!");
3472  }
3473 
3474  SetWindowParameters(a, b, c);
3475 }
3476 
3477 void HypreSmoother::SetWindowParameters(double a, double b, double c)
3478 {
3479  window_params[0] = a;
3480  window_params[1] = b;
3481  window_params[2] = c;
3482 }
3483 
3485 {
3486  A = const_cast<HypreParMatrix *>(dynamic_cast<const HypreParMatrix *>(&op));
3487  if (A == NULL)
3488  {
3489  mfem_error("HypreSmoother::SetOperator : not HypreParMatrix!");
3490  }
3491 
3492  height = A->Height();
3493  width = A->Width();
3494 
3495  auxX.Delete(); auxB.Delete();
3496  if (B) { delete B; }
3497  if (X) { delete X; }
3498  if (V) { delete V; }
3499  if (Z) { delete Z; }
3500  if (l1_norms)
3501  {
3502  mfem_hypre_TFree(l1_norms);
3503  }
3504  delete X0;
3505  delete X1;
3506 
3507  X1 = X0 = Z = V = B = X = NULL;
3508  auxB.Reset(); auxX.Reset();
3509 
3510  if (type >= 1 && type <= 4)
3511  {
3512  hypre_ParCSRComputeL1Norms(*A, type, NULL, &l1_norms);
3513  // The above call will set the hypre_error_flag when it encounters zero
3514  // rows in A.
3515  }
3516  else if (type == 5)
3517  {
3518  l1_norms = mfem_hypre_CTAlloc(double, height);
3519  Vector ones(height), diag(l1_norms, height);
3520  ones = 1.0;
3521  A->Mult(ones, diag);
3522  }
3523  else
3524  {
3525  l1_norms = NULL;
3526  }
3527  if (l1_norms && pos_l1_norms)
3528  {
3529 #if defined(HYPRE_USING_GPU)
3530  double *d_l1_norms = l1_norms; // avoid *this capture
3531  MFEM_GPU_FORALL(i, height,
3532  {
3533  d_l1_norms[i] = std::abs(d_l1_norms[i]);
3534  });
3535 #else
3536  for (int i = 0; i < height; i++)
3537  {
3538  l1_norms[i] = std::abs(l1_norms[i]);
3539  }
3540 #endif
3541  }
3542 
3543  if (type == 16)
3544  {
3545  poly_scale = 1;
3546  if (eig_est_cg_iter > 0)
3547  {
3548  hypre_ParCSRMaxEigEstimateCG(*A, poly_scale, eig_est_cg_iter,
3550  }
3551  else
3552  {
3553 #if MFEM_HYPRE_VERSION <= 22200
3554  min_eig_est = 0;
3555  hypre_ParCSRMaxEigEstimate(*A, poly_scale, &max_eig_est);
3556 #else
3557  hypre_ParCSRMaxEigEstimate(*A, poly_scale, &max_eig_est, &min_eig_est);
3558 #endif
3559  }
3560  Z = new HypreParVector(*A);
3561  }
3562  else if (type == 1001 || type == 1002)
3563  {
3564  poly_scale = 0;
3565  if (eig_est_cg_iter > 0)
3566  {
3567  hypre_ParCSRMaxEigEstimateCG(*A, poly_scale, eig_est_cg_iter,
3569  }
3570  else
3571  {
3572 #if MFEM_HYPRE_VERSION <= 22200
3573  min_eig_est = 0;
3574  hypre_ParCSRMaxEigEstimate(*A, poly_scale, &max_eig_est);
3575 #else
3576  hypre_ParCSRMaxEigEstimate(*A, poly_scale, &max_eig_est, &min_eig_est);
3577 #endif
3578  }
3579 
3580  // The Taubin and FIR polynomials are defined on [0, 2]
3581  max_eig_est /= 2;
3582 
3583  // Compute window function, Chebyshev coefficients, and allocate temps.
3584  if (type == 1002)
3585  {
3586  // Temporaries for Chebyshev recursive evaluation
3587  Z = new HypreParVector(*A);
3588  X0 = new HypreParVector(*A);
3589  X1 = new HypreParVector(*A);
3590 
3592  }
3593  }
3594 }
3595 
3597 {
3598  if (fir_coeffs)
3599  {
3600  delete [] fir_coeffs;
3601  }
3602 
3603  fir_coeffs = new double[poly_order+1];
3604 
3605  double* window_coeffs = new double[poly_order+1];
3606  double* cheby_coeffs = new double[poly_order+1];
3607 
3608  double a = window_params[0];
3609  double b = window_params[1];
3610  double c = window_params[2];
3611  for (int i = 0; i <= poly_order; i++)
3612  {
3613  double t = (i*M_PI)/(poly_order+1);
3614  window_coeffs[i] = a + b*cos(t) +c*cos(2*t);
3615  }
3616 
3617  double k_pb = poly_fraction*max_eig;
3618  double theta_pb = acos(1.0 -0.5*k_pb);
3619  double sigma = 0.0;
3620  cheby_coeffs[0] = (theta_pb +sigma)/M_PI;
3621  for (int i = 1; i <= poly_order; i++)
3622  {
3623  double t = i*(theta_pb+sigma);
3624  cheby_coeffs[i] = 2.0*sin(t)/(i*M_PI);
3625  }
3626 
3627  for (int i = 0; i <= poly_order; i++)
3628  {
3629  fir_coeffs[i] = window_coeffs[i]*cheby_coeffs[i];
3630  }
3631 
3632  delete[] window_coeffs;
3633  delete[] cheby_coeffs;
3634 }
3635 
3637 {
3638  if (A == NULL)
3639  {
3640  mfem_error("HypreSmoother::Mult (...) : HypreParMatrix A is missing");
3641  return;
3642  }
3643 
3644  // TODO: figure out where each function needs A, b, and x ...
3645 
3646  b.HypreRead();
3647  if (!iterative_mode)
3648  {
3649  x.HypreWrite();
3650  if (type == 0 && relax_times == 1)
3651  {
3652  // Note: hypre_ParCSRDiagScale() is not exposed in older versions
3653  HYPRE_ParCSRDiagScale(NULL, *A, b, x);
3654  if (relax_weight != 1.0)
3655  {
3656  hypre_ParVectorScale(relax_weight, x);
3657  }
3658  return;
3659  }
3660  hypre_ParVectorSetConstantValues(x, 0.0);
3661  }
3662  else
3663  {
3664  x.HypreReadWrite();
3665  }
3666 
3667  if (V == NULL)
3668  {
3669  V = new HypreParVector(*A);
3670  }
3671 
3672  if (type == 1001)
3673  {
3674  for (int sweep = 0; sweep < relax_times; sweep++)
3675  {
3677  max_eig_est,
3678  x, *V);
3679  }
3680  }
3681  else if (type == 1002)
3682  {
3683  for (int sweep = 0; sweep < relax_times; sweep++)
3684  {
3685  ParCSRRelax_FIR(*A, b,
3686  max_eig_est,
3687  poly_order,
3688  fir_coeffs,
3689  x,
3690  *X0, *X1, *V, *Z);
3691  }
3692  }
3693  else
3694  {
3695  int hypre_type = type;
3696  // hypre doesn't have lumped Jacobi, so treat the action as l1-Jacobi
3697  if (type == 5) { hypre_type = 1; }
3698 
3699  if (Z == NULL)
3700  {
3701  hypre_ParCSRRelax(*A, b, hypre_type,
3704  x, *V, NULL);
3705  }
3706  else
3707  {
3708  hypre_ParCSRRelax(*A, b, hypre_type,
3711  x, *V, *Z);
3712  }
3713  }
3714 }
3715 
3716 void HypreSmoother::Mult(const Vector &b, Vector &x) const
3717 {
3718  MFEM_ASSERT(b.Size() == NumCols(), "");
3719  MFEM_ASSERT(x.Size() == NumRows(), "");
3720 
3721  if (A == NULL)
3722  {
3723  mfem_error("HypreSmoother::Mult (...) : HypreParMatrix A is missing");
3724  return;
3725  }
3726 
3727  if (B == NULL)
3728  {
3729  B = new HypreParVector(A->GetComm(),
3730  A -> GetGlobalNumRows(),
3731  nullptr,
3732  A -> GetRowStarts());
3733  X = new HypreParVector(A->GetComm(),
3734  A -> GetGlobalNumCols(),
3735  nullptr,
3736  A -> GetColStarts());
3737  }
3738 
3739  const bool bshallow = CanShallowCopy(b.GetMemory(), GetHypreMemoryClass());
3740  const bool xshallow = CanShallowCopy(x.GetMemory(), GetHypreMemoryClass());
3741 
3742  if (bshallow)
3743  {
3744  B->WrapMemoryRead(b.GetMemory());
3745  }
3746  else
3747  {
3748  if (auxB.Empty()) { auxB.New(NumCols(), GetHypreMemoryType()); }
3749  auxB.CopyFrom(b.GetMemory(), auxB.Capacity()); // Deep copy
3750  B->WrapMemoryRead(auxB);
3751  }
3752 
3753  if (xshallow)
3754  {
3756  else { X->WrapMemoryWrite(x.GetMemory()); }
3757  }
3758  else
3759  {
3760  if (auxX.Empty()) { auxX.New(NumRows(), GetHypreMemoryType()); }
3761  if (iterative_mode)
3762  {
3763  auxX.CopyFrom(x.GetMemory(), x.Size()); // Deep copy
3765  }
3766  else
3767  {
3769  }
3770  }
3771 
3772  Mult(*B, *X);
3773 
3774  if (!xshallow) { x = *X; } // Deep copy
3775 }
3776 
3778 {
3779  if (A_is_symmetric || type == 0 || type == 1 || type == 5)
3780  {
3781  Mult(b, x);
3782  return;
3783  }
3784  mfem_error("HypreSmoother::MultTranspose (...) : undefined!\n");
3785 }
3786 
3788 {
3789  auxX.Delete(); auxB.Delete();
3790  if (B) { delete B; }
3791  if (X) { delete X; }
3792  if (V) { delete V; }
3793  if (Z) { delete Z; }
3794  if (l1_norms)
3795  {
3796  mfem_hypre_TFree(l1_norms);
3797  }
3798  if (fir_coeffs)
3799  {
3800  delete [] fir_coeffs;
3801  }
3802  if (X0) { delete X0; }
3803  if (X1) { delete X1; }
3804 }
3805 
3806 
3808 {
3809  A = NULL;
3810  setup_called = 0;
3811  B = X = NULL;
3812  auxB.Reset();
3813  auxX.Reset();
3815 }
3816 
3818  : Solver(A_->Height(), A_->Width())
3819 {
3820  A = A_;
3821  setup_called = 0;
3822  B = X = NULL;
3823  auxB.Reset();
3824  auxX.Reset();
3826 }
3827 
3828 bool HypreSolver::WrapVectors(const Vector &b, Vector &x) const
3829 {
3830  MFEM_ASSERT(b.Size() == NumCols(), "");
3831  MFEM_ASSERT(x.Size() == NumRows(), "");
3832 
3833  MFEM_VERIFY(A != NULL, "HypreParMatrix A is missing");
3834 
3835  if (B == NULL)
3836  {
3837  B = new HypreParVector(A->GetComm(), A->GetGlobalNumRows(),
3838  nullptr, A->GetRowStarts());
3839  X = new HypreParVector(A->GetComm(), A->GetGlobalNumCols(),
3840  nullptr, A->GetColStarts());
3841  }
3842 
3843  const bool bshallow = CanShallowCopy(b.GetMemory(), GetHypreMemoryClass());
3844  const bool xshallow = CanShallowCopy(x.GetMemory(), GetHypreMemoryClass());
3845 
3846  if (bshallow)
3847  {
3848  B->WrapMemoryRead(b.GetMemory());
3849  }
3850  else
3851  {
3852  if (auxB.Empty()) { auxB.New(NumCols(), GetHypreMemoryType()); }
3853  auxB.CopyFrom(b.GetMemory(), auxB.Capacity()); // Deep copy
3854  B->WrapMemoryRead(auxB);
3855  }
3856 
3857  if (xshallow)
3858  {
3860  else { X->WrapMemoryWrite(x.GetMemory()); }
3861  }
3862  else
3863  {
3864  if (auxX.Empty()) { auxX.New(NumRows(), GetHypreMemoryType()); }
3865  if (iterative_mode)
3866  {
3867  auxX.CopyFrom(x.GetMemory(), x.Size()); // Deep copy
3869  }
3870  else
3871  {
3873  }
3874  }
3875 
3876  return xshallow;
3877 }
3878 
3880 {
3881  if (setup_called) { return; }
3882 
3883  MFEM_VERIFY(A != NULL, "HypreParMatrix A is missing");
3884 
3885  HYPRE_Int err_flag = SetupFcn()(*this, *A, b, x);
3887  {
3888  if (err_flag)
3889  { MFEM_WARNING("Error during setup! Error code: " << err_flag); }
3890  }
3891  else if (error_mode == ABORT_HYPRE_ERRORS)
3892  {
3893  MFEM_VERIFY(!err_flag, "Error during setup! Error code: " << err_flag);
3894  }
3895  hypre_error_flag = 0;
3896  setup_called = 1;
3897 }
3898 
3899 void HypreSolver::Setup(const Vector &b, Vector &x) const
3900 {
3901  const bool x_shallow = WrapVectors(b, x);
3902  Setup(*B, *X);
3903  if (!x_shallow) { x = *X; } // Deep copy if shallow copy is impossible
3904 }
3905 
3907 {
3908  HYPRE_Int err_flag;
3909  if (A == NULL)
3910  {
3911  mfem_error("HypreSolver::Mult (...) : HypreParMatrix A is missing");
3912  return;
3913  }
3914 
3915  if (!iterative_mode)
3916  {
3917  x.HypreWrite();
3918  hypre_ParVectorSetConstantValues(x, 0.0);
3919  }
3920 
3921  b.HypreRead();
3922  x.HypreReadWrite();
3923 
3924  Setup(b, x);
3925 
3926  err_flag = SolveFcn()(*this, *A, b, x);
3928  {
3929  if (err_flag)
3930  { MFEM_WARNING("Error during solve! Error code: " << err_flag); }
3931  }
3932  else if (error_mode == ABORT_HYPRE_ERRORS)
3933  {
3934  MFEM_VERIFY(!err_flag, "Error during solve! Error code: " << err_flag);
3935  }
3936  hypre_error_flag = 0;
3937 }
3938 
3939 void HypreSolver::Mult(const Vector &b, Vector &x) const
3940 {
3941  const bool x_shallow = WrapVectors(b, x);
3942  Mult(*B, *X);
3943  if (!x_shallow) { x = *X; } // Deep copy if shallow copy is impossible
3944 }
3945 
3947 {
3948  if (B) { delete B; }
3949  if (X) { delete X; }
3950  auxB.Delete();
3951  auxX.Delete();
3952 }
3953 
3954 
3955 HyprePCG::HyprePCG(MPI_Comm comm) : precond(NULL)
3956 {
3957  iterative_mode = true;
3958 
3959  HYPRE_ParCSRPCGCreate(comm, &pcg_solver);
3960 }
3961 
3962 HyprePCG::HyprePCG(const HypreParMatrix &A_) : HypreSolver(&A_), precond(NULL)
3963 {
3964  MPI_Comm comm;
3965 
3966  iterative_mode = true;
3967 
3968  HYPRE_ParCSRMatrixGetComm(*A, &comm);
3969 
3970  HYPRE_ParCSRPCGCreate(comm, &pcg_solver);
3971 }
3972 
3974 {
3975  const HypreParMatrix *new_A = dynamic_cast<const HypreParMatrix *>(&op);
3976  MFEM_VERIFY(new_A, "new Operator must be a HypreParMatrix!");
3977 
3978  // update base classes: Operator, Solver, HypreSolver
3979  height = new_A->Height();
3980  width = new_A->Width();
3981  A = const_cast<HypreParMatrix *>(new_A);
3982  if (precond)
3983  {
3984  precond->SetOperator(*A);
3985  this->SetPreconditioner(*precond);
3986  }
3987  setup_called = 0;
3988  delete X;
3989  delete B;
3990  B = X = NULL;
3991  auxB.Delete(); auxB.Reset();
3992  auxX.Delete(); auxX.Reset();
3993 }
3994 
3995 void HyprePCG::SetTol(double tol)
3996 {
3997  HYPRE_PCGSetTol(pcg_solver, tol);
3998 }
3999 
4000 void HyprePCG::SetAbsTol(double atol)
4001 {
4002  HYPRE_PCGSetAbsoluteTol(pcg_solver, atol);
4003 }
4004 
4005 void HyprePCG::SetMaxIter(int max_iter)
4006 {
4007  HYPRE_PCGSetMaxIter(pcg_solver, max_iter);
4008 }
4009 
4010 void HyprePCG::SetLogging(int logging)
4011 {
4012  HYPRE_PCGSetLogging(pcg_solver, logging);
4013 }
4014 
4015 void HyprePCG::SetPrintLevel(int print_lvl)
4016 {
4017  HYPRE_ParCSRPCGSetPrintLevel(pcg_solver, print_lvl);
4018 }
4019 
4021 {
4022  precond = &precond_;
4023 
4024  HYPRE_ParCSRPCGSetPrecond(pcg_solver,
4025  precond_.SolveFcn(),
4026  precond_.SetupFcn(),
4027  precond_);
4028 }
4029 
4030 void HyprePCG::SetResidualConvergenceOptions(int res_frequency, double rtol)
4031 {
4032  HYPRE_PCGSetTwoNorm(pcg_solver, 1);
4033  if (res_frequency > 0)
4034  {
4035  HYPRE_PCGSetRecomputeResidualP(pcg_solver, res_frequency);
4036  }
4037  if (rtol > 0.0)
4038  {
4039  HYPRE_PCGSetResidualTol(pcg_solver, rtol);
4040  }
4041 }
4042 
4044 {
4045  int myid;
4046  HYPRE_Int time_index = 0;
4047  HYPRE_Int num_iterations;
4048  double final_res_norm;
4049  MPI_Comm comm;
4050  HYPRE_Int print_level;
4051 
4052  HYPRE_PCGGetPrintLevel(pcg_solver, &print_level);
4053  HYPRE_ParCSRPCGSetPrintLevel(pcg_solver, print_level%3);
4054 
4055  HYPRE_ParCSRMatrixGetComm(*A, &comm);
4056 
4057  if (!iterative_mode)
4058  {
4059  x.HypreWrite();
4060  hypre_ParVectorSetConstantValues(x, 0.0);
4061  }
4062 
4063  b.HypreRead();
4064  x.HypreReadWrite();
4065 
4066  if (!setup_called)
4067  {
4068  if (print_level > 0 && print_level < 3)
4069  {
4070  time_index = hypre_InitializeTiming("PCG Setup");
4071  hypre_BeginTiming(time_index);
4072  }
4073 
4074  HYPRE_ParCSRPCGSetup(pcg_solver, *A, b, x);
4075  setup_called = 1;
4076 
4077  if (print_level > 0 && print_level < 3)
4078  {
4079  hypre_EndTiming(time_index);
4080  hypre_PrintTiming("Setup phase times", comm);
4081  hypre_FinalizeTiming(time_index);
4082  hypre_ClearTiming();
4083  }
4084  }
4085 
4086  if (print_level > 0 && print_level < 3)
4087  {
4088  time_index = hypre_InitializeTiming("PCG Solve");
4089  hypre_BeginTiming(time_index);
4090  }
4091 
4092  HYPRE_ParCSRPCGSolve(pcg_solver, *A, b, x);
4093 
4094  if (print_level > 0)
4095  {
4096  if (print_level < 3)
4097  {
4098  hypre_EndTiming(time_index);
4099  hypre_PrintTiming("Solve phase times", comm);
4100  hypre_FinalizeTiming(time_index);
4101  hypre_ClearTiming();
4102  }
4103 
4104  HYPRE_ParCSRPCGGetNumIterations(pcg_solver, &num_iterations);
4105  HYPRE_ParCSRPCGGetFinalRelativeResidualNorm(pcg_solver,
4106  &final_res_norm);
4107 
4108  MPI_Comm_rank(comm, &myid);
4109 
4110  if (myid == 0)
4111  {
4112  mfem::out << "PCG Iterations = " << num_iterations << endl
4113  << "Final PCG Relative Residual Norm = " << final_res_norm
4114  << endl;
4115  }
4116  }
4117  HYPRE_ParCSRPCGSetPrintLevel(pcg_solver, print_level);
4118 }
4119 
4121 {
4122  HYPRE_ParCSRPCGDestroy(pcg_solver);
4123 }
4124 
4125 
4126 HypreGMRES::HypreGMRES(MPI_Comm comm) : precond(NULL)
4127 {
4128  iterative_mode = true;
4129 
4130  HYPRE_ParCSRGMRESCreate(comm, &gmres_solver);
4131  SetDefaultOptions();
4132 }
4133 
4135  : HypreSolver(&A_), precond(NULL)
4136 {
4137  MPI_Comm comm;
4138 
4139  iterative_mode = true;
4140 
4141  HYPRE_ParCSRMatrixGetComm(*A, &comm);
4142 
4143  HYPRE_ParCSRGMRESCreate(comm, &gmres_solver);
4144  SetDefaultOptions();
4145 }
4146 
4147 void HypreGMRES::SetDefaultOptions()
4148 {
4149  int k_dim = 50;
4150  int max_iter = 100;
4151  double tol = 1e-6;
4152 
4153  HYPRE_ParCSRGMRESSetKDim(gmres_solver, k_dim);
4154  HYPRE_ParCSRGMRESSetMaxIter(gmres_solver, max_iter);
4155  HYPRE_ParCSRGMRESSetTol(gmres_solver, tol);
4156 }
4157 
4159 {
4160  const HypreParMatrix *new_A = dynamic_cast<const HypreParMatrix *>(&op);
4161  MFEM_VERIFY(new_A, "new Operator must be a HypreParMatrix!");
4162 
4163  // update base classes: Operator, Solver, HypreSolver
4164  height = new_A->Height();
4165  width = new_A->Width();
4166  A = const_cast<HypreParMatrix *>(new_A);
4167  if (precond)
4168  {
4169  precond->SetOperator(*A);
4170  this->SetPreconditioner(*precond);
4171  }
4172  setup_called = 0;
4173  delete X;
4174  delete B;
4175  B = X = NULL;
4176  auxB.Delete(); auxB.Reset();
4177  auxX.Delete(); auxX.Reset();
4178 }
4179 
4180 void HypreGMRES::SetTol(double tol)
4181 {
4182  HYPRE_GMRESSetTol(gmres_solver, tol);
4183 }
4184 
4185 void HypreGMRES::SetAbsTol(double tol)
4186 {
4187  HYPRE_GMRESSetAbsoluteTol(gmres_solver, tol);
4188 }
4189 
4190 void HypreGMRES::SetMaxIter(int max_iter)
4191 {
4192  HYPRE_GMRESSetMaxIter(gmres_solver, max_iter);
4193 }
4194 
4195 void HypreGMRES::SetKDim(int k_dim)
4196 {
4197  HYPRE_GMRESSetKDim(gmres_solver, k_dim);
4198 }
4199 
4200 void HypreGMRES::SetLogging(int logging)
4201 {
4202  HYPRE_GMRESSetLogging(gmres_solver, logging);
4203 }
4204 
4205 void HypreGMRES::SetPrintLevel(int print_lvl)
4206 {
4207  HYPRE_GMRESSetPrintLevel(gmres_solver, print_lvl);
4208 }
4209 
4211 {
4212  precond = &precond_;
4213 
4214  HYPRE_ParCSRGMRESSetPrecond(gmres_solver,
4215  precond_.SolveFcn(),
4216  precond_.SetupFcn(),
4217  precond_);
4218 }
4219 
4221 {
4222  int myid;
4223  HYPRE_Int time_index = 0;
4224  HYPRE_Int num_iterations;
4225  double final_res_norm;
4226  MPI_Comm comm;
4227  HYPRE_Int print_level;
4228 
4229  HYPRE_GMRESGetPrintLevel(gmres_solver, &print_level);
4230 
4231  HYPRE_ParCSRMatrixGetComm(*A, &comm);
4232 
4233  if (!iterative_mode)
4234  {
4235  x.HypreWrite();
4236  hypre_ParVectorSetConstantValues(x, 0.0);
4237  }
4238 
4239  b.HypreRead();
4240  x.HypreReadWrite();
4241 
4242  if (!setup_called)
4243  {
4244  if (print_level > 0)
4245  {
4246  time_index = hypre_InitializeTiming("GMRES Setup");
4247  hypre_BeginTiming(time_index);
4248  }
4249 
4250  HYPRE_ParCSRGMRESSetup(gmres_solver, *A, b, x);
4251  setup_called = 1;
4252 
4253  if (print_level > 0)
4254  {
4255  hypre_EndTiming(time_index);
4256  hypre_PrintTiming("Setup phase times", comm);
4257  hypre_FinalizeTiming(time_index);
4258  hypre_ClearTiming();
4259  }
4260  }
4261 
4262  if (print_level > 0)
4263  {
4264  time_index = hypre_InitializeTiming("GMRES Solve");
4265  hypre_BeginTiming(time_index);
4266  }
4267 
4268  HYPRE_ParCSRGMRESSolve(gmres_solver, *A, b, x);
4269 
4270  if (print_level > 0)
4271  {
4272  hypre_EndTiming(time_index);
4273  hypre_PrintTiming("Solve phase times", comm);
4274  hypre_FinalizeTiming(time_index);
4275  hypre_ClearTiming();
4276 
4277  HYPRE_ParCSRGMRESGetNumIterations(gmres_solver, &num_iterations);
4278  HYPRE_ParCSRGMRESGetFinalRelativeResidualNorm(gmres_solver,
4279  &final_res_norm);
4280 
4281  MPI_Comm_rank(comm, &myid);
4282 
4283  if (myid == 0)
4284  {
4285  mfem::out << "GMRES Iterations = " << num_iterations << endl
4286  << "Final GMRES Relative Residual Norm = " << final_res_norm
4287  << endl;
4288  }
4289  }
4290 }
4291 
4293 {
4294  HYPRE_ParCSRGMRESDestroy(gmres_solver);
4295 }
4296 
4297 
4298 HypreFGMRES::HypreFGMRES(MPI_Comm comm) : precond(NULL)
4299 {
4300  iterative_mode = true;
4301 
4302  HYPRE_ParCSRFlexGMRESCreate(comm, &fgmres_solver);
4303  SetDefaultOptions();
4304 }
4305 
4307  : HypreSolver(&A_), precond(NULL)
4308 {
4309  MPI_Comm comm;
4310 
4311  iterative_mode = true;
4312 
4313  HYPRE_ParCSRMatrixGetComm(*A, &comm);
4314 
4315  HYPRE_ParCSRFlexGMRESCreate(comm, &fgmres_solver);
4316  SetDefaultOptions();
4317 }
4318 
4319 void HypreFGMRES::SetDefaultOptions()
4320 {
4321  int k_dim = 50;
4322  int max_iter = 100;
4323  double tol = 1e-6;
4324 
4325  HYPRE_ParCSRFlexGMRESSetKDim(fgmres_solver, k_dim);
4326  HYPRE_ParCSRFlexGMRESSetMaxIter(fgmres_solver, max_iter);
4327  HYPRE_ParCSRFlexGMRESSetTol(fgmres_solver, tol);
4328 }
4329 
4331 {
4332  const HypreParMatrix *new_A = dynamic_cast<const HypreParMatrix *>(&op);
4333  MFEM_VERIFY(new_A, "new Operator must be a HypreParMatrix!");
4334 
4335  // update base classes: Operator, Solver, HypreSolver
4336  height = new_A->Height();
4337  width = new_A->Width();
4338  A = const_cast<HypreParMatrix *>(new_A);
4339  if (precond)
4340  {
4341  precond->SetOperator(*A);
4342  this->SetPreconditioner(*precond);
4343  }
4344  setup_called = 0;
4345  delete X;
4346  delete B;
4347  B = X = NULL;
4348  auxB.Delete(); auxB.Reset();
4349  auxX.Delete(); auxX.Reset();
4350 }
4351 
4352 void HypreFGMRES::SetTol(double tol)
4353 {
4354  HYPRE_ParCSRFlexGMRESSetTol(fgmres_solver, tol);
4355 }
4356 
4357 void HypreFGMRES::SetMaxIter(int max_iter)
4358 {
4359  HYPRE_ParCSRFlexGMRESSetMaxIter(fgmres_solver, max_iter);
4360 }
4361 
4362 void HypreFGMRES::SetKDim(int k_dim)
4363 {
4364  HYPRE_ParCSRFlexGMRESSetKDim(fgmres_solver, k_dim);
4365 }
4366 
4367 void HypreFGMRES::SetLogging(int logging)
4368 {
4369  HYPRE_ParCSRFlexGMRESSetLogging(fgmres_solver, logging);
4370 }
4371 
4372 void HypreFGMRES::SetPrintLevel(int print_lvl)
4373 {
4374  HYPRE_ParCSRFlexGMRESSetPrintLevel(fgmres_solver, print_lvl);
4375 }
4376 
4378 {
4379  precond = &precond_;
4380  HYPRE_ParCSRFlexGMRESSetPrecond(fgmres_solver,
4381  precond_.SolveFcn(),
4382  precond_.SetupFcn(),
4383  precond_);
4384 }
4385 
4387 {
4388  int myid;
4389  HYPRE_Int time_index = 0;
4390  HYPRE_Int num_iterations;
4391  double final_res_norm;
4392  MPI_Comm comm;
4393  HYPRE_Int print_level;
4394 
4395  HYPRE_FlexGMRESGetPrintLevel(fgmres_solver, &print_level);
4396 
4397  HYPRE_ParCSRMatrixGetComm(*A, &comm);
4398 
4399  if (!iterative_mode)
4400  {
4401  x.HypreWrite();
4402  hypre_ParVectorSetConstantValues(x, 0.0);
4403  }
4404 
4405  b.HypreRead();
4406  x.HypreReadWrite();
4407 
4408  if (!setup_called)
4409  {
4410  if (print_level > 0)
4411  {
4412  time_index = hypre_InitializeTiming("FGMRES Setup");
4413  hypre_BeginTiming(time_index);
4414  }
4415 
4416  HYPRE_ParCSRFlexGMRESSetup(fgmres_solver, *A, b, x);
4417  setup_called = 1;
4418 
4419  if (print_level > 0)
4420  {
4421  hypre_EndTiming(time_index);
4422  hypre_PrintTiming("Setup phase times", comm);
4423  hypre_FinalizeTiming(time_index);
4424  hypre_ClearTiming();
4425  }
4426  }
4427 
4428  if (print_level > 0)
4429  {
4430  time_index = hypre_InitializeTiming("FGMRES Solve");
4431  hypre_BeginTiming(time_index);
4432  }
4433 
4434  HYPRE_ParCSRFlexGMRESSolve(fgmres_solver, *A, b, x);
4435 
4436  if (print_level > 0)
4437  {
4438  hypre_EndTiming(time_index);
4439  hypre_PrintTiming("Solve phase times", comm);
4440  hypre_FinalizeTiming(time_index);
4441  hypre_ClearTiming();
4442 
4443  HYPRE_ParCSRFlexGMRESGetNumIterations(fgmres_solver, &num_iterations);
4444  HYPRE_ParCSRFlexGMRESGetFinalRelativeResidualNorm(fgmres_solver,
4445  &final_res_norm);
4446 
4447  MPI_Comm_rank(comm, &myid);
4448 
4449  if (myid == 0)
4450  {
4451  mfem::out << "FGMRES Iterations = " << num_iterations << endl
4452  << "Final FGMRES Relative Residual Norm = " << final_res_norm
4453  << endl;
4454  }
4455  }
4456 }
4457 
4459 {
4460  HYPRE_ParCSRFlexGMRESDestroy(fgmres_solver);
4461 }
4462 
4463 
4465 {
4466  const HypreParMatrix *new_A = dynamic_cast<const HypreParMatrix *>(&op);
4467  MFEM_VERIFY(new_A, "new Operator must be a HypreParMatrix!");
4468 
4469  // update base classes: Operator, Solver, HypreSolver
4470  height = new_A->Height();
4471  width = new_A->Width();
4472  A = const_cast<HypreParMatrix *>(new_A);
4473  setup_called = 0;
4474  delete X;
4475  delete B;
4476  B = X = NULL;
4477  auxB.Delete(); auxB.Reset();
4478  auxX.Delete(); auxX.Reset();
4479 }
4480 
4481 
4483 {
4484  HYPRE_ParaSailsCreate(comm, &sai_precond);
4485  SetDefaultOptions();
4486 }
4487 
4489 {
4490  MPI_Comm comm;
4491 
4492  HYPRE_ParCSRMatrixGetComm(A, &comm);
4493 
4494  HYPRE_ParaSailsCreate(comm, &sai_precond);
4495  SetDefaultOptions();
4496 }
4497 
4498 void HypreParaSails::SetDefaultOptions()
4499 {
4500  int sai_max_levels = 1;
4501  double sai_threshold = 0.1;
4502  double sai_filter = 0.1;
4503  int sai_sym = 0;
4504  double sai_loadbal = 0.0;
4505  int sai_reuse = 0;
4506  int sai_logging = 1;
4507 
4508  HYPRE_ParaSailsSetParams(sai_precond, sai_threshold, sai_max_levels);
4509  HYPRE_ParaSailsSetFilter(sai_precond, sai_filter);
4510  HYPRE_ParaSailsSetSym(sai_precond, sai_sym);
4511  HYPRE_ParaSailsSetLoadbal(sai_precond, sai_loadbal);
4512  HYPRE_ParaSailsSetReuse(sai_precond, sai_reuse);
4513  HYPRE_ParaSailsSetLogging(sai_precond, sai_logging);
4514 }
4515 
4516 void HypreParaSails::ResetSAIPrecond(MPI_Comm comm)
4517 {
4518  HYPRE_Int sai_max_levels;
4519  HYPRE_Real sai_threshold;
4520  HYPRE_Real sai_filter;
4521  HYPRE_Int sai_sym;
4522  HYPRE_Real sai_loadbal;
4523  HYPRE_Int sai_reuse;
4524  HYPRE_Int sai_logging;
4525 
4526  // hypre_ParAMGData *amg_data = (hypre_ParAMGData *)sai_precond;
4527  HYPRE_ParaSailsGetNlevels(sai_precond, &sai_max_levels);
4528  HYPRE_ParaSailsGetThresh(sai_precond, &sai_threshold);
4529  HYPRE_ParaSailsGetFilter(sai_precond, &sai_filter);
4530  HYPRE_ParaSailsGetSym(sai_precond, &sai_sym);
4531  HYPRE_ParaSailsGetLoadbal(sai_precond, &sai_loadbal);
4532  HYPRE_ParaSailsGetReuse(sai_precond, &sai_reuse);
4533  HYPRE_ParaSailsGetLogging(sai_precond, &sai_logging);
4534 
4535  HYPRE_ParaSailsDestroy(sai_precond);
4536  HYPRE_ParaSailsCreate(comm, &sai_precond);
4537 
4538  HYPRE_ParaSailsSetParams(sai_precond, sai_threshold, sai_max_levels);
4539  HYPRE_ParaSailsSetFilter(sai_precond, sai_filter);
4540  HYPRE_ParaSailsSetSym(sai_precond, sai_sym);
4541  HYPRE_ParaSailsSetLoadbal(sai_precond, sai_loadbal);
4542  HYPRE_ParaSailsSetReuse(sai_precond, sai_reuse);
4543  HYPRE_ParaSailsSetLogging(sai_precond, sai_logging);
4544 }
4545 
4547 {
4548  const HypreParMatrix *new_A = dynamic_cast<const HypreParMatrix *>(&op);
4549  MFEM_VERIFY(new_A, "new Operator must be a HypreParMatrix!");
4550 
4551  if (A)
4552  {
4553  MPI_Comm comm;
4554  HYPRE_ParCSRMatrixGetComm(*A, &comm);
4555  ResetSAIPrecond(comm);
4556  }
4557 
4558  // update base classes: Operator, Solver, HypreSolver
4559  height = new_A->Height();
4560  width = new_A->Width();
4561  A = const_cast<HypreParMatrix *>(new_A);
4562  setup_called = 0;
4563  delete X;
4564  delete B;
4565  B = X = NULL;
4566  auxB.Delete(); auxB.Reset();
4567  auxX.Delete(); auxX.Reset();
4568 }
4569 
4570 void HypreParaSails::SetParams(double threshold, int max_levels)
4571 {
4572  HYPRE_ParaSailsSetParams(sai_precond, threshold, max_levels);
4573 }
4574 
4575 void HypreParaSails::SetFilter(double filter)
4576 {
4577  HYPRE_ParaSailsSetFilter(sai_precond, filter);
4578 }
4579 
4580 void HypreParaSails::SetLoadBal(double loadbal)
4581 {
4582  HYPRE_ParaSailsSetLoadbal(sai_precond, loadbal);
4583 }
4584 
4586 {
4587  HYPRE_ParaSailsSetReuse(sai_precond, reuse);
4588 }
4589 
4591 {
4592  HYPRE_ParaSailsSetLogging(sai_precond, logging);
4593 }
4594 
4596 {
4597  HYPRE_ParaSailsSetSym(sai_precond, sym);
4598 }
4599 
4601 {
4602  HYPRE_ParaSailsDestroy(sai_precond);
4603 }
4604 
4605 
4607 {
4608  HYPRE_EuclidCreate(comm, &euc_precond);
4609  SetDefaultOptions();
4610 }
4611 
4613 {
4614  MPI_Comm comm;
4615 
4616  HYPRE_ParCSRMatrixGetComm(A, &comm);
4617 
4618  HYPRE_EuclidCreate(comm, &euc_precond);
4619  SetDefaultOptions();
4620 }
4621 
4622 void HypreEuclid::SetDefaultOptions()
4623 {
4624  int euc_level = 1; // We use ILU(1)
4625  int euc_stats = 0; // No logging
4626  int euc_mem = 0; // No memory logging
4627  int euc_bj = 0; // 1: Use Block Jacobi
4628  int euc_ro_sc = 0; // 1: Use Row scaling
4629 
4630  HYPRE_EuclidSetLevel(euc_precond, euc_level);
4631  HYPRE_EuclidSetStats(euc_precond, euc_stats);
4632  HYPRE_EuclidSetMem(euc_precond, euc_mem);
4633  HYPRE_EuclidSetBJ(euc_precond, euc_bj);
4634  HYPRE_EuclidSetRowScale(euc_precond, euc_ro_sc);
4635 }
4636 
4637 void HypreEuclid::SetLevel(int level)
4638 {
4639  HYPRE_EuclidSetLevel(euc_precond, level);
4640 }
4641 
4642 void HypreEuclid::SetStats(int stats)
4643 {
4644  HYPRE_EuclidSetStats(euc_precond, stats);
4645 }
4646 
4648 {
4649  HYPRE_EuclidSetMem(euc_precond, mem);
4650 }
4651 
4653 {
4654  HYPRE_EuclidSetBJ(euc_precond, bj);
4655 }
4656 
4657 void HypreEuclid::SetRowScale(int row_scale)
4658 {
4659  HYPRE_EuclidSetRowScale(euc_precond, row_scale);
4660 }
4661 
4662 void HypreEuclid::ResetEuclidPrecond(MPI_Comm comm)
4663 {
4664  // Euclid does not seem to offer access to its current configuration, so we
4665  // simply reset it to its default options.
4666  HYPRE_EuclidDestroy(euc_precond);
4667  HYPRE_EuclidCreate(comm, &euc_precond);
4668 
4669  SetDefaultOptions();
4670 }
4671 
4673 {
4674  const HypreParMatrix *new_A = dynamic_cast<const HypreParMatrix *>(&op);
4675  MFEM_VERIFY(new_A, "new Operator must be a HypreParMatrix!");
4676 
4677  if (A)
4678  {
4679  MPI_Comm comm;
4680  HYPRE_ParCSRMatrixGetComm(*new_A, &comm);
4681  ResetEuclidPrecond(comm);
4682  }
4683 
4684  // update base classes: Operator, Solver, HypreSolver
4685  height = new_A->Height();
4686  width = new_A->Width();
4687  A = const_cast<HypreParMatrix *>(new_A);
4688  setup_called = 0;
4689  delete X;
4690  delete B;
4691  B = X = NULL;
4692  auxB.Delete(); auxB.Reset();
4693  auxX.Delete(); auxX.Reset();
4694 }
4695 
4697 {
4698  HYPRE_EuclidDestroy(euc_precond);
4699 }
4700 
4701 
4702 #if MFEM_HYPRE_VERSION >= 21900
4704 {
4705  HYPRE_ILUCreate(&ilu_precond);
4706  SetDefaultOptions();
4707 }
4708 
4709 void HypreILU::SetDefaultOptions()
4710 {
4711  // The type of incomplete LU used locally and globally (see class doc)
4712  HYPRE_Int ilu_type = 0; // ILU(k) locally and block Jacobi globally
4713  HYPRE_ILUSetType(ilu_precond, ilu_type);
4714 
4715  // Maximum iterations; 1 iter for preconditioning
4716  HYPRE_Int max_iter = 1;
4717  HYPRE_ILUSetMaxIter(ilu_precond, max_iter);
4718 
4719  // The tolerance when used as a smoother; set to 0.0 for preconditioner
4720  HYPRE_Real tol = 0.0;
4721  HYPRE_ILUSetTol(ilu_precond, tol);
4722 
4723  // Fill level for ILU(k)
4724  HYPRE_Int lev_fill = 1;
4725  HYPRE_ILUSetLevelOfFill(ilu_precond, lev_fill);
4726 
4727  // Local reordering scheme; 0 = no reordering, 1 = reverse Cuthill-McKee
4728  HYPRE_Int reorder_type = 1;
4729  HYPRE_ILUSetLocalReordering(ilu_precond, reorder_type);
4730 
4731  // Information print level; 0 = none, 1 = setup, 2 = solve, 3 = setup+solve
4732  HYPRE_Int print_level = 0;
4733  HYPRE_ILUSetPrintLevel(ilu_precond, print_level);
4734 }
4735 
4736 void HypreILU::ResetILUPrecond()
4737 {
4738  if (ilu_precond)
4739  {
4740  HYPRE_ILUDestroy(ilu_precond);
4741  }
4742  HYPRE_ILUCreate(&ilu_precond);
4743  SetDefaultOptions();
4744 }
4745 
4746 void HypreILU::SetLevelOfFill(HYPRE_Int lev_fill)
4747 {
4748  HYPRE_ILUSetLevelOfFill(ilu_precond, lev_fill);
4749 }
4750 
4751 void HypreILU::SetType(HYPRE_Int ilu_type)
4752 {
4753  HYPRE_ILUSetType(ilu_precond, ilu_type);
4754 }
4755 
4756 void HypreILU::SetMaxIter(HYPRE_Int max_iter)
4757 {
4758  HYPRE_ILUSetMaxIter(ilu_precond, max_iter);
4759 }
4760 
4761 void HypreILU::SetTol(HYPRE_Real tol)
4762 {
4763  HYPRE_ILUSetTol(ilu_precond, tol);
4764 }
4765 
4766 void HypreILU::SetLocalReordering(HYPRE_Int reorder_type)
4767 {
4768  HYPRE_ILUSetLocalReordering(ilu_precond, reorder_type);
4769 }
4770 
4771 void HypreILU::SetPrintLevel(HYPRE_Int print_level)
4772 {
4773  HYPRE_ILUSetPrintLevel(ilu_precond, print_level);
4774 }
4775 
4777 {
4778  const HypreParMatrix *new_A = dynamic_cast<const HypreParMatrix *>(&op);
4779  MFEM_VERIFY(new_A, "new Operator must be a HypreParMatrix!");
4780 
4781  if (A) { ResetILUPrecond(); }
4782 
4783  // update base classes: Operator, Solver, HypreSolver
4784  height = new_A->Height();
4785  width = new_A->Width();
4786  A = const_cast<HypreParMatrix *>(new_A);
4787  setup_called = 0;
4788  delete X;
4789  delete B;
4790  B = X = NULL;
4791  auxB.Delete(); auxB.Reset();
4792  auxX.Delete(); auxX.Reset();
4793 }
4794 
4796 {
4797  HYPRE_ILUDestroy(ilu_precond);
4798 }
4799 #endif
4800 
4801 
4803 {
4804  HYPRE_BoomerAMGCreate(&amg_precond);
4805  SetDefaultOptions();
4806 }
4807 
4809 {
4810  HYPRE_BoomerAMGCreate(&amg_precond);
4811  SetDefaultOptions();
4812 }
4813 
4814 void HypreBoomerAMG::SetDefaultOptions()
4815 {
4816 #if !defined(HYPRE_USING_GPU)
4817  // AMG coarsening options:
4818  int coarsen_type = 10; // 10 = HMIS, 8 = PMIS, 6 = Falgout, 0 = CLJP
4819  int agg_levels = 1; // number of aggressive coarsening levels
4820  double theta = 0.25; // strength threshold: 0.25, 0.5, 0.8
4821 
4822  // AMG interpolation options:
4823  int interp_type = 6; // 6 = extended+i, 0 = classical
4824  int Pmax = 4; // max number of elements per row in P
4825 
4826  // AMG relaxation options:
4827  int relax_type = 8; // 8 = l1-GS, 6 = symm. GS, 3 = GS, 18 = l1-Jacobi
4828  int relax_sweeps = 1; // relaxation sweeps on each level
4829 
4830  // Additional options:
4831  int print_level = 1; // print AMG iterations? 1 = no, 2 = yes
4832  int max_levels = 25; // max number of levels in AMG hierarchy
4833 #else
4834  // AMG coarsening options:
4835  int coarsen_type = 8; // 10 = HMIS, 8 = PMIS, 6 = Falgout, 0 = CLJP
4836  int agg_levels = 0; // number of aggressive coarsening levels
4837  double theta = 0.25; // strength threshold: 0.25, 0.5, 0.8
4838 
4839  // AMG interpolation options:
4840  int interp_type = 6; // 6 = extended+i, or 18 = extended+e
4841  int Pmax = 4; // max number of elements per row in P
4842 
4843  // AMG relaxation options:
4844  int relax_type = 18; // 18 = l1-Jacobi, or 16 = Chebyshev
4845  int relax_sweeps = 1; // relaxation sweeps on each level
4846 
4847  // Additional options:
4848  int print_level = 1; // print AMG iterations? 1 = no, 2 = yes
4849  int max_levels = 25; // max number of levels in AMG hierarchy
4850 #endif
4851 
4852  HYPRE_BoomerAMGSetCoarsenType(amg_precond, coarsen_type);
4853  HYPRE_BoomerAMGSetAggNumLevels(amg_precond, agg_levels);
4854  HYPRE_BoomerAMGSetRelaxType(amg_precond, relax_type);
4855  // default in hypre is 1.0 with some exceptions, e.g. for relax_type = 7
4856  // HYPRE_BoomerAMGSetRelaxWt(amg_precond, 1.0);
4857  HYPRE_BoomerAMGSetNumSweeps(amg_precond, relax_sweeps);
4858  HYPRE_BoomerAMGSetStrongThreshold(amg_precond, theta);
4859  HYPRE_BoomerAMGSetInterpType(amg_precond, interp_type);
4860  HYPRE_BoomerAMGSetPMaxElmts(amg_precond, Pmax);
4861  HYPRE_BoomerAMGSetPrintLevel(amg_precond, print_level);
4862  HYPRE_BoomerAMGSetMaxLevels(amg_precond, max_levels);
4863 
4864  // Use as a preconditioner (one V-cycle, zero tolerance)
4865  HYPRE_BoomerAMGSetMaxIter(amg_precond, 1);
4866  HYPRE_BoomerAMGSetTol(amg_precond, 0.0);
4867 }
4868 
4869 void HypreBoomerAMG::ResetAMGPrecond()
4870 {
4871  HYPRE_Int coarsen_type;
4872  HYPRE_Int agg_levels;
4873  HYPRE_Int relax_type;
4874  HYPRE_Int relax_sweeps;
4875  HYPRE_Real theta;
4876  HYPRE_Int interp_type;
4877  HYPRE_Int Pmax;
4878  HYPRE_Int print_level;
4879  HYPRE_Int max_levels;
4880  HYPRE_Int dim;
4881  HYPRE_Int nrbms = rbms.Size();
4882  HYPRE_Int nodal;
4883  HYPRE_Int nodal_diag;
4884  HYPRE_Int relax_coarse;
4885  HYPRE_Int interp_vec_variant;
4886  HYPRE_Int q_max;
4887  HYPRE_Int smooth_interp_vectors;
4888  HYPRE_Int interp_refine;
4889 
4890  hypre_ParAMGData *amg_data = (hypre_ParAMGData *)amg_precond;
4891 
4892  // read options from amg_precond
4893  HYPRE_BoomerAMGGetCoarsenType(amg_precond, &coarsen_type);
4894  agg_levels = hypre_ParAMGDataAggNumLevels(amg_data);
4895  relax_type = hypre_ParAMGDataUserRelaxType(amg_data);
4896  relax_sweeps = hypre_ParAMGDataUserNumSweeps(amg_data);
4897  HYPRE_BoomerAMGGetStrongThreshold(amg_precond, &theta);
4898  hypre_BoomerAMGGetInterpType(amg_precond, &interp_type);
4899  HYPRE_BoomerAMGGetPMaxElmts(amg_precond, &Pmax);
4900  HYPRE_BoomerAMGGetPrintLevel(amg_precond, &print_level);
4901  HYPRE_BoomerAMGGetMaxLevels(amg_precond, &max_levels);
4902  HYPRE_BoomerAMGGetNumFunctions(amg_precond, &dim);
4903  if (nrbms) // elasticity solver options
4904  {
4905  nodal = hypre_ParAMGDataNodal(amg_data);
4906  nodal_diag = hypre_ParAMGDataNodalDiag(amg_data);
4907  HYPRE_BoomerAMGGetCycleRelaxType(amg_precond, &relax_coarse, 3);
4908  interp_vec_variant = hypre_ParAMGInterpVecVariant(amg_data);
4909  q_max = hypre_ParAMGInterpVecQMax(amg_data);
4910  smooth_interp_vectors = hypre_ParAMGSmoothInterpVectors(amg_data);
4911  interp_refine = hypre_ParAMGInterpRefine(amg_data);
4912  }
4913 
4914  HYPRE_BoomerAMGDestroy(amg_precond);
4915  HYPRE_BoomerAMGCreate(&amg_precond);
4916 
4917  HYPRE_BoomerAMGSetCoarsenType(amg_precond, coarsen_type);
4918  HYPRE_BoomerAMGSetAggNumLevels(amg_precond, agg_levels);
4919  HYPRE_BoomerAMGSetRelaxType(amg_precond, relax_type);
4920  HYPRE_BoomerAMGSetNumSweeps(amg_precond, relax_sweeps);
4921  HYPRE_BoomerAMGSetMaxLevels(amg_precond, max_levels);
4922  HYPRE_BoomerAMGSetTol(amg_precond, 0.0);
4923  HYPRE_BoomerAMGSetMaxIter(amg_precond, 1); // one V-cycle
4924  HYPRE_BoomerAMGSetStrongThreshold(amg_precond, theta);
4925  HYPRE_BoomerAMGSetInterpType(amg_precond, interp_type);
4926  HYPRE_BoomerAMGSetPMaxElmts(amg_precond, Pmax);
4927  HYPRE_BoomerAMGSetPrintLevel(amg_precond, print_level);
4928  HYPRE_BoomerAMGSetNumFunctions(amg_precond, dim);
4929  if (nrbms)
4930  {
4931  HYPRE_BoomerAMGSetNodal(amg_precond, nodal);
4932  HYPRE_BoomerAMGSetNodalDiag(amg_precond, nodal_diag);
4933  HYPRE_BoomerAMGSetCycleRelaxType(amg_precond, relax_coarse, 3);
4934  HYPRE_BoomerAMGSetInterpVecVariant(amg_precond, interp_vec_variant);
4935  HYPRE_BoomerAMGSetInterpVecQMax(amg_precond, q_max);
4936  HYPRE_BoomerAMGSetSmoothInterpVectors(amg_precond, smooth_interp_vectors);
4937  HYPRE_BoomerAMGSetInterpRefine(amg_precond, interp_refine);
4938  RecomputeRBMs();
4939  HYPRE_BoomerAMGSetInterpVectors(amg_precond, rbms.Size(), rbms.GetData());
4940  }
4941 }
4942 
4944 {
4945  const HypreParMatrix *new_A = dynamic_cast<const HypreParMatrix *>(&op);
4946  MFEM_VERIFY(new_A, "new Operator must be a HypreParMatrix!");
4947 
4948  if (A) { ResetAMGPrecond(); }
4949 
4950  // update base classes: Operator, Solver, HypreSolver
4951  height = new_A->Height();
4952  width = new_A->Width();
4953  A = const_cast<HypreParMatrix *>(new_A);
4954  setup_called = 0;
4955  delete X;
4956  delete B;
4957  B = X = NULL;
4958  auxB.Delete(); auxB.Reset();
4959  auxX.Delete(); auxX.Reset();
4960 }
4961 
4962 void HypreBoomerAMG::SetSystemsOptions(int dim, bool order_bynodes)
4963 {
4964  HYPRE_BoomerAMGSetNumFunctions(amg_precond, dim);
4965 
4966  // The default "system" ordering in hypre is Ordering::byVDIM. When we are
4967  // using Ordering::byNODES, we have to specify the ordering explicitly with
4968  // HYPRE_BoomerAMGSetDofFunc as in the following code.
4969  if (order_bynodes)
4970  {
4971  // Generate DofFunc mapping on the host
4972  HYPRE_Int *h_mapping = mfem_hypre_CTAlloc_host(HYPRE_Int, height);
4973  int h_nnodes = height / dim; // nodes owned in linear algebra (not fem)
4974  MFEM_VERIFY(height % dim == 0, "Ordering does not work as claimed!");
4975  int k = 0;
4976  for (int i = 0; i < dim; ++i)
4977  {
4978  for (int j = 0; j < h_nnodes; ++j)
4979  {
4980  h_mapping[k++] = i;
4981  }
4982  }
4983 
4984  // After the addition of hypre_IntArray, mapping is assumed
4985  // to be a device pointer. Previously, it was assumed to be
4986  // a host pointer.
4987 #if defined(hypre_IntArrayData) && defined(HYPRE_USING_GPU)
4988  HYPRE_Int *mapping = mfem_hypre_CTAlloc(HYPRE_Int, height);
4989  hypre_TMemcpy(mapping, h_mapping, HYPRE_Int, height,
4990  HYPRE_MEMORY_DEVICE, HYPRE_MEMORY_HOST);
4991  mfem_hypre_TFree_host(h_mapping);
4992 #else
4993  HYPRE_Int *mapping = h_mapping;
4994 #endif
4995 
4996  // hypre actually deletes the mapping pointer in HYPRE_BoomerAMGDestroy,
4997  // so we don't need to track it
4998  HYPRE_BoomerAMGSetDofFunc(amg_precond, mapping);
4999  }
5000 
5001  // More robust options with respect to convergence
5002  HYPRE_BoomerAMGSetAggNumLevels(amg_precond, 0);
5003  HYPRE_BoomerAMGSetStrongThreshold(amg_precond, 0.5);
5004 }
5005 
5006 // Rotational rigid-body mode functions, used in SetElasticityOptions()
5007 static void func_rxy(const Vector &x, Vector &y)
5008 {
5009  y = 0.0; y(0) = x(1); y(1) = -x(0);
5010 }
5011 static void func_ryz(const Vector &x, Vector &y)
5012 {
5013  y = 0.0; y(1) = x(2); y(2) = -x(1);
5014 }
5015 static void func_rzx(const Vector &x, Vector &y)
5016 {
5017  y = 0.0; y(2) = x(0); y(0) = -x(2);
5018 }
5019 
5020 void HypreBoomerAMG::RecomputeRBMs()
5021 {
5022  int nrbms;
5023  Array<HypreParVector*> gf_rbms;
5024  int dim = fespace->GetParMesh()->Dimension();
5025 
5026  for (int i = 0; i < rbms.Size(); i++)
5027  {
5028  HYPRE_ParVectorDestroy(rbms[i]);
5029  }
5030 
5031  if (dim == 2)
5032  {
5033  nrbms = 1;
5034 
5035  VectorFunctionCoefficient coeff_rxy(2, func_rxy);
5036 
5037  ParGridFunction rbms_rxy(fespace);
5038  rbms_rxy.ProjectCoefficient(coeff_rxy);
5039 
5040  rbms.SetSize(nrbms);
5041  gf_rbms.SetSize(nrbms);
5042  gf_rbms[0] = fespace->NewTrueDofVector();
5043  rbms_rxy.GetTrueDofs(*gf_rbms[0]);
5044  }
5045  else if (dim == 3)
5046  {
5047  nrbms = 3;
5048 
5049  VectorFunctionCoefficient coeff_rxy(3, func_rxy);
5050  VectorFunctionCoefficient coeff_ryz(3, func_ryz);
5051  VectorFunctionCoefficient coeff_rzx(3, func_rzx);
5052 
5053  ParGridFunction rbms_rxy(fespace);
5054  ParGridFunction rbms_ryz(fespace);
5055  ParGridFunction rbms_rzx(fespace);
5056  rbms_rxy.ProjectCoefficient(coeff_rxy);
5057  rbms_ryz.ProjectCoefficient(coeff_ryz);
5058  rbms_rzx.ProjectCoefficient(coeff_rzx);
5059 
5060  rbms.SetSize(nrbms);
5061  gf_rbms.SetSize(nrbms);
5062  gf_rbms[0] = fespace->NewTrueDofVector();
5063  gf_rbms[1] = fespace->NewTrueDofVector();
5064  gf_rbms[2] = fespace->NewTrueDofVector();
5065  rbms_rxy.GetTrueDofs(*gf_rbms[0]);
5066  rbms_ryz.GetTrueDofs(*gf_rbms[1]);
5067  rbms_rzx.GetTrueDofs(*gf_rbms[2]);
5068  }
5069  else
5070  {
5071  nrbms = 0;
5072  rbms.SetSize(nrbms);
5073  }
5074 
5075  // Transfer the RBMs from the ParGridFunction to the HYPRE_ParVector objects
5076  for (int i = 0; i < nrbms; i++)
5077  {
5078  rbms[i] = gf_rbms[i]->StealParVector();
5079  delete gf_rbms[i];
5080  }
5081 }
5082 
5084 {
5085 #ifdef HYPRE_USING_GPU
5086  MFEM_ABORT("this method is not supported in hypre built with GPU support");
5087 #endif
5088 
5089  // Save the finite element space to support multiple calls to SetOperator()
5090  this->fespace = fespace_;
5091 
5092  // Make sure the systems AMG options are set
5093  int dim = fespace_->GetParMesh()->Dimension();
5094  SetSystemsOptions(dim);
5095 
5096  // Nodal coarsening options (nodal coarsening is required for this solver)
5097  // See hypre's new_ij driver and the paper for descriptions.
5098  int nodal = 4; // strength reduction norm: 1, 3 or 4
5099  int nodal_diag = 1; // diagonal in strength matrix: 0, 1 or 2
5100  int relax_coarse = 8; // smoother on the coarsest grid: 8, 99 or 29
5101 
5102  // Elasticity interpolation options
5103  int interp_vec_variant = 2; // 1 = GM-1, 2 = GM-2, 3 = LN
5104  int q_max = 4; // max elements per row for each Q
5105  int smooth_interp_vectors = 1; // smooth the rigid-body modes?
5106 
5107  // Optionally pre-process the interpolation matrix through iterative weight
5108  // refinement (this is generally applicable for any system)
5109  int interp_refine = 1;
5110 
5111  HYPRE_BoomerAMGSetNodal(amg_precond, nodal);
5112  HYPRE_BoomerAMGSetNodalDiag(amg_precond, nodal_diag);
5113  HYPRE_BoomerAMGSetCycleRelaxType(amg_precond, relax_coarse, 3);
5114  HYPRE_BoomerAMGSetInterpVecVariant(amg_precond, interp_vec_variant);
5115  HYPRE_BoomerAMGSetInterpVecQMax(amg_precond, q_max);
5116  HYPRE_BoomerAMGSetSmoothInterpVectors(amg_precond, smooth_interp_vectors);
5117  HYPRE_BoomerAMGSetInterpRefine(amg_precond, interp_refine);
5118 
5119  RecomputeRBMs();
5120  HYPRE_BoomerAMGSetInterpVectors(amg_precond, rbms.Size(), rbms.GetData());
5121 
5122  // The above BoomerAMG options may result in singular matrices on the coarse
5123  // grids, which are handled correctly in hypre's Solve method, but can produce
5124  // hypre errors in the Setup (specifically in the l1 row norm computation).
5125  // See the documentation of SetErrorMode() for more details.
5127 }
5128 
5129 #if MFEM_HYPRE_VERSION >= 21800
5130 
5132  const std::string &prerelax,
5133  const std::string &postrelax)
5134 {
5135  // Hypre parameters
5136  int Sabs = 0;
5137  int interp_type = 100;
5138  int relax_type = 10;
5139  int coarsen_type = 6;
5140  double strength_tolC = 0.1;
5141  double strength_tolR = 0.01;
5142  double filter_tolR = 0.0;
5143  double filterA_tol = 0.0;
5144 
5145  // Set relaxation on specified grid points
5146  int ns_down, ns_up, ns_coarse;
5147  if (distanceR > 0)
5148  {
5149  ns_down = prerelax.length();
5150  ns_up = postrelax.length();
5151  ns_coarse = 1;
5152 
5153  // Array to store relaxation scheme and pass to Hypre
5154  HYPRE_Int **grid_relax_points = mfem_hypre_TAlloc(HYPRE_Int*, 4);
5155  grid_relax_points[0] = NULL;
5156  grid_relax_points[1] = mfem_hypre_TAlloc(HYPRE_Int, ns_down);
5157  grid_relax_points[2] = mfem_hypre_TAlloc(HYPRE_Int, ns_up);
5158  grid_relax_points[3] = mfem_hypre_TAlloc(HYPRE_Int, 1);
5159  grid_relax_points[3][0] = 0;
5160 
5161  // set down relax scheme
5162  for (int i = 0; i<ns_down; i++)
5163  {
5164  if (prerelax[i] == 'F')
5165  {
5166  grid_relax_points[1][i] = -1;
5167  }
5168  else if (prerelax[i] == 'C')
5169  {
5170  grid_relax_points[1][i] = 1;
5171  }
5172  else if (prerelax[i] == 'A')
5173  {
5174  grid_relax_points[1][i] = 0;
5175  }
5176  }
5177 
5178  // set up relax scheme
5179  for (int i = 0; i<ns_up; i++)
5180  {
5181  if (postrelax[i] == 'F')
5182  {
5183  grid_relax_points[2][i] = -1;
5184  }
5185  else if (postrelax[i] == 'C')
5186  {
5187  grid_relax_points[2][i] = 1;
5188  }
5189  else if (postrelax[i] == 'A')
5190  {
5191  grid_relax_points[2][i] = 0;
5192  }
5193  }
5194 
5195  HYPRE_BoomerAMGSetRestriction(amg_precond, distanceR);
5196 
5197  HYPRE_BoomerAMGSetGridRelaxPoints(amg_precond, grid_relax_points);
5198 
5199  HYPRE_BoomerAMGSetInterpType(amg_precond, interp_type);
5200  }
5201 
5202  if (Sabs)
5203  {
5204  HYPRE_BoomerAMGSetSabs(amg_precond, Sabs);
5205  }
5206 
5207  HYPRE_BoomerAMGSetCoarsenType(amg_precond, coarsen_type);
5208 
5209  // does not support aggressive coarsening
5210  HYPRE_BoomerAMGSetAggNumLevels(amg_precond, 0);
5211 
5212  HYPRE_BoomerAMGSetStrongThreshold(amg_precond, strength_tolC);
5213 
5214  if (distanceR > 0)
5215  {
5216  HYPRE_BoomerAMGSetStrongThresholdR(amg_precond, strength_tolR);
5217  HYPRE_BoomerAMGSetFilterThresholdR(amg_precond, filter_tolR);
5218  }
5219 
5220  if (relax_type > -1)
5221  {
5222  HYPRE_BoomerAMGSetRelaxType(amg_precond, relax_type);
5223  }
5224 
5225  if (distanceR > 0)
5226  {
5227  HYPRE_BoomerAMGSetCycleNumSweeps(amg_precond, ns_coarse, 3);
5228  HYPRE_BoomerAMGSetCycleNumSweeps(amg_precond, ns_down, 1);
5229  HYPRE_BoomerAMGSetCycleNumSweeps(amg_precond, ns_up, 2);
5230 
5231  HYPRE_BoomerAMGSetADropTol(amg_precond, filterA_tol);
5232  // type = -1: drop based on row inf-norm
5233  HYPRE_BoomerAMGSetADropType(amg_precond, -1);
5234  }
5235 }
5236 
5237 #endif
5238 
5240 {
5241  for (int i = 0; i < rbms.Size(); i++)
5242  {
5243  HYPRE_ParVectorDestroy(rbms[i]);
5244  }
5245 
5246  HYPRE_BoomerAMGDestroy(amg_precond);
5247 }
5248 
5250 {
5251  Init(edge_fespace);
5252 }
5253 
5255  : HypreSolver(&A)
5256 {
5257  Init(edge_fespace);
5258 }
5259 
5262  : HypreSolver(&A),
5263  x(x_),
5264  y(y_),
5265  z(z_),
5266  G(G_),
5267  Pi(NULL),
5268  Pix(NULL),
5269  Piy(NULL),
5270  Piz(NULL)
5271 {
5272  MFEM_ASSERT(G != NULL, "");
5273  MFEM_ASSERT(x != NULL, "");
5274  MFEM_ASSERT(y != NULL, "");
5275  int sdim = (z == NULL) ? 2 : 3;
5276  int cycle_type = 13;
5277 
5278  MakeSolver(sdim, cycle_type);
5279 
5280  HYPRE_ParVector pz = z ? static_cast<HYPRE_ParVector>(*z) : NULL;
5281  HYPRE_AMSSetCoordinateVectors(ams, *x, *y, pz);
5282  HYPRE_AMSSetDiscreteGradient(ams, *G);
5283 }
5284 
5285 void HypreAMS::MakeSolver(int sdim, int cycle_type)
5286 {
5287  int rlx_sweeps = 1;
5288  double rlx_weight = 1.0;
5289  double rlx_omega = 1.0;
5290 #if !defined(HYPRE_USING_GPU)
5291  int amg_coarsen_type = 10;
5292  int amg_agg_levels = 1;
5293  int amg_rlx_type = 8;
5294  int rlx_type = 2;
5295  double theta = 0.25;
5296  int amg_interp_type = 6;
5297  int amg_Pmax = 4;
5298 #else
5299  int amg_coarsen_type = 8;
5300  int amg_agg_levels = 0;
5301  int amg_rlx_type = 18;
5302  int rlx_type = 1;
5303  double theta = 0.25;
5304  int amg_interp_type = 6;
5305  int amg_Pmax = 4;
5306 #endif
5307 
5308  space_dim = sdim;
5309  ams_cycle_type = cycle_type;
5310  HYPRE_AMSCreate(&ams);
5311 
5312  HYPRE_AMSSetDimension(ams, sdim); // 2D H(div) and 3D H(curl) problems
5313  HYPRE_AMSSetTol(ams, 0.0);
5314  HYPRE_AMSSetMaxIter(ams, 1); // use as a preconditioner
5315  HYPRE_AMSSetCycleType(ams, cycle_type);
5316  HYPRE_AMSSetPrintLevel(ams, 1);
5317 
5318  // set additional AMS options
5319  HYPRE_AMSSetSmoothingOptions(ams, rlx_type, rlx_sweeps, rlx_weight, rlx_omega);
5320  HYPRE_AMSSetAlphaAMGOptions(ams, amg_coarsen_type, amg_agg_levels, amg_rlx_type,
5321  theta, amg_interp_type, amg_Pmax);
5322  HYPRE_AMSSetBetaAMGOptions(ams, amg_coarsen_type, amg_agg_levels, amg_rlx_type,
5323  theta, amg_interp_type, amg_Pmax);
5324 
5325  HYPRE_AMSSetAlphaAMGCoarseRelaxType(ams, amg_rlx_type);
5326  HYPRE_AMSSetBetaAMGCoarseRelaxType(ams, amg_rlx_type);
5327 
5328  // The AMS preconditioner may sometimes require inverting singular matrices
5329  // with BoomerAMG, which are handled correctly in hypre's Solve method, but
5330  // can produce hypre errors in the Setup (specifically in the l1 row norm
5331  // computation). See the documentation of SetErrorMode() for more details.
5333 }
5334 
5335 void HypreAMS::MakeGradientAndInterpolation(
5336  ParFiniteElementSpace *edge_fespace, int cycle_type)
5337 {
5338  int dim = edge_fespace->GetMesh()->Dimension();
5339  int sdim = edge_fespace->GetMesh()->SpaceDimension();
5340  const FiniteElementCollection *edge_fec = edge_fespace->FEColl();
5341 
5342  bool trace_space, rt_trace_space;
5343  ND_Trace_FECollection *nd_tr_fec = NULL;
5344  trace_space = dynamic_cast<const ND_Trace_FECollection*>(edge_fec);
5345  rt_trace_space = dynamic_cast<const RT_Trace_FECollection*>(edge_fec);
5346  trace_space = trace_space || rt_trace_space;
5347 
5348  int p = 1;
5349  if (edge_fespace->GetNE() > 0)
5350  {
5351  MFEM_VERIFY(!edge_fespace->IsVariableOrder(), "");
5352  if (trace_space)
5353  {
5354  p = edge_fespace->GetFaceOrder(0);
5355  if (dim == 2) { p++; }
5356  }
5357  else
5358  {
5359  p = edge_fespace->GetElementOrder(0);
5360  }
5361  }
5362 
5363  ParMesh *pmesh = edge_fespace->GetParMesh();
5364  if (rt_trace_space)
5365  {
5366  nd_tr_fec = new ND_Trace_FECollection(p, dim);
5367  edge_fespace = new ParFiniteElementSpace(pmesh, nd_tr_fec);
5368  }
5369 
5370  // define the nodal linear finite element space associated with edge_fespace
5371  FiniteElementCollection *vert_fec;
5372  if (trace_space)
5373  {
5374  vert_fec = new H1_Trace_FECollection(p, dim);
5375  }
5376  else
5377  {
5378  vert_fec = new H1_FECollection(p, dim);
5379  }
5380  ParFiniteElementSpace *vert_fespace = new ParFiniteElementSpace(pmesh,
5381  vert_fec);
5382 
5383  // generate and set the vertex coordinates
5384  if (p == 1 && pmesh->GetNodes() == NULL)
5385  {
5386  ParGridFunction x_coord(vert_fespace);
5387  ParGridFunction y_coord(vert_fespace);
5388  ParGridFunction z_coord(vert_fespace);
5389  double *coord;
5390  for (int i = 0; i < pmesh->GetNV(); i++)
5391  {
5392  coord = pmesh -> GetVertex(i);
5393  x_coord(i) = coord[0];
5394  if (sdim >= 2) { y_coord(i) = coord[1]; }
5395  if (sdim == 3) { z_coord(i) = coord[2]; }
5396  }
5397  x = x_coord.ParallelProject();
5398  y = NULL;
5399  z = NULL;
5400  x->HypreReadWrite();
5401 
5402  if (sdim >= 2)
5403  {
5404  y = y_coord.ParallelProject();
5405  y->HypreReadWrite();
5406  }
5407  if (sdim == 3)
5408  {
5409  z = z_coord.ParallelProject();
5410  z->HypreReadWrite();
5411  }
5412 
5413  HYPRE_AMSSetCoordinateVectors(ams,
5414  x ? (HYPRE_ParVector)(*x) : NULL,
5415  y ? (HYPRE_ParVector)(*y) : NULL,
5416  z ? (HYPRE_ParVector)(*z) : NULL);
5417  }
5418  else
5419  {
5420  x = NULL;
5421  y = NULL;
5422  z = NULL;
5423  }
5424 
5425  // generate and set the discrete gradient
5426  ParDiscreteLinearOperator *grad;
5427  grad = new ParDiscreteLinearOperator(vert_fespace, edge_fespace);
5428  if (trace_space)
5429  {
5430  grad->AddTraceFaceInterpolator(new GradientInterpolator);
5431  }
5432  else
5433  {
5434  grad->AddDomainInterpolator(new GradientInterpolator);
5435  }
5436  grad->Assemble();
5437  grad->Finalize();
5438  G = grad->ParallelAssemble();
5439  HYPRE_AMSSetDiscreteGradient(ams, *G);
5440  delete grad;
5441 
5442  // generate and set the Nedelec interpolation matrices
5443  Pi = Pix = Piy = Piz = NULL;
5444  if (p > 1 || pmesh->GetNodes() != NULL)
5445  {
5446  ParFiniteElementSpace *vert_fespace_d
5447  = new ParFiniteElementSpace(pmesh, vert_fec, sdim, Ordering::byVDIM);
5448 
5449  ParDiscreteLinearOperator *id_ND;
5450  id_ND = new ParDiscreteLinearOperator(vert_fespace_d, edge_fespace);
5451  if (trace_space)
5452  {
5453  id_ND->AddTraceFaceInterpolator(new IdentityInterpolator);
5454  }
5455  else
5456  {
5457  id_ND->AddDomainInterpolator(new IdentityInterpolator);
5458  }
5459  id_ND->Assemble();
5460  id_ND->Finalize();
5461 
5462  if (cycle_type < 10)
5463  {
5464  Pi = id_ND->ParallelAssemble();
5465  }
5466  else
5467  {
5468  Array2D<HypreParMatrix *> Pi_blocks;
5469  id_ND->GetParBlocks(Pi_blocks);
5470  Pix = Pi_blocks(0,0);
5471  if (sdim >= 2) { Piy = Pi_blocks(0,1); }
5472  if (sdim == 3) { Piz = Pi_blocks(0,2); }
5473  }
5474 
5475  delete id_ND;
5476 
5477  HYPRE_ParCSRMatrix HY_Pi = (Pi) ? (HYPRE_ParCSRMatrix) *Pi : NULL;
5478  HYPRE_ParCSRMatrix HY_Pix = (Pix) ? (HYPRE_ParCSRMatrix) *Pix : NULL;
5479  HYPRE_ParCSRMatrix HY_Piy = (Piy) ? (HYPRE_ParCSRMatrix) *Piy : NULL;
5480  HYPRE_ParCSRMatrix HY_Piz = (Piz) ? (HYPRE_ParCSRMatrix) *Piz : NULL;
5481  HYPRE_AMSSetInterpolations(ams, HY_Pi, HY_Pix, HY_Piy, HY_Piz);
5482 
5483  delete vert_fespace_d;
5484  }
5485 
5486  delete vert_fespace;
5487  delete vert_fec;
5488 
5489  if (rt_trace_space)
5490  {
5491  delete edge_fespace;
5492  delete nd_tr_fec;
5493  }
5494 }
5495 
5496 void HypreAMS::Init(ParFiniteElementSpace *edge_fespace)
5497 {
5498  int cycle_type = 13;
5499  int sdim = edge_fespace->GetMesh()->SpaceDimension();
5500  MakeSolver(sdim, cycle_type);
5501  MakeGradientAndInterpolation(edge_fespace, cycle_type);
5502 }
5503 
5504 void HypreAMS::ResetAMSPrecond()
5505 {
5506 #if MFEM_HYPRE_VERSION >= 22600
5507  /* Read options from ams */
5508  auto *ams_data = (hypre_AMSData *)ams;
5509 
5510  /* Space dimension */
5511  HYPRE_Int dim = hypre_AMSDataDimension(ams_data);
5512 
5513  /* Vertex space data */
5514  hypre_ParCSRMatrix *hy_G = hypre_AMSDataDiscreteGradient(ams_data);
5515 
5516  HYPRE_Int beta_is_zero = hypre_AMSDataBetaIsZero(ams_data);
5517 
5518  /* Vector vertex space data */
5519  hypre_ParCSRMatrix *hy_Pi hypre_AMSDataPiInterpolation(ams_data);
5520  hypre_ParCSRMatrix *hy_Pix = ams_data->Pix;
5521  hypre_ParCSRMatrix *hy_Piy = ams_data->Piy;
5522  hypre_ParCSRMatrix *hy_Piz = ams_data->Piz;
5523  HYPRE_Int owns_Pi = hypre_AMSDataOwnsPiInterpolation(ams_data);
5524  if (owns_Pi)
5525  {
5526  ams_data->owns_Pi = 0; // we're stealing Pi
5527  }
5528 
5529  /* Coordinates of the vertices */
5530  hypre_ParVector *hy_x = hypre_AMSDataVertexCoordinateX(ams_data);
5531  hypre_ParVector *hy_y = hypre_AMSDataVertexCoordinateY(ams_data);
5532  hypre_ParVector *hy_z = hypre_AMSDataVertexCoordinateZ(ams_data);
5533 
5534  /* Solver options */
5535  HYPRE_Int maxit = hypre_AMSDataMaxIter(ams_data);
5536  HYPRE_Real tol = hypre_AMSDataTol(ams_data);
5537  HYPRE_Int cycle_type = hypre_AMSDataCycleType(ams_data);
5538  HYPRE_Int ams_print_level = hypre_AMSDataPrintLevel(ams_data);
5539 
5540  /* Smoothing and AMG options */
5541  HYPRE_Int A_relax_type = hypre_AMSDataARelaxType(ams_data);
5542  HYPRE_Int A_relax_times = hypre_AMSDataARelaxTimes(ams_data);
5543  HYPRE_Real A_relax_weight = hypre_AMSDataARelaxWeight(ams_data);
5544  HYPRE_Real A_omega = hypre_AMSDataAOmega(ams_data);
5545  HYPRE_Int A_cheby_order = hypre_AMSDataAChebyOrder(ams_data);
5546  HYPRE_Real A_cheby_fraction = hypre_AMSDataAChebyFraction(ams_data);
5547 
5548  HYPRE_Int B_Pi_coarsen_type = hypre_AMSDataPoissonAlphaAMGCoarsenType(ams_data);
5549  HYPRE_Int B_Pi_agg_levels = hypre_AMSDataPoissonAlphaAMGAggLevels(ams_data);
5550  HYPRE_Int B_Pi_relax_type = hypre_AMSDataPoissonAlphaAMGRelaxType(ams_data);
5551  HYPRE_Int B_Pi_coarse_relax_type = ams_data->B_Pi_coarse_relax_type;
5552  HYPRE_Real B_Pi_theta = hypre_AMSDataPoissonAlphaAMGStrengthThreshold(ams_data);
5553  HYPRE_Int B_Pi_interp_type = ams_data->B_Pi_interp_type;
5554  HYPRE_Int B_Pi_Pmax = ams_data->B_Pi_Pmax;
5555 
5556  HYPRE_Int B_G_coarsen_type = hypre_AMSDataPoissonBetaAMGCoarsenType(ams_data);
5557  HYPRE_Int B_G_agg_levels = hypre_AMSDataPoissonBetaAMGAggLevels(ams_data);
5558  HYPRE_Int B_G_relax_type = hypre_AMSDataPoissonBetaAMGRelaxType(ams_data);
5559  HYPRE_Int B_G_coarse_relax_type = ams_data->B_G_coarse_relax_type;
5560  HYPRE_Real B_G_theta = hypre_AMSDataPoissonBetaAMGStrengthThreshold(ams_data);
5561  HYPRE_Int B_G_interp_type = ams_data->B_G_interp_type;
5562  HYPRE_Int B_G_Pmax = ams_data->B_G_Pmax;
5563 
5564  HYPRE_AMSDestroy(ams);
5565  HYPRE_AMSCreate(&ams);
5566  ams_data = (hypre_AMSData *)ams;
5567 
5568  HYPRE_AMSSetDimension(ams, dim); // 2D H(div) and 3D H(curl) problems
5569  HYPRE_AMSSetTol(ams, tol);
5570  HYPRE_AMSSetMaxIter(ams, maxit); // use as a preconditioner
5571  HYPRE_AMSSetCycleType(ams, cycle_type);
5572  HYPRE_AMSSetPrintLevel(ams, ams_print_level);
5573 
5574  HYPRE_AMSSetCoordinateVectors(ams, hy_x, hy_y, hy_z);
5575 
5576  HYPRE_AMSSetDiscreteGradient(ams, hy_G);
5577  HYPRE_AMSSetCoordinateVectors(ams, hy_x, hy_y, hy_z);
5578  HYPRE_AMSSetInterpolations(ams, hy_Pi, hy_Pix, hy_Piy, hy_Piz);
5579  ams_data->owns_Pi = owns_Pi;
5580 
5581  // set additional AMS options
5582  HYPRE_AMSSetSmoothingOptions(ams, A_relax_type, A_relax_times, A_relax_weight,
5583  A_omega);
5584 
5585  hypre_AMSDataAChebyOrder(ams_data) = A_cheby_order;
5586  hypre_AMSDataAChebyFraction(ams_data) = A_cheby_fraction;
5587 
5588  HYPRE_AMSSetAlphaAMGOptions(ams, B_Pi_coarsen_type, B_Pi_agg_levels,
5589  B_Pi_relax_type,
5590  B_Pi_theta, B_Pi_interp_type, B_Pi_Pmax);
5591  HYPRE_AMSSetBetaAMGOptions(ams, B_G_coarsen_type, B_G_agg_levels,
5592  B_G_relax_type,
5593  B_G_theta, B_G_interp_type, B_G_Pmax);
5594 
5595  HYPRE_AMSSetAlphaAMGCoarseRelaxType(ams, B_Pi_coarse_relax_type);
5596  HYPRE_AMSSetBetaAMGCoarseRelaxType(ams, B_G_coarse_relax_type);
5597 
5598  ams_data->beta_is_zero = beta_is_zero;
5599 
5600 #else
5601  HYPRE_AMSDestroy(ams);
5602 
5603  MakeSolver(space_dim, ams_cycle_type);
5604 
5605  HYPRE_AMSSetPrintLevel(ams, print_level);
5606  if (singular) { HYPRE_AMSSetBetaPoissonMatrix(ams, NULL); }
5607 
5608  HYPRE_AMSSetDiscreteGradient(ams, *G);
5609  if (x != nullptr)
5610  {
5611  HYPRE_AMSSetCoordinateVectors(ams,
5612  x ? (HYPRE_ParVector)(*x) : nullptr,
5613  y ? (HYPRE_ParVector)(*y) : nullptr,
5614  z ? (HYPRE_ParVector)(*z) : nullptr);
5615  }
5616  else
5617  {
5618  HYPRE_AMSSetInterpolations(ams,
5619  Pi ? (HYPRE_ParCSRMatrix) *Pi : nullptr,
5620  Pix ? (HYPRE_ParCSRMatrix) *Pix : nullptr,
5621  Piy ? (HYPRE_ParCSRMatrix) *Piy : nullptr,
5622  Piz ? (HYPRE_ParCSRMatrix) *Piz : nullptr);
5623  }
5624 #endif
5625 }
5626 
5628 {
5629  const HypreParMatrix *new_A = dynamic_cast<const HypreParMatrix *>(&op);
5630  MFEM_VERIFY(new_A, "new Operator must be a HypreParMatrix!");
5631 
5632  if (A) { ResetAMSPrecond(); }
5633 
5634  // update base classes: Operator, Solver, HypreSolver
5635  height = new_A->Height();
5636  width = new_A->Width();
5637  A = const_cast<HypreParMatrix *>(new_A);
5638 
5639  setup_called = 0;
5640  delete X;
5641  delete B;
5642  B = X = NULL;
5643  auxB.Delete(); auxB.Reset();
5644  auxX.Delete(); auxX.Reset();
5645 }
5646 
5648 {
5649  HYPRE_AMSDestroy(ams);
5650 
5651  delete x;
5652  delete y;
5653  delete z;
5654 
5655  delete G;
5656  delete Pi;
5657  delete Pix;
5658  delete Piy;
5659  delete Piz;
5660 }
5661 
5662 void HypreAMS::SetPrintLevel(int print_lvl)
5663 {
5664  HYPRE_AMSSetPrintLevel(ams, print_lvl);
5665  print_level = print_lvl;
5666 }
5667 
5669 {
5670  Init(face_fespace);
5671 }
5672 
5674  : HypreSolver(&A)
5675 {
5676  Init(face_fespace);
5677 }
5678 
5680  const HypreParMatrix &A, HypreParMatrix *C_, HypreParMatrix *G_,
5682  : HypreSolver(&A),
5683  x(x_), y(y_), z(z_),
5684  G(G_), C(C_),
5685  ND_Pi(NULL), ND_Pix(NULL), ND_Piy(NULL), ND_Piz(NULL),
5686  RT_Pi(NULL), RT_Pix(NULL), RT_Piy(NULL), RT_Piz(NULL)
5687 {
5688  MFEM_ASSERT(C != NULL, "");
5689  MFEM_ASSERT(G != NULL, "");
5690  MFEM_ASSERT(x != NULL, "");
5691  MFEM_ASSERT(y != NULL, "");
5692  MFEM_ASSERT(z != NULL, "");
5693 
5694  MakeSolver();
5695 
5696  HYPRE_ADSSetCoordinateVectors(ads, *x, *y, *z);
5697  HYPRE_ADSSetDiscreteCurl(ads, *C);
5698  HYPRE_ADSSetDiscreteGradient(ads, *G);
5699 }
5700 
5701 void HypreADS::MakeSolver()
5702 {
5703  int rlx_sweeps = 1;
5704  double rlx_weight = 1.0;
5705  double rlx_omega = 1.0;
5706 #if !defined(HYPRE_USING_GPU)
5707  int rlx_type = 2;
5708  int amg_coarsen_type = 10;
5709  int amg_agg_levels = 1;
5710  int amg_rlx_type = 8;
5711  double theta = 0.25;
5712  int amg_interp_type = 6;
5713  int amg_Pmax = 4;
5714 #else
5715  int rlx_type = 1;
5716  int amg_coarsen_type = 8;
5717  int amg_agg_levels = 0;
5718  int amg_rlx_type = 18;
5719  double theta = 0.25;
5720  int amg_interp_type = 6;
5721  int amg_Pmax = 4;
5722 #endif
5723 
5724  HYPRE_ADSCreate(&ads);
5725 
5726  HYPRE_ADSSetTol(ads, 0.0);
5727  HYPRE_ADSSetMaxIter(ads, 1); // use as a preconditioner
5728  HYPRE_ADSSetCycleType(ads, cycle_type);
5729  HYPRE_ADSSetPrintLevel(ads, 1);
5730 
5731  // set additional ADS options
5732  HYPRE_ADSSetSmoothingOptions(ads, rlx_type, rlx_sweeps, rlx_weight, rlx_omega);
5733  HYPRE_ADSSetAMGOptions(ads, amg_coarsen_type, amg_agg_levels, amg_rlx_type,
5734  theta, amg_interp_type, amg_Pmax);
5735  HYPRE_ADSSetAMSOptions(ads, ams_cycle_type, amg_coarsen_type, amg_agg_levels,
5736  amg_rlx_type, theta, amg_interp_type, amg_Pmax);
5737 
5738  // The ADS preconditioner requires inverting singular matrices with BoomerAMG,
5739  // which are handled correctly in hypre's Solve method, but can produce hypre
5740  // errors in the Setup (specifically in the l1 row norm computation). See the
5741  // documentation of SetErrorMode() for more details.
5743 }
5744 
5745 void HypreADS::MakeDiscreteMatrices(ParFiniteElementSpace *face_fespace)
5746 {
5747  const FiniteElementCollection *face_fec = face_fespace->FEColl();
5748  bool trace_space =
5749  (dynamic_cast<const RT_Trace_FECollection*>(face_fec) != NULL);
5750  int p = 1;
5751  if (face_fespace->GetNE() > 0)
5752  {
5753  MFEM_VERIFY(!face_fespace->IsVariableOrder(), "");
5754  if (trace_space)
5755  {
5756  p = face_fespace->GetFaceOrder(0) + 1;
5757  }
5758  else
5759  {
5760  p = face_fespace->GetElementOrder(0);
5761  }
5762  }
5763 
5764  // define the nodal and edge finite element spaces associated with face_fespace
5765  ParMesh *pmesh = (ParMesh *) face_fespace->GetMesh();
5766  FiniteElementCollection *vert_fec, *edge_fec;
5767  if (trace_space)
5768  {
5769  vert_fec = new H1_Trace_FECollection(p, 3);
5770  edge_fec = new ND_Trace_FECollection(p, 3);
5771  }
5772  else
5773  {
5774  vert_fec = new H1_FECollection(p, 3);
5775  edge_fec = new ND_FECollection(p, 3);
5776  }
5777 
5778  ParFiniteElementSpace *vert_fespace = new ParFiniteElementSpace(pmesh,
5779  vert_fec);
5780  ParFiniteElementSpace *edge_fespace = new ParFiniteElementSpace(pmesh,
5781  edge_fec);
5782 
5783  // generate and set the vertex coordinates
5784  if (p == 1 && pmesh->GetNodes() == NULL)
5785  {
5786  ParGridFunction x_coord(vert_fespace);
5787  ParGridFunction y_coord(vert_fespace);
5788  ParGridFunction z_coord(vert_fespace);
5789  double *coord;
5790  for (int i = 0; i < pmesh->GetNV(); i++)
5791  {
5792  coord = pmesh -> GetVertex(i);
5793  x_coord(i) = coord[0];
5794  y_coord(i) = coord[1];
5795  z_coord(i) = coord[2];
5796  }
5797  x = x_coord.ParallelProject();
5798  y = y_coord.ParallelProject();
5799  z = z_coord.ParallelProject();
5800  x->HypreReadWrite();
5801  y->HypreReadWrite();
5802  z->HypreReadWrite();
5803  HYPRE_ADSSetCoordinateVectors(ads, *x, *y, *z);
5804  }
5805  else
5806  {
5807  x = NULL;
5808  y = NULL;
5809  z = NULL;
5810  }
5811 
5812  // generate and set the discrete curl
5813  ParDiscreteLinearOperator *curl;
5814  curl = new ParDiscreteLinearOperator(edge_fespace, face_fespace);
5815  if (trace_space)
5816  {
5817  curl->AddTraceFaceInterpolator(new CurlInterpolator);
5818  }
5819  else
5820  {
5821  curl->AddDomainInterpolator(new CurlInterpolator);
5822  }
5823  curl->Assemble();
5824  curl->Finalize();
5825  C = curl->ParallelAssemble();
5826  C->CopyColStarts(); // since we'll delete edge_fespace
5827  HYPRE_ADSSetDiscreteCurl(ads, *C);
5828  delete curl;
5829 
5830  // generate and set the discrete gradient
5831  ParDiscreteLinearOperator *grad;
5832  grad = new ParDiscreteLinearOperator(vert_fespace, edge_fespace);
5833  if (trace_space)
5834  {
5835  grad->AddTraceFaceInterpolator(new GradientInterpolator);
5836  }
5837  else
5838  {
5839  grad->AddDomainInterpolator(new GradientInterpolator);
5840  }
5841  grad->Assemble();
5842  grad->Finalize();
5843  G = grad->ParallelAssemble();
5844  G->CopyColStarts(); // since we'll delete vert_fespace
5845  G->CopyRowStarts(); // since we'll delete edge_fespace
5846  HYPRE_ADSSetDiscreteGradient(ads, *G);
5847  delete grad;
5848 
5849  // generate and set the Nedelec and Raviart-Thomas interpolation matrices
5850  RT_Pi = RT_Pix = RT_Piy = RT_Piz = NULL;
5851  ND_Pi = ND_Pix = ND_Piy = ND_Piz = NULL;
5852  if (p > 1 || pmesh->GetNodes() != NULL)
5853  {
5854  ParFiniteElementSpace *vert_fespace_d
5855  = new ParFiniteElementSpace(pmesh, vert_fec, 3, Ordering::byVDIM);
5856 
5857  ParDiscreteLinearOperator *id_ND;
5858  id_ND = new ParDiscreteLinearOperator(vert_fespace_d, edge_fespace);
5859  if (trace_space)
5860  {
5861  id_ND->AddTraceFaceInterpolator(new IdentityInterpolator);
5862  }
5863  else
5864  {
5865  id_ND->AddDomainInterpolator(new IdentityInterpolator);
5866  }
5867  id_ND->Assemble();
5868  id_ND->Finalize();
5869 
5870  if (ams_cycle_type < 10)
5871  {
5872  ND_Pi = id_ND->ParallelAssemble();
5873  ND_Pi->CopyColStarts(); // since we'll delete vert_fespace_d
5874  ND_Pi->CopyRowStarts(); // since we'll delete edge_fespace
5875  }
5876  else
5877  {
5878  Array2D<HypreParMatrix *> ND_Pi_blocks;
5879  id_ND->GetParBlocks(ND_Pi_blocks);
5880  ND_Pix = ND_Pi_blocks(0,0);
5881  ND_Piy = ND_Pi_blocks(0,1);
5882  ND_Piz = ND_Pi_blocks(0,2);
5883  }
5884 
5885  delete id_ND;
5886 
5887  ParDiscreteLinearOperator *id_RT;
5888  id_RT = new ParDiscreteLinearOperator(vert_fespace_d, face_fespace);
5889  if (trace_space)
5890  {
5891  id_RT->AddTraceFaceInterpolator(new NormalInterpolator);
5892  }
5893  else
5894  {
5895  id_RT->AddDomainInterpolator(new IdentityInterpolator);
5896  }
5897  id_RT->Assemble();
5898  id_RT->Finalize();
5899 
5900  if (cycle_type < 10)
5901  {
5902  RT_Pi = id_RT->ParallelAssemble();
5903  RT_Pi->CopyColStarts(); // since we'll delete vert_fespace_d
5904  }
5905  else
5906  {
5907  Array2D<HypreParMatrix *> RT_Pi_blocks;
5908  id_RT->GetParBlocks(RT_Pi_blocks);
5909  RT_Pix = RT_Pi_blocks(0,0);
5910  RT_Piy = RT_Pi_blocks(0,1);
5911  RT_Piz = RT_Pi_blocks(0,2);
5912  }
5913 
5914  delete id_RT;
5915 
5916  HYPRE_ParCSRMatrix HY_RT_Pi, HY_RT_Pix, HY_RT_Piy, HY_RT_Piz;
5917  HY_RT_Pi = (RT_Pi) ? (HYPRE_ParCSRMatrix) *RT_Pi : NULL;
5918  HY_RT_Pix = (RT_Pix) ? (HYPRE_ParCSRMatrix) *RT_Pix : NULL;
5919  HY_RT_Piy = (RT_Piy) ? (HYPRE_ParCSRMatrix) *RT_Piy : NULL;
5920  HY_RT_Piz = (RT_Piz) ? (HYPRE_ParCSRMatrix) *RT_Piz : NULL;
5921  HYPRE_ParCSRMatrix HY_ND_Pi, HY_ND_Pix, HY_ND_Piy, HY_ND_Piz;
5922  HY_ND_Pi = (ND_Pi) ? (HYPRE_ParCSRMatrix) *ND_Pi : NULL;
5923  HY_ND_Pix = (ND_Pix) ? (HYPRE_ParCSRMatrix) *ND_Pix : NULL;
5924  HY_ND_Piy = (ND_Piy) ? (HYPRE_ParCSRMatrix) *ND_Piy : NULL;
5925  HY_ND_Piz = (ND_Piz) ? (HYPRE_ParCSRMatrix) *ND_Piz : NULL;
5926  HYPRE_ADSSetInterpolations(ads,
5927  HY_RT_Pi, HY_RT_Pix, HY_RT_Piy, HY_RT_Piz,
5928  HY_ND_Pi, HY_ND_Pix, HY_ND_Piy, HY_ND_Piz);
5929 
5930  delete vert_fespace_d;
5931  }
5932 
5933  delete vert_fec;
5934  delete vert_fespace;
5935  delete edge_fec;
5936  delete edge_fespace;
5937 }
5938 
5939 void HypreADS::Init(ParFiniteElementSpace *face_fespace)
5940 {
5941  MakeSolver();
5942  MakeDiscreteMatrices(face_fespace);
5943 }
5944 
5945 void HypreADS::ResetADSPrecond()
5946 {
5947  HYPRE_ADSDestroy(ads);
5948 
5949  MakeSolver();
5950 
5951  HYPRE_ADSSetPrintLevel(ads, print_level);
5952 
5953  HYPRE_ADSSetDiscreteCurl(ads, *C);
5954  HYPRE_ADSSetDiscreteGradient(ads, *G);
5955  if (x != nullptr)
5956  {
5957  MFEM_VERIFY(x && y && z, "");
5958  HYPRE_ADSSetCoordinateVectors(ads, *x, *y, *z);
5959  }
5960  else
5961  {
5962  HYPRE_ParCSRMatrix HY_RT_Pi, HY_RT_Pix, HY_RT_Piy, HY_RT_Piz;
5963  HY_RT_Pi = (RT_Pi) ? (HYPRE_ParCSRMatrix) *RT_Pi : NULL;
5964  HY_RT_Pix = (RT_Pix) ? (HYPRE_ParCSRMatrix) *RT_Pix : NULL;
5965  HY_RT_Piy = (RT_Piy) ? (HYPRE_ParCSRMatrix) *RT_Piy : NULL;
5966  HY_RT_Piz = (RT_Piz) ? (HYPRE_ParCSRMatrix) *RT_Piz : NULL;
5967  HYPRE_ParCSRMatrix HY_ND_Pi, HY_ND_Pix, HY_ND_Piy, HY_ND_Piz;
5968  HY_ND_Pi = (ND_Pi) ? (HYPRE_ParCSRMatrix) *ND_Pi : NULL;
5969  HY_ND_Pix = (ND_Pix) ? (HYPRE_ParCSRMatrix) *ND_Pix : NULL;
5970  HY_ND_Piy = (ND_Piy) ? (HYPRE_ParCSRMatrix) *ND_Piy : NULL;
5971  HY_ND_Piz = (ND_Piz) ? (HYPRE_ParCSRMatrix) *ND_Piz : NULL;
5972  HYPRE_ADSSetInterpolations(ads,
5973  HY_RT_Pi, HY_RT_Pix, HY_RT_Piy, HY_RT_Piz,
5974  HY_ND_Pi, HY_ND_Pix, HY_ND_Piy, HY_ND_Piz);
5975  }
5976 }
5977 
5979 {
5980  const HypreParMatrix *new_A = dynamic_cast<const HypreParMatrix *>(&op);
5981  MFEM_VERIFY(new_A, "new Operator must be a HypreParMatrix!");
5982 
5983  if (A) { ResetADSPrecond(); }
5984 
5985  // update base classes: Operator, Solver, HypreSolver
5986  height = new_A->Height();
5987  width = new_A->Width();
5988  A = const_cast<HypreParMatrix *>(new_A);
5989 
5990  setup_called = 0;
5991  delete X;
5992  delete B;
5993  B = X = NULL;
5994  auxB.Delete(); auxB.Reset();
5995  auxX.Delete(); auxX.Reset();
5996 }
5997 
5999 {
6000  HYPRE_ADSDestroy(ads);
6001 
6002  delete x;
6003  delete y;
6004  delete z;
6005 
6006  delete G;
6007  delete C;
6008 
6009  delete RT_Pi;
6010  delete RT_Pix;
6011  delete RT_Piy;
6012  delete RT_Piz;
6013 
6014  delete ND_Pi;
6015  delete ND_Pix;
6016  delete ND_Piy;
6017  delete ND_Piz;
6018 }
6019 
6020 void HypreADS::SetPrintLevel(int print_lvl)
6021 {
6022  HYPRE_ADSSetPrintLevel(ads, print_lvl);
6023  print_level = print_lvl;
6024 }
6025 
6026 HypreLOBPCG::HypreMultiVector::HypreMultiVector(int n, HypreParVector & v,
6027  mv_InterfaceInterpreter & interpreter)
6028  : hpv(NULL),
6029  nv(n)
6030 {
6031  mv_ptr = mv_MultiVectorCreateFromSampleVector(&interpreter, nv,
6032  (HYPRE_ParVector)v);
6033 
6034  HYPRE_ParVector* vecs = NULL;
6035  {
6036  mv_TempMultiVector* tmp =
6037  (mv_TempMultiVector*)mv_MultiVectorGetData(mv_ptr);
6038  vecs = (HYPRE_ParVector*)(tmp -> vector);
6039  }
6040 
6041  hpv = new HypreParVector*[nv];
6042  for (int i=0; i<nv; i++)
6043  {
6044  hpv[i] = new HypreParVector(vecs[i]);
6045  }
6046 }
6047 
6048 HypreLOBPCG::HypreMultiVector::~HypreMultiVector()
6049 {
6050  if ( hpv != NULL )
6051  {
6052  for (int i=0; i<nv; i++)
6053  {
6054  delete hpv[i];
6055  }
6056  delete [] hpv;
6057  }
6058 
6059  mv_MultiVectorDestroy(mv_ptr);
6060 }
6061 
6062 void
6063 HypreLOBPCG::HypreMultiVector::Randomize(HYPRE_Int seed_)
6064 {
6065  mv_MultiVectorSetRandom(mv_ptr, seed_);
6066 }
6067 
6068 HypreParVector &
6069 HypreLOBPCG::HypreMultiVector::GetVector(unsigned int i)
6070 {
6071  MFEM_ASSERT((int)i < nv, "index out of range");
6072 
6073  return ( *hpv[i] );
6074 }
6075 
6076 HypreParVector **
6077 HypreLOBPCG::HypreMultiVector::StealVectors()
6078 {
6079  HypreParVector ** hpv_ret = hpv;
6080 
6081  hpv = NULL;
6082 
6083  mv_TempMultiVector * mv_tmp =
6084  (mv_TempMultiVector*)mv_MultiVectorGetData(mv_ptr);
6085 
6086  mv_tmp->ownsVectors = 0;
6087 
6088  for (int i=0; i<nv; i++)
6089  {
6090  hpv_ret[i]->SetOwnership(1);
6091  }
6092 
6093  return hpv_ret;
6094 }
6095 
6097  : comm(c),
6098  myid(0),
6099  numProcs(1),
6100  nev(10),
6101  seed(75),
6102  glbSize(-1),
6103  part(NULL),
6104  multi_vec(NULL),
6105  x(NULL),
6106  subSpaceProj(NULL)
6107 {
6108  MPI_Comm_size(comm,&numProcs);
6109  MPI_Comm_rank(comm,&myid);
6110 
6111  HYPRE_ParCSRSetupInterpreter(&interpreter);
6112  HYPRE_ParCSRSetupMatvec(&matvec_fn);
6113  HYPRE_LOBPCGCreate(&interpreter, &matvec_fn, &lobpcg_solver);
6114 }
6115 
6117 {
6118  delete multi_vec;
6119  delete x;
6120  delete [] part;
6121 
6122  HYPRE_LOBPCGDestroy(lobpcg_solver);
6123 }
6124 
6125 void
6127 {
6128  HYPRE_LOBPCGSetTol(lobpcg_solver, tol);
6129 }
6130 
6131 void
6132 HypreLOBPCG::SetRelTol(double rel_tol)
6133 {
6134 #if MFEM_HYPRE_VERSION >= 21101
6135  HYPRE_LOBPCGSetRTol(lobpcg_solver, rel_tol);
6136 #else
6137  MFEM_ABORT("This method requires HYPRE version >= 2.11.1");
6138 #endif
6139 }
6140 
6141 void
6143 {
6144  HYPRE_LOBPCGSetMaxIter(lobpcg_solver, max_iter);
6145 }
6146 
6147 void
6149 {
6150  if (myid == 0)
6151  {
6152  HYPRE_LOBPCGSetPrintLevel(lobpcg_solver, logging);
6153  }
6154 }
6155 
6156 void
6158 {
6159  HYPRE_LOBPCGSetPrecondUsageMode(lobpcg_solver, pcg_mode);
6160 }
6161 
6162 void
6164 {
6165  HYPRE_LOBPCGSetPrecond(lobpcg_solver,
6166  (HYPRE_PtrToSolverFcn)this->PrecondSolve,
6167  (HYPRE_PtrToSolverFcn)this->PrecondSetup,
6168  (HYPRE_Solver)&precond);
6169 }
6170 
6171 void
6173 {
6174  HYPRE_BigInt locSize = A.Width();
6175 
6176  if (HYPRE_AssumedPartitionCheck())
6177  {
6178  part = new HYPRE_BigInt[2];
6179 
6180  MPI_Scan(&locSize, &part[1], 1, HYPRE_MPI_BIG_INT, MPI_SUM, comm);
6181 
6182  part[0] = part[1] - locSize;
6183 
6184  MPI_Allreduce(&locSize, &glbSize, 1, HYPRE_MPI_BIG_INT, MPI_SUM, comm);
6185  }
6186  else
6187  {
6188  part = new HYPRE_BigInt[numProcs+1];
6189 
6190  MPI_Allgather(&locSize, 1, HYPRE_MPI_BIG_INT,
6191  &part[1], 1, HYPRE_MPI_BIG_INT, comm);
6192 
6193  part[0] = 0;
6194  for (int i=0; i<numProcs; i++)
6195  {
6196  part[i+1] += part[i];
6197  }
6198 
6199  glbSize = part[numProcs];
6200  }
6201 
6202  if ( x != NULL )
6203  {
6204  delete x;
6205  }
6206 
6207  // Create a distributed vector without a data array.
6208  const bool is_device_ptr = true;
6209  x = new HypreParVector(comm,glbSize,NULL,part,is_device_ptr);
6210 
6211  matvec_fn.MatvecCreate = this->OperatorMatvecCreate;
6212  matvec_fn.Matvec = this->OperatorMatvec;
6213  matvec_fn.MatvecDestroy = this->OperatorMatvecDestroy;
6214 
6215  HYPRE_LOBPCGSetup(lobpcg_solver,(HYPRE_Matrix)&A,NULL,NULL);
6216 }
6217 
6218 void
6220 {
6221  matvec_fn.MatvecCreate = this->OperatorMatvecCreate;
6222  matvec_fn.Matvec = this->OperatorMatvec;
6223  matvec_fn.MatvecDestroy = this->OperatorMatvecDestroy;
6224 
6225  HYPRE_LOBPCGSetupB(lobpcg_solver,(HYPRE_Matrix)&M,NULL);
6226 }
6227 
6228 void
6230 {
6231  // Initialize eigenvalues array with marker values
6232  eigs.SetSize(nev);
6233 
6234  for (int i=0; i<nev; i++)
6235  {
6236  eigs[i] = eigenvalues[i];
6237  }
6238 }
6239 
6240 const HypreParVector &
6241 HypreLOBPCG::GetEigenvector(unsigned int i) const
6242 {
6243  return multi_vec->GetVector(i);
6244 }
6245 
6246 void
6248 {
6249  // Initialize HypreMultiVector object if necessary
6250  if ( multi_vec == NULL )
6251  {
6252  MFEM_ASSERT(x != NULL, "In HypreLOBPCG::SetInitialVectors()");
6253 
6254  multi_vec = new HypreMultiVector(nev, *x, interpreter);
6255  }
6256 
6257  // Copy the vectors provided
6258  for (int i=0; i < min(num_vecs,nev); i++)
6259  {
6260  multi_vec->GetVector(i) = *vecs[i];
6261  }
6262 
6263  // Randomize any remaining vectors
6264  for (int i=min(num_vecs,nev); i < nev; i++)
6265  {
6266  multi_vec->GetVector(i).Randomize(seed);
6267  }
6268 
6269  // Ensure all vectors are in the proper subspace
6270  if ( subSpaceProj != NULL )
6271  {
6272  HypreParVector y(*x);
6273  y = multi_vec->GetVector(0);
6274 
6275  for (int i=1; i<nev; i++)
6276  {
6277  subSpaceProj->Mult(multi_vec->GetVector(i),
6278  multi_vec->GetVector(i-1));
6279  }
6280  subSpaceProj->Mult(y,
6281  multi_vec->GetVector(nev-1));
6282  }
6283 }
6284 
6285 void
6287 {
6288  // Initialize HypreMultiVector object if necessary
6289  if ( multi_vec == NULL )
6290  {
6291  MFEM_ASSERT(x != NULL, "In HypreLOBPCG::Solve()");
6292 
6293  multi_vec = new HypreMultiVector(nev, *x, interpreter);
6294  multi_vec->Randomize(seed);
6295 
6296  if ( subSpaceProj != NULL )
6297  {
6298  HypreParVector y(*x);
6299  y = multi_vec->GetVector(0);
6300 
6301  for (int i=1; i<nev; i++)
6302  {
6303  subSpaceProj->Mult(multi_vec->GetVector(i),
6304  multi_vec->GetVector(i-1));
6305  }
6306  subSpaceProj->Mult(y, multi_vec->GetVector(nev-1));
6307  }
6308  }
6309 
6310  eigenvalues.SetSize(nev);
6311  eigenvalues = NAN;
6312 
6313  // Perform eigenmode calculation
6314  //
6315  // The eigenvalues are computed in ascending order (internally the
6316  // order is determined by the LAPACK routine 'dsydv'.)
6317  HYPRE_LOBPCGSolve(lobpcg_solver, NULL, *multi_vec, eigenvalues);
6318 }
6319 
6320 void *
6321 HypreLOBPCG::OperatorMatvecCreate( void *A,
6322  void *x )
6323 {
6324  void *matvec_data;
6325 
6326  matvec_data = NULL;
6327 
6328  return ( matvec_data );
6329 }
6330 
6331 HYPRE_Int
6332 HypreLOBPCG::OperatorMatvec( void *matvec_data,
6333  HYPRE_Complex alpha,
6334  void *A,
6335  void *x,
6336  HYPRE_Complex beta,
6337  void *y )
6338 {
6339  MFEM_VERIFY(alpha == 1.0 && beta == 0.0, "values not supported");
6340 
6341  Operator *Aop = (Operator*)A;
6342 
6343  hypre_ParVector * xPar = (hypre_ParVector *)x;
6344  hypre_ParVector * yPar = (hypre_ParVector *)y;
6345 
6346  HypreParVector xVec(xPar);
6347  HypreParVector yVec(yPar);
6348 
6349  Aop->Mult( xVec, yVec );
6350 
6351  // Move data back to hypre's device memory location in case the above Mult
6352  // operation moved it to host.
6353  yVec.HypreReadWrite();
6354 
6355  return 0;
6356 }
6357 
6358 HYPRE_Int
6359 HypreLOBPCG::OperatorMatvecDestroy( void *matvec_data )
6360 {
6361  return 0;
6362 }
6363 
6364 HYPRE_Int
6365 HypreLOBPCG::PrecondSolve(void *solver,
6366  void *A,
6367  void *b,
6368  void *x)
6369 {
6370  Solver *PC = (Solver*)solver;
6371 
6372  hypre_ParVector * bPar = (hypre_ParVector *)b;
6373  hypre_ParVector * xPar = (hypre_ParVector *)x;
6374 
6375  HypreParVector bVec(bPar);
6376  HypreParVector xVec(xPar);
6377 
6378  PC->Mult( bVec, xVec );
6379 
6380  // Move data back to hypre's device memory location in case the above Mult
6381  // operation moved it to host.
6382  xVec.HypreReadWrite();
6383 
6384  return 0;
6385 }
6386 
6387 HYPRE_Int
6388 HypreLOBPCG::PrecondSetup(void *solver,
6389  void *A,
6390  void *b,
6391  void *x)
6392 {
6393  return 0;
6394 }
6395 
6396 HypreAME::HypreAME(MPI_Comm comm)
6397  : myid(0),
6398  numProcs(1),
6399  nev(10),
6400  setT(false),
6401  ams_precond(NULL),
6402  eigenvalues(NULL),
6403  multi_vec(NULL),
6404  eigenvectors(NULL)
6405 {
6406  MPI_Comm_size(comm,&numProcs);
6407  MPI_Comm_rank(comm,&myid);
6408 
6409  HYPRE_AMECreate(&ame_solver);
6410  HYPRE_AMESetPrintLevel(ame_solver, 0);
6411 }
6412 
6414 {
6415  if ( multi_vec )
6416  {
6417  mfem_hypre_TFree_host(multi_vec);
6418  }
6419 
6420  if ( eigenvectors )
6421  {
6422  for (int i=0; i<nev; i++)
6423  {
6424  delete eigenvectors[i];
6425  }
6426  }
6427  delete [] eigenvectors;
6428 
6429  if ( eigenvalues )
6430  {
6431  mfem_hypre_TFree_host(eigenvalues);
6432  }
6433 
6434  HYPRE_AMEDestroy(ame_solver);
6435 }
6436 
6437 void
6439 {
6440  nev = num_eigs;
6441 
6442  HYPRE_AMESetBlockSize(ame_solver, nev);
6443 }
6444 
6445 void
6446 HypreAME::SetTol(double tol)
6447 {
6448  HYPRE_AMESetTol(ame_solver, tol);
6449 }
6450 
6451 void
6452 HypreAME::SetRelTol(double rel_tol)
6453 {
6454 #if MFEM_HYPRE_VERSION >= 21101
6455  HYPRE_AMESetRTol(ame_solver, rel_tol);
6456 #else
6457  MFEM_ABORT("This method requires HYPRE version >= 2.11.1");
6458 #endif
6459 }
6460 
6461 void
6463 {
6464  HYPRE_AMESetMaxIter(ame_solver, max_iter);
6465 }
6466 
6467 void
6469 {
6470  if (myid == 0)
6471  {
6472  HYPRE_AMESetPrintLevel(ame_solver, logging);
6473  }
6474 }
6475 
6476 void
6478 {
6479  ams_precond = &precond;
6480 }
6481 
6482 void
6484 {
6485  if ( !setT )
6486  {
6487  HYPRE_Solver ams_precond_ptr = (HYPRE_Solver)*ams_precond;
6488 
6489  ams_precond->SetupFcn()(*ams_precond,A,NULL,NULL);
6490 
6491  HYPRE_AMESetAMSSolver(ame_solver, ams_precond_ptr);
6492  }
6493 
6494  HYPRE_AMESetup(ame_solver);
6495 }
6496 
6497 void
6499 {
6500  HYPRE_ParCSRMatrix parcsr_M = M;
6501  HYPRE_AMESetMassMatrix(ame_solver,(HYPRE_ParCSRMatrix)parcsr_M);
6502 }
6503 
6504 void
6506 {
6507  HYPRE_AMESolve(ame_solver);
6508 
6509  // Grab a pointer to the eigenvalues from AME
6510  HYPRE_AMEGetEigenvalues(ame_solver,&eigenvalues);
6511 
6512  // Grad a pointer to the eigenvectors from AME
6513  HYPRE_AMEGetEigenvectors(ame_solver,&multi_vec);
6514 }
6515 
6516 void
6518 {
6519  // Initialize eigenvalues array with marker values
6520  eigs.SetSize(nev); eigs = -1.0;
6521 
6522  // Copy eigenvalues to eigs array
6523  for (int i=0; i<nev; i++)
6524  {
6525  eigs[i] = eigenvalues[i];
6526  }
6527 }
6528 
6529 void
6530 HypreAME::createDummyVectors() const
6531 {
6532  eigenvectors = new HypreParVector*[nev];
6533  for (int i=0; i<nev; i++)
6534  {
6535  eigenvectors[i] = new HypreParVector(multi_vec[i]);
6536  eigenvectors[i]->SetOwnership(1);
6537  }
6538 }
6539 
6540 const HypreParVector &
6541 HypreAME::GetEigenvector(unsigned int i) const
6542 {
6543  if ( eigenvectors == NULL )
6544  {
6545  this->createDummyVectors();
6546  }
6547 
6548  return *eigenvectors[i];
6549 }
6550 
6551 HypreParVector **
6553 {
6554  if ( eigenvectors == NULL )
6555  {
6556  this->createDummyVectors();
6557  }
6558 
6559  // Set the local pointers to NULL so that they won't be deleted later
6560  HypreParVector ** vecs = eigenvectors;
6561  eigenvectors = NULL;
6562  multi_vec = NULL;
6563 
6564  return vecs;
6565 }
6566 
6567 }
6568 
6569 #endif
MemoryType GetHypreMemoryType()
The MemoryType used by MFEM when allocating arrays for Hypre objects.
Definition: hypre.hpp:149
virtual ~HypreBoomerAMG()
Definition: hypre.cpp:5239
void SetAbsTol(double atol)
Definition: hypre.cpp:4000
Hash function for data sequences.
Definition: hash.hpp:455
Memory< double > auxX
Definition: hypre.hpp:972
void SetTol(double tol)
Definition: hypre.cpp:3995
void SetPreconditioner(HypreSolver &precond)
Set the hypre solver to be used as a preconditioner.
Definition: hypre.cpp:4210
void SetType(HYPRE_Int ilu_type)
Definition: hypre.cpp:4751
int Size() const
Return the logical size of the array.
Definition: array.hpp:138
void EliminateRowsCols(const Array< int > &rows_cols, const HypreParVector &X, HypreParVector &B)
Definition: hypre.cpp:2272
HypreADS(ParFiniteElementSpace *face_fespace)
Definition: hypre.cpp:5668
virtual int NumNonZeroElems() const
Returns the number of the nonzero elements in the matrix.
Definition: sparsemat.cpp:1584
static constexpr Type default_type
Definition: hypre.hpp:1034
MPI_Comm GetComm() const
MPI communicator.
Definition: hypre.hpp:534
HypreEuclid(MPI_Comm comm)
Definition: hypre.cpp:4606
virtual void SetOperator(const Operator &op)
Set/update the solver for the given operator.
Definition: hypre.cpp:4330
Vector * GlobalVector() const
Returns the global vector in each processor.
Definition: hypre.cpp:264
void EliminateBC(const HypreParMatrix &A, const HypreParMatrix &Ae, const Array< int > &ess_dof_list, const Vector &X, Vector &B)
Eliminate essential BC specified by ess_dof_list from the solution X to the r.h.s. B.
Definition: hypre.cpp:3263
HypreParVector * NewTrueDofVector()
Definition: pfespace.hpp:331
HypreParVector * X0
FIR Filter Temporary Vectors.
Definition: hypre.hpp:976
void Clear()
Clear the contents of the SparseMatrix.
Definition: sparsemat.hpp:209
double min_eig_est
Minimal eigenvalue estimate for polynomial smoothing.
Definition: hypre.hpp:1008
Memory< T > & GetMemory()
Return a reference to the Memory object used by the Array.
Definition: array.hpp:117
int NumCols() const
Definition: array.hpp:373
void MakeDataOwner() const
Set the Vector data (host pointer) ownership flag.
Definition: vector.hpp:188
Memory< double > data
Definition: vector.hpp:64
static MemoryClass GetHostMemoryClass()
Get the current Host MemoryClass. This is the MemoryClass used by most MFEM host Memory objects...
Definition: device.hpp:268
Memory< HYPRE_Int > I