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