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