1228 diagOwner = HypreCsrToMem(A->diag, host_mt,
true, mem_diag);
1229 offdOwner = HypreCsrToMem(A->offd, host_mt,
true, mem_offd);
1232 hypre_CSRMatrixSetRownnz(A->diag);
1233 hypre_CSRMatrixSetRownnz(A->offd);
1252 if (HYPRE_AssumedPartitionCheck())
1255 my_col_start = cols[0];
1256 my_col_end = cols[1];
1261 MPI_Comm_rank(comm, &myid);
1262 MPI_Comm_size(comm, &part_size);
1264 my_col_start = cols[myid];
1265 my_col_end = cols[myid+1];
1272 row_starts = col_starts = mfem_hypre_TAlloc_host(
HYPRE_BigInt, part_size);
1273 for (
int i = 0; i < part_size; i++)
1275 row_starts[i] = rows[i];
1280 row_starts = mfem_hypre_TAlloc_host(
HYPRE_BigInt, part_size);
1281 col_starts = mfem_hypre_TAlloc_host(
HYPRE_BigInt, part_size);
1282 for (
int i = 0; i < part_size; i++)
1284 row_starts[i] = rows[i];
1285 col_starts[i] = cols[i];
1291 HYPRE_Int diag_nnz = 0, offd_nnz = 0, offd_num_cols = 0;
1292 map<HYPRE_BigInt, HYPRE_Int> offd_map;
1293 for (HYPRE_Int j = 0, loc_nnz = I[nrows]; j < loc_nnz; j++)
1296 if (my_col_start <= glob_col && glob_col < my_col_end)
1302 offd_map.insert(pair<const HYPRE_BigInt, HYPRE_Int>(glob_col, -1));
1307 for (
auto it = offd_map.begin(); it != offd_map.end(); ++it)
1309 it->second = offd_num_cols++;
1313 A = hypre_ParCSRMatrixCreate(comm, glob_nrows, glob_ncols,
1314 row_starts, col_starts, offd_num_cols,
1315 diag_nnz, offd_nnz);
1316 hypre_ParCSRMatrixInitialize(A);
1322 HYPRE_Int *diag_i, *diag_j, *offd_i, *offd_j;
1324 real_t *diag_data, *offd_data;
1325 diag_i = A->diag->i;
1326 diag_j = A->diag->j;
1327 diag_data = A->diag->data;
1328 offd_i = A->offd->i;
1329 offd_j = A->offd->j;
1330 offd_data = A->offd->data;
1331 offd_col_map = A->col_map_offd;
1333 diag_nnz = offd_nnz = 0;
1334 for (HYPRE_Int i = 0, j = 0; i < nrows; i++)
1336 diag_i[i] = diag_nnz;
1337 offd_i[i] = offd_nnz;
1338 for (HYPRE_Int j_end = I[i+1]; j < j_end; j++)
1341 if (my_col_start <= glob_col && glob_col < my_col_end)
1343 diag_j[diag_nnz] = glob_col - my_col_start;
1344 diag_data[diag_nnz] = data[j];
1349 offd_j[offd_nnz] = offd_map[glob_col];
1350 offd_data[offd_nnz] = data[j];
1355 diag_i[nrows] = diag_nnz;
1356 offd_i[nrows] = offd_nnz;
1357 for (
auto it = offd_map.begin(); it != offd_map.end(); ++it)
1359 offd_col_map[it->second] = it->first;
1362 hypre_ParCSRMatrixSetNumNonzeros(A);
1364 if (row_starts == col_starts)
1366 hypre_CSRMatrixReorder(hypre_ParCSRMatrixDiag(A));
1368#if MFEM_HYPRE_VERSION > 22200
1369 mfem_hypre_TFree_host(row_starts);
1372 mfem_hypre_TFree_host(col_starts);
1375 hypre_MatvecCommPkgCreate(A);
1385 hypre_ParCSRMatrix *Ph =
static_cast<hypre_ParCSRMatrix *
>(P);
1390 A = hypre_ParCSRMatrixCompleteClone(Ph);
1392 hypre_ParCSRMatrixCopy(Ph, A, 1);
1400 hypre_ParCSRMatrixSetNumNonzeros(A);
1402 hypre_MatvecCommPkgCreate(A);
1431 MFEM_ASSERT(diagOwner < 0 && offdOwner < 0 && colMapOwner == -1,
"");
1432 MFEM_ASSERT(diagOwner == offdOwner,
"");
1433 MFEM_ASSERT(ParCSROwner,
"");
1434 hypre_ParCSRMatrix *R = A;
1435#ifdef HYPRE_USING_GPU
1442 ParCSROwner =
false;
1470 colMapOwner = colmap;
1475#if MFEM_HYPRE_VERSION <= 22200
1476 if (!A || hypre_ParCSRMatrixOwnsRowStarts(A) ||
1477 (hypre_ParCSRMatrixRowStarts(A) == hypre_ParCSRMatrixColStarts(A) &&
1478 hypre_ParCSRMatrixOwnsColStarts(A)))
1483 int row_starts_size;
1484 if (HYPRE_AssumedPartitionCheck())
1486 row_starts_size = 2;
1490 MPI_Comm_size(hypre_ParCSRMatrixComm(A), &row_starts_size);
1494 HYPRE_BigInt *old_row_starts = hypre_ParCSRMatrixRowStarts(A);
1497 for (
int i = 0; i < row_starts_size; i++)
1499 new_row_starts[i] = old_row_starts[i];
1502 hypre_ParCSRMatrixRowStarts(A) = new_row_starts;
1503 hypre_ParCSRMatrixOwnsRowStarts(A) = 1;
1505 if (hypre_ParCSRMatrixColStarts(A) == old_row_starts)
1507 hypre_ParCSRMatrixColStarts(A) = new_row_starts;
1508 hypre_ParCSRMatrixOwnsColStarts(A) = 0;
1515#if MFEM_HYPRE_VERSION <= 22200
1516 if (!A || hypre_ParCSRMatrixOwnsColStarts(A) ||
1517 (hypre_ParCSRMatrixRowStarts(A) == hypre_ParCSRMatrixColStarts(A) &&
1518 hypre_ParCSRMatrixOwnsRowStarts(A)))
1523 int col_starts_size;
1524 if (HYPRE_AssumedPartitionCheck())
1526 col_starts_size = 2;
1530 MPI_Comm_size(hypre_ParCSRMatrixComm(A), &col_starts_size);
1534 HYPRE_BigInt *old_col_starts = hypre_ParCSRMatrixColStarts(A);
1537 for (
int i = 0; i < col_starts_size; i++)
1539 new_col_starts[i] = old_col_starts[i];
1542 hypre_ParCSRMatrixColStarts(A) = new_col_starts;
1544 if (hypre_ParCSRMatrixRowStarts(A) == old_col_starts)
1546 hypre_ParCSRMatrixRowStarts(A) = new_col_starts;
1547 hypre_ParCSRMatrixOwnsRowStarts(A) = 1;
1548 hypre_ParCSRMatrixOwnsColStarts(A) = 0;
1552 hypre_ParCSRMatrixOwnsColStarts(A) = 1;
1559 const int size =
Height();
1564 MemoryClass hypre_mc = (hypre_ml == HYPRE_MEMORY_HOST) ?
1567#if MFEM_HYPRE_VERSION >= 21800
1568 MFEM_VERIFY(A->diag->memory_location == hypre_ml,
1569 "unexpected HypreParMatrix memory location!");
1571 const HYPRE_Int *A_diag_i = A->diag->i;
1572 const real_t *A_diag_d = A->diag->data;
1574 const HYPRE_Int *A_diag_j = A->diag->j;
1578 diag_hd[i] = A_diag_d[A_diag_i[i]];
1580 A_diag_j[A_diag_i[i]] == i,
1581 "The first entry in each row must be the diagonal one!");
1585static void MakeSparseMatrixWrapper(
int nrows,
int ncols,
1586 HYPRE_Int *I, HYPRE_Int *J,
real_t *data,
1590 SparseMatrix tmp(I, J, data, nrows, ncols,
false,
false,
false);
1593 for (
int i = 0; i <= nrows; i++)
1595 mI[i] = internal::to_int(I[i]);
1597 const int nnz = mI[nrows];
1598 int *mJ = Memory<int>(nnz);
1599 for (
int j = 0; j < nnz; j++)
1601 mJ[j] = internal::to_int(J[j]);
1603 SparseMatrix tmp(mI, mJ, data, nrows, ncols,
true,
false,
false);
1608static void MakeWrapper(
const hypre_CSRMatrix *mat,
1609 const MemoryIJData &mem,
1610 SparseMatrix &wrapper)
1612 const int nrows = internal::to_int(hypre_CSRMatrixNumRows(mat));
1613 const int ncols = internal::to_int(hypre_CSRMatrixNumCols(mat));
1614 const int nnz = internal::to_int(mat->num_nonzeros);
1618 MakeSparseMatrixWrapper(nrows, ncols,
1619 const_cast<HYPRE_Int*
>(I),
1620 const_cast<HYPRE_Int*
>(J),
1621 const_cast<real_t*
>(data),
1627 MakeWrapper(A->diag, mem_diag, diag);
1632 MakeWrapper(A->offd, mem_offd, offd);
1633 cmap = A->col_map_offd;
1639 hypre_CSRMatrix *hypre_merged = hypre_MergeDiagAndOffd(A);
1643#if MFEM_HYPRE_VERSION >= 21600
1644 hypre_CSRMatrixBigJtoJ(hypre_merged);
1646 MakeSparseMatrixWrapper(
1647 internal::to_int(hypre_merged->num_rows),
1648 internal::to_int(hypre_merged->num_cols),
1655 merged = merged_tmp;
1657 hypre_CSRMatrixDestroy(hypre_merged);
1661 bool interleaved_rows,
1662 bool interleaved_cols)
const
1667 hypre_ParCSRMatrix **hypre_blocks =
new hypre_ParCSRMatrix*[nr * nc];
1669 internal::hypre_ParCSRMatrixSplit(A, nr, nc, hypre_blocks,
1670 interleaved_rows, interleaved_cols);
1673 for (
int i = 0; i < nr; i++)
1675 for (
int j = 0; j < nc; j++)
1681 delete [] hypre_blocks;
1686 hypre_ParCSRMatrix * At;
1687 hypre_ParCSRMatrixTranspose(A, &At, 1);
1688 hypre_ParCSRMatrixSetNumNonzeros(At);
1690 if (!hypre_ParCSRMatrixCommPkg(At)) { hypre_MatvecCommPkgCreate(At); }
1696 hypre_CSRMatrixReorder(hypre_ParCSRMatrixDiag(At));
1702#if MFEM_HYPRE_VERSION >= 21800
1712 hypre_MatvecCommPkgCreate(A);
1715 hypre_ParCSRMatrix *submat;
1718 int local_num_vars = hypre_CSRMatrixNumRows(hypre_ParCSRMatrixDiag(A));
1721#ifdef hypre_IntArrayData
1723 hypre_IntArray *CF_marker;
1725 CF_marker = hypre_IntArrayCreate(local_num_vars);
1726 hypre_IntArrayInitialize_v2(CF_marker, HYPRE_MEMORY_HOST);
1727 hypre_IntArraySetConstantValues(CF_marker, 1);
1732 for (
int j=0; j<indices.
Size(); j++)
1734 if (indices[j] > local_num_vars)
1736 MFEM_WARNING(
"WARNING : " << indices[j] <<
" > " << local_num_vars);
1738#ifdef hypre_IntArrayData
1739 hypre_IntArrayData(CF_marker)[indices[j]] = -1;
1741 CF_marker[indices[j]] = -1;
1746#if (MFEM_HYPRE_VERSION > 22300) || (MFEM_HYPRE_VERSION == 22300 && HYPRE_DEVELOP_NUMBER >=8)
1749 hypre_BoomerAMGCoarseParms(MPI_COMM_WORLD, local_num_vars, 1, NULL,
1750 CF_marker, NULL, cpts_global);
1753 hypre_BoomerAMGCoarseParms(MPI_COMM_WORLD, local_num_vars, 1, NULL,
1754 CF_marker, NULL, &cpts_global);
1758#ifdef hypre_IntArrayData
1759 hypre_ParCSRMatrixExtractSubmatrixFC(A, hypre_IntArrayData(CF_marker),
1760 cpts_global,
"FF", &submat,
1763 hypre_ParCSRMatrixExtractSubmatrixFC(A, CF_marker, cpts_global,
1764 "FF", &submat, threshold);
1767#if (MFEM_HYPRE_VERSION <= 22300) && !(MFEM_HYPRE_VERSION == 22300 && HYPRE_DEVELOP_NUMBER >=8)
1768 mfem_hypre_TFree(cpts_global);
1770#ifdef hypre_IntArrayData
1771 hypre_IntArrayDestroy(CF_marker);
1782#if (MFEM_HYPRE_VERSION == 22500 && HYPRE_DEVELOP_NUMBER >= 1) || \
1783 (MFEM_HYPRE_VERSION > 22500)
1784#ifdef HYPRE_USING_GPU
1787 hypre_ParCSRMatrixLocalTranspose(A);
1795#if (MFEM_HYPRE_VERSION == 22500 && HYPRE_DEVELOP_NUMBER >= 1) || \
1796 (MFEM_HYPRE_VERSION > 22500)
1797#ifdef HYPRE_USING_GPU
1802 hypre_CSRMatrixDestroy(A->diagT);
1807 hypre_CSRMatrixDestroy(A->offdT);
1820 return hypre_ParCSRMatrixMatvec(
a, A, x,
b, y);
1825 MFEM_ASSERT(x.
Size() ==
Width(),
"invalid x.Size() = " << x.
Size()
1826 <<
", expected size = " <<
Width());
1827 MFEM_ASSERT(y.
Size() ==
Height(),
"invalid y.Size() = " << y.
Size()
1828 <<
", expected size = " <<
Height());
1875 hypre_ParCSRMatrixMatvec(
a, A, *X,
b, *Y);
1877 if (!yshallow) { y = *Y; }
1883 MFEM_ASSERT(x.
Size() ==
Height(),
"invalid x.Size() = " << x.
Size()
1884 <<
", expected size = " <<
Height());
1885 MFEM_ASSERT(y.
Size() ==
Width(),
"invalid y.Size() = " << y.
Size()
1886 <<
", expected size = " <<
Width());
1939 hypre_ParCSRMatrixMatvecT(
a, A, *Y,
b, *X);
1941 if (!yshallow) { y = *X; }
1947 return hypre_ParCSRMatrixMatvec(
a, A, (hypre_ParVector *) x,
b,
1948 (hypre_ParVector *) y);
1957 return hypre_ParCSRMatrixMatvecT(
a, A, x,
b, y);
1963 MFEM_ASSERT(x.
Size() ==
Width(),
"invalid x.Size() = " << x.
Size()
1964 <<
", expected size = " <<
Width());
1965 MFEM_ASSERT(y.
Size() ==
Height(),
"invalid y.Size() = " << y.
Size()
1966 <<
", expected size = " <<
Height());
1972 internal::hypre_ParCSRMatrixAbsMatvec(A,
a,
const_cast<real_t*
>(x_data),
1980 MFEM_ASSERT(x.
Size() ==
Height(),
"invalid x.Size() = " << x.
Size()
1981 <<
", expected size = " <<
Height());
1982 MFEM_ASSERT(y.
Size() ==
Width(),
"invalid y.Size() = " << y.
Size()
1983 <<
", expected size = " <<
Width());
1989 internal::hypre_ParCSRMatrixAbsMatvecT(A,
a,
const_cast<real_t*
>(x_data),
1997 const bool assumed_partition = HYPRE_AssumedPartitionCheck();
1998 const bool row_starts_given = (row_starts != NULL);
1999 if (!row_starts_given)
2001 row_starts = hypre_ParCSRMatrixRowStarts(A);
2002 MFEM_VERIFY(D.
Height() == hypre_CSRMatrixNumRows(A->diag),
2003 "the matrix D is NOT compatible with the row starts of"
2004 " this HypreParMatrix, row_starts must be given.");
2009 if (assumed_partition)
2015 MPI_Comm_rank(
GetComm(), &offset);
2017 int local_num_rows = row_starts[offset+1]-row_starts[offset];
2018 MFEM_VERIFY(local_num_rows == D.
Height(),
"the number of rows in D is "
2019 " not compatible with the given row_starts");
2026 if (assumed_partition)
2029 if (row_starts_given)
2031 global_num_rows = row_starts[2];
2038 global_num_rows = hypre_ParCSRMatrixGlobalNumRows(A);
2043 MPI_Comm_size(
GetComm(), &part_size);
2044 global_num_rows = row_starts[part_size];
2048 HYPRE_BigInt *col_starts = hypre_ParCSRMatrixColStarts(A);
2054 GetOffd(A_offd, col_map_offd);
2064 DuplicateAs<HYPRE_BigInt>(row_starts, part_size,
false);
2066 (row_starts == col_starts ? new_row_starts :
2067 DuplicateAs<HYPRE_BigInt>(col_starts, part_size,
false));
2069 DuplicateAs<HYPRE_BigInt>(col_map_offd, A_offd.
Width());
2073 const bool own_diag_offd =
true;
2078 global_num_rows, hypre_ParCSRMatrixGlobalNumCols(A),
2079 new_row_starts, new_col_starts,
2080 DA_diag, DA_offd, new_col_map_offd,
2083#if MFEM_HYPRE_VERSION <= 22200
2085 hypre_ParCSRMatrixSetRowStartsOwner(DA->A, 1);
2086 hypre_ParCSRMatrixSetColStartsOwner(DA->A, 1);
2088 mfem_hypre_TFree_host(new_row_starts);
2089 mfem_hypre_TFree_host(new_col_starts);
2091 DA->colMapOwner = 1;
2098 if (hypre_CSRMatrixNumRows(A->diag) != hypre_CSRMatrixNumRows(A->offd))
2103 if (hypre_CSRMatrixNumRows(A->diag) != diag.
Size())
2105 mfem_error(
"Note the Vector diag is not of compatible dimensions with A\n");
2112 real_t *Adiag_data = hypre_CSRMatrixData(A->diag);
2113 HYPRE_Int *Adiag_i = hypre_CSRMatrixI(A->diag);
2115 real_t *Aoffd_data = hypre_CSRMatrixData(A->offd);
2116 HYPRE_Int *Aoffd_i = hypre_CSRMatrixI(A->offd);
2119 for (
int i(0); i < size; ++i)
2122 for (jj = Adiag_i[i]; jj < Adiag_i[i+1]; ++jj)
2124 Adiag_data[jj] *= val;
2126 for (jj = Aoffd_i[i]; jj < Aoffd_i[i+1]; ++jj)
2128 Aoffd_data[jj] *= val;
2137 if (hypre_CSRMatrixNumRows(A->diag) != hypre_CSRMatrixNumRows(A->offd))
2142 if (hypre_CSRMatrixNumRows(A->diag) != diag.
Size())
2144 mfem_error(
"Note the Vector diag is not of compatible dimensions with A\n");
2151 real_t *Adiag_data = hypre_CSRMatrixData(A->diag);
2152 HYPRE_Int *Adiag_i = hypre_CSRMatrixI(A->diag);
2155 real_t *Aoffd_data = hypre_CSRMatrixData(A->offd);
2156 HYPRE_Int *Aoffd_i = hypre_CSRMatrixI(A->offd);
2159 for (
int i(0); i < size; ++i)
2164 mfem_error(
"HypreParMatrix::InvDiagScale : Division by 0");
2168 for (jj = Adiag_i[i]; jj < Adiag_i[i+1]; ++jj)
2170 Adiag_data[jj] *= val;
2172 for (jj = Aoffd_i[i]; jj < Aoffd_i[i+1]; ++jj)
2174 Aoffd_data[jj] *= val;
2183 if (hypre_CSRMatrixNumRows(A->diag) != hypre_CSRMatrixNumRows(A->offd))
2190 HYPRE_Int size=hypre_CSRMatrixNumRows(A->diag);
2193 real_t *Adiag_data = hypre_CSRMatrixData(A->diag);
2194 HYPRE_Int *Adiag_i = hypre_CSRMatrixI(A->diag);
2195 for (jj = 0; jj < Adiag_i[size]; ++jj)
2197 Adiag_data[jj] *=
s;
2200 real_t *Aoffd_data = hypre_CSRMatrixData(A->offd);
2201 HYPRE_Int *Aoffd_i = hypre_CSRMatrixI(A->offd);
2202 for (jj = 0; jj < Aoffd_i[size]; ++jj)
2204 Aoffd_data[jj] *=
s;
2210static void get_sorted_rows_cols(
const Array<int> &rows_cols,
2216 for (
int i = 0; i < rows_cols.
Size(); i++)
2218 hypre_sorted[i] = rows_cols[i];
2219 if (i && rows_cols[i-1] > rows_cols[i]) { sorted =
false; }
2221 if (!sorted) { hypre_sorted.
Sort(); }
2229 hypre_CSRMatrix * csr_A;
2230 hypre_CSRMatrix * csr_A_wo_z;
2231 hypre_ParCSRMatrix * parcsr_A_ptr;
2236 comm = hypre_ParCSRMatrixComm(A);
2238 ierr += hypre_ParCSRMatrixGetLocalRange(A,
2239 &row_start,&row_end,
2240 &col_start,&col_end );
2242 row_starts = hypre_ParCSRMatrixRowStarts(A);
2243 col_starts = hypre_ParCSRMatrixColStarts(A);
2245#if MFEM_HYPRE_VERSION <= 22200
2246 bool old_owns_row = hypre_ParCSRMatrixOwnsRowStarts(A);
2247 bool old_owns_col = hypre_ParCSRMatrixOwnsColStarts(A);
2249 HYPRE_BigInt global_num_rows = hypre_ParCSRMatrixGlobalNumRows(A);
2250 HYPRE_BigInt global_num_cols = hypre_ParCSRMatrixGlobalNumCols(A);
2251 parcsr_A_ptr = hypre_ParCSRMatrixCreate(comm, global_num_rows,
2253 row_starts, col_starts,
2255#if MFEM_HYPRE_VERSION <= 22200
2256 hypre_ParCSRMatrixOwnsRowStarts(parcsr_A_ptr) = old_owns_row;
2257 hypre_ParCSRMatrixOwnsColStarts(parcsr_A_ptr) = old_owns_col;
2258 hypre_ParCSRMatrixOwnsRowStarts(A) = 0;
2259 hypre_ParCSRMatrixOwnsColStarts(A) = 0;
2262 csr_A = hypre_MergeDiagAndOffd(A);
2268 csr_A_wo_z = hypre_CSRMatrixDeleteZeros(csr_A,threshold);
2272 if (csr_A_wo_z == NULL)
2278 ierr += hypre_CSRMatrixDestroy(csr_A);
2284 ierr += GenerateDiagAndOffd(csr_A_wo_z,parcsr_A_ptr,
2287 ierr += hypre_CSRMatrixDestroy(csr_A_wo_z);
2289 MFEM_VERIFY(ierr == 0,
"");
2293 hypre_ParCSRMatrixSetNumNonzeros(A);
2295#if MFEM_HYPRE_VERSION <= 22200
2296 if (row_starts == col_starts)
2298 if ((row_starts[0] == col_starts[0]) &&
2299 (row_starts[1] == col_starts[1]))
2302 hypre_CSRMatrixReorder(hypre_ParCSRMatrixDiag(A));
2304 if (!hypre_ParCSRMatrixCommPkg(A)) { hypre_MatvecCommPkgCreate(A); }
2311 HYPRE_Int old_err = hypre_error_flag;
2312 hypre_error_flag = 0;
2314#if MFEM_HYPRE_VERSION < 21400
2319 HYPRE_Int *diag_I = A->diag->i, *offd_I = A->offd->i;
2320 real_t *diag_d = A->diag->data, *offd_d = A->offd->data;
2321 HYPRE_Int local_num_rows = A->diag->num_rows;
2322 real_t max_l2_row_norm = 0.0;
2324 for (HYPRE_Int r = 0; r < local_num_rows; r++)
2329 l2_row_norm = std::hypot(l2_row_norm, row.
Norml2());
2330 max_l2_row_norm = std::max(max_l2_row_norm, l2_row_norm);
2332 real_t loc_max_l2_row_norm = max_l2_row_norm;
2333 MPI_Allreduce(&loc_max_l2_row_norm, &max_l2_row_norm, 1,
2336 threshold = tol * max_l2_row_norm;
2341#elif MFEM_HYPRE_VERSION < 21800
2343 HYPRE_Int err_flag = hypre_ParCSRMatrixDropSmallEntries(A, tol);
2344 MFEM_VERIFY(!err_flag,
"error encountered: error code = " << err_flag);
2348 HYPRE_Int err_flag = hypre_ParCSRMatrixDropSmallEntries(A, tol, 2);
2349 MFEM_VERIFY(!err_flag,
"error encountered: error code = " << err_flag);
2353 hypre_error_flag = old_err;
2361 get_sorted_rows_cols(rows_cols, rc_sorted);
2363 internal::hypre_ParCSRMatrixEliminateAXB(
2370 get_sorted_rows_cols(rows_cols, rc_sorted);
2372 hypre_ParCSRMatrix* Ae;
2374 internal::hypre_ParCSRMatrixEliminateAAe(
2384 get_sorted_rows_cols(cols, rc_sorted);
2386 hypre_ParCSRMatrix* Ae;
2388 internal::hypre_ParCSRMatrixEliminateAAe(
2389 A, &Ae, rc_sorted.
Size(), rc_sorted.
GetData(), 1);
2397 if (rows.
Size() > 0)
2400 get_sorted_rows_cols(rows, r_sorted);
2402 internal::hypre_ParCSRMatrixEliminateRows(A, r_sorted.
Size(),
2413 Ae.
Mult(-1.0, x, 1.0,
b);
2417 if (ess_dof_list.
Size() == 0) {
return; }
2420 hypre_CSRMatrix *A_diag = hypre_ParCSRMatrixDiag(A);
2421 real_t *data = hypre_CSRMatrixData(A_diag);
2422 HYPRE_Int *I = hypre_CSRMatrixI(A_diag);
2424 HYPRE_Int *J = hypre_CSRMatrixJ(A_diag);
2425 hypre_CSRMatrix *A_offd = hypre_ParCSRMatrixOffd(A);
2426 HYPRE_Int *I_offd = hypre_CSRMatrixI(A_offd);
2427 real_t *data_offd = hypre_CSRMatrixData(A_offd);
2434 for (
int i = 0; i < ess_dof_list.
Size(); i++)
2436 int r = ess_dof_list[i];
2437 b(r) = data[I[r]] * x(r);
2439 MFEM_ASSERT(I[r] < I[r+1],
"empty row found!");
2445 MFEM_ABORT(
"the diagonal entry must be the first entry in the row!");
2447 for (
int j = I[r]+1; j < I[r+1]; j++)
2451 MFEM_ABORT(
"all off-diagonal entries must be zero!");
2454 for (
int j = I_offd[r]; j < I_offd[r+1]; j++)
2456 if (data_offd[j] != 0.0)
2458 MFEM_ABORT(
"all off-diagonal entries must be zero!");
2469 hypre_ParCSRMatrix *A_hypre = *
this;
2472 hypre_CSRMatrix *diag = hypre_ParCSRMatrixDiag(A_hypre);
2473 hypre_CSRMatrix *offd = hypre_ParCSRMatrixOffd(A_hypre);
2475 HYPRE_Int diag_nrows = hypre_CSRMatrixNumRows(diag);
2476 HYPRE_Int offd_ncols = hypre_CSRMatrixNumCols(offd);
2478 const int n_ess_dofs = ess_dofs.
Size();
2479 const auto ess_dofs_d = ess_dofs.
GetMemory().Read(
2484 hypre_ParCSRCommHandle *comm_handle;
2485 HYPRE_Int *int_buf_data, *eliminate_row, *eliminate_col;
2487 eliminate_row = mfem_hypre_CTAlloc(HYPRE_Int, diag_nrows);
2488 eliminate_col = mfem_hypre_CTAlloc(HYPRE_Int, offd_ncols);
2491 hypre_ParCSRCommPkg *comm_pkg = hypre_ParCSRMatrixCommPkg(A_hypre);
2494 hypre_MatvecCommPkgCreate(A_hypre);
2495 comm_pkg = hypre_ParCSRMatrixCommPkg(A_hypre);
2501 eliminate_row[i] = 0;
2505 eliminate_row[ess_dofs_d[i]] = 1;
2511 HYPRE_Int num_sends = hypre_ParCSRCommPkgNumSends(comm_pkg);
2512 HYPRE_Int int_buf_sz = hypre_ParCSRCommPkgSendMapStart(comm_pkg, num_sends);
2513 int_buf_data = mfem_hypre_CTAlloc(HYPRE_Int, int_buf_sz);
2515 HYPRE_Int *send_map_elmts;
2516#if defined(HYPRE_USING_GPU)
2519 hypre_ParCSRCommPkgCopySendMapElmtsToDevice(comm_pkg);
2520 send_map_elmts = hypre_ParCSRCommPkgDeviceSendMapElmts(comm_pkg);
2525 send_map_elmts = hypre_ParCSRCommPkgSendMapElmts(comm_pkg);
2529 int k = send_map_elmts[i];
2530 int_buf_data[i] = eliminate_row[k];
2533#if defined(HYPRE_USING_GPU)
2537 comm_handle = hypre_ParCSRCommHandleCreate_v2(
2538 11, comm_pkg, HYPRE_MEMORY_DEVICE, int_buf_data,
2539 HYPRE_MEMORY_DEVICE, eliminate_col);
2544 comm_handle = hypre_ParCSRCommHandleCreate(
2545 11, comm_pkg, int_buf_data, eliminate_col );
2551 const auto I = diag->i;
2552 const auto J = diag->j;
2553 auto data = diag->data;
2557 const int idof = ess_dofs_d[i];
2558 for (
auto j=I[idof]; j<I[idof+1]; ++j)
2560 const auto jdof = J[j];
2576 for (
auto k=I[jdof]; k<I[jdof+1]; ++k)
2591 const auto I = offd->i;
2592 auto data = offd->data;
2595 const int idof = ess_dofs_d[i];
2596 for (
auto j=I[idof]; j<I[idof+1]; ++j)
2604 hypre_ParCSRCommHandleDestroy(comm_handle);
2605 mfem_hypre_TFree(int_buf_data);
2606 mfem_hypre_TFree(eliminate_row);
2610 const int nrows_offd = hypre_CSRMatrixNumRows(offd);
2611 const auto I = offd->i;
2612 const auto J = offd->j;
2613 auto data = offd->data;
2616 for (
auto j=I[i]; j<I[i+1]; ++j)
2618 data[j] *= 1 - eliminate_col[J[j]];
2623 mfem_hypre_TFree(eliminate_col);
2627 HYPRE_Int offj)
const
2630 hypre_ParCSRMatrixPrintIJ(A,offi,offj,fname);
2634void HypreParMatrix::Read(MPI_Comm comm,
const char *fname)
2639 HYPRE_Int base_i, base_j;
2640 hypre_ParCSRMatrixReadIJ(comm, fname, &base_i, &base_j, &A);
2641 hypre_ParCSRMatrixSetNumNonzeros(A);
2643 if (!hypre_ParCSRMatrixCommPkg(A)) { hypre_MatvecCommPkgCreate(A); }
2654 HYPRE_IJMatrix A_ij;
2655 HYPRE_IJMatrixRead(fname, comm, 5555, &A_ij);
2657 HYPRE_ParCSRMatrix A_parcsr;
2658 HYPRE_IJMatrixGetObject(A_ij, (
void**) &A_parcsr);
2660 A = (hypre_ParCSRMatrix*)A_parcsr;
2662 hypre_ParCSRMatrixSetNumNonzeros(A);
2664 if (!hypre_ParCSRMatrixCommPkg(A)) { hypre_MatvecCommPkgCreate(A); }
2672 hypre_ParCSRCommPkg *comm_pkg = A->comm_pkg;
2673 MPI_Comm comm = A->comm;
2675 const int tag = 46801;
2677 MPI_Comm_rank(comm, &myid);
2678 MPI_Comm_size(comm, &nproc);
2682 MPI_Recv(&c, 1, MPI_CHAR, myid-1, tag, comm, MPI_STATUS_IGNORE);
2686 os <<
"\nHypreParMatrix: hypre_ParCSRCommPkg:\n";
2688 os <<
"Rank " << myid <<
":\n"
2689 " number of sends = " << comm_pkg->num_sends <<
2690 " (" <<
sizeof(
real_t)*comm_pkg->send_map_starts[comm_pkg->num_sends] <<
2692 " number of recvs = " << comm_pkg->num_recvs <<
2693 " (" <<
sizeof(
real_t)*comm_pkg->recv_vec_starts[comm_pkg->num_recvs] <<
2695 if (myid != nproc-1)
2698 MPI_Send(&c, 1, MPI_CHAR, myid+1, tag, comm);
2711 os <<
"global number of rows : " << A->global_num_rows <<
'\n'
2712 <<
"global number of columns : " << A->global_num_cols <<
'\n'
2713 <<
"first row index : " << A->first_row_index <<
'\n'
2714 <<
" last row index : " << A->last_row_index <<
'\n'
2715 <<
"first col diag : " << A->first_col_diag <<
'\n'
2716 <<
" last col diag : " << A->last_col_diag <<
'\n'
2717 <<
"number of nonzeros : " << A->num_nonzeros <<
'\n';
2719 hypre_CSRMatrix *csr = A->diag;
2720 const char *csr_name =
"diag";
2721 for (
int m = 0; m < 2; m++)
2723 auto csr_nnz = csr->i[csr->num_rows];
2724 os << csr_name <<
" num rows : " << csr->num_rows <<
'\n'
2725 << csr_name <<
" num cols : " << csr->num_cols <<
'\n'
2726 << csr_name <<
" num nnz : " << csr->num_nonzeros <<
'\n'
2727 << csr_name <<
" i last : " << csr_nnz
2728 << (csr_nnz == csr->num_nonzeros ?
2729 " [good]" :
" [** BAD **]") <<
'\n';
2731 os << csr_name <<
" i hash : " << hf.
GetHash() <<
'\n';
2732 os << csr_name <<
" j hash : ";
2733 if (csr->j ==
nullptr)
2742#if MFEM_HYPRE_VERSION >= 21600
2743 os << csr_name <<
" big j hash : ";
2744 if (csr->big_j ==
nullptr)
2754 os << csr_name <<
" data hash : ";
2755 if (csr->data ==
nullptr)
2769 hf.
AppendInts(A->col_map_offd, A->offd->num_cols);
2770 os <<
"col map offd hash : " << hf.
GetHash() <<
'\n';
2775 HYPRE_BigInt *A_col_map_offd = hypre_ParCSRMatrixColMapOffd(A);
2776 int size = hypre_CSRMatrixNumCols(hypre_ParCSRMatrixOffd(A));
2780void HypreParMatrix::Destroy()
2782 if ( X != NULL ) {
delete X; }
2783 if ( Y != NULL ) {
delete Y; }
2787 if (A == NULL) {
return; }
2789#ifdef HYPRE_USING_GPU
2790 if (
HypreUsingGPU() && ParCSROwner && (diagOwner < 0 || offdOwner < 0))
2798 MFEM_VERIFY(!(diagOwner < 0 && offdOwner < 0) || diagOwner == offdOwner,
2801 MemoryClass mc = (diagOwner == -1 || offdOwner == -1) ?
2803 Write(mc, diagOwner < 0, offdOwner <0);
2812 hypre_CSRMatrixI(A->diag) = NULL;
2813 hypre_CSRMatrixJ(A->diag) = NULL;
2814 hypre_CSRMatrixData(A->diag) = NULL;
2821 hypre_CSRMatrixI(A->offd) = NULL;
2822 hypre_CSRMatrixJ(A->offd) = NULL;
2823 hypre_CSRMatrixData(A->offd) = NULL;
2825 if (colMapOwner >= 0)
2827 if (colMapOwner & 1)
2831 hypre_ParCSRMatrixColMapOffd(A) = NULL;
2836 hypre_ParCSRMatrixDestroy(A);
2845 MFEM_CONTRACT_VAR(own_j);
2846 MFEM_ASSERT(own_i == own_j,
"Inconsistent ownership");
2860#if MFEM_HYPRE_VERSION >= 21800
2869 hypre_ParCSRMatrix *C_hypre;
2870 hypre_ParcsrBdiagInvScal(*A, blocksize, &C_hypre);
2871 hypre_ParCSRMatrixDropSmallEntries(C_hypre, 1e-15, 1);
2881 hypre_ParVector *d_hypre;
2882 hypre_ParvecBdiagInvScal(b_Hypre, blocksize, &d_hypre, *A);
2890#if MFEM_HYPRE_VERSION < 21400
2895 hypre_ParCSRMatrix *C_hypre =
2898 MFEM_VERIFY(C_hypre,
"error in hypre_ParCSRMatrixAdd");
2900 if (!hypre_ParCSRMatrixCommPkg(C_hypre)) { hypre_MatvecCommPkgCreate(C_hypre); }
2911 hypre_ParCSRMatrix * C = internal::hypre_ParCSRMatrixAdd(*A,*B);
2913 if (!hypre_ParCSRMatrixCommPkg(C)) { hypre_MatvecCommPkgCreate(C); }
2923 hypre_ParCSRMatrix *C;
2924#if MFEM_HYPRE_VERSION <= 22000
2927 hypre_ParCSRMatrixAdd(
alpha, A,
beta, B, &C);
2929 if (!hypre_ParCSRMatrixCommPkg(C)) { hypre_MatvecCommPkgCreate(C); }
2931 return new HypreParMatrix(C);
2934HypreParMatrix *
ParAdd(
const HypreParMatrix *A,
const HypreParMatrix *B)
2936 hypre_ParCSRMatrix *C;
2937#if MFEM_HYPRE_VERSION <= 22000
2938 hypre_ParcsrAdd(1.0, *A, 1.0, *B, &C);
2940 hypre_ParCSRMatrixAdd(1.0, *A, 1.0, *B, &C);
2942 if (!hypre_ParCSRMatrixCommPkg(C)) { hypre_MatvecCommPkgCreate(C); }
2944 return new HypreParMatrix(C);
2952 hypre_ParCSRMatrix * ab;
2953#ifdef HYPRE_USING_GPU
2956 ab = hypre_ParCSRMatMat(*A, *B);
2961 ab = hypre_ParMatmul(*A,*B);
2963 hypre_ParCSRMatrixSetNumNonzeros(ab);
2965 if (!hypre_ParCSRMatrixCommPkg(ab)) { hypre_MatvecCommPkgCreate(ab); }
2977 hypre_ParCSRMatrix * rap;
2979#ifdef HYPRE_USING_GPU
2988 hypre_ParCSRMatrix *Q = hypre_ParCSRMatMat(*A,*P);
2989 const bool keepTranspose =
false;
2990 rap = hypre_ParCSRTMatMatKT(*P,Q,keepTranspose);
2991 hypre_ParCSRMatrixDestroy(Q);
2999#if MFEM_HYPRE_VERSION <= 22200
3000 HYPRE_Int P_owns_its_col_starts =
3001 hypre_ParCSRMatrixOwnsColStarts((hypre_ParCSRMatrix*)(*P));
3004 hypre_BoomerAMGBuildCoarseOperator(*P,*A,*P,&rap);
3006#if MFEM_HYPRE_VERSION <= 22200
3009 hypre_ParCSRMatrixSetRowStartsOwner(rap,0);
3010 hypre_ParCSRMatrixSetColStartsOwner(rap,0);
3011 if (P_owns_its_col_starts)
3013 hypre_ParCSRMatrixSetColStartsOwner(*P, 1);
3018 hypre_ParCSRMatrixSetNumNonzeros(rap);
3027 hypre_ParCSRMatrix * rap;
3029#ifdef HYPRE_USING_GPU
3032 hypre_ParCSRMatrix *Q = hypre_ParCSRMatMat(*A,*P);
3033 rap = hypre_ParCSRTMatMat(*Rt,Q);
3034 hypre_ParCSRMatrixDestroy(Q);
3039#if MFEM_HYPRE_VERSION <= 22200
3040 HYPRE_Int P_owns_its_col_starts =
3041 hypre_ParCSRMatrixOwnsColStarts((hypre_ParCSRMatrix*)(*P));
3042 HYPRE_Int Rt_owns_its_col_starts =
3043 hypre_ParCSRMatrixOwnsColStarts((hypre_ParCSRMatrix*)(*Rt));
3046 hypre_BoomerAMGBuildCoarseOperator(*Rt,*A,*P,&rap);
3048#if MFEM_HYPRE_VERSION <= 22200
3051 hypre_ParCSRMatrixSetRowStartsOwner(rap,0);
3052 hypre_ParCSRMatrixSetColStartsOwner(rap,0);
3053 if (P_owns_its_col_starts)
3055 hypre_ParCSRMatrixSetColStartsOwner(*P, 1);
3057 if (Rt_owns_its_col_starts)
3059 hypre_ParCSRMatrixSetColStartsOwner(*Rt, 1);
3064 hypre_ParCSRMatrixSetNumNonzeros(rap);
3073 const int num_loc,
const Array<int> &offsets,
3074 std::vector<int> &all_num_loc,
const int numBlocks,
3075 std::vector<std::vector<HYPRE_BigInt>> &blockProcOffsets,
3076 std::vector<HYPRE_BigInt> &procOffsets,
3077 std::vector<std::vector<int>> &procBlockOffsets,
3080 std::vector<std::vector<int>> all_block_num_loc(numBlocks);
3082 MPI_Allgather(&num_loc, 1, MPI_INT, all_num_loc.data(), 1, MPI_INT, comm);
3084 for (
int j = 0; j < numBlocks; ++j)
3086 all_block_num_loc[j].resize(nprocs);
3087 blockProcOffsets[j].resize(nprocs);
3089 const int blockNumRows = offsets[j + 1] - offsets[j];
3090 MPI_Allgather(&blockNumRows, 1, MPI_INT, all_block_num_loc[j].data(), 1,
3092 blockProcOffsets[j][0] = 0;
3093 for (
int i = 0; i < nprocs - 1; ++i)
3095 blockProcOffsets[j][i + 1] = blockProcOffsets[j][i]
3096 + all_block_num_loc[j][i];
3103 for (
int i = 0; i < nprocs; ++i)
3105 globalNum += all_num_loc[i];
3106 MFEM_VERIFY(globalNum >= 0,
"overflow in global size");
3109 firstLocal += all_num_loc[i];
3114 procOffsets[i + 1] = procOffsets[i] + all_num_loc[i];
3117 procBlockOffsets[i].resize(numBlocks);
3118 procBlockOffsets[i][0] = 0;
3119 for (
int j = 1; j < numBlocks; ++j)
3121 procBlockOffsets[i][j] = procBlockOffsets[i][j - 1]
3122 + all_block_num_loc[j - 1][i];
3130 const int numBlockRows = blocks.
NumRows();
3131 const int numBlockCols = blocks.
NumCols();
3133 MFEM_VERIFY(numBlockRows > 0 &&
3134 numBlockCols > 0,
"Invalid input to HypreParMatrixFromBlocks");
3136 if (blockCoeff != NULL)
3138 MFEM_VERIFY(numBlockRows == blockCoeff->
NumRows() &&
3139 numBlockCols == blockCoeff->
NumCols(),
3140 "Invalid input to HypreParMatrixFromBlocks");
3146 int nonNullBlockRow0 = -1;
3147 for (
int j=0; j<numBlockCols; ++j)
3149 if (blocks(0,j) != NULL)
3151 nonNullBlockRow0 = j;
3156 MFEM_VERIFY(nonNullBlockRow0 >= 0,
"Null row of blocks");
3157 MPI_Comm comm = blocks(0,nonNullBlockRow0)->GetComm();
3162 for (
int i=0; i<numBlockRows; ++i)
3164 for (
int j=0; j<numBlockCols; ++j)
3166 if (blocks(i,j) != NULL)
3168 const int nrows = blocks(i,j)->
NumRows();
3169 const int ncols = blocks(i,j)->
NumCols();
3171 if (rowOffsets[i+1] == 0)
3173 rowOffsets[i+1] = nrows;
3177 MFEM_VERIFY(rowOffsets[i+1] == nrows,
3178 "Inconsistent blocks in HypreParMatrixFromBlocks");
3181 if (colOffsets[j+1] == 0)
3183 colOffsets[j+1] = ncols;
3187 MFEM_VERIFY(colOffsets[j+1] == ncols,
3188 "Inconsistent blocks in HypreParMatrixFromBlocks");
3192 rowOffsets[i+1] += rowOffsets[i];
3195 for (
int j=0; j<numBlockCols; ++j)
3197 colOffsets[j+1] += colOffsets[j];
3200 const int num_loc_rows = rowOffsets[numBlockRows];
3201 const int num_loc_cols = colOffsets[numBlockCols];
3204 MPI_Comm_rank(comm, &rank);
3205 MPI_Comm_size(comm, &nprocs);
3207 std::vector<int> all_num_loc_rows(nprocs);
3208 std::vector<int> all_num_loc_cols(nprocs);
3209 std::vector<HYPRE_BigInt> procRowOffsets(nprocs);
3210 std::vector<HYPRE_BigInt> procColOffsets(nprocs);
3211 std::vector<std::vector<HYPRE_BigInt>> blockRowProcOffsets(numBlockRows);
3212 std::vector<std::vector<HYPRE_BigInt>> blockColProcOffsets(numBlockCols);
3213 std::vector<std::vector<int>> procBlockRowOffsets(nprocs);
3214 std::vector<std::vector<int>> procBlockColOffsets(nprocs);
3216 HYPRE_BigInt first_loc_row, glob_nrows, first_loc_col, glob_ncols;
3218 all_num_loc_rows, numBlockRows, blockRowProcOffsets,
3219 procRowOffsets, procBlockRowOffsets, first_loc_row,
3223 all_num_loc_cols, numBlockCols, blockColProcOffsets,
3224 procColOffsets, procBlockColOffsets, first_loc_col,
3227 std::vector<int> opI(num_loc_rows + 1);
3228 std::vector<int> cnt(num_loc_rows);
3230 for (
int i = 0; i < num_loc_rows; ++i)
3236 opI[num_loc_rows] = 0;
3241 for (
int i = 0; i < numBlockRows; ++i)
3243 for (
int j = 0; j < numBlockCols; ++j)
3245 if (blocks(i, j) == NULL)
3247 csr_blocks(i, j) = NULL;
3251 blocks(i, j)->HostRead();
3252 csr_blocks(i, j) = hypre_MergeDiagAndOffd(*blocks(i, j));
3253 blocks(i, j)->HypreRead();
3255 for (
int k = 0; k < csr_blocks(i, j)->num_rows; ++k)
3257 opI[rowOffsets[i] + k + 1] +=
3258 csr_blocks(i, j)->i[k + 1] - csr_blocks(i, j)->i[k];
3265 for (
int i = 0; i < num_loc_rows; ++i)
3267 opI[i + 1] += opI[i];
3270 const int nnz = opI[num_loc_rows];
3272 std::vector<HYPRE_BigInt> opJ(nnz);
3273 std::vector<real_t> data(nnz);
3276 for (
int i = 0; i < numBlockRows; ++i)
3278 for (
int j = 0; j < numBlockCols; ++j)
3280 if (csr_blocks(i, j) != NULL)
3282 const int nrows = csr_blocks(i, j)->num_rows;
3283 const real_t cij = blockCoeff ? (*blockCoeff)(i, j) : 1.0;
3284#if MFEM_HYPRE_VERSION >= 21600
3285 const bool usingBigJ = (csr_blocks(i, j)->big_j != NULL);
3288 for (
int k = 0; k < nrows; ++k)
3290 const int rowg = rowOffsets[i] + k;
3291 const int nnz_k = csr_blocks(i,j)->i[k+1]-csr_blocks(i,j)->i[k];
3292 const int osk = csr_blocks(i, j)->i[k];
3294 for (
int l = 0; l < nnz_k; ++l)
3297#if MFEM_HYPRE_VERSION >= 21600
3298 const HYPRE_Int bcol = usingBigJ ?
3299 csr_blocks(i, j)->big_j[osk + l] :
3300 csr_blocks(i, j)->j[osk + l];
3302 const HYPRE_Int bcol = csr_blocks(i, j)->j[osk + l];
3306 const auto &offs = blockColProcOffsets[j];
3307 const int bcolproc =
3308 std::upper_bound(offs.begin() + 1, offs.end(), bcol)
3311 opJ[opI[rowg] + cnt[rowg]] = procColOffsets[bcolproc] +
3312 procBlockColOffsets[bcolproc][j]
3314 - blockColProcOffsets[j][bcolproc];
3315 data[opI[rowg] + cnt[rowg]] = cij * csr_blocks(i, j)->data[osk + l];
3323 for (
int i = 0; i < numBlockRows; ++i)
3325 for (
int j = 0; j < numBlockCols; ++j)
3327 if (csr_blocks(i, j) != NULL)
3329 hypre_CSRMatrixDestroy(csr_blocks(i, j));
3334 MFEM_VERIFY(HYPRE_AssumedPartitionCheck(),
3335 "only 'assumed partition' mode is supported");
3337 std::vector<HYPRE_BigInt> rowStarts2(2);
3338 rowStarts2[0] = first_loc_row;
3339 rowStarts2[1] = first_loc_row + all_num_loc_rows[rank];
3341 int square = std::equal(all_num_loc_rows.begin(), all_num_loc_rows.end(),
3342 all_num_loc_cols.begin());
3345 return new HypreParMatrix(comm, num_loc_rows, glob_nrows, glob_ncols,
3346 opI.data(), opJ.data(),
3353 std::vector<HYPRE_BigInt> colStarts2(2);
3354 colStarts2[0] = first_loc_col;
3355 colStarts2[1] = first_loc_col + all_num_loc_cols[rank];
3357 return new HypreParMatrix(comm, num_loc_rows, glob_nrows, glob_ncols,
3358 opI.data(), opJ.data(),
3369 for (
int i = 0; i < blocks.
NumRows(); ++i)
3371 for (
int j = 0; j < blocks.
NumCols(); ++j)
3373 constBlocks(i, j) = blocks(i, j);
3399 hypre_CSRMatrix *A_diag = hypre_ParCSRMatrixDiag(A);
3400 HYPRE_Int num_rows = hypre_CSRMatrixNumRows(A_diag);
3402 real_t *u_data = hypre_VectorData(hypre_ParVectorLocalVector(
u));
3403 real_t *r_data = hypre_VectorData(hypre_ParVectorLocalVector(r));
3405 for (
int i = 0; i < N; i++)
3408 hypre_ParVectorCopy(
f, r);
3409 hypre_ParCSRMatrixMatvec(-1.0, A,
u, 1.0, r);
3412 (0 == (i % 2)) ? coef = lambda : coef =
mu;
3414 for (HYPRE_Int j = 0; j < num_rows; j++)
3416 u_data[j] += coef*r_data[j] / max_eig;
3432 hypre_ParVector *x0,
3433 hypre_ParVector *x1,
3434 hypre_ParVector *x2,
3435 hypre_ParVector *x3)
3438 hypre_CSRMatrix *A_diag = hypre_ParCSRMatrixDiag(A);
3439 HYPRE_Int num_rows = hypre_CSRMatrixNumRows(A_diag);
3441 real_t *u_data = hypre_VectorData(hypre_ParVectorLocalVector(
u));
3443 real_t *x0_data = hypre_VectorData(hypre_ParVectorLocalVector(x0));
3444 real_t *x1_data = hypre_VectorData(hypre_ParVectorLocalVector(x1));
3445 real_t *x2_data = hypre_VectorData(hypre_ParVectorLocalVector(x2));
3446 real_t *x3_data = hypre_VectorData(hypre_ParVectorLocalVector(x3));
3448 hypre_ParVectorCopy(
u, x0);
3451 hypre_ParVectorCopy(
f, x1);
3452 hypre_ParCSRMatrixMatvec(-1.0, A, x0, 1.0, x1);
3454 for (HYPRE_Int i = 0; i < num_rows; i++)
3456 x1_data[i] /= -max_eig;
3460 for (HYPRE_Int i = 0; i < num_rows; i++)
3462 x1_data[i] = x0_data[i] -x1_data[i];
3466 for (HYPRE_Int i = 0; i < num_rows; i++)
3468 x3_data[i] = fir_coeffs[0]*x0_data[i] +fir_coeffs[1]*x1_data[i];
3471 for (
int n = 2; n <= poly_order; n++)
3474 hypre_ParVectorCopy(
f, x2);
3475 hypre_ParCSRMatrixMatvec(-1.0, A, x1, 1.0, x2);
3477 for (HYPRE_Int i = 0; i < num_rows; i++)
3479 x2_data[i] /= -max_eig;
3487 for (HYPRE_Int i = 0; i < num_rows; i++)
3489 x2_data[i] = (x1_data[i]-x0_data[i]) +(x1_data[i]-2*x2_data[i]);
3490 x3_data[i] += fir_coeffs[n]*x2_data[i];
3491 x0_data[i] = x1_data[i];
3492 x1_data[i] = x2_data[i];
3496 for (HYPRE_Int i = 0; i < num_rows; i++)
3498 u_data[i] = x3_data[i];
3519 B =
X =
V =
Z = NULL;
3527 int relax_times_,
real_t relax_weight_,
3528 real_t omega_,
int poly_order_,
3529 real_t poly_fraction_,
int eig_est_cg_iter_)
3541 B =
X =
V =
Z = NULL;
3552 type =
static_cast<int>(type_);
3563 int eig_est_cg_iter_)
3581 if (!strcmp(name,
"Rectangular")) {
a = 1.0,
b = 0.0, c = 0.0; }
3582 if (!strcmp(name,
"Hanning")) {
a = 0.5,
b = 0.5, c = 0.0; }
3583 if (!strcmp(name,
"Hamming")) {
a = 0.54,
b = 0.46, c = 0.0; }
3584 if (!strcmp(name,
"Blackman")) {
a = 0.42,
b = 0.50, c = 0.08; }
3587 mfem_error(
"HypreSmoother::SetWindowByName : name not recognized!");
3605 mfem_error(
"HypreSmoother::SetOperator : not HypreParMatrix!");
3612 if (
B) {
delete B; }
3613 if (
X) {
delete X; }
3614 if (
V) {
delete V; }
3615 if (
Z) {
delete Z; }
3637 A->
Mult(ones, diag);
3648 d_l1_norms[i] = std::abs(d_l1_norms[i]);
3652#if MFEM_HYPRE_VERSION < 22100
3664#elif defined(HYPRE_USING_GPU)
3694#if MFEM_HYPRE_VERSION <= 22200
3703 else if (
type == 1001 ||
type == 1002)
3713#if MFEM_HYPRE_VERSION <= 22200
3755 window_coeffs[i] =
a +
b*cos(
t) +c*cos(2*
t);
3759 real_t theta_pb = acos(1.0 -0.5*k_pb);
3761 cheby_coeffs[0] = (theta_pb +
sigma)/M_PI;
3765 cheby_coeffs[i] = 2.0*sin(
t)/(i*M_PI);
3770 fir_coeffs[i] = window_coeffs[i]*cheby_coeffs[i];
3773 delete[] window_coeffs;
3774 delete[] cheby_coeffs;
3781 mfem_error(
"HypreSmoother::Mult (...) : HypreParMatrix A is missing");
3794 HYPRE_ParCSRDiagScale(NULL, *
A,
b, x);
3801 hypre_ParVectorSetConstantValues(x, 0.0);
3822 else if (
type == 1002)
3836 int hypre_type =
type;
3838 if (
type == 5) { hypre_type = 1; }
3842 hypre_ParCSRRelax(*
A,
b, hypre_type,
3849 hypre_ParCSRRelax(*
A,
b, hypre_type,
3859 MFEM_ASSERT(
b.Size() ==
NumCols(),
"");
3864 mfem_error(
"HypreSmoother::Mult (...) : HypreParMatrix A is missing");
3871 A -> GetGlobalNumRows(),
3873 A -> GetRowStarts());
3875 A -> GetGlobalNumCols(),
3877 A -> GetColStarts());
3915 if (!xshallow) { x = *
X; }
3925 mfem_error(
"HypreSmoother::MultTranspose (...) : undefined!\n");
3931 if (
B) {
delete B; }
3932 if (
X) {
delete X; }
3933 if (
V) {
delete V; }
3934 if (
Z) {
delete Z; }
3943 if (
X0) {
delete X0; }
3944 if (
X1) {
delete X1; }
3959 :
Solver(A_->Height(), A_->Width())
3971 MFEM_ASSERT(
b.Size() ==
NumCols(),
"");
3974 MFEM_VERIFY(
A != NULL,
"HypreParMatrix A is missing");
4024 MFEM_VERIFY(
A != NULL,
"HypreParMatrix A is missing");
4026 HYPRE_Int err_flag =
SetupFcn()(*
this, *
A,
b, x);
4030 { MFEM_WARNING(
"Error during setup! Error code: " << err_flag); }
4034 MFEM_VERIFY(!err_flag,
"Error during setup! Error code: " << err_flag);
4036 hypre_error_flag = 0;
4044 if (!x_shallow) { x = *
X; }
4052 mfem_error(
"HypreSolver::Mult (...) : HypreParMatrix A is missing");
4059 hypre_ParVectorSetConstantValues(x, 0.0);
4071 { MFEM_WARNING(
"Error during solve! Error code: " << err_flag); }
4075 MFEM_VERIFY(!err_flag,
"Error during solve! Error code: " << err_flag);
4077 hypre_error_flag = 0;
4084 if (!x_shallow) { x = *
X; }
4089 if (
B) {
delete B; }
4090 if (
X) {
delete X; }
4100 HYPRE_ParCSRPCGCreate(comm, &pcg_solver);
4109 HYPRE_ParCSRMatrixGetComm(*
A, &comm);
4111 HYPRE_ParCSRPCGCreate(comm, &pcg_solver);
4117 MFEM_VERIFY(new_A,
"new Operator must be a HypreParMatrix!");
4138 HYPRE_PCGSetTol(pcg_solver, tol);
4143 HYPRE_PCGSetAbsoluteTol(pcg_solver, atol);
4148 HYPRE_PCGSetMaxIter(pcg_solver, max_iter);
4153 HYPRE_PCGSetLogging(pcg_solver, logging);
4158 HYPRE_ParCSRPCGSetPrintLevel(pcg_solver, print_lvl);
4163 precond = &precond_;
4165 HYPRE_ParCSRPCGSetPrecond(pcg_solver,
4173 HYPRE_PCGSetTwoNorm(pcg_solver, 1);
4174 if (res_frequency > 0)
4176 HYPRE_PCGSetRecomputeResidualP(pcg_solver, res_frequency);
4180 HYPRE_PCGSetResidualTol(pcg_solver, rtol);
4187 HYPRE_Int time_index = 0;
4188 HYPRE_Int num_iterations;
4191 HYPRE_Int print_level;
4193 HYPRE_PCGGetPrintLevel(pcg_solver, &print_level);
4194 HYPRE_ParCSRPCGSetPrintLevel(pcg_solver, print_level%3);
4196 HYPRE_ParCSRMatrixGetComm(*
A, &comm);
4201 hypre_ParVectorSetConstantValues(x, 0.0);
4209 if (print_level > 0 && print_level < 3)
4211 time_index = hypre_InitializeTiming(
"PCG Setup");
4212 hypre_BeginTiming(time_index);
4215 HYPRE_ParCSRPCGSetup(pcg_solver, *
A,
b, x);
4218 if (print_level > 0 && print_level < 3)
4220 hypre_EndTiming(time_index);
4221 hypre_PrintTiming(
"Setup phase times", comm);
4222 hypre_FinalizeTiming(time_index);
4223 hypre_ClearTiming();
4227 if (print_level > 0 && print_level < 3)
4229 time_index = hypre_InitializeTiming(
"PCG Solve");
4230 hypre_BeginTiming(time_index);
4233 HYPRE_ParCSRPCGSolve(pcg_solver, *
A,
b, x);
4235 if (print_level > 0)
4237 if (print_level < 3)
4239 hypre_EndTiming(time_index);
4240 hypre_PrintTiming(
"Solve phase times", comm);
4241 hypre_FinalizeTiming(time_index);
4242 hypre_ClearTiming();
4245 HYPRE_ParCSRPCGGetNumIterations(pcg_solver, &num_iterations);
4246 HYPRE_ParCSRPCGGetFinalRelativeResidualNorm(pcg_solver,
4249 MPI_Comm_rank(comm, &myid);
4253 mfem::out <<
"PCG Iterations = " << num_iterations << endl
4254 <<
"Final PCG Relative Residual Norm = " << final_res_norm
4258 HYPRE_ParCSRPCGSetPrintLevel(pcg_solver, print_level);
4263 HYPRE_ParCSRPCGDestroy(pcg_solver);
4271 HYPRE_ParCSRGMRESCreate(comm, &gmres_solver);
4272 SetDefaultOptions();
4282 HYPRE_ParCSRMatrixGetComm(*
A, &comm);
4284 HYPRE_ParCSRGMRESCreate(comm, &gmres_solver);
4285 SetDefaultOptions();
4288void HypreGMRES::SetDefaultOptions()
4294 HYPRE_ParCSRGMRESSetKDim(gmres_solver, k_dim);
4295 HYPRE_ParCSRGMRESSetMaxIter(gmres_solver, max_iter);
4296 HYPRE_ParCSRGMRESSetTol(gmres_solver, tol);
4302 MFEM_VERIFY(new_A,
"new Operator must be a HypreParMatrix!");
4323 HYPRE_GMRESSetTol(gmres_solver, tol);
4328 HYPRE_GMRESSetAbsoluteTol(gmres_solver, tol);
4333 HYPRE_GMRESSetMaxIter(gmres_solver, max_iter);
4338 HYPRE_GMRESSetKDim(gmres_solver, k_dim);
4343 HYPRE_GMRESSetLogging(gmres_solver, logging);
4348 HYPRE_GMRESSetPrintLevel(gmres_solver, print_lvl);
4353 precond = &precond_;
4355 HYPRE_ParCSRGMRESSetPrecond(gmres_solver,
4364 HYPRE_Int time_index = 0;
4365 HYPRE_Int num_iterations;
4368 HYPRE_Int print_level;
4370 HYPRE_GMRESGetPrintLevel(gmres_solver, &print_level);
4372 HYPRE_ParCSRMatrixGetComm(*
A, &comm);
4377 hypre_ParVectorSetConstantValues(x, 0.0);
4385 if (print_level > 0)
4387 time_index = hypre_InitializeTiming(
"GMRES Setup");
4388 hypre_BeginTiming(time_index);
4391 HYPRE_ParCSRGMRESSetup(gmres_solver, *
A,
b, x);
4394 if (print_level > 0)
4396 hypre_EndTiming(time_index);
4397 hypre_PrintTiming(
"Setup phase times", comm);
4398 hypre_FinalizeTiming(time_index);
4399 hypre_ClearTiming();
4403 if (print_level > 0)
4405 time_index = hypre_InitializeTiming(
"GMRES Solve");
4406 hypre_BeginTiming(time_index);
4409 HYPRE_ParCSRGMRESSolve(gmres_solver, *
A,
b, x);
4411 if (print_level > 0)
4413 hypre_EndTiming(time_index);
4414 hypre_PrintTiming(
"Solve phase times", comm);
4415 hypre_FinalizeTiming(time_index);
4416 hypre_ClearTiming();
4418 HYPRE_ParCSRGMRESGetNumIterations(gmres_solver, &num_iterations);
4419 HYPRE_ParCSRGMRESGetFinalRelativeResidualNorm(gmres_solver,
4422 MPI_Comm_rank(comm, &myid);
4426 mfem::out <<
"GMRES Iterations = " << num_iterations << endl
4427 <<
"Final GMRES Relative Residual Norm = " << final_res_norm
4435 HYPRE_ParCSRGMRESDestroy(gmres_solver);
4443 HYPRE_ParCSRFlexGMRESCreate(comm, &fgmres_solver);
4444 SetDefaultOptions();
4454 HYPRE_ParCSRMatrixGetComm(*
A, &comm);
4456 HYPRE_ParCSRFlexGMRESCreate(comm, &fgmres_solver);
4457 SetDefaultOptions();
4460void HypreFGMRES::SetDefaultOptions()
4466 HYPRE_ParCSRFlexGMRESSetKDim(fgmres_solver, k_dim);
4467 HYPRE_ParCSRFlexGMRESSetMaxIter(fgmres_solver, max_iter);
4468 HYPRE_ParCSRFlexGMRESSetTol(fgmres_solver, tol);
4474 MFEM_VERIFY(new_A,
"new Operator must be a HypreParMatrix!");
4495 HYPRE_ParCSRFlexGMRESSetTol(fgmres_solver, tol);
4500 HYPRE_ParCSRFlexGMRESSetMaxIter(fgmres_solver, max_iter);
4505 HYPRE_ParCSRFlexGMRESSetKDim(fgmres_solver, k_dim);
4510 HYPRE_ParCSRFlexGMRESSetLogging(fgmres_solver, logging);
4515 HYPRE_ParCSRFlexGMRESSetPrintLevel(fgmres_solver, print_lvl);
4520 precond = &precond_;
4521 HYPRE_ParCSRFlexGMRESSetPrecond(fgmres_solver,
4530 HYPRE_Int time_index = 0;
4531 HYPRE_Int num_iterations;
4534 HYPRE_Int print_level;
4536 HYPRE_FlexGMRESGetPrintLevel(fgmres_solver, &print_level);
4538 HYPRE_ParCSRMatrixGetComm(*
A, &comm);
4543 hypre_ParVectorSetConstantValues(x, 0.0);
4551 if (print_level > 0)
4553 time_index = hypre_InitializeTiming(
"FGMRES Setup");
4554 hypre_BeginTiming(time_index);
4557 HYPRE_ParCSRFlexGMRESSetup(fgmres_solver, *
A,
b, x);
4560 if (print_level > 0)
4562 hypre_EndTiming(time_index);
4563 hypre_PrintTiming(
"Setup phase times", comm);
4564 hypre_FinalizeTiming(time_index);
4565 hypre_ClearTiming();
4569 if (print_level > 0)
4571 time_index = hypre_InitializeTiming(
"FGMRES Solve");
4572 hypre_BeginTiming(time_index);
4575 HYPRE_ParCSRFlexGMRESSolve(fgmres_solver, *
A,
b, x);
4577 if (print_level > 0)
4579 hypre_EndTiming(time_index);
4580 hypre_PrintTiming(
"Solve phase times", comm);
4581 hypre_FinalizeTiming(time_index);
4582 hypre_ClearTiming();
4584 HYPRE_ParCSRFlexGMRESGetNumIterations(fgmres_solver, &num_iterations);
4585 HYPRE_ParCSRFlexGMRESGetFinalRelativeResidualNorm(fgmres_solver,
4588 MPI_Comm_rank(comm, &myid);
4592 mfem::out <<
"FGMRES Iterations = " << num_iterations << endl
4593 <<
"Final FGMRES Relative Residual Norm = " << final_res_norm
4601 HYPRE_ParCSRFlexGMRESDestroy(fgmres_solver);
4608 MFEM_VERIFY(new_A,
"new Operator must be a HypreParMatrix!");
4625 HYPRE_ParaSailsCreate(comm, &sai_precond);
4626 SetDefaultOptions();
4633 HYPRE_ParCSRMatrixGetComm(
A, &comm);
4635 HYPRE_ParaSailsCreate(comm, &sai_precond);
4636 SetDefaultOptions();
4639void HypreParaSails::SetDefaultOptions()
4641 int sai_max_levels = 1;
4642 real_t sai_threshold = 0.1;
4645 real_t sai_loadbal = 0.0;
4647 int sai_logging = 1;
4649 HYPRE_ParaSailsSetParams(sai_precond, sai_threshold, sai_max_levels);
4650 HYPRE_ParaSailsSetFilter(sai_precond, sai_filter);
4651 HYPRE_ParaSailsSetSym(sai_precond, sai_sym);
4652 HYPRE_ParaSailsSetLoadbal(sai_precond, sai_loadbal);
4653 HYPRE_ParaSailsSetReuse(sai_precond, sai_reuse);
4654 HYPRE_ParaSailsSetLogging(sai_precond, sai_logging);
4657void HypreParaSails::ResetSAIPrecond(MPI_Comm comm)
4659 HYPRE_Int sai_max_levels;
4660 HYPRE_Real sai_threshold;
4661 HYPRE_Real sai_filter;
4663 HYPRE_Real sai_loadbal;
4664 HYPRE_Int sai_reuse;
4665 HYPRE_Int sai_logging;
4668 HYPRE_ParaSailsGetNlevels(sai_precond, &sai_max_levels);
4669 HYPRE_ParaSailsGetThresh(sai_precond, &sai_threshold);
4670 HYPRE_ParaSailsGetFilter(sai_precond, &sai_filter);
4671 HYPRE_ParaSailsGetSym(sai_precond, &sai_sym);
4672 HYPRE_ParaSailsGetLoadbal(sai_precond, &sai_loadbal);
4673 HYPRE_ParaSailsGetReuse(sai_precond, &sai_reuse);
4674 HYPRE_ParaSailsGetLogging(sai_precond, &sai_logging);
4676 HYPRE_ParaSailsDestroy(sai_precond);
4677 HYPRE_ParaSailsCreate(comm, &sai_precond);
4679 HYPRE_ParaSailsSetParams(sai_precond, sai_threshold, sai_max_levels);
4680 HYPRE_ParaSailsSetFilter(sai_precond, sai_filter);
4681 HYPRE_ParaSailsSetSym(sai_precond, sai_sym);
4682 HYPRE_ParaSailsSetLoadbal(sai_precond, sai_loadbal);
4683 HYPRE_ParaSailsSetReuse(sai_precond, sai_reuse);
4684 HYPRE_ParaSailsSetLogging(sai_precond, sai_logging);
4690 MFEM_VERIFY(new_A,
"new Operator must be a HypreParMatrix!");
4695 HYPRE_ParCSRMatrixGetComm(*
A, &comm);
4696 ResetSAIPrecond(comm);
4713 HYPRE_ParaSailsSetParams(sai_precond, threshold, max_levels);
4718 HYPRE_ParaSailsSetFilter(sai_precond, filter);
4723 HYPRE_ParaSailsSetSym(sai_precond, sym);
4728 HYPRE_ParaSailsSetLoadbal(sai_precond, loadbal);
4733 HYPRE_ParaSailsSetReuse(sai_precond, reuse);
4738 HYPRE_ParaSailsSetLogging(sai_precond, logging);
4743 HYPRE_ParaSailsDestroy(sai_precond);
4749 HYPRE_EuclidCreate(comm, &euc_precond);
4750 SetDefaultOptions();
4757 HYPRE_ParCSRMatrixGetComm(
A, &comm);
4759 HYPRE_EuclidCreate(comm, &euc_precond);
4760 SetDefaultOptions();
4763void HypreEuclid::SetDefaultOptions()
4771 HYPRE_EuclidSetLevel(euc_precond, euc_level);
4772 HYPRE_EuclidSetStats(euc_precond, euc_stats);
4773 HYPRE_EuclidSetMem(euc_precond, euc_mem);
4774 HYPRE_EuclidSetBJ(euc_precond, euc_bj);
4775 HYPRE_EuclidSetRowScale(euc_precond, euc_ro_sc);
4780 HYPRE_EuclidSetLevel(euc_precond, level);
4785 HYPRE_EuclidSetStats(euc_precond, stats);
4790 HYPRE_EuclidSetMem(euc_precond, mem);
4795 HYPRE_EuclidSetBJ(euc_precond, bj);
4800 HYPRE_EuclidSetRowScale(euc_precond, row_scale);
4803void HypreEuclid::ResetEuclidPrecond(MPI_Comm comm)
4807 HYPRE_EuclidDestroy(euc_precond);
4808 HYPRE_EuclidCreate(comm, &euc_precond);
4810 SetDefaultOptions();
4816 MFEM_VERIFY(new_A,
"new Operator must be a HypreParMatrix!");
4821 HYPRE_ParCSRMatrixGetComm(*new_A, &comm);
4822 ResetEuclidPrecond(comm);
4839 HYPRE_EuclidDestroy(euc_precond);
4843#if MFEM_HYPRE_VERSION >= 21900
4846 HYPRE_ILUCreate(&ilu_precond);
4847 SetDefaultOptions();
4850void HypreILU::SetDefaultOptions()
4853 HYPRE_Int ilu_type = 0;
4854 HYPRE_ILUSetType(ilu_precond, ilu_type);
4857 HYPRE_Int max_iter = 1;
4858 HYPRE_ILUSetMaxIter(ilu_precond, max_iter);
4861 HYPRE_Real tol = 0.0;
4862 HYPRE_ILUSetTol(ilu_precond, tol);
4865 HYPRE_Int lev_fill = 1;
4866 HYPRE_ILUSetLevelOfFill(ilu_precond, lev_fill);
4869 HYPRE_Int reorder_type = 1;
4870 HYPRE_ILUSetLocalReordering(ilu_precond, reorder_type);
4873 HYPRE_Int print_level = 0;
4874 HYPRE_ILUSetPrintLevel(ilu_precond, print_level);
4877void HypreILU::ResetILUPrecond()
4881 HYPRE_ILUDestroy(ilu_precond);
4883 HYPRE_ILUCreate(&ilu_precond);
4884 SetDefaultOptions();
4889 HYPRE_ILUSetLevelOfFill(ilu_precond, lev_fill);
4894 HYPRE_ILUSetType(ilu_precond, ilu_type);
4899 HYPRE_ILUSetMaxIter(ilu_precond, max_iter);
4904 HYPRE_ILUSetTol(ilu_precond, tol);
4909 HYPRE_ILUSetLocalReordering(ilu_precond, reorder_type);
4914 HYPRE_ILUSetPrintLevel(ilu_precond, print_level);
4920 MFEM_VERIFY(new_A,
"new Operator must be a HypreParMatrix!");
4922 if (
A) { ResetILUPrecond(); }
4938 HYPRE_ILUDestroy(ilu_precond);
4945 HYPRE_BoomerAMGCreate(&amg_precond);
4946 SetDefaultOptions();
4951 HYPRE_BoomerAMGCreate(&amg_precond);
4952 SetDefaultOptions();
4955void HypreBoomerAMG::SetDefaultOptions()
4958 int coarsen_type, agg_levels, interp_type, Pmax, relax_type, relax_sweeps,
4959 print_level, max_levels;
5001 HYPRE_BoomerAMGSetCoarsenType(amg_precond, coarsen_type);
5002 HYPRE_BoomerAMGSetAggNumLevels(amg_precond, agg_levels);
5003 HYPRE_BoomerAMGSetRelaxType(amg_precond, relax_type);
5006 HYPRE_BoomerAMGSetNumSweeps(amg_precond, relax_sweeps);
5007 HYPRE_BoomerAMGSetStrongThreshold(amg_precond, theta);
5008 HYPRE_BoomerAMGSetInterpType(amg_precond, interp_type);
5009 HYPRE_BoomerAMGSetPMaxElmts(amg_precond, Pmax);
5010 HYPRE_BoomerAMGSetPrintLevel(amg_precond, print_level);
5011 HYPRE_BoomerAMGSetMaxLevels(amg_precond, max_levels);
5014 HYPRE_BoomerAMGSetMaxIter(amg_precond, 1);
5015 HYPRE_BoomerAMGSetTol(amg_precond, 0.0);
5018void HypreBoomerAMG::ResetAMGPrecond()
5020 HYPRE_Int coarsen_type;
5021 HYPRE_Int agg_levels;
5022 HYPRE_Int relax_type;
5023 HYPRE_Int relax_sweeps;
5025 HYPRE_Int interp_type;
5027 HYPRE_Int print_level;
5028 HYPRE_Int max_levels;
5030 HYPRE_Int nrbms = rbms.
Size();
5032 HYPRE_Int nodal_diag;
5033 HYPRE_Int relax_coarse;
5034 HYPRE_Int interp_vec_variant;
5036 HYPRE_Int smooth_interp_vectors;
5037 HYPRE_Int interp_refine;
5039 hypre_ParAMGData *amg_data = (hypre_ParAMGData *)amg_precond;
5042 HYPRE_BoomerAMGGetCoarsenType(amg_precond, &coarsen_type);
5043 agg_levels = hypre_ParAMGDataAggNumLevels(amg_data);
5044 relax_type = hypre_ParAMGDataUserRelaxType(amg_data);
5045 relax_sweeps = hypre_ParAMGDataUserNumSweeps(amg_data);
5046 HYPRE_BoomerAMGGetStrongThreshold(amg_precond, &theta);
5047 hypre_BoomerAMGGetInterpType(amg_precond, &interp_type);
5048 HYPRE_BoomerAMGGetPMaxElmts(amg_precond, &Pmax);
5049 HYPRE_BoomerAMGGetPrintLevel(amg_precond, &print_level);
5050 HYPRE_BoomerAMGGetMaxLevels(amg_precond, &max_levels);
5051 HYPRE_BoomerAMGGetNumFunctions(amg_precond, &
dim);
5054 nodal = hypre_ParAMGDataNodal(amg_data);
5055 nodal_diag = hypre_ParAMGDataNodalDiag(amg_data);
5056 HYPRE_BoomerAMGGetCycleRelaxType(amg_precond, &relax_coarse, 3);
5057 interp_vec_variant = hypre_ParAMGInterpVecVariant(amg_data);
5058 q_max = hypre_ParAMGInterpVecQMax(amg_data);
5059 smooth_interp_vectors = hypre_ParAMGSmoothInterpVectors(amg_data);
5060 interp_refine = hypre_ParAMGInterpRefine(amg_data);
5063 HYPRE_BoomerAMGDestroy(amg_precond);
5064 HYPRE_BoomerAMGCreate(&amg_precond);
5066 HYPRE_BoomerAMGSetCoarsenType(amg_precond, coarsen_type);
5067 HYPRE_BoomerAMGSetAggNumLevels(amg_precond, agg_levels);
5068 HYPRE_BoomerAMGSetRelaxType(amg_precond, relax_type);
5069 HYPRE_BoomerAMGSetNumSweeps(amg_precond, relax_sweeps);
5070 HYPRE_BoomerAMGSetMaxLevels(amg_precond, max_levels);
5071 HYPRE_BoomerAMGSetTol(amg_precond, 0.0);
5072 HYPRE_BoomerAMGSetMaxIter(amg_precond, 1);
5073 HYPRE_BoomerAMGSetStrongThreshold(amg_precond, theta);
5074 HYPRE_BoomerAMGSetInterpType(amg_precond, interp_type);
5075 HYPRE_BoomerAMGSetPMaxElmts(amg_precond, Pmax);
5076 HYPRE_BoomerAMGSetPrintLevel(amg_precond, print_level);
5077 HYPRE_BoomerAMGSetNumFunctions(amg_precond,
dim);
5080 HYPRE_BoomerAMGSetNodal(amg_precond, nodal);
5081 HYPRE_BoomerAMGSetNodalDiag(amg_precond, nodal_diag);
5082 HYPRE_BoomerAMGSetCycleRelaxType(amg_precond, relax_coarse, 3);
5083 HYPRE_BoomerAMGSetInterpVecVariant(amg_precond, interp_vec_variant);
5084 HYPRE_BoomerAMGSetInterpVecQMax(amg_precond, q_max);
5085 HYPRE_BoomerAMGSetSmoothInterpVectors(amg_precond, smooth_interp_vectors);
5086 HYPRE_BoomerAMGSetInterpRefine(amg_precond, interp_refine);
5088 HYPRE_BoomerAMGSetInterpVectors(amg_precond, rbms.
Size(), rbms.
GetData());
5095 MFEM_VERIFY(new_A,
"new Operator must be a HypreParMatrix!");
5097 if (
A) { ResetAMGPrecond(); }
5113 HYPRE_BoomerAMGSetNumFunctions(amg_precond,
dim);
5121 HYPRE_Int *h_mapping = mfem_hypre_CTAlloc_host(HYPRE_Int,
height);
5123 MFEM_VERIFY(
height %
dim == 0,
"Ordering does not work as claimed!");
5125 for (
int i = 0; i <
dim; ++i)
5127 for (
int j = 0; j < h_nnodes; ++j)
5136 HYPRE_Int *mapping =
nullptr;
5137#if defined(hypre_IntArrayData) && defined(HYPRE_USING_GPU)
5140 mapping = mfem_hypre_CTAlloc(HYPRE_Int,
height);
5141 hypre_TMemcpy(mapping, h_mapping, HYPRE_Int,
height,
5142 HYPRE_MEMORY_DEVICE, HYPRE_MEMORY_HOST);
5143 mfem_hypre_TFree_host(h_mapping);
5148 mapping = h_mapping;
5153 HYPRE_BoomerAMGSetDofFunc(amg_precond, mapping);
5157 HYPRE_BoomerAMGSetAggNumLevels(amg_precond, 0);
5158 HYPRE_BoomerAMGSetStrongThreshold(amg_precond, 0.5);
5164 y = 0.0; y(0) = x(1); y(1) = -x(0);
5166static void func_ryz(
const Vector &x, Vector &y)
5168 y = 0.0; y(1) = x(2); y(2) = -x(1);
5170static void func_rzx(
const Vector &x, Vector &y)
5172 y = 0.0; y(2) = x(0); y(0) = -x(2);
5175void HypreBoomerAMG::RecomputeRBMs()
5178 Array<HypreParVector*> gf_rbms;
5181 for (
int i = 0; i < rbms.
Size(); i++)
5183 HYPRE_ParVectorDestroy(rbms[i]);
5190 VectorFunctionCoefficient coeff_rxy(2, func_rxy);
5192 ParGridFunction rbms_rxy(fespace);
5193 rbms_rxy.ProjectCoefficient(coeff_rxy);
5196 gf_rbms.SetSize(nrbms);
5198 rbms_rxy.GetTrueDofs(*gf_rbms[0]);
5204 VectorFunctionCoefficient coeff_rxy(3, func_rxy);
5205 VectorFunctionCoefficient coeff_ryz(3, func_ryz);
5206 VectorFunctionCoefficient coeff_rzx(3, func_rzx);
5208 ParGridFunction rbms_rxy(fespace);
5209 ParGridFunction rbms_ryz(fespace);
5210 ParGridFunction rbms_rzx(fespace);
5211 rbms_rxy.ProjectCoefficient(coeff_rxy);
5212 rbms_ryz.ProjectCoefficient(coeff_ryz);
5213 rbms_rzx.ProjectCoefficient(coeff_rzx);
5216 gf_rbms.SetSize(nrbms);
5220 rbms_rxy.GetTrueDofs(*gf_rbms[0]);
5221 rbms_ryz.GetTrueDofs(*gf_rbms[1]);
5222 rbms_rzx.GetTrueDofs(*gf_rbms[2]);
5231 for (
int i = 0; i < nrbms; i++)
5233 rbms[i] = gf_rbms[i]->StealParVector();
5239 bool interp_refine_)
5241#ifdef HYPRE_USING_GPU
5244 MFEM_ABORT(
"this method is not supported in hypre built with GPU support");
5249 this->fespace = fespace_;
5259 int relax_coarse = 8;
5262 int interp_vec_variant = 2;
5264 int smooth_interp_vectors = 1;
5268 int interp_refine = interp_refine_;
5270 HYPRE_BoomerAMGSetNodal(amg_precond, nodal);
5271 HYPRE_BoomerAMGSetNodalDiag(amg_precond, nodal_diag);
5272 HYPRE_BoomerAMGSetCycleRelaxType(amg_precond, relax_coarse, 3);
5273 HYPRE_BoomerAMGSetInterpVecVariant(amg_precond, interp_vec_variant);
5274 HYPRE_BoomerAMGSetInterpVecQMax(amg_precond, q_max);
5275 HYPRE_BoomerAMGSetSmoothInterpVectors(amg_precond, smooth_interp_vectors);
5276 HYPRE_BoomerAMGSetInterpRefine(amg_precond, interp_refine);
5279 HYPRE_BoomerAMGSetInterpVectors(amg_precond, rbms.
Size(), rbms.
GetData());
5288#if MFEM_HYPRE_VERSION >= 21800
5291 const std::string &prerelax,
5292 const std::string &postrelax)
5296 int interp_type = 100;
5297 int relax_type = 10;
5298 int coarsen_type = 6;
5299 real_t strength_tolC = 0.1;
5300 real_t strength_tolR = 0.01;
5301 real_t filter_tolR = 0.0;
5302 real_t filterA_tol = 0.0;
5305 int ns_down = 0, ns_up = 0, ns_coarse;
5308 ns_down = prerelax.length();
5309 ns_up = postrelax.length();
5313 HYPRE_Int **grid_relax_points = mfem_hypre_TAlloc(HYPRE_Int*, 4);
5314 grid_relax_points[0] = NULL;
5315 grid_relax_points[1] = mfem_hypre_TAlloc(HYPRE_Int, ns_down);
5316 grid_relax_points[2] = mfem_hypre_TAlloc(HYPRE_Int, ns_up);
5317 grid_relax_points[3] = mfem_hypre_TAlloc(HYPRE_Int, 1);
5318 grid_relax_points[3][0] = 0;
5321 for (
int i = 0; i<ns_down; i++)
5323 if (prerelax[i] ==
'F')
5325 grid_relax_points[1][i] = -1;
5327 else if (prerelax[i] ==
'C')
5329 grid_relax_points[1][i] = 1;
5331 else if (prerelax[i] ==
'A')
5333 grid_relax_points[1][i] = 0;
5338 for (
int i = 0; i<ns_up; i++)
5340 if (postrelax[i] ==
'F')
5342 grid_relax_points[2][i] = -1;
5344 else if (postrelax[i] ==
'C')
5346 grid_relax_points[2][i] = 1;
5348 else if (postrelax[i] ==
'A')
5350 grid_relax_points[2][i] = 0;
5354 HYPRE_BoomerAMGSetRestriction(amg_precond, distanceR);
5356 HYPRE_BoomerAMGSetGridRelaxPoints(amg_precond, grid_relax_points);
5358 HYPRE_BoomerAMGSetInterpType(amg_precond, interp_type);
5363 HYPRE_BoomerAMGSetSabs(amg_precond, Sabs);
5366 HYPRE_BoomerAMGSetCoarsenType(amg_precond, coarsen_type);
5369 HYPRE_BoomerAMGSetAggNumLevels(amg_precond, 0);
5371 HYPRE_BoomerAMGSetStrongThreshold(amg_precond, strength_tolC);
5375 HYPRE_BoomerAMGSetStrongThresholdR(amg_precond, strength_tolR);
5376 HYPRE_BoomerAMGSetFilterThresholdR(amg_precond, filter_tolR);
5379 if (relax_type > -1)
5381 HYPRE_BoomerAMGSetRelaxType(amg_precond, relax_type);
5386 HYPRE_BoomerAMGSetCycleNumSweeps(amg_precond, ns_coarse, 3);
5387 HYPRE_BoomerAMGSetCycleNumSweeps(amg_precond, ns_down, 1);
5388 HYPRE_BoomerAMGSetCycleNumSweeps(amg_precond, ns_up, 2);
5390 HYPRE_BoomerAMGSetADropTol(amg_precond, filterA_tol);
5392 HYPRE_BoomerAMGSetADropType(amg_precond, -1);
5400 for (
int i = 0; i < rbms.
Size(); i++)
5402 HYPRE_ParVectorDestroy(rbms[i]);
5405 HYPRE_BoomerAMGDestroy(amg_precond);
5431 MFEM_ASSERT(G != NULL,
"");
5432 MFEM_ASSERT(x != NULL,
"");
5433 MFEM_ASSERT(y != NULL,
"");
5434 int sdim = (z == NULL) ? 2 : 3;
5435 int cycle_type = 13;
5436 MakeSolver(sdim, cycle_type);
5438 HYPRE_ParVector pz = z ?
static_cast<HYPRE_ParVector
>(*z) : NULL;
5439 HYPRE_AMSSetCoordinateVectors(ams, *x, *y, pz);
5440 HYPRE_AMSSetDiscreteGradient(ams, *G);
5448 int cycle_type = 13;
5453 trace_space = trace_space || rt_trace_space;
5459 "HypreAMS does not support variable order spaces");
5466 MakeSolver(std::max(sdim, vdim), cycle_type);
5467 MakeGradientAndInterpolation(edge_fespace, cycle_type);
5471 delete edge_fespace;
5476void HypreAMS::MakeSolver(
int sdim,
int cycle_type)
5482 int amg_coarsen_type = hypre_gpu ? 8 : 10;
5483 int amg_agg_levels = hypre_gpu ? 0 : 1;
5484 int amg_rlx_type = hypre_gpu ? 18 : 8;
5485 int rlx_type = hypre_gpu ? 1: 2;
5487 int amg_interp_type = 6;
5491 ams_cycle_type = cycle_type;
5492 HYPRE_AMSCreate(&ams);
5494 HYPRE_AMSSetDimension(ams, sdim);
5495 HYPRE_AMSSetTol(ams, 0.0);
5496 HYPRE_AMSSetMaxIter(ams, 1);
5497 HYPRE_AMSSetCycleType(ams, cycle_type);
5498 HYPRE_AMSSetPrintLevel(ams, 1);
5501 HYPRE_AMSSetSmoothingOptions(ams, rlx_type, rlx_sweeps, rlx_weight, rlx_omega);
5502 HYPRE_AMSSetAlphaAMGOptions(ams, amg_coarsen_type, amg_agg_levels, amg_rlx_type,
5503 theta, amg_interp_type, amg_Pmax);
5504 HYPRE_AMSSetBetaAMGOptions(ams, amg_coarsen_type, amg_agg_levels, amg_rlx_type,
5505 theta, amg_interp_type, amg_Pmax);
5507 HYPRE_AMSSetAlphaAMGCoarseRelaxType(ams, amg_rlx_type);
5508 HYPRE_AMSSetBetaAMGCoarseRelaxType(ams, amg_rlx_type);
5517void HypreAMS::MakeGradientAndInterpolation(
5518 ParFiniteElementSpace *edge_fespace,
int cycle_type)
5520 const FiniteElementCollection *edge_fec = edge_fespace->FEColl();
5521 bool trace_space =
dynamic_cast<const ND_Trace_FECollection *
>(edge_fec);
5523 ParMesh *pmesh = edge_fespace->GetParMesh();
5525 int sdim = pmesh->SpaceDimension();
5526 int vdim = edge_fespace->FEColl()->GetRangeDim(
dim - trace_space);
5529 MFEM_VERIFY(!edge_fespace->IsVariableOrder(),
5530 "HypreAMS does not support variable order spaces");
5531 int p = edge_fec->GetOrder() + (
dim - trace_space == 1 ? 1 : 0);
5534 FiniteElementCollection *vert_fec;
5537 vert_fec =
new H1_Trace_FECollection(
p,
dim);
5541 vert_fec =
new H1_FECollection(
p,
dim);
5543 ParFiniteElementSpace *vert_fespace =
new ParFiniteElementSpace(pmesh,
5547 ParDiscreteLinearOperator *grad;
5548 grad =
new ParDiscreteLinearOperator(vert_fespace, edge_fespace);
5551 grad->AddTraceFaceInterpolator(
new GradientInterpolator);
5555 grad->AddDomainInterpolator(
new GradientInterpolator);
5559 G = grad->ParallelAssemble();
5560 HYPRE_AMSSetDiscreteGradient(ams, *G);
5565 Pi = Pix = Piy = Piz = NULL;
5566 if (
p == 1 && pmesh->GetNodes() == NULL && vdim <= sdim)
5568 ParGridFunction x_coord(vert_fespace);
5569 ParGridFunction y_coord(vert_fespace);
5570 ParGridFunction z_coord(vert_fespace);
5572 for (
int i = 0; i < pmesh->GetNV(); i++)
5574 coord = pmesh -> GetVertex(i);
5575 x_coord(i) = coord[0];
5576 if (sdim >= 2) { y_coord(i) = coord[1]; }
5577 if (sdim == 3) { z_coord(i) = coord[2]; }
5579 x = x_coord.ParallelProject();
5586 y = y_coord.ParallelProject();
5591 z = z_coord.ParallelProject();
5595 HYPRE_AMSSetCoordinateVectors(ams,
5596 x ? (HYPRE_ParVector)(*x) : NULL,
5597 y ? (HYPRE_ParVector)(*y) : NULL,
5598 z ? (HYPRE_ParVector)(*z) : NULL);
5602 ParFiniteElementSpace *vert_fespace_d =
5603 new ParFiniteElementSpace(pmesh, vert_fec, std::max(sdim, vdim),
5606 ParDiscreteLinearOperator *id_ND;
5607 id_ND =
new ParDiscreteLinearOperator(vert_fespace_d, edge_fespace);
5610 id_ND->AddTraceFaceInterpolator(
new IdentityInterpolator);
5614 id_ND->AddDomainInterpolator(
new IdentityInterpolator);
5619 if (cycle_type < 10)
5621 Pi = id_ND->ParallelAssemble();
5625 Array2D<HypreParMatrix *> Pi_blocks;
5626 id_ND->GetParBlocks(Pi_blocks);
5627 Pix = Pi_blocks(0,0);
5628 if (std::max(sdim, vdim) >= 2) { Piy = Pi_blocks(0,1); }
5629 if (std::max(sdim, vdim) == 3) { Piz = Pi_blocks(0,2); }
5634 HYPRE_ParCSRMatrix HY_Pi = (Pi) ? (HYPRE_ParCSRMatrix) *Pi : NULL;
5635 HYPRE_ParCSRMatrix HY_Pix = (Pix) ? (HYPRE_ParCSRMatrix) *Pix : NULL;
5636 HYPRE_ParCSRMatrix HY_Piy = (Piy) ? (HYPRE_ParCSRMatrix) *Piy : NULL;
5637 HYPRE_ParCSRMatrix HY_Piz = (Piz) ? (HYPRE_ParCSRMatrix) *Piz : NULL;
5638 HYPRE_AMSSetInterpolations(ams, HY_Pi, HY_Pix, HY_Piy, HY_Piz);
5640 delete vert_fespace_d;
5643 delete vert_fespace;
5647void HypreAMS::ResetAMSPrecond()
5649#if MFEM_HYPRE_VERSION >= 22600
5651 auto *ams_data = (hypre_AMSData *)ams;
5654 HYPRE_Int
dim = hypre_AMSDataDimension(ams_data);
5657 hypre_ParCSRMatrix *hy_G = hypre_AMSDataDiscreteGradient(ams_data);
5659 HYPRE_Int beta_is_zero = hypre_AMSDataBetaIsZero(ams_data);
5662 hypre_ParCSRMatrix *hy_Pi hypre_AMSDataPiInterpolation(ams_data);
5663 hypre_ParCSRMatrix *hy_Pix = ams_data->Pix;
5664 hypre_ParCSRMatrix *hy_Piy = ams_data->Piy;
5665 hypre_ParCSRMatrix *hy_Piz = ams_data->Piz;
5666 HYPRE_Int owns_Pi = hypre_AMSDataOwnsPiInterpolation(ams_data);
5669 ams_data->owns_Pi = 0;
5673 hypre_ParVector *hy_x = hypre_AMSDataVertexCoordinateX(ams_data);
5674 hypre_ParVector *hy_y = hypre_AMSDataVertexCoordinateY(ams_data);
5675 hypre_ParVector *hy_z = hypre_AMSDataVertexCoordinateZ(ams_data);
5678 HYPRE_Int maxit = hypre_AMSDataMaxIter(ams_data);
5679 HYPRE_Real tol = hypre_AMSDataTol(ams_data);
5680 HYPRE_Int cycle_type = hypre_AMSDataCycleType(ams_data);
5681 HYPRE_Int ams_print_level = hypre_AMSDataPrintLevel(ams_data);
5684 HYPRE_Int A_relax_type = hypre_AMSDataARelaxType(ams_data);
5685 HYPRE_Int A_relax_times = hypre_AMSDataARelaxTimes(ams_data);
5686 HYPRE_Real A_relax_weight = hypre_AMSDataARelaxWeight(ams_data);
5687 HYPRE_Real A_omega = hypre_AMSDataAOmega(ams_data);
5688 HYPRE_Int A_cheby_order = hypre_AMSDataAChebyOrder(ams_data);
5689 HYPRE_Real A_cheby_fraction = hypre_AMSDataAChebyFraction(ams_data);
5691 HYPRE_Int B_Pi_coarsen_type = hypre_AMSDataPoissonAlphaAMGCoarsenType(ams_data);
5692 HYPRE_Int B_Pi_agg_levels = hypre_AMSDataPoissonAlphaAMGAggLevels(ams_data);
5693 HYPRE_Int B_Pi_relax_type = hypre_AMSDataPoissonAlphaAMGRelaxType(ams_data);
5694 HYPRE_Int B_Pi_coarse_relax_type = ams_data->B_Pi_coarse_relax_type;
5695 HYPRE_Real B_Pi_theta = hypre_AMSDataPoissonAlphaAMGStrengthThreshold(ams_data);
5696 HYPRE_Int B_Pi_interp_type = ams_data->B_Pi_interp_type;
5697 HYPRE_Int B_Pi_Pmax = ams_data->B_Pi_Pmax;
5699 HYPRE_Int B_G_coarsen_type = hypre_AMSDataPoissonBetaAMGCoarsenType(ams_data);
5700 HYPRE_Int B_G_agg_levels = hypre_AMSDataPoissonBetaAMGAggLevels(ams_data);
5701 HYPRE_Int B_G_relax_type = hypre_AMSDataPoissonBetaAMGRelaxType(ams_data);
5702 HYPRE_Int B_G_coarse_relax_type = ams_data->B_G_coarse_relax_type;
5703 HYPRE_Real B_G_theta = hypre_AMSDataPoissonBetaAMGStrengthThreshold(ams_data);
5704 HYPRE_Int B_G_interp_type = ams_data->B_G_interp_type;
5705 HYPRE_Int B_G_Pmax = ams_data->B_G_Pmax;
5707 HYPRE_AMSDestroy(ams);
5708 HYPRE_AMSCreate(&ams);
5709 ams_data = (hypre_AMSData *)ams;
5711 HYPRE_AMSSetDimension(ams,
dim);
5712 HYPRE_AMSSetTol(ams, tol);
5713 HYPRE_AMSSetMaxIter(ams, maxit);
5714 HYPRE_AMSSetCycleType(ams, cycle_type);
5715 HYPRE_AMSSetPrintLevel(ams, ams_print_level);
5717 HYPRE_AMSSetCoordinateVectors(ams, hy_x, hy_y, hy_z);
5719 HYPRE_AMSSetDiscreteGradient(ams, hy_G);
5720 HYPRE_AMSSetCoordinateVectors(ams, hy_x, hy_y, hy_z);
5721 HYPRE_AMSSetInterpolations(ams, hy_Pi, hy_Pix, hy_Piy, hy_Piz);
5722 ams_data->owns_Pi = owns_Pi;
5725 HYPRE_AMSSetSmoothingOptions(ams, A_relax_type, A_relax_times, A_relax_weight,
5728 hypre_AMSDataAChebyOrder(ams_data) = A_cheby_order;
5729 hypre_AMSDataAChebyFraction(ams_data) = A_cheby_fraction;
5731 HYPRE_AMSSetAlphaAMGOptions(ams, B_Pi_coarsen_type, B_Pi_agg_levels,
5733 B_Pi_theta, B_Pi_interp_type, B_Pi_Pmax);
5734 HYPRE_AMSSetBetaAMGOptions(ams, B_G_coarsen_type, B_G_agg_levels,
5736 B_G_theta, B_G_interp_type, B_G_Pmax);
5738 HYPRE_AMSSetAlphaAMGCoarseRelaxType(ams, B_Pi_coarse_relax_type);
5739 HYPRE_AMSSetBetaAMGCoarseRelaxType(ams, B_G_coarse_relax_type);
5741 ams_data->beta_is_zero = beta_is_zero;
5744 HYPRE_AMSDestroy(ams);
5746 MakeSolver(space_dim, ams_cycle_type);
5748 HYPRE_AMSSetPrintLevel(ams, print_level);
5749 if (singular) { HYPRE_AMSSetBetaPoissonMatrix(ams, NULL); }
5751 HYPRE_AMSSetDiscreteGradient(ams, *G);
5754 HYPRE_AMSSetCoordinateVectors(ams,
5755 x ? (HYPRE_ParVector)(*x) : nullptr,
5756 y ? (HYPRE_ParVector)(*y) : nullptr,
5757 z ? (HYPRE_ParVector)(*z) : nullptr);
5761 HYPRE_AMSSetInterpolations(ams,
5762 Pi ? (HYPRE_ParCSRMatrix) *Pi : nullptr,
5763 Pix ? (HYPRE_ParCSRMatrix) *Pix : nullptr,
5764 Piy ? (HYPRE_ParCSRMatrix) *Piy : nullptr,
5765 Piz ? (HYPRE_ParCSRMatrix) *Piz : nullptr);
5773 MFEM_VERIFY(new_A,
"new Operator must be a HypreParMatrix!");
5775 if (
A) { ResetAMSPrecond(); }
5792 HYPRE_AMSDestroy(ams);
5807 HYPRE_AMSSetPrintLevel(ams, print_lvl);
5808 print_level = print_lvl;
5826 x(x_), y(y_), z(z_),
5828 ND_Pi(NULL), ND_Pix(NULL), ND_Piy(NULL), ND_Piz(NULL),
5829 RT_Pi(NULL), RT_Pix(NULL), RT_Piy(NULL), RT_Piz(NULL)
5831 MFEM_ASSERT(C != NULL,
"");
5832 MFEM_ASSERT(G != NULL,
"");
5833 MFEM_ASSERT(x != NULL,
"");
5834 MFEM_ASSERT(y != NULL,
"");
5835 MFEM_ASSERT(z != NULL,
"");
5839 HYPRE_ADSSetCoordinateVectors(ads, *x, *y, *z);
5840 HYPRE_ADSSetDiscreteCurl(ads, *C);
5841 HYPRE_ADSSetDiscreteGradient(ads, *G);
5844void HypreADS::MakeSolver()
5850 int rlx_type = hypre_gpu ? 1 : 2;
5851 int amg_coarsen_type = hypre_gpu ? 8 : 10;
5852 int amg_agg_levels = hypre_gpu ? 0 : 1;
5853 int amg_rlx_type = hypre_gpu ? 18 : 8;
5855 int amg_interp_type = 6;
5858 HYPRE_ADSCreate(&ads);
5860 HYPRE_ADSSetTol(ads, 0.0);
5861 HYPRE_ADSSetMaxIter(ads, 1);
5862 HYPRE_ADSSetCycleType(ads, cycle_type);
5863 HYPRE_ADSSetPrintLevel(ads, 1);
5866 HYPRE_ADSSetSmoothingOptions(ads, rlx_type, rlx_sweeps, rlx_weight, rlx_omega);
5867 HYPRE_ADSSetAMGOptions(ads, amg_coarsen_type, amg_agg_levels, amg_rlx_type,
5868 theta, amg_interp_type, amg_Pmax);
5869 HYPRE_ADSSetAMSOptions(ads, ams_cycle_type, amg_coarsen_type, amg_agg_levels,
5870 amg_rlx_type, theta, amg_interp_type, amg_Pmax);
5879void HypreADS::MakeDiscreteMatrices(ParFiniteElementSpace *face_fespace)
5881 const FiniteElementCollection *face_fec = face_fespace->FEColl();