MFEM  v3.4
Finite element discretization library
hypre.hpp
Go to the documentation of this file.
1 // Copyright (c) 2010, Lawrence Livermore National Security, LLC. Produced at
2 // the Lawrence Livermore National Laboratory. LLNL-CODE-443211. All Rights
3 // reserved. See file COPYRIGHT for details.
4 //
5 // This file is part of the MFEM library. For more information and source code
6 // availability see http://mfem.org.
7 //
8 // MFEM is free software; you can redistribute it and/or modify it under the
9 // terms of the GNU Lesser General Public License (as published by the Free
10 // Software Foundation) version 2.1 dated February 1999.
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 #include "sparsemat.hpp"
36 #include "hypre_parcsr.hpp"
37 #ifdef MFEM_USE_SUNDIALS
38 #include <nvector/nvector_parhyp.h>
39 #endif
40 
41 namespace mfem
42 {
43 
44 class ParFiniteElementSpace;
45 class HypreParMatrix;
46 
47 namespace internal
48 {
49 
50 template <typename int_type>
51 inline int to_int(int_type i)
52 {
53  MFEM_ASSERT(int_type(int(i)) == i, "overflow converting int_type to int");
54  return int(i);
55 }
56 
57 // Specialization for to_int(int)
58 template <> inline int to_int(int i) { return i; }
59 
60 // Convert a HYPRE_Int to int
61 #ifdef HYPRE_BIGINT
62 template <>
63 inline int to_int(HYPRE_Int i)
64 {
65  MFEM_ASSERT(HYPRE_Int(int(i)) == i, "overflow converting HYPRE_Int to int");
66  return int(i);
67 }
68 #endif
69 
70 }
71 
72 /// Wrapper for hypre's parallel vector class
73 class HypreParVector : public Vector
74 {
75 private:
76  int own_ParVector;
77 
78  /// The actual object
79  hypre_ParVector *x;
80 
81  friend class HypreParMatrix;
82 
83  // Set Vector::data and Vector::size from *x
84  inline void _SetDataAndSize_();
85 
86 public:
87  /** Creates vector with given global size and partitioning of the columns.
88  Processor P owns columns [col[P],col[P+1]) */
89  HypreParVector(MPI_Comm comm, HYPRE_Int glob_size, HYPRE_Int *col);
90  /** Creates vector with given global size, partitioning of the columns,
91  and data. The data must be allocated and destroyed outside.
92  If _data is NULL, a dummy vector without a valid data array will
93  be created. */
94  HypreParVector(MPI_Comm comm, HYPRE_Int glob_size, double *_data,
95  HYPRE_Int *col);
96  /// Creates vector compatible with y
98  /// Creates vector compatible with (i.e. in the domain of) A or A^T
99  explicit HypreParVector(const HypreParMatrix &A, int transpose = 0);
100  /// Creates vector wrapping y
101  explicit HypreParVector(HYPRE_ParVector y);
102  /// Create a true dof parallel vector on a given ParFiniteElementSpace
103  explicit HypreParVector(ParFiniteElementSpace *pfes);
104 
105  /// MPI communicator
106  MPI_Comm GetComm() { return x->comm; }
107 
108  /// Returns the row partitioning
109  inline HYPRE_Int *Partitioning() { return x->partitioning; }
110 
111  /// Returns the global number of rows
112  inline HYPRE_Int GlobalSize() { return x->global_size; }
113 
114  /// Typecasting to hypre's hypre_ParVector*
115  operator hypre_ParVector*() const { return x; }
116 #ifndef HYPRE_PAR_VECTOR_STRUCT
117  /// Typecasting to hypre's HYPRE_ParVector, a.k.a. void *
118  operator HYPRE_ParVector() const { return (HYPRE_ParVector) x; }
119 #endif
120  /// Changes the ownership of the the vector
121  hypre_ParVector *StealParVector() { own_ParVector = 0; return x; }
122 
123  /// Sets ownership of the internal hypre_ParVector
124  void SetOwnership(int own) { own_ParVector = own; }
125 
126  /// Gets ownership of the internal hypre_ParVector
127  int GetOwnership() const { return own_ParVector; }
128 
129  /// Returns the global vector in each processor
130  Vector* GlobalVector() const;
131 
132  /// Set constant values
133  HypreParVector& operator= (double d);
134  /// Define '=' for hypre vectors.
136 
137  /** Sets the data of the Vector and the hypre_ParVector to _data.
138  Must be used only for HypreParVectors that do not own the data,
139  e.g. created with the constructor:
140  HypreParVector(MPI_Comm comm, HYPRE_Int glob_size, double *_data,
141  HYPRE_Int *col). */
142  void SetData(double *_data);
143 
144  /// Set random values
145  HYPRE_Int Randomize(HYPRE_Int seed);
146 
147  /// Prints the locally owned rows in parallel
148  void Print(const char *fname) const;
149 
150  /// Calls hypre's destroy function
151  ~HypreParVector();
152 
153 #ifdef MFEM_USE_SUNDIALS
154  /// Return a new wrapper SUNDIALS N_Vector of type SUNDIALS_NVEC_PARHYP.
155  /** The returned N_Vector must be destroyed by the caller. */
156  virtual N_Vector ToNVector() { return N_VMake_ParHyp(x); }
157 
158  /** @brief Update an existing wrapper SUNDIALS N_Vector of type
159  SUNDIALS_NVEC_PARHYP to point to this Vector. */
160  virtual void ToNVector(N_Vector &nv);
161 #endif
162 };
163 
164 /// Returns the inner product of x and y
165 double InnerProduct(HypreParVector &x, HypreParVector &y);
166 double InnerProduct(HypreParVector *x, HypreParVector *y);
167 
168 
169 /** @brief Compute the l_p norm of the Vector which is split without overlap
170  across the given communicator. */
171 double ParNormlp(const Vector &vec, double p, MPI_Comm comm);
172 
173 
174 /// Wrapper for hypre's ParCSR matrix class
175 class HypreParMatrix : public Operator
176 {
177 private:
178  /// The actual object
179  hypre_ParCSRMatrix *A;
180 
181  /// Auxiliary vectors for typecasting
182  mutable HypreParVector *X, *Y;
183 
184  // Flags indicating ownership of A->diag->{i,j,data}, A->offd->{i,j,data},
185  // and A->col_map_offd.
186  // The possible values for diagOwner are:
187  // -1: no special treatment of A->diag (default)
188  // 0: prevent hypre from destroying A->diag->{i,j,data}
189  // 1: same as 0, plus take ownership of A->diag->{i,j}
190  // 2: same as 0, plus take ownership of A->diag->data
191  // 3: same as 0, plus take ownership of A->diag->{i,j,data}
192  // The same values and rules apply to offdOwner and A->offd.
193  // The possible values for colMapOwner are:
194  // -1: no special treatment of A->col_map_offd (default)
195  // 0: prevent hypre from destroying A->col_map_offd
196  // 1: same as 0, plus take ownership of A->col_map_offd
197  // All owned arrays are destroyed with 'delete []'.
198  signed char diagOwner, offdOwner, colMapOwner;
199 
200  // Does the object own the pointer A?
201  signed char ParCSROwner;
202 
203  // Initialize with defaults. Does not initialize inherited members.
204  void Init();
205 
206  // Delete all owned data. Does not perform re-initialization with defaults.
207  void Destroy();
208 
209  // Copy (shallow or deep, based on HYPRE_BIGINT) the I and J arrays from csr
210  // to hypre_csr. Shallow copy the data. Return the appropriate ownership
211  // flag.
212  static char CopyCSR(SparseMatrix *csr, hypre_CSRMatrix *hypre_csr);
213  // Copy (shallow or deep, based on HYPRE_BIGINT) the I and J arrays from
214  // bool_csr to hypre_csr. Allocate the data array and set it to all ones.
215  // Return the appropriate ownership flag.
216  static char CopyBoolCSR(Table *bool_csr, hypre_CSRMatrix *hypre_csr);
217 
218  // Copy the j array of a hypre_CSRMatrix to the given J array, converting
219  // the indices from HYPRE_Int to int.
220  static void CopyCSR_J(hypre_CSRMatrix *hypre_csr, int *J);
221 
222 public:
223  /// An empty matrix to be used as a reference to an existing matrix
224  HypreParMatrix();
225 
226  /// Converts hypre's format to HypreParMatrix
227  /** If @a owner is false, ownership of @a a is not transferred */
228  explicit HypreParMatrix(hypre_ParCSRMatrix *a, bool owner = true)
229  {
230  Init();
231  A = a;
232  if (!owner) { ParCSROwner = 0; }
233  height = GetNumRows();
234  width = GetNumCols();
235  }
236 
237  /** Creates block-diagonal square parallel matrix. Diagonal is given by diag
238  which must be in CSR format (finalized). The new HypreParMatrix does not
239  take ownership of any of the input arrays.
240  @warning The ordering of the columns in each row in @a *diag may be
241  changed by this constructor to ensure that the first entry in each row is
242  the diagonal one. This is expected by most hypre functions. */
243  HypreParMatrix(MPI_Comm comm, HYPRE_Int glob_size, HYPRE_Int *row_starts,
244  SparseMatrix *diag); // constructor with 4 arguments, v1
245 
246  /** Creates block-diagonal rectangular parallel matrix. Diagonal is given by
247  diag which must be in CSR format (finalized). The new HypreParMatrix does
248  not take ownership of any of the input arrays. */
249  HypreParMatrix(MPI_Comm comm, HYPRE_Int global_num_rows,
250  HYPRE_Int global_num_cols, HYPRE_Int *row_starts,
251  HYPRE_Int *col_starts,
252  SparseMatrix *diag); // constructor with 6 arguments, v1
253 
254  /** Creates general (rectangular) parallel matrix. The new HypreParMatrix
255  does not take ownership of any of the input arrays. */
256  HypreParMatrix(MPI_Comm comm, HYPRE_Int global_num_rows,
257  HYPRE_Int global_num_cols, HYPRE_Int *row_starts,
258  HYPRE_Int *col_starts, SparseMatrix *diag, SparseMatrix *offd,
259  HYPRE_Int *cmap); // constructor with 8 arguments
260 
261  /** Creates general (rectangular) parallel matrix. The new HypreParMatrix
262  takes ownership of all input arrays, except col_starts and row_starts. */
263  HypreParMatrix(MPI_Comm comm,
264  HYPRE_Int global_num_rows, HYPRE_Int global_num_cols,
265  HYPRE_Int *row_starts, HYPRE_Int *col_starts,
266  HYPRE_Int *diag_i, HYPRE_Int *diag_j, double *diag_data,
267  HYPRE_Int *offd_i, HYPRE_Int *offd_j, double *offd_data,
268  HYPRE_Int offd_num_cols,
269  HYPRE_Int *offd_col_map); // constructor with 13 arguments
270 
271  /// Creates a parallel matrix from SparseMatrix on processor 0.
272  HypreParMatrix(MPI_Comm comm, HYPRE_Int *row_starts, HYPRE_Int *col_starts,
273  SparseMatrix *a); // constructor with 4 arguments, v2
274 
275  /** Creates boolean block-diagonal rectangular parallel matrix. The new
276  HypreParMatrix does not take ownership of any of the input arrays. */
277  HypreParMatrix(MPI_Comm comm, HYPRE_Int global_num_rows,
278  HYPRE_Int global_num_cols, HYPRE_Int *row_starts,
279  HYPRE_Int *col_starts,
280  Table *diag); // constructor with 6 arguments, v2
281 
282  /** Creates boolean rectangular parallel matrix. The new HypreParMatrix takes
283  ownership of the arrays i_diag, j_diag, i_offd, j_offd, and cmap; does
284  not take ownership of the arrays row and col. */
285  HypreParMatrix(MPI_Comm comm, int id, int np, HYPRE_Int *row, HYPRE_Int *col,
286  HYPRE_Int *i_diag, HYPRE_Int *j_diag, HYPRE_Int *i_offd,
287  HYPRE_Int *j_offd, HYPRE_Int *cmap,
288  HYPRE_Int cmap_size); // constructor with 11 arguments
289 
290  /** Creates a general parallel matrix from a local CSR matrix on each
291  processor described by the I, J and data arrays. The local matrix should
292  be of size (local) nrows by (global) glob_ncols. The new parallel matrix
293  contains copies of all input arrays (so they can be deleted). */
294  HypreParMatrix(MPI_Comm comm, int nrows, HYPRE_Int glob_nrows,
295  HYPRE_Int glob_ncols, int *I, HYPRE_Int *J,
296  double *data, HYPRE_Int *rows,
297  HYPRE_Int *cols); // constructor with 9 arguments
298 
299  /// Make this HypreParMatrix a reference to 'master'
300  void MakeRef(const HypreParMatrix &master);
301 
302  /// MPI communicator
303  MPI_Comm GetComm() const { return A->comm; }
304 
305  /// Typecasting to hypre's hypre_ParCSRMatrix*
306  operator hypre_ParCSRMatrix*() const { return A; }
307 #ifndef HYPRE_PAR_CSR_MATRIX_STRUCT
308  /// Typecasting to hypre's HYPRE_ParCSRMatrix, a.k.a. void *
309  operator HYPRE_ParCSRMatrix() { return (HYPRE_ParCSRMatrix) A; }
310 #endif
311  /// Changes the ownership of the the matrix
312  hypre_ParCSRMatrix* StealData();
313 
314  /// Explicitly set the three ownership flags, see docs for diagOwner etc.
315  void SetOwnerFlags(signed char diag, signed char offd, signed char colmap)
316  { diagOwner = diag, offdOwner = offd, colMapOwner = colmap; }
317 
318  /// Get diag ownership flag
319  signed char OwnsDiag() const { return diagOwner; }
320  /// Get offd ownership flag
321  signed char OwnsOffd() const { return offdOwner; }
322  /// Get colmap ownership flag
323  signed char OwnsColMap() const { return colMapOwner; }
324 
325  /** If the HypreParMatrix does not own the row-starts array, make a copy of
326  it that the HypreParMatrix will own. If the col-starts array is the same
327  as the row-starts array, col-starts is also replaced. */
328  void CopyRowStarts();
329  /** If the HypreParMatrix does not own the col-starts array, make a copy of
330  it that the HypreParMatrix will own. If the row-starts array is the same
331  as the col-starts array, row-starts is also replaced. */
332  void CopyColStarts();
333 
334  /// Returns the global number of nonzeros
335  inline HYPRE_Int NNZ() const { return A->num_nonzeros; }
336  /// Returns the row partitioning
337  inline HYPRE_Int *RowPart() { return A->row_starts; }
338  /// Returns the column partitioning
339  inline HYPRE_Int *ColPart() { return A->col_starts; }
340  /// Returns the row partitioning (const version)
341  inline const HYPRE_Int *RowPart() const { return A->row_starts; }
342  /// Returns the column partitioning (const version)
343  inline const HYPRE_Int *ColPart() const { return A->col_starts; }
344  /// Returns the global number of rows
345  inline HYPRE_Int M() const { return A->global_num_rows; }
346  /// Returns the global number of columns
347  inline HYPRE_Int N() const { return A->global_num_cols; }
348 
349  /// Get the local diagonal of the matrix.
350  void GetDiag(Vector &diag) const;
351  /// Get the local diagonal block. NOTE: 'diag' will not own any data.
352  void GetDiag(SparseMatrix &diag) const;
353  /// Get the local off-diagonal block. NOTE: 'offd' will not own any data.
354  void GetOffd(SparseMatrix &offd, HYPRE_Int* &cmap) const;
355 
356  /** Split the matrix into M x N equally sized blocks of parallel matrices.
357  The size of 'blocks' must already be set to M x N. */
358  void GetBlocks(Array2D<HypreParMatrix*> &blocks,
359  bool interleaved_rows = false,
360  bool interleaved_cols = false) const;
361 
362  /// Returns the transpose of *this
363  HypreParMatrix * Transpose() const;
364 
365  /// Returns the number of rows in the diagonal block of the ParCSRMatrix
366  int GetNumRows() const
367  {
368  return internal::to_int(
369  hypre_CSRMatrixNumRows(hypre_ParCSRMatrixDiag(A)));
370  }
371 
372  /// Returns the number of columns in the diagonal block of the ParCSRMatrix
373  int GetNumCols() const
374  {
375  return internal::to_int(
376  hypre_CSRMatrixNumCols(hypre_ParCSRMatrixDiag(A)));
377  }
378 
379  HYPRE_Int GetGlobalNumRows() const
380  { return hypre_ParCSRMatrixGlobalNumRows(A); }
381 
382  HYPRE_Int GetGlobalNumCols() const
383  { return hypre_ParCSRMatrixGlobalNumCols(A); }
384 
385  HYPRE_Int *GetRowStarts() const { return hypre_ParCSRMatrixRowStarts(A); }
386 
387  HYPRE_Int *GetColStarts() const { return hypre_ParCSRMatrixColStarts(A); }
388 
389  /// Computes y = alpha * A * x + beta * y
390  HYPRE_Int Mult(HypreParVector &x, HypreParVector &y,
391  double alpha = 1.0, double beta = 0.0);
392  /// Computes y = alpha * A * x + beta * y
393  HYPRE_Int Mult(HYPRE_ParVector x, HYPRE_ParVector y,
394  double alpha = 1.0, double beta = 0.0);
395  /// Computes y = alpha * A^t * x + beta * y
397  double alpha = 1.0, double beta = 0.0);
398 
399  void Mult(double a, const Vector &x, double b, Vector &y) const;
400  void MultTranspose(double a, const Vector &x, double b, Vector &y) const;
401 
402  virtual void Mult(const Vector &x, Vector &y) const
403  { Mult(1.0, x, 0.0, y); }
404  virtual void MultTranspose(const Vector &x, Vector &y) const
405  { MultTranspose(1.0, x, 0.0, y); }
406 
407  /** The "Boolean" analog of y = alpha * A * x + beta * y, where elements in
408  the sparsity pattern of the matrix are treated as "true". */
409  void BooleanMult(int alpha, int *x, int beta, int *y)
410  {
411  internal::hypre_ParCSRMatrixBooleanMatvec(A, alpha, x, beta, y);
412  }
413 
414  /** The "Boolean" analog of y = alpha * A^T * x + beta * y, where elements in
415  the sparsity pattern of the matrix are treated as "true". */
416  void BooleanMultTranspose(int alpha, int *x, int beta, int *y)
417  {
418  internal::hypre_ParCSRMatrixBooleanMatvecT(A, alpha, x, beta, y);
419  }
420 
421  /// Initialize all entries with value.
422  HypreParMatrix &operator=(double value)
423  { internal::hypre_ParCSRMatrixSetConstantValues(A, value); return *this; }
424 
425  /** Perform the operation `*this += B`, assuming that both matrices use the
426  same row and column partitions and the same col_map_offd arrays. We also
427  assume that the sparsity pattern of `*this` contains that of `B`. */
428  HypreParMatrix &operator+=(const HypreParMatrix &B) { return Add(1.0, B); }
429 
430  /** Perform the operation `*this += beta*B`, assuming that both matrices use
431  the same row and column partitions and the same col_map_offd arrays. We
432  also assume that the sparsity pattern of `*this` contains that of `B`. */
433  HypreParMatrix &Add(const double beta, const HypreParMatrix &B)
434  {
435  MFEM_VERIFY(internal::hypre_ParCSRMatrixSum(A, beta, B.A) == 0,
436  "error in hypre_ParCSRMatrixSum");
437  return *this;
438  }
439 
440  /** @brief Multiply the HypreParMatrix on the left by a block-diagonal
441  parallel matrix @a D and return the result as a new HypreParMatrix. */
442  /** If @a D has a different number of rows than @a A (this matrix), @a D's
443  row starts array needs to be given (as returned by the methods
444  GetDofOffsets/GetTrueDofOffsets of ParFiniteElementSpace). The new
445  matrix @a D*A uses copies of the row- and column-starts arrays, so "this"
446  matrix and @a row_starts can be deleted.
447  @note This operation is local and does not require communication. */
449  HYPRE_Int* row_starts = NULL) const;
450 
451  /// Scale the local row i by s(i).
452  void ScaleRows(const Vector & s);
453  /// Scale the local row i by 1./s(i)
454  void InvScaleRows(const Vector & s);
455  /// Scale all entries by s: A_scaled = s*A.
456  void operator*=(double s);
457 
458  /// Remove values smaller in absolute value than some threshold
459  void Threshold(double threshold = 0.0);
460 
461  /// If a row contains only zeros, set its diagonal to 1.
462  void EliminateZeroRows() { hypre_ParCSRMatrixFixZeroRows(A); }
463 
464  /** Eliminate rows and columns from the matrix, and rows from the vector B.
465  Modify B with the BC values in X. */
466  void EliminateRowsCols(const Array<int> &rows_cols, const HypreParVector &X,
467  HypreParVector &B);
468 
469  /** Eliminate rows and columns from the matrix and store the eliminated
470  elements in a new matrix Ae (returned), so that the modified matrix and
471  Ae sum to the original matrix. */
472  HypreParMatrix* EliminateRowsCols(const Array<int> &rows_cols);
473 
474  /// Prints the locally owned rows in parallel
475  void Print(const char *fname, HYPRE_Int offi = 0, HYPRE_Int offj = 0);
476  /// Reads the matrix from a file
477  void Read(MPI_Comm comm, const char *fname);
478  /// Read a matrix saved as a HYPRE_IJMatrix
479  void Read_IJMatrix(MPI_Comm comm, const char *fname);
480 
481  /// Print information about the hypre_ParCSRCommPkg of the HypreParMatrix.
482  void PrintCommPkg(std::ostream &out = mfem::out) const;
483 
484  /// Calls hypre's destroy function
485  virtual ~HypreParMatrix() { Destroy(); }
486 
487  Type GetType() const { return Hypre_ParCSR; }
488 };
489 
490 /** @brief Return a new matrix `C = alpha*A + beta*B`, assuming that both `A`
491  and `B` use the same row and column partitions and the same `col_map_offd`
492  arrays. */
493 HypreParMatrix *Add(double alpha, const HypreParMatrix &A,
494  double beta, const HypreParMatrix &B);
495 
496 /// Returns the matrix A * B
497 HypreParMatrix * ParMult(const HypreParMatrix *A, const HypreParMatrix *B);
498 /// Returns the matrix A + B
499 /** It is assumed that both matrices use the same row and column partitions and
500  the same col_map_offd arrays. */
501 HypreParMatrix * ParAdd(const HypreParMatrix *A, const HypreParMatrix *B);
502 
503 /// Returns the matrix P^t * A * P
504 HypreParMatrix * RAP(const HypreParMatrix *A, const HypreParMatrix *P);
505 /// Returns the matrix Rt^t * A * P
506 HypreParMatrix * RAP(const HypreParMatrix * Rt, const HypreParMatrix *A,
507  const HypreParMatrix *P);
508 
509 /** Eliminate essential BC specified by 'ess_dof_list' from the solution X to
510  the r.h.s. B. Here A is a matrix with eliminated BC, while Ae is such that
511  (A+Ae) is the original (Neumann) matrix before elimination. */
512 void EliminateBC(HypreParMatrix &A, HypreParMatrix &Ae,
513  const Array<int> &ess_dof_list, const Vector &X, Vector &B);
514 
515 
516 /// Parallel smoothers in hypre
517 class HypreSmoother : public Solver
518 {
519 protected:
520  /// The linear system matrix
522  /// Right-hand side and solution vectors
523  mutable HypreParVector *B, *X;
524  /// Temporary vectors
525  mutable HypreParVector *V, *Z;
526  /// FIR Filter Temporary Vectors
527  mutable HypreParVector *X0, *X1;
528 
529  /** Smoother type from hypre_ParCSRRelax() in ams.c plus extensions, see the
530  enumeration Type below. */
531  int type;
532  /// Number of relaxation sweeps
534  /// Damping coefficient (usually <= 1)
535  double relax_weight;
536  /// SOR parameter (usually in (0,2))
537  double omega;
538  /// Order of the smoothing polynomial
540  /// Fraction of spectrum to smooth for polynomial relaxation
542  /// Apply the polynomial smoother to A or D^{-1/2} A D^{-1/2}
544 
545  /// Taubin's lambda-mu method parameters
546  double lambda;
547  double mu;
549 
550  /// l1 norms of the rows of A
551  double *l1_norms;
552  /// If set, take absolute values of the computed l1_norms
554  /// Maximal eigenvalue estimate for polynomial smoothing
555  double max_eig_est;
556  /// Minimal eigenvalue estimate for polynomial smoothing
557  double min_eig_est;
558  /// Parameters for windowing function of FIR filter
559  double window_params[3];
560 
561  /// Combined coefficients for windowing and Chebyshev polynomials.
562  double* fir_coeffs;
563 
564 public:
565  /** Hypre smoother types:
566  0 = Jacobi
567  1 = l1-scaled Jacobi
568  2 = l1-scaled block Gauss-Seidel/SSOR
569  4 = truncated l1-scaled block Gauss-Seidel/SSOR
570  5 = lumped Jacobi
571  6 = Gauss-Seidel
572  16 = Chebyshev
573  1001 = Taubin polynomial smoother
574  1002 = FIR polynomial smoother. */
575  enum Type { Jacobi = 0, l1Jacobi = 1, l1GS = 2, l1GStr = 4, lumpedJacobi = 5,
576  GS = 6, Chebyshev = 16, Taubin = 1001, FIR = 1002
577  };
578 
579  HypreSmoother();
580 
582  int relax_times = 1, double relax_weight = 1.0,
583  double omega = 1.0, int poly_order = 2,
584  double poly_fraction = .3);
585 
586  /// Set the relaxation type and number of sweeps
588  /// Set SOR-related parameters
589  void SetSOROptions(double relax_weight, double omega);
590  /// Set parameters for polynomial smoothing
591  void SetPolyOptions(int poly_order, double poly_fraction);
592  /// Set parameters for Taubin's lambda-mu method
593  void SetTaubinOptions(double lambda, double mu, int iter);
594 
595  /// Convenience function for setting canonical windowing parameters
596  void SetWindowByName(const char* window_name);
597  /// Set parameters for windowing function for FIR smoother.
598  void SetWindowParameters(double a, double b, double c);
599  /// Compute window and Chebyshev coefficients for given polynomial order.
600  void SetFIRCoefficients(double max_eig);
601 
602  /// After computing l1-norms, replace them with their absolute values.
603  /** By default, the l1-norms take their sign from the corresponding diagonal
604  entries in the associated matrix. */
605  void SetPositiveDiagonal(bool pos = true) { pos_l1_norms = pos; }
606 
607  /** Set/update the associated operator. Must be called after setting the
608  HypreSmoother type and options. */
609  virtual void SetOperator(const Operator &op);
610 
611  /// Relax the linear system Ax=b
612  virtual void Mult(const HypreParVector &b, HypreParVector &x) const;
613  virtual void Mult(const Vector &b, Vector &x) const;
614 
615  virtual ~HypreSmoother();
616 };
617 
618 
619 /// Abstract class for hypre's solvers and preconditioners
620 class HypreSolver : public Solver
621 {
622 protected:
623  /// The linear system matrix
625 
626  /// Right-hand side and solution vector
627  mutable HypreParVector *B, *X;
628 
629  /// Was hypre's Setup function called already?
630  mutable int setup_called;
631 
632 public:
633  HypreSolver();
634 
636 
637  /// Typecast to HYPRE_Solver -- return the solver
638  virtual operator HYPRE_Solver() const = 0;
639 
640  /// hypre's internal Setup function
641  virtual HYPRE_PtrToParSolverFcn SetupFcn() const = 0;
642  /// hypre's internal Solve function
643  virtual HYPRE_PtrToParSolverFcn SolveFcn() const = 0;
644 
645  virtual void SetOperator(const Operator &op)
646  { mfem_error("HypreSolvers do not support SetOperator!"); }
647 
648  /// Solve the linear system Ax=b
649  virtual void Mult(const HypreParVector &b, HypreParVector &x) const;
650  virtual void Mult(const Vector &b, Vector &x) const;
651 
652  virtual ~HypreSolver();
653 };
654 
655 /// PCG solver in hypre
656 class HyprePCG : public HypreSolver
657 {
658 private:
659  HYPRE_Solver pcg_solver;
660 
661 public:
663 
664  void SetTol(double tol);
665  void SetMaxIter(int max_iter);
666  void SetLogging(int logging);
667  void SetPrintLevel(int print_lvl);
668 
669  /// Set the hypre solver to be used as a preconditioner
670  void SetPreconditioner(HypreSolver &precond);
671 
672  /** Use the L2 norm of the residual for measuring PCG convergence, plus
673  (optionally) 1) periodically recompute true residuals from scratch; and
674  2) enable residual-based stopping criteria. */
675  void SetResidualConvergenceOptions(int res_frequency=-1, double rtol=0.0);
676 
677  /// non-hypre setting
679 
680  void GetNumIterations(int &num_iterations)
681  {
682  HYPRE_Int num_it;
683  HYPRE_ParCSRPCGGetNumIterations(pcg_solver, &num_it);
684  num_iterations = internal::to_int(num_it);
685  }
686 
687  /// The typecast to HYPRE_Solver returns the internal pcg_solver
688  virtual operator HYPRE_Solver() const { return pcg_solver; }
689 
690  /// PCG Setup function
691  virtual HYPRE_PtrToParSolverFcn SetupFcn() const
692  { return (HYPRE_PtrToParSolverFcn) HYPRE_ParCSRPCGSetup; }
693  /// PCG Solve function
694  virtual HYPRE_PtrToParSolverFcn SolveFcn() const
695  { return (HYPRE_PtrToParSolverFcn) HYPRE_ParCSRPCGSolve; }
696 
697  /// Solve Ax=b with hypre's PCG
698  virtual void Mult(const HypreParVector &b, HypreParVector &x) const;
699  using HypreSolver::Mult;
700 
701  virtual ~HyprePCG();
702 };
703 
704 /// GMRES solver in hypre
705 class HypreGMRES : public HypreSolver
706 {
707 private:
708  HYPRE_Solver gmres_solver;
709 
710 public:
712 
713  void SetTol(double tol);
714  void SetMaxIter(int max_iter);
715  void SetKDim(int dim);
716  void SetLogging(int logging);
717  void SetPrintLevel(int print_lvl);
718 
719  /// Set the hypre solver to be used as a preconditioner
720  void SetPreconditioner(HypreSolver &precond);
721 
722  /// non-hypre setting
724 
725  /// The typecast to HYPRE_Solver returns the internal gmres_solver
726  virtual operator HYPRE_Solver() const { return gmres_solver; }
727 
728  /// GMRES Setup function
729  virtual HYPRE_PtrToParSolverFcn SetupFcn() const
730  { return (HYPRE_PtrToParSolverFcn) HYPRE_ParCSRGMRESSetup; }
731  /// GMRES Solve function
732  virtual HYPRE_PtrToParSolverFcn SolveFcn() const
733  { return (HYPRE_PtrToParSolverFcn) HYPRE_ParCSRGMRESSolve; }
734 
735  /// Solve Ax=b with hypre's GMRES
736  virtual void Mult (const HypreParVector &b, HypreParVector &x) const;
737  using HypreSolver::Mult;
738 
739  virtual ~HypreGMRES();
740 };
741 
742 /// The identity operator as a hypre solver
744 {
745 public:
746  virtual operator HYPRE_Solver() const { return NULL; }
747 
748  virtual HYPRE_PtrToParSolverFcn SetupFcn() const
749  { return (HYPRE_PtrToParSolverFcn) hypre_ParKrylovIdentitySetup; }
750  virtual HYPRE_PtrToParSolverFcn SolveFcn() const
751  { return (HYPRE_PtrToParSolverFcn) hypre_ParKrylovIdentity; }
752 
753  virtual ~HypreIdentity() { }
754 };
755 
756 /// Jacobi preconditioner in hypre
758 {
759 public:
762  virtual operator HYPRE_Solver() const { return NULL; }
763 
764  virtual HYPRE_PtrToParSolverFcn SetupFcn() const
765  { return (HYPRE_PtrToParSolverFcn) HYPRE_ParCSRDiagScaleSetup; }
766  virtual HYPRE_PtrToParSolverFcn SolveFcn() const
767  { return (HYPRE_PtrToParSolverFcn) HYPRE_ParCSRDiagScale; }
768 
769  HypreParMatrix* GetData() { return A; }
770  virtual ~HypreDiagScale() { }
771 };
772 
773 /// The ParaSails preconditioner in hypre
775 {
776 private:
777  HYPRE_Solver sai_precond;
778 
779 public:
781 
782  void SetSymmetry(int sym);
783 
784  /// The typecast to HYPRE_Solver returns the internal sai_precond
785  virtual operator HYPRE_Solver() const { return sai_precond; }
786 
787  virtual HYPRE_PtrToParSolverFcn SetupFcn() const
788  { return (HYPRE_PtrToParSolverFcn) HYPRE_ParaSailsSetup; }
789  virtual HYPRE_PtrToParSolverFcn SolveFcn() const
790  { return (HYPRE_PtrToParSolverFcn) HYPRE_ParaSailsSolve; }
791 
792  virtual ~HypreParaSails();
793 };
794 
795 /// The BoomerAMG solver in hypre
797 {
798 private:
799  HYPRE_Solver amg_precond;
800 
801  /// Rigid body modes
803 
804  /// Finite element space for elasticity problems, see SetElasticityOptions()
805  ParFiniteElementSpace *fespace;
806 
807  /// Recompute the rigid-body modes vectors (in the rbms array)
808  void RecomputeRBMs();
809 
810  /// Default, generally robust, BoomerAMG options
811  void SetDefaultOptions();
812 
813  // If amg_precond is NULL, allocates it and sets default options.
814  // Otherwise saves the options from amg_precond, destroys it, allocates a new
815  // one, and sets its options to the saved values.
816  void ResetAMGPrecond();
817 
818 public:
819  HypreBoomerAMG();
820 
822 
823  virtual void SetOperator(const Operator &op);
824 
825  /** More robust options for systems, such as elasticity. Note that BoomerAMG
826  assumes Ordering::byVDIM in the finite element space used to generate the
827  matrix A. */
828  void SetSystemsOptions(int dim);
829 
830  /** A special elasticity version of BoomerAMG that takes advantage of
831  geometric rigid body modes and could perform better on some problems, see
832  "Improving algebraic multigrid interpolation operators for linear
833  elasticity problems", Baker, Kolev, Yang, NLAA 2009, DOI:10.1002/nla.688.
834  As with SetSystemsOptions(), this solver assumes Ordering::byVDIM. */
836 
837  void SetPrintLevel(int print_level)
838  { HYPRE_BoomerAMGSetPrintLevel(amg_precond, print_level); }
839 
840  /// The typecast to HYPRE_Solver returns the internal amg_precond
841  virtual operator HYPRE_Solver() const { return amg_precond; }
842 
843  virtual HYPRE_PtrToParSolverFcn SetupFcn() const
844  { return (HYPRE_PtrToParSolverFcn) HYPRE_BoomerAMGSetup; }
845  virtual HYPRE_PtrToParSolverFcn SolveFcn() const
846  { return (HYPRE_PtrToParSolverFcn) HYPRE_BoomerAMGSolve; }
847 
848  virtual ~HypreBoomerAMG();
849 };
850 
851 /// Compute the discrete gradient matrix between the nodal linear and ND1 spaces
852 HypreParMatrix* DiscreteGrad(ParFiniteElementSpace *edge_fespace,
853  ParFiniteElementSpace *vert_fespace);
854 /// Compute the discrete curl matrix between the ND1 and RT0 spaces
855 HypreParMatrix* DiscreteCurl(ParFiniteElementSpace *face_fespace,
856  ParFiniteElementSpace *edge_fespace);
857 
858 /// The Auxiliary-space Maxwell Solver in hypre
859 class HypreAMS : public HypreSolver
860 {
861 private:
862  HYPRE_Solver ams;
863 
864  /// Vertex coordinates
865  HypreParVector *x, *y, *z;
866  /// Discrete gradient matrix
867  HypreParMatrix *G;
868  /// Nedelec interpolation matrix and its components
869  HypreParMatrix *Pi, *Pix, *Piy, *Piz;
870 
871 public:
873 
874  void SetPrintLevel(int print_lvl);
875 
876  /// Set this option when solving a curl-curl problem with zero mass term
877  void SetSingularProblem() { HYPRE_AMSSetBetaPoissonMatrix(ams, NULL); }
878 
879  /// The typecast to HYPRE_Solver returns the internal ams object
880  virtual operator HYPRE_Solver() const { return ams; }
881 
882  virtual HYPRE_PtrToParSolverFcn SetupFcn() const
883  { return (HYPRE_PtrToParSolverFcn) HYPRE_AMSSetup; }
884  virtual HYPRE_PtrToParSolverFcn SolveFcn() const
885  { return (HYPRE_PtrToParSolverFcn) HYPRE_AMSSolve; }
886 
887  virtual ~HypreAMS();
888 };
889 
890 /// The Auxiliary-space Divergence Solver in hypre
891 class HypreADS : public HypreSolver
892 {
893 private:
894  HYPRE_Solver ads;
895 
896  /// Vertex coordinates
897  HypreParVector *x, *y, *z;
898  /// Discrete gradient matrix
899  HypreParMatrix *G;
900  /// Discrete curl matrix
901  HypreParMatrix *C;
902  /// Nedelec interpolation matrix and its components
903  HypreParMatrix *ND_Pi, *ND_Pix, *ND_Piy, *ND_Piz;
904  /// Raviart-Thomas interpolation matrix and its components
905  HypreParMatrix *RT_Pi, *RT_Pix, *RT_Piy, *RT_Piz;
906 
907 public:
909 
910  void SetPrintLevel(int print_lvl);
911 
912  /// The typecast to HYPRE_Solver returns the internal ads object
913  virtual operator HYPRE_Solver() const { return ads; }
914 
915  virtual HYPRE_PtrToParSolverFcn SetupFcn() const
916  { return (HYPRE_PtrToParSolverFcn) HYPRE_ADSSetup; }
917  virtual HYPRE_PtrToParSolverFcn SolveFcn() const
918  { return (HYPRE_PtrToParSolverFcn) HYPRE_ADSSolve; }
919 
920  virtual ~HypreADS();
921 };
922 
923 /** LOBPCG eigenvalue solver in hypre
924 
925  The Locally Optimal Block Preconditioned Conjugate Gradient (LOBPCG)
926  eigenvalue solver is designed to find the lowest eigenmodes of the
927  generalized eigenvalue problem:
928  A x = lambda M x
929  where A is symmetric, potentially indefinite and M is symmetric positive
930  definite. The eigenvectors are M-orthonormal, meaning that
931  x^T M x = 1 and x^T M y = 0,
932  if x and y are distinct eigenvectors. The matrix M is optional and is
933  assumed to be the identity if left unset.
934 
935  The efficiency of LOBPCG relies on the availability of a suitable
936  preconditioner for the matrix A. The preconditioner is supplied through the
937  SetPreconditioner() method. It should be noted that the operator used with
938  the preconditioner need not be A itself.
939 
940  For more information regarding LOBPCG see "Block Locally Optimal
941  Preconditioned Eigenvalue Xolvers (BLOPEX) in Hypre and PETSc" by
942  A. Knyazev, M. Argentati, I. Lashuk, and E. Ovtchinnikov, SISC, 29(5),
943  2224-2239, 2007.
944 */
946 {
947 private:
948  MPI_Comm comm;
949  int myid;
950  int numProcs;
951  int nev; // Number of desired eigenmodes
952  int seed; // Random seed used for initial vectors
953 
954  HYPRE_Int glbSize; // Global number of DoFs in the linear system
955  HYPRE_Int * part; // Row partitioning of the linear system
956 
957  // Pointer to HYPRE's solver struct
958  HYPRE_Solver lobpcg_solver;
959 
960  // Interface for matrix storage type
961  mv_InterfaceInterpreter interpreter;
962 
963  // Interface for setting up and performing matrix-vector products
964  HYPRE_MatvecFunctions matvec_fn;
965 
966  // Eigenvalues
967  Array<double> eigenvalues;
968 
969  // Forward declaration
970  class HypreMultiVector;
971 
972  // MultiVector to store eigenvectors
973  HypreMultiVector * multi_vec;
974 
975  // Empty vectors used to setup the matrices and preconditioner
976  HypreParVector * x;
977 
978  // An optional operator which projects vectors into a desired subspace
979  Operator * subSpaceProj;
980 
981  /// Internal class to represent a set of eigenvectors
982  class HypreMultiVector
983  {
984  private:
985  // Pointer to hypre's multi-vector object
986  mv_MultiVectorPtr mv_ptr;
987 
988  // Wrappers for each member of the multivector
989  HypreParVector ** hpv;
990 
991  // Number of vectors in the multivector
992  int nv;
993 
994  public:
995  HypreMultiVector(int n, HypreParVector & v,
996  mv_InterfaceInterpreter & interpreter);
997  ~HypreMultiVector();
998 
999  /// Set random values
1000  void Randomize(HYPRE_Int seed);
1001 
1002  /// Extract a single HypreParVector object
1003  HypreParVector & GetVector(unsigned int i);
1004 
1005  /// Transfers ownership of data to returned array of vectors
1006  HypreParVector ** StealVectors();
1007 
1008  operator mv_MultiVectorPtr() const { return mv_ptr; }
1009 
1010  mv_MultiVectorPtr & GetMultiVector() { return mv_ptr; }
1011  };
1012 
1013  static void * OperatorMatvecCreate( void *A, void *x );
1014  static HYPRE_Int OperatorMatvec( void *matvec_data,
1015  HYPRE_Complex alpha,
1016  void *A,
1017  void *x,
1018  HYPRE_Complex beta,
1019  void *y );
1020  static HYPRE_Int OperatorMatvecDestroy( void *matvec_data );
1021 
1022  static HYPRE_Int PrecondSolve(void *solver,
1023  void *A,
1024  void *b,
1025  void *x);
1026  static HYPRE_Int PrecondSetup(void *solver,
1027  void *A,
1028  void *b,
1029  void *x);
1030 
1031 public:
1032  HypreLOBPCG(MPI_Comm comm);
1033  ~HypreLOBPCG();
1034 
1035  void SetTol(double tol);
1036  void SetRelTol(double rel_tol);
1037  void SetMaxIter(int max_iter);
1038  void SetPrintLevel(int logging);
1039  void SetNumModes(int num_eigs) { nev = num_eigs; }
1040  void SetPrecondUsageMode(int pcg_mode);
1041  void SetRandomSeed(int s) { seed = s; }
1042  void SetInitialVectors(int num_vecs, HypreParVector ** vecs);
1043 
1044  // The following four methods support general operators
1045  void SetPreconditioner(Solver & precond);
1046  void SetOperator(Operator & A);
1047  void SetMassMatrix(Operator & M);
1048  void SetSubSpaceProjector(Operator & proj) { subSpaceProj = &proj; }
1049 
1050  /// Solve the eigenproblem
1051  void Solve();
1052 
1053  /// Collect the converged eigenvalues
1054  void GetEigenvalues(Array<double> & eigenvalues);
1055 
1056  /// Extract a single eigenvector
1057  HypreParVector & GetEigenvector(unsigned int i);
1058 
1059  /// Transfer ownership of the converged eigenvectors
1060  HypreParVector ** StealEigenvectors() { return multi_vec->StealVectors(); }
1061 };
1062 
1063 /** AME eigenvalue solver in hypre
1064 
1065  The Auxiliary space Maxwell Eigensolver (AME) is designed to find
1066  the lowest eigenmodes of the generalized eigenvalue problem:
1067  Curl Curl x = lambda M x
1068  where the Curl Curl operator is discretized using Nedelec finite element
1069  basis functions. Properties of this discretization are essential to
1070  eliminating the large null space of the Curl Curl operator.
1071 
1072  This eigensolver relies upon the LOBPCG eigensolver internally. It is also
1073  expected that the preconditioner supplied to this method will be the
1074  HypreAMS preconditioner defined above.
1075 
1076  As with LOBPCG, the operator set in the preconditioner need not be the same
1077  as A. This flexibility may be useful in solving eigenproblems which bare a
1078  strong resemblance to the Curl Curl problems for which AME is designed.
1079 
1080  Unlike LOBPCG, this eigensolver requires that the mass matrix be set.
1081  It is possible to circumvent this by passing an identity operator as the
1082  mass matrix but it seems unlikely that this would be useful so it is not the
1083  default behavior.
1084 */
1086 {
1087 private:
1088  int myid;
1089  int numProcs;
1090  int nev; // Number of desired eigenmodes
1091  bool setT;
1092 
1093  // Pointer to HYPRE's AME solver struct
1094  HYPRE_Solver ame_solver;
1095 
1096  // Pointer to HYPRE's AMS solver struct
1097  HypreSolver * ams_precond;
1098 
1099  // Eigenvalues
1100  HYPRE_Real * eigenvalues;
1101 
1102  // MultiVector to store eigenvectors
1103  HYPRE_ParVector * multi_vec;
1104 
1105  HypreParVector ** eigenvectors;
1106 
1107  void createDummyVectors();
1108 
1109 public:
1110  HypreAME(MPI_Comm comm);
1111  ~HypreAME();
1112 
1113  void SetTol(double tol);
1114  void SetRelTol(double rel_tol);
1115  void SetMaxIter(int max_iter);
1116  void SetPrintLevel(int logging);
1117  void SetNumModes(int num_eigs);
1118 
1119  // The following four methods support operators of type HypreParMatrix.
1120  void SetPreconditioner(HypreSolver & precond);
1121  void SetOperator(HypreParMatrix & A);
1122  void SetMassMatrix(HypreParMatrix & M);
1123 
1124  /// Solve the eigenproblem
1125  void Solve();
1126 
1127  /// Collect the converged eigenvalues
1128  void GetEigenvalues(Array<double> & eigenvalues);
1129 
1130  /// Extract a single eigenvector
1131  HypreParVector & GetEigenvector(unsigned int i);
1132 
1133  /// Transfer ownership of the converged eigenvectors
1135 };
1136 
1137 }
1138 
1139 #endif // MFEM_USE_MPI
1140 
1141 #endif
virtual HYPRE_PtrToParSolverFcn SetupFcn() const
hypre&#39;s internal Setup function
Definition: hypre.hpp:882
void GetOffd(SparseMatrix &offd, HYPRE_Int *&cmap) const
Get the local off-diagonal block. NOTE: &#39;offd&#39; will not own any data.
Definition: hypre.cpp:959
virtual ~HypreBoomerAMG()
Definition: hypre.cpp:2702
void SetTol(double tol)
Definition: hypre.cpp:2170
void SetPreconditioner(HypreSolver &precond)
Set the hypre solver to be used as a preconditioner.
Definition: hypre.cpp:2333
virtual void Mult(const HypreParVector &b, HypreParVector &x) const
Solve Ax=b with hypre&#39;s GMRES.
Definition: hypre.cpp:2341
void EliminateRowsCols(const Array< int > &rows_cols, const HypreParVector &X, HypreParVector &B)
Definition: hypre.cpp:1368
The Auxiliary-space Maxwell Solver in hypre.
Definition: hypre.hpp:859
HypreParVector * X0
FIR Filter Temporary Vectors.
Definition: hypre.hpp:527
HypreParVector & GetEigenvector(unsigned int i)
Extract a single eigenvector.
Definition: hypre.cpp:3371
double min_eig_est
Minimal eigenvalue estimate for polynomial smoothing.
Definition: hypre.hpp:557
virtual HYPRE_PtrToParSolverFcn SolveFcn() const
hypre&#39;s internal Solve function
Definition: hypre.hpp:845
HYPRE_Int GetGlobalNumCols() const
Definition: hypre.hpp:382
The Auxiliary-space Divergence Solver in hypre.
Definition: hypre.hpp:891
void SetPrintLevel(int print_lvl)
Definition: hypre.cpp:3152
void Print(const char *fname, HYPRE_Int offi=0, HYPRE_Int offj=0)
Prints the locally owned rows in parallel.
Definition: hypre.cpp:1391
MPI_Comm GetComm()
MPI communicator.
Definition: hypre.hpp:106
HypreParMatrix * RAP(const HypreParMatrix *A, const HypreParMatrix *P)
Returns the matrix P^t * A * P.
Definition: hypre.cpp:1558
int setup_called
Was hypre&#39;s Setup function called already?
Definition: hypre.hpp:630
virtual void Mult(const Vector &x, Vector &y) const
Operator application: y=A(x).
Definition: hypre.hpp:402
HYPRE_Int * ColPart()
Returns the column partitioning.
Definition: hypre.hpp:339
HYPRE_Int * GetRowStarts() const
Definition: hypre.hpp:385
void MakeRef(const HypreParMatrix &master)
Make this HypreParMatrix a reference to &#39;master&#39;.
Definition: hypre.cpp:821
void Read_IJMatrix(MPI_Comm comm, const char *fname)
Read a matrix saved as a HYPRE_IJMatrix.
Definition: hypre.cpp:1411
virtual HYPRE_PtrToParSolverFcn SolveFcn() const
GMRES Solve function.
Definition: hypre.hpp:732
HYPRE_Int MultTranspose(HypreParVector &x, HypreParVector &y, double alpha=1.0, double beta=0.0)
Computes y = alpha * A^t * x + beta * y.
Definition: hypre.cpp:1075
HypreParMatrix & Add(const double beta, const HypreParMatrix &B)
Definition: hypre.hpp:433
HypreParVector * B
Right-hand side and solution vector.
Definition: hypre.hpp:627
HypreDiagScale(HypreParMatrix &A)
Definition: hypre.hpp:761
double window_params[3]
Parameters for windowing function of FIR filter.
Definition: hypre.hpp:559
virtual void Mult(const HypreParVector &b, HypreParVector &x) const
Solve Ax=b with hypre&#39;s PCG.
Definition: hypre.cpp:2211
virtual HYPRE_PtrToParSolverFcn SetupFcn() const
PCG Setup function.
Definition: hypre.hpp:691
void SetInitialVectors(int num_vecs, HypreParVector **vecs)
Definition: hypre.cpp:3377
HypreADS(HypreParMatrix &A, ParFiniteElementSpace *face_fespace)
Definition: hypre.cpp:2912
void SetWindowByName(const char *window_name)
Convenience function for setting canonical windowing parameters.
Definition: hypre.cpp:1842
void GetEigenvalues(Array< double > &eigenvalues)
Collect the converged eigenvalues.
Definition: hypre.cpp:3638
void SetPreconditioner(HypreSolver &precond)
Definition: hypre.cpp:3604
virtual N_Vector ToNVector()
Return a new wrapper SUNDIALS N_Vector of type SUNDIALS_NVEC_PARHYP.
Definition: hypre.hpp:156
void SetTol(double tol)
Definition: hypre.cpp:3573
int GetNumCols() const
Returns the number of columns in the diagonal block of the ParCSRMatrix.
Definition: hypre.hpp:373
HYPRE_Int M() const
Returns the global number of rows.
Definition: hypre.hpp:345
void BooleanMultTranspose(int alpha, int *x, int beta, int *y)
Definition: hypre.hpp:416
virtual ~HypreGMRES()
Definition: hypre.cpp:2409
void SetPreconditioner(Solver &precond)
Definition: hypre.cpp:3294
void SetMassMatrix(Operator &M)
Definition: hypre.cpp:3349
Abstract parallel finite element space.
Definition: pfespace.hpp:28
void SetPrintLevel(int logging)
Definition: hypre.cpp:3279
virtual HYPRE_PtrToParSolverFcn SetupFcn() const
hypre&#39;s internal Setup function
Definition: hypre.hpp:915
bool iterative_mode
If true, use the second argument of Mult() as an initial guess.
Definition: operator.hpp:272
virtual ~HypreAMS()
Definition: hypre.cpp:2892
void SetTol(double tol)
Definition: hypre.cpp:2308
virtual HYPRE_PtrToParSolverFcn SolveFcn() const
hypre&#39;s internal Solve function
Definition: hypre.hpp:917
void SetOperator(HypreParMatrix &A)
Definition: hypre.cpp:3610
HypreParVector ** StealEigenvectors()
Transfer ownership of the converged eigenvectors.
Definition: hypre.hpp:1060
virtual void Mult(const HypreParVector &b, HypreParVector &x) const
Relax the linear system Ax=b.
Definition: hypre.cpp:1981
HypreParMatrix & operator+=(const HypreParMatrix &B)
Definition: hypre.hpp:428
bool pos_l1_norms
If set, take absolute values of the computed l1_norms.
Definition: hypre.hpp:553
virtual HYPRE_PtrToParSolverFcn SetupFcn() const
hypre&#39;s internal Setup function
Definition: hypre.hpp:787
virtual void Mult(const HypreParVector &b, HypreParVector &x) const
Solve the linear system Ax=b.
Definition: hypre.cpp:2105
HYPRE_Int Mult(HypreParVector &x, HypreParVector &y, double alpha=1.0, double beta=0.0)
Computes y = alpha * A * x + beta * y.
Definition: hypre.cpp:1005
void SetPrintLevel(int print_lvl)
Definition: hypre.cpp:2328
virtual HYPRE_PtrToParSolverFcn SolveFcn() const
hypre&#39;s internal Solve function
Definition: hypre.hpp:766
void ScaleRows(const Vector &s)
Scale the local row i by s(i).
Definition: hypre.cpp:1178
double poly_fraction
Fraction of spectrum to smooth for polynomial relaxation.
Definition: hypre.hpp:541
HypreParMatrix * ParMult(const HypreParMatrix *A, const HypreParMatrix *B)
Returns the matrix A * B.
Definition: hypre.cpp:1538
int GetOwnership() const
Gets ownership of the internal hypre_ParVector.
Definition: hypre.hpp:127
void SetPrintLevel(int print_lvl)
Definition: hypre.cpp:2185
const HYPRE_Int * ColPart() const
Returns the column partitioning (const version)
Definition: hypre.hpp:343
HypreParMatrix * DiscreteGrad(ParFiniteElementSpace *edge_fespace, ParFiniteElementSpace *vert_fespace)
Compute the discrete gradient matrix between the nodal linear and ND1 spaces.
void CopyRowStarts()
Definition: hypre.cpp:845
int poly_scale
Apply the polynomial smoother to A or D^{-1/2} A D^{-1/2}.
Definition: hypre.hpp:543
virtual HYPRE_PtrToParSolverFcn SetupFcn() const
GMRES Setup function.
Definition: hypre.hpp:729
void SetTol(double tol)
Definition: hypre.cpp:3257
HypreParVector * X1
Definition: hypre.hpp:527
The identity operator as a hypre solver.
Definition: hypre.hpp:743
void SetWindowParameters(double a, double b, double c)
Set parameters for windowing function for FIR smoother.
Definition: hypre.cpp:1857
void EliminateZeroRows()
If a row contains only zeros, set its diagonal to 1.
Definition: hypre.hpp:462
void SetPositiveDiagonal(bool pos=true)
After computing l1-norms, replace them with their absolute values.
Definition: hypre.hpp:605
void EliminateBC(HypreParMatrix &A, HypreParMatrix &Ae, const Array< int > &ess_dof_list, const Vector &X, Vector &B)
Definition: hypre.cpp:1612
void SetSystemsOptions(int dim)
Definition: hypre.cpp:2584
HypreLOBPCG(MPI_Comm comm)
Definition: hypre.cpp:3227
virtual HYPRE_PtrToParSolverFcn SolveFcn() const
PCG Solve function.
Definition: hypre.hpp:694
void Add(const DenseMatrix &A, const DenseMatrix &B, double alpha, DenseMatrix &C)
C = A + alpha*B.
Definition: densemat.cpp:2942
HYPRE_Int NNZ() const
Returns the global number of nonzeros.
Definition: hypre.hpp:335
hypre_ParVector * StealParVector()
Changes the ownership of the the vector.
Definition: hypre.hpp:121
The BoomerAMG solver in hypre.
Definition: hypre.hpp:796
void SetSymmetry(int sym)
Definition: hypre.cpp:2438
int dim
Definition: ex3.cpp:47
virtual HYPRE_PtrToParSolverFcn SolveFcn() const
hypre&#39;s internal Solve function
Definition: hypre.hpp:750
virtual ~HypreParaSails()
Definition: hypre.cpp:2443
The ParaSails preconditioner in hypre.
Definition: hypre.hpp:774
virtual HYPRE_PtrToParSolverFcn SetupFcn() const
hypre&#39;s internal Setup function
Definition: hypre.hpp:843
Type GetType() const
Definition: hypre.hpp:487
void SetMaxIter(int max_iter)
Definition: hypre.cpp:3273
double * l1_norms
l1 norms of the rows of A
Definition: hypre.hpp:551
Data type sparse matrix.
Definition: sparsemat.hpp:38
void SetOwnerFlags(signed char diag, signed char offd, signed char colmap)
Explicitly set the three ownership flags, see docs for diagOwner etc.
Definition: hypre.hpp:315
void SetLogging(int logging)
Definition: hypre.cpp:2180
void SetZeroInintialIterate()
non-hypre setting
Definition: hypre.hpp:678
Jacobi preconditioner in hypre.
Definition: hypre.hpp:757
virtual ~HypreSolver()
Definition: hypre.cpp:2152
HyprePCG(HypreParMatrix &_A)
Definition: hypre.cpp:2159
virtual HYPRE_PtrToParSolverFcn SolveFcn() const =0
hypre&#39;s internal Solve function
void SetRelTol(double rel_tol)
Definition: hypre.cpp:3579
void SetKDim(int dim)
Definition: hypre.cpp:2318
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:146
void SetResidualConvergenceOptions(int res_frequency=-1, double rtol=0.0)
Definition: hypre.cpp:2198
int to_int(const std::string &str)
Definition: text.hpp:70
void SetPrintLevel(int print_level)
Definition: hypre.hpp:837
virtual HYPRE_PtrToParSolverFcn SetupFcn() const
hypre&#39;s internal Setup function
Definition: hypre.hpp:764
HypreParVector ** StealEigenvectors()
Transfer ownership of the converged eigenvectors.
Definition: hypre.cpp:3685
HypreParMatrix * DiscreteCurl(ParFiniteElementSpace *face_fespace, ParFiniteElementSpace *edge_fespace)
Compute the discrete curl matrix between the ND1 and RT0 spaces.
HYPRE_Int Randomize(HYPRE_Int seed)
Set random values.
Definition: hypre.cpp:213
double relax_weight
Damping coefficient (usually <= 1)
Definition: hypre.hpp:535
Parallel smoothers in hypre.
Definition: hypre.hpp:517
void SetElasticityOptions(ParFiniteElementSpace *fespace)
Definition: hypre.cpp:2666
void GetDiag(Vector &diag) const
Get the local diagonal of the matrix.
Definition: hypre.cpp:923
HypreParMatrix * LeftDiagMult(const SparseMatrix &D, HYPRE_Int *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:1081
void BooleanMult(int alpha, int *x, int beta, int *y)
Definition: hypre.hpp:409
HypreParMatrix * A
The linear system matrix.
Definition: hypre.hpp:521
HypreAMS(HypreParMatrix &A, ParFiniteElementSpace *edge_fespace)
Definition: hypre.cpp:2713
HYPRE_Int N() const
Returns the global number of columns.
Definition: hypre.hpp:347
Dynamic 2D array using row-major layout.
Definition: array.hpp:289
virtual void SetOperator(const Operator &op)
Definition: hypre.cpp:1864
virtual HYPRE_PtrToParSolverFcn SetupFcn() const =0
hypre&#39;s internal Setup function
void SetLogging(int logging)
Definition: hypre.cpp:2323
virtual ~HypreADS()
Definition: hypre.cpp:3130
HypreParMatrix()
An empty matrix to be used as a reference to an existing matrix.
Definition: hypre.cpp:304
void SetMaxIter(int max_iter)
Definition: hypre.cpp:2175
HypreParMatrix * GetData()
Definition: hypre.hpp:769
HYPRE_Int * RowPart()
Returns the row partitioning.
Definition: hypre.hpp:337
virtual HYPRE_PtrToParSolverFcn SolveFcn() const
hypre&#39;s internal Solve function
Definition: hypre.hpp:884
void SetSOROptions(double relax_weight, double omega)
Set SOR-related parameters.
Definition: hypre.cpp:1822
MPI_Comm GetComm() const
MPI communicator.
Definition: hypre.hpp:303
Wrapper for hypre&#39;s parallel vector class.
Definition: hypre.hpp:73
HypreParVector * X
Definition: hypre.hpp:627
Type
Enumeration defining IDs for some classes derived from Operator.
Definition: operator.hpp:124
void GetEigenvalues(Array< double > &eigenvalues)
Collect the converged eigenvalues.
Definition: hypre.cpp:3359
PCG solver in hypre.
Definition: hypre.hpp:656
GMRES solver in hypre.
Definition: hypre.hpp:705
void SetPrintLevel(int print_lvl)
Definition: hypre.cpp:2907
void PrintCommPkg(std::ostream &out=mfem::out) const
Print information about the hypre_ParCSRCommPkg of the HypreParMatrix.
Definition: hypre.cpp:1432
HypreParVector * X
Definition: hypre.hpp:523
signed char OwnsColMap() const
Get colmap ownership flag.
Definition: hypre.hpp:323
void Solve()
Solve the eigenproblem.
Definition: hypre.cpp:3632
void SetNumModes(int num_eigs)
Definition: hypre.cpp:3565
HypreParMatrix(hypre_ParCSRMatrix *a, bool owner=true)
Converts hypre&#39;s format to HypreParMatrix.
Definition: hypre.hpp:228
virtual ~HypreDiagScale()
Definition: hypre.hpp:770
virtual void SetOperator(const Operator &op)
Set/update the solver for the given operator.
Definition: hypre.hpp:645
HypreParMatrix * ParAdd(const HypreParMatrix *A, const HypreParMatrix *B)
Returns the matrix A + B.
Definition: hypre.cpp:1549
double max_eig_est
Maximal eigenvalue estimate for polynomial smoothing.
Definition: hypre.hpp:555
void SetRelTol(double rel_tol)
Definition: hypre.cpp:3263
HypreAME(MPI_Comm comm)
Definition: hypre.cpp:3523
HypreParMatrix * Transpose() const
Returns the transpose of *this.
Definition: hypre.cpp:987
virtual ~HypreIdentity()
Definition: hypre.hpp:753
void GetBlocks(Array2D< HypreParMatrix *> &blocks, bool interleaved_rows=false, bool interleaved_cols=false) const
Definition: hypre.cpp:965
int height
Dimension of the output / number of rows in the matrix.
Definition: operator.hpp:24
void SetMaxIter(int max_iter)
Definition: hypre.cpp:3589
void SetFIRCoefficients(double max_eig)
Compute window and Chebyshev coefficients for given polynomial order.
Definition: hypre.cpp:1941
void SetRandomSeed(int s)
Definition: hypre.hpp:1041
double * fir_coeffs
Combined coefficients for windowing and Chebyshev polynomials.
Definition: hypre.hpp:562
void Threshold(double threshold=0.0)
Remove values smaller in absolute value than some threshold.
Definition: hypre.cpp:1292
double InnerProduct(HypreParVector *x, HypreParVector *y)
Definition: hypre.cpp:252
HYPRE_Int GlobalSize()
Returns the global number of rows.
Definition: hypre.hpp:112
HYPRE_Int * GetColStarts() const
Definition: hypre.hpp:387
virtual ~HypreSmoother()
Definition: hypre.cpp:2071
HypreParVector(MPI_Comm comm, HYPRE_Int glob_size, HYPRE_Int *col)
Definition: hypre.cpp:98
void SetSubSpaceProjector(Operator &proj)
Definition: hypre.hpp:1048
void Print(const char *fname) const
Prints the locally owned rows in parallel.
Definition: hypre.cpp:218
virtual void SetOperator(const Operator &op)
Set/update the solver for the given operator.
Definition: hypre.cpp:2567
void SetData(double *_data)
Definition: hypre.cpp:208
void InvScaleRows(const Vector &s)
Scale the local row i by 1./s(i)
Definition: hypre.cpp:1213
void SetTaubinOptions(double lambda, double mu, int iter)
Set parameters for Taubin&#39;s lambda-mu method.
Definition: hypre.cpp:1834
HypreParVector * B
Right-hand side and solution vectors.
Definition: hypre.hpp:523
void SetMassMatrix(HypreParMatrix &M)
Definition: hypre.cpp:3625
virtual HYPRE_PtrToParSolverFcn SolveFcn() const
hypre&#39;s internal Solve function
Definition: hypre.hpp:789
void SetZeroInintialIterate()
non-hypre setting
Definition: hypre.hpp:723
virtual ~HyprePCG()
Definition: hypre.cpp:2284
void SetOperator(Operator &A)
Definition: hypre.cpp:3303
void SetPreconditioner(HypreSolver &precond)
Set the hypre solver to be used as a preconditioner.
Definition: hypre.cpp:2190
void SetPolyOptions(int poly_order, double poly_fraction)
Set parameters for polynomial smoothing.
Definition: hypre.cpp:1828
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:404
void GetNumIterations(int &num_iterations)
Definition: hypre.hpp:680
void CopyColStarts()
Definition: hypre.cpp:882
int relax_times
Number of relaxation sweeps.
Definition: hypre.hpp:533
void operator*=(double s)
Scale all entries by s: A_scaled = s*A.
Definition: hypre.cpp:1254
void SetPrintLevel(int logging)
Definition: hypre.cpp:3595
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:263
void SetOwnership(int own)
Sets ownership of the internal hypre_ParVector.
Definition: hypre.hpp:124
Abstract class for hypre&#39;s solvers and preconditioners.
Definition: hypre.hpp:620
HypreParVector & operator=(double d)
Set constant values.
Definition: hypre.cpp:186
const double alpha
Definition: ex15.cpp:337
const HYPRE_Int * RowPart() const
Returns the row partitioning (const version)
Definition: hypre.hpp:341
Vector data type.
Definition: vector.hpp:48
signed char OwnsDiag() const
Get diag ownership flag.
Definition: hypre.hpp:319
HypreParMatrix & operator=(double value)
Initialize all entries with value.
Definition: hypre.hpp:422
HypreGMRES(HypreParMatrix &_A)
Definition: hypre.cpp:2290
virtual HYPRE_PtrToParSolverFcn SetupFcn() const
hypre&#39;s internal Setup function
Definition: hypre.hpp:748
hypre_ParCSRMatrix * StealData()
Changes the ownership of the the matrix.
Definition: hypre.cpp:831
void SetPrecondUsageMode(int pcg_mode)
Definition: hypre.cpp:3288
virtual ~HypreParMatrix()
Calls hypre&#39;s destroy function.
Definition: hypre.hpp:485
HypreParMatrix * A
The linear system matrix.
Definition: hypre.hpp:624
ID for class HypreParMatrix.
Definition: operator.hpp:128
void SetMaxIter(int max_iter)
Definition: hypre.cpp:2313
Vector * GlobalVector() const
Returns the global vector in each processor.
Definition: hypre.cpp:176
HypreParVector & GetEigenvector(unsigned int i)
Extract a single eigenvector.
Definition: hypre.cpp:3674
double omega
SOR parameter (usually in (0,2))
Definition: hypre.hpp:537
signed char OwnsOffd() const
Get offd ownership flag.
Definition: hypre.hpp:321
Base class for solvers.
Definition: operator.hpp:268
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:64
Abstract operator.
Definition: operator.hpp:21
Wrapper for hypre&#39;s ParCSR matrix class.
Definition: hypre.hpp:175
~HypreParVector()
Calls hypre&#39;s destroy function.
Definition: hypre.cpp:223
void SetNumModes(int num_eigs)
Definition: hypre.hpp:1039
HypreParVector * Z
Definition: hypre.hpp:525
void SetSingularProblem()
Set this option when solving a curl-curl problem with zero mass term.
Definition: hypre.hpp:877
void Solve()
Solve the eigenproblem.
Definition: hypre.cpp:3416
void Read(MPI_Comm comm, const char *fname)
Reads the matrix from a file.
Definition: hypre.cpp:1396
HypreParVector * V
Temporary vectors.
Definition: hypre.hpp:525
void SetType(HypreSmoother::Type type, int relax_times=1)
Set the relaxation type and number of sweeps.
Definition: hypre.cpp:1816
HypreParaSails(HypreParMatrix &A)
Definition: hypre.cpp:2415
HYPRE_Int GetGlobalNumRows() const
Definition: hypre.hpp:379
int width
Dimension of the input / number of columns in the matrix.
Definition: operator.hpp:25
int poly_order
Order of the smoothing polynomial.
Definition: hypre.hpp:539
HYPRE_Int * Partitioning()
Returns the row partitioning.
Definition: hypre.hpp:109
int GetNumRows() const
Returns the number of rows in the diagonal block of the ParCSRMatrix.
Definition: hypre.hpp:366
double lambda
Taubin&#39;s lambda-mu method parameters.
Definition: hypre.hpp:546