MFEM  v4.3.0
Finite element discretization library
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Pages
hypre.hpp
Go to the documentation of this file.
1 // Copyright (c) 2010-2021, 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 #ifndef MFEM_HYPRE
13 #define MFEM_HYPRE
14 
15 #include "../config/config.hpp"
16 
17 #ifdef MFEM_USE_MPI
18 
19 #include <mpi.h>
20 
21 // Enable internal hypre timing routines
22 #define HYPRE_TIMING
23 
24 // hypre header files
25 #include "seq_mv.h"
26 #include "_hypre_parcsr_mv.h"
27 #include "_hypre_parcsr_ls.h"
28 #include "temp_multivector.h"
29 #include "../general/globals.hpp"
30 
31 #ifdef HYPRE_COMPLEX
32 #error "MFEM does not work with HYPRE's complex numbers support"
33 #endif
34 
35 #if defined(HYPRE_USING_CUDA) && !defined(MFEM_USE_CUDA)
36 #error "MFEM_USE_CUDA=YES is required when HYPRE is built with CUDA!"
37 #endif
38 
39 #include "sparsemat.hpp"
40 #include "hypre_parcsr.hpp"
41 
42 namespace mfem
43 {
44 
45 class ParFiniteElementSpace;
46 class HypreParMatrix;
47 
48 namespace internal
49 {
50 
51 template <typename int_type>
52 inline int to_int(int_type i)
53 {
54  MFEM_ASSERT(int_type(int(i)) == i, "overflow converting int_type to int");
55  return int(i);
56 }
57 
58 // Specialization for to_int(int)
59 template <> inline int to_int(int i) { return i; }
60 
61 // Convert a HYPRE_Int to int
62 #ifdef HYPRE_BIGINT
63 template <>
64 inline int to_int(HYPRE_Int i)
65 {
66  MFEM_ASSERT(HYPRE_Int(int(i)) == i, "overflow converting HYPRE_Int to int");
67  return int(i);
68 }
69 #endif
70 
71 } // namespace internal
72 
73 
74 /// The MemoryClass used by Hypre objects.
75 inline constexpr MemoryClass GetHypreMemoryClass()
76 {
77 #ifndef HYPRE_USING_CUDA
78  return MemoryClass::HOST;
79 #elif defined(HYPRE_USING_UNIFIED_MEMORY)
80  return MemoryClass::MANAGED;
81 #else
82  return MemoryClass::DEVICE;
83 #endif
84 }
85 
86 /// The MemoryType used by MFEM when allocating arrays for Hypre objects.
88 {
89 #ifndef HYPRE_USING_CUDA
91 #elif defined(HYPRE_USING_UNIFIED_MEMORY)
92  return MemoryType::MANAGED;
93 #else
94  return MemoryType::DEVICE;
95 #endif
96 }
97 
98 /// Wrapper for hypre's parallel vector class
99 class HypreParVector : public Vector
100 {
101 private:
102  int own_ParVector;
103 
104  /// The actual object
105  hypre_ParVector *x;
106 
107  friend class HypreParMatrix;
108 
109  // Set Vector::data and Vector::size from *x
110  inline void _SetDataAndSize_();
111 
112 public:
113 
114  /// Default constructor, no underlying @a hypre_ParVector is created.
116  {
117  own_ParVector = false;
118  x = NULL;
119  }
120 
121  /** @brief Creates vector with given global size and parallel partitioning of
122  the rows/columns given by @a col. */
123  /** @anchor hypre_partitioning_descr
124  The partitioning is defined in one of two ways depending on the
125  configuration of HYPRE:
126  1. If HYPRE_AssumedPartitionCheck() returns true (the default),
127  then col is of length 2 and the local processor owns columns
128  [col[0],col[1]).
129  2. If HYPRE_AssumedPartitionCheck() returns false, then col is of
130  length (number of processors + 1) and processor P owns columns
131  [col[P],col[P+1]) i.e. each processor has a copy of the same col
132  array. */
133  HypreParVector(MPI_Comm comm, HYPRE_BigInt glob_size, HYPRE_BigInt *col);
134  /** @brief Creates vector with given global size, partitioning of the
135  columns, and data. */
136  /** The data must be allocated and destroyed outside. If @a data_ is NULL, a
137  dummy vector without a valid data array will be created. See @ref
138  hypre_partitioning_descr "here" for a description of the @a col array.
139 
140  If @a is_device_ptr is true, the pointer @a data_ is assumed to be
141  allocated in the memory location HYPRE_MEMORY_DEVICE. */
142  HypreParVector(MPI_Comm comm, HYPRE_BigInt glob_size, double *data_,
143  HYPRE_BigInt *col, bool is_device_ptr = false);
144  /// Creates vector compatible with y
145  HypreParVector(const HypreParVector &y);
146  /// Creates vector compatible with (i.e. in the domain of) A or A^T
147  explicit HypreParVector(const HypreParMatrix &A, int transpose = 0);
148  /// Creates vector wrapping y
149  explicit HypreParVector(HYPRE_ParVector y);
150  /// Create a true dof parallel vector on a given ParFiniteElementSpace
151  explicit HypreParVector(ParFiniteElementSpace *pfes);
152 
153  /// MPI communicator
154  MPI_Comm GetComm() const { return x->comm; }
155 
156  /// Converts hypre's format to HypreParVector
157  void WrapHypreParVector(hypre_ParVector *y, bool owner=true);
158 
159  /// Returns the parallel row/column partitioning
160  /** See @ref hypre_partitioning_descr "here" for a description of the
161  partitioning array. */
162  inline const HYPRE_BigInt *Partitioning() const { return x->partitioning; }
163 
164  /// @brief Returns a non-const pointer to the parallel row/column
165  /// partitioning.
166  /// Deprecated in favor of HypreParVector::Partitioning() const.
167  MFEM_DEPRECATED
168  inline HYPRE_BigInt *Partitioning() { return x->partitioning; }
169 
170  /// Returns the global number of rows
171  inline HYPRE_BigInt GlobalSize() const { return x->global_size; }
172 
173  /// Typecasting to hypre's hypre_ParVector*
174  operator hypre_ParVector*() const { return x; }
175 #ifndef HYPRE_PAR_VECTOR_STRUCT
176  /// Typecasting to hypre's HYPRE_ParVector, a.k.a. void *
177  operator HYPRE_ParVector() const { return (HYPRE_ParVector) x; }
178 #endif
179  /// Changes the ownership of the vector
180  hypre_ParVector *StealParVector() { own_ParVector = 0; return x; }
181 
182  /// Sets ownership of the internal hypre_ParVector
183  void SetOwnership(int own) { own_ParVector = own; }
184 
185  /// Gets ownership of the internal hypre_ParVector
186  int GetOwnership() const { return own_ParVector; }
187 
188  /// Returns the global vector in each processor
189  Vector* GlobalVector() const;
190 
191  /// Set constant values
192  HypreParVector& operator= (double d);
193  /// Define '=' for hypre vectors.
195 
196  using Vector::Read;
197 
198  /// Sets the data of the Vector and the hypre_ParVector to @a data_.
199  /** Must be used only for HypreParVector%s that do not own the data,
200  e.g. created with the constructor:
201  HypreParVector(MPI_Comm, HYPRE_BigInt, double *, HYPRE_BigInt *). */
202  void SetData(double *data_);
203 
204  /** @brief Prepare the HypreParVector for read access in hypre's device
205  memory space, HYPRE_MEMORY_DEVICE. */
206  void HypreRead() const;
207 
208  /** @brief Prepare the HypreParVector for read and write access in hypre's
209  device memory space, HYPRE_MEMORY_DEVICE. */
210  void HypreReadWrite();
211 
212  /** @brief Prepare the HypreParVector for write access in hypre's device
213  memory space, HYPRE_MEMORY_DEVICE. */
214  void HypreWrite();
215 
216  /** @brief Replace the HypreParVector's data with the given Memory, @a mem,
217  and prepare the vector for read access in hypre's device memory space,
218  HYPRE_MEMORY_DEVICE. */
219  /** This method must be used with HypreParVector%s that do not own the data,
220  e.g. created with the constructor:
221  HypreParVector(MPI_Comm, HYPRE_BigInt, double *, HYPRE_BigInt *).
222 
223  The Memory @a mem must be accessible with the hypre MemoryClass defined
224  by GetHypreMemoryClass(). */
225  void WrapMemoryRead(const Memory<double> &mem);
226 
227  /** @brief Replace the HypreParVector's data with the given Memory, @a mem,
228  and prepare the vector for read and write access in hypre's device memory
229  space, HYPRE_MEMORY_DEVICE. */
230  /** This method must be used with HypreParVector%s that do not own the data,
231  e.g. created with the constructor:
232  HypreParVector(MPI_Comm, HYPRE_BigInt, double *, HYPRE_BigInt *).
233 
234  The Memory @a mem must be accessible with the hypre MemoryClass defined
235  by GetHypreMemoryClass(). */
237 
238  /** @brief Replace the HypreParVector's data with the given Memory, @a mem,
239  and prepare the vector for write access in hypre's device memory space,
240  HYPRE_MEMORY_DEVICE. */
241  /** This method must be used with HypreParVector%s that do not own the data,
242  e.g. created with the constructor:
243  HypreParVector(MPI_Comm, HYPRE_BigInt, double *, HYPRE_BigInt *).
244 
245  The Memory @a mem must be accessible with the hypre MemoryClass defined
246  by GetHypreMemoryClass(). */
247  void WrapMemoryWrite(Memory<double> &mem);
248 
249  /// Set random values
250  HYPRE_Int Randomize(HYPRE_Int seed);
251 
252  /// Prints the locally owned rows in parallel
253  void Print(const char *fname) const;
254 
255  /// Calls hypre's destroy function
256  ~HypreParVector();
257 
258 #ifdef MFEM_USE_SUNDIALS
259  /// (DEPRECATED) Return a new wrapper SUNDIALS N_Vector of type SUNDIALS_NVEC_PARALLEL.
260  /** @deprecated The returned N_Vector must be destroyed by the caller. */
261  MFEM_DEPRECATED virtual N_Vector ToNVector();
262  using Vector::ToNVector;
263 #endif
264 };
265 
266 /// Returns the inner product of x and y
267 double InnerProduct(HypreParVector &x, HypreParVector &y);
268 double InnerProduct(HypreParVector *x, HypreParVector *y);
269 
270 
271 /** @brief Compute the l_p norm of the Vector which is split without overlap
272  across the given communicator. */
273 double ParNormlp(const Vector &vec, double p, MPI_Comm comm);
274 
275 
276 /// Wrapper for hypre's ParCSR matrix class
277 class HypreParMatrix : public Operator
278 {
279 private:
280  /// The actual object
281  hypre_ParCSRMatrix *A;
282 
283  /// Auxiliary vectors for typecasting
284  mutable HypreParVector *X, *Y;
285  /** @brief Auxiliary buffers for the case when the input or output arrays in
286  methods like Mult(double, const Vector &, double, Vector &) need to be
287  deep copied in order to be used by hypre. */
288  mutable Memory<double> auxX, auxY;
289 
290  // Flags indicating ownership of A->diag->{i,j,data}, A->offd->{i,j,data},
291  // and A->col_map_offd.
292  // The possible values for diagOwner are:
293  // -1: no special treatment of A->diag (default)
294  // when hypre is built with CUDA support, A->diag owns the "host"
295  // pointers (according to A->diag->owns_data)
296  // -2: used when hypre is built with CUDA support, A->diag owns the "hypre"
297  // pointers (according to A->diag->owns_data)
298  // 0: prevent hypre from destroying A->diag->{i,j,data}
299  // 1: same as 0, plus own the "host" A->diag->{i,j}
300  // 2: same as 0, plus own the "host" A->diag->data
301  // 3: same as 0, plus own the "host" A->diag->{i,j,data}
302  // The same values and rules apply to offdOwner and A->offd.
303  // The possible values for colMapOwner are:
304  // -1: no special treatment of A->col_map_offd (default)
305  // 0: prevent hypre from destroying A->col_map_offd
306  // 1: same as 0, plus take ownership of A->col_map_offd
307  // All owned arrays are destroyed with 'delete []'.
308  signed char diagOwner, offdOwner, colMapOwner;
309 
310  // Does the object own the pointer A?
311  signed char ParCSROwner;
312 
313  MemoryIJData mem_diag, mem_offd;
314 
315  // Initialize with defaults. Does not initialize inherited members.
316  void Init();
317 
318  // Delete all owned data. Does not perform re-initialization with defaults.
319  void Destroy();
320 
321  void Read(MemoryClass mc) const;
322  void ReadWrite(MemoryClass mc);
323  // The Boolean flags are used in Destroy().
324  void Write(MemoryClass mc, bool set_diag = true, bool set_offd = true);
325 
326  // Copy (shallow/deep, based on HYPRE_BIGINT) the I and J arrays from csr to
327  // hypre_csr. Shallow copy the data. Return the appropriate ownership flag.
328  // The CSR arrays are wrapped in the mem_csr struct which is used to move
329  // these arrays to device, if necessary.
330  static signed char CopyCSR(SparseMatrix *csr,
331  MemoryIJData &mem_csr,
332  hypre_CSRMatrix *hypre_csr,
333  bool mem_owner);
334  // Copy (shallow or deep, based on HYPRE_BIGINT) the I and J arrays from
335  // bool_csr to hypre_csr. Allocate the data array and set it to all ones.
336  // Return the appropriate ownership flag. The CSR arrays are wrapped in the
337  // mem_csr struct which is used to move these arrays to device, if necessary.
338  static signed char CopyBoolCSR(Table *bool_csr,
339  MemoryIJData &mem_csr,
340  hypre_CSRMatrix *hypre_csr);
341 
342  // Copy the j array of a hypre_CSRMatrix to the given J array, converting
343  // the indices from HYPRE_Int/HYPRE_BigInt to int.
344  static void CopyCSR_J(hypre_CSRMatrix *hypre_csr, int *J);
345 
346  // Wrap the data from h_mat into mem with the given ownership flag.
347  // If the new Memory arrays in mem are not suitable to be accessed via
348  // GetHypreMemoryClass(), then mem will be re-allocated using the memory type
349  // returned by GetHypreMemoryType(), the data will be deep copied, and h_mat
350  // will be updated with the new pointers.
351  static signed char HypreCsrToMem(hypre_CSRMatrix *h_mat, MemoryType h_mat_mt,
352  bool own_ija, MemoryIJData &mem);
353 
354 public:
355  /// An empty matrix to be used as a reference to an existing matrix
356  HypreParMatrix();
357 
358  /// Converts hypre's format to HypreParMatrix
359  /** If @a owner is false, ownership of @a a is not transferred */
360  void WrapHypreParCSRMatrix(hypre_ParCSRMatrix *a, bool owner = true);
361 
362  /// Converts hypre's format to HypreParMatrix
363  /** If @a owner is false, ownership of @a a is not transferred */
364  explicit HypreParMatrix(hypre_ParCSRMatrix *a, bool owner = true)
365  {
366  Init();
367  WrapHypreParCSRMatrix(a, owner);
368  }
369 
370  /// Creates block-diagonal square parallel matrix.
371  /** Diagonal is given by @a diag which must be in CSR format (finalized). The
372  new HypreParMatrix does not take ownership of any of the input arrays.
373  See @ref hypre_partitioning_descr "here" for a description of the row
374  partitioning array @a row_starts.
375 
376  @warning The ordering of the columns in each row in @a *diag may be
377  changed by this constructor to ensure that the first entry in each row is
378  the diagonal one. This is expected by most hypre functions. */
379  HypreParMatrix(MPI_Comm comm, HYPRE_BigInt glob_size,
380  HYPRE_BigInt *row_starts,
381  SparseMatrix *diag); // constructor with 4 arguments, v1
382 
383  /// Creates block-diagonal rectangular parallel matrix.
384  /** Diagonal is given by @a diag which must be in CSR format (finalized). The
385  new HypreParMatrix does not take ownership of any of the input arrays.
386  See @ref hypre_partitioning_descr "here" for a description of the
387  partitioning arrays @a row_starts and @a col_starts. */
388  HypreParMatrix(MPI_Comm comm, HYPRE_BigInt global_num_rows,
389  HYPRE_BigInt global_num_cols, HYPRE_BigInt *row_starts,
390  HYPRE_BigInt *col_starts,
391  SparseMatrix *diag); // constructor with 6 arguments, v1
392 
393  /// Creates general (rectangular) parallel matrix.
394  /** The new HypreParMatrix does not take ownership of any of the input
395  arrays, if @a own_diag_offd is false (default). If @a own_diag_offd is
396  true, ownership of @a diag and @a offd is transferred to the
397  HypreParMatrix.
398 
399  See @ref hypre_partitioning_descr "here" for a description of the
400  partitioning arrays @a row_starts and @a col_starts. */
401  HypreParMatrix(MPI_Comm comm, HYPRE_BigInt global_num_rows,
402  HYPRE_BigInt global_num_cols, HYPRE_BigInt *row_starts,
403  HYPRE_BigInt *col_starts, SparseMatrix *diag,
404  SparseMatrix *offd, HYPRE_BigInt *cmap,
405  bool own_diag_offd = false); // constructor with 8+1 arguments
406 
407  /// Creates general (rectangular) parallel matrix.
408  /** The new HypreParMatrix takes ownership of all input arrays, except
409  @a col_starts and @a row_starts. See @ref hypre_partitioning_descr "here"
410  for a description of the partitioning arrays @a row_starts and @a
411  col_starts.
412 
413  If @a hypre_arrays is false, all arrays (except @a row_starts and
414  @a col_starts) are assumed to be allocated according to the MemoryType
415  returned by Device::GetHostMemoryType(). If @a hypre_arrays is true, then
416  the same arrays are assumed to be allocated by hypre as host arrays. */
417  HypreParMatrix(MPI_Comm comm,
418  HYPRE_BigInt global_num_rows, HYPRE_BigInt global_num_cols,
419  HYPRE_BigInt *row_starts, HYPRE_BigInt *col_starts,
420  HYPRE_Int *diag_i, HYPRE_Int *diag_j, double *diag_data,
421  HYPRE_Int *offd_i, HYPRE_Int *offd_j, double *offd_data,
422  HYPRE_Int offd_num_cols,
423  HYPRE_BigInt *offd_col_map,
424  bool hypre_arrays = false); // constructor with 13+1 arguments
425 
426  /// Creates a parallel matrix from SparseMatrix on processor 0.
427  /** See @ref hypre_partitioning_descr "here" for a description of the
428  partitioning arrays @a row_starts and @a col_starts. */
429  HypreParMatrix(MPI_Comm comm, HYPRE_BigInt *row_starts,
430  HYPRE_BigInt *col_starts,
431  SparseMatrix *a); // constructor with 4 arguments, v2
432 
433  /// Creates boolean block-diagonal rectangular parallel matrix.
434  /** The new HypreParMatrix does not take ownership of any of the input
435  arrays. See @ref hypre_partitioning_descr "here" for a description of the
436  partitioning arrays @a row_starts and @a col_starts. */
437  HypreParMatrix(MPI_Comm comm, HYPRE_BigInt global_num_rows,
438  HYPRE_BigInt global_num_cols, HYPRE_BigInt *row_starts,
439  HYPRE_BigInt *col_starts,
440  Table *diag); // constructor with 6 arguments, v2
441 
442  /// Creates boolean rectangular parallel matrix.
443  /** The new HypreParMatrix takes ownership of the arrays @a i_diag,
444  @a j_diag, @a i_offd, @a j_offd, and @a cmap; does not take ownership of
445  the arrays @a row and @a col. See @ref hypre_partitioning_descr "here"
446  for a description of the partitioning arrays @a row and @a col. */
447  HypreParMatrix(MPI_Comm comm, int id, int np, HYPRE_BigInt *row,
448  HYPRE_BigInt *col,
449  HYPRE_Int *i_diag, HYPRE_Int *j_diag, HYPRE_Int *i_offd,
450  HYPRE_Int *j_offd, HYPRE_BigInt *cmap,
451  HYPRE_Int cmap_size); // constructor with 11 arguments
452 
453  /** @brief Creates a general parallel matrix from a local CSR matrix on each
454  processor described by the @a I, @a J and @a data arrays. */
455  /** The local matrix should be of size (local) @a nrows by (global)
456  @a glob_ncols. The new parallel matrix contains copies of all input
457  arrays (so they can be deleted). See @ref hypre_partitioning_descr "here"
458  for a description of the partitioning arrays @a rows and @a cols. */
459  HypreParMatrix(MPI_Comm comm, int nrows, HYPRE_BigInt glob_nrows,
460  HYPRE_BigInt glob_ncols, int *I, HYPRE_BigInt *J,
461  double *data, HYPRE_BigInt *rows,
462  HYPRE_BigInt *cols); // constructor with 9 arguments
463 
464  /** @brief Copy constructor for a ParCSR matrix which creates a deep copy of
465  structure and data from @a P. */
466  HypreParMatrix(const HypreParMatrix &P);
467 
468  /// Make this HypreParMatrix a reference to 'master'
469  void MakeRef(const HypreParMatrix &master);
470 
471  /// MPI communicator
472  MPI_Comm GetComm() const { return A->comm; }
473 
474  /// Typecasting to hypre's hypre_ParCSRMatrix*
475  operator hypre_ParCSRMatrix*() const { return A; }
476 #ifndef HYPRE_PAR_CSR_MATRIX_STRUCT
477  /// Typecasting to hypre's HYPRE_ParCSRMatrix, a.k.a. void *
478  operator HYPRE_ParCSRMatrix() { return (HYPRE_ParCSRMatrix) A; }
479 #endif
480  /// Changes the ownership of the matrix
481  hypre_ParCSRMatrix* StealData();
482 
483  /// Explicitly set the three ownership flags, see docs for diagOwner etc.
484  void SetOwnerFlags(signed char diag, signed char offd, signed char colmap);
485 
486  /// Get diag ownership flag
487  signed char OwnsDiag() const { return diagOwner; }
488  /// Get offd ownership flag
489  signed char OwnsOffd() const { return offdOwner; }
490  /// Get colmap ownership flag
491  signed char OwnsColMap() const { return colMapOwner; }
492 
493  /** If the HypreParMatrix does not own the row-starts array, make a copy of
494  it that the HypreParMatrix will own. If the col-starts array is the same
495  as the row-starts array, col-starts is also replaced. */
496  void CopyRowStarts();
497  /** If the HypreParMatrix does not own the col-starts array, make a copy of
498  it that the HypreParMatrix will own. If the row-starts array is the same
499  as the col-starts array, row-starts is also replaced. */
500  void CopyColStarts();
501 
502  /// Returns the global number of nonzeros
503  inline HYPRE_BigInt NNZ() const { return A->num_nonzeros; }
504  /// Returns the row partitioning
505  /** See @ref hypre_partitioning_descr "here" for a description of the
506  partitioning array. */
507  inline HYPRE_BigInt *RowPart() { return A->row_starts; }
508  /// Returns the column partitioning
509  /** See @ref hypre_partitioning_descr "here" for a description of the
510  partitioning array. */
511  inline HYPRE_BigInt *ColPart() { return A->col_starts; }
512  /// Returns the row partitioning (const version)
513  /** See @ref hypre_partitioning_descr "here" for a description of the
514  partitioning array. */
515  inline const HYPRE_BigInt *RowPart() const { return A->row_starts; }
516  /// Returns the column partitioning (const version)
517  /** See @ref hypre_partitioning_descr "here" for a description of the
518  partitioning array. */
519  inline const HYPRE_BigInt *ColPart() const { return A->col_starts; }
520  /// Returns the global number of rows
521  inline HYPRE_BigInt M() const { return A->global_num_rows; }
522  /// Returns the global number of columns
523  inline HYPRE_BigInt N() const { return A->global_num_cols; }
524 
525  /// Get the local diagonal of the matrix.
526  void GetDiag(Vector &diag) const;
527  /// Get the local diagonal block. NOTE: 'diag' will not own any data.
528  void GetDiag(SparseMatrix &diag) const;
529  /// Get the local off-diagonal block. NOTE: 'offd' will not own any data.
530  void GetOffd(SparseMatrix &offd, HYPRE_BigInt* &cmap) const;
531  /** @brief Get a single SparseMatrix containing all rows from this processor,
532  merged from the diagonal and off-diagonal blocks stored by the
533  HypreParMatrix. */
534  /** @note The number of columns in the SparseMatrix will be the global number
535  of columns in the parallel matrix, so using this method may result in an
536  integer overflow in the column indices. */
537  void MergeDiagAndOffd(SparseMatrix &merged);
538 
539  /// Return the diagonal of the matrix (Operator interface).
540  virtual void AssembleDiagonal(Vector &diag) const { GetDiag(diag); }
541 
542  /** Split the matrix into M x N equally sized blocks of parallel matrices.
543  The size of 'blocks' must already be set to M x N. */
544  void GetBlocks(Array2D<HypreParMatrix*> &blocks,
545  bool interleaved_rows = false,
546  bool interleaved_cols = false) const;
547 
548  /// Returns the transpose of *this
549  HypreParMatrix * Transpose() const;
550 
551  /** Returns principle submatrix given by array of indices of connections
552  with relative size > @a threshold in *this. */
553 #if MFEM_HYPRE_VERSION >= 21800
555  double threshold=0.0) const;
556 #endif
557 
558  /// Returns the number of rows in the diagonal block of the ParCSRMatrix
559  int GetNumRows() const
560  {
561  return internal::to_int(
562  hypre_CSRMatrixNumRows(hypre_ParCSRMatrixDiag(A)));
563  }
564 
565  /// Returns the number of columns in the diagonal block of the ParCSRMatrix
566  int GetNumCols() const
567  {
568  return internal::to_int(
569  hypre_CSRMatrixNumCols(hypre_ParCSRMatrixDiag(A)));
570  }
571 
572  /// Return the global number of rows
574  { return hypre_ParCSRMatrixGlobalNumRows(A); }
575 
576  /// Return the global number of columns
578  { return hypre_ParCSRMatrixGlobalNumCols(A); }
579 
580  /// Return the parallel row partitioning array.
581  /** See @ref hypre_partitioning_descr "here" for a description of the
582  partitioning array. */
583  HYPRE_BigInt *GetRowStarts() const { return hypre_ParCSRMatrixRowStarts(A); }
584 
585  /// Return the parallel column partitioning array.
586  /** See @ref hypre_partitioning_descr "here" for a description of the
587  partitioning array. */
588  HYPRE_BigInt *GetColStarts() const { return hypre_ParCSRMatrixColStarts(A); }
589 
590  virtual MemoryClass GetMemoryClass() const { return GetHypreMemoryClass(); }
591 
592  /// Computes y = alpha * A * x + beta * y
593  HYPRE_Int Mult(HypreParVector &x, HypreParVector &y,
594  double alpha = 1.0, double beta = 0.0) const;
595  /// Computes y = alpha * A * x + beta * y
596  HYPRE_Int Mult(HYPRE_ParVector x, HYPRE_ParVector y,
597  double alpha = 1.0, double beta = 0.0) const;
598  /// Computes y = alpha * A^t * x + beta * y
600  double alpha = 1.0, double beta = 0.0) const;
601 
602  void Mult(double a, const Vector &x, double b, Vector &y) const;
603  void MultTranspose(double a, const Vector &x, double b, Vector &y) const;
604 
605  virtual void Mult(const Vector &x, Vector &y) const
606  { Mult(1.0, x, 0.0, y); }
607  virtual void MultTranspose(const Vector &x, Vector &y) const
608  { MultTranspose(1.0, x, 0.0, y); }
609 
610  /** @brief Computes y = a * |A| * x + b * y, using entry-wise absolute values
611  of the matrix A. */
612  void AbsMult(double a, const Vector &x, double b, Vector &y) const;
613 
614  /** @brief Computes y = a * |At| * x + b * y, using entry-wise absolute
615  values of the transpose of the matrix A. */
616  void AbsMultTranspose(double a, const Vector &x, double b, Vector &y) const;
617 
618  /** @brief The "Boolean" analog of y = alpha * A * x + beta * y, where
619  elements in the sparsity pattern of the matrix are treated as "true". */
620  void BooleanMult(int alpha, const int *x, int beta, int *y)
621  {
622  HostRead();
623  internal::hypre_ParCSRMatrixBooleanMatvec(A, alpha, const_cast<int*>(x),
624  beta, y);
625  HypreRead();
626  }
627 
628  /** @brief The "Boolean" analog of y = alpha * A^T * x + beta * y, where
629  elements in the sparsity pattern of the matrix are treated as "true". */
630  void BooleanMultTranspose(int alpha, const int *x, int beta, int *y)
631  {
632  HostRead();
633  internal::hypre_ParCSRMatrixBooleanMatvecT(A, alpha, const_cast<int*>(x),
634  beta, y);
635  HypreRead();
636  }
637 
638  /// Initialize all entries with value.
639  HypreParMatrix &operator=(double value)
640  {
641 #if MFEM_HYPRE_VERSION < 22200
642  internal::hypre_ParCSRMatrixSetConstantValues(A, value);
643 #else
644  hypre_ParCSRMatrixSetConstantValues(A, value);
645 #endif
646  return *this;
647  }
648 
649  /** Perform the operation `*this += B`, assuming that both matrices use the
650  same row and column partitions and the same col_map_offd arrays, or B has
651  an empty off-diagonal block. We also assume that the sparsity pattern of
652  `*this` contains that of `B`. */
653  HypreParMatrix &operator+=(const HypreParMatrix &B) { return Add(1.0, B); }
654 
655  /** Perform the operation `*this += beta*B`, assuming that both matrices use
656  the same row and column partitions and the same col_map_offd arrays, or
657  B has an empty off-diagonal block. We also assume that the sparsity
658  pattern of `*this` contains that of `B`. For a more general case consider
659  the stand-alone function ParAdd described below. */
660  HypreParMatrix &Add(const double beta, const HypreParMatrix &B)
661  {
662  MFEM_VERIFY(internal::hypre_ParCSRMatrixSum(A, beta, B.A) == 0,
663  "error in hypre_ParCSRMatrixSum");
664  return *this;
665  }
666 
667  /** @brief Multiply the HypreParMatrix on the left by a block-diagonal
668  parallel matrix @a D and return the result as a new HypreParMatrix. */
669  /** If @a D has a different number of rows than @a A (this matrix), @a D's
670  row starts array needs to be given (as returned by the methods
671  GetDofOffsets/GetTrueDofOffsets of ParFiniteElementSpace). The new
672  matrix @a D*A uses copies of the row- and column-starts arrays, so "this"
673  matrix and @a row_starts can be deleted.
674  @note This operation is local and does not require communication. */
676  HYPRE_BigInt* row_starts = NULL) const;
677 
678  /// Scale the local row i by s(i).
679  void ScaleRows(const Vector & s);
680  /// Scale the local row i by 1./s(i)
681  void InvScaleRows(const Vector & s);
682  /// Scale all entries by s: A_scaled = s*A.
683  void operator*=(double s);
684 
685  /// Remove values smaller in absolute value than some threshold
686  void Threshold(double threshold = 0.0);
687 
688  /** @brief Wrapper for hypre_ParCSRMatrixDropSmallEntries in different
689  versions of hypre. Drop off-diagonal entries that are smaller than
690  tol * l2 norm of its row */
691  /** For HYPRE versions < 2.14, this method just calls Threshold() with
692  threshold = tol * max(l2 row norm). */
693  void DropSmallEntries(double tol);
694 
695  /// If a row contains only zeros, set its diagonal to 1.
696  void EliminateZeroRows() { hypre_ParCSRMatrixFixZeroRows(A); }
697 
698  /** Eliminate rows and columns from the matrix, and rows from the vector B.
699  Modify B with the BC values in X. */
700  void EliminateRowsCols(const Array<int> &rows_cols, const HypreParVector &X,
701  HypreParVector &B);
702 
703  /** Eliminate rows and columns from the matrix and store the eliminated
704  elements in a new matrix Ae (returned), so that the modified matrix and
705  Ae sum to the original matrix. */
706  HypreParMatrix* EliminateRowsCols(const Array<int> &rows_cols);
707 
708  /** Eliminate columns from the matrix and store the eliminated elements in a
709  new matrix Ae (returned) so that the modified matrix and Ae sum to the
710  original matrix. */
712 
713  /// Eliminate rows from the diagonal and off-diagonal blocks of the matrix.
714  void EliminateRows(const Array<int> &rows);
715 
716  /** @brief Eliminate essential BC specified by @a ess_dof_list from the
717  solution @a X to the r.h.s. @a B. */
718  /** This matrix is the matrix with eliminated BC, while @a Ae is such that
719  (A+Ae) is the original (Neumann) matrix before elimination. */
720  void EliminateBC(const HypreParMatrix &Ae, const Array<int> &ess_dof_list,
721  const Vector &X, Vector &B) const;
722 
723  /// Update the internal hypre_ParCSRMatrix object, A, to be on host.
724  /** After this call A's diagonal and off-diagonal should not be modified
725  until after a suitable call to {Host,Hypre}{Write,ReadWrite}. */
726  void HostRead() const { Read(Device::GetHostMemoryClass()); }
727 
728  /// Update the internal hypre_ParCSRMatrix object, A, to be on host.
729  /** After this call A's diagonal and off-diagonal can be modified on host
730  and subsequent calls to Hypre{Read,Write,ReadWrite} will require a deep
731  copy of the data if hypre is built with device support. */
732  void HostReadWrite() { ReadWrite(Device::GetHostMemoryClass()); }
733 
734  /// Update the internal hypre_ParCSRMatrix object, A, to be on host.
735  /** Similar to HostReadWrite(), except that the data will never be copied
736  from device to host to ensure host contains the correct current data. */
738 
739  /** @brief Update the internal hypre_ParCSRMatrix object, A, to be in hypre
740  memory space. */
741  /** After this call A's diagonal and off-diagonal should not be modified
742  until after a suitable call to {Host,Hypre}{Write,ReadWrite}. */
743  void HypreRead() const { Read(GetHypreMemoryClass()); }
744 
745  /** @brief Update the internal hypre_ParCSRMatrix object, A, to be in hypre
746  memory space. */
747  /** After this call A's diagonal and off-diagonal can be modified in hypre
748  memory space and subsequent calls to Host{Read,Write,ReadWrite} will
749  require a deep copy of the data if hypre is built with device support. */
750  void HypreReadWrite() { ReadWrite(GetHypreMemoryClass()); }
751 
752  /** @brief Update the internal hypre_ParCSRMatrix object, A, to be in hypre
753  memory space. */
754  /** Similar to HostReadWrite(), except that the data will never be copied
755  from host to hypre memory space to ensure the latter contains the correct
756  current data. */
757  void HypreWrite() { Write(GetHypreMemoryClass()); }
758 
759  /// Prints the locally owned rows in parallel
760  void Print(const char *fname, HYPRE_Int offi = 0, HYPRE_Int offj = 0) const;
761  /// Reads the matrix from a file
762  void Read(MPI_Comm comm, const char *fname);
763  /// Read a matrix saved as a HYPRE_IJMatrix
764  void Read_IJMatrix(MPI_Comm comm, const char *fname);
765 
766  /// Print information about the hypre_ParCSRCommPkg of the HypreParMatrix.
767  void PrintCommPkg(std::ostream &out = mfem::out) const;
768 
769  /** @brief Print sizes and hashes for all data arrays of the HypreParMatrix
770  from the local MPI rank. */
771  /** This is a compact text representation of the local data of the
772  HypreParMatrix that can be used to compare matrices from different runs
773  without the need to save the whole matrix. */
774  void PrintHash(std::ostream &out) const;
775 
776  /// Calls hypre's destroy function
777  virtual ~HypreParMatrix() { Destroy(); }
778 
779  Type GetType() const { return Hypre_ParCSR; }
780 };
781 
782 #if MFEM_HYPRE_VERSION >= 21800
783 
785 {
786  MATRIX_ONLY,
787  RHS_ONLY,
789 };
790 
791 /** Constructs and applies block diagonal inverse of HypreParMatrix.
792  The enum @a job specifies whether the matrix or the RHS should be
793  scaled (or both). */
794 void BlockInverseScale(const HypreParMatrix *A, HypreParMatrix *C,
795  const Vector *b, HypreParVector *d,
796  int blocksize, BlockInverseScaleJob job);
797 #endif
798 
799 /** @brief Return a new matrix `C = alpha*A + beta*B`, assuming that both `A`
800  and `B` use the same row and column partitions and the same `col_map_offd`
801  arrays. */
802 HypreParMatrix *Add(double alpha, const HypreParMatrix &A,
803  double beta, const HypreParMatrix &B);
804 
805 /** Returns the matrix @a A * @a B. Returned matrix does not necessarily own
806  row or column starts unless the bool @a own_matrix is set to true. */
807 HypreParMatrix * ParMult(const HypreParMatrix *A, const HypreParMatrix *B,
808  bool own_matrix = false);
809 /// Returns the matrix A + B
810 /** It is assumed that both matrices use the same row and column partitions and
811  the same col_map_offd arrays. */
812 HypreParMatrix * ParAdd(const HypreParMatrix *A, const HypreParMatrix *B);
813 
814 /// Returns the matrix P^t * A * P
815 HypreParMatrix * RAP(const HypreParMatrix *A, const HypreParMatrix *P);
816 /// Returns the matrix Rt^t * A * P
817 HypreParMatrix * RAP(const HypreParMatrix * Rt, const HypreParMatrix *A,
818  const HypreParMatrix *P);
819 
820 /// Returns a merged hypre matrix constructed from hypre matrix blocks.
821 /** It is assumed that all block matrices use the same communicator, and the
822  block sizes are consistent in rows and columns. Rows and columns are
823  renumbered but not redistributed in parallel, e.g. the block rows owned by
824  each process remain on that process in the resulting matrix. Some blocks can
825  be NULL. Each block and the entire system can be rectangular. Scalability to
826  extremely large processor counts is limited by global MPI communication, see
827  GatherBlockOffsetData() in hypre.cpp. */
828 HypreParMatrix * HypreParMatrixFromBlocks(Array2D<HypreParMatrix*> &blocks,
829  Array2D<double> *blockCoeff=NULL);
830 
831 /** @brief Eliminate essential BC specified by @a ess_dof_list from the solution
832  @a X to the r.h.s. @a B. */
833 /** Here @a A is a matrix with eliminated BC, while @a Ae is such that (A+Ae) is
834  the original (Neumann) matrix before elimination. */
835 void EliminateBC(const HypreParMatrix &A, const HypreParMatrix &Ae,
836  const Array<int> &ess_dof_list, const Vector &X, Vector &B);
837 
838 
839 /// Parallel smoothers in hypre
840 class HypreSmoother : public Solver
841 {
842 protected:
843  /// The linear system matrix
845  /// Right-hand side and solution vectors
846  mutable HypreParVector *B, *X;
847  /** @brief Auxiliary buffers for the case when the input or output arrays in
848  methods like Mult(const Vector &, Vector &) need to be deep copied in
849  order to be used by hypre. */
851  /// Temporary vectors
852  mutable HypreParVector *V, *Z;
853  /// FIR Filter Temporary Vectors
854  mutable HypreParVector *X0, *X1;
855 
856  /** Smoother type from hypre_ParCSRRelax() in ams.c plus extensions, see the
857  enumeration Type below. */
858  int type;
859  /// Number of relaxation sweeps
861  /// Damping coefficient (usually <= 1)
862  double relax_weight;
863  /// SOR parameter (usually in (0,2))
864  double omega;
865  /// Order of the smoothing polynomial
867  /// Fraction of spectrum to smooth for polynomial relaxation
869  /// Apply the polynomial smoother to A or D^{-1/2} A D^{-1/2}
871 
872  /// Taubin's lambda-mu method parameters
873  double lambda;
874  double mu;
876 
877  /// l1 norms of the rows of A
878  double *l1_norms;
879  /// If set, take absolute values of the computed l1_norms
881  /// Number of CG iterations to determine eigenvalue estimates
883  /// Maximal eigenvalue estimate for polynomial smoothing
884  double max_eig_est;
885  /// Minimal eigenvalue estimate for polynomial smoothing
886  double min_eig_est;
887  /// Parameters for windowing function of FIR filter
888  double window_params[3];
889 
890  /// Combined coefficients for windowing and Chebyshev polynomials.
891  double* fir_coeffs;
892 
893  /// A flag that indicates whether the linear system matrix A is symmetric
895 
896 public:
897  /** Hypre smoother types:
898  0 = Jacobi
899  1 = l1-scaled Jacobi
900  2 = l1-scaled block Gauss-Seidel/SSOR
901  4 = truncated l1-scaled block Gauss-Seidel/SSOR
902  5 = lumped Jacobi
903  6 = Gauss-Seidel
904  10 = On-processor forward solve for matrix w/ triangular structure
905  16 = Chebyshev
906  1001 = Taubin polynomial smoother
907  1002 = FIR polynomial smoother. */
908  enum Type { Jacobi = 0, l1Jacobi = 1, l1GS = 2, l1GStr = 4, lumpedJacobi = 5,
909  GS = 6, OPFS = 10, Chebyshev = 16, Taubin = 1001, FIR = 1002
910  };
911 #ifndef HYPRE_USING_CUDA
912  static constexpr Type default_type = l1GS;
913 #else
914  static constexpr Type default_type = l1Jacobi;
915 #endif
916 
917  HypreSmoother();
918 
920  int relax_times = 1, double relax_weight = 1.0,
921  double omega = 1.0, int poly_order = 2,
922  double poly_fraction = .3, int eig_est_cg_iter = 10);
923 
924  /// Set the relaxation type and number of sweeps
926  /// Set SOR-related parameters
927  void SetSOROptions(double relax_weight, double omega);
928  /// Set parameters for polynomial smoothing
929  /** By default, 10 iterations of CG are used to estimate the eigenvalues.
930  Setting eig_est_cg_iter = 0 uses hypre's hypre_ParCSRMaxEigEstimate() instead. */
931  void SetPolyOptions(int poly_order, double poly_fraction,
932  int eig_est_cg_iter = 10);
933  /// Set parameters for Taubin's lambda-mu method
934  void SetTaubinOptions(double lambda, double mu, int iter);
935 
936  /// Convenience function for setting canonical windowing parameters
937  void SetWindowByName(const char* window_name);
938  /// Set parameters for windowing function for FIR smoother.
939  void SetWindowParameters(double a, double b, double c);
940  /// Compute window and Chebyshev coefficients for given polynomial order.
941  void SetFIRCoefficients(double max_eig);
942 
943  /// After computing l1-norms, replace them with their absolute values.
944  /** By default, the l1-norms take their sign from the corresponding diagonal
945  entries in the associated matrix. */
946  void SetPositiveDiagonal(bool pos = true) { pos_l1_norms = pos; }
947 
948  /** Explicitly indicate whether the linear system matrix A is symmetric. If A
949  is symmetric, the smoother will also be symmetric. In this case, calling
950  MultTranspose will be redirected to Mult. (This is also done if the
951  smoother is diagonal.) By default, A is assumed to be nonsymmetric. */
952  void SetOperatorSymmetry(bool is_sym) { A_is_symmetric = is_sym; }
953 
954  /** Set/update the associated operator. Must be called after setting the
955  HypreSmoother type and options. */
956  virtual void SetOperator(const Operator &op);
957 
958  /// Relax the linear system Ax=b
959  virtual void Mult(const HypreParVector &b, HypreParVector &x) const;
960  virtual void Mult(const Vector &b, Vector &x) const;
961 
962  /// Apply transpose of the smoother to relax the linear system Ax=b
963  virtual void MultTranspose(const Vector &b, Vector &x) const;
964 
965  virtual ~HypreSmoother();
966 };
967 
968 
969 /// Abstract class for hypre's solvers and preconditioners
970 class HypreSolver : public Solver
971 {
972 public:
973  /// How to treat errors returned by hypre function calls.
975  {
976  IGNORE_HYPRE_ERRORS, ///< Ignore hypre errors (see e.g. HypreADS)
977  WARN_HYPRE_ERRORS, ///< Issue warnings on hypre errors
978  ABORT_HYPRE_ERRORS ///< Abort on hypre errors (default in base class)
979  };
980 
981 protected:
982  /// The linear system matrix
984 
985  /// Right-hand side and solution vector
986  mutable HypreParVector *B, *X;
987 
989 
990  /// Was hypre's Setup function called already?
991  mutable int setup_called;
992 
993  /// How to treat hypre errors.
995 
996 public:
997  HypreSolver();
998 
999  HypreSolver(const HypreParMatrix *A_);
1000 
1001  /// Typecast to HYPRE_Solver -- return the solver
1002  virtual operator HYPRE_Solver() const = 0;
1003 
1004  /// hypre's internal Setup function
1005  virtual HYPRE_PtrToParSolverFcn SetupFcn() const = 0;
1006  /// hypre's internal Solve function
1007  virtual HYPRE_PtrToParSolverFcn SolveFcn() const = 0;
1008 
1009  virtual void SetOperator(const Operator &op)
1010  { mfem_error("HypreSolvers do not support SetOperator!"); }
1011 
1012  virtual MemoryClass GetMemoryClass() const { return GetHypreMemoryClass(); }
1013 
1014  /// Solve the linear system Ax=b
1015  virtual void Mult(const HypreParVector &b, HypreParVector &x) const;
1016  virtual void Mult(const Vector &b, Vector &x) const;
1017 
1018  /** @brief Set the behavior for treating hypre errors, see the ErrorMode
1019  enum. The default mode in the base class is ABORT_HYPRE_ERRORS. */
1020  /** Currently, there are three cases in derived classes where the error flag
1021  is set to IGNORE_HYPRE_ERRORS:
1022  * in the method HypreBoomerAMG::SetElasticityOptions(), and
1023  * in the constructor of classes HypreAMS and HypreADS.
1024  The reason for this is that a nonzero hypre error is returned) when
1025  hypre_ParCSRComputeL1Norms() encounters zero row in a matrix, which is
1026  expected in some cases with the above solvers. */
1027  void SetErrorMode(ErrorMode err_mode) const { error_mode = err_mode; }
1028 
1029  virtual ~HypreSolver();
1030 };
1031 
1032 
1033 #if MFEM_HYPRE_VERSION >= 21800
1034 /** Preconditioner for HypreParMatrices that are triangular in some ordering.
1035  Finds correct ordering and performs forward substitution on processor
1036  as approximate inverse. Exact on one processor. */
1038 {
1039 public:
1041  explicit HypreTriSolve(const HypreParMatrix &A) : HypreSolver(&A) { }
1042  virtual operator HYPRE_Solver() const { return NULL; }
1043 
1044  virtual HYPRE_PtrToParSolverFcn SetupFcn() const
1045  { return (HYPRE_PtrToParSolverFcn) HYPRE_ParCSROnProcTriSetup; }
1046  virtual HYPRE_PtrToParSolverFcn SolveFcn() const
1047  { return (HYPRE_PtrToParSolverFcn) HYPRE_ParCSROnProcTriSolve; }
1048 
1049  const HypreParMatrix* GetData() const { return A; }
1050 
1051  /// Deprecated. Use HypreTriSolve::GetData() const instead.
1052  MFEM_DEPRECATED HypreParMatrix* GetData()
1053  { return const_cast<HypreParMatrix*>(A); }
1054 
1055  virtual ~HypreTriSolve() { }
1056 };
1057 #endif
1058 
1059 /// PCG solver in hypre
1060 class HyprePCG : public HypreSolver
1061 {
1062 private:
1063  HYPRE_Solver pcg_solver;
1064 
1065  HypreSolver * precond;
1066 
1067 public:
1068  HyprePCG(MPI_Comm comm);
1069 
1070  HyprePCG(const HypreParMatrix &A_);
1071 
1072  virtual void SetOperator(const Operator &op);
1073 
1074  void SetTol(double tol);
1075  void SetAbsTol(double atol);
1076  void SetMaxIter(int max_iter);
1077  void SetLogging(int logging);
1078  void SetPrintLevel(int print_lvl);
1079 
1080  /// Set the hypre solver to be used as a preconditioner
1081  void SetPreconditioner(HypreSolver &precond);
1082 
1083  /** Use the L2 norm of the residual for measuring PCG convergence, plus
1084  (optionally) 1) periodically recompute true residuals from scratch; and
1085  2) enable residual-based stopping criteria. */
1086  void SetResidualConvergenceOptions(int res_frequency=-1, double rtol=0.0);
1087 
1088  /// deprecated: use SetZeroInitialIterate()
1089  MFEM_DEPRECATED void SetZeroInintialIterate() { iterative_mode = false; }
1090 
1091  /// non-hypre setting
1093 
1094  void GetNumIterations(int &num_iterations) const
1095  {
1096  HYPRE_Int num_it;
1097  HYPRE_ParCSRPCGGetNumIterations(pcg_solver, &num_it);
1098  num_iterations = internal::to_int(num_it);
1099  }
1100 
1101  /// The typecast to HYPRE_Solver returns the internal pcg_solver
1102  virtual operator HYPRE_Solver() const { return pcg_solver; }
1103 
1104  /// PCG Setup function
1105  virtual HYPRE_PtrToParSolverFcn SetupFcn() const
1106  { return (HYPRE_PtrToParSolverFcn) HYPRE_ParCSRPCGSetup; }
1107  /// PCG Solve function
1108  virtual HYPRE_PtrToParSolverFcn SolveFcn() const
1109  { return (HYPRE_PtrToParSolverFcn) HYPRE_ParCSRPCGSolve; }
1110 
1111  /// Solve Ax=b with hypre's PCG
1112  virtual void Mult(const HypreParVector &b, HypreParVector &x) const;
1113  using HypreSolver::Mult;
1114 
1115  virtual ~HyprePCG();
1116 };
1117 
1118 /// GMRES solver in hypre
1119 class HypreGMRES : public HypreSolver
1120 {
1121 private:
1122  HYPRE_Solver gmres_solver;
1123 
1124  HypreSolver * precond;
1125 
1126  /// Default, generally robust, GMRES options
1127  void SetDefaultOptions();
1128 
1129 public:
1130  HypreGMRES(MPI_Comm comm);
1131 
1132  HypreGMRES(const HypreParMatrix &A_);
1133 
1134  virtual void SetOperator(const Operator &op);
1135 
1136  void SetTol(double tol);
1137  void SetAbsTol(double tol);
1138  void SetMaxIter(int max_iter);
1139  void SetKDim(int dim);
1140  void SetLogging(int logging);
1141  void SetPrintLevel(int print_lvl);
1142 
1143  /// Set the hypre solver to be used as a preconditioner
1144  void SetPreconditioner(HypreSolver &precond);
1145 
1146  /// deprecated: use SetZeroInitialIterate()
1147  MFEM_DEPRECATED void SetZeroInintialIterate() { iterative_mode = false; }
1148 
1149  /// non-hypre setting
1151 
1152  /// The typecast to HYPRE_Solver returns the internal gmres_solver
1153  virtual operator HYPRE_Solver() const { return gmres_solver; }
1154 
1155  /// GMRES Setup function
1156  virtual HYPRE_PtrToParSolverFcn SetupFcn() const
1157  { return (HYPRE_PtrToParSolverFcn) HYPRE_ParCSRGMRESSetup; }
1158  /// GMRES Solve function
1159  virtual HYPRE_PtrToParSolverFcn SolveFcn() const
1160  { return (HYPRE_PtrToParSolverFcn) HYPRE_ParCSRGMRESSolve; }
1161 
1162  /// Solve Ax=b with hypre's GMRES
1163  virtual void Mult (const HypreParVector &b, HypreParVector &x) const;
1164  using HypreSolver::Mult;
1165 
1166  virtual ~HypreGMRES();
1167 };
1168 
1169 /// Flexible GMRES solver in hypre
1170 class HypreFGMRES : public HypreSolver
1171 {
1172 private:
1173  HYPRE_Solver fgmres_solver;
1174 
1175  HypreSolver * precond;
1176 
1177  /// Default, generally robust, FGMRES options
1178  void SetDefaultOptions();
1179 
1180 public:
1181  HypreFGMRES(MPI_Comm comm);
1182 
1183  HypreFGMRES(const HypreParMatrix &A_);
1184 
1185  virtual void SetOperator(const Operator &op);
1186 
1187  void SetTol(double tol);
1188  void SetMaxIter(int max_iter);
1189  void SetKDim(int dim);
1190  void SetLogging(int logging);
1191  void SetPrintLevel(int print_lvl);
1192 
1193  /// Set the hypre solver to be used as a preconditioner
1194  void SetPreconditioner(HypreSolver &precond);
1195 
1196  /// deprecated: use SetZeroInitialIterate()
1197  MFEM_DEPRECATED void SetZeroInintialIterate() { iterative_mode = false; }
1198 
1199  /// non-hypre setting
1201 
1202  /// The typecast to HYPRE_Solver returns the internal fgmres_solver
1203  virtual operator HYPRE_Solver() const { return fgmres_solver; }
1204 
1205  /// FGMRES Setup function
1206  virtual HYPRE_PtrToParSolverFcn SetupFcn() const
1207  { return (HYPRE_PtrToParSolverFcn) HYPRE_ParCSRFlexGMRESSetup; }
1208  /// FGMRES Solve function
1209  virtual HYPRE_PtrToParSolverFcn SolveFcn() const
1210  { return (HYPRE_PtrToParSolverFcn) HYPRE_ParCSRFlexGMRESSolve; }
1211 
1212  /// Solve Ax=b with hypre's FGMRES
1213  virtual void Mult (const HypreParVector &b, HypreParVector &x) const;
1214  using HypreSolver::Mult;
1215 
1216  virtual ~HypreFGMRES();
1217 };
1218 
1219 /// The identity operator as a hypre solver
1221 {
1222 public:
1223  virtual operator HYPRE_Solver() const { return NULL; }
1224 
1225  virtual HYPRE_PtrToParSolverFcn SetupFcn() const
1226  { return (HYPRE_PtrToParSolverFcn) hypre_ParKrylovIdentitySetup; }
1227  virtual HYPRE_PtrToParSolverFcn SolveFcn() const
1228  { return (HYPRE_PtrToParSolverFcn) hypre_ParKrylovIdentity; }
1229 
1230  virtual ~HypreIdentity() { }
1231 };
1232 
1233 /// Jacobi preconditioner in hypre
1235 {
1236 public:
1238  explicit HypreDiagScale(const HypreParMatrix &A) : HypreSolver(&A) { }
1239  virtual operator HYPRE_Solver() const { return NULL; }
1240 
1241  virtual void SetOperator(const Operator &op);
1242 
1243  virtual HYPRE_PtrToParSolverFcn SetupFcn() const
1244  { return (HYPRE_PtrToParSolverFcn) HYPRE_ParCSRDiagScaleSetup; }
1245  virtual HYPRE_PtrToParSolverFcn SolveFcn() const
1246  { return (HYPRE_PtrToParSolverFcn) HYPRE_ParCSRDiagScale; }
1247 
1248  const HypreParMatrix* GetData() const { return A; }
1249 
1250  /// Deprecated. Use HypreDiagScale::GetData() const instead.
1251  MFEM_DEPRECATED HypreParMatrix* GetData()
1252  { return const_cast<HypreParMatrix*>(A); }
1253 
1254  virtual ~HypreDiagScale() { }
1255 };
1256 
1257 /// The ParaSails preconditioner in hypre
1259 {
1260 private:
1261  HYPRE_Solver sai_precond;
1262 
1263  /// Default, generally robust, ParaSails options
1264  void SetDefaultOptions();
1265 
1266  // If sai_precond is NULL, this method allocates it and sets default options.
1267  // Otherwise the method saves the options from sai_precond, destroys it,
1268  // allocates a new object, and sets its options to the saved values.
1269  void ResetSAIPrecond(MPI_Comm comm);
1270 
1271 public:
1272  HypreParaSails(MPI_Comm comm);
1273 
1275 
1276  virtual void SetOperator(const Operator &op);
1277 
1278  void SetSymmetry(int sym);
1279 
1280  /// The typecast to HYPRE_Solver returns the internal sai_precond
1281  virtual operator HYPRE_Solver() const { return sai_precond; }
1282 
1283  virtual HYPRE_PtrToParSolverFcn SetupFcn() const
1284  { return (HYPRE_PtrToParSolverFcn) HYPRE_ParaSailsSetup; }
1285  virtual HYPRE_PtrToParSolverFcn SolveFcn() const
1286  { return (HYPRE_PtrToParSolverFcn) HYPRE_ParaSailsSolve; }
1287 
1288  virtual ~HypreParaSails();
1289 };
1290 
1291 /** The Euclid preconditioner in Hypre
1292 
1293  Euclid implements the Parallel Incomplete LU factorization technique. For
1294  more information see:
1295 
1296  "A Scalable Parallel Algorithm for Incomplete Factor Preconditioning" by
1297  David Hysom and Alex Pothen, https://doi.org/10.1137/S1064827500376193
1298 */
1299 class HypreEuclid : public HypreSolver
1300 {
1301 private:
1302  HYPRE_Solver euc_precond;
1303 
1304  /// Default, generally robust, Euclid options
1305  void SetDefaultOptions();
1306 
1307  // If euc_precond is NULL, this method allocates it and sets default options.
1308  // Otherwise the method saves the options from euc_precond, destroys it,
1309  // allocates a new object, and sets its options to the saved values.
1310  void ResetEuclidPrecond(MPI_Comm comm);
1311 
1312 public:
1313  HypreEuclid(MPI_Comm comm);
1314 
1315  HypreEuclid(const HypreParMatrix &A);
1316 
1317  virtual void SetOperator(const Operator &op);
1318 
1319  /// The typecast to HYPRE_Solver returns the internal euc_precond
1320  virtual operator HYPRE_Solver() const { return euc_precond; }
1321 
1322  virtual HYPRE_PtrToParSolverFcn SetupFcn() const
1323  { return (HYPRE_PtrToParSolverFcn) HYPRE_EuclidSetup; }
1324  virtual HYPRE_PtrToParSolverFcn SolveFcn() const
1325  { return (HYPRE_PtrToParSolverFcn) HYPRE_EuclidSolve; }
1326 
1327  virtual ~HypreEuclid();
1328 };
1329 
1330 #if MFEM_HYPRE_VERSION >= 21900
1331 /**
1332 @brief Wrapper for Hypre's native parallel ILU preconditioner.
1333 
1334 The default ILU factorization type is ILU(k). If you need to change this, or
1335 any other option, you can use the HYPRE_Solver method to cast the object for use
1336 with Hypre's native functions. For example, if want to use natural ordering
1337 rather than RCM reordering, you can use the following approach:
1338 
1339 @code
1340 mfem::HypreILU ilu();
1341 int reorder_type = 0;
1342 HYPRE_ILUSetLocalReordering(ilu, reorder_type);
1343 @endcode
1344 */
1345 class HypreILU : public HypreSolver
1346 {
1347 private:
1348  HYPRE_Solver ilu_precond;
1349 
1350  /// Set the ILU default options
1351  void SetDefaultOptions();
1352 
1353  /** Reset the ILU preconditioner.
1354  @note If ilu_precond is NULL, this method allocates; otherwise it destroys
1355  ilu_precond and allocates a new object. In both cases the default options
1356  are set. */
1357  void ResetILUPrecond();
1358 
1359 public:
1360  /// Constructor; sets the default options
1361  HypreILU();
1362 
1363  virtual ~HypreILU();
1364 
1365  /// Set the fill level for ILU(k); the default is k=1.
1366  void SetLevelOfFill(HYPRE_Int lev_fill);
1367 
1368  /// Set the print level: 0 = none, 1 = setup, 2 = solve, 3 = setup+solve
1369  void SetPrintLevel(HYPRE_Int print_level);
1370 
1371  /// The typecast to HYPRE_Solver returns the internal ilu_precond
1372  virtual operator HYPRE_Solver() const { return ilu_precond; }
1373 
1374  virtual void SetOperator(const Operator &op);
1375 
1376  /// ILU Setup function
1377  virtual HYPRE_PtrToParSolverFcn SetupFcn() const
1378  { return (HYPRE_PtrToParSolverFcn) HYPRE_ILUSetup; }
1379 
1380  /// ILU Solve function
1381  virtual HYPRE_PtrToParSolverFcn SolveFcn() const
1382  { return (HYPRE_PtrToParSolverFcn) HYPRE_ILUSolve; }
1383 };
1384 #endif
1385 
1386 /// The BoomerAMG solver in hypre
1388 {
1389 private:
1390  HYPRE_Solver amg_precond;
1391 
1392  /// Rigid body modes
1394 
1395  /// Finite element space for elasticity problems, see SetElasticityOptions()
1396  ParFiniteElementSpace *fespace;
1397 
1398  /// Recompute the rigid-body modes vectors (in the rbms array)
1399  void RecomputeRBMs();
1400 
1401  /// Default, generally robust, BoomerAMG options
1402  void SetDefaultOptions();
1403 
1404  // If amg_precond is NULL, allocates it and sets default options.
1405  // Otherwise saves the options from amg_precond, destroys it, allocates a new
1406  // one, and sets its options to the saved values.
1407  void ResetAMGPrecond();
1408 
1409 public:
1410  HypreBoomerAMG();
1411 
1413 
1414  virtual void SetOperator(const Operator &op);
1415 
1416  /** More robust options for systems, such as elasticity. */
1417  void SetSystemsOptions(int dim, bool order_bynodes=false);
1418 
1419  /** A special elasticity version of BoomerAMG that takes advantage of
1420  geometric rigid body modes and could perform better on some problems, see
1421  "Improving algebraic multigrid interpolation operators for linear
1422  elasticity problems", Baker, Kolev, Yang, NLAA 2009, DOI:10.1002/nla.688.
1423  This solver assumes Ordering::byVDIM in the FiniteElementSpace used to
1424  construct A. */
1426 
1427 #if MFEM_HYPRE_VERSION >= 21800
1428  /** Hypre parameters to use AIR AMG solve for advection-dominated problems.
1429  See "Nonsymmetric Algebraic Multigrid Based on Local Approximate Ideal
1430  Restriction (AIR)," Manteuffel, Ruge, Southworth, SISC (2018),
1431  DOI:/10.1137/17M1144350. Options: "distanceR" -> distance of neighbor
1432  DOFs for the restriction operator; options include 1, 2, and 15 (1.5).
1433  Strings "prerelax" and "postrelax" indicate points to relax on:
1434  F = F-points, C = C-points, A = all points. E.g., FFC -> relax on
1435  F-points, relax again on F-points, then relax on C-points. */
1436  void SetAdvectiveOptions(int distance=15, const std::string &prerelax="",
1437  const std::string &postrelax="FFC");
1438 
1439  /// Expert option - consult hypre documentation/team
1440  void SetStrongThresholdR(double strengthR)
1441  { HYPRE_BoomerAMGSetStrongThresholdR(amg_precond, strengthR); }
1442 
1443  /// Expert option - consult hypre documentation/team
1444  void SetFilterThresholdR(double filterR)
1445  { HYPRE_BoomerAMGSetFilterThresholdR(amg_precond, filterR); }
1446 
1447  /// Expert option - consult hypre documentation/team
1448  void SetRestriction(int restrict_type)
1449  { HYPRE_BoomerAMGSetRestriction(amg_precond, restrict_type); }
1450 
1451  /// Expert option - consult hypre documentation/team
1453  { HYPRE_BoomerAMGSetIsTriangular(amg_precond, 1); }
1454 
1455  /// Expert option - consult hypre documentation/team
1456  void SetGMRESSwitchR(int gmres_switch)
1457  { HYPRE_BoomerAMGSetGMRESSwitchR(amg_precond, gmres_switch); }
1458 
1459  /// Expert option - consult hypre documentation/team
1460  void SetCycleNumSweeps(int prerelax, int postrelax)
1461  {
1462  HYPRE_BoomerAMGSetCycleNumSweeps(amg_precond, prerelax, 1);
1463  HYPRE_BoomerAMGSetCycleNumSweeps(amg_precond, postrelax, 2);
1464  }
1465 #endif
1466 
1467  void SetPrintLevel(int print_level)
1468  { HYPRE_BoomerAMGSetPrintLevel(amg_precond, print_level); }
1469 
1470  void SetMaxIter(int max_iter)
1471  { HYPRE_BoomerAMGSetMaxIter(amg_precond, max_iter); }
1472 
1473  /// Expert option - consult hypre documentation/team
1474  void SetMaxLevels(int max_levels)
1475  { HYPRE_BoomerAMGSetMaxLevels(amg_precond, max_levels); }
1476 
1477  /// Expert option - consult hypre documentation/team
1478  void SetTol(double tol)
1479  { HYPRE_BoomerAMGSetTol(amg_precond, tol); }
1480 
1481  /// Expert option - consult hypre documentation/team
1482  void SetStrengthThresh(double strength)
1483  { HYPRE_BoomerAMGSetStrongThreshold(amg_precond, strength); }
1484 
1485  /// Expert option - consult hypre documentation/team
1486  void SetInterpolation(int interp_type)
1487  { HYPRE_BoomerAMGSetInterpType(amg_precond, interp_type); }
1488 
1489  /// Expert option - consult hypre documentation/team
1490  void SetCoarsening(int coarsen_type)
1491  { HYPRE_BoomerAMGSetCoarsenType(amg_precond, coarsen_type); }
1492 
1493  /// Expert option - consult hypre documentation/team
1494  void SetRelaxType(int relax_type)
1495  { HYPRE_BoomerAMGSetRelaxType(amg_precond, relax_type); }
1496 
1497  /// Expert option - consult hypre documentation/team
1498  void SetCycleType(int cycle_type)
1499  { HYPRE_BoomerAMGSetCycleType(amg_precond, cycle_type); }
1500 
1501  void GetNumIterations(int &num_iterations) const
1502  {
1503  HYPRE_Int num_it;
1504  HYPRE_BoomerAMGGetNumIterations(amg_precond, &num_it);
1505  num_iterations = internal::to_int(num_it);
1506  }
1507 
1508  /// Expert option - consult hypre documentation/team
1509  void SetNodal(int blocksize)
1510  {
1511  HYPRE_BoomerAMGSetNumFunctions(amg_precond, blocksize);
1512  HYPRE_BoomerAMGSetNodal(amg_precond, 1);
1513  }
1514 
1515  /// Expert option - consult hypre documentation/team
1516  void SetAggressiveCoarsening(int num_levels)
1517  { HYPRE_BoomerAMGSetAggNumLevels(amg_precond, num_levels); }
1518 
1519  /// The typecast to HYPRE_Solver returns the internal amg_precond
1520  virtual operator HYPRE_Solver() const { return amg_precond; }
1521 
1522  virtual HYPRE_PtrToParSolverFcn SetupFcn() const
1523  { return (HYPRE_PtrToParSolverFcn) HYPRE_BoomerAMGSetup; }
1524  virtual HYPRE_PtrToParSolverFcn SolveFcn() const
1525  { return (HYPRE_PtrToParSolverFcn) HYPRE_BoomerAMGSolve; }
1526 
1527  using HypreSolver::Mult;
1528 
1529  virtual ~HypreBoomerAMG();
1530 };
1531 
1532 /// Compute the discrete gradient matrix between the nodal linear and ND1 spaces
1533 HypreParMatrix* DiscreteGrad(ParFiniteElementSpace *edge_fespace,
1534  ParFiniteElementSpace *vert_fespace);
1535 /// Compute the discrete curl matrix between the ND1 and RT0 spaces
1536 HypreParMatrix* DiscreteCurl(ParFiniteElementSpace *face_fespace,
1537  ParFiniteElementSpace *edge_fespace);
1538 
1539 /// The Auxiliary-space Maxwell Solver in hypre
1540 class HypreAMS : public HypreSolver
1541 {
1542 private:
1543  /// Constuct AMS solver from finite element space
1544  void Init(ParFiniteElementSpace *edge_space);
1545 
1546  HYPRE_Solver ams;
1547 
1548  /// Vertex coordinates
1549  HypreParVector *x, *y, *z;
1550  /// Discrete gradient matrix
1551  HypreParMatrix *G;
1552  /// Nedelec interpolation matrix and its components
1553  HypreParMatrix *Pi, *Pix, *Piy, *Piz;
1554 
1555 public:
1556  HypreAMS(ParFiniteElementSpace *edge_fespace);
1557 
1558  HypreAMS(const HypreParMatrix &A, ParFiniteElementSpace *edge_fespace);
1559 
1560  virtual void SetOperator(const Operator &op);
1561 
1562  void SetPrintLevel(int print_lvl);
1563 
1564  /// Set this option when solving a curl-curl problem with zero mass term
1565  void SetSingularProblem() { HYPRE_AMSSetBetaPoissonMatrix(ams, NULL); }
1566 
1567  /// The typecast to HYPRE_Solver returns the internal ams object
1568  virtual operator HYPRE_Solver() const { return ams; }
1569 
1570  virtual HYPRE_PtrToParSolverFcn SetupFcn() const
1571  { return (HYPRE_PtrToParSolverFcn) HYPRE_AMSSetup; }
1572  virtual HYPRE_PtrToParSolverFcn SolveFcn() const
1573  { return (HYPRE_PtrToParSolverFcn) HYPRE_AMSSolve; }
1574 
1575  virtual ~HypreAMS();
1576 };
1577 
1578 /// The Auxiliary-space Divergence Solver in hypre
1579 class HypreADS : public HypreSolver
1580 {
1581 private:
1582  /// Constuct ADS solver from finite element space
1583  void Init(ParFiniteElementSpace *face_fespace);
1584 
1585  HYPRE_Solver ads;
1586 
1587  /// Vertex coordinates
1588  HypreParVector *x, *y, *z;
1589  /// Discrete gradient matrix
1590  HypreParMatrix *G;
1591  /// Discrete curl matrix
1592  HypreParMatrix *C;
1593  /// Nedelec interpolation matrix and its components
1594  HypreParMatrix *ND_Pi, *ND_Pix, *ND_Piy, *ND_Piz;
1595  /// Raviart-Thomas interpolation matrix and its components
1596  HypreParMatrix *RT_Pi, *RT_Pix, *RT_Piy, *RT_Piz;
1597 
1598 public:
1599  HypreADS(ParFiniteElementSpace *face_fespace);
1600 
1601  HypreADS(const HypreParMatrix &A, ParFiniteElementSpace *face_fespace);
1602 
1603  virtual void SetOperator(const Operator &op);
1604 
1605  void SetPrintLevel(int print_lvl);
1606 
1607  /// The typecast to HYPRE_Solver returns the internal ads object
1608  virtual operator HYPRE_Solver() const { return ads; }
1609 
1610  virtual HYPRE_PtrToParSolverFcn SetupFcn() const
1611  { return (HYPRE_PtrToParSolverFcn) HYPRE_ADSSetup; }
1612  virtual HYPRE_PtrToParSolverFcn SolveFcn() const
1613  { return (HYPRE_PtrToParSolverFcn) HYPRE_ADSSolve; }
1614 
1615  virtual ~HypreADS();
1616 };
1617 
1618 /** LOBPCG eigenvalue solver in hypre
1619 
1620  The Locally Optimal Block Preconditioned Conjugate Gradient (LOBPCG)
1621  eigenvalue solver is designed to find the lowest eigenmodes of the
1622  generalized eigenvalue problem:
1623  A x = lambda M x
1624  where A is symmetric, potentially indefinite and M is symmetric positive
1625  definite. The eigenvectors are M-orthonormal, meaning that
1626  x^T M x = 1 and x^T M y = 0,
1627  if x and y are distinct eigenvectors. The matrix M is optional and is
1628  assumed to be the identity if left unset.
1629 
1630  The efficiency of LOBPCG relies on the availability of a suitable
1631  preconditioner for the matrix A. The preconditioner is supplied through the
1632  SetPreconditioner() method. It should be noted that the operator used with
1633  the preconditioner need not be A itself.
1634 
1635  For more information regarding LOBPCG see "Block Locally Optimal
1636  Preconditioned Eigenvalue Xolvers (BLOPEX) in Hypre and PETSc" by
1637  A. Knyazev, M. Argentati, I. Lashuk, and E. Ovtchinnikov, SISC, 29(5),
1638  2224-2239, 2007.
1639 */
1641 {
1642 private:
1643  MPI_Comm comm;
1644  int myid;
1645  int numProcs;
1646  int nev; // Number of desired eigenmodes
1647  int seed; // Random seed used for initial vectors
1648 
1649  HYPRE_BigInt glbSize; // Global number of DoFs in the linear system
1650  HYPRE_BigInt * part; // Row partitioning of the linear system
1651 
1652  // Pointer to HYPRE's solver struct
1653  HYPRE_Solver lobpcg_solver;
1654 
1655  // Interface for matrix storage type
1656  mv_InterfaceInterpreter interpreter;
1657 
1658  // Interface for setting up and performing matrix-vector products
1659  HYPRE_MatvecFunctions matvec_fn;
1660 
1661  // Eigenvalues
1662  Array<double> eigenvalues;
1663 
1664  // Forward declaration
1665  class HypreMultiVector;
1666 
1667  // MultiVector to store eigenvectors
1668  HypreMultiVector * multi_vec;
1669 
1670  // Empty vectors used to setup the matrices and preconditioner
1671  HypreParVector * x;
1672 
1673  // An optional operator which projects vectors into a desired subspace
1674  Operator * subSpaceProj;
1675 
1676  /// Internal class to represent a set of eigenvectors
1677  class HypreMultiVector
1678  {
1679  private:
1680  // Pointer to hypre's multi-vector object
1681  mv_MultiVectorPtr mv_ptr;
1682 
1683  // Wrappers for each member of the multivector
1684  HypreParVector ** hpv;
1685 
1686  // Number of vectors in the multivector
1687  int nv;
1688 
1689  public:
1690  HypreMultiVector(int n, HypreParVector & v,
1691  mv_InterfaceInterpreter & interpreter);
1692  ~HypreMultiVector();
1693 
1694  /// Set random values
1695  void Randomize(HYPRE_Int seed);
1696 
1697  /// Extract a single HypreParVector object
1698  HypreParVector & GetVector(unsigned int i);
1699 
1700  /// Transfers ownership of data to returned array of vectors
1701  HypreParVector ** StealVectors();
1702 
1703  operator mv_MultiVectorPtr() const { return mv_ptr; }
1704 
1705  mv_MultiVectorPtr & GetMultiVector() { return mv_ptr; }
1706  };
1707 
1708  static void * OperatorMatvecCreate( void *A, void *x );
1709  static HYPRE_Int OperatorMatvec( void *matvec_data,
1710  HYPRE_Complex alpha,
1711  void *A,
1712  void *x,
1713  HYPRE_Complex beta,
1714  void *y );
1715  static HYPRE_Int OperatorMatvecDestroy( void *matvec_data );
1716 
1717  static HYPRE_Int PrecondSolve(void *solver,
1718  void *A,
1719  void *b,
1720  void *x);
1721  static HYPRE_Int PrecondSetup(void *solver,
1722  void *A,
1723  void *b,
1724  void *x);
1725 
1726 public:
1727  HypreLOBPCG(MPI_Comm comm);
1728  ~HypreLOBPCG();
1729 
1730  void SetTol(double tol);
1731  void SetRelTol(double rel_tol);
1732  void SetMaxIter(int max_iter);
1733  void SetPrintLevel(int logging);
1734  void SetNumModes(int num_eigs) { nev = num_eigs; }
1735  void SetPrecondUsageMode(int pcg_mode);
1736  void SetRandomSeed(int s) { seed = s; }
1737  void SetInitialVectors(int num_vecs, HypreParVector ** vecs);
1738 
1739  // The following four methods support general operators
1740  void SetPreconditioner(Solver & precond);
1741  void SetOperator(Operator & A);
1742  void SetMassMatrix(Operator & M);
1743  void SetSubSpaceProjector(Operator & proj) { subSpaceProj = &proj; }
1744 
1745  /// Solve the eigenproblem
1746  void Solve();
1747 
1748  /// Collect the converged eigenvalues
1749  void GetEigenvalues(Array<double> & eigenvalues) const;
1750 
1751  /// Extract a single eigenvector
1752  const HypreParVector & GetEigenvector(unsigned int i) const;
1753 
1754  /// Transfer ownership of the converged eigenvectors
1755  HypreParVector ** StealEigenvectors() { return multi_vec->StealVectors(); }
1756 };
1757 
1758 /** AME eigenvalue solver in hypre
1759 
1760  The Auxiliary space Maxwell Eigensolver (AME) is designed to find
1761  the lowest eigenmodes of the generalized eigenvalue problem:
1762  Curl Curl x = lambda M x
1763  where the Curl Curl operator is discretized using Nedelec finite element
1764  basis functions. Properties of this discretization are essential to
1765  eliminating the large null space of the Curl Curl operator.
1766 
1767  This eigensolver relies upon the LOBPCG eigensolver internally. It is also
1768  expected that the preconditioner supplied to this method will be the
1769  HypreAMS preconditioner defined above.
1770 
1771  As with LOBPCG, the operator set in the preconditioner need not be the same
1772  as A. This flexibility may be useful in solving eigenproblems which bare a
1773  strong resemblance to the Curl Curl problems for which AME is designed.
1774 
1775  Unlike LOBPCG, this eigensolver requires that the mass matrix be set.
1776  It is possible to circumvent this by passing an identity operator as the
1777  mass matrix but it seems unlikely that this would be useful so it is not the
1778  default behavior.
1779 */
1781 {
1782 private:
1783  int myid;
1784  int numProcs;
1785  int nev; // Number of desired eigenmodes
1786  bool setT;
1787 
1788  // Pointer to HYPRE's AME solver struct
1789  HYPRE_Solver ame_solver;
1790 
1791  // Pointer to HYPRE's AMS solver struct
1792  HypreSolver * ams_precond;
1793 
1794  // Eigenvalues
1795  HYPRE_Real * eigenvalues;
1796 
1797  // MultiVector to store eigenvectors
1798  HYPRE_ParVector * multi_vec;
1799 
1800  // HypreParVector wrappers to contain eigenvectors
1801  mutable HypreParVector ** eigenvectors;
1802 
1803  void createDummyVectors() const;
1804 
1805 public:
1806  HypreAME(MPI_Comm comm);
1807  ~HypreAME();
1808 
1809  void SetTol(double tol);
1810  void SetRelTol(double rel_tol);
1811  void SetMaxIter(int max_iter);
1812  void SetPrintLevel(int logging);
1813  void SetNumModes(int num_eigs);
1814 
1815  // The following four methods support operators of type HypreParMatrix.
1816  void SetPreconditioner(HypreSolver & precond);
1817  void SetOperator(const HypreParMatrix & A);
1818  void SetMassMatrix(const HypreParMatrix & M);
1819 
1820  /// Solve the eigenproblem
1821  void Solve();
1822 
1823  /// Collect the converged eigenvalues
1824  void GetEigenvalues(Array<double> & eigenvalues) const;
1825 
1826  /// Extract a single eigenvector
1827  const HypreParVector & GetEigenvector(unsigned int i) const;
1828 
1829  /// Transfer ownership of the converged eigenvectors
1831 };
1832 
1833 }
1834 
1835 #endif // MFEM_USE_MPI
1836 
1837 #endif
Type GetType() const
Definition: hypre.hpp:779
MemoryType GetHypreMemoryType()
The MemoryType used by MFEM when allocating arrays for Hypre objects.
Definition: hypre.hpp:87
virtual ~HypreBoomerAMG()
Definition: hypre.cpp:4725
void SetAbsTol(double atol)
Definition: hypre.cpp:3574
Memory< double > auxX
Definition: hypre.hpp:850
void SetTol(double tol)
Definition: hypre.cpp:3569
void SetPreconditioner(HypreSolver &precond)
Set the hypre solver to be used as a preconditioner.
Definition: hypre.cpp:3784
void EliminateRowsCols(const Array< int > &rows_cols, const HypreParVector &X, HypreParVector &B)
Definition: hypre.cpp:2052
HypreADS(ParFiniteElementSpace *face_fespace)
Definition: hypre.cpp:4986
const HypreParMatrix * GetData() const
Definition: hypre.hpp:1049
static constexpr Type default_type
Definition: hypre.hpp:912
MPI_Comm GetComm() const
MPI communicator.
Definition: hypre.hpp:472
HypreEuclid(MPI_Comm comm)
Definition: hypre.cpp:4155
void SetErrorMode(ErrorMode err_mode) const
Set the behavior for treating hypre errors, see the ErrorMode enum. The default mode in the base clas...
Definition: hypre.hpp:1027
virtual void SetOperator(const Operator &op)
Set/update the solver for the given operator.
Definition: hypre.cpp:3904
Vector * GlobalVector() const
Returns the global vector in each processor.
Definition: hypre.cpp:182
void EliminateBC(const HypreParMatrix &A, const HypreParMatrix &Ae, const Array< int > &ess_dof_list, const Vector &X, Vector &B)
Eliminate essential BC specified by ess_dof_list from the solution X to the r.h.s. B.
Definition: hypre.cpp:2857
The Auxiliary-space Maxwell Solver in hypre.
Definition: hypre.hpp:1540
HypreParVector * X0
FIR Filter Temporary Vectors.
Definition: hypre.hpp:854
void SetStrengthThresh(double strength)
Expert option - consult hypre documentation/team.
Definition: hypre.hpp:1482
double min_eig_est
Minimal eigenvalue estimate for polynomial smoothing.
Definition: hypre.hpp:886
virtual HYPRE_PtrToParSolverFcn SolveFcn() const
hypre&#39;s internal Solve function
Definition: hypre.hpp:1227
static MemoryClass GetHostMemoryClass()
Get the current Host MemoryClass. This is the MemoryClass used by most MFEM host Memory objects...
Definition: device.hpp:268
ErrorMode
How to treat errors returned by hypre function calls.
Definition: hypre.hpp:974
The Auxiliary-space Divergence Solver in hypre.
Definition: hypre.hpp:1579
virtual HYPRE_PtrToParSolverFcn SolveFcn() const
FGMRES Solve function.
Definition: hypre.hpp:1209
void SetPrintLevel(int print_lvl)
Definition: hypre.cpp:5273
MFEM_DEPRECATED HypreParMatrix * GetData()
Deprecated. Use HypreTriSolve::GetData() const instead.
Definition: hypre.hpp:1052
void SetMassMatrix(const HypreParMatrix &M)
Definition: hypre.cpp:5747
HypreParMatrix * RAP(const HypreParMatrix *A, const HypreParMatrix *P)
Returns the matrix P^t * A * P.
Definition: hypre.cpp:2485
virtual ~HypreEuclid()
Definition: hypre.cpp:4220
void WrapMemoryWrite(Memory< double > &mem)
Replace the HypreParVector&#39;s data with the given Memory, mem, and prepare the vector for write access...
Definition: hypre.cpp:274
int setup_called
Was hypre&#39;s Setup function called already?
Definition: hypre.hpp:991
Device memory; using CUDA or HIP *Malloc and *Free.
void HypreRead() const
Prepare the HypreParVector for read access in hypre&#39;s device memory space, HYPRE_MEMORY_DEVICE.
Definition: hypre.cpp:217
void EliminateBC(const HypreParMatrix &Ae, const Array< int > &ess_dof_list, const Vector &X, Vector &B) const
Eliminate essential BC specified by ess_dof_list from the solution X to the r.h.s. B.
Definition: hypre.cpp:2104
void HypreRead() const
Update the internal hypre_ParCSRMatrix object, A, to be in hypre memory space.
Definition: hypre.hpp:743
HYPRE_BigInt GlobalSize() const
Returns the global number of rows.
Definition: hypre.hpp:171
virtual void SetOperator(const Operator &op)
Set/update the solver for the given operator.
Definition: hypre.cpp:3547
void SetPreconditioner(HypreSolver &precond)
Set the hypre solver to be used as a preconditioner.
Definition: hypre.cpp:3951
void MakeRef(const HypreParMatrix &master)
Make this HypreParMatrix a reference to &#39;master&#39;.
Definition: hypre.cpp:1213
void Read_IJMatrix(MPI_Comm comm, const char *fname)
Read a matrix saved as a HYPRE_IJMatrix.
Definition: hypre.cpp:2185
HypreParMatrix & Add(const double beta, const HypreParMatrix &B)
Definition: hypre.hpp:660
HypreParVector * B
Right-hand side and solution vector.
Definition: hypre.hpp:986
Wrapper for Hypre&#39;s native parallel ILU preconditioner.
Definition: hypre.hpp:1345
void GetNumIterations(int &num_iterations) const
Definition: hypre.hpp:1094
const HypreParMatrix * GetData() const
Definition: hypre.hpp:1248
double window_params[3]
Parameters for windowing function of FIR filter.
Definition: hypre.hpp:888
virtual void Mult(const Vector &x, Vector &y) const
Operator application: y=A(x).
Definition: hypre.hpp:605
virtual HYPRE_PtrToParSolverFcn SolveFcn() const
hypre&#39;s internal Solve function
Definition: hypre.hpp:1572
void SetInitialVectors(int num_vecs, HypreParVector **vecs)
Definition: hypre.cpp:5499
virtual void SetOperator(const Operator &op)
Set/update the solver for the given operator.
Definition: hypre.cpp:3732
void SetWindowByName(const char *window_name)
Convenience function for setting canonical windowing parameters.
Definition: hypre.cpp:3056
HypreTriSolve(const HypreParMatrix &A)
Definition: hypre.hpp:1041
Issue warnings on hypre errors.
Definition: hypre.hpp:977
void SetPreconditioner(HypreSolver &precond)
Definition: hypre.cpp:5726
HyprePCG(MPI_Comm comm)
Definition: hypre.cpp:3529
void SetTol(double tol)
Definition: hypre.cpp:5695
HypreILU()
Constructor; sets the default options.
Definition: hypre.cpp:4227
void SetRestriction(int restrict_type)
Expert option - consult hypre documentation/team.
Definition: hypre.hpp:1448
HypreAMS(ParFiniteElementSpace *edge_fespace)
Definition: hypre.cpp:4735
void SetKDim(int dim)
Definition: hypre.cpp:3936
virtual ~HypreGMRES()
Definition: hypre.cpp:3866
void SetCycleNumSweeps(int prerelax, int postrelax)
Expert option - consult hypre documentation/team.
Definition: hypre.hpp:1460
void SetPrintLevel(int print_lvl)
Definition: hypre.cpp:3946
virtual void SetOperator(const Operator &op)
Set/update the solver for the given operator.
Definition: hypre.cpp:5233
void SetPreconditioner(Solver &precond)
Definition: hypre.cpp:5415
MFEM_DEPRECATED void SetZeroInintialIterate()
deprecated: use SetZeroInitialIterate()
Definition: hypre.hpp:1147
void SetMassMatrix(Operator &M)
Definition: hypre.cpp:5471
Abstract parallel finite element space.
Definition: pfespace.hpp:28
HypreParVector()
Default constructor, no underlying hypre_ParVector is created.
Definition: hypre.hpp:115
void SetPrintLevel(int logging)
Definition: hypre.cpp:5400
bool iterative_mode
If true, use the second argument of Mult() as an initial guess.
Definition: operator.hpp:652
virtual ~HypreAMS()
Definition: hypre.cpp:4966
void SetTol(double tol)
Definition: hypre.cpp:3754
virtual ~HypreILU()
Definition: hypre.cpp:4299
void SetAdvectiveOptions(int distance=15, const std::string &prerelax="", const std::string &postrelax="FFC")
Definition: hypre.cpp:4617
virtual HYPRE_PtrToParSolverFcn SolveFcn() const
GMRES Solve function.
Definition: hypre.hpp:1159
void BlockInverseScale(const HypreParMatrix *A, HypreParMatrix *C, const Vector *b, HypreParVector *d, int blocksize, BlockInverseScaleJob job)
Definition: hypre.cpp:2378
HypreParVector ** StealEigenvectors()
Transfer ownership of the converged eigenvectors.
Definition: hypre.hpp:1755
void HostWrite()
Update the internal hypre_ParCSRMatrix object, A, to be on host.
Definition: hypre.hpp:737
HypreParMatrix & operator+=(const HypreParMatrix &B)
Definition: hypre.hpp:653
bool pos_l1_norms
If set, take absolute values of the computed l1_norms.
Definition: hypre.hpp:880
HypreParMatrix * ParMult(const HypreParMatrix *A, const HypreParMatrix *B, bool own_matrix)
Definition: hypre.cpp:2464
void SetPrintLevel(int print_lvl)
Definition: hypre.cpp:3779
void SetZeroInitialIterate()
non-hypre setting
Definition: hypre.hpp:1200
void SetStrongThresholdR(double strengthR)
Expert option - consult hypre documentation/team.
Definition: hypre.hpp:1440
MPI_Comm GetComm() const
MPI communicator.
Definition: hypre.hpp:154
void ScaleRows(const Vector &s)
Scale the local row i by s(i).
Definition: hypre.cpp:1802
virtual MFEM_DEPRECATED N_Vector ToNVector()
(DEPRECATED) Return a new wrapper SUNDIALS N_Vector of type SUNDIALS_NVEC_PARALLEL.
Definition: hypre.cpp:308
double poly_fraction
Fraction of spectrum to smooth for polynomial relaxation.
Definition: hypre.hpp:868
void SetPrintLevel(int print_lvl)
Definition: hypre.cpp:3589
virtual void SetOperator(const Operator &op)
Set/update the solver for the given operator.
Definition: hypre.cpp:4280
HypreParMatrix * DiscreteGrad(ParFiniteElementSpace *edge_fespace, ParFiniteElementSpace *vert_fespace)
Compute the discrete gradient matrix between the nodal linear and ND1 spaces.
int poly_scale
Apply the polynomial smoother to A or D^{-1/2} A D^{-1/2}.
Definition: hypre.hpp:870
void MergeDiagAndOffd(SparseMatrix &merged)
Get a single SparseMatrix containing all rows from this processor, merged from the diagonal and off-d...
Definition: hypre.cpp:1423
void GetEigenvalues(Array< double > &eigenvalues) const
Collect the converged eigenvalues.
Definition: hypre.cpp:5481
void BooleanMultTranspose(int alpha, const int *x, int beta, int *y)
The &quot;Boolean&quot; analog of y = alpha * A^T * x + beta * y, where elements in the sparsity pattern of the...
Definition: hypre.hpp:630
void SetTol(double tol)
Definition: hypre.cpp:5378
HypreParVector * X1
Definition: hypre.hpp:854
HypreFGMRES(MPI_Comm comm)
Definition: hypre.cpp:3872
HYPRE_BigInt N() const
Returns the global number of columns.
Definition: hypre.hpp:523
virtual void Mult(const HypreParVector &b, HypreParVector &x) const
Solve Ax=b with hypre&#39;s FGMRES.
Definition: hypre.cpp:3960
HYPRE_BigInt M() const
Returns the global number of rows.
Definition: hypre.hpp:521
The identity operator as a hypre solver.
Definition: hypre.hpp:1220
void SetWindowParameters(double a, double b, double c)
Set parameters for windowing function for FIR smoother.
Definition: hypre.cpp:3071
void EliminateZeroRows()
If a row contains only zeros, set its diagonal to 1.
Definition: hypre.hpp:696
virtual ~HypreTriSolve()
Definition: hypre.hpp:1055
HypreGMRES(MPI_Comm comm)
Definition: hypre.cpp:3700
void GetOffd(SparseMatrix &offd, HYPRE_BigInt *&cmap) const
Get the local off-diagonal block. NOTE: &#39;offd&#39; will not own any data.
Definition: hypre.cpp:1417
virtual void SetOperator(const Operator &op)
Set/update the solver for the given operator.
Definition: hypre.cpp:4120
virtual void Mult(const HypreParVector &b, HypreParVector &x) const
Solve Ax=b with hypre&#39;s GMRES.
Definition: hypre.cpp:3794
virtual void AssembleDiagonal(Vector &diag) const
Return the diagonal of the matrix (Operator interface).
Definition: hypre.hpp:540
void SetPositiveDiagonal(bool pos=true)
After computing l1-norms, replace them with their absolute values.
Definition: hypre.hpp:946
void Print(const char *fname) const
Prints the locally owned rows in parallel.
Definition: hypre.cpp:293
void HypreWrite()
Update the internal hypre_ParCSRMatrix object, A, to be in hypre memory space.
Definition: hypre.hpp:757
void AbsMultTranspose(double a, const Vector &x, double b, Vector &y) const
Computes y = a * |At| * x + b * y, using entry-wise absolute values of the transpose of the matrix A...
Definition: hypre.cpp:1688
HypreLOBPCG(MPI_Comm comm)
Definition: hypre.cpp:5348
HYPRE_BigInt * GetRowStarts() const
Return the parallel row partitioning array.
Definition: hypre.hpp:583
void Add(const DenseMatrix &A, const DenseMatrix &B, double alpha, DenseMatrix &C)
C = A + alpha*B.
Definition: densemat.cpp:1931
MFEM_DEPRECATED HYPRE_BigInt * Partitioning()
Returns a non-const pointer to the parallel row/column partitioning. Deprecated in favor of HypreParV...
Definition: hypre.hpp:168
virtual void Mult(const HypreParVector &b, HypreParVector &x) const
Relax the linear system Ax=b.
Definition: hypre.cpp:3222
hypre_ParVector * StealParVector()
Changes the ownership of the vector.
Definition: hypre.hpp:180
The BoomerAMG solver in hypre.
Definition: hypre.hpp:1387
virtual void SetOperator(const Operator &op)
Set/update the solver for the given operator.
Definition: hypre.cpp:4038
void SetRelaxType(int relax_type)
Expert option - consult hypre documentation/team.
Definition: hypre.hpp:1494
void SetSymmetry(int sym)
Definition: hypre.cpp:4144
virtual MemoryClass GetMemoryClass() const
Return the MemoryClass preferred by the Operator.
Definition: hypre.hpp:590
virtual void MultTranspose(const Vector &x, Vector &y) const
Action of the transpose operator: y=A^t(x). The default behavior in class Operator is to generate an ...
Definition: hypre.hpp:607
virtual ~HypreParaSails()
Definition: hypre.cpp:4149
The ParaSails preconditioner in hypre.
Definition: hypre.hpp:1258
void SetMaxIter(int max_iter)
Definition: hypre.cpp:5394
double * l1_norms
l1 norms of the rows of A
Definition: hypre.hpp:878
Data type sparse matrix.
Definition: sparsemat.hpp:41
BlockInverseScaleJob
Definition: hypre.hpp:784
void SetOwnerFlags(signed char diag, signed char offd, signed char colmap)
Explicitly set the three ownership flags, see docs for diagOwner etc.
Definition: hypre.cpp:1250
void PrintHash(std::ostream &out) const
Print sizes and hashes for all data arrays of the HypreParMatrix from the local MPI rank...
Definition: hypre.cpp:2243
void SetLogging(int logging)
Definition: hypre.cpp:3584
Jacobi preconditioner in hypre.
Definition: hypre.hpp:1234
virtual ~HypreSolver()
Definition: hypre.cpp:3520
virtual HYPRE_PtrToParSolverFcn SolveFcn() const =0
hypre&#39;s internal Solve function
void SetRelTol(double rel_tol)
Definition: hypre.cpp:5701
void SetKDim(int dim)
Definition: hypre.cpp:3769
void mfem_error(const char *msg)
Function called when an error is encountered. Used by the macros MFEM_ABORT, MFEM_ASSERT, MFEM_VERIFY.
Definition: error.cpp:153
void SetResidualConvergenceOptions(int res_frequency=-1, double rtol=0.0)
Definition: hypre.cpp:3604
double b
Definition: lissajous.cpp:42
int to_int(const std::string &str)
Convert a string to an int.
Definition: text.hpp:62
Abort on hypre errors (default in base class)
Definition: hypre.hpp:978
void SetLogging(int logging)
Definition: hypre.cpp:3941
void SetPrintLevel(int print_level)
Definition: hypre.hpp:1467
virtual HYPRE_PtrToParSolverFcn SetupFcn() const
FGMRES Setup function.
Definition: hypre.hpp:1206
void SetOperatorSymmetry(bool is_sym)
Definition: hypre.hpp:952
void GetBlocks(Array2D< HypreParMatrix * > &blocks, bool interleaved_rows=false, bool interleaved_cols=false) const
Definition: hypre.cpp:1447
virtual void Mult(const HypreParVector &b, HypreParVector &x) const
Solve the linear system Ax=b.
Definition: hypre.cpp:3414
HypreParaSails(MPI_Comm comm)
Definition: hypre.cpp:4056
void HypreWrite()
Prepare the HypreParVector for write access in hypre&#39;s device memory space, HYPRE_MEMORY_DEVICE.
Definition: hypre.cpp:236
void SetTol(double tol)
Definition: hypre.cpp:3926
void WrapHypreParVector(hypre_ParVector *y, bool owner=true)
Converts hypre&#39;s format to HypreParVector.
Definition: hypre.cpp:173
void SetTol(double tol)
Expert option - consult hypre documentation/team.
Definition: hypre.hpp:1478
virtual HYPRE_PtrToParSolverFcn SolveFcn() const
PCG Solve function.
Definition: hypre.hpp:1108
virtual void SetOperator(const Operator &op)
Set/update the solver for the given operator.
Definition: hypre.cpp:4948
virtual MemoryClass GetMemoryClass() const
Return the MemoryClass preferred by the Operator.
Definition: hypre.hpp:1012
HypreParVector ** StealEigenvectors()
Transfer ownership of the converged eigenvectors.
Definition: hypre.cpp:5801
const HypreParMatrix * A
The linear system matrix.
Definition: hypre.hpp:983
HypreParMatrix * DiscreteCurl(ParFiniteElementSpace *face_fespace, ParFiniteElementSpace *edge_fespace)
Compute the discrete curl matrix between the ND1 and RT0 spaces.
HypreParMatrix * Transpose() const
Returns the transpose of *this.
Definition: hypre.cpp:1471
HYPRE_Int Randomize(HYPRE_Int seed)
Set random values.
Definition: hypre.cpp:288
void SetFilterThresholdR(double filterR)
Expert option - consult hypre documentation/team.
Definition: hypre.hpp:1444
Ignore hypre errors (see e.g. HypreADS)
Definition: hypre.hpp:976
double relax_weight
Damping coefficient (usually &lt;= 1)
Definition: hypre.hpp:862
Parallel smoothers in hypre.
Definition: hypre.hpp:840
void SetElasticityOptions(ParFiniteElementSpace *fespace)
Definition: hypre.cpp:4569
Memory< double > auxB
Definition: hypre.hpp:988
void GetDiag(Vector &diag) const
Get the local diagonal of the matrix.
Definition: hypre.cpp:1344
virtual HYPRE_PtrToParSolverFcn SetupFcn() const
hypre&#39;s internal Setup function
Definition: hypre.hpp:1522
virtual HYPRE_PtrToParSolverFcn SolveFcn() const
ILU Solve function.
Definition: hypre.hpp:1381
virtual HYPRE_PtrToParSolverFcn SetupFcn() const
hypre&#39;s internal Setup function
Definition: hypre.hpp:1225
void SetZeroInitialIterate()
non-hypre setting
Definition: hypre.hpp:1150
void SetLevelOfFill(HYPRE_Int lev_fill)
Set the fill level for ILU(k); the default is k=1.
Definition: hypre.cpp:4270
HypreParMatrix * A
The linear system matrix.
Definition: hypre.hpp:844
void SetCoarsening(int coarsen_type)
Expert option - consult hypre documentation/team.
Definition: hypre.hpp:1490
void SetGMRESSwitchR(int gmres_switch)
Expert option - consult hypre documentation/team.
Definition: hypre.hpp:1456
virtual HYPRE_PtrToParSolverFcn SolveFcn() const
hypre&#39;s internal Solve function
Definition: hypre.hpp:1245
Dynamic 2D array using row-major layout.
Definition: array.hpp:347
virtual void SetOperator(const Operator &op)
Definition: hypre.cpp:3078
virtual HYPRE_PtrToParSolverFcn SetupFcn() const =0
hypre&#39;s internal Setup function
void SetLogging(int logging)
Definition: hypre.cpp:3774
virtual ~HypreADS()
Definition: hypre.cpp:5251
HYPRE_BigInt GetGlobalNumCols() const
Return the global number of columns.
Definition: hypre.hpp:577
HypreParMatrix()
An empty matrix to be used as a reference to an existing matrix.
Definition: hypre.cpp:504
virtual HYPRE_PtrToParSolverFcn SetupFcn() const
hypre&#39;s internal Setup function
Definition: hypre.hpp:1243
bool A_is_symmetric
A flag that indicates whether the linear system matrix A is symmetric.
Definition: hypre.hpp:894
void SetMaxIter(int max_iter)
Definition: hypre.cpp:3579
void SetMaxIter(int max_iter)
Definition: hypre.hpp:1470
void SetMaxIter(int max_iter)
Definition: hypre.cpp:3931
void SetPolyOptions(int poly_order, double poly_fraction, int eig_est_cg_iter=10)
Set parameters for polynomial smoothing.
Definition: hypre.cpp:3040
HypreDiagScale(const HypreParMatrix &A)
Definition: hypre.hpp:1238
void SetSOROptions(double relax_weight, double omega)
Set SOR-related parameters.
Definition: hypre.cpp:3034
void BooleanMult(int alpha, const int *x, int beta, int *y)
The &quot;Boolean&quot; analog of y = alpha * A * x + beta * y, where elements in the sparsity pattern of the m...
Definition: hypre.hpp:620
void SetInterpolation(int interp_type)
Expert option - consult hypre documentation/team.
Definition: hypre.hpp:1486
Wrapper for hypre&#39;s parallel vector class.
Definition: hypre.hpp:99
HypreParMatrix * EliminateCols(const Array< int > &cols)
Definition: hypre.cpp:2077
HypreParVector * X
Definition: hypre.hpp:986
virtual HYPRE_PtrToParSolverFcn SetupFcn() const
hypre&#39;s internal Setup function
Definition: hypre.hpp:1610
Type
Enumeration defining IDs for some classes derived from Operator.
Definition: operator.hpp:252
static MemoryType GetHostMemoryType()
Get the current Host MemoryType. This is the MemoryType used by most MFEM classes when allocating mem...
Definition: device.hpp:264
PCG solver in hypre.
Definition: hypre.hpp:1060
GMRES solver in hypre.
Definition: hypre.hpp:1119
signed char OwnsOffd() const
Get offd ownership flag.
Definition: hypre.hpp:489
void HostRead() const
Update the internal hypre_ParCSRMatrix object, A, to be on host.
Definition: hypre.hpp:726
MemoryType
Memory types supported by MFEM.
Definition: mem_manager.hpp:31
signed char OwnsColMap() const
Get colmap ownership flag.
Definition: hypre.hpp:491
void SetPrintLevel(int print_lvl)
Definition: hypre.cpp:4981
Memory< double > auxX
Definition: hypre.hpp:988
HypreParVector * X
Definition: hypre.hpp:846
void Solve()
Solve the eigenproblem.
Definition: hypre.cpp:5754
HYPRE_Int HYPRE_BigInt
void SetData(double *data_)
Sets the data of the Vector and the hypre_ParVector to data_.
Definition: hypre.cpp:211
void SetNumModes(int num_eigs)
Definition: hypre.cpp:5687
virtual HYPRE_PtrToParSolverFcn SolveFcn() const
hypre&#39;s internal Solve function
Definition: hypre.hpp:1285
HypreParMatrix(hypre_ParCSRMatrix *a, bool owner=true)
Converts hypre&#39;s format to HypreParMatrix.
Definition: hypre.hpp:364
virtual ~HypreDiagScale()
Definition: hypre.hpp:1254
const HypreParVector & GetEigenvector(unsigned int i) const
Extract a single eigenvector.
Definition: hypre.cpp:5790
int GetNumRows() const
Returns the number of rows in the diagonal block of the ParCSRMatrix.
Definition: hypre.hpp:559
virtual HYPRE_PtrToParSolverFcn SolveFcn() const
hypre&#39;s internal Solve function
Definition: hypre.hpp:1046
const HYPRE_BigInt * ColPart() const
Returns the column partitioning (const version)
Definition: hypre.hpp:519
virtual void SetOperator(const Operator &op)
Set/update the solver for the given operator.
Definition: hypre.hpp:1009
HypreParMatrix * ParAdd(const HypreParMatrix *A, const HypreParMatrix *B)
Returns the matrix A + B.
Definition: hypre.cpp:2425
double max_eig_est
Maximal eigenvalue estimate for polynomial smoothing.
Definition: hypre.hpp:884
A class to initialize the size of a Tensor.
Definition: dtensor.hpp:54
void SetRelTol(double rel_tol)
Definition: hypre.cpp:5384
HypreAME(MPI_Comm comm)
Definition: hypre.cpp:5645
HYPRE_Int MultTranspose(HypreParVector &x, HypreParVector &y, double alpha=1.0, double beta=0.0) const
Computes y = alpha * A^t * x + beta * y.
Definition: hypre.cpp:1663
void GetNumIterations(int &num_iterations) const
Definition: hypre.hpp:1501
virtual HYPRE_PtrToParSolverFcn SetupFcn() const
hypre&#39;s internal Setup function
Definition: hypre.hpp:1283
virtual ~HypreIdentity()
Definition: hypre.hpp:1230
void SetMaxLevels(int max_levels)
Expert option - consult hypre documentation/team.
Definition: hypre.hpp:1474
void SetMaxIter(int max_iter)
Definition: hypre.cpp:5711
void SetFIRCoefficients(double max_eig)
Compute window and Chebyshev coefficients for given polynomial order.
Definition: hypre.cpp:3182
void SetRandomSeed(int s)
Definition: hypre.hpp:1736
double a
Definition: lissajous.cpp:41
void EliminateRows(const Array< int > &rows)
Eliminate rows from the diagonal and off-diagonal blocks of the matrix.
Definition: hypre.cpp:2091
signed char OwnsDiag() const
Get diag ownership flag.
Definition: hypre.hpp:487
double * fir_coeffs
Combined coefficients for windowing and Chebyshev polynomials.
Definition: hypre.hpp:891
void SetIsTriangular()
Expert option - consult hypre documentation/team.
Definition: hypre.hpp:1452
void Threshold(double threshold=0.0)
Remove values smaller in absolute value than some threshold.
Definition: hypre.cpp:1930
double InnerProduct(HypreParVector *x, HypreParVector *y)
Definition: hypre.cpp:316
void Print(const char *fname, HYPRE_Int offi=0, HYPRE_Int offj=0) const
Prints the locally owned rows in parallel.
Definition: hypre.cpp:2162
virtual ~HypreSmoother()
Definition: hypre.cpp:3373
void SetPrintLevel(HYPRE_Int print_level)
Set the print level: 0 = none, 1 = setup, 2 = solve, 3 = setup+solve.
Definition: hypre.cpp:4275
Host memory; using new[] and delete[].
void SetSubSpaceProjector(Operator &proj)
Definition: hypre.hpp:1743
virtual MFEM_DEPRECATED N_Vector ToNVector()
(DEPRECATED) Return a new wrapper SUNDIALS N_Vector of type SUNDIALS_NVEC_SERIAL. ...
Definition: vector.hpp:455
constexpr MemoryClass GetHypreMemoryClass()
The MemoryClass used by Hypre objects.
Definition: hypre.hpp:75
HYPRE_BigInt GetGlobalNumRows() const
Return the global number of rows.
Definition: hypre.hpp:573
HYPRE_BigInt * ColPart()
Returns the column partitioning.
Definition: hypre.hpp:511
void HypreReadWrite()
Prepare the HypreParVector for read and write access in hypre&#39;s device memory space, HYPRE_MEMORY_DEVICE.
Definition: hypre.cpp:227
virtual void SetOperator(const Operator &op)
Set/update the solver for the given operator.
Definition: hypre.cpp:4447
void InvScaleRows(const Vector &s)
Scale the local row i by 1./s(i)
Definition: hypre.cpp:1841
MFEM_DEPRECATED void SetZeroInintialIterate()
deprecated: use SetZeroInitialIterate()
Definition: hypre.hpp:1089
void SetTaubinOptions(double lambda, double mu, int iter)
Set parameters for Taubin&#39;s lambda-mu method.
Definition: hypre.cpp:3048
HypreParVector * B
Right-hand side and solution vectors.
Definition: hypre.hpp:846
HYPRE_BigInt * RowPart()
Returns the row partitioning.
Definition: hypre.hpp:507
virtual void SetOperator(const Operator &op)
Set/update the solver for the given operator.
Definition: hypre.cpp:4196
virtual void MultTranspose(const Vector &b, Vector &x) const
Apply transpose of the smoother to relax the linear system Ax=b.
Definition: hypre.cpp:3363
int dim
Definition: ex24.cpp:53
virtual ~HyprePCG()
Definition: hypre.cpp:3694
int GetNumCols() const
Returns the number of columns in the diagonal block of the ParCSRMatrix.
Definition: hypre.hpp:566
void SetZeroInitialIterate()
non-hypre setting
Definition: hypre.hpp:1092
void SetOperator(Operator &A)
Definition: hypre.cpp:5424
void SetPreconditioner(HypreSolver &precond)
Set the hypre solver to be used as a preconditioner.
Definition: hypre.cpp:3594
virtual HYPRE_PtrToParSolverFcn SetupFcn() const
hypre&#39;s internal Setup function
Definition: hypre.hpp:1044
virtual HYPRE_PtrToParSolverFcn SolveFcn() const
hypre&#39;s internal Solve function
Definition: hypre.hpp:1324
int relax_times
Number of relaxation sweeps.
Definition: hypre.hpp:860
void GetEigenvalues(Array< double > &eigenvalues) const
Collect the converged eigenvalues.
Definition: hypre.cpp:5766
void SetCycleType(int cycle_type)
Expert option - consult hypre documentation/team.
Definition: hypre.hpp:1498
void operator*=(double s)
Scale all entries by s: A_scaled = s*A.
Definition: hypre.cpp:1887
void SetPrintLevel(int logging)
Definition: hypre.cpp:5717
double ParNormlp(const Vector &vec, double p, MPI_Comm comm)
Compute the l_p norm of the Vector which is split without overlap across the given communicator...
Definition: hypre.cpp:327
int eig_est_cg_iter
Number of CG iterations to determine eigenvalue estimates.
Definition: hypre.hpp:882
void SetOwnership(int own)
Sets ownership of the internal hypre_ParVector.
Definition: hypre.hpp:183
Abstract class for hypre&#39;s solvers and preconditioners.
Definition: hypre.hpp:970
Flexible GMRES solver in hypre.
Definition: hypre.hpp:1170
ErrorMode error_mode
How to treat hypre errors.
Definition: hypre.hpp:994
HYPRE_BigInt NNZ() const
Returns the global number of nonzeros.
Definition: hypre.hpp:503
const HypreParVector & GetEigenvector(unsigned int i) const
Extract a single eigenvector.
Definition: hypre.cpp:5493
HYPRE_BigInt * GetColStarts() const
Return the parallel column partitioning array.
Definition: hypre.hpp:588
virtual HYPRE_PtrToParSolverFcn SetupFcn() const
ILU Setup function.
Definition: hypre.hpp:1377
const HYPRE_BigInt * Partitioning() const
Returns the parallel row/column partitioning.
Definition: hypre.hpp:162
HypreParVector & operator=(double d)
Set constant values.
Definition: hypre.cpp:192
const double alpha
Definition: ex15.cpp:369
Vector data type.
Definition: vector.hpp:60
void SetAbsTol(double tol)
Definition: hypre.cpp:3759
void PrintCommPkg(std::ostream &out=mfem::out) const
Print information about the hypre_ParCSRCommPkg of the HypreParMatrix.
Definition: hypre.cpp:2206
HypreParMatrix & operator=(double value)
Initialize all entries with value.
Definition: hypre.hpp:639
HypreParMatrix * ExtractSubmatrix(const Array< int > &indices, double threshold=0.0) const
Definition: hypre.cpp:1490
hypre_ParCSRMatrix * StealData()
Changes the ownership of the matrix.
Definition: hypre.cpp:1231
void SetPrecondUsageMode(int pcg_mode)
Definition: hypre.cpp:5409
void SetOperator(const HypreParMatrix &A)
Definition: hypre.cpp:5732
void SetSystemsOptions(int dim, bool order_bynodes=false)
Definition: hypre.cpp:4466
virtual ~HypreParMatrix()
Calls hypre&#39;s destroy function.
Definition: hypre.hpp:777
ID for class HypreParMatrix.
Definition: operator.hpp:256
HYPRE_Int Mult(HypreParVector &x, HypreParVector &y, double alpha=1.0, double beta=0.0) const
Computes y = alpha * A * x + beta * y.
Definition: hypre.cpp:1529
void DropSmallEntries(double tol)
Wrapper for hypre_ParCSRMatrixDropSmallEntries in different versions of hypre. Drop off-diagonal entr...
Definition: hypre.cpp:2006
void SetAggressiveCoarsening(int num_levels)
Expert option - consult hypre documentation/team.
Definition: hypre.hpp:1516
void SetMaxIter(int max_iter)
Definition: hypre.cpp:3764
RefCoord s[3]
const HYPRE_BigInt * RowPart() const
Returns the row partitioning (const version)
Definition: hypre.hpp:515
double omega
SOR parameter (usually in (0,2))
Definition: hypre.hpp:864
HypreParMatrix * LeftDiagMult(const SparseMatrix &D, HYPRE_BigInt *row_starts=NULL) const
Multiply the HypreParMatrix on the left by a block-diagonal parallel matrix D and return the result a...
Definition: hypre.cpp:1705
Base class for solvers.
Definition: operator.hpp:648
virtual HYPRE_PtrToParSolverFcn SetupFcn() const
hypre&#39;s internal Setup function
Definition: hypre.hpp:1570
void SetNodal(int blocksize)
Expert option - consult hypre documentation/team.
Definition: hypre.hpp:1509
int GetOwnership() const
Gets ownership of the internal hypre_ParVector.
Definition: hypre.hpp:186
virtual HYPRE_PtrToParSolverFcn SetupFcn() const
PCG Setup function.
Definition: hypre.hpp:1105
void WrapHypreParCSRMatrix(hypre_ParCSRMatrix *a, bool owner=true)
Converts hypre&#39;s format to HypreParMatrix.
Definition: hypre.cpp:510
OutStream out(std::cout)
Global stream used by the library for standard output. Initially it uses the same std::streambuf as s...
Definition: globals.hpp:66
Abstract operator.
Definition: operator.hpp:24
Wrapper for hypre&#39;s ParCSR matrix class.
Definition: hypre.hpp:277
virtual HYPRE_PtrToParSolverFcn SetupFcn() const
hypre&#39;s internal Setup function
Definition: hypre.hpp:1322
void HostReadWrite()
Update the internal hypre_ParCSRMatrix object, A, to be on host.
Definition: hypre.hpp:732
void AbsMult(double a, const Vector &x, double b, Vector &y) const
Computes y = a * |A| * x + b * y, using entry-wise absolute values of the matrix A.
Definition: hypre.cpp:1671
Memory< double > auxB
Auxiliary buffers for the case when the input or output arrays in methods like Mult(const Vector &amp;...
Definition: hypre.hpp:850
virtual HYPRE_PtrToParSolverFcn SolveFcn() const
hypre&#39;s internal Solve function
Definition: hypre.hpp:1524
~HypreParVector()
Calls hypre&#39;s destroy function.
Definition: hypre.cpp:298
void SetNumModes(int num_eigs)
Definition: hypre.hpp:1734
virtual void Mult(const HypreParVector &b, HypreParVector &x) const
Solve Ax=b with hypre&#39;s PCG.
Definition: hypre.cpp:3617
HypreParMatrix * HypreParMatrixFromBlocks(Array2D< HypreParMatrix * > &blocks, Array2D< double > *blockCoeff)
Returns a merged hypre matrix constructed from hypre matrix blocks.
Definition: hypre.cpp:2626
HypreParVector * Z
Definition: hypre.hpp:852
void SetSingularProblem()
Set this option when solving a curl-curl problem with zero mass term.
Definition: hypre.hpp:1565
void Solve()
Solve the eigenproblem.
Definition: hypre.cpp:5538
virtual ~HypreFGMRES()
Definition: hypre.cpp:4032
HypreParVector * V
Temporary vectors.
Definition: hypre.hpp:852
void SetType(HypreSmoother::Type type, int relax_times=1)
Set the relaxation type and number of sweeps.
Definition: hypre.cpp:3028
virtual const double * Read(bool on_dev=true) const
Shortcut for mfem::Read(vec.GetMemory(), vec.Size(), on_dev).
Definition: vector.hpp:426
MemoryClass
Memory classes identify sets of memory types.
Definition: mem_manager.hpp:73
void WrapMemoryReadWrite(Memory< double > &mem)
Replace the HypreParVector&#39;s data with the given Memory, mem, and prepare the vector for read and wri...
Definition: hypre.cpp:260
virtual HYPRE_PtrToParSolverFcn SetupFcn() const
GMRES Setup function.
Definition: hypre.hpp:1156
int poly_order
Order of the smoothing polynomial.
Definition: hypre.hpp:866
double lambda
Taubin&#39;s lambda-mu method parameters.
Definition: hypre.hpp:873
void WrapMemoryRead(const Memory< double > &mem)
Replace the HypreParVector&#39;s data with the given Memory, mem, and prepare the vector for read access ...
Definition: hypre.cpp:245
MFEM_DEPRECATED void SetZeroInintialIterate()
deprecated: use SetZeroInitialIterate()
Definition: hypre.hpp:1197
MFEM_DEPRECATED HypreParMatrix * GetData()
Deprecated. Use HypreDiagScale::GetData() const instead.
Definition: hypre.hpp:1251
virtual HYPRE_PtrToParSolverFcn SolveFcn() const
hypre&#39;s internal Solve function
Definition: hypre.hpp:1612
void HypreReadWrite()
Update the internal hypre_ParCSRMatrix object, A, to be in hypre memory space.
Definition: hypre.hpp:750