MFEM  v3.0
 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.googlecode.com.
7 //
8 // MFEM is free software; you can redistribute it and/or modify it under the
9 // terms of the GNU Lesser General Public License (as published by the Free
10 // Software Foundation) version 2.1 dated February 1999.
11 
12 #include "../config/config.hpp"
13 
14 #ifdef MFEM_USE_MPI
15 
16 #include "linalg.hpp"
17 #include "../fem/fem.hpp"
18 
19 #include <fstream>
20 #include <iomanip>
21 #include <cmath>
22 #include <cstdlib>
23 
24 using namespace std;
25 
26 namespace mfem
27 {
28 
29 HypreParVector::HypreParVector(MPI_Comm comm, int glob_size,
30  int *col) : Vector()
31 {
32  x = hypre_ParVectorCreate(comm,glob_size,col);
33  hypre_ParVectorInitialize(x);
34  hypre_ParVectorSetPartitioningOwner(x,0);
35  // The data will be destroyed by hypre (this is the default)
36  hypre_ParVectorSetDataOwner(x,1);
37  hypre_SeqVectorSetDataOwner(hypre_ParVectorLocalVector(x),1);
38  SetDataAndSize(hypre_VectorData(hypre_ParVectorLocalVector(x)),
39  hypre_VectorSize(hypre_ParVectorLocalVector(x)));
40  own_ParVector = 1;
41 }
42 
43 HypreParVector::HypreParVector(MPI_Comm comm, int glob_size,
44  double *_data, int *col) : Vector()
45 {
46  x = hypre_ParVectorCreate(comm,glob_size,col);
47  hypre_ParVectorSetDataOwner(x,1); // owns the seq vector
48  hypre_SeqVectorSetDataOwner(hypre_ParVectorLocalVector(x),0);
49  hypre_ParVectorSetPartitioningOwner(x,0);
50  hypre_VectorData(hypre_ParVectorLocalVector(x)) = _data;
51  // If hypre_ParVectorLocalVector(x) is non-NULL, hypre_ParVectorInitialize(x)
52  // does not allocate memory!
53  hypre_ParVectorInitialize(x);
54  SetDataAndSize(hypre_VectorData(hypre_ParVectorLocalVector(x)),
55  hypre_VectorSize(hypre_ParVectorLocalVector(x)));
56  own_ParVector = 1;
57 }
58 
60 {
61  x = hypre_ParVectorCreate(y.x -> comm, y.x -> global_size,
62  y.x -> partitioning);
63  hypre_ParVectorInitialize(x);
64  hypre_ParVectorSetPartitioningOwner(x,0);
65  hypre_ParVectorSetDataOwner(x,1);
66  hypre_SeqVectorSetDataOwner(hypre_ParVectorLocalVector(x),1);
67  SetDataAndSize(hypre_VectorData(hypre_ParVectorLocalVector(x)),
68  hypre_VectorSize(hypre_ParVectorLocalVector(x)));
69  own_ParVector = 1;
70 }
71 
73 {
74  if (!tr)
75  x = hypre_ParVectorInDomainOf(A);
76  else
77  x = hypre_ParVectorInRangeOf(A);
78  SetDataAndSize(hypre_VectorData(hypre_ParVectorLocalVector(x)),
79  hypre_VectorSize(hypre_ParVectorLocalVector(x)));
80  own_ParVector = 1;
81 }
82 
83 HypreParVector::HypreParVector(HYPRE_ParVector y) : Vector()
84 {
85  x = (hypre_ParVector *) y;
86  SetDataAndSize(hypre_VectorData(hypre_ParVectorLocalVector(x)),
87  hypre_VectorSize(hypre_ParVectorLocalVector(x)));
88  own_ParVector = 0;
89 }
90 
92 {
93  x = hypre_ParVectorCreate(pfes->GetComm(), pfes->GlobalTrueVSize(),
94  pfes->GetTrueDofOffsets());
95  hypre_ParVectorInitialize(x);
96  hypre_ParVectorSetPartitioningOwner(x,0);
97  // The data will be destroyed by hypre (this is the default)
98  hypre_ParVectorSetDataOwner(x,1);
99  hypre_SeqVectorSetDataOwner(hypre_ParVectorLocalVector(x),1);
100  SetDataAndSize(hypre_VectorData(hypre_ParVectorLocalVector(x)),
101  hypre_VectorSize(hypre_ParVectorLocalVector(x)));
102  own_ParVector = 1;
103 }
104 
105 HypreParVector::operator hypre_ParVector*() const
106 {
107  return x;
108 }
109 
110 #ifndef HYPRE_PAR_VECTOR_STRUCT
111 HypreParVector::operator HYPRE_ParVector() const
112 {
113  return (HYPRE_ParVector) x;
114 }
115 #endif
116 
118 {
119  hypre_Vector *hv = hypre_ParVectorToVectorAll(*this);
120  Vector *v = new Vector(hv->data, hv->size);
121  v->MakeDataOwner();
122  hypre_SeqVectorSetDataOwner(hv,0);
123  hypre_SeqVectorDestroy(hv);
124  return v;
125 }
126 
128 {
129  hypre_ParVectorSetConstantValues(x,d);
130  return *this;
131 }
132 
134 {
135 #ifdef MFEM_DEBUG
136  if (size != y.Size())
137  mfem_error("HypreParVector::operator=");
138 #endif
139 
140  for (int i = 0; i < size; i++)
141  data[i] = y.data[i];
142  return *this;
143 }
144 
145 void HypreParVector::SetData(double *_data)
146 {
147  Vector::data = hypre_VectorData(hypre_ParVectorLocalVector(x)) = _data;
148 }
149 
151 {
152  return hypre_ParVectorSetRandomValues(x,seed);
153 }
154 
155 void HypreParVector::Print(const char *fname)
156 {
157  hypre_ParVectorPrint(x,fname);
158 }
159 
161 {
162  if (own_ParVector)
163  hypre_ParVectorDestroy(x);
164 }
165 
166 
168 {
169  return hypre_ParVectorInnerProd(*x, *y);
170 }
171 
173 {
174  return hypre_ParVectorInnerProd(x, y);
175 }
176 
177 
178 HypreParMatrix::HypreParMatrix(MPI_Comm comm, int glob_size, int *row_starts,
179  SparseMatrix *diag)
180  : Operator(diag->Height(), diag->Width())
181 {
182  A = hypre_ParCSRMatrixCreate(comm, glob_size, glob_size, row_starts,
183  row_starts, 0, diag->NumNonZeroElems(), 0);
184  hypre_ParCSRMatrixSetDataOwner(A,0);
185  hypre_ParCSRMatrixSetRowStartsOwner(A,0);
186  hypre_ParCSRMatrixSetColStartsOwner(A,0);
187 
188  hypre_CSRMatrixSetDataOwner(A->diag,0);
189  hypre_CSRMatrixI(A->diag) = diag->GetI();
190  hypre_CSRMatrixJ(A->diag) = diag->GetJ();
191  hypre_CSRMatrixData(A->diag) = diag->GetData();
192  hypre_CSRMatrixSetRownnz(A->diag);
193 
194  hypre_CSRMatrixSetDataOwner(A->offd,1);
195  hypre_CSRMatrixI(A->offd) = hypre_CTAlloc(int, diag->Height()+1);
196 
197  /* Don't need to call these, since they allocate memory only
198  if it was not already allocated */
199  // hypre_CSRMatrixInitialize(A->diag);
200  // hypre_ParCSRMatrixInitialize(A);
201 
202  hypre_ParCSRMatrixSetNumNonzeros(A);
203 
204  /* Make sure that the first entry in each row is the diagonal one. */
205  hypre_CSRMatrixReorder(hypre_ParCSRMatrixDiag(A));
206 
207  hypre_MatvecCommPkgCreate(A);
208 
209  CommPkg = NULL;
210  X = Y = NULL;
211 }
212 
213 
215  int global_num_rows, int global_num_cols,
216  int *row_starts, int *col_starts,
217  SparseMatrix *diag)
218  : Operator(diag->Height(), diag->Width())
219 {
220  A = hypre_ParCSRMatrixCreate(comm, global_num_rows, global_num_cols,
221  row_starts, col_starts,
222  0, diag->NumNonZeroElems(), 0);
223  hypre_ParCSRMatrixSetDataOwner(A,0);
224  hypre_ParCSRMatrixSetRowStartsOwner(A,0);
225  hypre_ParCSRMatrixSetColStartsOwner(A,0);
226 
227  hypre_CSRMatrixSetDataOwner(A->diag,0);
228  hypre_CSRMatrixI(A->diag) = diag->GetI();
229  hypre_CSRMatrixJ(A->diag) = diag->GetJ();
230  hypre_CSRMatrixData(A->diag) = diag->GetData();
231  hypre_CSRMatrixSetRownnz(A->diag);
232 
233  hypre_CSRMatrixSetDataOwner(A->offd,1);
234  hypre_CSRMatrixI(A->offd) = hypre_CTAlloc(int, diag->Height()+1);
235 
236  hypre_ParCSRMatrixSetNumNonzeros(A);
237 
238  /* Make sure that the first entry in each row is the diagonal one. */
239  if (row_starts == col_starts)
240  hypre_CSRMatrixReorder(hypre_ParCSRMatrixDiag(A));
241 
242  hypre_MatvecCommPkgCreate(A);
243 
244  CommPkg = NULL;
245  X = Y = NULL;
246 }
247 
249  int global_num_rows, int global_num_cols,
250  int *row_starts, int *col_starts,
251  SparseMatrix *diag, SparseMatrix *offd,
252  int *cmap)
253  : Operator(diag->Height(), diag->Width())
254 {
255  A = hypre_ParCSRMatrixCreate(comm, global_num_rows, global_num_cols,
256  row_starts, col_starts,
257  offd->Width(), diag->NumNonZeroElems(),
258  offd->NumNonZeroElems());
259  hypre_ParCSRMatrixSetDataOwner(A,0);
260  hypre_ParCSRMatrixSetRowStartsOwner(A,0);
261  hypre_ParCSRMatrixSetColStartsOwner(A,0);
262 
263  hypre_CSRMatrixSetDataOwner(A->diag,0);
264  hypre_CSRMatrixI(A->diag) = diag->GetI();
265  hypre_CSRMatrixJ(A->diag) = diag->GetJ();
266  hypre_CSRMatrixData(A->diag) = diag->GetData();
267  hypre_CSRMatrixSetRownnz(A->diag);
268 
269  hypre_CSRMatrixSetDataOwner(A->offd,0);
270  hypre_CSRMatrixI(A->offd) = offd->GetI();
271  hypre_CSRMatrixJ(A->offd) = offd->GetJ();
272  hypre_CSRMatrixData(A->offd) = offd->GetData();
273  hypre_CSRMatrixSetRownnz(A->offd);
274 
275  hypre_ParCSRMatrixColMapOffd(A) = cmap;
276 
277  hypre_ParCSRMatrixSetNumNonzeros(A);
278 
279  /* Make sure that the first entry in each row is the diagonal one. */
280  if (row_starts == col_starts)
281  hypre_CSRMatrixReorder(hypre_ParCSRMatrixDiag(A));
282 
283  hypre_MatvecCommPkgCreate(A);
284 
285  CommPkg = NULL;
286  X = Y = NULL;
287 }
288 
289 HypreParMatrix::HypreParMatrix(MPI_Comm comm, int *row_starts, int *col_starts,
290  SparseMatrix *sm_a)
291 {
292 #ifdef MFEM_DEBUG
293  if (sm_a == NULL)
294  mfem_error("HypreParMatrix::HypreParMatrix: sm_a==NULL");
295 #endif
296 
297  hypre_CSRMatrix *csr_a;
298 
299  csr_a = hypre_CSRMatrixCreate(sm_a -> Height(), sm_a -> Width(),
300  sm_a -> NumNonZeroElems());
301 
302  hypre_CSRMatrixSetDataOwner(csr_a,0);
303  hypre_CSRMatrixI(csr_a) = sm_a -> GetI();
304  hypre_CSRMatrixJ(csr_a) = sm_a -> GetJ();
305  hypre_CSRMatrixData(csr_a) = sm_a -> GetData();
306  hypre_CSRMatrixSetRownnz(csr_a);
307 
308  A = hypre_CSRMatrixToParCSRMatrix(comm, csr_a, row_starts, col_starts);
309 
310  CommPkg = NULL;
311  X = Y = NULL;
312 
313  height = GetNumRows();
314  width = GetNumCols();
315 
316  /* Make sure that the first entry in each row is the diagonal one. */
317  if (row_starts == col_starts)
318  hypre_CSRMatrixReorder(hypre_ParCSRMatrixDiag(A));
319 
320  hypre_MatvecCommPkgCreate(A);
321 }
322 
324  int global_num_rows, int global_num_cols,
325  int *row_starts, int *col_starts, Table *diag)
326 {
327  int nnz = diag->Size_of_connections();
328  A = hypre_ParCSRMatrixCreate(comm, global_num_rows, global_num_cols,
329  row_starts, col_starts, 0, nnz, 0);
330  hypre_ParCSRMatrixSetDataOwner(A,1);
331  hypre_ParCSRMatrixSetRowStartsOwner(A,0);
332  hypre_ParCSRMatrixSetColStartsOwner(A,0);
333 
334  hypre_CSRMatrixSetDataOwner(A->diag,1);
335  hypre_CSRMatrixI(A->diag) = diag->GetI();
336  hypre_CSRMatrixJ(A->diag) = diag->GetJ();
337 
338  hypre_CSRMatrixData(A->diag) = hypre_TAlloc(double, nnz);
339  for (int k = 0; k < nnz; k++)
340  (hypre_CSRMatrixData(A->diag))[k] = 1.0;
341 
342  hypre_CSRMatrixSetRownnz(A->diag);
343 
344  hypre_CSRMatrixSetDataOwner(A->offd,1);
345  hypre_CSRMatrixI(A->offd) = hypre_CTAlloc(int, diag->Size()+1);
346 
347  hypre_ParCSRMatrixSetNumNonzeros(A);
348 
349  /* Make sure that the first entry in each row is the diagonal one. */
350  if (row_starts == col_starts)
351  hypre_CSRMatrixReorder(hypre_ParCSRMatrixDiag(A));
352 
353  hypre_MatvecCommPkgCreate(A);
354 
355  CommPkg = NULL;
356  X = Y = NULL;
357 
358  height = GetNumRows();
359  width = GetNumCols();
360 }
361 
362 HypreParMatrix::HypreParMatrix(MPI_Comm comm, int id, int np,
363  int *row, int *col,
364  int *i_diag, int *j_diag,
365  int *i_offd, int *j_offd,
366  int *cmap, int cmap_size)
367 {
368  int diag_col, offd_col;
369 
370  if (HYPRE_AssumedPartitionCheck())
371  {
372  diag_col = i_diag[row[1]-row[0]];
373  offd_col = i_offd[row[1]-row[0]];
374 
375  A = hypre_ParCSRMatrixCreate(comm, row[2], col[2], row, col,
376  cmap_size, diag_col, offd_col);
377  }
378  else
379  {
380  diag_col = i_diag[row[id+1]-row[id]];
381  offd_col = i_offd[row[id+1]-row[id]];
382 
383  A = hypre_ParCSRMatrixCreate(comm, row[np], col[np], row, col,
384  cmap_size, diag_col, offd_col);
385  }
386 
387  hypre_ParCSRMatrixSetDataOwner(A,1);
388  hypre_ParCSRMatrixSetRowStartsOwner(A,0);
389  hypre_ParCSRMatrixSetColStartsOwner(A,0);
390 
391  int i;
392 
393  double *a_diag = hypre_TAlloc(double, diag_col);
394  for (i = 0; i < diag_col; i++)
395  a_diag[i] = 1.0;
396 
397  double *a_offd = hypre_TAlloc(double, offd_col);
398  for (i = 0; i < offd_col; i++)
399  a_offd[i] = 1.0;
400 
401  hypre_CSRMatrixSetDataOwner(A->diag,1);
402  hypre_CSRMatrixI(A->diag) = i_diag;
403  hypre_CSRMatrixJ(A->diag) = j_diag;
404  hypre_CSRMatrixData(A->diag) = a_diag;
405  hypre_CSRMatrixSetRownnz(A->diag);
406 
407  hypre_CSRMatrixSetDataOwner(A->offd,1);
408  hypre_CSRMatrixI(A->offd) = i_offd;
409  hypre_CSRMatrixJ(A->offd) = j_offd;
410  hypre_CSRMatrixData(A->offd) = a_offd;
411  hypre_CSRMatrixSetRownnz(A->offd);
412 
413  hypre_ParCSRMatrixColMapOffd(A) = cmap;
414 
415  hypre_ParCSRMatrixSetNumNonzeros(A);
416 
417  /* Make sure that the first entry in each row is the diagonal one. */
418  if (row == col)
419  hypre_CSRMatrixReorder(hypre_ParCSRMatrixDiag(A));
420 
421  hypre_MatvecCommPkgCreate(A);
422 
423  CommPkg = NULL;
424  X = Y = NULL;
425 
426  height = GetNumRows();
427  width = GetNumCols();
428 }
429 
430 HypreParMatrix::HypreParMatrix(MPI_Comm comm, int nrows, int glob_nrows,
431  int glob_ncols, int *I, int *J, double *data,
432  int *rows, int *cols)
433 {
434  // construct the local CSR matrix
435  int nnz = I[nrows];
436  hypre_CSRMatrix *local = hypre_CSRMatrixCreate(nrows, glob_ncols, nnz);
437  hypre_CSRMatrixI(local) = I;
438  hypre_CSRMatrixJ(local) = J;
439  hypre_CSRMatrixData(local) = data;
440  hypre_CSRMatrixRownnz(local) = NULL;
441  hypre_CSRMatrixOwnsData(local) = 1;
442  hypre_CSRMatrixNumRownnz(local) = nrows;
443 
444  int part_size, myid;
445  if (HYPRE_AssumedPartitionCheck())
446  {
447  myid = 0;
448  part_size = 2;
449  }
450  else
451  {
452  MPI_Comm_rank(comm, &myid);
453  MPI_Comm_size(comm, &part_size);
454  part_size++;
455  }
456 
457  // copy in the row and column partitionings
458  int *row_starts = hypre_TAlloc(int, part_size);
459  int *col_starts = hypre_TAlloc(int, part_size);
460  for (int i = 0; i < part_size; i++)
461  {
462  row_starts[i] = rows[i];
463  col_starts[i] = cols[i];
464  }
465 
466  // construct the global ParCSR matrix
467  A = hypre_ParCSRMatrixCreate(comm, glob_nrows, glob_ncols,
468  row_starts, col_starts, 0, 0, 0);
469  hypre_ParCSRMatrixOwnsRowStarts(A) = 1;
470  hypre_ParCSRMatrixOwnsColStarts(A) = 1;
471  GenerateDiagAndOffd(local, A, col_starts[myid], col_starts[myid+1]-1);
472  hypre_ParCSRMatrixSetNumNonzeros(A);
473  /* Make sure that the first entry in each row is the diagonal one. */
474  if (rows == cols)
475  hypre_CSRMatrixReorder(hypre_ParCSRMatrixDiag(A));
476  hypre_MatvecCommPkgCreate(A);
477 
478  // delete the local CSR matrix
479  hypre_CSRMatrixI(local) = NULL;
480  hypre_CSRMatrixJ(local) = NULL;
481  hypre_CSRMatrixData(local) = NULL;
482  hypre_CSRMatrixDestroy(local);
483 
484  CommPkg = NULL;
485  X = Y = NULL;
486  height = GetNumRows();
487  width = GetNumCols();
488 }
489 
490 void HypreParMatrix::SetCommPkg(hypre_ParCSRCommPkg *comm_pkg)
491 {
492  CommPkg = comm_pkg;
493 
494  if (hypre_ParCSRMatrixCommPkg(A))
495  hypre_MatvecCommPkgDestroy(hypre_ParCSRMatrixCommPkg(A));
496 
497  hypre_ParCSRMatrixCommPkg(A) = comm_pkg;
498 }
499 
501 {
502 #ifdef MFEM_DEBUG
503  if (CommPkg == NULL || CommPkg != hypre_ParCSRMatrixCommPkg(A))
504  mfem_error("\nHypreParMatrix::CheckCommPkg()");
505 #endif
506 }
507 
509 {
510  if (CommPkg == NULL)
511  return;
512  hypre_TFree(CommPkg->send_procs);
513  hypre_TFree(CommPkg->send_map_starts);
514  hypre_TFree(CommPkg->send_map_elmts);
515  hypre_TFree(CommPkg->recv_procs);
516  hypre_TFree(CommPkg->recv_vec_starts);
517  if (CommPkg->send_mpi_types)
518  hypre_TFree(CommPkg->send_mpi_types);
519  if (CommPkg->recv_mpi_types)
520  hypre_TFree(CommPkg->recv_mpi_types);
521  if (hypre_ParCSRMatrixCommPkg(A) == CommPkg)
522  hypre_ParCSRMatrixCommPkg(A) = NULL;
523  delete CommPkg;
524  CommPkg = NULL;
525 }
526 
527 HypreParMatrix::operator hypre_ParCSRMatrix*()
528 {
529  return (this) ? A : NULL;
530 }
531 
532 #ifndef HYPRE_PAR_CSR_MATRIX_STRUCT
533 HypreParMatrix::operator HYPRE_ParCSRMatrix()
534 {
535  return (this) ? (HYPRE_ParCSRMatrix) A : (HYPRE_ParCSRMatrix) NULL;
536 }
537 #endif
538 
539 hypre_ParCSRMatrix* HypreParMatrix::StealData()
540 {
541  hypre_ParCSRMatrix *R = A;
542  A = NULL;
543  return R;
544 }
545 
547 {
548  int size=hypre_CSRMatrixNumRows(A->diag);
549  diag.SetSize(size);
550  for (int j = 0; j < size; j++)
551  {
552  diag(j) = A->diag->data[A->diag->i[j]];
553 #ifdef MFEM_DEBUG
554  if (A->diag->j[A->diag->i[j]] != j)
555  mfem_error("HypreParMatrix::GetDiag");
556 #endif
557  }
558 }
559 
560 
562 {
563  hypre_ParCSRMatrix * At;
564  hypre_ParCSRMatrixTranspose(A, &At, 1);
565  hypre_ParCSRMatrixSetNumNonzeros(At);
566 
567  hypre_MatvecCommPkgCreate(At);
568 
569  return new HypreParMatrix(At);
570 }
571 
573  double a, double b)
574 {
575  return hypre_ParCSRMatrixMatvec(a, A, x, b, y);
576 }
577 
578 void HypreParMatrix::Mult(double a, const Vector &x, double b, Vector &y) const
579 {
580  if (X == NULL)
581  {
582  X = new HypreParVector(A->comm,
584  x.GetData(),
585  GetColStarts());
586  Y = new HypreParVector(A->comm,
588  y.GetData(),
589  GetRowStarts());
590  }
591  else
592  {
593  X->SetData(x.GetData());
594  Y->SetData(y.GetData());
595  }
596 
597  hypre_ParCSRMatrixMatvec(a, A, *X, b, *Y);
598 }
599 
600 void HypreParMatrix::MultTranspose(double a, const Vector &x,
601  double b, Vector &y) const
602 {
603  // Note: x has the dimensions of Y (height), and
604  // y has the dimensions of X (width)
605  if (X == NULL)
606  {
607  X = new HypreParVector(A->comm,
609  y.GetData(),
610  GetColStarts());
611  Y = new HypreParVector(A->comm,
613  x.GetData(),
614  GetRowStarts());
615  }
616  else
617  {
618  X->SetData(y.GetData());
619  Y->SetData(x.GetData());
620  }
621 
622  hypre_ParCSRMatrixMatvecT(a, A, *Y, b, *X);
623 }
624 
625 int HypreParMatrix::Mult(HYPRE_ParVector x, HYPRE_ParVector y,
626  double a, double b)
627 {
628  return hypre_ParCSRMatrixMatvec(a,A,(hypre_ParVector *)x,b,(hypre_ParVector *)y);
629 }
630 
632  double a, double b)
633 {
634  return hypre_ParCSRMatrixMatvecT(a,A,x,b,y);
635 }
636 
638 {
639 
640  if(hypre_CSRMatrixNumRows(A->diag) != hypre_CSRMatrixNumRows(A->offd))
641  mfem_error("Row does not match");
642 
643  if(hypre_CSRMatrixNumRows(A->diag) != diag.Size())
644  mfem_error("Note the Vector diag is not of compatible dimensions with A\n");
645 
646  int size=hypre_CSRMatrixNumRows(A->diag);
647  double *Adiag_data = hypre_CSRMatrixData(A->diag);
648  HYPRE_Int *Adiag_i = hypre_CSRMatrixI(A->diag);
649 
650 
651  double *Aoffd_data = hypre_CSRMatrixData(A->offd);
652  HYPRE_Int *Aoffd_i = hypre_CSRMatrixI(A->offd);
653  double val;
654  int jj;
655  for (int i(0); i < size; ++i)
656  {
657  val = diag[i];
658  for (jj = Adiag_i[i]; jj < Adiag_i[i+1]; ++jj)
659  Adiag_data[jj] *= val;
660  for (jj = Aoffd_i[i]; jj < Aoffd_i[i+1]; ++jj)
661  Aoffd_data[jj] *= val;
662  }
663 }
664 
666 {
667 
668  if(hypre_CSRMatrixNumRows(A->diag) != hypre_CSRMatrixNumRows(A->offd))
669  mfem_error("Row does not match");
670 
671  if(hypre_CSRMatrixNumRows(A->diag) != diag.Size())
672  mfem_error("Note the Vector diag is not of compatible dimensions with A\n");
673 
674  int size=hypre_CSRMatrixNumRows(A->diag);
675  double *Adiag_data = hypre_CSRMatrixData(A->diag);
676  HYPRE_Int *Adiag_i = hypre_CSRMatrixI(A->diag);
677 
678 
679  double *Aoffd_data = hypre_CSRMatrixData(A->offd);
680  HYPRE_Int *Aoffd_i = hypre_CSRMatrixI(A->offd);
681  double val;
682  int jj;
683  for (int i(0); i < size; ++i)
684  {
685 #ifdef MFEM_DEBUG
686  if(0.0 == diag(i))
687  mfem_error("HypreParMatrix::InvDiagScale : Division by 0");
688 #endif
689  val = 1./diag(i);
690  for (jj = Adiag_i[i]; jj < Adiag_i[i+1]; ++jj)
691  Adiag_data[jj] *= val;
692  for (jj = Aoffd_i[i]; jj < Aoffd_i[i+1]; ++jj)
693  Aoffd_data[jj] *= val;
694  }
695 }
696 
698 {
699  if (hypre_CSRMatrixNumRows(A->diag) != hypre_CSRMatrixNumRows(A->offd))
700  mfem_error("Row does not match");
701 
702  int size=hypre_CSRMatrixNumRows(A->diag);
703  int jj;
704 
705  double *Adiag_data = hypre_CSRMatrixData(A->diag);
706  HYPRE_Int *Adiag_i = hypre_CSRMatrixI(A->diag);
707  for (jj = 0; jj < Adiag_i[size]; ++jj)
708  Adiag_data[jj] *= s;
709 
710  double *Aoffd_data = hypre_CSRMatrixData(A->offd);
711  HYPRE_Int *Aoffd_i = hypre_CSRMatrixI(A->offd);
712  for (jj = 0; jj < Aoffd_i[size]; ++jj)
713  Aoffd_data[jj] *= s;
714 }
715 
716 void HypreParMatrix::Print(const char *fname, int offi, int offj)
717 {
718  hypre_ParCSRMatrixPrintIJ(A,offi,offj,fname);
719 }
720 
721 void HypreParMatrix::Read(MPI_Comm comm, const char *fname)
722 {
723  if (A) hypre_ParCSRMatrixDestroy(A);
724  int io,jo;
725  hypre_ParCSRMatrixReadIJ(comm, fname, &io, &jo, &A);
726  hypre_ParCSRMatrixSetNumNonzeros(A);
727 
728  hypre_MatvecCommPkgCreate(A);
729 }
730 
732 {
733  DestroyCommPkg();
734 
735  if (A)
736  {
737  if (hypre_ParCSRMatrixCommPkg(A))
738  hypre_MatvecCommPkgDestroy(hypre_ParCSRMatrixCommPkg(A));
739  hypre_ParCSRMatrixCommPkg(A) = NULL;
740 
741  if (hypre_CSRMatrixOwnsData(A->diag))
742  {
743  hypre_CSRMatrixDestroy(A->diag);
744  A->diag = NULL;
745  }
746  else
747  {
748  if (hypre_CSRMatrixRownnz(A->diag))
749  hypre_TFree(hypre_CSRMatrixRownnz(A->diag));
750  hypre_TFree(A->diag);
751  A->diag = NULL;
752  }
753 
754  if (hypre_CSRMatrixOwnsData(A->offd))
755  {
756  hypre_CSRMatrixDestroy(A->offd);
757  A->offd = NULL;
758  }
759  else
760  {
761  if (hypre_CSRMatrixRownnz(A->offd))
762  hypre_TFree(hypre_CSRMatrixRownnz(A->offd));
763  }
764  hypre_ParCSRMatrixDestroy(A);
765  }
766 
767  delete X;
768  delete Y;
769 }
770 
772 {
773  hypre_ParCSRMatrix * ab;
774  ab = hypre_ParMatmul(*A,*B);
775 
776  hypre_MatvecCommPkgCreate(ab);
777 
778  return new HypreParMatrix(ab);
779 }
780 
782 {
783  int P_owns_its_col_starts =
784  hypre_ParCSRMatrixOwnsColStarts((hypre_ParCSRMatrix*)(*P));
785  hypre_ParCSRMatrix * rap;
786  hypre_BoomerAMGBuildCoarseOperator(*P,*A,*P,&rap);
787 
788  hypre_ParCSRMatrixSetNumNonzeros(rap);
789  // hypre_MatvecCommPkgCreate(rap);
790  if (!P_owns_its_col_starts)
791  {
792  /* Warning: hypre_BoomerAMGBuildCoarseOperator steals the col_starts
793  from P (even if it does not own them)! */
794  hypre_ParCSRMatrixSetRowStartsOwner(rap,0);
795  hypre_ParCSRMatrixSetColStartsOwner(rap,0);
796  }
797  return new HypreParMatrix(rap);
798 }
799 
801 {
802  int P_owns_its_col_starts =
803  hypre_ParCSRMatrixOwnsColStarts((hypre_ParCSRMatrix*)(*P));
804  int Rt_owns_its_col_starts =
805  hypre_ParCSRMatrixOwnsColStarts((hypre_ParCSRMatrix*)(*Rt));
806 
807  hypre_ParCSRMatrix * rap;
808  hypre_BoomerAMGBuildCoarseOperator(*Rt,*A,*P,&rap);
809 
810  hypre_ParCSRMatrixSetNumNonzeros(rap);
811  // hypre_MatvecCommPkgCreate(rap);
812  if (!P_owns_its_col_starts)
813  {
814  /* Warning: hypre_BoomerAMGBuildCoarseOperator steals the col_starts
815  from P (even if it does not own them)! */
816  hypre_ParCSRMatrixSetColStartsOwner(rap,0);
817  }
818  if (!Rt_owns_its_col_starts)
819  {
820  /* Warning: hypre_BoomerAMGBuildCoarseOperator steals the col_starts
821  from P (even if it does not own them)! */
822  hypre_ParCSRMatrixSetRowStartsOwner(rap,0);
823  }
824  return new HypreParMatrix(rap);
825 }
826 
828  Array<int> &ess_dof_list,
830 {
831  // b -= Ae*x
832  Ae.Mult(x, b, -1.0, 1.0);
833 
834  hypre_CSRMatrix *A_diag = hypre_ParCSRMatrixDiag((hypre_ParCSRMatrix *)A);
835  double *data = hypre_CSRMatrixData(A_diag);
836  int *I = hypre_CSRMatrixI(A_diag);
837 #ifdef MFEM_DEBUG
838  int *J = hypre_CSRMatrixJ(A_diag);
839  int *I_offd =
840  hypre_CSRMatrixI(hypre_ParCSRMatrixOffd((hypre_ParCSRMatrix *)A));
841 #endif
842 
843  for (int i = 0; i < ess_dof_list.Size(); i++)
844  {
845  int r = ess_dof_list[i];
846  b(r) = data[I[r]] * x(r);
847 #ifdef MFEM_DEBUG
848  // Check that in the rows specified by the ess_dof_list, the matrix A has
849  // only one entry -- the diagonal.
850  if (I[r+1] != I[r]+1 || J[I[r]] != r || I_offd[r] != I_offd[r+1])
851  mfem_error("EliminateBC (hypre.cpp)");
852 #endif
853  }
854 }
855 
856 // Taubin or "lambda-mu" scheme, which alternates between positive and
857 // negative step sizes to approximate low-pass filter effect.
858 
859 int ParCSRRelax_Taubin(hypre_ParCSRMatrix *A, // matrix to relax with
860  hypre_ParVector *f, // right-hand side
861  double lambda,
862  double mu,
863  int N,
864  double max_eig,
865  hypre_ParVector *u, // initial/updated approximation
866  hypre_ParVector *r // another temp vector
867  )
868 {
869  hypre_CSRMatrix *A_diag = hypre_ParCSRMatrixDiag(A);
870  int num_rows = hypre_CSRMatrixNumRows(A_diag);
871 
872  double *u_data = hypre_VectorData(hypre_ParVectorLocalVector(u));
873  double *r_data = hypre_VectorData(hypre_ParVectorLocalVector(r));
874 
875  for (int i = 0; i < N; i++)
876  {
877  // get residual: r = f - A*u
878  hypre_ParVectorCopy(f, r);
879  hypre_ParCSRMatrixMatvec(-1.0, A, u, 1.0, r);
880 
881  double coef;
882  (0 == (i % 2)) ? coef = lambda : coef = mu;
883 
884  for (int i = 0; i < num_rows; i++)
885  {
886  u_data[i] += coef*r_data[i] / max_eig;
887  }
888  }
889 
890  return 0;
891 }
892 
893 // FIR scheme, which uses Chebyshev polynomials and a window function
894 // to approximate a low-pass step filter.
895 
896 int ParCSRRelax_FIR(hypre_ParCSRMatrix *A, // matrix to relax with
897  hypre_ParVector *f, // right-hand side
898  double max_eig,
899  int poly_order,
900  double* fir_coeffs,
901  hypre_ParVector *u, // initial/updated approximation
902  hypre_ParVector *x0, // temporaries
903  hypre_ParVector *x1,
904  hypre_ParVector *x2,
905  hypre_ParVector *x3)
906 
907 {
908  hypre_CSRMatrix *A_diag = hypre_ParCSRMatrixDiag(A);
909  int num_rows = hypre_CSRMatrixNumRows(A_diag);
910 
911  double *u_data = hypre_VectorData(hypre_ParVectorLocalVector(u));
912 
913  double *x0_data = hypre_VectorData(hypre_ParVectorLocalVector(x0));
914  double *x1_data = hypre_VectorData(hypre_ParVectorLocalVector(x1));
915  double *x2_data = hypre_VectorData(hypre_ParVectorLocalVector(x2));
916  double *x3_data = hypre_VectorData(hypre_ParVectorLocalVector(x3));
917 
918  hypre_ParVectorCopy(u, x0);
919 
920  // x1 = f -A*x0/max_eig
921  hypre_ParVectorCopy(f, x1);
922  hypre_ParCSRMatrixMatvec(-1.0, A, x0, 1.0, x1);
923 
924  for (int i = 0; i < num_rows; i++)
925  {
926  x1_data[i] /= -max_eig;
927  }
928 
929  // x1 = x0 -x1
930  for (int i = 0; i < num_rows; i++)
931  {
932  x1_data[i] = x0_data[i] -x1_data[i];
933  }
934 
935  // x3 = f0*x0 +f1*x1
936  for (int i = 0; i < num_rows; i++)
937  {
938  x3_data[i] = fir_coeffs[0]*x0_data[i] +fir_coeffs[1]*x1_data[i];
939  }
940 
941  for (int n = 2; n <= poly_order; n++)
942  {
943  // x2 = f - A*x1/max_eig
944  hypre_ParVectorCopy(f, x2);
945  hypre_ParCSRMatrixMatvec(-1.0, A, x1, 1.0, x2);
946 
947  for (int i = 0; i < num_rows; i++)
948  {
949  x2_data[i] /= -max_eig;
950  }
951 
952  // x2 = (x1-x0) +(x1-2*x2)
953  // x3 = x3 +f[n]*x2
954  // x0 = x1
955  // x1 = x2
956 
957  for (int i = 0; i < num_rows; i++)
958  {
959  x2_data[i] = (x1_data[i]-x0_data[i]) +(x1_data[i]-2*x2_data[i]);
960  x3_data[i] += fir_coeffs[n]*x2_data[i];
961  x0_data[i] = x1_data[i];
962  x1_data[i] = x2_data[i];
963  }
964  }
965 
966  for (int i = 0; i < num_rows; i++)
967  {
968  u_data[i] = x3_data[i];
969  }
970 
971  return 0;
972 }
973 
975 {
976  type = 2;
977  relax_times = 1;
978  relax_weight = 1.0;
979  omega = 1.0;
980  poly_order = 2;
981  poly_fraction = .3;
982  lambda = 0.5;
983  mu = -0.5;
984  taubin_iter = 40;
985 
986  l1_norms = NULL;
987  B = X = V = Z = NULL;
988  X0 = X1 = NULL;
989  fir_coeffs = NULL;
990 }
991 
993  int _relax_times, double _relax_weight, double _omega,
994  int _poly_order, double _poly_fraction)
995 {
996  type = _type;
997  relax_times = _relax_times;
998  relax_weight = _relax_weight;
999  omega = _omega;
1000  poly_order = _poly_order;
1001  poly_fraction = _poly_fraction;
1002 
1003  l1_norms = NULL;
1004  B = X = V = Z = NULL;
1005  X0 = X1 = NULL;
1006  fir_coeffs = NULL;
1007 
1008  SetOperator(_A);
1009 }
1010 
1011 void HypreSmoother::SetType(HypreSmoother::Type _type, int _relax_times)
1012 {
1013  type = static_cast<int>(_type);
1014  relax_times = _relax_times;
1015 }
1016 
1017 void HypreSmoother::SetSOROptions(double _relax_weight, double _omega)
1018 {
1019  relax_weight = _relax_weight;
1020  omega = _omega;
1021 }
1022 
1023 void HypreSmoother::SetPolyOptions(int _poly_order, double _poly_fraction)
1024 {
1025  poly_order = _poly_order;
1026  poly_fraction = _poly_fraction;
1027 }
1028 
1029 void HypreSmoother::SetTaubinOptions(double _lambda, double _mu,
1030  int _taubin_iter)
1031 {
1032  lambda = _lambda;
1033  mu = _mu;
1034  taubin_iter = _taubin_iter;
1035 }
1036 
1037 void HypreSmoother::SetWindowByName(const char* name)
1038 {
1039  double a = -1, b, c;
1040  if (!strcmp(name,"Rectangular")) a = 1.0, b = 0.0, c = 0.0;
1041  if (!strcmp(name,"Hanning")) a = 0.5, b = 0.5, c = 0.0;
1042  if (!strcmp(name,"Hamming")) a = 0.54, b = 0.46, c = 0.0;
1043  if (!strcmp(name,"Blackman")) a = 0.42, b = 0.50, c = 0.08;
1044  if (a < 0)
1045  mfem_error("HypreSmoother::SetWindowByName : name not recognized!");
1046 
1047  SetWindowParameters(a, b, c);
1048 }
1049 
1050 void HypreSmoother::SetWindowParameters(double a, double b, double c)
1051 {
1052  window_params[0] = a;
1053  window_params[1] = b;
1054  window_params[2] = c;
1055 }
1056 
1058 {
1059  A = const_cast<HypreParMatrix *>(dynamic_cast<const HypreParMatrix *>(&op));
1060  if (A == NULL)
1061  mfem_error("HypreSmoother::SetOperator : not HypreParMatrix!");
1062 
1063  height = A->Height();
1064  width = A->Width();
1065 
1066  if (B) delete B;
1067  if (X) delete X;
1068  if (V) delete V;
1069  if (Z) delete Z;
1070  if (l1_norms)
1071  hypre_TFree(l1_norms);
1072  delete X0;
1073  delete X1;
1074 
1075  X1 = X0 = Z = V = B = X = NULL;
1076 
1077  if (type >= 1 && type <= 4)
1078  {
1079  hypre_ParCSRComputeL1Norms(*A, type, NULL, &l1_norms);
1080  }
1081  else if (type == 5)
1082  {
1083  l1_norms = hypre_CTAlloc(double, height);
1084  Vector ones(height), diag(l1_norms, height);
1085  ones = 1.0;
1086  A->Mult(ones, diag);
1087  type = 1;
1088  }
1089  else
1090  l1_norms = NULL;
1091 
1092  if (type == 16)
1093  {
1094  poly_scale = 1;
1095  hypre_ParCSRMaxEigEstimateCG(*A, poly_scale, 10,
1097  Z = new HypreParVector(*A);
1098  }
1099  else if (type == 1001 || type == 1002)
1100  {
1101  poly_scale = 0;
1102  hypre_ParCSRMaxEigEstimateCG(*A, poly_scale, 10,
1104 
1105  // The Taubin and FIR polynomials are defined on [0, 2]
1106  max_eig_est /= 2;
1107 
1108  // Compute window function, Chebyshev coefficients, and allocate temps.
1109  if (type == 1002)
1110  {
1111  // Temporaries for Chebyshev recursive evaluation
1112  Z = new HypreParVector(*A);
1113  X0 = new HypreParVector(*A);
1114  X1 = new HypreParVector(*A);
1115 
1117  }
1118  }
1119 }
1120 
1122 {
1123  if (fir_coeffs)
1124  delete [] fir_coeffs;
1125 
1126  fir_coeffs = new double[poly_order+1];
1127 
1128  double* window_coeffs = new double[poly_order+1];
1129  double* cheby_coeffs = new double[poly_order+1];
1130 
1131  double a = window_params[0];
1132  double b = window_params[1];
1133  double c = window_params[2];
1134  for (int i = 0; i <= poly_order; i++)
1135  {
1136  double t = (i*M_PI)/(poly_order+1);
1137  window_coeffs[i] = a + b*cos(t) +c*cos(2*t);
1138  }
1139 
1140  double k_pb = poly_fraction*max_eig;
1141  double theta_pb = acos(1.0 -0.5*k_pb);
1142  double sigma = 0.0;
1143  cheby_coeffs[0] = (theta_pb +sigma)/M_PI;
1144  for (int i = 1; i <= poly_order; i++)
1145  {
1146  double t = i*(theta_pb+sigma);
1147  cheby_coeffs[i] = 2.0*sin(t)/(i*M_PI);
1148  }
1149 
1150  for (int i = 0; i <= poly_order; i++)
1151  {
1152  fir_coeffs[i] = window_coeffs[i]*cheby_coeffs[i];
1153  }
1154 
1155  delete[] window_coeffs;
1156  delete[] cheby_coeffs;
1157 }
1158 
1160 {
1161  if (A == NULL)
1162  {
1163  mfem_error("HypreSmoother::Mult (...) : HypreParMatrix A is missing");
1164  return;
1165  }
1166 
1167  if (!iterative_mode)
1168  {
1169  if (type == 0 && relax_times == 1)
1170  {
1171  HYPRE_ParCSRDiagScale(NULL, *A, b, x);
1172  if (relax_weight != 1.0)
1173  x *= relax_weight;
1174  return;
1175  }
1176  x = 0.0;
1177  }
1178 
1179  if (V == NULL)
1180  V = new HypreParVector(*A);
1181 
1182  if (type == 1001)
1183  {
1184  for (int sweep = 0; sweep < relax_times; sweep++)
1185  {
1187  max_eig_est,
1188  x, *V);
1189  }
1190  }
1191  else if (type == 1002)
1192  {
1193  for (int sweep = 0; sweep < relax_times; sweep++)
1194  {
1195  ParCSRRelax_FIR(*A, b,
1196  max_eig_est,
1197  poly_order,
1198  fir_coeffs,
1199  x,
1200  *X0, *X1, *V, *Z);
1201  }
1202  }
1203  else
1204  {
1205  if (Z == NULL)
1206  hypre_ParCSRRelax(*A, b, type,
1209  x, *V, NULL);
1210  else
1211  hypre_ParCSRRelax(*A, b, type,
1214  x, *V, *Z);
1215  }
1216 }
1217 
1218 void HypreSmoother::Mult(const Vector &b, Vector &x) const
1219 {
1220  if (A == NULL)
1221  {
1222  mfem_error("HypreSmoother::Mult (...) : HypreParMatrix A is missing");
1223  return;
1224  }
1225  if (B == NULL)
1226  {
1227  B = new HypreParVector(A->GetComm(),
1228  A -> GetGlobalNumRows(),
1229  b.GetData(),
1230  A -> GetRowStarts());
1231  X = new HypreParVector(A->GetComm(),
1232  A -> GetGlobalNumCols(),
1233  x.GetData(),
1234  A -> GetColStarts());
1235  }
1236  else
1237  {
1238  B -> SetData(b.GetData());
1239  X -> SetData(x.GetData());
1240  }
1241 
1242  Mult(*B, *X);
1243 }
1244 
1246 {
1247  if (B) delete B;
1248  if (X) delete X;
1249  if (V) delete V;
1250  if (Z) delete Z;
1251  if (l1_norms)
1252  hypre_TFree(l1_norms);
1253  if (fir_coeffs)
1254  delete [] fir_coeffs;
1255  if (X0) delete X0;
1256  if (X1) delete X1;
1257 }
1258 
1259 
1261 {
1262  A = NULL;
1263  setup_called = 0;
1264  B = X = NULL;
1265 }
1266 
1268  : Solver(_A->Height(), _A->Width())
1269 {
1270  A = _A;
1271  setup_called = 0;
1272  B = X = NULL;
1273 }
1274 
1276 {
1277  if (A == NULL)
1278  {
1279  mfem_error("HypreSolver::Mult (...) : HypreParMatrix A is missing");
1280  return;
1281  }
1282  if (!setup_called)
1283  {
1284  SetupFcn()(*this, *A, b, x);
1285  setup_called = 1;
1286  }
1287 
1288  if (!iterative_mode)
1289  x = 0.0;
1290  SolveFcn()(*this, *A, b, x);
1291 }
1292 
1293 void HypreSolver::Mult(const Vector &b, Vector &x) const
1294 {
1295  if (A == NULL)
1296  {
1297  mfem_error("HypreSolver::Mult (...) : HypreParMatrix A is missing");
1298  return;
1299  }
1300  if (B == NULL)
1301  {
1302  B = new HypreParVector(A->GetComm(),
1303  A -> GetGlobalNumRows(),
1304  b.GetData(),
1305  A -> GetRowStarts());
1306  X = new HypreParVector(A->GetComm(),
1307  A -> GetGlobalNumCols(),
1308  x.GetData(),
1309  A -> GetColStarts());
1310  }
1311  else
1312  {
1313  B -> SetData(b.GetData());
1314  X -> SetData(x.GetData());
1315  }
1316 
1317  Mult(*B, *X);
1318 }
1319 
1321 {
1322  if (B) delete B;
1323  if (X) delete X;
1324 }
1325 
1326 
1328 {
1329  MPI_Comm comm;
1330 
1331  print_level = 0;
1332  iterative_mode = true;
1333 
1334  HYPRE_ParCSRMatrixGetComm(*A, &comm);
1335 
1336  HYPRE_ParCSRPCGCreate(comm, &pcg_solver);
1337 }
1338 
1339 void HyprePCG::SetTol(double tol)
1340 {
1341  HYPRE_ParCSRPCGSetTol(pcg_solver, tol);
1342 }
1343 
1344 void HyprePCG::SetMaxIter(int max_iter)
1345 {
1346  HYPRE_ParCSRPCGSetMaxIter(pcg_solver, max_iter);
1347 }
1348 
1349 void HyprePCG::SetLogging(int logging)
1350 {
1351  HYPRE_ParCSRPCGSetLogging(pcg_solver, logging);
1352 }
1353 
1354 void HyprePCG::SetPrintLevel(int print_lvl)
1355 {
1356  print_level = print_lvl;
1357  HYPRE_ParCSRPCGSetPrintLevel(pcg_solver, print_level);
1358 }
1359 
1361 {
1362  HYPRE_ParCSRPCGSetPrecond(pcg_solver,
1363  precond.SolveFcn(),
1364  precond.SetupFcn(),
1365  precond);
1366 }
1367 
1368 void HyprePCG::SetResidualConvergenceOptions(int res_frequency, double rtol)
1369 {
1370  HYPRE_PCGSetTwoNorm(pcg_solver, 1);
1371  if (res_frequency > 0)
1372  HYPRE_PCGSetRecomputeResidualP(pcg_solver, res_frequency);
1373  if (rtol > 0.0)
1374  HYPRE_PCGSetResidualTol(pcg_solver, rtol);
1375 }
1376 
1378 {
1379  int myid;
1380  int time_index = 0;
1381  int num_iterations;
1382  double final_res_norm;
1383  MPI_Comm comm;
1384 
1385  HYPRE_ParCSRMatrixGetComm(*A, &comm);
1386 
1387  if (!setup_called)
1388  {
1389  if (print_level > 0)
1390  {
1391  time_index = hypre_InitializeTiming("PCG Setup");
1392  hypre_BeginTiming(time_index);
1393  }
1394 
1395  HYPRE_ParCSRPCGSetup(pcg_solver, *A, b, x);
1396  setup_called = 1;
1397 
1398  if (print_level > 0)
1399  {
1400  hypre_EndTiming(time_index);
1401  hypre_PrintTiming("Setup phase times", comm);
1402  hypre_FinalizeTiming(time_index);
1403  hypre_ClearTiming();
1404  }
1405  }
1406 
1407  if (print_level > 0)
1408  {
1409  time_index = hypre_InitializeTiming("PCG Solve");
1410  hypre_BeginTiming(time_index);
1411  }
1412 
1413  if (!iterative_mode)
1414  x = 0.0;
1415 
1416  HYPRE_ParCSRPCGSolve(pcg_solver, *A, b, x);
1417 
1418  if (print_level > 0)
1419  {
1420  hypre_EndTiming(time_index);
1421  hypre_PrintTiming("Solve phase times", comm);
1422  hypre_FinalizeTiming(time_index);
1423  hypre_ClearTiming();
1424 
1425  HYPRE_ParCSRPCGGetNumIterations(pcg_solver, &num_iterations);
1426  HYPRE_ParCSRPCGGetFinalRelativeResidualNorm(pcg_solver,
1427  &final_res_norm);
1428 
1429  MPI_Comm_rank(comm, &myid);
1430 
1431  if (myid == 0)
1432  {
1433  cout << "PCG Iterations = " << num_iterations << endl
1434  << "Final PCG Relative Residual Norm = " << final_res_norm
1435  << endl;
1436  }
1437  }
1438 }
1439 
1441 {
1442  HYPRE_ParCSRPCGDestroy(pcg_solver);
1443 }
1444 
1445 
1447 {
1448  MPI_Comm comm;
1449 
1450  int k_dim = 50;
1451  int max_iter = 100;
1452  double tol = 1e-6;
1453 
1454  print_level = 0;
1455  iterative_mode = true;
1456 
1457  HYPRE_ParCSRMatrixGetComm(*A, &comm);
1458 
1459  HYPRE_ParCSRGMRESCreate(comm, &gmres_solver);
1460  HYPRE_ParCSRGMRESSetKDim(gmres_solver, k_dim);
1461  HYPRE_ParCSRGMRESSetMaxIter(gmres_solver, max_iter);
1462  HYPRE_ParCSRGMRESSetTol(gmres_solver, tol);
1463 }
1464 
1465 void HypreGMRES::SetTol(double tol)
1466 {
1467  HYPRE_ParCSRGMRESSetTol(gmres_solver, tol);
1468 }
1469 
1470 void HypreGMRES::SetMaxIter(int max_iter)
1471 {
1472  HYPRE_ParCSRGMRESSetMaxIter(gmres_solver, max_iter);
1473 }
1474 
1475 void HypreGMRES::SetKDim(int k_dim)
1476 {
1477  HYPRE_ParCSRGMRESSetKDim(gmres_solver, k_dim);
1478 }
1479 
1480 void HypreGMRES::SetLogging(int logging)
1481 {
1482  HYPRE_ParCSRGMRESSetLogging(gmres_solver, logging);
1483 }
1484 
1485 void HypreGMRES::SetPrintLevel(int print_lvl)
1486 {
1487  print_level = print_lvl;
1488  HYPRE_ParCSRGMRESSetPrintLevel(gmres_solver, print_level);
1489 }
1490 
1492 {
1493  HYPRE_ParCSRGMRESSetPrecond(gmres_solver,
1494  precond.SolveFcn(),
1495  precond.SetupFcn(),
1496  precond);
1497 }
1498 
1500 {
1501  int myid;
1502  int time_index = 0;
1503  int num_iterations;
1504  double final_res_norm;
1505  MPI_Comm comm;
1506 
1507  HYPRE_ParCSRMatrixGetComm(*A, &comm);
1508 
1509  if (!setup_called)
1510  {
1511  if (print_level > 0)
1512  {
1513  time_index = hypre_InitializeTiming("GMRES Setup");
1514  hypre_BeginTiming(time_index);
1515  }
1516 
1517  HYPRE_ParCSRGMRESSetup(gmres_solver, *A, b, x);
1518  setup_called = 1;
1519 
1520  if (print_level > 0)
1521  {
1522  hypre_EndTiming(time_index);
1523  hypre_PrintTiming("Setup phase times", comm);
1524  hypre_FinalizeTiming(time_index);
1525  hypre_ClearTiming();
1526  }
1527  }
1528 
1529  if (print_level > 0)
1530  {
1531  time_index = hypre_InitializeTiming("GMRES Solve");
1532  hypre_BeginTiming(time_index);
1533  }
1534 
1535  if (!iterative_mode)
1536  x = 0.0;
1537 
1538  HYPRE_ParCSRGMRESSolve(gmres_solver, *A, b, x);
1539 
1540  if (print_level > 0)
1541  {
1542  hypre_EndTiming(time_index);
1543  hypre_PrintTiming("Solve phase times", comm);
1544  hypre_FinalizeTiming(time_index);
1545  hypre_ClearTiming();
1546 
1547  HYPRE_ParCSRGMRESGetNumIterations(gmres_solver, &num_iterations);
1548  HYPRE_ParCSRGMRESGetFinalRelativeResidualNorm(gmres_solver,
1549  &final_res_norm);
1550 
1551  MPI_Comm_rank(comm, &myid);
1552 
1553  if (myid == 0)
1554  {
1555  cout << "GMRES Iterations = " << num_iterations << endl
1556  << "Final GMRES Relative Residual Norm = " << final_res_norm
1557  << endl;
1558  }
1559  }
1560 }
1561 
1563 {
1564  HYPRE_ParCSRGMRESDestroy(gmres_solver);
1565 }
1566 
1567 
1569 {
1570  MPI_Comm comm;
1571 
1572  int sai_max_levels = 1;
1573  double sai_threshold = 0.1;
1574  double sai_filter = 0.1;
1575  int sai_sym = 0;
1576  double sai_loadbal = 0.0;
1577  int sai_reuse = 0;
1578  int sai_logging = 1;
1579 
1580  HYPRE_ParCSRMatrixGetComm(A, &comm);
1581 
1582  HYPRE_ParaSailsCreate(comm, &sai_precond);
1583  HYPRE_ParaSailsSetParams(sai_precond, sai_threshold, sai_max_levels);
1584  HYPRE_ParaSailsSetFilter(sai_precond, sai_filter);
1585  HYPRE_ParaSailsSetSym(sai_precond, sai_sym);
1586  HYPRE_ParaSailsSetLoadbal(sai_precond, sai_loadbal);
1587  HYPRE_ParaSailsSetReuse(sai_precond, sai_reuse);
1588  HYPRE_ParaSailsSetLogging(sai_precond, sai_logging);
1589 }
1590 
1592 {
1593  HYPRE_ParaSailsSetSym(sai_precond, sym);
1594 }
1595 
1597 {
1598  HYPRE_ParaSailsDestroy(sai_precond);
1599 }
1600 
1601 
1603 {
1604  int coarsen_type = 10;
1605  int agg_levels = 1;
1606  int relax_type = 8;
1607  int relax_sweeps = 1;
1608  double theta = 0.25;
1609  int interp_type = 6;
1610  int Pmax = 4;
1611  int print_level = 1;
1612 
1613  HYPRE_BoomerAMGCreate(&amg_precond);
1614 
1615  HYPRE_BoomerAMGSetCoarsenType(amg_precond, coarsen_type);
1616  HYPRE_BoomerAMGSetAggNumLevels(amg_precond, agg_levels);
1617  HYPRE_BoomerAMGSetRelaxType(amg_precond, relax_type);
1618  HYPRE_BoomerAMGSetNumSweeps(amg_precond, relax_sweeps);
1619  HYPRE_BoomerAMGSetMaxLevels(amg_precond, 25);
1620  HYPRE_BoomerAMGSetTol(amg_precond, 0.0);
1621  HYPRE_BoomerAMGSetMaxIter(amg_precond, 1); // one V-cycle
1622  HYPRE_BoomerAMGSetStrongThreshold(amg_precond, theta);
1623  HYPRE_BoomerAMGSetInterpType(amg_precond, interp_type);
1624  HYPRE_BoomerAMGSetPMaxElmts(amg_precond, Pmax);
1625  HYPRE_BoomerAMGSetPrintLevel(amg_precond, print_level);
1626 }
1627 
1629 {
1630  HYPRE_BoomerAMGSetNumFunctions(amg_precond, dim);
1631 
1632  // More robust options with respect to convergence
1633  HYPRE_BoomerAMGSetAggNumLevels(amg_precond, 0);
1634  HYPRE_BoomerAMGSetStrongThreshold(amg_precond, 0.5);
1635 }
1636 
1638 {
1639  HYPRE_BoomerAMGDestroy(amg_precond);
1640 }
1641 
1642 
1644  : HypreSolver(&A)
1645 {
1646  int cycle_type = 13;
1647  int rlx_type = 2;
1648  int rlx_sweeps = 1;
1649  double rlx_weight = 1.0;
1650  double rlx_omega = 1.0;
1651  int amg_coarsen_type = 10;
1652  int amg_agg_levels = 1;
1653  int amg_rlx_type = 8;
1654  double theta = 0.25;
1655  int amg_interp_type = 6;
1656  int amg_Pmax = 4;
1657 
1658  int p = 1;
1659  if (edge_fespace->GetNE() > 0)
1660  p = edge_fespace->GetOrder(0);
1661  int dim = edge_fespace->GetMesh()->Dimension();
1662 
1663  HYPRE_AMSCreate(&ams);
1664 
1665  HYPRE_AMSSetDimension(ams, dim); // 2D H(div) and 3D H(curl) problems
1666  HYPRE_AMSSetTol(ams, 0.0);
1667  HYPRE_AMSSetMaxIter(ams, 1); // use as a preconditioner
1668  HYPRE_AMSSetCycleType(ams, cycle_type);
1669  HYPRE_AMSSetPrintLevel(ams, 1);
1670 
1671  // define the nodal linear finite element space associated with edge_fespace
1672  ParMesh *pmesh = (ParMesh *) edge_fespace->GetMesh();
1673  FiniteElementCollection *vert_fec = new H1_FECollection(p, dim);
1674  ParFiniteElementSpace *vert_fespace = new ParFiniteElementSpace(pmesh, vert_fec);
1675 
1676  // generate and set the vertex coordinates
1677  if (p == 1)
1678  {
1679  ParGridFunction x_coord(vert_fespace);
1680  ParGridFunction y_coord(vert_fespace);
1681  ParGridFunction z_coord(vert_fespace);
1682  double *coord;
1683  for (int i = 0; i < pmesh->GetNV(); i++)
1684  {
1685  coord = pmesh -> GetVertex(i);
1686  x_coord(i) = coord[0];
1687  y_coord(i) = coord[1];
1688  if (dim == 3)
1689  z_coord(i) = coord[2];
1690  }
1691  x = x_coord.ParallelAverage();
1692  y = y_coord.ParallelAverage();
1693  if (dim == 2)
1694  {
1695  z = NULL;
1696  HYPRE_AMSSetCoordinateVectors(ams, *x, *y, NULL);
1697  }
1698  else
1699  {
1700  z = z_coord.ParallelAverage();
1701  HYPRE_AMSSetCoordinateVectors(ams, *x, *y, *z);
1702  }
1703  }
1704  else
1705  {
1706  x = NULL;
1707  y = NULL;
1708  z = NULL;
1709  }
1710 
1711  // generate and set the discrete gradient
1713  grad = new ParDiscreteLinearOperator(vert_fespace, edge_fespace);
1715  grad->Assemble();
1716  grad->Finalize();
1717  G = grad->ParallelAssemble();
1718  HYPRE_AMSSetDiscreteGradient(ams, *G);
1719  delete grad;
1720 
1721  // generate and set the Nedelec interpolation matrices
1722  Pi = Pix = Piy = Piz = NULL;
1723  if (p > 1)
1724  {
1725  ParFiniteElementSpace *vert_fespace_d;
1726  if (cycle_type < 10)
1727  vert_fespace_d = new ParFiniteElementSpace(pmesh, vert_fec, dim,
1729  else
1730  vert_fespace_d = new ParFiniteElementSpace(pmesh, vert_fec, dim,
1732 
1734  id_ND = new ParDiscreteLinearOperator(vert_fespace_d, edge_fespace);
1736  id_ND->Assemble();
1737 
1738  if (cycle_type < 10)
1739  {
1740  id_ND->Finalize();
1741  Pi = id_ND->ParallelAssemble();
1742  }
1743  else
1744  {
1745  Array2D<HypreParMatrix *> Pi_blocks;
1746  id_ND->GetParBlocks(Pi_blocks);
1747  Pix = Pi_blocks(0,0);
1748  Piy = Pi_blocks(0,1);
1749  if (dim == 3)
1750  Piz = Pi_blocks(0,2);
1751  }
1752 
1753  delete id_ND;
1754 
1755  HYPRE_AMSSetInterpolations(ams, *Pi, *Pix, *Piy, *Piz);
1756 
1757  delete vert_fespace_d;
1758  }
1759 
1760  delete vert_fespace;
1761  delete vert_fec;
1762 
1763  // set additional AMS options
1764  HYPRE_AMSSetSmoothingOptions(ams, rlx_type, rlx_sweeps, rlx_weight, rlx_omega);
1765  HYPRE_AMSSetAlphaAMGOptions(ams, amg_coarsen_type, amg_agg_levels, amg_rlx_type,
1766  theta, amg_interp_type, amg_Pmax);
1767  HYPRE_AMSSetBetaAMGOptions(ams, amg_coarsen_type, amg_agg_levels, amg_rlx_type,
1768  theta, amg_interp_type, amg_Pmax);
1769 }
1770 
1772 {
1773  HYPRE_AMSDestroy(ams);
1774 
1775  delete x;
1776  delete y;
1777  delete z;
1778 
1779  delete G;
1780  delete Pi;
1781  delete Pix;
1782  delete Piy;
1783  delete Piz;
1784 }
1785 
1787  : HypreSolver(&A)
1788 {
1789  int cycle_type = 11;
1790  int rlx_type = 2;
1791  int rlx_sweeps = 1;
1792  double rlx_weight = 1.0;
1793  double rlx_omega = 1.0;
1794  int amg_coarsen_type = 10;
1795  int amg_agg_levels = 1;
1796  int amg_rlx_type = 6;
1797  double theta = 0.25;
1798  int amg_interp_type = 6;
1799  int amg_Pmax = 4;
1800  int ams_cycle_type = 14;
1801 
1802  int p = 1;
1803  if (face_fespace->GetNE() > 0)
1804  p = face_fespace->GetOrder(0);
1805 
1806  HYPRE_ADSCreate(&ads);
1807 
1808  HYPRE_ADSSetTol(ads, 0.0);
1809  HYPRE_ADSSetMaxIter(ads, 1); // use as a preconditioner
1810  HYPRE_ADSSetCycleType(ads, cycle_type);
1811  HYPRE_ADSSetPrintLevel(ads, 1);
1812 
1813  // define the nodal and edge finite element spaces associated with face_fespace
1814  ParMesh *pmesh = (ParMesh *) face_fespace->GetMesh();
1815  FiniteElementCollection *vert_fec = new H1_FECollection(p, 3);
1816  ParFiniteElementSpace *vert_fespace = new ParFiniteElementSpace(pmesh, vert_fec);
1817  FiniteElementCollection *edge_fec = new ND_FECollection(p, 3);
1818  ParFiniteElementSpace *edge_fespace = new ParFiniteElementSpace(pmesh, edge_fec);
1819 
1820  // generate and set the vertex coordinates
1821  if (p == 1)
1822  {
1823  ParGridFunction x_coord(vert_fespace);
1824  ParGridFunction y_coord(vert_fespace);
1825  ParGridFunction z_coord(vert_fespace);
1826  double *coord;
1827  for (int i = 0; i < pmesh->GetNV(); i++)
1828  {
1829  coord = pmesh -> GetVertex(i);
1830  x_coord(i) = coord[0];
1831  y_coord(i) = coord[1];
1832  z_coord(i) = coord[2];
1833  }
1834  x = x_coord.ParallelAverage();
1835  y = y_coord.ParallelAverage();
1836  z = z_coord.ParallelAverage();
1837  HYPRE_ADSSetCoordinateVectors(ads, *x, *y, *z);
1838  }
1839  else
1840  {
1841  x = NULL;
1842  y = NULL;
1843  z = NULL;
1844  }
1845 
1846  // generate and set the discrete curl
1848  curl = new ParDiscreteLinearOperator(edge_fespace, face_fespace);
1850  curl->Assemble();
1851  curl->Finalize();
1852  C = curl->ParallelAssemble();
1853  HYPRE_ADSSetDiscreteCurl(ads, *C);
1854  delete curl;
1855 
1856  // generate and set the discrete gradient
1858  grad = new ParDiscreteLinearOperator(vert_fespace, edge_fespace);
1860  grad->Assemble();
1861  grad->Finalize();
1862  G = grad->ParallelAssemble();
1863  HYPRE_ADSSetDiscreteGradient(ads, *G);
1864  delete grad;
1865 
1866  // generate and set the Nedelec and Raviart-Thomas interpolation matrices
1867  RT_Pi = RT_Pix = RT_Piy = RT_Piz = NULL;
1868  ND_Pi = ND_Pix = ND_Piy = ND_Piz = NULL;
1869  if (p > 1)
1870  {
1871  ParFiniteElementSpace *vert_fespace_d;
1872 
1873  if (ams_cycle_type < 10)
1874  vert_fespace_d = new ParFiniteElementSpace(pmesh, vert_fec, 3,
1876  else
1877  vert_fespace_d = new ParFiniteElementSpace(pmesh, vert_fec, 3,
1879 
1881  id_ND = new ParDiscreteLinearOperator(vert_fespace_d, edge_fespace);
1883  id_ND->Assemble();
1884 
1885  if (ams_cycle_type < 10)
1886  {
1887  id_ND->Finalize();
1888  ND_Pi = id_ND->ParallelAssemble();
1889  }
1890  else
1891  {
1892  Array2D<HypreParMatrix *> ND_Pi_blocks;
1893  id_ND->GetParBlocks(ND_Pi_blocks);
1894  ND_Pix = ND_Pi_blocks(0,0);
1895  ND_Piy = ND_Pi_blocks(0,1);
1896  ND_Piz = ND_Pi_blocks(0,2);
1897  }
1898 
1899  delete id_ND;
1900 
1901  if (cycle_type < 10 && ams_cycle_type > 10)
1902  {
1903  delete vert_fespace_d;
1904  vert_fespace_d = new ParFiniteElementSpace(pmesh, vert_fec, 3,
1906  }
1907  else if (cycle_type > 10 && ams_cycle_type < 10)
1908  {
1909  delete vert_fespace_d;
1910  vert_fespace_d = new ParFiniteElementSpace(pmesh, vert_fec, 3,
1912  }
1913 
1915  id_RT = new ParDiscreteLinearOperator(vert_fespace_d, face_fespace);
1917  id_RT->Assemble();
1918 
1919  if (cycle_type < 10)
1920  {
1921  id_RT->Finalize();
1922  RT_Pi = id_RT->ParallelAssemble();
1923  }
1924  else
1925  {
1926  Array2D<HypreParMatrix *> RT_Pi_blocks;
1927  id_RT->GetParBlocks(RT_Pi_blocks);
1928  RT_Pix = RT_Pi_blocks(0,0);
1929  RT_Piy = RT_Pi_blocks(0,1);
1930  RT_Piz = RT_Pi_blocks(0,2);
1931  }
1932 
1933  delete id_RT;
1934 
1935  HYPRE_ADSSetInterpolations(ads,
1936  *RT_Pi, *RT_Pix, *RT_Piy, *RT_Piz,
1937  *ND_Pi, *ND_Pix, *ND_Piy, *ND_Piz);
1938 
1939  delete vert_fespace_d;
1940  }
1941 
1942  delete vert_fec;
1943  delete vert_fespace;
1944  delete edge_fec;
1945  delete edge_fespace;
1946 
1947  // set additional ADS options
1948  HYPRE_ADSSetSmoothingOptions(ads, rlx_type, rlx_sweeps, rlx_weight, rlx_omega);
1949  HYPRE_ADSSetAMGOptions(ads, amg_coarsen_type, amg_agg_levels, amg_rlx_type,
1950  theta, amg_interp_type, amg_Pmax);
1951  HYPRE_ADSSetAMSOptions(ads, ams_cycle_type, amg_coarsen_type, amg_agg_levels,
1952  amg_rlx_type, theta, amg_interp_type, amg_Pmax);
1953 }
1954 
1956 {
1957  HYPRE_ADSDestroy(ads);
1958 
1959  delete x;
1960  delete y;
1961  delete z;
1962 
1963  delete G;
1964  delete C;
1965 
1966  delete RT_Pi;
1967  delete RT_Pix;
1968  delete RT_Piy;
1969  delete RT_Piz;
1970 
1971  delete ND_Pi;
1972  delete ND_Pix;
1973  delete ND_Piy;
1974  delete ND_Piz;
1975 }
1976 
1977 }
1978 
1979 #endif
virtual ~HypreBoomerAMG()
Definition: hypre.cpp:1637
void SetTol(double tol)
Definition: hypre.cpp:1339
void SetPreconditioner(HypreSolver &precond)
Set the hypre solver to be used as a preconditioner.
Definition: hypre.cpp:1491
int Size() const
Logical size of the array.
Definition: array.hpp:108
virtual int NumNonZeroElems() const
Returns the number of the nonzero elements in the matrix.
Definition: sparsemat.cpp:693
void DestroyCommPkg()
Definition: hypre.cpp:508
Vector()
Default constructor for Vector. Sets size = 0 and data = NULL.
Definition: vector.hpp:39
int * GetJ()
Definition: table.hpp:88
HypreParVector * X0
FIR Filter Temporary Vectors.
Definition: hypre.hpp:260
double min_eig_est
Minimal eigenvalue estimate for polynomial smoothing.
Definition: hypre.hpp:288
int setup_called
Was hypre&#39;s Setup function called already?
Definition: hypre.hpp:355
void SetSize(int s)
Resizes the vector if the new size is different.
Definition: vector.hpp:248
HypreParVector * B
Right-hand side and solution vector.
Definition: hypre.hpp:352
double window_params[3]
Paramters for windowing function of FIR filter.
Definition: hypre.hpp:290
void MakeDataOwner()
Definition: vector.hpp:70
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:1786
void SetWindowByName(const char *window_name)
Convenience function for setting canonical windowing parameters.
Definition: hypre.cpp:1037
HypreParMatrix(hypre_ParCSRMatrix *a)
Converts hypre&#39;s format to HypreParMatrix.
Definition: hypre.hpp:117
int Size() const
Returns the size of the vector.
Definition: vector.hpp:76
virtual ~HypreGMRES()
Definition: hypre.cpp:1562
Abstract parallel finite element space.
Definition: pfespace.hpp:28
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:1771
void SetTol(double tol)
Definition: hypre.cpp:1465
void Print(const char *fname, int offi=0, int offj=0)
Prints the locally owned rows in parallel.
Definition: hypre.cpp:716
double * GetData() const
Definition: vector.hpp:80
void SetPrintLevel(int print_lvl)
Definition: hypre.cpp:1485
void ScaleRows(const Vector &s)
Scale the local row i by s(i).
Definition: hypre.cpp:637
int Mult(HypreParVector &x, HypreParVector &y, double alpha=1.0, double beta=0.0)
Computes y = alpha * A * x + beta * y.
Definition: hypre.cpp:572
int Size_of_connections() const
Definition: table.hpp:72
double poly_fraction
Fraction of spectrum to smooth for polynomial relaxation.
Definition: hypre.hpp:274
void SetPrintLevel(int print_lvl)
Definition: hypre.cpp:1354
int poly_scale
Apply the polynomial smoother to A or D^{-1/2} A D^{-1/2}.
Definition: hypre.hpp:276
void AddDomainInterpolator(DiscreteInterpolator *di)
HypreParVector(MPI_Comm comm, int glob_size, int *col)
Definition: hypre.cpp:29
HypreParVector * X1
Definition: hypre.hpp:260
void SetWindowParameters(double a, double b, double c)
Set parameters for windowing function for FIR smoother.
Definition: hypre.cpp:1050
virtual void Mult(const HypreParVector &b, HypreParVector &x) const
Solve Ax=b with hypre&#39;s GMRES.
Definition: hypre.cpp:1499
void SetSystemsOptions(int dim)
Definition: hypre.cpp:1628
virtual void Mult(const HypreParVector &b, HypreParVector &x) const
Relax the linear system Ax=b.
Definition: hypre.cpp:1159
int GetNE() const
Returns number of elements in the mesh.
Definition: fespace.hpp:169
void SetSymmetry(int sym)
Definition: hypre.cpp:1591
virtual ~HypreParaSails()
Definition: hypre.cpp:1596
void Print(const char *fname)
Prints the locally owned rows in parallel.
Definition: hypre.cpp:155
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:631
MPI_Comm GetComm()
MPI communicator.
Definition: hypre.hpp:156
double * l1_norms
l1 norms of the rows of A
Definition: hypre.hpp:284
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:1349
Mesh * GetMesh() const
Returns the mesh.
Definition: fespace.hpp:132
virtual ~HypreSolver()
Definition: hypre.cpp:1320
HyprePCG(HypreParMatrix &_A)
Definition: hypre.cpp:1327
virtual HYPRE_PtrToParSolverFcn SolveFcn() const =0
hypre&#39;s internal Solve function
void SetKDim(int dim)
Definition: hypre.cpp:1475
void SetResidualConvergenceOptions(int res_frequency=-1, double rtol=0.0)
Definition: hypre.cpp:1368
virtual void Mult(const HypreParVector &b, HypreParVector &x) const
Solve the linear system Ax=b.
Definition: hypre.cpp:1275
HypreParMatrix * RAP(HypreParMatrix *A, HypreParMatrix *P)
Returns the matrix P^t * A * P.
Definition: hypre.cpp:781
int Randomize(int seed)
Set random values.
Definition: hypre.cpp:150
double relax_weight
Damping coefficient (usually &lt;= 1)
Definition: hypre.hpp:268
void EliminateBC(HypreParMatrix &A, HypreParMatrix &Ae, Array< int > &ess_dof_list, HypreParVector &x, HypreParVector &b)
Definition: hypre.cpp:827
int Size() const
Returns the number of TYPE I elements.
Definition: table.hpp:66
int Dimension() const
Definition: mesh.hpp:417
HypreParMatrix * A
The linear system matrix.
Definition: hypre.hpp:254
double * GetData() const
Return element data.
Definition: sparsemat.hpp:101
HypreAMS(HypreParMatrix &A, ParFiniteElementSpace *edge_fespace)
Definition: hypre.cpp:1643
int * GetI() const
Return the array I.
Definition: sparsemat.hpp:97
virtual void SetOperator(const Operator &op)
Definition: hypre.cpp:1057
virtual HYPRE_PtrToParSolverFcn SetupFcn() const =0
hypre&#39;s internal Setup function
void SetLogging(int logging)
Definition: hypre.cpp:1480
virtual ~HypreADS()
Definition: hypre.cpp:1955
void SetMaxIter(int max_iter)
Definition: hypre.cpp:1344
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:896
void SetSOROptions(double relax_weight, double omega)
Set SOR-related parameters.
Definition: hypre.cpp:1017
Wrapper for hypre&#39;s parallel vector class.
Definition: hypre.hpp:38
HypreParVector * X
Definition: hypre.hpp:352
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:859
int * GetColStarts() const
Definition: hypre.hpp:197
Vector * GlobalVector()
Returns the global vector in each processor.
Definition: hypre.cpp:117
HypreParVector * X
Definition: hypre.hpp:256
void mfem_error(const char *msg)
Definition: error.cpp:23
int GetOrder(int i) const
Returns the order of the i&#39;th finite element.
Definition: fespace.cpp:23
int GetGlobalNumCols() const
Definition: hypre.hpp:193
int GetNumRows() const
Returns the number of rows in the diagonal block of the ParCSRMatrix.
Definition: hypre.hpp:184
double max_eig_est
Maximal eigenvalue estimate for polynomial smoothing.
Definition: hypre.hpp:286
void SetDataAndSize(double *d, int s)
Definition: vector.hpp:64
HypreParMatrix * Transpose()
Returns the transpose of *this.
Definition: hypre.cpp:561
void SetFIRCoefficients(double max_eig)
Compute window and Chebyshev coefficients for given polynomial order.
Definition: hypre.cpp:1121
int GetNV() const
Definition: mesh.hpp:393
double * fir_coeffs
Combined coefficients for windowing and Chebyshev polynomials.
Definition: hypre.hpp:293
double InnerProduct(HypreParVector *x, HypreParVector *y)
Definition: hypre.cpp:167
virtual ~HypreSmoother()
Definition: hypre.cpp:1245
int * GetRowStarts() const
Definition: hypre.hpp:195
void SetData(double *_data)
Definition: hypre.cpp:145
void InvScaleRows(const Vector &s)
Scale the local row i by 1./s(i)
Definition: hypre.cpp:665
void SetTaubinOptions(double lambda, double mu, int iter)
Set parameters for Taubin&#39;s lambda-mu method.
Definition: hypre.cpp:1029
HypreParVector * B
Right-hand side and solution vectors.
Definition: hypre.hpp:256
void SetCommPkg(hypre_ParCSRCommPkg *comm_pkg)
Definition: hypre.cpp:490
virtual ~HyprePCG()
Definition: hypre.cpp:1440
int GetNumCols() const
Returns the number of columns in the diagonal block of the ParCSRMatrix.
Definition: hypre.hpp:188
void SetPreconditioner(HypreSolver &precond)
Set the hypre solver to be used as a preconditioner.
Definition: hypre.cpp:1360
void SetPolyOptions(int poly_order, double poly_fraction)
Set parameters for polynomial smoothing.
Definition: hypre.cpp:1023
int relax_times
Number of relaxation sweeps.
Definition: hypre.hpp:266
void operator*=(double s)
Scale all entries by s: A_scaled = s*A.
Definition: hypre.cpp:697
Abstract class for hypre&#39;s solvers and preconditioners.
Definition: hypre.hpp:345
int GetGlobalNumRows() const
Definition: hypre.hpp:191
HypreParVector & operator=(double d)
Set constant values.
Definition: hypre.cpp:127
Arbitrary order H(curl)-conforming Nedelec finite elements.
Definition: fe_coll.hpp:158
Vector data type.
Definition: vector.hpp:29
void GetDiag(Vector &diag)
Get the diagonal of the matrix.
Definition: hypre.cpp:546
HypreGMRES(HypreParMatrix &_A)
Definition: hypre.cpp:1446
virtual ~HypreParMatrix()
Calls hypre&#39;s destroy function.
Definition: hypre.cpp:731
void GetParBlocks(Array2D< HypreParMatrix * > &blocks) const
hypre_ParCSRMatrix * StealData()
Changes the ownership of the the matrix.
Definition: hypre.cpp:539
HypreParMatrix * A
The linear system matrix.
Definition: hypre.hpp:349
int * GetI()
Definition: table.hpp:87
Arbitrary order H1-conforming (continuous) finite elements.
Definition: fe_coll.hpp:52
void SetMaxIter(int max_iter)
Definition: hypre.cpp:1470
double omega
SOR parameter (usually in (0,2))
Definition: hypre.hpp:270
Base class for solvers.
Definition: operator.hpp:102
Class for parallel grid function.
Definition: pgridfunc.hpp:31
double * data
Definition: vector.hpp:34
HypreParMatrix * ParallelAssemble(SparseMatrix *m)
Abstract operator.
Definition: operator.hpp:21
Wrapper for hypre&#39;s ParCSR matrix class.
Definition: hypre.hpp:103
~HypreParVector()
Calls hypre&#39;s destroy function.
Definition: hypre.cpp:160
virtual void Assemble(int skip_zeros=1)
int * GetJ() const
Return the array J.
Definition: sparsemat.hpp:99
virtual void Mult(const HypreParVector &b, HypreParVector &x) const
Solve Ax=b with hypre&#39;s PCG.
Definition: hypre.cpp:1377
Class for parallel meshes.
Definition: pmesh.hpp:27
HypreParVector * Z
Definition: hypre.hpp:258
void Read(MPI_Comm comm, const char *fname)
Reads the matrix from a file.
Definition: hypre.cpp:721
void ParallelAverage(Vector &tv) const
Returns the vector averaged on the true dofs.
Definition: pgridfunc.cpp:110
HypreParVector * V
Temporary vectors.
Definition: hypre.hpp:258
void SetType(HypreSmoother::Type type, int relax_times=1)
Set the relaxation type and number of sweeps.
Definition: hypre.cpp:1011
HypreParaSails(HypreParMatrix &A)
Definition: hypre.cpp:1568
int poly_order
Order of the smoothing polynomial.
Definition: hypre.hpp:272
double lambda
Taubin&#39;s lambda-mu method parameters.
Definition: hypre.hpp:279
HypreParMatrix * ParMult(HypreParMatrix *A, HypreParMatrix *B)
Returns the matrix A * B.
Definition: hypre.cpp:771
HypreBoomerAMG(HypreParMatrix &A)
Definition: hypre.cpp:1602