12 #include "../config/config.hpp"
17 #include "../fem/fem.hpp"
30 template<
typename TargetT,
typename SourceT>
31 static TargetT *DuplicateAs(
const SourceT *array,
int size,
32 bool cplusplus =
true)
34 TargetT *target_array = cplusplus ?
new TargetT[size]
35 : hypre_TAlloc(TargetT, size);
36 for (
int i = 0; i < size; i++)
38 target_array[i] = array[i];
43 inline void HypreParVector::_SetDataAndSize_()
45 SetDataAndSize(hypre_VectorData(hypre_ParVectorLocalVector(x)),
47 hypre_VectorSize(hypre_ParVectorLocalVector(x))));
50 HypreParVector::HypreParVector(MPI_Comm comm, HYPRE_Int glob_size,
53 x = hypre_ParVectorCreate(comm,glob_size,col);
54 hypre_ParVectorInitialize(x);
55 hypre_ParVectorSetPartitioningOwner(x,0);
57 hypre_ParVectorSetDataOwner(x,1);
58 hypre_SeqVectorSetDataOwner(hypre_ParVectorLocalVector(x),1);
64 double *_data, HYPRE_Int *col) :
Vector()
66 x = hypre_ParVectorCreate(comm,glob_size,col);
67 hypre_ParVectorSetDataOwner(x,1);
68 hypre_SeqVectorSetDataOwner(hypre_ParVectorLocalVector(x),0);
69 hypre_ParVectorSetPartitioningOwner(x,0);
71 hypre_VectorData(hypre_ParVectorLocalVector(x)) = &tmp;
74 hypre_ParVectorInitialize(x);
76 hypre_VectorData(hypre_ParVectorLocalVector(x)) = _data;
83 x = hypre_ParVectorCreate(y.x -> comm, y.x -> global_size,
85 hypre_ParVectorInitialize(x);
86 hypre_ParVectorSetPartitioningOwner(x,0);
87 hypre_ParVectorSetDataOwner(x,1);
88 hypre_SeqVectorSetDataOwner(hypre_ParVectorLocalVector(x),1);
98 x = hypre_ParVectorInDomainOf(const_cast<HypreParMatrix&>(A));
102 x = hypre_ParVectorInRangeOf(const_cast<HypreParMatrix&>(A));
110 x = (hypre_ParVector *) y;
119 hypre_ParVectorInitialize(x);
120 hypre_ParVectorSetPartitioningOwner(x,0);
122 hypre_ParVectorSetDataOwner(x,1);
123 hypre_SeqVectorSetDataOwner(hypre_ParVectorLocalVector(x),1);
128 HypreParVector::operator hypre_ParVector*()
const
133 #ifndef HYPRE_PAR_VECTOR_STRUCT
134 HypreParVector::operator HYPRE_ParVector()
const
136 return (HYPRE_ParVector) x;
142 hypre_Vector *hv = hypre_ParVectorToVectorAll(*
this);
143 Vector *v =
new Vector(hv->data, internal::to_int(hv->size));
145 hypre_SeqVectorSetDataOwner(hv,0);
146 hypre_SeqVectorDestroy(hv);
152 hypre_ParVectorSetConstantValues(x,d);
159 if (size != y.
Size())
165 for (
int i = 0; i <
size; i++)
174 Vector::data = hypre_VectorData(hypre_ParVectorLocalVector(x)) = _data;
179 return hypre_ParVectorSetRandomValues(x,seed);
184 hypre_ParVectorPrint(x,fname);
191 hypre_ParVectorDestroy(x);
198 return hypre_ParVectorInnerProd(*x, *y);
203 return hypre_ParVectorInnerProd(x, y);
212 double loc_norm = vec.
Norml1();
213 MPI_Allreduce(&loc_norm, &norm, 1, MPI_DOUBLE, MPI_SUM, comm);
217 double loc_norm = vec*vec;
218 MPI_Allreduce(&loc_norm, &norm, 1, MPI_DOUBLE, MPI_SUM, comm);
221 if (p < std::numeric_limits<double>::infinity())
224 for (
int i = 0; i < vec.
Size(); i++)
226 sum += pow(fabs(vec(i)), p);
228 MPI_Allreduce(&sum, &norm, 1, MPI_DOUBLE, MPI_SUM, comm);
229 norm = pow(norm, 1.0/p);
234 MPI_Allreduce(&loc_norm, &norm, 1, MPI_DOUBLE, MPI_MAX, comm);
240 void HypreParMatrix::Init()
244 diagOwner = offdOwner = colMapOwner = -1;
254 char HypreParMatrix::CopyCSR(
SparseMatrix *csr, hypre_CSRMatrix *hypre_csr)
256 hypre_CSRMatrixData(hypre_csr) = csr->
GetData();
258 hypre_CSRMatrixI(hypre_csr) = csr->
GetI();
259 hypre_CSRMatrixJ(hypre_csr) = csr->
GetJ();
263 hypre_CSRMatrixI(hypre_csr) =
264 DuplicateAs<HYPRE_Int>(csr->
GetI(), csr->
Height()+1);
265 hypre_CSRMatrixJ(hypre_csr) =
272 char HypreParMatrix::CopyBoolCSR(Table *bool_csr, hypre_CSRMatrix *hypre_csr)
274 int nnz = bool_csr->Size_of_connections();
275 double *data =
new double[nnz];
276 for (
int i = 0; i < nnz; i++)
280 hypre_CSRMatrixData(hypre_csr) = data;
282 hypre_CSRMatrixI(hypre_csr) = bool_csr->GetI();
283 hypre_CSRMatrixJ(hypre_csr) = bool_csr->GetJ();
287 hypre_CSRMatrixI(hypre_csr) =
288 DuplicateAs<HYPRE_Int>(bool_csr->GetI(), bool_csr->Size()+1);
289 hypre_CSRMatrixJ(hypre_csr) =
290 DuplicateAs<HYPRE_Int>(bool_csr->GetJ(), nnz);
296 void HypreParMatrix::CopyCSR_J(hypre_CSRMatrix *hypre_csr,
int *J)
298 HYPRE_Int nnz = hypre_CSRMatrixNumNonzeros(hypre_csr);
299 for (HYPRE_Int j = 0; j < nnz; j++)
301 J[j] = int(hypre_CSRMatrixJ(hypre_csr)[j]);
308 :
Operator(diag->Height(), diag->Width())
311 A = hypre_ParCSRMatrixCreate(comm, glob_size, glob_size, row_starts,
313 hypre_ParCSRMatrixSetDataOwner(A,1);
314 hypre_ParCSRMatrixSetRowStartsOwner(A,0);
315 hypre_ParCSRMatrixSetColStartsOwner(A,0);
317 hypre_CSRMatrixSetDataOwner(A->diag,0);
318 diagOwner = CopyCSR(diag, A->diag);
319 hypre_CSRMatrixSetRownnz(A->diag);
321 hypre_CSRMatrixSetDataOwner(A->offd,1);
322 hypre_CSRMatrixI(A->offd) = hypre_CTAlloc(HYPRE_Int, diag->
Height()+1);
329 hypre_ParCSRMatrixSetNumNonzeros(A);
332 hypre_CSRMatrixReorder(hypre_ParCSRMatrixDiag(A));
334 CopyCSR_J(A->diag, diag->
GetJ());
337 hypre_MatvecCommPkgCreate(A);
342 HYPRE_Int global_num_rows,
343 HYPRE_Int global_num_cols,
344 HYPRE_Int *row_starts, HYPRE_Int *col_starts,
346 :
Operator(diag->Height(), diag->Width())
349 A = hypre_ParCSRMatrixCreate(comm, global_num_rows, global_num_cols,
350 row_starts, col_starts,
352 hypre_ParCSRMatrixSetDataOwner(A,1);
353 hypre_ParCSRMatrixSetRowStartsOwner(A,0);
354 hypre_ParCSRMatrixSetColStartsOwner(A,0);
356 hypre_CSRMatrixSetDataOwner(A->diag,0);
357 diagOwner = CopyCSR(diag, A->diag);
358 hypre_CSRMatrixSetRownnz(A->diag);
360 hypre_CSRMatrixSetDataOwner(A->offd,1);
361 hypre_CSRMatrixI(A->offd) = hypre_CTAlloc(HYPRE_Int, diag->
Height()+1);
363 hypre_ParCSRMatrixSetNumNonzeros(A);
366 if (row_starts == col_starts)
368 hypre_CSRMatrixReorder(hypre_ParCSRMatrixDiag(A));
370 CopyCSR_J(A->diag, diag->
GetJ());
374 hypre_MatvecCommPkgCreate(A);
379 HYPRE_Int global_num_rows,
380 HYPRE_Int global_num_cols,
381 HYPRE_Int *row_starts, HYPRE_Int *col_starts,
384 :
Operator(diag->Height(), diag->Width())
387 A = hypre_ParCSRMatrixCreate(comm, global_num_rows, global_num_cols,
388 row_starts, col_starts,
391 hypre_ParCSRMatrixSetDataOwner(A,1);
392 hypre_ParCSRMatrixSetRowStartsOwner(A,0);
393 hypre_ParCSRMatrixSetColStartsOwner(A,0);
395 hypre_CSRMatrixSetDataOwner(A->diag,0);
396 diagOwner = CopyCSR(diag, A->diag);
397 hypre_CSRMatrixSetRownnz(A->diag);
399 hypre_CSRMatrixSetDataOwner(A->offd,0);
400 offdOwner = CopyCSR(offd, A->offd);
401 hypre_CSRMatrixSetRownnz(A->offd);
403 hypre_ParCSRMatrixColMapOffd(A) = cmap;
407 hypre_ParCSRMatrixSetNumNonzeros(A);
410 if (row_starts == col_starts)
412 hypre_CSRMatrixReorder(hypre_ParCSRMatrixDiag(A));
414 CopyCSR_J(A->diag, diag->
GetJ());
418 hypre_MatvecCommPkgCreate(A);
424 HYPRE_Int global_num_rows, HYPRE_Int global_num_cols,
425 HYPRE_Int *row_starts, HYPRE_Int *col_starts,
426 HYPRE_Int *diag_i, HYPRE_Int *diag_j,
double *diag_data,
427 HYPRE_Int *offd_i, HYPRE_Int *offd_j,
double *offd_data,
428 HYPRE_Int offd_num_cols, HYPRE_Int *offd_col_map)
431 A = hypre_ParCSRMatrixCreate(comm, global_num_rows, global_num_cols,
432 row_starts, col_starts, offd_num_cols, 0, 0);
433 hypre_ParCSRMatrixSetDataOwner(A,1);
434 hypre_ParCSRMatrixSetRowStartsOwner(A,0);
435 hypre_ParCSRMatrixSetColStartsOwner(A,0);
437 HYPRE_Int local_num_rows = hypre_CSRMatrixNumRows(A->diag);
439 hypre_CSRMatrixSetDataOwner(A->diag,0);
440 hypre_CSRMatrixI(A->diag) = diag_i;
441 hypre_CSRMatrixJ(A->diag) = diag_j;
442 hypre_CSRMatrixData(A->diag) = diag_data;
443 hypre_CSRMatrixNumNonzeros(A->diag) = diag_i[local_num_rows];
444 hypre_CSRMatrixSetRownnz(A->diag);
448 hypre_CSRMatrixSetDataOwner(A->offd,0);
449 hypre_CSRMatrixI(A->offd) = offd_i;
450 hypre_CSRMatrixJ(A->offd) = offd_j;
451 hypre_CSRMatrixData(A->offd) = offd_data;
452 hypre_CSRMatrixNumNonzeros(A->offd) = offd_i[local_num_rows];
453 hypre_CSRMatrixSetRownnz(A->offd);
457 hypre_ParCSRMatrixColMapOffd(A) = offd_col_map;
461 hypre_ParCSRMatrixSetNumNonzeros(A);
464 if (row_starts == col_starts)
466 hypre_CSRMatrixReorder(hypre_ParCSRMatrixDiag(A));
469 hypre_MatvecCommPkgCreate(A);
477 HYPRE_Int *row_starts, HYPRE_Int *col_starts,
480 MFEM_ASSERT(sm_a != NULL,
"invalid input");
481 MFEM_VERIFY(!HYPRE_AssumedPartitionCheck(),
482 "this method can not be used with assumed partition");
486 hypre_CSRMatrix *csr_a;
487 csr_a = hypre_CSRMatrixCreate(sm_a ->
Height(), sm_a ->
Width(),
488 sm_a -> NumNonZeroElems());
490 hypre_CSRMatrixSetDataOwner(csr_a,0);
491 CopyCSR(sm_a, csr_a);
492 hypre_CSRMatrixSetRownnz(csr_a);
494 A = hypre_CSRMatrixToParCSRMatrix(comm, csr_a, row_starts, col_starts);
497 delete [] hypre_CSRMatrixI(csr_a);
498 delete [] hypre_CSRMatrixJ(csr_a);
500 hypre_CSRMatrixI(csr_a) = NULL;
501 hypre_CSRMatrixDestroy(csr_a);
507 if (row_starts == col_starts)
509 hypre_CSRMatrixReorder(hypre_ParCSRMatrixDiag(A));
512 hypre_MatvecCommPkgCreate(A);
517 HYPRE_Int global_num_rows,
518 HYPRE_Int global_num_cols,
519 HYPRE_Int *row_starts, HYPRE_Int *col_starts,
524 A = hypre_ParCSRMatrixCreate(comm, global_num_rows, global_num_cols,
525 row_starts, col_starts, 0, nnz, 0);
526 hypre_ParCSRMatrixSetDataOwner(A,1);
527 hypre_ParCSRMatrixSetRowStartsOwner(A,0);
528 hypre_ParCSRMatrixSetColStartsOwner(A,0);
530 hypre_CSRMatrixSetDataOwner(A->diag,0);
531 diagOwner = CopyBoolCSR(diag, A->diag);
532 hypre_CSRMatrixSetRownnz(A->diag);
534 hypre_CSRMatrixSetDataOwner(A->offd,1);
535 hypre_CSRMatrixI(A->offd) = hypre_CTAlloc(HYPRE_Int, diag->
Size()+1);
537 hypre_ParCSRMatrixSetNumNonzeros(A);
540 if (row_starts == col_starts)
542 hypre_CSRMatrixReorder(hypre_ParCSRMatrixDiag(A));
544 CopyCSR_J(A->diag, diag->
GetJ());
548 hypre_MatvecCommPkgCreate(A);
556 HYPRE_Int *row, HYPRE_Int *col,
557 HYPRE_Int *i_diag, HYPRE_Int *j_diag,
558 HYPRE_Int *i_offd, HYPRE_Int *j_offd,
559 HYPRE_Int *cmap, HYPRE_Int cmap_size)
561 HYPRE_Int diag_nnz, offd_nnz;
564 if (HYPRE_AssumedPartitionCheck())
566 diag_nnz = i_diag[row[1]-row[0]];
567 offd_nnz = i_offd[row[1]-row[0]];
569 A = hypre_ParCSRMatrixCreate(comm, row[2], col[2], row, col,
570 cmap_size, diag_nnz, offd_nnz);
574 diag_nnz = i_diag[row[
id+1]-row[id]];
575 offd_nnz = i_offd[row[
id+1]-row[id]];
577 A = hypre_ParCSRMatrixCreate(comm, row[np], col[np], row, col,
578 cmap_size, diag_nnz, offd_nnz);
581 hypre_ParCSRMatrixSetDataOwner(A,1);
582 hypre_ParCSRMatrixSetRowStartsOwner(A,0);
583 hypre_ParCSRMatrixSetColStartsOwner(A,0);
587 double *a_diag =
new double[diag_nnz];
588 for (i = 0; i < diag_nnz; i++)
593 double *a_offd =
new double[offd_nnz];
594 for (i = 0; i < offd_nnz; i++)
599 hypre_CSRMatrixSetDataOwner(A->diag,0);
600 hypre_CSRMatrixI(A->diag) = i_diag;
601 hypre_CSRMatrixJ(A->diag) = j_diag;
602 hypre_CSRMatrixData(A->diag) = a_diag;
603 hypre_CSRMatrixSetRownnz(A->diag);
607 hypre_CSRMatrixSetDataOwner(A->offd,0);
608 hypre_CSRMatrixI(A->offd) = i_offd;
609 hypre_CSRMatrixJ(A->offd) = j_offd;
610 hypre_CSRMatrixData(A->offd) = a_offd;
611 hypre_CSRMatrixSetRownnz(A->offd);
615 hypre_ParCSRMatrixColMapOffd(A) = cmap;
619 hypre_ParCSRMatrixSetNumNonzeros(A);
624 hypre_CSRMatrixReorder(hypre_ParCSRMatrixDiag(A));
627 hypre_MatvecCommPkgCreate(A);
636 HYPRE_Int glob_ncols,
int *I, HYPRE_Int *J,
637 double *data, HYPRE_Int *rows, HYPRE_Int *cols)
643 HYPRE_Int my_col_start, my_col_end;
644 if (HYPRE_AssumedPartitionCheck())
647 my_col_start = cols[0];
648 my_col_end = cols[1];
653 MPI_Comm_rank(comm, &myid);
654 MPI_Comm_size(comm, &part_size);
656 my_col_start = cols[myid];
657 my_col_end = cols[myid+1];
661 HYPRE_Int *row_starts, *col_starts;
664 row_starts = col_starts = hypre_TAlloc(HYPRE_Int, part_size);
665 for (
int i = 0; i < part_size; i++)
667 row_starts[i] = rows[i];
672 row_starts = hypre_TAlloc(HYPRE_Int, part_size);
673 col_starts = hypre_TAlloc(HYPRE_Int, part_size);
674 for (
int i = 0; i < part_size; i++)
676 row_starts[i] = rows[i];
677 col_starts[i] = cols[i];
683 HYPRE_Int diag_nnz = 0, offd_nnz = 0, offd_num_cols = 0;
684 map<HYPRE_Int, HYPRE_Int> offd_map;
685 for (HYPRE_Int j = 0, loc_nnz = I[nrows]; j < loc_nnz; j++)
687 HYPRE_Int glob_col = J[j];
688 if (my_col_start <= glob_col && glob_col < my_col_end)
694 offd_map.insert(pair<const HYPRE_Int, HYPRE_Int>(glob_col, -1));
699 for (map<HYPRE_Int, HYPRE_Int>::iterator it = offd_map.begin();
700 it != offd_map.end(); ++it)
702 it->second = offd_num_cols++;
706 A = hypre_ParCSRMatrixCreate(comm, glob_nrows, glob_ncols,
707 row_starts, col_starts, offd_num_cols,
709 hypre_ParCSRMatrixInitialize(A);
711 HYPRE_Int *diag_i, *diag_j, *offd_i, *offd_j, *offd_col_map;
712 double *diag_data, *offd_data;
715 diag_data = A->diag->data;
718 offd_data = A->offd->data;
719 offd_col_map = A->col_map_offd;
721 diag_nnz = offd_nnz = 0;
722 for (HYPRE_Int i = 0, j = 0; i < nrows; i++)
724 diag_i[i] = diag_nnz;
725 offd_i[i] = offd_nnz;
726 for (HYPRE_Int j_end = I[i+1]; j < j_end; j++)
728 HYPRE_Int glob_col = J[j];
729 if (my_col_start <= glob_col && glob_col < my_col_end)
731 diag_j[diag_nnz] = glob_col - my_col_start;
732 diag_data[diag_nnz] = data[j];
737 offd_j[offd_nnz] = offd_map[glob_col];
738 offd_data[offd_nnz] = data[j];
743 diag_i[nrows] = diag_nnz;
744 offd_i[nrows] = offd_nnz;
745 for (map<HYPRE_Int, HYPRE_Int>::iterator it = offd_map.begin();
746 it != offd_map.end(); ++it)
748 offd_col_map[it->second] = it->first;
751 hypre_ParCSRMatrixSetNumNonzeros(A);
753 if (row_starts == col_starts)
755 hypre_CSRMatrixReorder(hypre_ParCSRMatrixDiag(A));
757 hypre_MatvecCommPkgCreate(A);
778 MFEM_ASSERT(diagOwner == -1 && offdOwner == -1 && colMapOwner == -1,
"");
779 MFEM_ASSERT(ParCSROwner,
"");
780 hypre_ParCSRMatrix *R = A;
789 if (!A || hypre_ParCSRMatrixOwnsRowStarts(A) ||
790 (hypre_ParCSRMatrixRowStarts(A) == hypre_ParCSRMatrixColStarts(A) &&
791 hypre_ParCSRMatrixOwnsColStarts(A)))
797 if (HYPRE_AssumedPartitionCheck())
803 MPI_Comm_size(hypre_ParCSRMatrixComm(A), &row_starts_size);
807 HYPRE_Int *old_row_starts = hypre_ParCSRMatrixRowStarts(A);
808 HYPRE_Int *new_row_starts = hypre_CTAlloc(HYPRE_Int, row_starts_size);
809 for (
int i = 0; i < row_starts_size; i++)
811 new_row_starts[i] = old_row_starts[i];
814 hypre_ParCSRMatrixRowStarts(A) = new_row_starts;
815 hypre_ParCSRMatrixOwnsRowStarts(A) = 1;
817 if (hypre_ParCSRMatrixColStarts(A) == old_row_starts)
819 hypre_ParCSRMatrixColStarts(A) = new_row_starts;
820 hypre_ParCSRMatrixOwnsColStarts(A) = 0;
826 if (!A || hypre_ParCSRMatrixOwnsColStarts(A) ||
827 (hypre_ParCSRMatrixRowStarts(A) == hypre_ParCSRMatrixColStarts(A) &&
828 hypre_ParCSRMatrixOwnsRowStarts(A)))
834 if (HYPRE_AssumedPartitionCheck())
840 MPI_Comm_size(hypre_ParCSRMatrixComm(A), &col_starts_size);
844 HYPRE_Int *old_col_starts = hypre_ParCSRMatrixColStarts(A);
845 HYPRE_Int *new_col_starts = hypre_CTAlloc(HYPRE_Int, col_starts_size);
846 for (
int i = 0; i < col_starts_size; i++)
848 new_col_starts[i] = old_col_starts[i];
851 hypre_ParCSRMatrixColStarts(A) = new_col_starts;
853 if (hypre_ParCSRMatrixRowStarts(A) == old_col_starts)
855 hypre_ParCSRMatrixRowStarts(A) = new_col_starts;
856 hypre_ParCSRMatrixOwnsRowStarts(A) = 1;
857 hypre_ParCSRMatrixOwnsColStarts(A) = 0;
861 hypre_ParCSRMatrixOwnsColStarts(A) = 1;
869 for (
int j = 0; j < size; j++)
871 diag(j) = A->diag->data[A->diag->i[j]];
872 MFEM_ASSERT(A->diag->j[A->diag->i[j]] == j,
873 "the first entry in each row must be the diagonal one");
877 static void MakeWrapper(
const hypre_CSRMatrix *mat,
SparseMatrix &wrapper)
879 HYPRE_Int nr = hypre_CSRMatrixNumRows(mat);
880 HYPRE_Int nc = hypre_CSRMatrixNumCols(mat);
883 hypre_CSRMatrixJ(mat),
884 hypre_CSRMatrixData(mat),
885 nr, nc,
false,
false,
false);
887 HYPRE_Int nnz = hypre_CSRMatrixNumNonzeros(mat);
888 SparseMatrix tmp(DuplicateAs<int>(hypre_CSRMatrixI(mat), nr+1),
889 DuplicateAs<int>(hypre_CSRMatrixJ(mat), nnz),
890 hypre_CSRMatrixData(mat),
891 nr, nc,
true,
false,
false);
898 MakeWrapper(A->diag, diag);
903 MakeWrapper(A->offd, offd);
904 cmap = A->col_map_offd;
908 bool interleaved_rows,
909 bool interleaved_cols)
const
914 hypre_ParCSRMatrix **hypre_blocks =
new hypre_ParCSRMatrix*[nr * nc];
915 internal::hypre_ParCSRMatrixSplit(A, nr, nc, hypre_blocks,
916 interleaved_rows, interleaved_cols);
918 for (
int i = 0; i < nr; i++)
920 for (
int j = 0; j < nc; j++)
926 delete [] hypre_blocks;
931 hypre_ParCSRMatrix * At;
932 hypre_ParCSRMatrixTranspose(A, &At, 1);
933 hypre_ParCSRMatrixSetNumNonzeros(At);
935 hypre_MatvecCommPkgCreate(At);
943 return hypre_ParCSRMatrixMatvec(a, A, x, b, y);
948 MFEM_ASSERT(x.
Size() ==
Width(),
"invalid x.Size() = " << x.
Size()
949 <<
", expected size = " <<
Width());
950 MFEM_ASSERT(y.
Size() ==
Height(),
"invalid y.Size() = " << y.
Size()
951 <<
", expected size = " <<
Height());
970 hypre_ParCSRMatrixMatvec(a, A, *X, b, *Y);
974 double b,
Vector &y)
const
976 MFEM_ASSERT(x.
Size() ==
Height(),
"invalid x.Size() = " << x.
Size()
977 <<
", expected size = " <<
Height());
978 MFEM_ASSERT(y.
Size() ==
Width(),
"invalid y.Size() = " << y.
Size()
979 <<
", expected size = " <<
Width());
1000 hypre_ParCSRMatrixMatvecT(a, A, *Y, b, *X);
1006 return hypre_ParCSRMatrixMatvec(a, A, (hypre_ParVector *) x, b,
1007 (hypre_ParVector *) y);
1013 return hypre_ParCSRMatrixMatvecT(a, A, x, b, y);
1017 HYPRE_Int* row_starts)
const
1019 const bool assumed_partition = HYPRE_AssumedPartitionCheck();
1020 const bool same_rows = (D.
Height() == hypre_CSRMatrixNumRows(A->diag));
1023 if (assumed_partition)
1029 MPI_Comm_size(
GetComm(), &part_size);
1033 HYPRE_Int global_num_rows;
1036 row_starts = hypre_ParCSRMatrixRowStarts(A);
1037 global_num_rows = hypre_ParCSRMatrixGlobalNumRows(A);
1041 MFEM_VERIFY(row_starts != NULL,
"the number of rows in D and A is not "
1042 "the same; row_starts must be given (not NULL)");
1047 assumed_partition ? row_starts[2] : row_starts[part_size-1];
1050 HYPRE_Int *col_starts = hypre_ParCSRMatrixColStarts(A);
1051 HYPRE_Int *col_map_offd;
1056 GetOffd(A_offd, col_map_offd);
1064 global_num_rows, hypre_ParCSRMatrixGlobalNumCols(A),
1065 DuplicateAs<HYPRE_Int>(row_starts, part_size,
false),
1066 DuplicateAs<HYPRE_Int>(col_starts, part_size,
false),
1068 DuplicateAs<HYPRE_Int>(col_map_offd, A_offd.
Width()));
1073 #ifndef HYPRE_BIGINT
1084 hypre_ParCSRMatrixSetRowStartsOwner(DA->A, 1);
1085 hypre_ParCSRMatrixSetColStartsOwner(DA->A, 1);
1087 DA->diagOwner = DA->offdOwner = 3;
1088 DA->colMapOwner = 1;
1095 if (hypre_CSRMatrixNumRows(A->diag) != hypre_CSRMatrixNumRows(A->offd))
1100 if (hypre_CSRMatrixNumRows(A->diag) != diag.
Size())
1102 mfem_error(
"Note the Vector diag is not of compatible dimensions with A\n");
1106 double *Adiag_data = hypre_CSRMatrixData(A->diag);
1107 HYPRE_Int *Adiag_i = hypre_CSRMatrixI(A->diag);
1110 double *Aoffd_data = hypre_CSRMatrixData(A->offd);
1111 HYPRE_Int *Aoffd_i = hypre_CSRMatrixI(A->offd);
1114 for (
int i(0); i < size; ++i)
1117 for (jj = Adiag_i[i]; jj < Adiag_i[i+1]; ++jj)
1119 Adiag_data[jj] *= val;
1121 for (jj = Aoffd_i[i]; jj < Aoffd_i[i+1]; ++jj)
1123 Aoffd_data[jj] *= val;
1130 if (hypre_CSRMatrixNumRows(A->diag) != hypre_CSRMatrixNumRows(A->offd))
1135 if (hypre_CSRMatrixNumRows(A->diag) != diag.
Size())
1137 mfem_error(
"Note the Vector diag is not of compatible dimensions with A\n");
1141 double *Adiag_data = hypre_CSRMatrixData(A->diag);
1142 HYPRE_Int *Adiag_i = hypre_CSRMatrixI(A->diag);
1145 double *Aoffd_data = hypre_CSRMatrixData(A->offd);
1146 HYPRE_Int *Aoffd_i = hypre_CSRMatrixI(A->offd);
1149 for (
int i(0); i < size; ++i)
1154 mfem_error(
"HypreParMatrix::InvDiagScale : Division by 0");
1158 for (jj = Adiag_i[i]; jj < Adiag_i[i+1]; ++jj)
1160 Adiag_data[jj] *= val;
1162 for (jj = Aoffd_i[i]; jj < Aoffd_i[i+1]; ++jj)
1164 Aoffd_data[jj] *= val;
1171 if (hypre_CSRMatrixNumRows(A->diag) != hypre_CSRMatrixNumRows(A->offd))
1176 HYPRE_Int size=hypre_CSRMatrixNumRows(A->diag);
1179 double *Adiag_data = hypre_CSRMatrixData(A->diag);
1180 HYPRE_Int *Adiag_i = hypre_CSRMatrixI(A->diag);
1181 for (jj = 0; jj < Adiag_i[size]; ++jj)
1183 Adiag_data[jj] *= s;
1186 double *Aoffd_data = hypre_CSRMatrixData(A->offd);
1187 HYPRE_Int *Aoffd_i = hypre_CSRMatrixI(A->offd);
1188 for (jj = 0; jj < Aoffd_i[size]; ++jj)
1190 Aoffd_data[jj] *= s;
1194 static void get_sorted_rows_cols(
const Array<int> &rows_cols,
1199 for (
int i = 0; i < rows_cols.
Size(); i++)
1201 hypre_sorted[i] = rows_cols[i];
1202 if (i && rows_cols[i-1] > rows_cols[i]) { sorted =
false; }
1204 if (!sorted) { hypre_sorted.
Sort(); }
1213 hypre_CSRMatrix * csr_A;
1214 hypre_CSRMatrix * csr_A_wo_z;
1215 hypre_ParCSRMatrix * parcsr_A_ptr;
1216 int * row_starts = NULL;
int * col_starts = NULL;
1217 int row_start = -1;
int row_end = -1;
1218 int col_start = -1;
int col_end = -1;
1220 comm = hypre_ParCSRMatrixComm(A);
1222 MPI_Comm_size(comm, &num_procs);
1224 ierr += hypre_ParCSRMatrixGetLocalRange(A,
1225 &row_start,&row_end,
1226 &col_start,&col_end );
1228 row_starts = hypre_ParCSRMatrixRowStarts(A);
1229 col_starts = hypre_ParCSRMatrixColStarts(A);
1231 parcsr_A_ptr = hypre_ParCSRMatrixCreate(comm,row_starts[num_procs],
1232 col_starts[num_procs],row_starts,
1235 csr_A = hypre_MergeDiagAndOffd(A);
1237 csr_A_wo_z = hypre_CSRMatrixDeleteZeros(csr_A,threshold);
1241 if (csr_A_wo_z == NULL)
1247 ierr += hypre_CSRMatrixDestroy(csr_A);
1250 ierr += GenerateDiagAndOffd(csr_A_wo_z,parcsr_A_ptr,
1253 ierr += hypre_CSRMatrixDestroy(csr_A_wo_z);
1255 ierr += hypre_ParCSRMatrixDestroy(A);
1265 get_sorted_rows_cols(rows_cols, rc_sorted);
1267 internal::hypre_ParCSRMatrixEliminateAXB(
1274 get_sorted_rows_cols(rows_cols, rc_sorted);
1276 hypre_ParCSRMatrix* Ae;
1277 internal::hypre_ParCSRMatrixEliminateAAe(
1285 hypre_ParCSRMatrixPrintIJ(A,offi,offj,fname);
1293 HYPRE_Int base_i, base_j;
1294 hypre_ParCSRMatrixReadIJ(comm, fname, &base_i, &base_j, &A);
1295 hypre_ParCSRMatrixSetNumNonzeros(A);
1297 hypre_MatvecCommPkgCreate(A);
1303 void HypreParMatrix::Destroy()
1305 if ( X != NULL ) {
delete X; }
1306 if ( Y != NULL ) {
delete Y; }
1308 if (A == NULL) {
return; }
1314 delete [] hypre_CSRMatrixI(A->diag);
1315 delete [] hypre_CSRMatrixJ(A->diag);
1317 hypre_CSRMatrixI(A->diag) = NULL;
1318 hypre_CSRMatrixJ(A->diag) = NULL;
1321 delete [] hypre_CSRMatrixData(A->diag);
1323 hypre_CSRMatrixData(A->diag) = NULL;
1329 delete [] hypre_CSRMatrixI(A->offd);
1330 delete [] hypre_CSRMatrixJ(A->offd);
1332 hypre_CSRMatrixI(A->offd) = NULL;
1333 hypre_CSRMatrixJ(A->offd) = NULL;
1336 delete [] hypre_CSRMatrixData(A->offd);
1338 hypre_CSRMatrixData(A->offd) = NULL;
1340 if (colMapOwner >= 0)
1342 if (colMapOwner & 1)
1344 delete [] hypre_ParCSRMatrixColMapOffd(A);
1346 hypre_ParCSRMatrixColMapOffd(A) = NULL;
1351 hypre_ParCSRMatrixDestroy(A);
1357 hypre_ParCSRMatrix * ab;
1358 ab = hypre_ParMatmul(*A,*B);
1360 hypre_MatvecCommPkgCreate(ab);
1367 HYPRE_Int P_owns_its_col_starts =
1368 hypre_ParCSRMatrixOwnsColStarts((hypre_ParCSRMatrix*)(*P));
1370 hypre_ParCSRMatrix * rap;
1371 hypre_BoomerAMGBuildCoarseOperator(*P,*A,*P,&rap);
1372 hypre_ParCSRMatrixSetNumNonzeros(rap);
1377 hypre_ParCSRMatrixSetRowStartsOwner(rap,0);
1378 hypre_ParCSRMatrixSetColStartsOwner(rap,0);
1380 if (P_owns_its_col_starts)
1382 hypre_ParCSRMatrixSetColStartsOwner(*P, 1);
1390 HYPRE_Int P_owns_its_col_starts =
1391 hypre_ParCSRMatrixOwnsColStarts((hypre_ParCSRMatrix*)(*P));
1392 HYPRE_Int Rt_owns_its_col_starts =
1393 hypre_ParCSRMatrixOwnsColStarts((hypre_ParCSRMatrix*)(*Rt));
1395 hypre_ParCSRMatrix * rap;
1396 hypre_BoomerAMGBuildCoarseOperator(*Rt,*A,*P,&rap);
1398 hypre_ParCSRMatrixSetNumNonzeros(rap);
1400 if (!P_owns_its_col_starts)
1404 hypre_ParCSRMatrixSetColStartsOwner(rap,0);
1406 if (!Rt_owns_its_col_starts)
1410 hypre_ParCSRMatrixSetRowStartsOwner(rap,0);
1420 Ae.
Mult(-1.0, X, 1.0, B);
1422 hypre_CSRMatrix *A_diag = hypre_ParCSRMatrixDiag((hypre_ParCSRMatrix *)A);
1423 double *data = hypre_CSRMatrixData(A_diag);
1424 HYPRE_Int *I = hypre_CSRMatrixI(A_diag);
1426 HYPRE_Int *J = hypre_CSRMatrixJ(A_diag);
1427 hypre_CSRMatrix *A_offd = hypre_ParCSRMatrixOffd((hypre_ParCSRMatrix *)A);
1428 HYPRE_Int *I_offd = hypre_CSRMatrixI(A_offd);
1429 double *data_offd = hypre_CSRMatrixData(A_offd);
1432 for (
int i = 0; i < ess_dof_list.
Size(); i++)
1434 int r = ess_dof_list[i];
1435 B(r) = data[I[r]] * X(r);
1442 MFEM_ABORT(
"the diagonal entry must be the first entry in the row!");
1444 for (
int j = I[r]+1; j < I[r+1]; j++)
1448 MFEM_ABORT(
"all off-diagonal entries must be zero!");
1451 for (
int j = I_offd[r]; j < I_offd[r+1]; j++)
1453 if (data_offd[j] != 0.0)
1455 MFEM_ABORT(
"all off-diagonal entries must be zero!");
1475 hypre_CSRMatrix *A_diag = hypre_ParCSRMatrixDiag(A);
1476 HYPRE_Int num_rows = hypre_CSRMatrixNumRows(A_diag);
1478 double *u_data = hypre_VectorData(hypre_ParVectorLocalVector(u));
1479 double *r_data = hypre_VectorData(hypre_ParVectorLocalVector(r));
1481 for (
int i = 0; i < N; i++)
1484 hypre_ParVectorCopy(f, r);
1485 hypre_ParCSRMatrixMatvec(-1.0, A, u, 1.0, r);
1488 (0 == (i % 2)) ? coef = lambda : coef = mu;
1490 for (HYPRE_Int j = 0; j < num_rows; j++)
1492 u_data[j] += coef*r_data[j] / max_eig;
1508 hypre_ParVector *x0,
1509 hypre_ParVector *x1,
1510 hypre_ParVector *x2,
1511 hypre_ParVector *x3)
1514 hypre_CSRMatrix *A_diag = hypre_ParCSRMatrixDiag(A);
1515 HYPRE_Int num_rows = hypre_CSRMatrixNumRows(A_diag);
1517 double *u_data = hypre_VectorData(hypre_ParVectorLocalVector(u));
1519 double *x0_data = hypre_VectorData(hypre_ParVectorLocalVector(x0));
1520 double *x1_data = hypre_VectorData(hypre_ParVectorLocalVector(x1));
1521 double *x2_data = hypre_VectorData(hypre_ParVectorLocalVector(x2));
1522 double *x3_data = hypre_VectorData(hypre_ParVectorLocalVector(x3));
1524 hypre_ParVectorCopy(u, x0);
1527 hypre_ParVectorCopy(f, x1);
1528 hypre_ParCSRMatrixMatvec(-1.0, A, x0, 1.0, x1);
1530 for (HYPRE_Int i = 0; i < num_rows; i++)
1532 x1_data[i] /= -max_eig;
1536 for (HYPRE_Int i = 0; i < num_rows; i++)
1538 x1_data[i] = x0_data[i] -x1_data[i];
1542 for (HYPRE_Int i = 0; i < num_rows; i++)
1544 x3_data[i] = fir_coeffs[0]*x0_data[i] +fir_coeffs[1]*x1_data[i];
1547 for (
int n = 2; n <= poly_order; n++)
1550 hypre_ParVectorCopy(f, x2);
1551 hypre_ParCSRMatrixMatvec(-1.0, A, x1, 1.0, x2);
1553 for (HYPRE_Int i = 0; i < num_rows; i++)
1555 x2_data[i] /= -max_eig;
1563 for (HYPRE_Int i = 0; i < num_rows; i++)
1565 x2_data[i] = (x1_data[i]-x0_data[i]) +(x1_data[i]-2*x2_data[i]);
1566 x3_data[i] += fir_coeffs[n]*x2_data[i];
1567 x0_data[i] = x1_data[i];
1568 x1_data[i] = x2_data[i];
1572 for (HYPRE_Int i = 0; i < num_rows; i++)
1574 u_data[i] = x3_data[i];
1593 B =
X =
V =
Z = NULL;
1599 int _relax_times,
double _relax_weight,
double _omega,
1600 int _poly_order,
double _poly_fraction)
1610 B =
X =
V =
Z = NULL;
1619 type =
static_cast<int>(_type);
1645 double a = -1, b, c;
1646 if (!strcmp(name,
"Rectangular")) { a = 1.0, b = 0.0, c = 0.0; }
1647 if (!strcmp(name,
"Hanning")) { a = 0.5, b = 0.5, c = 0.0; }
1648 if (!strcmp(name,
"Hamming")) { a = 0.54, b = 0.46, c = 0.0; }
1649 if (!strcmp(name,
"Blackman")) { a = 0.42, b = 0.50, c = 0.08; }
1652 mfem_error(
"HypreSmoother::SetWindowByName : name not recognized!");
1670 mfem_error(
"HypreSmoother::SetOperator : not HypreParMatrix!");
1676 if (
B) {
delete B; }
1677 if (
X) {
delete X; }
1678 if (
V) {
delete V; }
1679 if (
Z) {
delete Z; }
1698 A->
Mult(ones, diag);
1713 else if (
type == 1001 ||
type == 1002)
1744 double* window_coeffs =
new double[
poly_order+1];
1745 double* cheby_coeffs =
new double[
poly_order+1];
1753 window_coeffs[i] = a + b*cos(t) +c*cos(2*t);
1757 double theta_pb = acos(1.0 -0.5*k_pb);
1759 cheby_coeffs[0] = (theta_pb +sigma)/M_PI;
1762 double t = i*(theta_pb+sigma);
1763 cheby_coeffs[i] = 2.0*sin(t)/(i*M_PI);
1768 fir_coeffs[i] = window_coeffs[i]*cheby_coeffs[i];
1771 delete[] window_coeffs;
1772 delete[] cheby_coeffs;
1779 mfem_error(
"HypreSmoother::Mult (...) : HypreParMatrix A is missing");
1787 HYPRE_ParCSRDiagScale(NULL, *
A, b, x);
1811 else if (
type == 1002)
1826 hypre_ParCSRRelax(*
A, b,
type,
1831 hypre_ParCSRRelax(*
A, b,
type,
1842 mfem_error(
"HypreSmoother::Mult (...) : HypreParMatrix A is missing");
1848 A -> GetGlobalNumRows(),
1850 A -> GetRowStarts());
1852 A -> GetGlobalNumCols(),
1854 A -> GetColStarts());
1867 if (
B) {
delete B; }
1868 if (
X) {
delete X; }
1869 if (
V) {
delete V; }
1870 if (
Z) {
delete Z; }
1879 if (
X0) {
delete X0; }
1880 if (
X1) {
delete X1; }
1892 :
Solver(_A->Height(), _A->Width())
1903 mfem_error(
"HypreSolver::Mult (...) : HypreParMatrix A is missing");
1923 mfem_error(
"HypreSolver::Mult (...) : HypreParMatrix A is missing");
1929 A -> GetGlobalNumRows(),
1931 A -> GetRowStarts());
1933 A -> GetGlobalNumCols(),
1935 A -> GetColStarts());
1948 if (
B) {
delete B; }
1949 if (
X) {
delete X; }
1959 HYPRE_ParCSRMatrixGetComm(*
A, &comm);
1961 HYPRE_ParCSRPCGCreate(comm, &pcg_solver);
1966 HYPRE_PCGSetTol(pcg_solver, tol);
1971 HYPRE_PCGSetMaxIter(pcg_solver, max_iter);
1976 HYPRE_PCGSetLogging(pcg_solver, logging);
1981 HYPRE_ParCSRPCGSetPrintLevel(pcg_solver, print_lvl);
1986 HYPRE_ParCSRPCGSetPrecond(pcg_solver,
1994 HYPRE_PCGSetTwoNorm(pcg_solver, 1);
1995 if (res_frequency > 0)
1997 HYPRE_PCGSetRecomputeResidualP(pcg_solver, res_frequency);
2001 HYPRE_PCGSetResidualTol(pcg_solver, rtol);
2008 HYPRE_Int time_index = 0;
2009 HYPRE_Int num_iterations;
2010 double final_res_norm;
2012 HYPRE_Int print_level;
2014 HYPRE_PCGGetPrintLevel(pcg_solver, &print_level);
2015 HYPRE_ParCSRPCGSetPrintLevel(pcg_solver, print_level%3);
2017 HYPRE_ParCSRMatrixGetComm(*
A, &comm);
2021 if (print_level > 0 && print_level < 3)
2023 time_index = hypre_InitializeTiming(
"PCG Setup");
2024 hypre_BeginTiming(time_index);
2027 HYPRE_ParCSRPCGSetup(pcg_solver, *
A, b, x);
2030 if (print_level > 0 && print_level < 3)
2032 hypre_EndTiming(time_index);
2033 hypre_PrintTiming(
"Setup phase times", comm);
2034 hypre_FinalizeTiming(time_index);
2035 hypre_ClearTiming();
2039 if (print_level > 0 && print_level < 3)
2041 time_index = hypre_InitializeTiming(
"PCG Solve");
2042 hypre_BeginTiming(time_index);
2050 HYPRE_ParCSRPCGSolve(pcg_solver, *
A, b, x);
2052 if (print_level > 0)
2054 if (print_level < 3)
2056 hypre_EndTiming(time_index);
2057 hypre_PrintTiming(
"Solve phase times", comm);
2058 hypre_FinalizeTiming(time_index);
2059 hypre_ClearTiming();
2062 HYPRE_ParCSRPCGGetNumIterations(pcg_solver, &num_iterations);
2063 HYPRE_ParCSRPCGGetFinalRelativeResidualNorm(pcg_solver,
2066 MPI_Comm_rank(comm, &myid);
2070 cout <<
"PCG Iterations = " << num_iterations << endl
2071 <<
"Final PCG Relative Residual Norm = " << final_res_norm
2075 HYPRE_ParCSRPCGSetPrintLevel(pcg_solver, print_level);
2080 HYPRE_ParCSRPCGDestroy(pcg_solver);
2094 HYPRE_ParCSRMatrixGetComm(*
A, &comm);
2096 HYPRE_ParCSRGMRESCreate(comm, &gmres_solver);
2097 HYPRE_ParCSRGMRESSetKDim(gmres_solver, k_dim);
2098 HYPRE_ParCSRGMRESSetMaxIter(gmres_solver, max_iter);
2099 HYPRE_ParCSRGMRESSetTol(gmres_solver, tol);
2104 HYPRE_GMRESSetTol(gmres_solver, tol);
2109 HYPRE_GMRESSetMaxIter(gmres_solver, max_iter);
2114 HYPRE_GMRESSetKDim(gmres_solver, k_dim);
2119 HYPRE_GMRESSetLogging(gmres_solver, logging);
2124 HYPRE_GMRESSetPrintLevel(gmres_solver, print_lvl);
2129 HYPRE_ParCSRGMRESSetPrecond(gmres_solver,
2138 HYPRE_Int time_index = 0;
2139 HYPRE_Int num_iterations;
2140 double final_res_norm;
2142 HYPRE_Int print_level;
2144 HYPRE_GMRESGetPrintLevel(gmres_solver, &print_level);
2146 HYPRE_ParCSRMatrixGetComm(*
A, &comm);
2150 if (print_level > 0)
2152 time_index = hypre_InitializeTiming(
"GMRES Setup");
2153 hypre_BeginTiming(time_index);
2156 HYPRE_ParCSRGMRESSetup(gmres_solver, *
A, b, x);
2159 if (print_level > 0)
2161 hypre_EndTiming(time_index);
2162 hypre_PrintTiming(
"Setup phase times", comm);
2163 hypre_FinalizeTiming(time_index);
2164 hypre_ClearTiming();
2168 if (print_level > 0)
2170 time_index = hypre_InitializeTiming(
"GMRES Solve");
2171 hypre_BeginTiming(time_index);
2179 HYPRE_ParCSRGMRESSolve(gmres_solver, *
A, b, x);
2181 if (print_level > 0)
2183 hypre_EndTiming(time_index);
2184 hypre_PrintTiming(
"Solve phase times", comm);
2185 hypre_FinalizeTiming(time_index);
2186 hypre_ClearTiming();
2188 HYPRE_ParCSRGMRESGetNumIterations(gmres_solver, &num_iterations);
2189 HYPRE_ParCSRGMRESGetFinalRelativeResidualNorm(gmres_solver,
2192 MPI_Comm_rank(comm, &myid);
2196 cout <<
"GMRES Iterations = " << num_iterations << endl
2197 <<
"Final GMRES Relative Residual Norm = " << final_res_norm
2205 HYPRE_ParCSRGMRESDestroy(gmres_solver);
2213 int sai_max_levels = 1;
2214 double sai_threshold = 0.1;
2215 double sai_filter = 0.1;
2217 double sai_loadbal = 0.0;
2219 int sai_logging = 1;
2221 HYPRE_ParCSRMatrixGetComm(A, &comm);
2223 HYPRE_ParaSailsCreate(comm, &sai_precond);
2224 HYPRE_ParaSailsSetParams(sai_precond, sai_threshold, sai_max_levels);
2225 HYPRE_ParaSailsSetFilter(sai_precond, sai_filter);
2226 HYPRE_ParaSailsSetSym(sai_precond, sai_sym);
2227 HYPRE_ParaSailsSetLoadbal(sai_precond, sai_loadbal);
2228 HYPRE_ParaSailsSetReuse(sai_precond, sai_reuse);
2229 HYPRE_ParaSailsSetLogging(sai_precond, sai_logging);
2234 HYPRE_ParaSailsSetSym(sai_precond, sym);
2239 HYPRE_ParaSailsDestroy(sai_precond);
2245 HYPRE_BoomerAMGCreate(&amg_precond);
2246 SetDefaultOptions();
2251 HYPRE_BoomerAMGCreate(&amg_precond);
2252 SetDefaultOptions();
2255 void HypreBoomerAMG::SetDefaultOptions()
2258 int coarsen_type = 10;
2260 double theta = 0.25;
2263 int interp_type = 6;
2268 int relax_sweeps = 1;
2271 int print_level = 1;
2272 int max_levels = 25;
2274 HYPRE_BoomerAMGSetCoarsenType(amg_precond, coarsen_type);
2275 HYPRE_BoomerAMGSetAggNumLevels(amg_precond, agg_levels);
2276 HYPRE_BoomerAMGSetRelaxType(amg_precond, relax_type);
2277 HYPRE_BoomerAMGSetNumSweeps(amg_precond, relax_sweeps);
2278 HYPRE_BoomerAMGSetStrongThreshold(amg_precond, theta);
2279 HYPRE_BoomerAMGSetInterpType(amg_precond, interp_type);
2280 HYPRE_BoomerAMGSetPMaxElmts(amg_precond, Pmax);
2281 HYPRE_BoomerAMGSetPrintLevel(amg_precond, print_level);
2282 HYPRE_BoomerAMGSetMaxLevels(amg_precond, max_levels);
2285 HYPRE_BoomerAMGSetMaxIter(amg_precond, 1);
2286 HYPRE_BoomerAMGSetTol(amg_precond, 0.0);
2289 void HypreBoomerAMG::ResetAMGPrecond()
2291 HYPRE_Int coarsen_type;
2292 HYPRE_Int agg_levels;
2293 HYPRE_Int relax_type;
2294 HYPRE_Int relax_sweeps;
2296 HYPRE_Int interp_type;
2298 HYPRE_Int print_level;
2300 HYPRE_Int nrbms = rbms.
Size();
2302 HYPRE_Int nodal_diag;
2303 HYPRE_Int relax_coarse;
2304 HYPRE_Int interp_vec_variant;
2306 HYPRE_Int smooth_interp_vectors;
2307 HYPRE_Int interp_refine;
2309 hypre_ParAMGData *amg_data = (hypre_ParAMGData *)amg_precond;
2312 HYPRE_BoomerAMGGetCoarsenType(amg_precond, &coarsen_type);
2313 agg_levels = hypre_ParAMGDataAggNumLevels(amg_data);
2314 relax_type = hypre_ParAMGDataUserRelaxType(amg_data);
2315 relax_sweeps = hypre_ParAMGDataUserNumSweeps(amg_data);
2316 HYPRE_BoomerAMGGetStrongThreshold(amg_precond, &theta);
2317 hypre_BoomerAMGGetInterpType(amg_precond, &interp_type);
2318 HYPRE_BoomerAMGGetPMaxElmts(amg_precond, &Pmax);
2319 HYPRE_BoomerAMGGetPrintLevel(amg_precond, &print_level);
2320 HYPRE_BoomerAMGGetNumFunctions(amg_precond, &dim);
2323 nodal = hypre_ParAMGDataNodal(amg_data);
2324 nodal_diag = hypre_ParAMGDataNodalDiag(amg_data);
2325 HYPRE_BoomerAMGGetCycleRelaxType(amg_precond, &relax_coarse, 3);
2326 interp_vec_variant = hypre_ParAMGInterpVecVariant(amg_data);
2327 q_max = hypre_ParAMGInterpVecQMax(amg_data);
2328 smooth_interp_vectors = hypre_ParAMGSmoothInterpVectors(amg_data);
2329 interp_refine = hypre_ParAMGInterpRefine(amg_data);
2332 HYPRE_BoomerAMGDestroy(amg_precond);
2333 HYPRE_BoomerAMGCreate(&amg_precond);
2335 HYPRE_BoomerAMGSetCoarsenType(amg_precond, coarsen_type);
2336 HYPRE_BoomerAMGSetAggNumLevels(amg_precond, agg_levels);
2337 HYPRE_BoomerAMGSetRelaxType(amg_precond, relax_type);
2338 HYPRE_BoomerAMGSetNumSweeps(amg_precond, relax_sweeps);
2339 HYPRE_BoomerAMGSetMaxLevels(amg_precond, 25);
2340 HYPRE_BoomerAMGSetTol(amg_precond, 0.0);
2341 HYPRE_BoomerAMGSetMaxIter(amg_precond, 1);
2342 HYPRE_BoomerAMGSetStrongThreshold(amg_precond, theta);
2343 HYPRE_BoomerAMGSetInterpType(amg_precond, interp_type);
2344 HYPRE_BoomerAMGSetPMaxElmts(amg_precond, Pmax);
2345 HYPRE_BoomerAMGSetPrintLevel(amg_precond, print_level);
2346 HYPRE_BoomerAMGSetNumFunctions(amg_precond, dim);
2349 HYPRE_BoomerAMGSetNodal(amg_precond, nodal);
2350 HYPRE_BoomerAMGSetNodalDiag(amg_precond, nodal_diag);
2351 HYPRE_BoomerAMGSetCycleRelaxType(amg_precond, relax_coarse, 3);
2352 HYPRE_BoomerAMGSetInterpVecVariant(amg_precond, interp_vec_variant);
2353 HYPRE_BoomerAMGSetInterpVecQMax(amg_precond, q_max);
2354 HYPRE_BoomerAMGSetSmoothInterpVectors(amg_precond, smooth_interp_vectors);
2355 HYPRE_BoomerAMGSetInterpRefine(amg_precond, interp_refine);
2357 HYPRE_BoomerAMGSetInterpVectors(amg_precond, rbms.
Size(), rbms.
GetData());
2364 MFEM_VERIFY(new_A,
"new Operator must be a HypreParMatrix!");
2366 if (
A) { ResetAMGPrecond(); }
2380 HYPRE_BoomerAMGSetNumFunctions(amg_precond, dim);
2383 HYPRE_BoomerAMGSetAggNumLevels(amg_precond, 0);
2384 HYPRE_BoomerAMGSetStrongThreshold(amg_precond, 0.5);
2390 y = 0.0; y(0) = x(1); y(1) = -x(0);
2392 static void func_ryz(
const Vector &x, Vector &y)
2394 y = 0.0; y(1) = x(2); y(2) = -x(1);
2396 static void func_rzx(
const Vector &x, Vector &y)
2398 y = 0.0; y(2) = x(0); y(0) = -x(2);
2401 void HypreBoomerAMG::RecomputeRBMs()
2404 Array<HypreParVector*> gf_rbms;
2407 for (
int i = 0; i < rbms.
Size(); i++)
2409 HYPRE_ParVectorDestroy(rbms[i]);
2416 VectorFunctionCoefficient coeff_rxy(2, func_rxy);
2418 ParGridFunction rbms_rxy(fespace);
2419 rbms_rxy.ProjectCoefficient(coeff_rxy);
2422 gf_rbms.SetSize(nrbms);
2423 gf_rbms[0] = rbms_rxy.ParallelAverage();
2429 VectorFunctionCoefficient coeff_rxy(3, func_rxy);
2430 VectorFunctionCoefficient coeff_ryz(3, func_ryz);
2431 VectorFunctionCoefficient coeff_rzx(3, func_rzx);
2433 ParGridFunction rbms_rxy(fespace);
2434 ParGridFunction rbms_ryz(fespace);
2435 ParGridFunction rbms_rzx(fespace);
2436 rbms_rxy.ProjectCoefficient(coeff_rxy);
2437 rbms_ryz.ProjectCoefficient(coeff_ryz);
2438 rbms_rzx.ProjectCoefficient(coeff_rzx);
2441 gf_rbms.SetSize(nrbms);
2442 gf_rbms[0] = rbms_rxy.ParallelAverage();
2443 gf_rbms[1] = rbms_ryz.ParallelAverage();
2444 gf_rbms[2] = rbms_rzx.ParallelAverage();
2453 for (
int i = 0; i < nrbms; i++)
2455 rbms[i] = gf_rbms[i]->StealParVector();
2463 this->fespace = fespace;
2473 int relax_coarse = 8;
2476 int interp_vec_variant = 2;
2478 int smooth_interp_vectors = 1;
2482 int interp_refine = 1;
2484 HYPRE_BoomerAMGSetNodal(amg_precond, nodal);
2485 HYPRE_BoomerAMGSetNodalDiag(amg_precond, nodal_diag);
2486 HYPRE_BoomerAMGSetCycleRelaxType(amg_precond, relax_coarse, 3);
2487 HYPRE_BoomerAMGSetInterpVecVariant(amg_precond, interp_vec_variant);
2488 HYPRE_BoomerAMGSetInterpVecQMax(amg_precond, q_max);
2489 HYPRE_BoomerAMGSetSmoothInterpVectors(amg_precond, smooth_interp_vectors);
2490 HYPRE_BoomerAMGSetInterpRefine(amg_precond, interp_refine);
2493 HYPRE_BoomerAMGSetInterpVectors(amg_precond, rbms.
Size(), rbms.
GetData());
2498 for (
int i = 0; i < rbms.
Size(); i++)
2500 HYPRE_ParVectorDestroy(rbms[i]);
2503 HYPRE_BoomerAMGDestroy(amg_precond);
2510 int cycle_type = 13;
2513 double rlx_weight = 1.0;
2514 double rlx_omega = 1.0;
2515 int amg_coarsen_type = 10;
2516 int amg_agg_levels = 1;
2517 int amg_rlx_type = 8;
2518 double theta = 0.25;
2519 int amg_interp_type = 6;
2526 bool trace_space, rt_trace_space;
2530 trace_space = trace_space || rt_trace_space;
2533 if (edge_fespace->
GetNE() > 0)
2538 if (dim == 2) { p++; }
2553 HYPRE_AMSCreate(&ams);
2555 HYPRE_AMSSetDimension(ams, sdim);
2556 HYPRE_AMSSetTol(ams, 0.0);
2557 HYPRE_AMSSetMaxIter(ams, 1);
2558 HYPRE_AMSSetCycleType(ams, cycle_type);
2559 HYPRE_AMSSetPrintLevel(ams, 1);
2581 for (
int i = 0; i < pmesh->
GetNV(); i++)
2583 coord = pmesh -> GetVertex(i);
2584 x_coord(i) = coord[0];
2585 y_coord(i) = coord[1];
2586 if (sdim == 3) { z_coord(i) = coord[2]; }
2593 HYPRE_AMSSetCoordinateVectors(ams, *x, *y, NULL);
2598 HYPRE_AMSSetCoordinateVectors(ams, *x, *y, *z);
2622 HYPRE_AMSSetDiscreteGradient(ams, *G);
2626 Pi = Pix = Piy = Piz = NULL;
2645 if (cycle_type < 10)
2653 Pix = Pi_blocks(0,0);
2654 Piy = Pi_blocks(0,1);
2655 if (sdim == 3) { Piz = Pi_blocks(0,2); }
2660 HYPRE_ParCSRMatrix HY_Pi = (Pi) ? (HYPRE_ParCSRMatrix) *Pi : NULL;
2661 HYPRE_ParCSRMatrix HY_Pix = (Pix) ? (HYPRE_ParCSRMatrix) *Pix : NULL;
2662 HYPRE_ParCSRMatrix HY_Piy = (Piy) ? (HYPRE_ParCSRMatrix) *Piy : NULL;
2663 HYPRE_ParCSRMatrix HY_Piz = (Piz) ? (HYPRE_ParCSRMatrix) *Piz : NULL;
2664 HYPRE_AMSSetInterpolations(ams, HY_Pi, HY_Pix, HY_Piy, HY_Piz);
2666 delete vert_fespace_d;
2669 delete vert_fespace;
2674 delete edge_fespace;
2679 HYPRE_AMSSetSmoothingOptions(ams, rlx_type, rlx_sweeps, rlx_weight, rlx_omega);
2680 HYPRE_AMSSetAlphaAMGOptions(ams, amg_coarsen_type, amg_agg_levels, amg_rlx_type,
2681 theta, amg_interp_type, amg_Pmax);
2682 HYPRE_AMSSetBetaAMGOptions(ams, amg_coarsen_type, amg_agg_levels, amg_rlx_type,
2683 theta, amg_interp_type, amg_Pmax);
2688 HYPRE_AMSDestroy(ams);
2703 HYPRE_AMSSetPrintLevel(ams, print_lvl);
2709 int cycle_type = 11;
2712 double rlx_weight = 1.0;
2713 double rlx_omega = 1.0;
2714 int amg_coarsen_type = 10;
2715 int amg_agg_levels = 1;
2716 int amg_rlx_type = 8;
2717 double theta = 0.25;
2718 int amg_interp_type = 6;
2720 int ams_cycle_type = 14;
2726 if (face_fespace->
GetNE() > 0)
2738 HYPRE_ADSCreate(&ads);
2740 HYPRE_ADSSetTol(ads, 0.0);
2741 HYPRE_ADSSetMaxIter(ads, 1);
2742 HYPRE_ADSSetCycleType(ads, cycle_type);
2743 HYPRE_ADSSetPrintLevel(ads, 1);
2771 for (
int i = 0; i < pmesh->
GetNV(); i++)
2773 coord = pmesh -> GetVertex(i);
2774 x_coord(i) = coord[0];
2775 y_coord(i) = coord[1];
2776 z_coord(i) = coord[2];
2781 HYPRE_ADSSetCoordinateVectors(ads, *x, *y, *z);
2805 HYPRE_ADSSetDiscreteCurl(ads, *C);
2824 HYPRE_ADSSetDiscreteGradient(ads, *G);
2828 RT_Pi = RT_Pix = RT_Piy = RT_Piz = NULL;
2829 ND_Pi = ND_Pix = ND_Piy = ND_Piz = NULL;
2848 if (ams_cycle_type < 10)
2858 ND_Pix = ND_Pi_blocks(0,0);
2859 ND_Piy = ND_Pi_blocks(0,1);
2860 ND_Piz = ND_Pi_blocks(0,2);
2878 if (cycle_type < 10)
2887 RT_Pix = RT_Pi_blocks(0,0);
2888 RT_Piy = RT_Pi_blocks(0,1);
2889 RT_Piz = RT_Pi_blocks(0,2);
2894 HYPRE_ParCSRMatrix HY_RT_Pi, HY_RT_Pix, HY_RT_Piy, HY_RT_Piz;
2895 HY_RT_Pi = (RT_Pi) ? (HYPRE_ParCSRMatrix) *RT_Pi : NULL;
2896 HY_RT_Pix = (RT_Pix) ? (HYPRE_ParCSRMatrix) *RT_Pix : NULL;
2897 HY_RT_Piy = (RT_Piy) ? (HYPRE_ParCSRMatrix) *RT_Piy : NULL;
2898 HY_RT_Piz = (RT_Piz) ? (HYPRE_ParCSRMatrix) *RT_Piz : NULL;
2899 HYPRE_ParCSRMatrix HY_ND_Pi, HY_ND_Pix, HY_ND_Piy, HY_ND_Piz;
2900 HY_ND_Pi = (ND_Pi) ? (HYPRE_ParCSRMatrix) *ND_Pi : NULL;
2901 HY_ND_Pix = (ND_Pix) ? (HYPRE_ParCSRMatrix) *ND_Pix : NULL;
2902 HY_ND_Piy = (ND_Piy) ? (HYPRE_ParCSRMatrix) *ND_Piy : NULL;
2903 HY_ND_Piz = (ND_Piz) ? (HYPRE_ParCSRMatrix) *ND_Piz : NULL;
2904 HYPRE_ADSSetInterpolations(ads,
2905 HY_RT_Pi, HY_RT_Pix, HY_RT_Piy, HY_RT_Piz,
2906 HY_ND_Pi, HY_ND_Pix, HY_ND_Piy, HY_ND_Piz);
2908 delete vert_fespace_d;
2912 delete vert_fespace;
2914 delete edge_fespace;
2917 HYPRE_ADSSetSmoothingOptions(ads, rlx_type, rlx_sweeps, rlx_weight, rlx_omega);
2918 HYPRE_ADSSetAMGOptions(ads, amg_coarsen_type, amg_agg_levels, amg_rlx_type,
2919 theta, amg_interp_type, amg_Pmax);
2920 HYPRE_ADSSetAMSOptions(ads, ams_cycle_type, amg_coarsen_type, amg_agg_levels,
2921 amg_rlx_type, theta, amg_interp_type, amg_Pmax);
2926 HYPRE_ADSDestroy(ads);
2948 HYPRE_ADSSetPrintLevel(ads, print_lvl);
2951 HypreLOBPCG::HypreMultiVector::HypreMultiVector(
int n,
HypreParVector & v,
2952 mv_InterfaceInterpreter & interpreter)
2956 mv_ptr = mv_MultiVectorCreateFromSampleVector(&interpreter, nv,
2957 (HYPRE_ParVector)v);
2959 HYPRE_ParVector* vecs = NULL;
2961 mv_TempMultiVector* tmp =
2962 (mv_TempMultiVector*)mv_MultiVectorGetData(mv_ptr);
2963 vecs = (HYPRE_ParVector*)(tmp -> vector);
2966 hpv =
new HypreParVector*[nv];
2967 for (
int i=0; i<nv; i++)
2969 hpv[i] =
new HypreParVector(vecs[i]);
2973 HypreLOBPCG::HypreMultiVector::~HypreMultiVector()
2977 for (
int i=0; i<nv; i++)
2984 mv_MultiVectorDestroy(mv_ptr);
2988 HypreLOBPCG::HypreMultiVector::Randomize(HYPRE_Int seed)
2990 mv_MultiVectorSetRandom(mv_ptr, seed);
2994 HypreLOBPCG::HypreMultiVector::GetVector(
unsigned int i)
2996 MFEM_ASSERT((
int)i < nv,
"index out of range");
3002 HypreLOBPCG::HypreMultiVector::StealVectors()
3004 HypreParVector ** hpv_ret = hpv;
3008 mv_TempMultiVector * mv_tmp =
3009 (mv_TempMultiVector*)mv_MultiVectorGetData(mv_ptr);
3011 mv_tmp->ownsVectors = 0;
3013 for (
int i=0; i<nv; i++)
3015 hpv_ret[i]->SetOwnership(1);
3032 MPI_Comm_size(comm,&numProcs);
3033 MPI_Comm_rank(comm,&myid);
3035 HYPRE_ParCSRSetupInterpreter(&interpreter);
3036 HYPRE_ParCSRSetupMatvec(&matvec_fn);
3037 HYPRE_LOBPCGCreate(&interpreter, &matvec_fn, &lobpcg_solver);
3046 HYPRE_LOBPCGDestroy(lobpcg_solver);
3052 HYPRE_LOBPCGSetTol(lobpcg_solver, tol);
3058 HYPRE_LOBPCGSetMaxIter(lobpcg_solver, max_iter);
3066 HYPRE_LOBPCGSetPrintLevel(lobpcg_solver, logging);
3073 HYPRE_LOBPCGSetPrecondUsageMode(lobpcg_solver, pcg_mode);
3079 HYPRE_LOBPCGSetPrecond(lobpcg_solver,
3080 (HYPRE_PtrToSolverFcn)this->PrecondSolve,
3081 (HYPRE_PtrToSolverFcn)this->PrecondSetup,
3082 (HYPRE_Solver)&precond);
3088 int locSize = A.
Width();
3090 if (HYPRE_AssumedPartitionCheck())
3092 part =
new HYPRE_Int[2];
3094 MPI_Scan(&locSize, &part[1], 1, HYPRE_MPI_INT, MPI_SUM, comm);
3096 part[0] = part[1] - locSize;
3098 MPI_Allreduce(&locSize, &glbSize, 1, HYPRE_MPI_INT, MPI_SUM, comm);
3102 part =
new HYPRE_Int[numProcs+1];
3104 MPI_Allgather(&locSize, 1, MPI_INT,
3105 &part[1], 1, HYPRE_MPI_INT, comm);
3108 for (
int i=0; i<numProcs; i++)
3110 part[i+1] += part[i];
3113 glbSize = part[numProcs];
3124 matvec_fn.MatvecCreate = this->OperatorMatvecCreate;
3125 matvec_fn.Matvec = this->OperatorMatvec;
3126 matvec_fn.MatvecDestroy = this->OperatorMatvecDestroy;
3128 HYPRE_LOBPCGSetup(lobpcg_solver,(HYPRE_Matrix)&A,NULL,NULL);
3134 matvec_fn.MatvecCreate = this->OperatorMatvecCreate;
3135 matvec_fn.Matvec = this->OperatorMatvec;
3136 matvec_fn.MatvecDestroy = this->OperatorMatvecDestroy;
3138 HYPRE_LOBPCGSetupB(lobpcg_solver,(HYPRE_Matrix)&M,NULL);
3147 for (
int i=0; i<nev; i++)
3149 eigs[i] = eigenvalues[i];
3156 return multi_vec->GetVector(i);
3163 if ( multi_vec == NULL )
3165 MFEM_ASSERT(x != NULL,
"In HypreLOBPCG::Solve()");
3167 multi_vec =
new HypreMultiVector(nev, *x, interpreter);
3168 multi_vec->Randomize(seed);
3178 HYPRE_LOBPCGSolve(lobpcg_solver, NULL, *multi_vec, eigenvalues);
3182 HypreLOBPCG::OperatorMatvecCreate(
void *A,
3189 return ( matvec_data );
3193 HypreLOBPCG::OperatorMatvec(
void *matvec_data,
3194 HYPRE_Complex
alpha,
3200 MFEM_VERIFY(alpha == 1.0 && beta == 0.0,
"values not supported");
3202 Operator *Aop = (Operator*)A;
3204 int width = Aop->Width();
3206 hypre_ParVector * xPar = (hypre_ParVector *)x;
3207 hypre_ParVector * yPar = (hypre_ParVector *)y;
3209 Vector xVec(xPar->local_vector->data, width);
3210 Vector yVec(yPar->local_vector->data, width);
3212 Aop->Mult( xVec, yVec );
3218 HypreLOBPCG::OperatorMatvecDestroy(
void *matvec_data )
3224 HypreLOBPCG::PrecondSolve(
void *solver,
3229 Solver *PC = (Solver*)solver;
3230 Operator *OP = (Operator*)A;
3232 int width = OP->Width();
3234 hypre_ParVector * bPar = (hypre_ParVector *)b;
3235 hypre_ParVector * xPar = (hypre_ParVector *)x;
3237 Vector bVec(bPar->local_vector->data, width);
3238 Vector xVec(xPar->local_vector->data, width);
3240 PC->Mult( bVec, xVec );
3246 HypreLOBPCG::PrecondSetup(
void *solver,
3264 MPI_Comm_size(comm,&numProcs);
3265 MPI_Comm_rank(comm,&myid);
3267 HYPRE_AMECreate(&ame_solver);
3268 HYPRE_AMESetPrintLevel(ame_solver, 0);
3275 hypre_TFree(multi_vec);
3280 for (
int i=0; i<nev; i++)
3282 delete eigenvectors[i];
3285 delete [] eigenvectors;
3289 hypre_TFree(eigenvalues);
3292 HYPRE_AMEDestroy(ame_solver);
3300 HYPRE_AMESetBlockSize(ame_solver, nev);
3306 HYPRE_AMESetTol(ame_solver, tol);
3312 HYPRE_AMESetMaxIter(ame_solver, max_iter);
3320 HYPRE_AMESetPrintLevel(ame_solver, logging);
3327 ams_precond = &precond;
3335 HYPRE_Solver ams_precond_ptr = (HYPRE_Solver)*ams_precond;
3337 ams_precond->
SetupFcn()(*ams_precond,A,NULL,NULL);
3339 HYPRE_AMESetAMSSolver(ame_solver, ams_precond_ptr);
3342 HYPRE_AMESetup(ame_solver);
3348 HYPRE_ParCSRMatrix parcsr_M = M;
3349 HYPRE_AMESetMassMatrix(ame_solver,(HYPRE_ParCSRMatrix)parcsr_M);
3355 HYPRE_AMESolve(ame_solver);
3362 eigs.
SetSize(nev); eigs = -1.0;
3364 if ( eigenvalues == NULL )
3367 HYPRE_AMEGetEigenvalues(ame_solver,&eigenvalues);
3371 for (
int i=0; i<nev; i++)
3373 eigs[i] = eigenvalues[i];
3378 HypreAME::createDummyVectors()
3380 if ( multi_vec == NULL )
3382 HYPRE_AMEGetEigenvectors(ame_solver,&multi_vec);
3385 eigenvectors =
new HypreParVector*[nev];
3386 for (
int i=0; i<nev; i++)
3388 eigenvectors[i] =
new HypreParVector(multi_vec[i]);
3397 if ( eigenvectors == NULL )
3399 this->createDummyVectors();
3402 return *eigenvectors[i];
3408 if ( eigenvectors == NULL )
3410 this->createDummyVectors();
3415 eigenvectors = NULL;
virtual ~HypreBoomerAMG()
void SetPreconditioner(HypreSolver &precond)
Set the hypre solver to be used as a preconditioner.
int Size() const
Logical size of the array.
void EliminateRowsCols(const Array< int > &rows_cols, const HypreParVector &X, HypreParVector &B)
virtual int NumNonZeroElems() const
Returns the number of the nonzero elements in the matrix.
Vector()
Default constructor for Vector. Sets size = 0 and data = NULL.
MPI_Comm GetComm() const
MPI communicator.
Vector * GlobalVector() const
Returns the global vector in each processor.
HypreParVector * X0
FIR Filter Temporary Vectors.
HypreParVector & GetEigenvector(unsigned int i)
Extract a single eigenvector.
double min_eig_est
Minimal eigenvalue estimate for polynomial smoothing.
void SetPrintLevel(int print_lvl)
void Print(const char *fname, HYPRE_Int offi=0, HYPRE_Int offj=0)
Prints the locally owned rows in parallel.
int setup_called
Was hypre's Setup function called already?
void MakeRef(const HypreParMatrix &master)
Make this HypreParMatrix a reference to 'master'.
HypreParMatrix * LeftDiagMult(const SparseMatrix &D, HYPRE_Int *row_starts=NULL) const
HYPRE_Int MultTranspose(HypreParVector &x, HypreParVector &y, double alpha=1.0, double beta=0.0)
Computes y = alpha * A^t * x + beta * y.
void SetSize(int s)
Resize the vector if the new size is different.
HypreParVector * B
Right-hand side and solution vector.
double window_params[3]
Parameters for windowing function of FIR filter.
void Mult(const Table &A, const Table &B, Table &C)
C = A * B (as boolean matrices)
int Width() const
Get the width (size of input) of the Operator. Synonym with NumCols.
HypreADS(HypreParMatrix &A, ParFiniteElementSpace *face_fespace)
void SetWindowByName(const char *window_name)
Convenience function for setting canonical windowing parameters.
int GetFaceOrder(int i) const
Returns the order of the i'th face finite element.
T * GetData()
Returns the data.
void GetEigenvalues(Array< double > &eigenvalues)
Collect the converged eigenvalues.
HYPRE_Int GlobalTrueVSize()
void SetPreconditioner(HypreSolver &precond)
int Size() const
Returns the size of the vector.
void SetPreconditioner(Solver &precond)
void SetMassMatrix(Operator &M)
Abstract parallel finite element space.
double Normlinf() const
Returns the l_infinity norm of the vector.
void SetPrintLevel(int logging)
bool iterative_mode
If true, use the second argument of Mult as an initial guess.
void SetOperator(HypreParMatrix &A)
void LoseData()
Lose the ownership of the graph (I, J) and data (A) arrays.
HYPRE_Int Mult(HypreParVector &x, HypreParVector &y, double alpha=1.0, double beta=0.0)
Computes y = alpha * A * x + beta * y.
void SetPrintLevel(int print_lvl)
void ScaleRows(const Vector &s)
Scale the local row i by s(i).
int Size_of_connections() const
double poly_fraction
Fraction of spectrum to smooth for polynomial relaxation.
void SetPrintLevel(int print_lvl)
int poly_scale
Apply the polynomial smoother to A or D^{-1/2} A D^{-1/2}.
void AddDomainInterpolator(DiscreteInterpolator *di)
HYPRE_Int * GetTrueDofOffsets()
void SetWindowParameters(double a, double b, double c)
Set parameters for windowing function for FIR smoother.
virtual void Mult(const HypreParVector &b, HypreParVector &x) const
Solve Ax=b with hypre's GMRES.
void Print(const char *fname) const
Prints the locally owned rows in parallel.
void EliminateBC(HypreParMatrix &A, HypreParMatrix &Ae, const Array< int > &ess_dof_list, const Vector &X, Vector &B)
void SetSystemsOptions(int dim)
HypreLOBPCG(MPI_Comm comm)
virtual void Mult(const HypreParVector &b, HypreParVector &x) const
Relax the linear system Ax=b.
int GetNE() const
Returns number of elements in the mesh.
HYPRE_Int GetGlobalNumRows() const
void SetSymmetry(int sym)
virtual ~HypreParaSails()
HYPRE_Int GetGlobalNumCols() const
void SetMaxIter(int max_iter)
double * l1_norms
l1 norms of the rows of A
int Height() const
Get the height (size of output) of the Operator. Synonym with NumRows.
void SetLogging(int logging)
Mesh * GetMesh() const
Returns the mesh.
HyprePCG(HypreParMatrix &_A)
virtual HYPRE_PtrToParSolverFcn SolveFcn() const =0
hypre's internal Solve function
void SetResidualConvergenceOptions(int res_frequency=-1, double rtol=0.0)
void GetBlocks(Array2D< HypreParMatrix * > &blocks, bool interleaved_rows=false, bool interleaved_cols=false) const
virtual void Mult(const HypreParVector &b, HypreParVector &x) const
Solve the linear system Ax=b.
HypreParMatrix * RAP(HypreParMatrix *A, HypreParMatrix *P)
Returns the matrix P^t * A * P.
HypreParVector ** StealEigenvectors()
Transfer ownership of the converged eigenvectors.
HYPRE_Int Randomize(HYPRE_Int seed)
Set random values.
double relax_weight
Damping coefficient (usually <= 1)
void SetElasticityOptions(ParFiniteElementSpace *fespace)
void GetDiag(Vector &diag) const
Get the local diagonal of the matrix.
void Sort()
Sorts the array. This requires operator< to be defined for T.
int Size() const
Returns the number of TYPE I elements.
HypreParMatrix * A
The linear system matrix.
double * GetData() const
Return element data.
HypreAMS(HypreParMatrix &A, ParFiniteElementSpace *edge_fespace)
int * GetI() const
Return the array I.
virtual void SetOperator(const Operator &op)
virtual HYPRE_PtrToParSolverFcn SetupFcn() const =0
hypre's internal Setup function
void SetLogging(int logging)
HypreParMatrix()
An empty matrix to be used as a reference to an existing matrix.
void SetMaxIter(int max_iter)
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)
HypreParMatrix * ParallelAssemble() const
Returns the matrix "assembled" on the true dofs.
void SetSOROptions(double relax_weight, double omega)
Set SOR-related parameters.
int SpaceDimension() const
Wrapper for hypre's parallel vector class.
double Norml1() const
Returns the l_1 norm of the vector.
int ParCSRRelax_Taubin(hypre_ParCSRMatrix *A, hypre_ParVector *f, double lambda, double mu, int N, double max_eig, hypre_ParVector *u, hypre_ParVector *r)
void GetEigenvalues(Array< double > &eigenvalues)
Collect the converged eigenvalues.
void SetPrintLevel(int print_lvl)
void mfem_error(const char *msg)
void SetSize(int nsize)
Change logical size of the array, keep existing entries.
void Solve()
Solve the eigenproblem.
void AddTraceFaceInterpolator(DiscreteInterpolator *di)
int GetOrder(int i) const
Returns the order of the i'th finite element.
void SetNumModes(int num_eigs)
int GetNumRows() const
Returns the number of rows in the diagonal block of the ParCSRMatrix.
double max_eig_est
Maximal eigenvalue estimate for polynomial smoothing.
HypreParMatrix * Transpose()
Returns the transpose of *this.
void SetMaxIter(int max_iter)
void SetFIRCoefficients(double max_eig)
Compute window and Chebyshev coefficients for given polynomial order.
double * fir_coeffs
Combined coefficients for windowing and Chebyshev polynomials.
void Threshold(double threshold=0.0)
Remove values smaller in absolute value than some threshold.
double InnerProduct(HypreParVector *x, HypreParVector *y)
HypreParVector(MPI_Comm comm, HYPRE_Int glob_size, HYPRE_Int *col)
virtual void SetOperator(const Operator &op)
Set/update the solver for the given operator.
void SetData(double *_data)
void InvScaleRows(const Vector &s)
Scale the local row i by 1./s(i)
void SetTaubinOptions(double lambda, double mu, int iter)
Set parameters for Taubin's lambda-mu method.
HypreParVector * B
Right-hand side and solution vectors.
void SetMassMatrix(HypreParMatrix &M)
int GetNumCols() const
Returns the number of columns in the diagonal block of the ParCSRMatrix.
void SetOperator(Operator &A)
void SetPreconditioner(HypreSolver &precond)
Set the hypre solver to be used as a preconditioner.
void SetPolyOptions(int poly_order, double poly_fraction)
Set parameters for polynomial smoothing.
int relax_times
Number of relaxation sweeps.
void operator*=(double s)
Scale all entries by s: A_scaled = s*A.
void SetPrintLevel(int logging)
double ParNormlp(const Vector &vec, double p, MPI_Comm comm)
Compute the l_p norm of the Vector which is split without overlap across the given communicator...
void SetOwnership(int own)
Sets ownership of the internal hypre_ParVector.
void SetDataOwner(bool owna)
Set the data ownership flag (A array).
Abstract class for hypre's solvers and preconditioners.
HypreParVector & operator=(double d)
Set constant values.
Arbitrary order H(curl)-conforming Nedelec finite elements.
const FiniteElementCollection * FEColl() const
HypreGMRES(HypreParMatrix &_A)
void GetParBlocks(Array2D< HypreParMatrix * > &blocks) const
hypre_ParCSRMatrix * StealData()
Changes the ownership of the the matrix.
void SetPrecondUsageMode(int pcg_mode)
HypreParMatrix * A
The linear system matrix.
Arbitrary order H1-conforming (continuous) finite elements.
void SetMaxIter(int max_iter)
void GetOffd(SparseMatrix &offd, HYPRE_Int *&cmap) const
Get the local off-diagonal block. NOTE: 'offd' will not own any data.
HYPRE_Int * GetColStarts() const
HypreParVector & GetEigenvector(unsigned int i)
Extract a single eigenvector.
double omega
SOR parameter (usually in (0,2))
HYPRE_Int * GetRowStarts() const
Class for parallel grid function.
Wrapper for hypre's ParCSR matrix class.
void Swap(SparseMatrix &other)
~HypreParVector()
Calls hypre's destroy function.
virtual void Assemble(int skip_zeros=1)
int * GetJ() const
Return the array J.
virtual void Mult(const HypreParVector &b, HypreParVector &x) const
Solve Ax=b with hypre's PCG.
Class for parallel meshes.
void Solve()
Solve the eigenproblem.
void Read(MPI_Comm comm, const char *fname)
Reads the matrix from a file.
HypreParVector * V
Temporary vectors.
void SetType(HypreSmoother::Type type, int relax_times=1)
Set the relaxation type and number of sweeps.
HypreParaSails(HypreParMatrix &A)
int poly_order
Order of the smoothing polynomial.
double lambda
Taubin's lambda-mu method parameters.
HypreParMatrix * ParMult(HypreParMatrix *A, HypreParMatrix *B)
Returns the matrix A * B.
void ParallelProject(Vector &tv) const
Returns the vector restricted to the true dofs.