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