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