MFEM  v3.4
Finite element discretization library
hypre.cpp
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 #include "../config/config.hpp"
13 
14 #ifdef MFEM_USE_MPI
15 
16 #include "linalg.hpp"
17 #include "../fem/fem.hpp"
18 
19 #include <fstream>
20 #include <iomanip>
21 #include <cmath>
22 #include <cstdlib>
23 
24 // Define macro wrappers for hypre_TAlloc, hypre_CTAlloc and hypre_TFree:
25 // mfem_hypre_TAlloc, mfem_hypre_CTAlloc, and mfem_hypre_TFree, respectively.
26 // Note: the same macros are defined in hypre_parcsr.cpp.
27 #if MFEM_HYPRE_VERSION < 21400
28 
29 #define mfem_hypre_TAlloc(type, size) hypre_TAlloc(type, size)
30 #define mfem_hypre_CTAlloc(type, size) hypre_CTAlloc(type, size)
31 #define mfem_hypre_TFree(ptr) hypre_TFree(ptr)
32 
33 #else // MFEM_HYPRE_VERSION >= 21400
34 
35 #define mfem_hypre_TAlloc(type, size) \
36  hypre_TAlloc(type, size, HYPRE_MEMORY_HOST)
37 #define mfem_hypre_CTAlloc(type, size) \
38  hypre_CTAlloc(type, size, HYPRE_MEMORY_HOST)
39 #define mfem_hypre_TFree(ptr) hypre_TFree(ptr, HYPRE_MEMORY_HOST)
40 
41 // Notes regarding allocation and deallocation of hypre objects in 2.14.0
42 //-----------------------------------------------------------------------
43 //
44 // 1. hypre_CSRMatrix: i, j, data, and rownnz use HYPRE_MEMORY_SHARED while the
45 // hypre_CSRMatrix structure uses HYPRE_MEMORY_HOST.
46 //
47 // Note: the function HYPRE_CSRMatrixCreate creates the i array using
48 // HYPRE_MEMORY_HOST!
49 // Note: the functions hypre_CSRMatrixAdd and hypre_CSRMatrixMultiply create
50 // C_i using HYPRE_MEMORY_HOST!
51 //
52 // 2. hypre_Vector: data uses HYPRE_MEMORY_SHARED while the hypre_Vector
53 // structure uses HYPRE_MEMORY_HOST.
54 //
55 // 3. hypre_ParVector: the structure hypre_ParVector uses HYPRE_MEMORY_HOST;
56 // partitioning uses HYPRE_MEMORY_HOST.
57 //
58 // 4. hypre_ParCSRMatrix: the structure hypre_ParCSRMatrix uses
59 // HYPRE_MEMORY_HOST; col_map_offd, row_starts, col_starts, rowindices,
60 // rowvalues also use HYPRE_MEMORY_HOST.
61 //
62 // Note: the function hypre_ParCSRMatrixToCSRMatrixAll allocates matrix_i
63 // using HYPRE_MEMORY_HOST!
64 //
65 // 5. The goal for the MFEM wrappers of hypre objects is to support only the
66 // standard hypre build case, i.e. when hypre is build without device support
67 // and all memory types correspond to host memory. In this case memory
68 // allocated with operator new can be used by hypre but (as usual) it must
69 // not be owned by hypre.
70 
71 #endif // #if MFEM_HYPRE_VERSION < 21400
72 
73 using namespace std;
74 
75 namespace mfem
76 {
77 
78 template<typename TargetT, typename SourceT>
79 static TargetT *DuplicateAs(const SourceT *array, int size,
80  bool cplusplus = true)
81 {
82  TargetT *target_array = cplusplus ? new TargetT[size]
83  /* */ : mfem_hypre_TAlloc(TargetT, size);
84  for (int i = 0; i < size; i++)
85  {
86  target_array[i] = array[i];
87  }
88  return target_array;
89 }
90 
91 inline void HypreParVector::_SetDataAndSize_()
92 {
93  SetDataAndSize(hypre_VectorData(hypre_ParVectorLocalVector(x)),
95  hypre_VectorSize(hypre_ParVectorLocalVector(x))));
96 }
97 
98 HypreParVector::HypreParVector(MPI_Comm comm, HYPRE_Int glob_size,
99  HYPRE_Int *col) : Vector()
100 {
101  x = hypre_ParVectorCreate(comm,glob_size,col);
102  hypre_ParVectorInitialize(x);
103  hypre_ParVectorSetPartitioningOwner(x,0);
104  // The data will be destroyed by hypre (this is the default)
105  hypre_ParVectorSetDataOwner(x,1);
106  hypre_SeqVectorSetDataOwner(hypre_ParVectorLocalVector(x),1);
107  _SetDataAndSize_();
108  own_ParVector = 1;
109 }
110 
111 HypreParVector::HypreParVector(MPI_Comm comm, HYPRE_Int glob_size,
112  double *_data, HYPRE_Int *col) : Vector()
113 {
114  x = hypre_ParVectorCreate(comm,glob_size,col);
115  hypre_ParVectorSetDataOwner(x,1); // owns the seq vector
116  hypre_SeqVectorSetDataOwner(hypre_ParVectorLocalVector(x),0);
117  hypre_ParVectorSetPartitioningOwner(x,0);
118  double tmp = 0.0;
119  hypre_VectorData(hypre_ParVectorLocalVector(x)) = &tmp;
120  // If hypre_ParVectorLocalVector(x) and &tmp are non-NULL,
121  // hypre_ParVectorInitialize(x) does not allocate memory!
122  hypre_ParVectorInitialize(x);
123  // Set the internal data array to the one passed in
124  hypre_VectorData(hypre_ParVectorLocalVector(x)) = _data;
125  _SetDataAndSize_();
126  own_ParVector = 1;
127 }
128 
130 {
131  x = hypre_ParVectorCreate(y.x -> comm, y.x -> global_size,
132  y.x -> partitioning);
133  hypre_ParVectorInitialize(x);
134  hypre_ParVectorSetPartitioningOwner(x,0);
135  hypre_ParVectorSetDataOwner(x,1);
136  hypre_SeqVectorSetDataOwner(hypre_ParVectorLocalVector(x),1);
137  _SetDataAndSize_();
138  own_ParVector = 1;
139 }
140 
142  int transpose) : Vector()
143 {
144  if (!transpose)
145  {
146  x = hypre_ParVectorInDomainOf(const_cast<HypreParMatrix&>(A));
147  }
148  else
149  {
150  x = hypre_ParVectorInRangeOf(const_cast<HypreParMatrix&>(A));
151  }
152  _SetDataAndSize_();
153  own_ParVector = 1;
154 }
155 
157 {
158  x = (hypre_ParVector *) y;
159  _SetDataAndSize_();
160  own_ParVector = 0;
161 }
162 
164 {
165  x = hypre_ParVectorCreate(pfes->GetComm(), pfes->GlobalTrueVSize(),
166  pfes->GetTrueDofOffsets());
167  hypre_ParVectorInitialize(x);
168  hypre_ParVectorSetPartitioningOwner(x,0);
169  // The data will be destroyed by hypre (this is the default)
170  hypre_ParVectorSetDataOwner(x,1);
171  hypre_SeqVectorSetDataOwner(hypre_ParVectorLocalVector(x),1);
172  _SetDataAndSize_();
173  own_ParVector = 1;
174 }
175 
177 {
178  hypre_Vector *hv = hypre_ParVectorToVectorAll(*this);
179  Vector *v = new Vector(hv->data, internal::to_int(hv->size));
180  v->MakeDataOwner();
181  hypre_SeqVectorSetDataOwner(hv,0);
182  hypre_SeqVectorDestroy(hv);
183  return v;
184 }
185 
187 {
188  hypre_ParVectorSetConstantValues(x,d);
189  return *this;
190 }
191 
193 {
194 #ifdef MFEM_DEBUG
195  if (size != y.Size())
196  {
197  mfem_error("HypreParVector::operator=");
198  }
199 #endif
200 
201  for (int i = 0; i < size; i++)
202  {
203  data[i] = y.data[i];
204  }
205  return *this;
206 }
207 
208 void HypreParVector::SetData(double *_data)
209 {
210  Vector::data = hypre_VectorData(hypre_ParVectorLocalVector(x)) = _data;
211 }
212 
213 HYPRE_Int HypreParVector::Randomize(HYPRE_Int seed)
214 {
215  return hypre_ParVectorSetRandomValues(x,seed);
216 }
217 
218 void HypreParVector::Print(const char *fname) const
219 {
220  hypre_ParVectorPrint(x,fname);
221 }
222 
224 {
225  if (own_ParVector)
226  {
227  hypre_ParVectorDestroy(x);
228  }
229 }
230 
231 #ifdef MFEM_USE_SUNDIALS
232 
233 #ifndef SUNFALSE
234 #define SUNFALSE FALSE
235 #endif
236 
237 void HypreParVector::ToNVector(N_Vector &nv)
238 {
239  MFEM_ASSERT(nv && N_VGetVectorID(nv) == SUNDIALS_NVEC_PARHYP,
240  "invalid N_Vector");
241  N_VectorContent_ParHyp nv_c = (N_VectorContent_ParHyp)(nv->content);
242  MFEM_ASSERT(nv_c->own_parvector == SUNFALSE, "invalid N_Vector");
243  nv_c->local_length = x->local_vector->size;
244  nv_c->global_length = x->global_size;
245  nv_c->comm = x->comm;
246  nv_c->x = x;
247 }
248 
249 #endif // MFEM_USE_SUNDIALS
250 
251 
253 {
254  return hypre_ParVectorInnerProd(*x, *y);
255 }
256 
258 {
259  return hypre_ParVectorInnerProd(x, y);
260 }
261 
262 
263 double ParNormlp(const Vector &vec, double p, MPI_Comm comm)
264 {
265  double norm = 0.0;
266  if (p == 1.0)
267  {
268  double loc_norm = vec.Norml1();
269  MPI_Allreduce(&loc_norm, &norm, 1, MPI_DOUBLE, MPI_SUM, comm);
270  }
271  if (p == 2.0)
272  {
273  double loc_norm = vec*vec;
274  MPI_Allreduce(&loc_norm, &norm, 1, MPI_DOUBLE, MPI_SUM, comm);
275  norm = sqrt(norm);
276  }
277  if (p < infinity())
278  {
279  double sum = 0.0;
280  for (int i = 0; i < vec.Size(); i++)
281  {
282  sum += pow(fabs(vec(i)), p);
283  }
284  MPI_Allreduce(&sum, &norm, 1, MPI_DOUBLE, MPI_SUM, comm);
285  norm = pow(norm, 1.0/p);
286  }
287  else
288  {
289  double loc_norm = vec.Normlinf();
290  MPI_Allreduce(&loc_norm, &norm, 1, MPI_DOUBLE, MPI_MAX, comm);
291  }
292  return norm;
293 }
294 
295 
296 void HypreParMatrix::Init()
297 {
298  A = NULL;
299  X = Y = NULL;
300  diagOwner = offdOwner = colMapOwner = -1;
301  ParCSROwner = 1;
302 }
303 
305 {
306  Init();
307  height = width = 0;
308 }
309 
310 char HypreParMatrix::CopyCSR(SparseMatrix *csr, hypre_CSRMatrix *hypre_csr)
311 {
312  hypre_CSRMatrixData(hypre_csr) = csr->GetData();
313 #ifndef HYPRE_BIGINT
314  hypre_CSRMatrixI(hypre_csr) = csr->GetI();
315  hypre_CSRMatrixJ(hypre_csr) = csr->GetJ();
316  // Prevent hypre from destroying hypre_csr->{i,j,data}
317  return 0;
318 #else
319  hypre_CSRMatrixI(hypre_csr) =
320  DuplicateAs<HYPRE_Int>(csr->GetI(), csr->Height()+1);
321  hypre_CSRMatrixJ(hypre_csr) =
322  DuplicateAs<HYPRE_Int>(csr->GetJ(), csr->NumNonZeroElems());
323  // Prevent hypre from destroying hypre_csr->{i,j,data}, own {i,j}
324  return 1;
325 #endif
326 }
327 
328 char HypreParMatrix::CopyBoolCSR(Table *bool_csr, hypre_CSRMatrix *hypre_csr)
329 {
330  int nnz = bool_csr->Size_of_connections();
331  double *data = new double[nnz];
332  for (int i = 0; i < nnz; i++)
333  {
334  data[i] = 1.0;
335  }
336  hypre_CSRMatrixData(hypre_csr) = data;
337 #ifndef HYPRE_BIGINT
338  hypre_CSRMatrixI(hypre_csr) = bool_csr->GetI();
339  hypre_CSRMatrixJ(hypre_csr) = bool_csr->GetJ();
340  // Prevent hypre from destroying hypre_csr->{i,j,data}, own {data}
341  return 2;
342 #else
343  hypre_CSRMatrixI(hypre_csr) =
344  DuplicateAs<HYPRE_Int>(bool_csr->GetI(), bool_csr->Size()+1);
345  hypre_CSRMatrixJ(hypre_csr) =
346  DuplicateAs<HYPRE_Int>(bool_csr->GetJ(), nnz);
347  // Prevent hypre from destroying hypre_csr->{i,j,data}, own {i,j,data}
348  return 3;
349 #endif
350 }
351 
352 void HypreParMatrix::CopyCSR_J(hypre_CSRMatrix *hypre_csr, int *J)
353 {
354  HYPRE_Int nnz = hypre_CSRMatrixNumNonzeros(hypre_csr);
355  for (HYPRE_Int j = 0; j < nnz; j++)
356  {
357  J[j] = int(hypre_CSRMatrixJ(hypre_csr)[j]);
358  }
359 }
360 
361 // Square block-diagonal constructor (4 arguments, v1)
362 HypreParMatrix::HypreParMatrix(MPI_Comm comm, HYPRE_Int glob_size,
363  HYPRE_Int *row_starts, SparseMatrix *diag)
364  : Operator(diag->Height(), diag->Width())
365 {
366  Init();
367  A = hypre_ParCSRMatrixCreate(comm, glob_size, glob_size, row_starts,
368  row_starts, 0, diag->NumNonZeroElems(), 0);
369  hypre_ParCSRMatrixSetDataOwner(A,1);
370  hypre_ParCSRMatrixSetRowStartsOwner(A,0);
371  hypre_ParCSRMatrixSetColStartsOwner(A,0);
372 
373  hypre_CSRMatrixSetDataOwner(A->diag,0);
374  diagOwner = CopyCSR(diag, A->diag);
375  hypre_CSRMatrixSetRownnz(A->diag);
376 
377  hypre_CSRMatrixSetDataOwner(A->offd,1);
378  hypre_CSRMatrixI(A->offd) = mfem_hypre_CTAlloc(HYPRE_Int, diag->Height()+1);
379 
380  /* Don't need to call these, since they allocate memory only
381  if it was not already allocated */
382  // hypre_CSRMatrixInitialize(A->diag);
383  // hypre_ParCSRMatrixInitialize(A);
384 
385  hypre_ParCSRMatrixSetNumNonzeros(A);
386 
387  /* Make sure that the first entry in each row is the diagonal one. */
388  hypre_CSRMatrixReorder(hypre_ParCSRMatrixDiag(A));
389 #ifdef HYPRE_BIGINT
390  CopyCSR_J(A->diag, diag->GetJ());
391 #endif
392 
393  hypre_MatvecCommPkgCreate(A);
394 }
395 
396 // Rectangular block-diagonal constructor (6 arguments, v1)
398  HYPRE_Int global_num_rows,
399  HYPRE_Int global_num_cols,
400  HYPRE_Int *row_starts, HYPRE_Int *col_starts,
401  SparseMatrix *diag)
402  : Operator(diag->Height(), diag->Width())
403 {
404  Init();
405  A = hypre_ParCSRMatrixCreate(comm, global_num_rows, global_num_cols,
406  row_starts, col_starts,
407  0, diag->NumNonZeroElems(), 0);
408  hypre_ParCSRMatrixSetDataOwner(A,1);
409  hypre_ParCSRMatrixSetRowStartsOwner(A,0);
410  hypre_ParCSRMatrixSetColStartsOwner(A,0);
411 
412  hypre_CSRMatrixSetDataOwner(A->diag,0);
413  diagOwner = CopyCSR(diag, A->diag);
414  hypre_CSRMatrixSetRownnz(A->diag);
415 
416  hypre_CSRMatrixSetDataOwner(A->offd,1);
417  hypre_CSRMatrixI(A->offd) = mfem_hypre_CTAlloc(HYPRE_Int, diag->Height()+1);
418 
419  hypre_ParCSRMatrixSetNumNonzeros(A);
420 
421  /* Make sure that the first entry in each row is the diagonal one. */
422  if (row_starts == col_starts)
423  {
424  hypre_CSRMatrixReorder(hypre_ParCSRMatrixDiag(A));
425 #ifdef HYPRE_BIGINT
426  CopyCSR_J(A->diag, diag->GetJ());
427 #endif
428  }
429 
430  hypre_MatvecCommPkgCreate(A);
431 }
432 
433 // General rectangular constructor with diagonal and off-diagonal (8 arguments)
435  HYPRE_Int global_num_rows,
436  HYPRE_Int global_num_cols,
437  HYPRE_Int *row_starts, HYPRE_Int *col_starts,
438  SparseMatrix *diag, SparseMatrix *offd,
439  HYPRE_Int *cmap)
440  : Operator(diag->Height(), diag->Width())
441 {
442  Init();
443  A = hypre_ParCSRMatrixCreate(comm, global_num_rows, global_num_cols,
444  row_starts, col_starts,
445  offd->Width(), diag->NumNonZeroElems(),
446  offd->NumNonZeroElems());
447  hypre_ParCSRMatrixSetDataOwner(A,1);
448  hypre_ParCSRMatrixSetRowStartsOwner(A,0);
449  hypre_ParCSRMatrixSetColStartsOwner(A,0);
450 
451  hypre_CSRMatrixSetDataOwner(A->diag,0);
452  diagOwner = CopyCSR(diag, A->diag);
453  hypre_CSRMatrixSetRownnz(A->diag);
454 
455  hypre_CSRMatrixSetDataOwner(A->offd,0);
456  offdOwner = CopyCSR(offd, A->offd);
457  hypre_CSRMatrixSetRownnz(A->offd);
458 
459  hypre_ParCSRMatrixColMapOffd(A) = cmap;
460  // Prevent hypre from destroying A->col_map_offd
461  colMapOwner = 0;
462 
463  hypre_ParCSRMatrixSetNumNonzeros(A);
464 
465  /* Make sure that the first entry in each row is the diagonal one. */
466  if (row_starts == col_starts)
467  {
468  hypre_CSRMatrixReorder(hypre_ParCSRMatrixDiag(A));
469 #ifdef HYPRE_BIGINT
470  CopyCSR_J(A->diag, diag->GetJ());
471 #endif
472  }
473 
474  hypre_MatvecCommPkgCreate(A);
475 }
476 
477 // General rectangular constructor with diagonal and off-diagonal (13 arguments)
479  MPI_Comm comm,
480  HYPRE_Int global_num_rows, HYPRE_Int global_num_cols,
481  HYPRE_Int *row_starts, HYPRE_Int *col_starts,
482  HYPRE_Int *diag_i, HYPRE_Int *diag_j, double *diag_data,
483  HYPRE_Int *offd_i, HYPRE_Int *offd_j, double *offd_data,
484  HYPRE_Int offd_num_cols, HYPRE_Int *offd_col_map)
485 {
486  Init();
487  A = hypre_ParCSRMatrixCreate(comm, global_num_rows, global_num_cols,
488  row_starts, col_starts, offd_num_cols, 0, 0);
489  hypre_ParCSRMatrixSetDataOwner(A,1);
490  hypre_ParCSRMatrixSetRowStartsOwner(A,0);
491  hypre_ParCSRMatrixSetColStartsOwner(A,0);
492 
493  HYPRE_Int local_num_rows = hypre_CSRMatrixNumRows(A->diag);
494 
495  hypre_CSRMatrixSetDataOwner(A->diag,0);
496  hypre_CSRMatrixI(A->diag) = diag_i;
497  hypre_CSRMatrixJ(A->diag) = diag_j;
498  hypre_CSRMatrixData(A->diag) = diag_data;
499  hypre_CSRMatrixNumNonzeros(A->diag) = diag_i[local_num_rows];
500  hypre_CSRMatrixSetRownnz(A->diag);
501  // Prevent hypre from destroying A->diag->{i,j,data}, own A->diag->{i,j,data}
502  diagOwner = 3;
503 
504  hypre_CSRMatrixSetDataOwner(A->offd,0);
505  hypre_CSRMatrixI(A->offd) = offd_i;
506  hypre_CSRMatrixJ(A->offd) = offd_j;
507  hypre_CSRMatrixData(A->offd) = offd_data;
508  hypre_CSRMatrixNumNonzeros(A->offd) = offd_i[local_num_rows];
509  hypre_CSRMatrixSetRownnz(A->offd);
510  // Prevent hypre from destroying A->offd->{i,j,data}, own A->offd->{i,j,data}
511  offdOwner = 3;
512 
513  hypre_ParCSRMatrixColMapOffd(A) = offd_col_map;
514  // Prevent hypre from destroying A->col_map_offd, own A->col_map_offd
515  colMapOwner = 1;
516 
517  hypre_ParCSRMatrixSetNumNonzeros(A);
518 
519  /* Make sure that the first entry in each row is the diagonal one. */
520  if (row_starts == col_starts)
521  {
522  hypre_CSRMatrixReorder(hypre_ParCSRMatrixDiag(A));
523  }
524 
525  hypre_MatvecCommPkgCreate(A);
526 
527  height = GetNumRows();
528  width = GetNumCols();
529 }
530 
531 // Constructor from a CSR matrix on rank 0 (4 arguments, v2)
533  HYPRE_Int *row_starts, HYPRE_Int *col_starts,
534  SparseMatrix *sm_a)
535 {
536  MFEM_ASSERT(sm_a != NULL, "invalid input");
537  MFEM_VERIFY(!HYPRE_AssumedPartitionCheck(),
538  "this method can not be used with assumed partition");
539 
540  Init();
541 
542  hypre_CSRMatrix *csr_a;
543  csr_a = hypre_CSRMatrixCreate(sm_a -> Height(), sm_a -> Width(),
544  sm_a -> NumNonZeroElems());
545 
546  hypre_CSRMatrixSetDataOwner(csr_a,0);
547  CopyCSR(sm_a, csr_a);
548  hypre_CSRMatrixSetRownnz(csr_a);
549 
550  A = hypre_CSRMatrixToParCSRMatrix(comm, csr_a, row_starts, col_starts);
551 
552 #ifdef HYPRE_BIGINT
553  delete [] hypre_CSRMatrixI(csr_a);
554  delete [] hypre_CSRMatrixJ(csr_a);
555 #endif
556  hypre_CSRMatrixI(csr_a) = NULL;
557  hypre_CSRMatrixDestroy(csr_a);
558 
559  height = GetNumRows();
560  width = GetNumCols();
561 
562  /* Make sure that the first entry in each row is the diagonal one. */
563  if (row_starts == col_starts)
564  {
565  hypre_CSRMatrixReorder(hypre_ParCSRMatrixDiag(A));
566  }
567 
568  hypre_MatvecCommPkgCreate(A);
569 }
570 
571 // Boolean, rectangular, block-diagonal constructor (6 arguments, v2)
573  HYPRE_Int global_num_rows,
574  HYPRE_Int global_num_cols,
575  HYPRE_Int *row_starts, HYPRE_Int *col_starts,
576  Table *diag)
577 {
578  Init();
579  int nnz = diag->Size_of_connections();
580  A = hypre_ParCSRMatrixCreate(comm, global_num_rows, global_num_cols,
581  row_starts, col_starts, 0, nnz, 0);
582  hypre_ParCSRMatrixSetDataOwner(A,1);
583  hypre_ParCSRMatrixSetRowStartsOwner(A,0);
584  hypre_ParCSRMatrixSetColStartsOwner(A,0);
585 
586  hypre_CSRMatrixSetDataOwner(A->diag,0);
587  diagOwner = CopyBoolCSR(diag, A->diag);
588  hypre_CSRMatrixSetRownnz(A->diag);
589 
590  hypre_CSRMatrixSetDataOwner(A->offd,1);
591  hypre_CSRMatrixI(A->offd) = mfem_hypre_CTAlloc(HYPRE_Int, diag->Size()+1);
592 
593  hypre_ParCSRMatrixSetNumNonzeros(A);
594 
595  /* Make sure that the first entry in each row is the diagonal one. */
596  if (row_starts == col_starts)
597  {
598  hypre_CSRMatrixReorder(hypre_ParCSRMatrixDiag(A));
599 #ifdef HYPRE_BIGINT
600  CopyCSR_J(A->diag, diag->GetJ());
601 #endif
602  }
603 
604  hypre_MatvecCommPkgCreate(A);
605 
606  height = GetNumRows();
607  width = GetNumCols();
608 }
609 
610 // Boolean, general rectangular constructor with diagonal and off-diagonal
611 // (11 arguments)
612 HypreParMatrix::HypreParMatrix(MPI_Comm comm, int id, int np,
613  HYPRE_Int *row, HYPRE_Int *col,
614  HYPRE_Int *i_diag, HYPRE_Int *j_diag,
615  HYPRE_Int *i_offd, HYPRE_Int *j_offd,
616  HYPRE_Int *cmap, HYPRE_Int cmap_size)
617 {
618  HYPRE_Int diag_nnz, offd_nnz;
619 
620  Init();
621  if (HYPRE_AssumedPartitionCheck())
622  {
623  diag_nnz = i_diag[row[1]-row[0]];
624  offd_nnz = i_offd[row[1]-row[0]];
625 
626  A = hypre_ParCSRMatrixCreate(comm, row[2], col[2], row, col,
627  cmap_size, diag_nnz, offd_nnz);
628  }
629  else
630  {
631  diag_nnz = i_diag[row[id+1]-row[id]];
632  offd_nnz = i_offd[row[id+1]-row[id]];
633 
634  A = hypre_ParCSRMatrixCreate(comm, row[np], col[np], row, col,
635  cmap_size, diag_nnz, offd_nnz);
636  }
637 
638  hypre_ParCSRMatrixSetDataOwner(A,1);
639  hypre_ParCSRMatrixSetRowStartsOwner(A,0);
640  hypre_ParCSRMatrixSetColStartsOwner(A,0);
641 
642  HYPRE_Int i;
643 
644  double *a_diag = new double[diag_nnz];
645  for (i = 0; i < diag_nnz; i++)
646  {
647  a_diag[i] = 1.0;
648  }
649 
650  double *a_offd = new double[offd_nnz];
651  for (i = 0; i < offd_nnz; i++)
652  {
653  a_offd[i] = 1.0;
654  }
655 
656  hypre_CSRMatrixSetDataOwner(A->diag,0);
657  hypre_CSRMatrixI(A->diag) = i_diag;
658  hypre_CSRMatrixJ(A->diag) = j_diag;
659  hypre_CSRMatrixData(A->diag) = a_diag;
660  hypre_CSRMatrixSetRownnz(A->diag);
661  // Prevent hypre from destroying A->diag->{i,j,data}, own A->diag->{i,j,data}
662  diagOwner = 3;
663 
664  hypre_CSRMatrixSetDataOwner(A->offd,0);
665  hypre_CSRMatrixI(A->offd) = i_offd;
666  hypre_CSRMatrixJ(A->offd) = j_offd;
667  hypre_CSRMatrixData(A->offd) = a_offd;
668  hypre_CSRMatrixSetRownnz(A->offd);
669  // Prevent hypre from destroying A->offd->{i,j,data}, own A->offd->{i,j,data}
670  offdOwner = 3;
671 
672  hypre_ParCSRMatrixColMapOffd(A) = cmap;
673  // Prevent hypre from destroying A->col_map_offd, own A->col_map_offd
674  colMapOwner = 1;
675 
676  hypre_ParCSRMatrixSetNumNonzeros(A);
677 
678  /* Make sure that the first entry in each row is the diagonal one. */
679  if (row == col)
680  {
681  hypre_CSRMatrixReorder(hypre_ParCSRMatrixDiag(A));
682  }
683 
684  hypre_MatvecCommPkgCreate(A);
685 
686  height = GetNumRows();
687  width = GetNumCols();
688 }
689 
690 // General rectangular constructor with diagonal and off-diagonal constructed
691 // from a CSR matrix that contains both diagonal and off-diagonal blocks
692 // (9 arguments)
693 HypreParMatrix::HypreParMatrix(MPI_Comm comm, int nrows, HYPRE_Int glob_nrows,
694  HYPRE_Int glob_ncols, int *I, HYPRE_Int *J,
695  double *data, HYPRE_Int *rows, HYPRE_Int *cols)
696 {
697  Init();
698 
699  // Determine partitioning size, and my column start and end
700  int part_size;
701  HYPRE_Int my_col_start, my_col_end; // my range: [my_col_start, my_col_end)
702  if (HYPRE_AssumedPartitionCheck())
703  {
704  part_size = 2;
705  my_col_start = cols[0];
706  my_col_end = cols[1];
707  }
708  else
709  {
710  int myid;
711  MPI_Comm_rank(comm, &myid);
712  MPI_Comm_size(comm, &part_size);
713  part_size++;
714  my_col_start = cols[myid];
715  my_col_end = cols[myid+1];
716  }
717 
718  // Copy in the row and column partitionings
719  HYPRE_Int *row_starts, *col_starts;
720  if (rows == cols)
721  {
722  row_starts = col_starts = mfem_hypre_TAlloc(HYPRE_Int, part_size);
723  for (int i = 0; i < part_size; i++)
724  {
725  row_starts[i] = rows[i];
726  }
727  }
728  else
729  {
730  row_starts = mfem_hypre_TAlloc(HYPRE_Int, part_size);
731  col_starts = mfem_hypre_TAlloc(HYPRE_Int, part_size);
732  for (int i = 0; i < part_size; i++)
733  {
734  row_starts[i] = rows[i];
735  col_starts[i] = cols[i];
736  }
737  }
738 
739  // Create a map for the off-diagonal indices - global to local. Count the
740  // number of diagonal and off-diagonal entries.
741  HYPRE_Int diag_nnz = 0, offd_nnz = 0, offd_num_cols = 0;
742  map<HYPRE_Int, HYPRE_Int> offd_map;
743  for (HYPRE_Int j = 0, loc_nnz = I[nrows]; j < loc_nnz; j++)
744  {
745  HYPRE_Int glob_col = J[j];
746  if (my_col_start <= glob_col && glob_col < my_col_end)
747  {
748  diag_nnz++;
749  }
750  else
751  {
752  offd_map.insert(pair<const HYPRE_Int, HYPRE_Int>(glob_col, -1));
753  offd_nnz++;
754  }
755  }
756  // count the number of columns in the off-diagonal and set the local indices
757  for (map<HYPRE_Int, HYPRE_Int>::iterator it = offd_map.begin();
758  it != offd_map.end(); ++it)
759  {
760  it->second = offd_num_cols++;
761  }
762 
763  // construct the global ParCSR matrix
764  A = hypre_ParCSRMatrixCreate(comm, glob_nrows, glob_ncols,
765  row_starts, col_starts, offd_num_cols,
766  diag_nnz, offd_nnz);
767  hypre_ParCSRMatrixInitialize(A);
768 
769  HYPRE_Int *diag_i, *diag_j, *offd_i, *offd_j, *offd_col_map;
770  double *diag_data, *offd_data;
771  diag_i = A->diag->i;
772  diag_j = A->diag->j;
773  diag_data = A->diag->data;
774  offd_i = A->offd->i;
775  offd_j = A->offd->j;
776  offd_data = A->offd->data;
777  offd_col_map = A->col_map_offd;
778 
779  diag_nnz = offd_nnz = 0;
780  for (HYPRE_Int i = 0, j = 0; i < nrows; i++)
781  {
782  diag_i[i] = diag_nnz;
783  offd_i[i] = offd_nnz;
784  for (HYPRE_Int j_end = I[i+1]; j < j_end; j++)
785  {
786  HYPRE_Int glob_col = J[j];
787  if (my_col_start <= glob_col && glob_col < my_col_end)
788  {
789  diag_j[diag_nnz] = glob_col - my_col_start;
790  diag_data[diag_nnz] = data[j];
791  diag_nnz++;
792  }
793  else
794  {
795  offd_j[offd_nnz] = offd_map[glob_col];
796  offd_data[offd_nnz] = data[j];
797  offd_nnz++;
798  }
799  }
800  }
801  diag_i[nrows] = diag_nnz;
802  offd_i[nrows] = offd_nnz;
803  for (map<HYPRE_Int, HYPRE_Int>::iterator it = offd_map.begin();
804  it != offd_map.end(); ++it)
805  {
806  offd_col_map[it->second] = it->first;
807  }
808 
809  hypre_ParCSRMatrixSetNumNonzeros(A);
810  /* Make sure that the first entry in each row is the diagonal one. */
811  if (row_starts == col_starts)
812  {
813  hypre_CSRMatrixReorder(hypre_ParCSRMatrixDiag(A));
814  }
815  hypre_MatvecCommPkgCreate(A);
816 
817  height = GetNumRows();
818  width = GetNumCols();
819 }
820 
822 {
823  Destroy();
824  Init();
825  A = master.A;
826  ParCSROwner = 0;
827  height = master.GetNumRows();
828  width = master.GetNumCols();
829 }
830 
831 hypre_ParCSRMatrix* HypreParMatrix::StealData()
832 {
833  // Only safe when (diagOwner == -1 && offdOwner == -1 && colMapOwner == -1)
834  // Otherwise, there may be memory leaks or hypre may destroy arrays allocated
835  // with operator new.
836  MFEM_ASSERT(diagOwner == -1 && offdOwner == -1 && colMapOwner == -1, "");
837  MFEM_ASSERT(ParCSROwner, "");
838  hypre_ParCSRMatrix *R = A;
839  A = NULL;
840  Destroy();
841  Init();
842  return R;
843 }
844 
846 {
847  if (!A || hypre_ParCSRMatrixOwnsRowStarts(A) ||
848  (hypre_ParCSRMatrixRowStarts(A) == hypre_ParCSRMatrixColStarts(A) &&
849  hypre_ParCSRMatrixOwnsColStarts(A)))
850  {
851  return;
852  }
853 
854  int row_starts_size;
855  if (HYPRE_AssumedPartitionCheck())
856  {
857  row_starts_size = 2;
858  }
859  else
860  {
861  MPI_Comm_size(hypre_ParCSRMatrixComm(A), &row_starts_size);
862  row_starts_size++; // num_proc + 1
863  }
864 
865  HYPRE_Int *old_row_starts = hypre_ParCSRMatrixRowStarts(A);
866  HYPRE_Int *new_row_starts = mfem_hypre_CTAlloc(HYPRE_Int, row_starts_size);
867  for (int i = 0; i < row_starts_size; i++)
868  {
869  new_row_starts[i] = old_row_starts[i];
870  }
871 
872  hypre_ParCSRMatrixRowStarts(A) = new_row_starts;
873  hypre_ParCSRMatrixOwnsRowStarts(A) = 1;
874 
875  if (hypre_ParCSRMatrixColStarts(A) == old_row_starts)
876  {
877  hypre_ParCSRMatrixColStarts(A) = new_row_starts;
878  hypre_ParCSRMatrixOwnsColStarts(A) = 0;
879  }
880 }
881 
883 {
884  if (!A || hypre_ParCSRMatrixOwnsColStarts(A) ||
885  (hypre_ParCSRMatrixRowStarts(A) == hypre_ParCSRMatrixColStarts(A) &&
886  hypre_ParCSRMatrixOwnsRowStarts(A)))
887  {
888  return;
889  }
890 
891  int col_starts_size;
892  if (HYPRE_AssumedPartitionCheck())
893  {
894  col_starts_size = 2;
895  }
896  else
897  {
898  MPI_Comm_size(hypre_ParCSRMatrixComm(A), &col_starts_size);
899  col_starts_size++; // num_proc + 1
900  }
901 
902  HYPRE_Int *old_col_starts = hypre_ParCSRMatrixColStarts(A);
903  HYPRE_Int *new_col_starts = mfem_hypre_CTAlloc(HYPRE_Int, col_starts_size);
904  for (int i = 0; i < col_starts_size; i++)
905  {
906  new_col_starts[i] = old_col_starts[i];
907  }
908 
909  hypre_ParCSRMatrixColStarts(A) = new_col_starts;
910 
911  if (hypre_ParCSRMatrixRowStarts(A) == old_col_starts)
912  {
913  hypre_ParCSRMatrixRowStarts(A) = new_col_starts;
914  hypre_ParCSRMatrixOwnsRowStarts(A) = 1;
915  hypre_ParCSRMatrixOwnsColStarts(A) = 0;
916  }
917  else
918  {
919  hypre_ParCSRMatrixOwnsColStarts(A) = 1;
920  }
921 }
922 
924 {
925  int size = Height();
926  diag.SetSize(size);
927  for (int j = 0; j < size; j++)
928  {
929  diag(j) = A->diag->data[A->diag->i[j]];
930  MFEM_ASSERT(A->diag->j[A->diag->i[j]] == j,
931  "the first entry in each row must be the diagonal one");
932  }
933 }
934 
935 static void MakeWrapper(const hypre_CSRMatrix *mat, SparseMatrix &wrapper)
936 {
937  HYPRE_Int nr = hypre_CSRMatrixNumRows(mat);
938  HYPRE_Int nc = hypre_CSRMatrixNumCols(mat);
939 #ifndef HYPRE_BIGINT
940  SparseMatrix tmp(hypre_CSRMatrixI(mat),
941  hypre_CSRMatrixJ(mat),
942  hypre_CSRMatrixData(mat),
943  nr, nc, false, false, false);
944 #else
945  HYPRE_Int nnz = hypre_CSRMatrixNumNonzeros(mat);
946  SparseMatrix tmp(DuplicateAs<int>(hypre_CSRMatrixI(mat), nr+1),
947  DuplicateAs<int>(hypre_CSRMatrixJ(mat), nnz),
948  hypre_CSRMatrixData(mat),
949  nr, nc, true, false, false);
950 #endif
951  wrapper.Swap(tmp);
952 }
953 
955 {
956  MakeWrapper(A->diag, diag);
957 }
958 
959 void HypreParMatrix::GetOffd(SparseMatrix &offd, HYPRE_Int* &cmap) const
960 {
961  MakeWrapper(A->offd, offd);
962  cmap = A->col_map_offd;
963 }
964 
966  bool interleaved_rows,
967  bool interleaved_cols) const
968 {
969  int nr = blocks.NumRows();
970  int nc = blocks.NumCols();
971 
972  hypre_ParCSRMatrix **hypre_blocks = new hypre_ParCSRMatrix*[nr * nc];
973  internal::hypre_ParCSRMatrixSplit(A, nr, nc, hypre_blocks,
974  interleaved_rows, interleaved_cols);
975 
976  for (int i = 0; i < nr; i++)
977  {
978  for (int j = 0; j < nc; j++)
979  {
980  blocks[i][j] = new HypreParMatrix(hypre_blocks[i*nc + j]);
981  }
982  }
983 
984  delete [] hypre_blocks;
985 }
986 
988 {
989  hypre_ParCSRMatrix * At;
990  hypre_ParCSRMatrixTranspose(A, &At, 1);
991  hypre_ParCSRMatrixSetNumNonzeros(At);
992 
993  hypre_MatvecCommPkgCreate(At);
994 
995  if ( M() == N() )
996  {
997  /* If the matrix is square, make sure that the first entry in each
998  row is the diagonal one. */
999  hypre_CSRMatrixReorder(hypre_ParCSRMatrixDiag(At));
1000  }
1001 
1002  return new HypreParMatrix(At);
1003 }
1004 
1006  double a, double b)
1007 {
1008  return hypre_ParCSRMatrixMatvec(a, A, x, b, y);
1009 }
1010 
1011 void HypreParMatrix::Mult(double a, const Vector &x, double b, Vector &y) const
1012 {
1013  MFEM_ASSERT(x.Size() == Width(), "invalid x.Size() = " << x.Size()
1014  << ", expected size = " << Width());
1015  MFEM_ASSERT(y.Size() == Height(), "invalid y.Size() = " << y.Size()
1016  << ", expected size = " << Height());
1017 
1018  if (X == NULL)
1019  {
1020  X = new HypreParVector(A->comm,
1021  GetGlobalNumCols(),
1022  x.GetData(),
1023  GetColStarts());
1024  Y = new HypreParVector(A->comm,
1025  GetGlobalNumRows(),
1026  y.GetData(),
1027  GetRowStarts());
1028  }
1029  else
1030  {
1031  X->SetData(x.GetData());
1032  Y->SetData(y.GetData());
1033  }
1034 
1035  hypre_ParCSRMatrixMatvec(a, A, *X, b, *Y);
1036 }
1037 
1038 void HypreParMatrix::MultTranspose(double a, const Vector &x,
1039  double b, Vector &y) const
1040 {
1041  MFEM_ASSERT(x.Size() == Height(), "invalid x.Size() = " << x.Size()
1042  << ", expected size = " << Height());
1043  MFEM_ASSERT(y.Size() == Width(), "invalid y.Size() = " << y.Size()
1044  << ", expected size = " << Width());
1045 
1046  // Note: x has the dimensions of Y (height), and
1047  // y has the dimensions of X (width)
1048  if (X == NULL)
1049  {
1050  X = new HypreParVector(A->comm,
1051  GetGlobalNumCols(),
1052  y.GetData(),
1053  GetColStarts());
1054  Y = new HypreParVector(A->comm,
1055  GetGlobalNumRows(),
1056  x.GetData(),
1057  GetRowStarts());
1058  }
1059  else
1060  {
1061  X->SetData(y.GetData());
1062  Y->SetData(x.GetData());
1063  }
1064 
1065  hypre_ParCSRMatrixMatvecT(a, A, *Y, b, *X);
1066 }
1067 
1068 HYPRE_Int HypreParMatrix::Mult(HYPRE_ParVector x, HYPRE_ParVector y,
1069  double a, double b)
1070 {
1071  return hypre_ParCSRMatrixMatvec(a, A, (hypre_ParVector *) x, b,
1072  (hypre_ParVector *) y);
1073 }
1074 
1076  double a, double b)
1077 {
1078  return hypre_ParCSRMatrixMatvecT(a, A, x, b, y);
1079 }
1080 
1082  HYPRE_Int* row_starts) const
1083 {
1084  const bool assumed_partition = HYPRE_AssumedPartitionCheck();
1085  const bool row_starts_given = (row_starts != NULL);
1086  if (!row_starts_given)
1087  {
1088  row_starts = hypre_ParCSRMatrixRowStarts(A);
1089  MFEM_VERIFY(D.Height() == hypre_CSRMatrixNumRows(A->diag),
1090  "the matrix D is NOT compatible with the row starts of"
1091  " this HypreParMatrix, row_starts must be given.");
1092  }
1093  else
1094  {
1095  int offset;
1096  if (assumed_partition)
1097  {
1098  offset = 0;
1099  }
1100  else
1101  {
1102  MPI_Comm_rank(GetComm(), &offset);
1103  }
1104  int local_num_rows = row_starts[offset+1]-row_starts[offset];
1105  MFEM_VERIFY(local_num_rows == D.Height(), "the number of rows in D is "
1106  " not compatible with the given row_starts");
1107  }
1108  // D.Width() will be checked for compatibility by the SparseMatrix
1109  // multiplication function, mfem::Mult(), called below.
1110 
1111  int part_size;
1112  HYPRE_Int global_num_rows;
1113  if (assumed_partition)
1114  {
1115  part_size = 2;
1116  if (row_starts_given)
1117  {
1118  global_num_rows = row_starts[2];
1119  // Here, we use row_starts[2], so row_starts must come from the
1120  // methods GetDofOffsets/GetTrueDofOffsets of ParFiniteElementSpace
1121  // (HYPRE's partitions have only 2 entries).
1122  }
1123  else
1124  {
1125  global_num_rows = hypre_ParCSRMatrixGlobalNumRows(A);
1126  }
1127  }
1128  else
1129  {
1130  MPI_Comm_size(GetComm(), &part_size);
1131  global_num_rows = row_starts[part_size];
1132  part_size++;
1133  }
1134 
1135  HYPRE_Int *col_starts = hypre_ParCSRMatrixColStarts(A);
1136  HYPRE_Int *col_map_offd;
1137 
1138  // get the diag and offd blocks as SparseMatrix wrappers
1139  SparseMatrix A_diag, A_offd;
1140  GetDiag(A_diag);
1141  GetOffd(A_offd, col_map_offd);
1142 
1143  // multiply the blocks with D and create a new HypreParMatrix
1144  SparseMatrix* DA_diag = mfem::Mult(D, A_diag);
1145  SparseMatrix* DA_offd = mfem::Mult(D, A_offd);
1146 
1147  HypreParMatrix* DA =
1148  new HypreParMatrix(GetComm(),
1149  global_num_rows, hypre_ParCSRMatrixGlobalNumCols(A),
1150  DuplicateAs<HYPRE_Int>(row_starts, part_size, false),
1151  DuplicateAs<HYPRE_Int>(col_starts, part_size, false),
1152  DA_diag, DA_offd,
1153  DuplicateAs<HYPRE_Int>(col_map_offd, A_offd.Width()));
1154 
1155  // When HYPRE_BIGINT is defined, we want DA_{diag,offd} to delete their I and
1156  // J arrays but not their data arrays; when HYPRE_BIGINT is not defined, we
1157  // don't want DA_{diag,offd} to delete anything.
1158 #ifndef HYPRE_BIGINT
1159  DA_diag->LoseData();
1160  DA_offd->LoseData();
1161 #else
1162  DA_diag->SetDataOwner(false);
1163  DA_offd->SetDataOwner(false);
1164 #endif
1165 
1166  delete DA_diag;
1167  delete DA_offd;
1168 
1169  hypre_ParCSRMatrixSetRowStartsOwner(DA->A, 1);
1170  hypre_ParCSRMatrixSetColStartsOwner(DA->A, 1);
1171 
1172  DA->diagOwner = DA->offdOwner = 3;
1173  DA->colMapOwner = 1;
1174 
1175  return DA;
1176 }
1177 
1179 {
1180  if (hypre_CSRMatrixNumRows(A->diag) != hypre_CSRMatrixNumRows(A->offd))
1181  {
1182  mfem_error("Row does not match");
1183  }
1184 
1185  if (hypre_CSRMatrixNumRows(A->diag) != diag.Size())
1186  {
1187  mfem_error("Note the Vector diag is not of compatible dimensions with A\n");
1188  }
1189 
1190  int size = Height();
1191  double *Adiag_data = hypre_CSRMatrixData(A->diag);
1192  HYPRE_Int *Adiag_i = hypre_CSRMatrixI(A->diag);
1193 
1194 
1195  double *Aoffd_data = hypre_CSRMatrixData(A->offd);
1196  HYPRE_Int *Aoffd_i = hypre_CSRMatrixI(A->offd);
1197  double val;
1198  HYPRE_Int jj;
1199  for (int i(0); i < size; ++i)
1200  {
1201  val = diag[i];
1202  for (jj = Adiag_i[i]; jj < Adiag_i[i+1]; ++jj)
1203  {
1204  Adiag_data[jj] *= val;
1205  }
1206  for (jj = Aoffd_i[i]; jj < Aoffd_i[i+1]; ++jj)
1207  {
1208  Aoffd_data[jj] *= val;
1209  }
1210  }
1211 }
1212 
1214 {
1215  if (hypre_CSRMatrixNumRows(A->diag) != hypre_CSRMatrixNumRows(A->offd))
1216  {
1217  mfem_error("Row does not match");
1218  }
1219 
1220  if (hypre_CSRMatrixNumRows(A->diag) != diag.Size())
1221  {
1222  mfem_error("Note the Vector diag is not of compatible dimensions with A\n");
1223  }
1224 
1225  int size = Height();
1226  double *Adiag_data = hypre_CSRMatrixData(A->diag);
1227  HYPRE_Int *Adiag_i = hypre_CSRMatrixI(A->diag);
1228 
1229 
1230  double *Aoffd_data = hypre_CSRMatrixData(A->offd);
1231  HYPRE_Int *Aoffd_i = hypre_CSRMatrixI(A->offd);
1232  double val;
1233  HYPRE_Int jj;
1234  for (int i(0); i < size; ++i)
1235  {
1236 #ifdef MFEM_DEBUG
1237  if (0.0 == diag(i))
1238  {
1239  mfem_error("HypreParMatrix::InvDiagScale : Division by 0");
1240  }
1241 #endif
1242  val = 1./diag(i);
1243  for (jj = Adiag_i[i]; jj < Adiag_i[i+1]; ++jj)
1244  {
1245  Adiag_data[jj] *= val;
1246  }
1247  for (jj = Aoffd_i[i]; jj < Aoffd_i[i+1]; ++jj)
1248  {
1249  Aoffd_data[jj] *= val;
1250  }
1251  }
1252 }
1253 
1255 {
1256  if (hypre_CSRMatrixNumRows(A->diag) != hypre_CSRMatrixNumRows(A->offd))
1257  {
1258  mfem_error("Row does not match");
1259  }
1260 
1261  HYPRE_Int size=hypre_CSRMatrixNumRows(A->diag);
1262  HYPRE_Int jj;
1263 
1264  double *Adiag_data = hypre_CSRMatrixData(A->diag);
1265  HYPRE_Int *Adiag_i = hypre_CSRMatrixI(A->diag);
1266  for (jj = 0; jj < Adiag_i[size]; ++jj)
1267  {
1268  Adiag_data[jj] *= s;
1269  }
1270 
1271  double *Aoffd_data = hypre_CSRMatrixData(A->offd);
1272  HYPRE_Int *Aoffd_i = hypre_CSRMatrixI(A->offd);
1273  for (jj = 0; jj < Aoffd_i[size]; ++jj)
1274  {
1275  Aoffd_data[jj] *= s;
1276  }
1277 }
1278 
1279 static void get_sorted_rows_cols(const Array<int> &rows_cols,
1280  Array<HYPRE_Int> &hypre_sorted)
1281 {
1282  hypre_sorted.SetSize(rows_cols.Size());
1283  bool sorted = true;
1284  for (int i = 0; i < rows_cols.Size(); i++)
1285  {
1286  hypre_sorted[i] = rows_cols[i];
1287  if (i && rows_cols[i-1] > rows_cols[i]) { sorted = false; }
1288  }
1289  if (!sorted) { hypre_sorted.Sort(); }
1290 }
1291 
1292 void HypreParMatrix::Threshold(double threshold)
1293 {
1294  int ierr = 0;
1295 
1296  MPI_Comm comm;
1297  hypre_CSRMatrix * csr_A;
1298  hypre_CSRMatrix * csr_A_wo_z;
1299  hypre_ParCSRMatrix * parcsr_A_ptr;
1300  HYPRE_Int * row_starts = NULL; HYPRE_Int * col_starts = NULL;
1301  HYPRE_Int row_start = -1; HYPRE_Int row_end = -1;
1302  HYPRE_Int col_start = -1; HYPRE_Int col_end = -1;
1303 
1304  comm = hypre_ParCSRMatrixComm(A);
1305 
1306  ierr += hypre_ParCSRMatrixGetLocalRange(A,
1307  &row_start,&row_end,
1308  &col_start,&col_end );
1309 
1310  row_starts = hypre_ParCSRMatrixRowStarts(A);
1311  col_starts = hypre_ParCSRMatrixColStarts(A);
1312 
1313  bool old_owns_row = hypre_ParCSRMatrixOwnsRowStarts(A);
1314  bool old_owns_col = hypre_ParCSRMatrixOwnsColStarts(A);
1315  HYPRE_Int global_num_rows = hypre_ParCSRMatrixGlobalNumRows(A);
1316  HYPRE_Int global_num_cols = hypre_ParCSRMatrixGlobalNumCols(A);
1317  parcsr_A_ptr = hypre_ParCSRMatrixCreate(comm, global_num_rows,
1318  global_num_cols,
1319  row_starts, col_starts,
1320  0, 0, 0);
1321  hypre_ParCSRMatrixOwnsRowStarts(parcsr_A_ptr) = old_owns_row;
1322  hypre_ParCSRMatrixOwnsColStarts(parcsr_A_ptr) = old_owns_col;
1323  hypre_ParCSRMatrixOwnsRowStarts(A) = 0;
1324  hypre_ParCSRMatrixOwnsColStarts(A) = 0;
1325 
1326  csr_A = hypre_MergeDiagAndOffd(A);
1327 
1328  // Free A, if owned
1329  Destroy();
1330  Init();
1331 
1332  csr_A_wo_z = hypre_CSRMatrixDeleteZeros(csr_A,threshold);
1333 
1334  /* hypre_CSRMatrixDeleteZeros will return a NULL pointer rather than a usable
1335  CSR matrix if it finds no non-zeros */
1336  if (csr_A_wo_z == NULL)
1337  {
1338  csr_A_wo_z = csr_A;
1339  }
1340  else
1341  {
1342  ierr += hypre_CSRMatrixDestroy(csr_A);
1343  }
1344 
1345  /* FIXME: GenerateDiagAndOffd() uses an int array of size equal to the number
1346  of columns in csr_A_wo_z which is the global number of columns in A. This
1347  does not scale well. */
1348  ierr += GenerateDiagAndOffd(csr_A_wo_z,parcsr_A_ptr,
1349  col_start,col_end);
1350 
1351  ierr += hypre_CSRMatrixDestroy(csr_A_wo_z);
1352 
1353  MFEM_VERIFY(ierr == 0, "");
1354 
1355  A = parcsr_A_ptr;
1356 
1357  hypre_ParCSRMatrixSetNumNonzeros(A);
1358  /* Make sure that the first entry in each row is the diagonal one. */
1359  if (row_starts == col_starts)
1360  {
1361  hypre_CSRMatrixReorder(hypre_ParCSRMatrixDiag(A));
1362  }
1363  hypre_MatvecCommPkgCreate(A);
1364  height = GetNumRows();
1365  width = GetNumCols();
1366 }
1367 
1369  const HypreParVector &X,
1370  HypreParVector &B)
1371 {
1372  Array<HYPRE_Int> rc_sorted;
1373  get_sorted_rows_cols(rows_cols, rc_sorted);
1374 
1375  internal::hypre_ParCSRMatrixEliminateAXB(
1376  A, rc_sorted.Size(), rc_sorted.GetData(), X, B);
1377 }
1378 
1380 {
1381  Array<HYPRE_Int> rc_sorted;
1382  get_sorted_rows_cols(rows_cols, rc_sorted);
1383 
1384  hypre_ParCSRMatrix* Ae;
1385  internal::hypre_ParCSRMatrixEliminateAAe(
1386  A, &Ae, rc_sorted.Size(), rc_sorted.GetData());
1387 
1388  return new HypreParMatrix(Ae);
1389 }
1390 
1391 void HypreParMatrix::Print(const char *fname, HYPRE_Int offi, HYPRE_Int offj)
1392 {
1393  hypre_ParCSRMatrixPrintIJ(A,offi,offj,fname);
1394 }
1395 
1396 void HypreParMatrix::Read(MPI_Comm comm, const char *fname)
1397 {
1398  Destroy();
1399  Init();
1400 
1401  HYPRE_Int base_i, base_j;
1402  hypre_ParCSRMatrixReadIJ(comm, fname, &base_i, &base_j, &A);
1403  hypre_ParCSRMatrixSetNumNonzeros(A);
1404 
1405  hypre_MatvecCommPkgCreate(A);
1406 
1407  height = GetNumRows();
1408  width = GetNumCols();
1409 }
1410 
1411 void HypreParMatrix::Read_IJMatrix(MPI_Comm comm, const char *fname)
1412 {
1413  Destroy();
1414  Init();
1415 
1416  HYPRE_IJMatrix A_ij;
1417  HYPRE_IJMatrixRead(fname, comm, 5555, &A_ij); // HYPRE_PARCSR = 5555
1418 
1419  HYPRE_ParCSRMatrix A_parcsr;
1420  HYPRE_IJMatrixGetObject(A_ij, (void**) &A_parcsr);
1421 
1422  A = (hypre_ParCSRMatrix*)A_parcsr;
1423 
1424  hypre_ParCSRMatrixSetNumNonzeros(A);
1425 
1426  hypre_MatvecCommPkgCreate(A);
1427 
1428  height = GetNumRows();
1429  width = GetNumCols();
1430 }
1431 
1432 void HypreParMatrix::PrintCommPkg(std::ostream &out) const
1433 {
1434  hypre_ParCSRCommPkg *comm_pkg = A->comm_pkg;
1435  MPI_Comm comm = A->comm;
1436  char c = '\0';
1437  const int tag = 46801;
1438  int myid, nproc;
1439  MPI_Comm_rank(comm, &myid);
1440  MPI_Comm_size(comm, &nproc);
1441 
1442  if (myid != 0)
1443  {
1444  MPI_Recv(&c, 1, MPI_CHAR, myid-1, tag, comm, MPI_STATUS_IGNORE);
1445  }
1446  else
1447  {
1448  out << "\nHypreParMatrix: hypre_ParCSRCommPkg:\n";
1449  }
1450  out << "Rank " << myid << ":\n"
1451  " number of sends = " << comm_pkg->num_sends <<
1452  " (" << sizeof(double)*comm_pkg->send_map_starts[comm_pkg->num_sends] <<
1453  " bytes)\n"
1454  " number of recvs = " << comm_pkg->num_recvs <<
1455  " (" << sizeof(double)*comm_pkg->recv_vec_starts[comm_pkg->num_recvs] <<
1456  " bytes)\n";
1457  if (myid != nproc-1)
1458  {
1459  out << std::flush;
1460  MPI_Send(&c, 1, MPI_CHAR, myid+1, tag, comm);
1461  }
1462  else
1463  {
1464  out << std::endl;
1465  }
1466  MPI_Barrier(comm);
1467 }
1468 
1469 void HypreParMatrix::Destroy()
1470 {
1471  if ( X != NULL ) { delete X; }
1472  if ( Y != NULL ) { delete Y; }
1473 
1474  if (A == NULL) { return; }
1475 
1476  if (diagOwner >= 0)
1477  {
1478  if (diagOwner & 1)
1479  {
1480  delete [] hypre_CSRMatrixI(A->diag);
1481  delete [] hypre_CSRMatrixJ(A->diag);
1482  }
1483  hypre_CSRMatrixI(A->diag) = NULL;
1484  hypre_CSRMatrixJ(A->diag) = NULL;
1485  if (diagOwner & 2)
1486  {
1487  delete [] hypre_CSRMatrixData(A->diag);
1488  }
1489  hypre_CSRMatrixData(A->diag) = NULL;
1490  }
1491  if (offdOwner >= 0)
1492  {
1493  if (offdOwner & 1)
1494  {
1495  delete [] hypre_CSRMatrixI(A->offd);
1496  delete [] hypre_CSRMatrixJ(A->offd);
1497  }
1498  hypre_CSRMatrixI(A->offd) = NULL;
1499  hypre_CSRMatrixJ(A->offd) = NULL;
1500  if (offdOwner & 2)
1501  {
1502  delete [] hypre_CSRMatrixData(A->offd);
1503  }
1504  hypre_CSRMatrixData(A->offd) = NULL;
1505  }
1506  if (colMapOwner >= 0)
1507  {
1508  if (colMapOwner & 1)
1509  {
1510  delete [] hypre_ParCSRMatrixColMapOffd(A);
1511  }
1512  hypre_ParCSRMatrixColMapOffd(A) = NULL;
1513  }
1514 
1515  if (ParCSROwner)
1516  {
1517  hypre_ParCSRMatrixDestroy(A);
1518  }
1519 }
1520 
1522  double beta, const HypreParMatrix &B)
1523 {
1524  hypre_ParCSRMatrix *C_hypre =
1525  internal::hypre_ParCSRMatrixAdd(const_cast<HypreParMatrix &>(A),
1526  const_cast<HypreParMatrix &>(B));
1527  MFEM_VERIFY(C_hypre, "error in hypre_ParCSRMatrixAdd");
1528 
1529  hypre_MatvecCommPkgCreate(C_hypre);
1530  HypreParMatrix *C = new HypreParMatrix(C_hypre);
1531  *C = 0.0;
1532  C->Add(alpha, A);
1533  C->Add(beta, B);
1534 
1535  return C;
1536 }
1537 
1539 {
1540  hypre_ParCSRMatrix * ab;
1541  ab = hypre_ParMatmul(*A,*B);
1542  hypre_ParCSRMatrixSetNumNonzeros(ab);
1543 
1544  hypre_MatvecCommPkgCreate(ab);
1545 
1546  return new HypreParMatrix(ab);
1547 }
1548 
1550 {
1551  hypre_ParCSRMatrix * C = internal::hypre_ParCSRMatrixAdd(*A,*B);
1552 
1553  hypre_MatvecCommPkgCreate(C);
1554 
1555  return new HypreParMatrix(C);
1556 }
1557 
1559 {
1560  HYPRE_Int P_owns_its_col_starts =
1561  hypre_ParCSRMatrixOwnsColStarts((hypre_ParCSRMatrix*)(*P));
1562 
1563  hypre_ParCSRMatrix * rap;
1564  hypre_BoomerAMGBuildCoarseOperator(*P,*A,*P,&rap);
1565  hypre_ParCSRMatrixSetNumNonzeros(rap);
1566  // hypre_MatvecCommPkgCreate(rap);
1567 
1568  /* Warning: hypre_BoomerAMGBuildCoarseOperator steals the col_starts
1569  from P (even if it does not own them)! */
1570  hypre_ParCSRMatrixSetRowStartsOwner(rap,0);
1571  hypre_ParCSRMatrixSetColStartsOwner(rap,0);
1572 
1573  if (P_owns_its_col_starts)
1574  {
1575  hypre_ParCSRMatrixSetColStartsOwner(*P, 1);
1576  }
1577 
1578  return new HypreParMatrix(rap);
1579 }
1580 
1582  const HypreParMatrix *P)
1583 {
1584  HYPRE_Int P_owns_its_col_starts =
1585  hypre_ParCSRMatrixOwnsColStarts((hypre_ParCSRMatrix*)(*P));
1586  HYPRE_Int Rt_owns_its_col_starts =
1587  hypre_ParCSRMatrixOwnsColStarts((hypre_ParCSRMatrix*)(*Rt));
1588 
1589  hypre_ParCSRMatrix * rap;
1590  hypre_BoomerAMGBuildCoarseOperator(*Rt,*A,*P,&rap);
1591 
1592  hypre_ParCSRMatrixSetNumNonzeros(rap);
1593  // hypre_MatvecCommPkgCreate(rap);
1594 
1595  /* Warning: hypre_BoomerAMGBuildCoarseOperator steals the col_starts
1596  from Rt and P (even if they do not own them)! */
1597  hypre_ParCSRMatrixSetRowStartsOwner(rap,0);
1598  hypre_ParCSRMatrixSetColStartsOwner(rap,0);
1599 
1600  if (P_owns_its_col_starts)
1601  {
1602  hypre_ParCSRMatrixSetColStartsOwner(*P, 1);
1603  }
1604  if (Rt_owns_its_col_starts)
1605  {
1606  hypre_ParCSRMatrixSetColStartsOwner(*Rt, 1);
1607  }
1608 
1609  return new HypreParMatrix(rap);
1610 }
1611 
1613  const Array<int> &ess_dof_list,
1614  const Vector &X, Vector &B)
1615 {
1616  // B -= Ae*X
1617  Ae.Mult(-1.0, X, 1.0, B);
1618 
1619  hypre_CSRMatrix *A_diag = hypre_ParCSRMatrixDiag((hypre_ParCSRMatrix *)A);
1620  double *data = hypre_CSRMatrixData(A_diag);
1621  HYPRE_Int *I = hypre_CSRMatrixI(A_diag);
1622 #ifdef MFEM_DEBUG
1623  HYPRE_Int *J = hypre_CSRMatrixJ(A_diag);
1624  hypre_CSRMatrix *A_offd = hypre_ParCSRMatrixOffd((hypre_ParCSRMatrix *)A);
1625  HYPRE_Int *I_offd = hypre_CSRMatrixI(A_offd);
1626  double *data_offd = hypre_CSRMatrixData(A_offd);
1627 #endif
1628 
1629  for (int i = 0; i < ess_dof_list.Size(); i++)
1630  {
1631  int r = ess_dof_list[i];
1632  B(r) = data[I[r]] * X(r);
1633 #ifdef MFEM_DEBUG
1634  // Check that in the rows specified by the ess_dof_list, the matrix A has
1635  // only one entry -- the diagonal.
1636  // if (I[r+1] != I[r]+1 || J[I[r]] != r || I_offd[r] != I_offd[r+1])
1637  if (J[I[r]] != r)
1638  {
1639  MFEM_ABORT("the diagonal entry must be the first entry in the row!");
1640  }
1641  for (int j = I[r]+1; j < I[r+1]; j++)
1642  {
1643  if (data[j] != 0.0)
1644  {
1645  MFEM_ABORT("all off-diagonal entries must be zero!");
1646  }
1647  }
1648  for (int j = I_offd[r]; j < I_offd[r+1]; j++)
1649  {
1650  if (data_offd[j] != 0.0)
1651  {
1652  MFEM_ABORT("all off-diagonal entries must be zero!");
1653  }
1654  }
1655 #endif
1656  }
1657 }
1658 
1659 // Taubin or "lambda-mu" scheme, which alternates between positive and
1660 // negative step sizes to approximate low-pass filter effect.
1661 
1662 int ParCSRRelax_Taubin(hypre_ParCSRMatrix *A, // matrix to relax with
1663  hypre_ParVector *f, // right-hand side
1664  double lambda,
1665  double mu,
1666  int N,
1667  double max_eig,
1668  hypre_ParVector *u, // initial/updated approximation
1669  hypre_ParVector *r // another temp vector
1670  )
1671 {
1672  hypre_CSRMatrix *A_diag = hypre_ParCSRMatrixDiag(A);
1673  HYPRE_Int num_rows = hypre_CSRMatrixNumRows(A_diag);
1674 
1675  double *u_data = hypre_VectorData(hypre_ParVectorLocalVector(u));
1676  double *r_data = hypre_VectorData(hypre_ParVectorLocalVector(r));
1677 
1678  for (int i = 0; i < N; i++)
1679  {
1680  // get residual: r = f - A*u
1681  hypre_ParVectorCopy(f, r);
1682  hypre_ParCSRMatrixMatvec(-1.0, A, u, 1.0, r);
1683 
1684  double coef;
1685  (0 == (i % 2)) ? coef = lambda : coef = mu;
1686 
1687  for (HYPRE_Int j = 0; j < num_rows; j++)
1688  {
1689  u_data[j] += coef*r_data[j] / max_eig;
1690  }
1691  }
1692 
1693  return 0;
1694 }
1695 
1696 // FIR scheme, which uses Chebyshev polynomials and a window function
1697 // to approximate a low-pass step filter.
1698 
1699 int ParCSRRelax_FIR(hypre_ParCSRMatrix *A, // matrix to relax with
1700  hypre_ParVector *f, // right-hand side
1701  double max_eig,
1702  int poly_order,
1703  double* fir_coeffs,
1704  hypre_ParVector *u, // initial/updated approximation
1705  hypre_ParVector *x0, // temporaries
1706  hypre_ParVector *x1,
1707  hypre_ParVector *x2,
1708  hypre_ParVector *x3)
1709 
1710 {
1711  hypre_CSRMatrix *A_diag = hypre_ParCSRMatrixDiag(A);
1712  HYPRE_Int num_rows = hypre_CSRMatrixNumRows(A_diag);
1713 
1714  double *u_data = hypre_VectorData(hypre_ParVectorLocalVector(u));
1715 
1716  double *x0_data = hypre_VectorData(hypre_ParVectorLocalVector(x0));
1717  double *x1_data = hypre_VectorData(hypre_ParVectorLocalVector(x1));
1718  double *x2_data = hypre_VectorData(hypre_ParVectorLocalVector(x2));
1719  double *x3_data = hypre_VectorData(hypre_ParVectorLocalVector(x3));
1720 
1721  hypre_ParVectorCopy(u, x0);
1722 
1723  // x1 = f -A*x0/max_eig
1724  hypre_ParVectorCopy(f, x1);
1725  hypre_ParCSRMatrixMatvec(-1.0, A, x0, 1.0, x1);
1726 
1727  for (HYPRE_Int i = 0; i < num_rows; i++)
1728  {
1729  x1_data[i] /= -max_eig;
1730  }
1731 
1732  // x1 = x0 -x1
1733  for (HYPRE_Int i = 0; i < num_rows; i++)
1734  {
1735  x1_data[i] = x0_data[i] -x1_data[i];
1736  }
1737 
1738  // x3 = f0*x0 +f1*x1
1739  for (HYPRE_Int i = 0; i < num_rows; i++)
1740  {
1741  x3_data[i] = fir_coeffs[0]*x0_data[i] +fir_coeffs[1]*x1_data[i];
1742  }
1743 
1744  for (int n = 2; n <= poly_order; n++)
1745  {
1746  // x2 = f - A*x1/max_eig
1747  hypre_ParVectorCopy(f, x2);
1748  hypre_ParCSRMatrixMatvec(-1.0, A, x1, 1.0, x2);
1749 
1750  for (HYPRE_Int i = 0; i < num_rows; i++)
1751  {
1752  x2_data[i] /= -max_eig;
1753  }
1754 
1755  // x2 = (x1-x0) +(x1-2*x2)
1756  // x3 = x3 +f[n]*x2
1757  // x0 = x1
1758  // x1 = x2
1759 
1760  for (HYPRE_Int i = 0; i < num_rows; i++)
1761  {
1762  x2_data[i] = (x1_data[i]-x0_data[i]) +(x1_data[i]-2*x2_data[i]);
1763  x3_data[i] += fir_coeffs[n]*x2_data[i];
1764  x0_data[i] = x1_data[i];
1765  x1_data[i] = x2_data[i];
1766  }
1767  }
1768 
1769  for (HYPRE_Int i = 0; i < num_rows; i++)
1770  {
1771  u_data[i] = x3_data[i];
1772  }
1773 
1774  return 0;
1775 }
1776 
1778 {
1779  type = 2;
1780  relax_times = 1;
1781  relax_weight = 1.0;
1782  omega = 1.0;
1783  poly_order = 2;
1784  poly_fraction = .3;
1785  lambda = 0.5;
1786  mu = -0.5;
1787  taubin_iter = 40;
1788 
1789  l1_norms = NULL;
1790  pos_l1_norms = false;
1791  B = X = V = Z = NULL;
1792  X0 = X1 = NULL;
1793  fir_coeffs = NULL;
1794 }
1795 
1797  int _relax_times, double _relax_weight, double _omega,
1798  int _poly_order, double _poly_fraction)
1799 {
1800  type = _type;
1801  relax_times = _relax_times;
1802  relax_weight = _relax_weight;
1803  omega = _omega;
1804  poly_order = _poly_order;
1805  poly_fraction = _poly_fraction;
1806 
1807  l1_norms = NULL;
1808  pos_l1_norms = false;
1809  B = X = V = Z = NULL;
1810  X0 = X1 = NULL;
1811  fir_coeffs = NULL;
1812 
1813  SetOperator(_A);
1814 }
1815 
1816 void HypreSmoother::SetType(HypreSmoother::Type _type, int _relax_times)
1817 {
1818  type = static_cast<int>(_type);
1819  relax_times = _relax_times;
1820 }
1821 
1822 void HypreSmoother::SetSOROptions(double _relax_weight, double _omega)
1823 {
1824  relax_weight = _relax_weight;
1825  omega = _omega;
1826 }
1827 
1828 void HypreSmoother::SetPolyOptions(int _poly_order, double _poly_fraction)
1829 {
1830  poly_order = _poly_order;
1831  poly_fraction = _poly_fraction;
1832 }
1833 
1834 void HypreSmoother::SetTaubinOptions(double _lambda, double _mu,
1835  int _taubin_iter)
1836 {
1837  lambda = _lambda;
1838  mu = _mu;
1839  taubin_iter = _taubin_iter;
1840 }
1841 
1842 void HypreSmoother::SetWindowByName(const char* name)
1843 {
1844  double a = -1, b, c;
1845  if (!strcmp(name,"Rectangular")) { a = 1.0, b = 0.0, c = 0.0; }
1846  if (!strcmp(name,"Hanning")) { a = 0.5, b = 0.5, c = 0.0; }
1847  if (!strcmp(name,"Hamming")) { a = 0.54, b = 0.46, c = 0.0; }
1848  if (!strcmp(name,"Blackman")) { a = 0.42, b = 0.50, c = 0.08; }
1849  if (a < 0)
1850  {
1851  mfem_error("HypreSmoother::SetWindowByName : name not recognized!");
1852  }
1853 
1854  SetWindowParameters(a, b, c);
1855 }
1856 
1857 void HypreSmoother::SetWindowParameters(double a, double b, double c)
1858 {
1859  window_params[0] = a;
1860  window_params[1] = b;
1861  window_params[2] = c;
1862 }
1863 
1865 {
1866  A = const_cast<HypreParMatrix *>(dynamic_cast<const HypreParMatrix *>(&op));
1867  if (A == NULL)
1868  {
1869  mfem_error("HypreSmoother::SetOperator : not HypreParMatrix!");
1870  }
1871 
1872  height = A->Height();
1873  width = A->Width();
1874 
1875  if (B) { delete B; }
1876  if (X) { delete X; }
1877  if (V) { delete V; }
1878  if (Z) { delete Z; }
1879  if (l1_norms)
1880  {
1881  mfem_hypre_TFree(l1_norms);
1882  }
1883  delete X0;
1884  delete X1;
1885 
1886  X1 = X0 = Z = V = B = X = NULL;
1887 
1888  if (type >= 1 && type <= 4)
1889  {
1890  hypre_ParCSRComputeL1Norms(*A, type, NULL, &l1_norms);
1891  }
1892  else if (type == 5)
1893  {
1894  l1_norms = mfem_hypre_CTAlloc(double, height);
1895  Vector ones(height), diag(l1_norms, height);
1896  ones = 1.0;
1897  A->Mult(ones, diag);
1898  type = 1;
1899  }
1900  else
1901  {
1902  l1_norms = NULL;
1903  }
1904  if (l1_norms && pos_l1_norms)
1905  {
1906  for (int i = 0; i < height; i++)
1907  {
1908  l1_norms[i] = std::abs(l1_norms[i]);
1909  }
1910  }
1911 
1912  if (type == 16)
1913  {
1914  poly_scale = 1;
1915  hypre_ParCSRMaxEigEstimateCG(*A, poly_scale, 10,
1917  Z = new HypreParVector(*A);
1918  }
1919  else if (type == 1001 || type == 1002)
1920  {
1921  poly_scale = 0;
1922  hypre_ParCSRMaxEigEstimateCG(*A, poly_scale, 10,
1924 
1925  // The Taubin and FIR polynomials are defined on [0, 2]
1926  max_eig_est /= 2;
1927 
1928  // Compute window function, Chebyshev coefficients, and allocate temps.
1929  if (type == 1002)
1930  {
1931  // Temporaries for Chebyshev recursive evaluation
1932  Z = new HypreParVector(*A);
1933  X0 = new HypreParVector(*A);
1934  X1 = new HypreParVector(*A);
1935 
1937  }
1938  }
1939 }
1940 
1942 {
1943  if (fir_coeffs)
1944  {
1945  delete [] fir_coeffs;
1946  }
1947 
1948  fir_coeffs = new double[poly_order+1];
1949 
1950  double* window_coeffs = new double[poly_order+1];
1951  double* cheby_coeffs = new double[poly_order+1];
1952 
1953  double a = window_params[0];
1954  double b = window_params[1];
1955  double c = window_params[2];
1956  for (int i = 0; i <= poly_order; i++)
1957  {
1958  double t = (i*M_PI)/(poly_order+1);
1959  window_coeffs[i] = a + b*cos(t) +c*cos(2*t);
1960  }
1961 
1962  double k_pb = poly_fraction*max_eig;
1963  double theta_pb = acos(1.0 -0.5*k_pb);
1964  double sigma = 0.0;
1965  cheby_coeffs[0] = (theta_pb +sigma)/M_PI;
1966  for (int i = 1; i <= poly_order; i++)
1967  {
1968  double t = i*(theta_pb+sigma);
1969  cheby_coeffs[i] = 2.0*sin(t)/(i*M_PI);
1970  }
1971 
1972  for (int i = 0; i <= poly_order; i++)
1973  {
1974  fir_coeffs[i] = window_coeffs[i]*cheby_coeffs[i];
1975  }
1976 
1977  delete[] window_coeffs;
1978  delete[] cheby_coeffs;
1979 }
1980 
1982 {
1983  if (A == NULL)
1984  {
1985  mfem_error("HypreSmoother::Mult (...) : HypreParMatrix A is missing");
1986  return;
1987  }
1988 
1989  if (!iterative_mode)
1990  {
1991  if (type == 0 && relax_times == 1)
1992  {
1993  HYPRE_ParCSRDiagScale(NULL, *A, b, x);
1994  if (relax_weight != 1.0)
1995  {
1996  x *= relax_weight;
1997  }
1998  return;
1999  }
2000  x = 0.0;
2001  }
2002 
2003  if (V == NULL)
2004  {
2005  V = new HypreParVector(*A);
2006  }
2007 
2008  if (type == 1001)
2009  {
2010  for (int sweep = 0; sweep < relax_times; sweep++)
2011  {
2013  max_eig_est,
2014  x, *V);
2015  }
2016  }
2017  else if (type == 1002)
2018  {
2019  for (int sweep = 0; sweep < relax_times; sweep++)
2020  {
2021  ParCSRRelax_FIR(*A, b,
2022  max_eig_est,
2023  poly_order,
2024  fir_coeffs,
2025  x,
2026  *X0, *X1, *V, *Z);
2027  }
2028  }
2029  else
2030  {
2031  if (Z == NULL)
2032  hypre_ParCSRRelax(*A, b, type,
2035  x, *V, NULL);
2036  else
2037  hypre_ParCSRRelax(*A, b, type,
2040  x, *V, *Z);
2041  }
2042 }
2043 
2044 void HypreSmoother::Mult(const Vector &b, Vector &x) const
2045 {
2046  if (A == NULL)
2047  {
2048  mfem_error("HypreSmoother::Mult (...) : HypreParMatrix A is missing");
2049  return;
2050  }
2051  if (B == NULL)
2052  {
2053  B = new HypreParVector(A->GetComm(),
2054  A -> GetGlobalNumRows(),
2055  b.GetData(),
2056  A -> GetRowStarts());
2057  X = new HypreParVector(A->GetComm(),
2058  A -> GetGlobalNumCols(),
2059  x.GetData(),
2060  A -> GetColStarts());
2061  }
2062  else
2063  {
2064  B -> SetData(b.GetData());
2065  X -> SetData(x.GetData());
2066  }
2067 
2068  Mult(*B, *X);
2069 }
2070 
2072 {
2073  if (B) { delete B; }
2074  if (X) { delete X; }
2075  if (V) { delete V; }
2076  if (Z) { delete Z; }
2077  if (l1_norms)
2078  {
2079  mfem_hypre_TFree(l1_norms);
2080  }
2081  if (fir_coeffs)
2082  {
2083  delete [] fir_coeffs;
2084  }
2085  if (X0) { delete X0; }
2086  if (X1) { delete X1; }
2087 }
2088 
2089 
2091 {
2092  A = NULL;
2093  setup_called = 0;
2094  B = X = NULL;
2095 }
2096 
2098  : Solver(_A->Height(), _A->Width())
2099 {
2100  A = _A;
2101  setup_called = 0;
2102  B = X = NULL;
2103 }
2104 
2106 {
2107  if (A == NULL)
2108  {
2109  mfem_error("HypreSolver::Mult (...) : HypreParMatrix A is missing");
2110  return;
2111  }
2112  if (!setup_called)
2113  {
2114  SetupFcn()(*this, *A, b, x);
2115  setup_called = 1;
2116  }
2117 
2118  if (!iterative_mode)
2119  {
2120  x = 0.0;
2121  }
2122  SolveFcn()(*this, *A, b, x);
2123 }
2124 
2125 void HypreSolver::Mult(const Vector &b, Vector &x) const
2126 {
2127  if (A == NULL)
2128  {
2129  mfem_error("HypreSolver::Mult (...) : HypreParMatrix A is missing");
2130  return;
2131  }
2132  if (B == NULL)
2133  {
2134  B = new HypreParVector(A->GetComm(),
2135  A -> GetGlobalNumRows(),
2136  b.GetData(),
2137  A -> GetRowStarts());
2138  X = new HypreParVector(A->GetComm(),
2139  A -> GetGlobalNumCols(),
2140  x.GetData(),
2141  A -> GetColStarts());
2142  }
2143  else
2144  {
2145  B -> SetData(b.GetData());
2146  X -> SetData(x.GetData());
2147  }
2148 
2149  Mult(*B, *X);
2150 }
2151 
2153 {
2154  if (B) { delete B; }
2155  if (X) { delete X; }
2156 }
2157 
2158 
2160 {
2161  MPI_Comm comm;
2162 
2163  iterative_mode = true;
2164 
2165  HYPRE_ParCSRMatrixGetComm(*A, &comm);
2166 
2167  HYPRE_ParCSRPCGCreate(comm, &pcg_solver);
2168 }
2169 
2170 void HyprePCG::SetTol(double tol)
2171 {
2172  HYPRE_PCGSetTol(pcg_solver, tol);
2173 }
2174 
2175 void HyprePCG::SetMaxIter(int max_iter)
2176 {
2177  HYPRE_PCGSetMaxIter(pcg_solver, max_iter);
2178 }
2179 
2180 void HyprePCG::SetLogging(int logging)
2181 {
2182  HYPRE_PCGSetLogging(pcg_solver, logging);
2183 }
2184 
2185 void HyprePCG::SetPrintLevel(int print_lvl)
2186 {
2187  HYPRE_ParCSRPCGSetPrintLevel(pcg_solver, print_lvl);
2188 }
2189 
2191 {
2192  HYPRE_ParCSRPCGSetPrecond(pcg_solver,
2193  precond.SolveFcn(),
2194  precond.SetupFcn(),
2195  precond);
2196 }
2197 
2198 void HyprePCG::SetResidualConvergenceOptions(int res_frequency, double rtol)
2199 {
2200  HYPRE_PCGSetTwoNorm(pcg_solver, 1);
2201  if (res_frequency > 0)
2202  {
2203  HYPRE_PCGSetRecomputeResidualP(pcg_solver, res_frequency);
2204  }
2205  if (rtol > 0.0)
2206  {
2207  HYPRE_PCGSetResidualTol(pcg_solver, rtol);
2208  }
2209 }
2210 
2212 {
2213  int myid;
2214  HYPRE_Int time_index = 0;
2215  HYPRE_Int num_iterations;
2216  double final_res_norm;
2217  MPI_Comm comm;
2218  HYPRE_Int print_level;
2219 
2220  HYPRE_PCGGetPrintLevel(pcg_solver, &print_level);
2221  HYPRE_ParCSRPCGSetPrintLevel(pcg_solver, print_level%3);
2222 
2223  HYPRE_ParCSRMatrixGetComm(*A, &comm);
2224 
2225  if (!setup_called)
2226  {
2227  if (print_level > 0 && print_level < 3)
2228  {
2229  time_index = hypre_InitializeTiming("PCG Setup");
2230  hypre_BeginTiming(time_index);
2231  }
2232 
2233  HYPRE_ParCSRPCGSetup(pcg_solver, *A, b, x);
2234  setup_called = 1;
2235 
2236  if (print_level > 0 && print_level < 3)
2237  {
2238  hypre_EndTiming(time_index);
2239  hypre_PrintTiming("Setup phase times", comm);
2240  hypre_FinalizeTiming(time_index);
2241  hypre_ClearTiming();
2242  }
2243  }
2244 
2245  if (print_level > 0 && print_level < 3)
2246  {
2247  time_index = hypre_InitializeTiming("PCG Solve");
2248  hypre_BeginTiming(time_index);
2249  }
2250 
2251  if (!iterative_mode)
2252  {
2253  x = 0.0;
2254  }
2255 
2256  HYPRE_ParCSRPCGSolve(pcg_solver, *A, b, x);
2257 
2258  if (print_level > 0)
2259  {
2260  if (print_level < 3)
2261  {
2262  hypre_EndTiming(time_index);
2263  hypre_PrintTiming("Solve phase times", comm);
2264  hypre_FinalizeTiming(time_index);
2265  hypre_ClearTiming();
2266  }
2267 
2268  HYPRE_ParCSRPCGGetNumIterations(pcg_solver, &num_iterations);
2269  HYPRE_ParCSRPCGGetFinalRelativeResidualNorm(pcg_solver,
2270  &final_res_norm);
2271 
2272  MPI_Comm_rank(comm, &myid);
2273 
2274  if (myid == 0)
2275  {
2276  mfem::out << "PCG Iterations = " << num_iterations << endl
2277  << "Final PCG Relative Residual Norm = " << final_res_norm
2278  << endl;
2279  }
2280  }
2281  HYPRE_ParCSRPCGSetPrintLevel(pcg_solver, print_level);
2282 }
2283 
2285 {
2286  HYPRE_ParCSRPCGDestroy(pcg_solver);
2287 }
2288 
2289 
2291 {
2292  MPI_Comm comm;
2293 
2294  int k_dim = 50;
2295  int max_iter = 100;
2296  double tol = 1e-6;
2297 
2298  iterative_mode = true;
2299 
2300  HYPRE_ParCSRMatrixGetComm(*A, &comm);
2301 
2302  HYPRE_ParCSRGMRESCreate(comm, &gmres_solver);
2303  HYPRE_ParCSRGMRESSetKDim(gmres_solver, k_dim);
2304  HYPRE_ParCSRGMRESSetMaxIter(gmres_solver, max_iter);
2305  HYPRE_ParCSRGMRESSetTol(gmres_solver, tol);
2306 }
2307 
2308 void HypreGMRES::SetTol(double tol)
2309 {
2310  HYPRE_GMRESSetTol(gmres_solver, tol);
2311 }
2312 
2313 void HypreGMRES::SetMaxIter(int max_iter)
2314 {
2315  HYPRE_GMRESSetMaxIter(gmres_solver, max_iter);
2316 }
2317 
2318 void HypreGMRES::SetKDim(int k_dim)
2319 {
2320  HYPRE_GMRESSetKDim(gmres_solver, k_dim);
2321 }
2322 
2323 void HypreGMRES::SetLogging(int logging)
2324 {
2325  HYPRE_GMRESSetLogging(gmres_solver, logging);
2326 }
2327 
2328 void HypreGMRES::SetPrintLevel(int print_lvl)
2329 {
2330  HYPRE_GMRESSetPrintLevel(gmres_solver, print_lvl);
2331 }
2332 
2334 {
2335  HYPRE_ParCSRGMRESSetPrecond(gmres_solver,
2336  precond.SolveFcn(),
2337  precond.SetupFcn(),
2338  precond);
2339 }
2340 
2342 {
2343  int myid;
2344  HYPRE_Int time_index = 0;
2345  HYPRE_Int num_iterations;
2346  double final_res_norm;
2347  MPI_Comm comm;
2348  HYPRE_Int print_level;
2349 
2350  HYPRE_GMRESGetPrintLevel(gmres_solver, &print_level);
2351 
2352  HYPRE_ParCSRMatrixGetComm(*A, &comm);
2353 
2354  if (!setup_called)
2355  {
2356  if (print_level > 0)
2357  {
2358  time_index = hypre_InitializeTiming("GMRES Setup");
2359  hypre_BeginTiming(time_index);
2360  }
2361 
2362  HYPRE_ParCSRGMRESSetup(gmres_solver, *A, b, x);
2363  setup_called = 1;
2364 
2365  if (print_level > 0)
2366  {
2367  hypre_EndTiming(time_index);
2368  hypre_PrintTiming("Setup phase times", comm);
2369  hypre_FinalizeTiming(time_index);
2370  hypre_ClearTiming();
2371  }
2372  }
2373 
2374  if (print_level > 0)
2375  {
2376  time_index = hypre_InitializeTiming("GMRES Solve");
2377  hypre_BeginTiming(time_index);
2378  }
2379 
2380  if (!iterative_mode)
2381  {
2382  x = 0.0;
2383  }
2384 
2385  HYPRE_ParCSRGMRESSolve(gmres_solver, *A, b, x);
2386 
2387  if (print_level > 0)
2388  {
2389  hypre_EndTiming(time_index);
2390  hypre_PrintTiming("Solve phase times", comm);
2391  hypre_FinalizeTiming(time_index);
2392  hypre_ClearTiming();
2393 
2394  HYPRE_ParCSRGMRESGetNumIterations(gmres_solver, &num_iterations);
2395  HYPRE_ParCSRGMRESGetFinalRelativeResidualNorm(gmres_solver,
2396  &final_res_norm);
2397 
2398  MPI_Comm_rank(comm, &myid);
2399 
2400  if (myid == 0)
2401  {
2402  mfem::out << "GMRES Iterations = " << num_iterations << endl
2403  << "Final GMRES Relative Residual Norm = " << final_res_norm
2404  << endl;
2405  }
2406  }
2407 }
2408 
2410 {
2411  HYPRE_ParCSRGMRESDestroy(gmres_solver);
2412 }
2413 
2414 
2416 {
2417  MPI_Comm comm;
2418 
2419  int sai_max_levels = 1;
2420  double sai_threshold = 0.1;
2421  double sai_filter = 0.1;
2422  int sai_sym = 0;
2423  double sai_loadbal = 0.0;
2424  int sai_reuse = 0;
2425  int sai_logging = 1;
2426 
2427  HYPRE_ParCSRMatrixGetComm(A, &comm);
2428 
2429  HYPRE_ParaSailsCreate(comm, &sai_precond);
2430  HYPRE_ParaSailsSetParams(sai_precond, sai_threshold, sai_max_levels);
2431  HYPRE_ParaSailsSetFilter(sai_precond, sai_filter);
2432  HYPRE_ParaSailsSetSym(sai_precond, sai_sym);
2433  HYPRE_ParaSailsSetLoadbal(sai_precond, sai_loadbal);
2434  HYPRE_ParaSailsSetReuse(sai_precond, sai_reuse);
2435  HYPRE_ParaSailsSetLogging(sai_precond, sai_logging);
2436 }
2437 
2439 {
2440  HYPRE_ParaSailsSetSym(sai_precond, sym);
2441 }
2442 
2444 {
2445  HYPRE_ParaSailsDestroy(sai_precond);
2446 }
2447 
2448 
2450 {
2451  HYPRE_BoomerAMGCreate(&amg_precond);
2452  SetDefaultOptions();
2453 }
2454 
2456 {
2457  HYPRE_BoomerAMGCreate(&amg_precond);
2458  SetDefaultOptions();
2459 }
2460 
2461 void HypreBoomerAMG::SetDefaultOptions()
2462 {
2463  // AMG coarsening options:
2464  int coarsen_type = 10; // 10 = HMIS, 8 = PMIS, 6 = Falgout, 0 = CLJP
2465  int agg_levels = 1; // number of aggressive coarsening levels
2466  double theta = 0.25; // strength threshold: 0.25, 0.5, 0.8
2467 
2468  // AMG interpolation options:
2469  int interp_type = 6; // 6 = extended+i, 0 = classical
2470  int Pmax = 4; // max number of elements per row in P
2471 
2472  // AMG relaxation options:
2473  int relax_type = 8; // 8 = l1-GS, 6 = symm. GS, 3 = GS, 18 = l1-Jacobi
2474  int relax_sweeps = 1; // relaxation sweeps on each level
2475 
2476  // Additional options:
2477  int print_level = 1; // print AMG iterations? 1 = no, 2 = yes
2478  int max_levels = 25; // max number of levels in AMG hierarchy
2479 
2480  HYPRE_BoomerAMGSetCoarsenType(amg_precond, coarsen_type);
2481  HYPRE_BoomerAMGSetAggNumLevels(amg_precond, agg_levels);
2482  HYPRE_BoomerAMGSetRelaxType(amg_precond, relax_type);
2483  HYPRE_BoomerAMGSetNumSweeps(amg_precond, relax_sweeps);
2484  HYPRE_BoomerAMGSetStrongThreshold(amg_precond, theta);
2485  HYPRE_BoomerAMGSetInterpType(amg_precond, interp_type);
2486  HYPRE_BoomerAMGSetPMaxElmts(amg_precond, Pmax);
2487  HYPRE_BoomerAMGSetPrintLevel(amg_precond, print_level);
2488  HYPRE_BoomerAMGSetMaxLevels(amg_precond, max_levels);
2489 
2490  // Use as a preconditioner (one V-cycle, zero tolerance)
2491  HYPRE_BoomerAMGSetMaxIter(amg_precond, 1);
2492  HYPRE_BoomerAMGSetTol(amg_precond, 0.0);
2493 }
2494 
2495 void HypreBoomerAMG::ResetAMGPrecond()
2496 {
2497  HYPRE_Int coarsen_type;
2498  HYPRE_Int agg_levels;
2499  HYPRE_Int relax_type;
2500  HYPRE_Int relax_sweeps;
2501  HYPRE_Real theta;
2502  HYPRE_Int interp_type;
2503  HYPRE_Int Pmax;
2504  HYPRE_Int print_level;
2505  HYPRE_Int dim;
2506  HYPRE_Int nrbms = rbms.Size();
2507  HYPRE_Int nodal;
2508  HYPRE_Int nodal_diag;
2509  HYPRE_Int relax_coarse;
2510  HYPRE_Int interp_vec_variant;
2511  HYPRE_Int q_max;
2512  HYPRE_Int smooth_interp_vectors;
2513  HYPRE_Int interp_refine;
2514 
2515  hypre_ParAMGData *amg_data = (hypre_ParAMGData *)amg_precond;
2516 
2517  // read options from amg_precond
2518  HYPRE_BoomerAMGGetCoarsenType(amg_precond, &coarsen_type);
2519  agg_levels = hypre_ParAMGDataAggNumLevels(amg_data);
2520  relax_type = hypre_ParAMGDataUserRelaxType(amg_data);
2521  relax_sweeps = hypre_ParAMGDataUserNumSweeps(amg_data);
2522  HYPRE_BoomerAMGGetStrongThreshold(amg_precond, &theta);
2523  hypre_BoomerAMGGetInterpType(amg_precond, &interp_type);
2524  HYPRE_BoomerAMGGetPMaxElmts(amg_precond, &Pmax);
2525  HYPRE_BoomerAMGGetPrintLevel(amg_precond, &print_level);
2526  HYPRE_BoomerAMGGetNumFunctions(amg_precond, &dim);
2527  if (nrbms) // elasticity solver options
2528  {
2529  nodal = hypre_ParAMGDataNodal(amg_data);
2530  nodal_diag = hypre_ParAMGDataNodalDiag(amg_data);
2531  HYPRE_BoomerAMGGetCycleRelaxType(amg_precond, &relax_coarse, 3);
2532  interp_vec_variant = hypre_ParAMGInterpVecVariant(amg_data);
2533  q_max = hypre_ParAMGInterpVecQMax(amg_data);
2534  smooth_interp_vectors = hypre_ParAMGSmoothInterpVectors(amg_data);
2535  interp_refine = hypre_ParAMGInterpRefine(amg_data);
2536  }
2537 
2538  HYPRE_BoomerAMGDestroy(amg_precond);
2539  HYPRE_BoomerAMGCreate(&amg_precond);
2540 
2541  HYPRE_BoomerAMGSetCoarsenType(amg_precond, coarsen_type);
2542  HYPRE_BoomerAMGSetAggNumLevels(amg_precond, agg_levels);
2543  HYPRE_BoomerAMGSetRelaxType(amg_precond, relax_type);
2544  HYPRE_BoomerAMGSetNumSweeps(amg_precond, relax_sweeps);
2545  HYPRE_BoomerAMGSetMaxLevels(amg_precond, 25);
2546  HYPRE_BoomerAMGSetTol(amg_precond, 0.0);
2547  HYPRE_BoomerAMGSetMaxIter(amg_precond, 1); // one V-cycle
2548  HYPRE_BoomerAMGSetStrongThreshold(amg_precond, theta);
2549  HYPRE_BoomerAMGSetInterpType(amg_precond, interp_type);
2550  HYPRE_BoomerAMGSetPMaxElmts(amg_precond, Pmax);
2551  HYPRE_BoomerAMGSetPrintLevel(amg_precond, print_level);
2552  HYPRE_BoomerAMGSetNumFunctions(amg_precond, dim);
2553  if (nrbms)
2554  {
2555  HYPRE_BoomerAMGSetNodal(amg_precond, nodal);
2556  HYPRE_BoomerAMGSetNodalDiag(amg_precond, nodal_diag);
2557  HYPRE_BoomerAMGSetCycleRelaxType(amg_precond, relax_coarse, 3);
2558  HYPRE_BoomerAMGSetInterpVecVariant(amg_precond, interp_vec_variant);
2559  HYPRE_BoomerAMGSetInterpVecQMax(amg_precond, q_max);
2560  HYPRE_BoomerAMGSetSmoothInterpVectors(amg_precond, smooth_interp_vectors);
2561  HYPRE_BoomerAMGSetInterpRefine(amg_precond, interp_refine);
2562  RecomputeRBMs();
2563  HYPRE_BoomerAMGSetInterpVectors(amg_precond, rbms.Size(), rbms.GetData());
2564  }
2565 }
2566 
2568 {
2569  const HypreParMatrix *new_A = dynamic_cast<const HypreParMatrix *>(&op);
2570  MFEM_VERIFY(new_A, "new Operator must be a HypreParMatrix!");
2571 
2572  if (A) { ResetAMGPrecond(); }
2573 
2574  // update base classes: Operator, Solver, HypreSolver
2575  height = new_A->Height();
2576  width = new_A->Width();
2577  A = const_cast<HypreParMatrix *>(new_A);
2578  setup_called = 0;
2579  delete X;
2580  delete B;
2581  B = X = NULL;
2582 }
2583 
2585 {
2586  HYPRE_BoomerAMGSetNumFunctions(amg_precond, dim);
2587 
2588  // More robust options with respect to convergence
2589  HYPRE_BoomerAMGSetAggNumLevels(amg_precond, 0);
2590  HYPRE_BoomerAMGSetStrongThreshold(amg_precond, 0.5);
2591 }
2592 
2593 // Rotational rigid-body mode functions, used in SetElasticityOptions()
2594 static void func_rxy(const Vector &x, Vector &y)
2595 {
2596  y = 0.0; y(0) = x(1); y(1) = -x(0);
2597 }
2598 static void func_ryz(const Vector &x, Vector &y)
2599 {
2600  y = 0.0; y(1) = x(2); y(2) = -x(1);
2601 }
2602 static void func_rzx(const Vector &x, Vector &y)
2603 {
2604  y = 0.0; y(2) = x(0); y(0) = -x(2);
2605 }
2606 
2607 void HypreBoomerAMG::RecomputeRBMs()
2608 {
2609  int nrbms;
2610  Array<HypreParVector*> gf_rbms;
2611  int dim = fespace->GetParMesh()->Dimension();
2612 
2613  for (int i = 0; i < rbms.Size(); i++)
2614  {
2615  HYPRE_ParVectorDestroy(rbms[i]);
2616  }
2617 
2618  if (dim == 2)
2619  {
2620  nrbms = 1;
2621 
2622  VectorFunctionCoefficient coeff_rxy(2, func_rxy);
2623 
2624  ParGridFunction rbms_rxy(fespace);
2625  rbms_rxy.ProjectCoefficient(coeff_rxy);
2626 
2627  rbms.SetSize(nrbms);
2628  gf_rbms.SetSize(nrbms);
2629  gf_rbms[0] = rbms_rxy.ParallelAverage();
2630  }
2631  else if (dim == 3)
2632  {
2633  nrbms = 3;
2634 
2635  VectorFunctionCoefficient coeff_rxy(3, func_rxy);
2636  VectorFunctionCoefficient coeff_ryz(3, func_ryz);
2637  VectorFunctionCoefficient coeff_rzx(3, func_rzx);
2638 
2639  ParGridFunction rbms_rxy(fespace);
2640  ParGridFunction rbms_ryz(fespace);
2641  ParGridFunction rbms_rzx(fespace);
2642  rbms_rxy.ProjectCoefficient(coeff_rxy);
2643  rbms_ryz.ProjectCoefficient(coeff_ryz);
2644  rbms_rzx.ProjectCoefficient(coeff_rzx);
2645 
2646  rbms.SetSize(nrbms);
2647  gf_rbms.SetSize(nrbms);
2648  gf_rbms[0] = rbms_rxy.ParallelAverage();
2649  gf_rbms[1] = rbms_ryz.ParallelAverage();
2650  gf_rbms[2] = rbms_rzx.ParallelAverage();
2651  }
2652  else
2653  {
2654  nrbms = 0;
2655  rbms.SetSize(nrbms);
2656  }
2657 
2658  // Transfer the RBMs from the ParGridFunction to the HYPRE_ParVector objects
2659  for (int i = 0; i < nrbms; i++)
2660  {
2661  rbms[i] = gf_rbms[i]->StealParVector();
2662  delete gf_rbms[i];
2663  }
2664 }
2665 
2667 {
2668  // Save the finite element space to support multiple calls to SetOperator()
2669  this->fespace = fespace;
2670 
2671  // Make sure the systems AMG options are set
2672  int dim = fespace->GetParMesh()->Dimension();
2674 
2675  // Nodal coarsening options (nodal coarsening is required for this solver)
2676  // See hypre's new_ij driver and the paper for descriptions.
2677  int nodal = 4; // strength reduction norm: 1, 3 or 4
2678  int nodal_diag = 1; // diagonal in strength matrix: 0, 1 or 2
2679  int relax_coarse = 8; // smoother on the coarsest grid: 8, 99 or 29
2680 
2681  // Elasticity interpolation options
2682  int interp_vec_variant = 2; // 1 = GM-1, 2 = GM-2, 3 = LN
2683  int q_max = 4; // max elements per row for each Q
2684  int smooth_interp_vectors = 1; // smooth the rigid-body modes?
2685 
2686  // Optionally pre-process the interpolation matrix through iterative weight
2687  // refinement (this is generally applicable for any system)
2688  int interp_refine = 1;
2689 
2690  HYPRE_BoomerAMGSetNodal(amg_precond, nodal);
2691  HYPRE_BoomerAMGSetNodalDiag(amg_precond, nodal_diag);
2692  HYPRE_BoomerAMGSetCycleRelaxType(amg_precond, relax_coarse, 3);
2693  HYPRE_BoomerAMGSetInterpVecVariant(amg_precond, interp_vec_variant);
2694  HYPRE_BoomerAMGSetInterpVecQMax(amg_precond, q_max);
2695  HYPRE_BoomerAMGSetSmoothInterpVectors(amg_precond, smooth_interp_vectors);
2696  HYPRE_BoomerAMGSetInterpRefine(amg_precond, interp_refine);
2697 
2698  RecomputeRBMs();
2699  HYPRE_BoomerAMGSetInterpVectors(amg_precond, rbms.Size(), rbms.GetData());
2700 }
2701 
2703 {
2704  for (int i = 0; i < rbms.Size(); i++)
2705  {
2706  HYPRE_ParVectorDestroy(rbms[i]);
2707  }
2708 
2709  HYPRE_BoomerAMGDestroy(amg_precond);
2710 }
2711 
2712 
2714  : HypreSolver(&A)
2715 {
2716  int cycle_type = 13;
2717  int rlx_type = 2;
2718  int rlx_sweeps = 1;
2719  double rlx_weight = 1.0;
2720  double rlx_omega = 1.0;
2721  int amg_coarsen_type = 10;
2722  int amg_agg_levels = 1;
2723  int amg_rlx_type = 8;
2724  double theta = 0.25;
2725  int amg_interp_type = 6;
2726  int amg_Pmax = 4;
2727 
2728  int dim = edge_fespace->GetMesh()->Dimension();
2729  int sdim = edge_fespace->GetMesh()->SpaceDimension();
2730  const FiniteElementCollection *edge_fec = edge_fespace->FEColl();
2731 
2732  bool trace_space, rt_trace_space;
2733  ND_Trace_FECollection *nd_tr_fec = NULL;
2734  trace_space = dynamic_cast<const ND_Trace_FECollection*>(edge_fec);
2735  rt_trace_space = dynamic_cast<const RT_Trace_FECollection*>(edge_fec);
2736  trace_space = trace_space || rt_trace_space;
2737 
2738  int p = 1;
2739  if (edge_fespace->GetNE() > 0)
2740  {
2741  if (trace_space)
2742  {
2743  p = edge_fespace->GetFaceOrder(0);
2744  if (dim == 2) { p++; }
2745  }
2746  else
2747  {
2748  p = edge_fespace->GetOrder(0);
2749  }
2750  }
2751 
2752  ParMesh *pmesh = edge_fespace->GetParMesh();
2753  if (rt_trace_space)
2754  {
2755  nd_tr_fec = new ND_Trace_FECollection(p, dim);
2756  edge_fespace = new ParFiniteElementSpace(pmesh, nd_tr_fec);
2757  }
2758 
2759  HYPRE_AMSCreate(&ams);
2760 
2761  HYPRE_AMSSetDimension(ams, sdim); // 2D H(div) and 3D H(curl) problems
2762  HYPRE_AMSSetTol(ams, 0.0);
2763  HYPRE_AMSSetMaxIter(ams, 1); // use as a preconditioner
2764  HYPRE_AMSSetCycleType(ams, cycle_type);
2765  HYPRE_AMSSetPrintLevel(ams, 1);
2766 
2767  // define the nodal linear finite element space associated with edge_fespace
2768  FiniteElementCollection *vert_fec;
2769  if (trace_space)
2770  {
2771  vert_fec = new H1_Trace_FECollection(p, dim);
2772  }
2773  else
2774  {
2775  vert_fec = new H1_FECollection(p, dim);
2776  }
2777  ParFiniteElementSpace *vert_fespace = new ParFiniteElementSpace(pmesh,
2778  vert_fec);
2779 
2780  // generate and set the vertex coordinates
2781  if (p == 1)
2782  {
2783  ParGridFunction x_coord(vert_fespace);
2784  ParGridFunction y_coord(vert_fespace);
2785  ParGridFunction z_coord(vert_fespace);
2786  double *coord;
2787  for (int i = 0; i < pmesh->GetNV(); i++)
2788  {
2789  coord = pmesh -> GetVertex(i);
2790  x_coord(i) = coord[0];
2791  y_coord(i) = coord[1];
2792  if (sdim == 3) { z_coord(i) = coord[2]; }
2793  }
2794  x = x_coord.ParallelProject();
2795  y = y_coord.ParallelProject();
2796  if (sdim == 2)
2797  {
2798  z = NULL;
2799  HYPRE_AMSSetCoordinateVectors(ams, *x, *y, NULL);
2800  }
2801  else
2802  {
2803  z = z_coord.ParallelProject();
2804  HYPRE_AMSSetCoordinateVectors(ams, *x, *y, *z);
2805  }
2806  }
2807  else
2808  {
2809  x = NULL;
2810  y = NULL;
2811  z = NULL;
2812  }
2813 
2814  // generate and set the discrete gradient
2816  grad = new ParDiscreteLinearOperator(vert_fespace, edge_fespace);
2817  if (trace_space)
2818  {
2820  }
2821  else
2822  {
2824  }
2825  grad->Assemble();
2826  grad->Finalize();
2827  G = grad->ParallelAssemble();
2828  HYPRE_AMSSetDiscreteGradient(ams, *G);
2829  delete grad;
2830 
2831  // generate and set the Nedelec interpolation matrices
2832  Pi = Pix = Piy = Piz = NULL;
2833  if (p > 1)
2834  {
2835  ParFiniteElementSpace *vert_fespace_d
2836  = new ParFiniteElementSpace(pmesh, vert_fec, sdim, Ordering::byVDIM);
2837 
2839  id_ND = new ParDiscreteLinearOperator(vert_fespace_d, edge_fespace);
2840  if (trace_space)
2841  {
2843  }
2844  else
2845  {
2847  }
2848  id_ND->Assemble();
2849  id_ND->Finalize();
2850 
2851  if (cycle_type < 10)
2852  {
2853  Pi = id_ND->ParallelAssemble();
2854  }
2855  else
2856  {
2857  Array2D<HypreParMatrix *> Pi_blocks;
2858  id_ND->GetParBlocks(Pi_blocks);
2859  Pix = Pi_blocks(0,0);
2860  Piy = Pi_blocks(0,1);
2861  if (sdim == 3) { Piz = Pi_blocks(0,2); }
2862  }
2863 
2864  delete id_ND;
2865 
2866  HYPRE_ParCSRMatrix HY_Pi = (Pi) ? (HYPRE_ParCSRMatrix) *Pi : NULL;
2867  HYPRE_ParCSRMatrix HY_Pix = (Pix) ? (HYPRE_ParCSRMatrix) *Pix : NULL;
2868  HYPRE_ParCSRMatrix HY_Piy = (Piy) ? (HYPRE_ParCSRMatrix) *Piy : NULL;
2869  HYPRE_ParCSRMatrix HY_Piz = (Piz) ? (HYPRE_ParCSRMatrix) *Piz : NULL;
2870  HYPRE_AMSSetInterpolations(ams, HY_Pi, HY_Pix, HY_Piy, HY_Piz);
2871 
2872  delete vert_fespace_d;
2873  }
2874 
2875  delete vert_fespace;
2876  delete vert_fec;
2877 
2878  if (rt_trace_space)
2879  {
2880  delete edge_fespace;
2881  delete nd_tr_fec;
2882  }
2883 
2884  // set additional AMS options
2885  HYPRE_AMSSetSmoothingOptions(ams, rlx_type, rlx_sweeps, rlx_weight, rlx_omega);
2886  HYPRE_AMSSetAlphaAMGOptions(ams, amg_coarsen_type, amg_agg_levels, amg_rlx_type,
2887  theta, amg_interp_type, amg_Pmax);
2888  HYPRE_AMSSetBetaAMGOptions(ams, amg_coarsen_type, amg_agg_levels, amg_rlx_type,
2889  theta, amg_interp_type, amg_Pmax);
2890 }
2891 
2893 {
2894  HYPRE_AMSDestroy(ams);
2895 
2896  delete x;
2897  delete y;
2898  delete z;
2899 
2900  delete G;
2901  delete Pi;
2902  delete Pix;
2903  delete Piy;
2904  delete Piz;
2905 }
2906 
2907 void HypreAMS::SetPrintLevel(int print_lvl)
2908 {
2909  HYPRE_AMSSetPrintLevel(ams, print_lvl);
2910 }
2911 
2913  : HypreSolver(&A)
2914 {
2915  int cycle_type = 11;
2916  int rlx_type = 2;
2917  int rlx_sweeps = 1;
2918  double rlx_weight = 1.0;
2919  double rlx_omega = 1.0;
2920  int amg_coarsen_type = 10;
2921  int amg_agg_levels = 1;
2922  int amg_rlx_type = 8;
2923  double theta = 0.25;
2924  int amg_interp_type = 6;
2925  int amg_Pmax = 4;
2926  int ams_cycle_type = 14;
2927 
2928  const FiniteElementCollection *face_fec = face_fespace->FEColl();
2929  bool trace_space =
2930  (dynamic_cast<const RT_Trace_FECollection*>(face_fec) != NULL);
2931  int p = 1;
2932  if (face_fespace->GetNE() > 0)
2933  {
2934  if (trace_space)
2935  {
2936  p = face_fespace->GetFaceOrder(0) + 1;
2937  }
2938  else
2939  {
2940  p = face_fespace->GetOrder(0);
2941  }
2942  }
2943 
2944  HYPRE_ADSCreate(&ads);
2945 
2946  HYPRE_ADSSetTol(ads, 0.0);
2947  HYPRE_ADSSetMaxIter(ads, 1); // use as a preconditioner
2948  HYPRE_ADSSetCycleType(ads, cycle_type);
2949  HYPRE_ADSSetPrintLevel(ads, 1);
2950 
2951  // define the nodal and edge finite element spaces associated with face_fespace
2952  ParMesh *pmesh = (ParMesh *) face_fespace->GetMesh();
2953  FiniteElementCollection *vert_fec, *edge_fec;
2954  if (trace_space)
2955  {
2956  vert_fec = new H1_Trace_FECollection(p, 3);
2957  edge_fec = new ND_Trace_FECollection(p, 3);
2958  }
2959  else
2960  {
2961  vert_fec = new H1_FECollection(p, 3);
2962  edge_fec = new ND_FECollection(p, 3);
2963  }
2964 
2965  ParFiniteElementSpace *vert_fespace = new ParFiniteElementSpace(pmesh,
2966  vert_fec);
2967  ParFiniteElementSpace *edge_fespace = new ParFiniteElementSpace(pmesh,
2968  edge_fec);
2969 
2970  // generate and set the vertex coordinates
2971  if (p == 1)
2972  {
2973  ParGridFunction x_coord(vert_fespace);
2974  ParGridFunction y_coord(vert_fespace);
2975  ParGridFunction z_coord(vert_fespace);
2976  double *coord;
2977  for (int i = 0; i < pmesh->GetNV(); i++)
2978  {
2979  coord = pmesh -> GetVertex(i);
2980  x_coord(i) = coord[0];
2981  y_coord(i) = coord[1];
2982  z_coord(i) = coord[2];
2983  }
2984  x = x_coord.ParallelProject();
2985  y = y_coord.ParallelProject();
2986  z = z_coord.ParallelProject();
2987  HYPRE_ADSSetCoordinateVectors(ads, *x, *y, *z);
2988  }
2989  else
2990  {
2991  x = NULL;
2992  y = NULL;
2993  z = NULL;
2994  }
2995 
2996  // generate and set the discrete curl
2998  curl = new ParDiscreteLinearOperator(edge_fespace, face_fespace);
2999  if (trace_space)
3000  {
3002  }
3003  else
3004  {
3006  }
3007  curl->Assemble();
3008  curl->Finalize();
3009  C = curl->ParallelAssemble();
3010  C->CopyColStarts(); // since we'll delete edge_fespace
3011  HYPRE_ADSSetDiscreteCurl(ads, *C);
3012  delete curl;
3013 
3014  // generate and set the discrete gradient
3016  grad = new ParDiscreteLinearOperator(vert_fespace, edge_fespace);
3017  if (trace_space)
3018  {
3020  }
3021  else
3022  {
3024  }
3025  grad->Assemble();
3026  grad->Finalize();
3027  G = grad->ParallelAssemble();
3028  G->CopyColStarts(); // since we'll delete vert_fespace
3029  G->CopyRowStarts(); // since we'll delete edge_fespace
3030  HYPRE_ADSSetDiscreteGradient(ads, *G);
3031  delete grad;
3032 
3033  // generate and set the Nedelec and Raviart-Thomas interpolation matrices
3034  RT_Pi = RT_Pix = RT_Piy = RT_Piz = NULL;
3035  ND_Pi = ND_Pix = ND_Piy = ND_Piz = NULL;
3036  if (p > 1)
3037  {
3038  ParFiniteElementSpace *vert_fespace_d
3039  = new ParFiniteElementSpace(pmesh, vert_fec, 3, Ordering::byVDIM);
3040 
3042  id_ND = new ParDiscreteLinearOperator(vert_fespace_d, edge_fespace);
3043  if (trace_space)
3044  {
3046  }
3047  else
3048  {
3050  }
3051  id_ND->Assemble();
3052  id_ND->Finalize();
3053 
3054  if (ams_cycle_type < 10)
3055  {
3056  ND_Pi = id_ND->ParallelAssemble();
3057  ND_Pi->CopyColStarts(); // since we'll delete vert_fespace_d
3058  ND_Pi->CopyRowStarts(); // since we'll delete edge_fespace
3059  }
3060  else
3061  {
3062  Array2D<HypreParMatrix *> ND_Pi_blocks;
3063  id_ND->GetParBlocks(ND_Pi_blocks);
3064  ND_Pix = ND_Pi_blocks(0,0);
3065  ND_Piy = ND_Pi_blocks(0,1);
3066  ND_Piz = ND_Pi_blocks(0,2);
3067  }
3068 
3069  delete id_ND;
3070 
3072  id_RT = new ParDiscreteLinearOperator(vert_fespace_d, face_fespace);
3073  if (trace_space)
3074  {
3076  }
3077  else
3078  {
3080  }
3081  id_RT->Assemble();
3082  id_RT->Finalize();
3083 
3084  if (cycle_type < 10)
3085  {
3086  RT_Pi = id_RT->ParallelAssemble();
3087  RT_Pi->CopyColStarts(); // since we'll delete vert_fespace_d
3088  }
3089  else
3090  {
3091  Array2D<HypreParMatrix *> RT_Pi_blocks;
3092  id_RT->GetParBlocks(RT_Pi_blocks);
3093  RT_Pix = RT_Pi_blocks(0,0);
3094  RT_Piy = RT_Pi_blocks(0,1);
3095  RT_Piz = RT_Pi_blocks(0,2);
3096  }
3097 
3098  delete id_RT;
3099 
3100  HYPRE_ParCSRMatrix HY_RT_Pi, HY_RT_Pix, HY_RT_Piy, HY_RT_Piz;
3101  HY_RT_Pi = (RT_Pi) ? (HYPRE_ParCSRMatrix) *RT_Pi : NULL;
3102  HY_RT_Pix = (RT_Pix) ? (HYPRE_ParCSRMatrix) *RT_Pix : NULL;
3103  HY_RT_Piy = (RT_Piy) ? (HYPRE_ParCSRMatrix) *RT_Piy : NULL;
3104  HY_RT_Piz = (RT_Piz) ? (HYPRE_ParCSRMatrix) *RT_Piz : NULL;
3105  HYPRE_ParCSRMatrix HY_ND_Pi, HY_ND_Pix, HY_ND_Piy, HY_ND_Piz;
3106  HY_ND_Pi = (ND_Pi) ? (HYPRE_ParCSRMatrix) *ND_Pi : NULL;
3107  HY_ND_Pix = (ND_Pix) ? (HYPRE_ParCSRMatrix) *ND_Pix : NULL;
3108  HY_ND_Piy = (ND_Piy) ? (HYPRE_ParCSRMatrix) *ND_Piy : NULL;
3109  HY_ND_Piz = (ND_Piz) ? (HYPRE_ParCSRMatrix) *ND_Piz : NULL;
3110  HYPRE_ADSSetInterpolations(ads,
3111  HY_RT_Pi, HY_RT_Pix, HY_RT_Piy, HY_RT_Piz,
3112  HY_ND_Pi, HY_ND_Pix, HY_ND_Piy, HY_ND_Piz);
3113 
3114  delete vert_fespace_d;
3115  }
3116 
3117  delete vert_fec;
3118  delete vert_fespace;
3119  delete edge_fec;
3120  delete edge_fespace;
3121 
3122  // set additional ADS options
3123  HYPRE_ADSSetSmoothingOptions(ads, rlx_type, rlx_sweeps, rlx_weight, rlx_omega);
3124  HYPRE_ADSSetAMGOptions(ads, amg_coarsen_type, amg_agg_levels, amg_rlx_type,
3125  theta, amg_interp_type, amg_Pmax);
3126  HYPRE_ADSSetAMSOptions(ads, ams_cycle_type, amg_coarsen_type, amg_agg_levels,
3127  amg_rlx_type, theta, amg_interp_type, amg_Pmax);
3128 }
3129 
3131 {
3132  HYPRE_ADSDestroy(ads);
3133 
3134  delete x;
3135  delete y;
3136  delete z;
3137 
3138  delete G;
3139  delete C;
3140 
3141  delete RT_Pi;
3142  delete RT_Pix;
3143  delete RT_Piy;
3144  delete RT_Piz;
3145 
3146  delete ND_Pi;
3147  delete ND_Pix;
3148  delete ND_Piy;
3149  delete ND_Piz;
3150 }
3151 
3152 void HypreADS::SetPrintLevel(int print_lvl)
3153 {
3154  HYPRE_ADSSetPrintLevel(ads, print_lvl);
3155 }
3156 
3157 HypreLOBPCG::HypreMultiVector::HypreMultiVector(int n, HypreParVector & v,
3158  mv_InterfaceInterpreter & interpreter)
3159  : hpv(NULL),
3160  nv(n)
3161 {
3162  mv_ptr = mv_MultiVectorCreateFromSampleVector(&interpreter, nv,
3163  (HYPRE_ParVector)v);
3164 
3165  HYPRE_ParVector* vecs = NULL;
3166  {
3167  mv_TempMultiVector* tmp =
3168  (mv_TempMultiVector*)mv_MultiVectorGetData(mv_ptr);
3169  vecs = (HYPRE_ParVector*)(tmp -> vector);
3170  }
3171 
3172  hpv = new HypreParVector*[nv];
3173  for (int i=0; i<nv; i++)
3174  {
3175  hpv[i] = new HypreParVector(vecs[i]);
3176  }
3177 }
3178 
3179 HypreLOBPCG::HypreMultiVector::~HypreMultiVector()
3180 {
3181  if ( hpv != NULL )
3182  {
3183  for (int i=0; i<nv; i++)
3184  {
3185  delete hpv[i];
3186  }
3187  delete [] hpv;
3188  }
3189 
3190  mv_MultiVectorDestroy(mv_ptr);
3191 }
3192 
3193 void
3194 HypreLOBPCG::HypreMultiVector::Randomize(HYPRE_Int seed)
3195 {
3196  mv_MultiVectorSetRandom(mv_ptr, seed);
3197 }
3198 
3199 HypreParVector &
3200 HypreLOBPCG::HypreMultiVector::GetVector(unsigned int i)
3201 {
3202  MFEM_ASSERT((int)i < nv, "index out of range");
3203 
3204  return ( *hpv[i] );
3205 }
3206 
3207 HypreParVector **
3208 HypreLOBPCG::HypreMultiVector::StealVectors()
3209 {
3210  HypreParVector ** hpv_ret = hpv;
3211 
3212  hpv = NULL;
3213 
3214  mv_TempMultiVector * mv_tmp =
3215  (mv_TempMultiVector*)mv_MultiVectorGetData(mv_ptr);
3216 
3217  mv_tmp->ownsVectors = 0;
3218 
3219  for (int i=0; i<nv; i++)
3220  {
3221  hpv_ret[i]->SetOwnership(1);
3222  }
3223 
3224  return hpv_ret;
3225 }
3226 
3228  : comm(c),
3229  myid(0),
3230  numProcs(1),
3231  nev(10),
3232  seed(75),
3233  glbSize(-1),
3234  part(NULL),
3235  multi_vec(NULL),
3236  x(NULL),
3237  subSpaceProj(NULL)
3238 {
3239  MPI_Comm_size(comm,&numProcs);
3240  MPI_Comm_rank(comm,&myid);
3241 
3242  HYPRE_ParCSRSetupInterpreter(&interpreter);
3243  HYPRE_ParCSRSetupMatvec(&matvec_fn);
3244  HYPRE_LOBPCGCreate(&interpreter, &matvec_fn, &lobpcg_solver);
3245 }
3246 
3248 {
3249  delete multi_vec;
3250  delete x;
3251  delete [] part;
3252 
3253  HYPRE_LOBPCGDestroy(lobpcg_solver);
3254 }
3255 
3256 void
3258 {
3259  HYPRE_LOBPCGSetTol(lobpcg_solver, tol);
3260 }
3261 
3262 void
3263 HypreLOBPCG::SetRelTol(double rel_tol)
3264 {
3265 #if MFEM_HYPRE_VERSION >= 21101
3266  HYPRE_LOBPCGSetRTol(lobpcg_solver, rel_tol);
3267 #else
3268  MFEM_ABORT("This method requires HYPRE version >= 2.11.1");
3269 #endif
3270 }
3271 
3272 void
3274 {
3275  HYPRE_LOBPCGSetMaxIter(lobpcg_solver, max_iter);
3276 }
3277 
3278 void
3280 {
3281  if (myid == 0)
3282  {
3283  HYPRE_LOBPCGSetPrintLevel(lobpcg_solver, logging);
3284  }
3285 }
3286 
3287 void
3289 {
3290  HYPRE_LOBPCGSetPrecondUsageMode(lobpcg_solver, pcg_mode);
3291 }
3292 
3293 void
3295 {
3296  HYPRE_LOBPCGSetPrecond(lobpcg_solver,
3297  (HYPRE_PtrToSolverFcn)this->PrecondSolve,
3298  (HYPRE_PtrToSolverFcn)this->PrecondSetup,
3299  (HYPRE_Solver)&precond);
3300 }
3301 
3302 void
3304 {
3305  int locSize = A.Width();
3306 
3307  if (HYPRE_AssumedPartitionCheck())
3308  {
3309  part = new HYPRE_Int[2];
3310 
3311  MPI_Scan(&locSize, &part[1], 1, HYPRE_MPI_INT, MPI_SUM, comm);
3312 
3313  part[0] = part[1] - locSize;
3314 
3315  MPI_Allreduce(&locSize, &glbSize, 1, HYPRE_MPI_INT, MPI_SUM, comm);
3316  }
3317  else
3318  {
3319  part = new HYPRE_Int[numProcs+1];
3320 
3321  MPI_Allgather(&locSize, 1, MPI_INT,
3322  &part[1], 1, HYPRE_MPI_INT, comm);
3323 
3324  part[0] = 0;
3325  for (int i=0; i<numProcs; i++)
3326  {
3327  part[i+1] += part[i];
3328  }
3329 
3330  glbSize = part[numProcs];
3331  }
3332 
3333  if ( x != NULL )
3334  {
3335  delete x;
3336  }
3337 
3338  // Create a distributed vector without a data array.
3339  x = new HypreParVector(comm,glbSize,NULL,part);
3340 
3341  matvec_fn.MatvecCreate = this->OperatorMatvecCreate;
3342  matvec_fn.Matvec = this->OperatorMatvec;
3343  matvec_fn.MatvecDestroy = this->OperatorMatvecDestroy;
3344 
3345  HYPRE_LOBPCGSetup(lobpcg_solver,(HYPRE_Matrix)&A,NULL,NULL);
3346 }
3347 
3348 void
3350 {
3351  matvec_fn.MatvecCreate = this->OperatorMatvecCreate;
3352  matvec_fn.Matvec = this->OperatorMatvec;
3353  matvec_fn.MatvecDestroy = this->OperatorMatvecDestroy;
3354 
3355  HYPRE_LOBPCGSetupB(lobpcg_solver,(HYPRE_Matrix)&M,NULL);
3356 }
3357 
3358 void
3360 {
3361  // Initialize eigenvalues array with marker values
3362  eigs.SetSize(nev);
3363 
3364  for (int i=0; i<nev; i++)
3365  {
3366  eigs[i] = eigenvalues[i];
3367  }
3368 }
3369 
3372 {
3373  return multi_vec->GetVector(i);
3374 }
3375 
3376 void
3378 {
3379  // Initialize HypreMultiVector object if necessary
3380  if ( multi_vec == NULL )
3381  {
3382  MFEM_ASSERT(x != NULL, "In HypreLOBPCG::SetInitialVectors()");
3383 
3384  multi_vec = new HypreMultiVector(nev, *x, interpreter);
3385  }
3386 
3387  // Copy the vectors provided
3388  for (int i=0; i < min(num_vecs,nev); i++)
3389  {
3390  multi_vec->GetVector(i) = *vecs[i];
3391  }
3392 
3393  // Randomize any remaining vectors
3394  for (int i=min(num_vecs,nev); i < nev; i++)
3395  {
3396  multi_vec->GetVector(i).Randomize(seed);
3397  }
3398 
3399  // Ensure all vectors are in the proper subspace
3400  if ( subSpaceProj != NULL )
3401  {
3402  HypreParVector y(*x);
3403  y = multi_vec->GetVector(0);
3404 
3405  for (int i=1; i<nev; i++)
3406  {
3407  subSpaceProj->Mult(multi_vec->GetVector(i),
3408  multi_vec->GetVector(i-1));
3409  }
3410  subSpaceProj->Mult(y,
3411  multi_vec->GetVector(nev-1));
3412  }
3413 }
3414 
3415 void
3417 {
3418  // Initialize HypreMultiVector object if necessary
3419  if ( multi_vec == NULL )
3420  {
3421  MFEM_ASSERT(x != NULL, "In HypreLOBPCG::Solve()");
3422 
3423  multi_vec = new HypreMultiVector(nev, *x, interpreter);
3424  multi_vec->Randomize(seed);
3425 
3426  if ( subSpaceProj != NULL )
3427  {
3428  HypreParVector y(*x);
3429  y = multi_vec->GetVector(0);
3430 
3431  for (int i=1; i<nev; i++)
3432  {
3433  subSpaceProj->Mult(multi_vec->GetVector(i),
3434  multi_vec->GetVector(i-1));
3435  }
3436  subSpaceProj->Mult(y, multi_vec->GetVector(nev-1));
3437  }
3438  }
3439 
3440  eigenvalues.SetSize(nev);
3441  eigenvalues = NAN;
3442 
3443  // Perform eigenmode calculation
3444  //
3445  // The eigenvalues are computed in ascending order (internally the
3446  // order is determined by the LAPACK routine 'dsydv'.)
3447  HYPRE_LOBPCGSolve(lobpcg_solver, NULL, *multi_vec, eigenvalues);
3448 }
3449 
3450 void *
3451 HypreLOBPCG::OperatorMatvecCreate( void *A,
3452  void *x )
3453 {
3454  void *matvec_data;
3455 
3456  matvec_data = NULL;
3457 
3458  return ( matvec_data );
3459 }
3460 
3461 HYPRE_Int
3462 HypreLOBPCG::OperatorMatvec( void *matvec_data,
3463  HYPRE_Complex alpha,
3464  void *A,
3465  void *x,
3466  HYPRE_Complex beta,
3467  void *y )
3468 {
3469  MFEM_VERIFY(alpha == 1.0 && beta == 0.0, "values not supported");
3470 
3471  Operator *Aop = (Operator*)A;
3472 
3473  int width = Aop->Width();
3474 
3475  hypre_ParVector * xPar = (hypre_ParVector *)x;
3476  hypre_ParVector * yPar = (hypre_ParVector *)y;
3477 
3478  Vector xVec(xPar->local_vector->data, width);
3479  Vector yVec(yPar->local_vector->data, width);
3480 
3481  Aop->Mult( xVec, yVec );
3482 
3483  return 0;
3484 }
3485 
3486 HYPRE_Int
3487 HypreLOBPCG::OperatorMatvecDestroy( void *matvec_data )
3488 {
3489  return 0;
3490 }
3491 
3492 HYPRE_Int
3493 HypreLOBPCG::PrecondSolve(void *solver,
3494  void *A,
3495  void *b,
3496  void *x)
3497 {
3498  Solver *PC = (Solver*)solver;
3499  Operator *OP = (Operator*)A;
3500 
3501  int width = OP->Width();
3502 
3503  hypre_ParVector * bPar = (hypre_ParVector *)b;
3504  hypre_ParVector * xPar = (hypre_ParVector *)x;
3505 
3506  Vector bVec(bPar->local_vector->data, width);
3507  Vector xVec(xPar->local_vector->data, width);
3508 
3509  PC->Mult( bVec, xVec );
3510 
3511  return 0;
3512 }
3513 
3514 HYPRE_Int
3515 HypreLOBPCG::PrecondSetup(void *solver,
3516  void *A,
3517  void *b,
3518  void *x)
3519 {
3520  return 0;
3521 }
3522 
3523 HypreAME::HypreAME(MPI_Comm comm)
3524  : myid(0),
3525  numProcs(1),
3526  nev(10),
3527  setT(false),
3528  ams_precond(NULL),
3529  eigenvalues(NULL),
3530  multi_vec(NULL),
3531  eigenvectors(NULL)
3532 {
3533  MPI_Comm_size(comm,&numProcs);
3534  MPI_Comm_rank(comm,&myid);
3535 
3536  HYPRE_AMECreate(&ame_solver);
3537  HYPRE_AMESetPrintLevel(ame_solver, 0);
3538 }
3539 
3541 {
3542  if ( multi_vec )
3543  {
3544  mfem_hypre_TFree(multi_vec);
3545  }
3546 
3547  if ( eigenvectors )
3548  {
3549  for (int i=0; i<nev; i++)
3550  {
3551  delete eigenvectors[i];
3552  }
3553  }
3554  delete [] eigenvectors;
3555 
3556  if ( eigenvalues )
3557  {
3558  mfem_hypre_TFree(eigenvalues);
3559  }
3560 
3561  HYPRE_AMEDestroy(ame_solver);
3562 }
3563 
3564 void
3566 {
3567  nev = num_eigs;
3568 
3569  HYPRE_AMESetBlockSize(ame_solver, nev);
3570 }
3571 
3572 void
3573 HypreAME::SetTol(double tol)
3574 {
3575  HYPRE_AMESetTol(ame_solver, tol);
3576 }
3577 
3578 void
3579 HypreAME::SetRelTol(double rel_tol)
3580 {
3581 #if MFEM_HYPRE_VERSION >= 21101
3582  HYPRE_AMESetRTol(ame_solver, rel_tol);
3583 #else
3584  MFEM_ABORT("This method requires HYPRE version >= 2.11.1");
3585 #endif
3586 }
3587 
3588 void
3590 {
3591  HYPRE_AMESetMaxIter(ame_solver, max_iter);
3592 }
3593 
3594 void
3596 {
3597  if (myid == 0)
3598  {
3599  HYPRE_AMESetPrintLevel(ame_solver, logging);
3600  }
3601 }
3602 
3603 void
3605 {
3606  ams_precond = &precond;
3607 }
3608 
3609 void
3611 {
3612  if ( !setT )
3613  {
3614  HYPRE_Solver ams_precond_ptr = (HYPRE_Solver)*ams_precond;
3615 
3616  ams_precond->SetupFcn()(*ams_precond,A,NULL,NULL);
3617 
3618  HYPRE_AMESetAMSSolver(ame_solver, ams_precond_ptr);
3619  }
3620 
3621  HYPRE_AMESetup(ame_solver);
3622 }
3623 
3624 void
3626 {
3627  HYPRE_ParCSRMatrix parcsr_M = M;
3628  HYPRE_AMESetMassMatrix(ame_solver,(HYPRE_ParCSRMatrix)parcsr_M);
3629 }
3630 
3631 void
3633 {
3634  HYPRE_AMESolve(ame_solver);
3635 }
3636 
3637 void
3639 {
3640  // Initialize eigenvalues array with marker values
3641  eigs.SetSize(nev); eigs = -1.0;
3642 
3643  if ( eigenvalues == NULL )
3644  {
3645  // Grab eigenvalues from AME
3646  HYPRE_AMEGetEigenvalues(ame_solver,&eigenvalues);
3647  }
3648 
3649  // Copy eigenvalues to eigs array
3650  for (int i=0; i<nev; i++)
3651  {
3652  eigs[i] = eigenvalues[i];
3653  }
3654 }
3655 
3656 void
3657 HypreAME::createDummyVectors()
3658 {
3659  if ( multi_vec == NULL )
3660  {
3661  HYPRE_AMEGetEigenvectors(ame_solver,&multi_vec);
3662  }
3663 
3664  eigenvectors = new HypreParVector*[nev];
3665  for (int i=0; i<nev; i++)
3666  {
3667  eigenvectors[i] = new HypreParVector(multi_vec[i]);
3668  eigenvectors[i]->SetOwnership(1);
3669  }
3670 
3671 }
3672 
3673 HypreParVector &
3675 {
3676  if ( eigenvectors == NULL )
3677  {
3678  this->createDummyVectors();
3679  }
3680 
3681  return *eigenvectors[i];
3682 }
3683 
3684 HypreParVector **
3686 {
3687  if ( eigenvectors == NULL )
3688  {
3689  this->createDummyVectors();
3690  }
3691 
3692  // Set the local pointers to NULL so that they won't be deleted later
3693  HypreParVector ** vecs = eigenvectors;
3694  eigenvectors = NULL;
3695  multi_vec = NULL;
3696 
3697  return vecs;
3698 }
3699 
3700 }
3701 
3702 #endif
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
virtual int NumNonZeroElems() const
Returns the number of the nonzero elements in the matrix.
Definition: sparsemat.cpp:983
Vector()
Default constructor for Vector. Sets size = 0 and data = NULL.
Definition: vector.hpp:58
int * GetI() const
Return the array I.
Definition: sparsemat.hpp:134
int * GetJ()
Definition: table.hpp:114
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
HYPRE_Int GetGlobalNumCols() const
Definition: hypre.hpp:382
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
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
int Dimension() const
Definition: mesh.hpp:645
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
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
void SetSize(int s)
Resize the vector to size s.
Definition: vector.hpp:328
HypreParMatrix & Add(const double beta, const HypreParMatrix &B)
Definition: hypre.hpp:433
HypreParVector * B
Right-hand side and solution vector.
Definition: hypre.hpp:627
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
void MakeDataOwner()
Definition: vector.hpp:114
void Mult(const Table &A, const Table &B, Table &C)
C = A * B (as boolean matrices)
Definition: table.cpp:478
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
int Width() const
Get the width (size of input) of the Operator. Synonym with NumCols().
Definition: operator.hpp:42
HypreParMatrix * ParallelAssemble() const
Returns the matrix "assembled" on the true dofs.
T * GetData()
Returns the data.
Definition: array.hpp:115
void GetEigenvalues(Array< double > &eigenvalues)
Collect the converged eigenvalues.
Definition: hypre.cpp:3638
int Size() const
Returns the size of the vector.
Definition: vector.hpp:120
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
virtual ~HypreGMRES()
Definition: hypre.cpp:2409
void SetPreconditioner(Solver &precond)
Definition: hypre.cpp:3294
void GetParBlocks(Array2D< HypreParMatrix *> &blocks) const
virtual void Mult(const Vector &x, Vector &y) const =0
Operator application: y=A(x).
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
bool iterative_mode
If true, use the second argument of Mult() as an initial guess.
Definition: operator.hpp:272
STL namespace.
virtual void Finalize(int skip_zeros=1)
Finalizes the matrix initialization.
virtual ~HypreAMS()
Definition: hypre.cpp:2892
void SetTol(double tol)
Definition: hypre.cpp:2308
void SetOperator(HypreParMatrix &A)
Definition: hypre.cpp:3610
virtual void Mult(const HypreParVector &b, HypreParVector &x) const
Relax the linear system Ax=b.
Definition: hypre.cpp:1981
void LoseData()
Lose the ownership of the graph (I, J) and data (A) arrays.
Definition: sparsemat.hpp:428
bool pos_l1_norms
If set, take absolute values of the computed l1_norms.
Definition: hypre.hpp:553
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
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
void SetPrintLevel(int print_lvl)
Definition: hypre.cpp:2185
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
void AddDomainInterpolator(DiscreteInterpolator *di)
void SetTol(double tol)
Definition: hypre.cpp:3257
HypreParVector * X1
Definition: hypre.hpp:527
void SetWindowParameters(double a, double b, double c)
Set parameters for windowing function for FIR smoother.
Definition: hypre.cpp:1857
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
int GetNV() const
Returns number of vertices. Vertices are only at the corners of elements, where you would expect them...
Definition: mesh.hpp:615
HypreLOBPCG(MPI_Comm comm)
Definition: hypre.cpp:3227
void Add(const DenseMatrix &A, const DenseMatrix &B, double alpha, DenseMatrix &C)
C = A + alpha*B.
Definition: densemat.cpp:2942
int NumRows() const
Definition: array.hpp:310
int dim
Definition: ex3.cpp:47
void SetSymmetry(int sym)
Definition: hypre.cpp:2438
virtual ~HypreParaSails()
Definition: hypre.cpp:2443
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 SetLogging(int logging)
Definition: hypre.cpp:2180
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
const FiniteElementCollection * FEColl() const
Definition: fespace.hpp:267
void SetRelTol(double rel_tol)
Definition: hypre.cpp:3579
HYPRE_Int GlobalTrueVSize() const
Definition: pfespace.hpp:247
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
HYPRE_Int * GetTrueDofOffsets() const
Definition: pfespace.hpp:244
int to_int(const std::string &str)
Definition: text.hpp:70
HypreParVector ** StealEigenvectors()
Transfer ownership of the converged eigenvectors.
Definition: hypre.cpp:3685
HYPRE_Int Randomize(HYPRE_Int seed)
Set random values.
Definition: hypre.cpp:213
double relax_weight
Damping coefficient (usually <= 1)
Definition: hypre.hpp:535
void SetElasticityOptions(ParFiniteElementSpace *fespace)
Definition: hypre.cpp:2666
void ParallelProject(Vector &tv) const
Returns the vector restricted to the true dofs.
Definition: pgridfunc.cpp:167
void Sort()
Sorts the array. This requires operator< to be defined for T.
Definition: array.hpp:237
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
int GetNE() const
Returns number of elements in the mesh.
Definition: fespace.hpp:277
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
double * GetData() const
Return a pointer to the beginning of the Vector data.
Definition: vector.hpp:129
int ParCSRRelax_FIR(hypre_ParCSRMatrix *A, hypre_ParVector *f, double max_eig, int poly_order, double *fir_coeffs, hypre_ParVector *u, hypre_ParVector *x0, hypre_ParVector *x1, hypre_ParVector *x2, hypre_ParVector *x3)
Definition: hypre.cpp:1699
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
int ParCSRRelax_Taubin(hypre_ParCSRMatrix *A, hypre_ParVector *f, double lambda, double mu, int N, double max_eig, hypre_ParVector *u, hypre_ParVector *r)
Definition: hypre.cpp:1662
Mesh * GetMesh() const
Returns the mesh.
Definition: fespace.hpp:224
void GetEigenvalues(Array< double > &eigenvalues)
Collect the converged eigenvalues.
Definition: hypre.cpp:3359
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
int GetOrder(int i) const
Returns the order of the i&#39;th finite element.
Definition: fespace.cpp:94
HypreParVector * X
Definition: hypre.hpp:523
void SetSize(int nsize)
Change logical size of the array, keep existing entries.
Definition: array.hpp:569
void Solve()
Solve the eigenproblem.
Definition: hypre.cpp:3632
void AddTraceFaceInterpolator(DiscreteInterpolator *di)
void SetNumModes(int num_eigs)
Definition: hypre.cpp:3565
int Height() const
Get the height (size of output) of the Operator. Synonym with NumRows().
Definition: operator.hpp:36
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
double Norml1() const
Returns the l_1 norm of the vector.
Definition: vector.cpp:713
int SpaceDimension() const
Definition: mesh.hpp:646
HypreParMatrix * Transpose() const
Returns the transpose of *this.
Definition: hypre.cpp:987
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
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 * 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 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
int * GetJ() const
Return the array J.
Definition: sparsemat.hpp:136
void InvScaleRows(const Vector &s)
Scale the local row i by 1./s(i)
Definition: hypre.cpp:1213
int Size() const
Returns the number of TYPE I elements.
Definition: table.hpp:92
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
int NumCols() const
Definition: array.hpp:311
void SetMassMatrix(HypreParMatrix &M)
Definition: hypre.cpp:3625
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
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 infinity()
Define a shortcut for std::numeric_limits<double>::infinity()
Definition: vector.hpp:42
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
void SetDataOwner(bool owna)
Set the data ownership flag (A array).
Definition: sparsemat.hpp:422
Abstract class for hypre&#39;s solvers and preconditioners.
Definition: hypre.hpp:620
int Size_of_connections() const
Definition: table.hpp:98
int Size() const
Logical size of the array.
Definition: array.hpp:133
HypreParVector & operator=(double d)
Set constant values.
Definition: hypre.cpp:186
MPI_Comm GetComm() const
Definition: pfespace.hpp:235
int GetFaceOrder(int i) const
Returns the order of the i&#39;th face finite element.
Definition: fespace.cpp:100
const double alpha
Definition: ex15.cpp:337
Arbitrary order H(curl)-conforming Nedelec finite elements.
Definition: fe_coll.hpp:230
Vector data type.
Definition: vector.hpp:48
HypreGMRES(HypreParMatrix &_A)
Definition: hypre.cpp:2290
hypre_ParCSRMatrix * StealData()
Changes the ownership of the the matrix.
Definition: hypre.cpp:831
void SetPrecondUsageMode(int pcg_mode)
Definition: hypre.cpp:3288
HypreParMatrix * A
The linear system matrix.
Definition: hypre.hpp:624
Arbitrary order H1-conforming (continuous) finite elements.
Definition: fe_coll.hpp:79
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
Base class for solvers.
Definition: operator.hpp:268
Class for parallel grid function.
Definition: pgridfunc.hpp:32
double * data
Definition: vector.hpp:53
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
void Swap(SparseMatrix &other)
Definition: sparsemat.cpp:3260
~HypreParVector()
Calls hypre&#39;s destroy function.
Definition: hypre.cpp:223
virtual void Assemble(int skip_zeros=1)
Class for parallel meshes.
Definition: pmesh.hpp:32
HypreParVector * Z
Definition: hypre.hpp:525
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
double Normlinf() const
Returns the l_infinity norm of the vector.
Definition: vector.cpp:703
void SetType(HypreSmoother::Type type, int relax_times=1)
Set the relaxation type and number of sweeps.
Definition: hypre.cpp:1816
double * GetData() const
Return element data, i.e. array A.
Definition: sparsemat.hpp:138
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
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
double sigma(const Vector &x)
Definition: maxwell.cpp:102