1 /* 2 Creates hypre ijmatrix from PETSc matrix 3 */ 4 5 #include <petscpkg_version.h> 6 #include <petsc/private/petschypre.h> 7 #include <petscmathypre.h> 8 #include <petsc/private/matimpl.h> 9 #include <petsc/private/deviceimpl.h> 10 #include <../src/mat/impls/hypre/mhypre.h> 11 #include <../src/mat/impls/aij/mpi/mpiaij.h> 12 #include <../src/vec/vec/impls/hypre/vhyp.h> 13 #include <HYPRE.h> 14 #include <HYPRE_utilities.h> 15 #include <_hypre_parcsr_ls.h> 16 #include <_hypre_sstruct_ls.h> 17 18 #if PETSC_PKG_HYPRE_VERSION_LT(2, 18, 0) 19 #define hypre_ParCSRMatrixClone(A, B) hypre_ParCSRMatrixCompleteClone(A) 20 #endif 21 22 static PetscErrorCode MatHYPRE_CreateFromMat(Mat, Mat_HYPRE *); 23 static PetscErrorCode MatHYPRE_IJMatrixPreallocate(Mat, Mat, HYPRE_IJMatrix); 24 static PetscErrorCode MatHYPRE_IJMatrixCopyIJ_MPIAIJ(Mat, HYPRE_IJMatrix); 25 static PetscErrorCode MatHYPRE_IJMatrixCopyIJ_SeqAIJ(Mat, HYPRE_IJMatrix); 26 static PetscErrorCode MatHYPRE_MultKernel_Private(Mat, HYPRE_Complex, Vec, HYPRE_Complex, Vec, PetscBool); 27 static PetscErrorCode MatSetValues_HYPRE(Mat, PetscInt, const PetscInt[], PetscInt, const PetscInt[], const PetscScalar[], InsertMode ins); 28 29 static PetscErrorCode MatHYPRE_IJMatrixPreallocate(Mat A_d, Mat A_o, HYPRE_IJMatrix ij) 30 { 31 PetscInt i, n_d, n_o; 32 const PetscInt *ia_d, *ia_o; 33 PetscBool done_d = PETSC_FALSE, done_o = PETSC_FALSE; 34 HYPRE_Int *nnz_d = NULL, *nnz_o = NULL; 35 36 PetscFunctionBegin; 37 if (A_d) { /* determine number of nonzero entries in local diagonal part */ 38 PetscCall(MatGetRowIJ(A_d, 0, PETSC_FALSE, PETSC_FALSE, &n_d, &ia_d, NULL, &done_d)); 39 if (done_d) { 40 PetscCall(PetscMalloc1(n_d, &nnz_d)); 41 for (i = 0; i < n_d; i++) nnz_d[i] = ia_d[i + 1] - ia_d[i]; 42 } 43 PetscCall(MatRestoreRowIJ(A_d, 0, PETSC_FALSE, PETSC_FALSE, NULL, &ia_d, NULL, &done_d)); 44 } 45 if (A_o) { /* determine number of nonzero entries in local off-diagonal part */ 46 PetscCall(MatGetRowIJ(A_o, 0, PETSC_FALSE, PETSC_FALSE, &n_o, &ia_o, NULL, &done_o)); 47 if (done_o) { 48 PetscCall(PetscMalloc1(n_o, &nnz_o)); 49 for (i = 0; i < n_o; i++) nnz_o[i] = ia_o[i + 1] - ia_o[i]; 50 } 51 PetscCall(MatRestoreRowIJ(A_o, 0, PETSC_FALSE, PETSC_FALSE, &n_o, &ia_o, NULL, &done_o)); 52 } 53 if (done_d) { /* set number of nonzeros in HYPRE IJ matrix */ 54 if (!done_o) { /* only diagonal part */ 55 PetscCall(PetscCalloc1(n_d, &nnz_o)); 56 } 57 #if PETSC_PKG_HYPRE_VERSION_GE(2, 16, 0) 58 { /* If we don't do this, the columns of the matrix will be all zeros! */ 59 hypre_AuxParCSRMatrix *aux_matrix; 60 aux_matrix = (hypre_AuxParCSRMatrix *)hypre_IJMatrixTranslator(ij); 61 hypre_AuxParCSRMatrixDestroy(aux_matrix); 62 hypre_IJMatrixTranslator(ij) = NULL; 63 PetscCallExternal(HYPRE_IJMatrixSetDiagOffdSizes, ij, nnz_d, nnz_o); 64 /* it seems they partially fixed it in 2.19.0 */ 65 #if PETSC_PKG_HYPRE_VERSION_LT(2, 19, 0) 66 aux_matrix = (hypre_AuxParCSRMatrix *)hypre_IJMatrixTranslator(ij); 67 hypre_AuxParCSRMatrixNeedAux(aux_matrix) = 1; 68 #endif 69 } 70 #else 71 PetscCallExternal(HYPRE_IJMatrixSetDiagOffdSizes, ij, nnz_d, nnz_o); 72 #endif 73 PetscCall(PetscFree(nnz_d)); 74 PetscCall(PetscFree(nnz_o)); 75 } 76 PetscFunctionReturn(PETSC_SUCCESS); 77 } 78 79 static PetscErrorCode MatHYPRE_CreateFromMat(Mat A, Mat_HYPRE *hA) 80 { 81 PetscInt rstart, rend, cstart, cend; 82 83 PetscFunctionBegin; 84 PetscCall(PetscLayoutSetUp(A->rmap)); 85 PetscCall(PetscLayoutSetUp(A->cmap)); 86 rstart = A->rmap->rstart; 87 rend = A->rmap->rend; 88 cstart = A->cmap->rstart; 89 cend = A->cmap->rend; 90 PetscHYPREInitialize(); 91 if (hA->ij) { 92 if (!hA->inner_free) hypre_IJMatrixObject(hA->ij) = NULL; 93 PetscCallExternal(HYPRE_IJMatrixDestroy, hA->ij); 94 } 95 PetscCallExternal(HYPRE_IJMatrixCreate, hA->comm, rstart, rend - 1, cstart, cend - 1, &hA->ij); 96 PetscCallExternal(HYPRE_IJMatrixSetObjectType, hA->ij, HYPRE_PARCSR); 97 { 98 PetscBool same; 99 Mat A_d, A_o; 100 const PetscInt *colmap; 101 PetscCall(PetscObjectBaseTypeCompare((PetscObject)A, MATMPIAIJ, &same)); 102 if (same) { 103 PetscCall(MatMPIAIJGetSeqAIJ(A, &A_d, &A_o, &colmap)); 104 PetscCall(MatHYPRE_IJMatrixPreallocate(A_d, A_o, hA->ij)); 105 PetscFunctionReturn(PETSC_SUCCESS); 106 } 107 PetscCall(PetscObjectBaseTypeCompare((PetscObject)A, MATMPIBAIJ, &same)); 108 if (same) { 109 PetscCall(MatMPIBAIJGetSeqBAIJ(A, &A_d, &A_o, &colmap)); 110 PetscCall(MatHYPRE_IJMatrixPreallocate(A_d, A_o, hA->ij)); 111 PetscFunctionReturn(PETSC_SUCCESS); 112 } 113 PetscCall(PetscObjectBaseTypeCompare((PetscObject)A, MATSEQAIJ, &same)); 114 if (same) { 115 PetscCall(MatHYPRE_IJMatrixPreallocate(A, NULL, hA->ij)); 116 PetscFunctionReturn(PETSC_SUCCESS); 117 } 118 PetscCall(PetscObjectBaseTypeCompare((PetscObject)A, MATSEQBAIJ, &same)); 119 if (same) { 120 PetscCall(MatHYPRE_IJMatrixPreallocate(A, NULL, hA->ij)); 121 PetscFunctionReturn(PETSC_SUCCESS); 122 } 123 } 124 PetscFunctionReturn(PETSC_SUCCESS); 125 } 126 127 static PetscErrorCode MatHYPRE_IJMatrixCopyIJ(Mat A, HYPRE_IJMatrix ij) 128 { 129 PetscBool flg; 130 131 PetscFunctionBegin; 132 #if PETSC_PKG_HYPRE_VERSION_LT(2, 19, 0) 133 PetscCallExternal(HYPRE_IJMatrixInitialize, ij); 134 #else 135 PetscCallExternal(HYPRE_IJMatrixInitialize_v2, ij, HYPRE_MEMORY_HOST); 136 #endif 137 PetscCall(PetscObjectBaseTypeCompare((PetscObject)A, MATMPIAIJ, &flg)); 138 if (flg) { 139 PetscCall(MatHYPRE_IJMatrixCopyIJ_MPIAIJ(A, ij)); 140 PetscFunctionReturn(PETSC_SUCCESS); 141 } 142 PetscCall(PetscObjectBaseTypeCompare((PetscObject)A, MATSEQAIJ, &flg)); 143 if (flg) { 144 PetscCall(MatHYPRE_IJMatrixCopyIJ_SeqAIJ(A, ij)); 145 PetscFunctionReturn(PETSC_SUCCESS); 146 } 147 PetscCheck(PETSC_FALSE, PetscObjectComm((PetscObject)A), PETSC_ERR_SUP, "No support for matrix type %s", ((PetscObject)A)->type_name); 148 } 149 150 static PetscErrorCode MatHYPRE_IJMatrixCopyIJ_SeqAIJ(Mat A, HYPRE_IJMatrix ij) 151 { 152 Mat_SeqAIJ *pdiag = (Mat_SeqAIJ *)A->data; 153 HYPRE_Int type; 154 hypre_ParCSRMatrix *par_matrix; 155 hypre_AuxParCSRMatrix *aux_matrix; 156 hypre_CSRMatrix *hdiag; 157 PetscBool sameint = (PetscBool)(sizeof(PetscInt) == sizeof(HYPRE_Int)); 158 159 PetscFunctionBegin; 160 PetscCallExternal(HYPRE_IJMatrixGetObjectType, ij, &type); 161 PetscCheck(type == HYPRE_PARCSR, PetscObjectComm((PetscObject)A), PETSC_ERR_SUP, "Only HYPRE_PARCSR is supported"); 162 PetscCallExternal(HYPRE_IJMatrixGetObject, ij, (void **)&par_matrix); 163 hdiag = hypre_ParCSRMatrixDiag(par_matrix); 164 /* 165 this is the Hack part where we monkey directly with the hypre datastructures 166 */ 167 if (sameint) { 168 PetscCall(PetscArraycpy(hdiag->i, pdiag->i, A->rmap->n + 1)); 169 PetscCall(PetscArraycpy(hdiag->j, pdiag->j, pdiag->nz)); 170 } else { 171 PetscInt i; 172 173 for (i = 0; i < A->rmap->n + 1; i++) hdiag->i[i] = (HYPRE_Int)pdiag->i[i]; 174 for (i = 0; i < pdiag->nz; i++) hdiag->j[i] = (HYPRE_Int)pdiag->j[i]; 175 } 176 177 aux_matrix = (hypre_AuxParCSRMatrix *)hypre_IJMatrixTranslator(ij); 178 hypre_AuxParCSRMatrixNeedAux(aux_matrix) = 0; 179 PetscFunctionReturn(PETSC_SUCCESS); 180 } 181 182 static PetscErrorCode MatHYPRE_IJMatrixCopyIJ_MPIAIJ(Mat A, HYPRE_IJMatrix ij) 183 { 184 Mat_MPIAIJ *pA = (Mat_MPIAIJ *)A->data; 185 Mat_SeqAIJ *pdiag, *poffd; 186 PetscInt i, *garray = pA->garray, *jj, cstart, *pjj; 187 HYPRE_Int *hjj, type; 188 hypre_ParCSRMatrix *par_matrix; 189 hypre_AuxParCSRMatrix *aux_matrix; 190 hypre_CSRMatrix *hdiag, *hoffd; 191 PetscBool sameint = (PetscBool)(sizeof(PetscInt) == sizeof(HYPRE_Int)); 192 193 PetscFunctionBegin; 194 pdiag = (Mat_SeqAIJ *)pA->A->data; 195 poffd = (Mat_SeqAIJ *)pA->B->data; 196 /* cstart is only valid for square MPIAIJ laid out in the usual way */ 197 PetscCall(MatGetOwnershipRange(A, &cstart, NULL)); 198 199 PetscCallExternal(HYPRE_IJMatrixGetObjectType, ij, &type); 200 PetscCheck(type == HYPRE_PARCSR, PetscObjectComm((PetscObject)A), PETSC_ERR_SUP, "Only HYPRE_PARCSR is supported"); 201 PetscCallExternal(HYPRE_IJMatrixGetObject, ij, (void **)&par_matrix); 202 hdiag = hypre_ParCSRMatrixDiag(par_matrix); 203 hoffd = hypre_ParCSRMatrixOffd(par_matrix); 204 205 if (sameint) { 206 PetscCall(PetscArraycpy(hdiag->i, pdiag->i, pA->A->rmap->n + 1)); 207 } else { 208 for (i = 0; i < pA->A->rmap->n + 1; i++) hdiag->i[i] = (HYPRE_Int)(pdiag->i[i]); 209 } 210 211 hjj = hdiag->j; 212 pjj = pdiag->j; 213 #if PETSC_PKG_HYPRE_VERSION_GE(2, 16, 0) 214 for (i = 0; i < pdiag->nz; i++) hjj[i] = pjj[i]; 215 #else 216 for (i = 0; i < pdiag->nz; i++) hjj[i] = cstart + pjj[i]; 217 #endif 218 if (sameint) { 219 PetscCall(PetscArraycpy(hoffd->i, poffd->i, pA->A->rmap->n + 1)); 220 } else { 221 for (i = 0; i < pA->A->rmap->n + 1; i++) hoffd->i[i] = (HYPRE_Int)(poffd->i[i]); 222 } 223 224 jj = (PetscInt *)hoffd->j; 225 #if PETSC_PKG_HYPRE_VERSION_GE(2, 16, 0) 226 PetscCallExternal(hypre_CSRMatrixBigInitialize, hoffd); 227 jj = (PetscInt *)hoffd->big_j; 228 #endif 229 pjj = poffd->j; 230 for (i = 0; i < poffd->nz; i++) jj[i] = garray[pjj[i]]; 231 232 aux_matrix = (hypre_AuxParCSRMatrix *)hypre_IJMatrixTranslator(ij); 233 hypre_AuxParCSRMatrixNeedAux(aux_matrix) = 0; 234 PetscFunctionReturn(PETSC_SUCCESS); 235 } 236 237 static PetscErrorCode MatConvert_HYPRE_IS(Mat A, MatType mtype, MatReuse reuse, Mat *B) 238 { 239 Mat_HYPRE *mhA = (Mat_HYPRE *)(A->data); 240 Mat lA; 241 ISLocalToGlobalMapping rl2g, cl2g; 242 IS is; 243 hypre_ParCSRMatrix *hA; 244 hypre_CSRMatrix *hdiag, *hoffd; 245 MPI_Comm comm; 246 HYPRE_Complex *hdd, *hod, *aa; 247 PetscScalar *data; 248 HYPRE_BigInt *col_map_offd; 249 HYPRE_Int *hdi, *hdj, *hoi, *hoj; 250 PetscInt *ii, *jj, *iptr, *jptr; 251 PetscInt cum, dr, dc, oc, str, stc, nnz, i, jd, jo, M, N; 252 HYPRE_Int type; 253 MatType lmattype = NULL; 254 PetscBool freeparcsr = PETSC_FALSE; 255 256 PetscFunctionBegin; 257 comm = PetscObjectComm((PetscObject)A); 258 PetscCallExternal(HYPRE_IJMatrixGetObjectType, mhA->ij, &type); 259 PetscCheck(type == HYPRE_PARCSR, comm, PETSC_ERR_SUP, "Only HYPRE_PARCSR is supported"); 260 PetscCallExternal(HYPRE_IJMatrixGetObject, mhA->ij, (void **)&hA); 261 #if defined(PETSC_HAVE_HYPRE_DEVICE) 262 if (HYPRE_MEMORY_DEVICE == hypre_IJMatrixMemoryLocation(mhA->ij)) { 263 /* Support by copying back on the host and copy to GPU 264 Kind of inefficient, but this is the best we can do now */ 265 #if defined(HYPRE_USING_HIP) 266 lmattype = MATSEQAIJHIPSPARSE; 267 #elif defined(HYPRE_USING_CUDA) 268 lmattype = MATSEQAIJCUSPARSE; 269 #endif 270 hA = hypre_ParCSRMatrixClone_v2(hA, 1, HYPRE_MEMORY_HOST); 271 freeparcsr = PETSC_TRUE; 272 } 273 #endif 274 M = hypre_ParCSRMatrixGlobalNumRows(hA); 275 N = hypre_ParCSRMatrixGlobalNumCols(hA); 276 str = hypre_ParCSRMatrixFirstRowIndex(hA); 277 stc = hypre_ParCSRMatrixFirstColDiag(hA); 278 hdiag = hypre_ParCSRMatrixDiag(hA); 279 hoffd = hypre_ParCSRMatrixOffd(hA); 280 dr = hypre_CSRMatrixNumRows(hdiag); 281 dc = hypre_CSRMatrixNumCols(hdiag); 282 nnz = hypre_CSRMatrixNumNonzeros(hdiag); 283 hdi = hypre_CSRMatrixI(hdiag); 284 hdj = hypre_CSRMatrixJ(hdiag); 285 hdd = hypre_CSRMatrixData(hdiag); 286 oc = hypre_CSRMatrixNumCols(hoffd); 287 nnz += hypre_CSRMatrixNumNonzeros(hoffd); 288 hoi = hypre_CSRMatrixI(hoffd); 289 hoj = hypre_CSRMatrixJ(hoffd); 290 hod = hypre_CSRMatrixData(hoffd); 291 if (reuse != MAT_REUSE_MATRIX) { 292 PetscInt *aux; 293 294 /* generate l2g maps for rows and cols */ 295 PetscCall(ISCreateStride(comm, dr, str, 1, &is)); 296 PetscCall(ISLocalToGlobalMappingCreateIS(is, &rl2g)); 297 PetscCall(ISDestroy(&is)); 298 col_map_offd = hypre_ParCSRMatrixColMapOffd(hA); 299 PetscCall(PetscMalloc1(dc + oc, &aux)); 300 for (i = 0; i < dc; i++) aux[i] = i + stc; 301 for (i = 0; i < oc; i++) aux[i + dc] = col_map_offd[i]; 302 PetscCall(ISCreateGeneral(comm, dc + oc, aux, PETSC_OWN_POINTER, &is)); 303 PetscCall(ISLocalToGlobalMappingCreateIS(is, &cl2g)); 304 PetscCall(ISDestroy(&is)); 305 /* create MATIS object */ 306 PetscCall(MatCreate(comm, B)); 307 PetscCall(MatSetSizes(*B, dr, dc, M, N)); 308 PetscCall(MatSetType(*B, MATIS)); 309 PetscCall(MatSetLocalToGlobalMapping(*B, rl2g, cl2g)); 310 PetscCall(ISLocalToGlobalMappingDestroy(&rl2g)); 311 PetscCall(ISLocalToGlobalMappingDestroy(&cl2g)); 312 313 /* allocate CSR for local matrix */ 314 PetscCall(PetscMalloc1(dr + 1, &iptr)); 315 PetscCall(PetscMalloc1(nnz, &jptr)); 316 PetscCall(PetscMalloc1(nnz, &data)); 317 } else { 318 PetscInt nr; 319 PetscBool done; 320 PetscCall(MatISGetLocalMat(*B, &lA)); 321 PetscCall(MatGetRowIJ(lA, 0, PETSC_FALSE, PETSC_FALSE, &nr, (const PetscInt **)&iptr, (const PetscInt **)&jptr, &done)); 322 PetscCheck(nr == dr, PETSC_COMM_SELF, PETSC_ERR_USER, "Cannot reuse mat: invalid number of rows in local mat! %" PetscInt_FMT " != %" PetscInt_FMT, nr, dr); 323 PetscCheck(iptr[nr] >= nnz, PETSC_COMM_SELF, PETSC_ERR_USER, "Cannot reuse mat: invalid number of nonzeros in local mat! reuse %" PetscInt_FMT " requested %" PetscInt_FMT, iptr[nr], nnz); 324 PetscCall(MatSeqAIJGetArrayWrite(lA, &data)); 325 } 326 /* merge local matrices */ 327 ii = iptr; 328 jj = jptr; 329 aa = (HYPRE_Complex *)data; /* this cast fixes the clang error when doing the assignments below: implicit conversion from 'HYPRE_Complex' (aka '_Complex double') to 'double' is not permitted in C++ */ 330 *ii = *(hdi++) + *(hoi++); 331 for (jd = 0, jo = 0, cum = 0; *ii < nnz; cum++) { 332 PetscScalar *aold = (PetscScalar *)aa; 333 PetscInt *jold = jj, nc = jd + jo; 334 for (; jd < *hdi; jd++) { 335 *jj++ = *hdj++; 336 *aa++ = *hdd++; 337 } 338 for (; jo < *hoi; jo++) { 339 *jj++ = *hoj++ + dc; 340 *aa++ = *hod++; 341 } 342 *(++ii) = *(hdi++) + *(hoi++); 343 PetscCall(PetscSortIntWithScalarArray(jd + jo - nc, jold, aold)); 344 } 345 for (; cum < dr; cum++) *(++ii) = nnz; 346 if (reuse != MAT_REUSE_MATRIX) { 347 Mat_SeqAIJ *a; 348 349 PetscCall(MatCreateSeqAIJWithArrays(PETSC_COMM_SELF, dr, dc + oc, iptr, jptr, data, &lA)); 350 /* hack SeqAIJ */ 351 a = (Mat_SeqAIJ *)(lA->data); 352 a->free_a = PETSC_TRUE; 353 a->free_ij = PETSC_TRUE; 354 if (lmattype) PetscCall(MatConvert(lA, lmattype, MAT_INPLACE_MATRIX, &lA)); 355 PetscCall(MatISSetLocalMat(*B, lA)); 356 PetscCall(MatDestroy(&lA)); 357 } else { 358 PetscCall(MatSeqAIJRestoreArrayWrite(lA, &data)); 359 } 360 PetscCall(MatAssemblyBegin(*B, MAT_FINAL_ASSEMBLY)); 361 PetscCall(MatAssemblyEnd(*B, MAT_FINAL_ASSEMBLY)); 362 if (reuse == MAT_INPLACE_MATRIX) PetscCall(MatHeaderReplace(A, B)); 363 if (freeparcsr) PetscCallExternal(hypre_ParCSRMatrixDestroy, hA); 364 PetscFunctionReturn(PETSC_SUCCESS); 365 } 366 367 static PetscErrorCode MatHYPRE_DestroyCOOMat(Mat mat) 368 { 369 Mat_HYPRE *hA = (Mat_HYPRE *)mat->data; 370 371 PetscFunctionBegin; 372 if (hA->cooMat) { /* If cooMat is present we need to destroy the column indices */ 373 PetscCall(MatDestroy(&hA->cooMat)); 374 if (hA->cooMatAttached) { 375 hypre_CSRMatrix *csr; 376 hypre_ParCSRMatrix *parcsr; 377 HYPRE_MemoryLocation mem; 378 379 PetscCallExternal(HYPRE_IJMatrixGetObject, hA->ij, (void **)&parcsr); 380 csr = hypre_ParCSRMatrixDiag(parcsr); 381 if (csr) { 382 mem = hypre_CSRMatrixMemoryLocation(csr); 383 PetscStackCallExternalVoid("hypre_TFree", hypre_TFree(hypre_CSRMatrixJ(csr), mem)); 384 PetscStackCallExternalVoid("hypre_TFree", hypre_TFree(hypre_CSRMatrixBigJ(csr), mem)); 385 } 386 csr = hypre_ParCSRMatrixOffd(parcsr); 387 if (csr) { 388 mem = hypre_CSRMatrixMemoryLocation(csr); 389 PetscStackCallExternalVoid("hypre_TFree", hypre_TFree(hypre_CSRMatrixJ(csr), mem)); 390 PetscStackCallExternalVoid("hypre_TFree", hypre_TFree(hypre_CSRMatrixBigJ(csr), mem)); 391 } 392 } 393 } 394 hA->cooMatAttached = PETSC_FALSE; 395 PetscFunctionReturn(PETSC_SUCCESS); 396 } 397 398 static PetscErrorCode MatHYPRE_CreateCOOMat(Mat mat) 399 { 400 MPI_Comm comm; 401 PetscMPIInt size; 402 PetscLayout rmap, cmap; 403 Mat_HYPRE *hmat = (Mat_HYPRE *)mat->data; 404 MatType matType = MATAIJ; /* default type of cooMat */ 405 406 PetscFunctionBegin; 407 /* Build an agent matrix cooMat with AIJ format 408 It has the same sparsity pattern as mat, and also shares the data array with mat. We use cooMat to do the COO work. 409 */ 410 PetscCall(PetscObjectGetComm((PetscObject)mat, &comm)); 411 PetscCallMPI(MPI_Comm_size(comm, &size)); 412 PetscCall(PetscLayoutSetUp(mat->rmap)); 413 PetscCall(PetscLayoutSetUp(mat->cmap)); 414 PetscCall(MatGetLayouts(mat, &rmap, &cmap)); 415 416 #if defined(PETSC_HAVE_HYPRE_DEVICE) 417 if (!mat->boundtocpu) { /* mat will be on device, so will cooMat */ 418 #if defined(HYPRE_USING_HIP) 419 matType = MATAIJHIPSPARSE; 420 #elif defined(HYPRE_USING_CUDA) 421 matType = MATAIJCUSPARSE; 422 #else 423 SETERRQ(comm, PETSC_ERR_SUP, "Do not know the HYPRE device"); 424 #endif 425 } 426 #endif 427 428 /* Do COO preallocation through cooMat */ 429 PetscCall(MatHYPRE_DestroyCOOMat(mat)); 430 PetscCall(MatCreate(comm, &hmat->cooMat)); 431 PetscCall(MatSetType(hmat->cooMat, matType)); 432 PetscCall(MatSetLayouts(hmat->cooMat, rmap, cmap)); 433 434 /* allocate local matrices if needed */ 435 PetscCall(MatMPIAIJSetPreallocation(hmat->cooMat, 0, NULL, 0, NULL)); 436 PetscFunctionReturn(PETSC_SUCCESS); 437 } 438 439 /* Attach cooMat data array to hypre matrix. 440 When AIJCUPMSPARSE will support raw device pointers and not THRUSTARRAY 441 we should swap the arrays: i.e., attach hypre matrix array to cooMat 442 This is because hypre should be in charge of handling the memory, 443 cooMat is only a way to reuse PETSc COO code. 444 attaching the memory will then be done at MatSetValuesCOO time and it will dynamically 445 support hypre matrix migrating to host. 446 */ 447 static PetscErrorCode MatHYPRE_AttachCOOMat(Mat mat) 448 { 449 Mat_HYPRE *hmat = (Mat_HYPRE *)mat->data; 450 hypre_CSRMatrix *diag, *offd; 451 hypre_ParCSRMatrix *parCSR; 452 HYPRE_MemoryLocation hmem = HYPRE_MEMORY_HOST; 453 PetscMemType pmem; 454 Mat A, B; 455 PetscScalar *a; 456 PetscMPIInt size; 457 MPI_Comm comm; 458 459 PetscFunctionBegin; 460 PetscCheck(hmat->cooMat, PetscObjectComm((PetscObject)mat), PETSC_ERR_PLIB, "HYPRE COO delegate matrix has not been created yet"); 461 if (hmat->cooMatAttached) PetscFunctionReturn(PETSC_SUCCESS); 462 PetscCheck(hmat->cooMat->preallocated, PetscObjectComm((PetscObject)mat), PETSC_ERR_PLIB, "HYPRE COO delegate matrix is not preallocated"); 463 PetscCall(PetscObjectSetName((PetscObject)hmat->cooMat, "_internal_COO_mat_for_hypre")); 464 PetscCall(PetscObjectGetComm((PetscObject)mat, &comm)); 465 PetscCallMPI(MPI_Comm_size(comm, &size)); 466 467 /* Alias cooMat's data array to IJMatrix's */ 468 PetscCallExternal(HYPRE_IJMatrixGetObject, hmat->ij, (void **)&parCSR); 469 diag = hypre_ParCSRMatrixDiag(parCSR); 470 offd = hypre_ParCSRMatrixOffd(parCSR); 471 472 A = (size == 1) ? hmat->cooMat : ((Mat_MPIAIJ *)hmat->cooMat->data)->A; 473 B = (size == 1) ? NULL : ((Mat_MPIAIJ *)hmat->cooMat->data)->B; 474 475 PetscCall(PetscObjectSetName((PetscObject)A, "_internal_COO_mat_for_hypre")); 476 hmem = hypre_CSRMatrixMemoryLocation(diag); 477 PetscCall(MatSeqAIJGetCSRAndMemType(A, NULL, NULL, &a, &pmem)); 478 PetscAssert((PetscMemTypeHost(pmem) && hmem == HYPRE_MEMORY_HOST) || (PetscMemTypeDevice(pmem) && hmem == HYPRE_MEMORY_DEVICE), comm, PETSC_ERR_PLIB, "PETSc and hypre's memory types mismatch"); 479 PetscStackCallExternalVoid("hypre_TFree", hypre_TFree(hypre_CSRMatrixData(diag), hmem)); 480 hypre_CSRMatrixData(diag) = (HYPRE_Complex *)a; 481 hypre_CSRMatrixOwnsData(diag) = 0; /* Take ownership of (j,a) away from hypre. As a result, we need to free them on our own */ 482 483 if (B) { 484 hmem = hypre_CSRMatrixMemoryLocation(offd); 485 PetscCall(MatSeqAIJGetCSRAndMemType(B, NULL, NULL, &a, &pmem)); 486 PetscAssert((PetscMemTypeHost(pmem) && hmem == HYPRE_MEMORY_HOST) || (PetscMemTypeDevice(pmem) && hmem == HYPRE_MEMORY_DEVICE), comm, PETSC_ERR_PLIB, "PETSc and hypre's memory types mismatch"); 487 PetscStackCallExternalVoid("hypre_TFree", hypre_TFree(hypre_CSRMatrixData(offd), hmem)); 488 hypre_CSRMatrixData(offd) = (HYPRE_Complex *)a; 489 hypre_CSRMatrixOwnsData(offd) = 0; 490 } 491 hmat->cooMatAttached = PETSC_TRUE; 492 PetscFunctionReturn(PETSC_SUCCESS); 493 } 494 495 static PetscErrorCode CSRtoCOO_Private(PetscInt n, const PetscInt ii[], const PetscInt jj[], PetscCount *ncoo, PetscInt **coo_i, PetscInt **coo_j) 496 { 497 PetscInt *cooi, *cooj; 498 499 PetscFunctionBegin; 500 *ncoo = ii[n]; 501 PetscCall(PetscMalloc2(*ncoo, &cooi, *ncoo, &cooj)); 502 for (PetscInt i = 0; i < n; i++) { 503 for (PetscInt j = ii[i]; j < ii[i + 1]; j++) cooi[j] = i; 504 } 505 PetscCall(PetscArraycpy(cooj, jj, *ncoo)); 506 *coo_i = cooi; 507 *coo_j = cooj; 508 PetscFunctionReturn(PETSC_SUCCESS); 509 } 510 511 static PetscErrorCode CSRtoCOO_HYPRE_Int_Private(PetscInt n, const HYPRE_Int ii[], const HYPRE_Int jj[], PetscCount *ncoo, PetscInt **coo_i, PetscInt **coo_j) 512 { 513 PetscInt *cooi, *cooj; 514 515 PetscFunctionBegin; 516 *ncoo = ii[n]; 517 PetscCall(PetscMalloc2(*ncoo, &cooi, *ncoo, &cooj)); 518 for (PetscInt i = 0; i < n; i++) { 519 for (HYPRE_Int j = ii[i]; j < ii[i + 1]; j++) cooi[j] = i; 520 } 521 for (PetscCount i = 0; i < *ncoo; i++) cooj[i] = jj[i]; 522 *coo_i = cooi; 523 *coo_j = cooj; 524 PetscFunctionReturn(PETSC_SUCCESS); 525 } 526 527 static PetscErrorCode MatSeqAIJGetCOO_Private(Mat A, PetscCount *ncoo, PetscInt **coo_i, PetscInt **coo_j) 528 { 529 PetscInt n; 530 const PetscInt *ii, *jj; 531 PetscBool done; 532 533 PetscFunctionBegin; 534 PetscCall(MatGetRowIJ(A, 0, PETSC_FALSE, PETSC_FALSE, &n, &ii, &jj, &done)); 535 PetscCheck(done, PetscObjectComm((PetscObject)A), PETSC_ERR_PLIB, "Failure for MatGetRowIJ"); 536 PetscCall(CSRtoCOO_Private(n, ii, jj, ncoo, coo_i, coo_j)); 537 PetscCall(MatRestoreRowIJ(A, 0, PETSC_FALSE, PETSC_FALSE, &n, &ii, &jj, &done)); 538 PetscCheck(done, PetscObjectComm((PetscObject)A), PETSC_ERR_PLIB, "Failure for MatRestoreRowIJ"); 539 PetscFunctionReturn(PETSC_SUCCESS); 540 } 541 542 static PetscErrorCode hypreCSRMatrixGetCOO_Private(hypre_CSRMatrix *A, PetscCount *ncoo, PetscInt **coo_i, PetscInt **coo_j) 543 { 544 PetscInt n = hypre_CSRMatrixNumRows(A); 545 HYPRE_Int *ii, *jj; 546 HYPRE_MemoryLocation mem = HYPRE_MEMORY_HOST; 547 548 PetscFunctionBegin; 549 #if defined(PETSC_HAVE_HYPRE_DEVICE) 550 mem = hypre_CSRMatrixMemoryLocation(A); 551 if (mem != HYPRE_MEMORY_HOST) { 552 PetscCount nnz = hypre_CSRMatrixNumNonzeros(A); 553 PetscCall(PetscMalloc2(n + 1, &ii, nnz, &jj)); 554 hypre_TMemcpy(ii, hypre_CSRMatrixI(A), HYPRE_Int, n + 1, HYPRE_MEMORY_HOST, mem); 555 hypre_TMemcpy(jj, hypre_CSRMatrixJ(A), HYPRE_Int, nnz, HYPRE_MEMORY_HOST, mem); 556 } else { 557 #else 558 { 559 #endif 560 ii = hypre_CSRMatrixI(A); 561 jj = hypre_CSRMatrixJ(A); 562 } 563 PetscCall(CSRtoCOO_HYPRE_Int_Private(n, ii, jj, ncoo, coo_i, coo_j)); 564 if (mem != HYPRE_MEMORY_HOST) PetscCall(PetscFree2(ii, jj)); 565 PetscFunctionReturn(PETSC_SUCCESS); 566 } 567 568 static PetscErrorCode MatSetValuesCOOFromCSRMatrix_Private(Mat A, hypre_CSRMatrix *H) 569 { 570 PetscBool iscpu = PETSC_TRUE; 571 PetscScalar *a; 572 HYPRE_MemoryLocation mem = HYPRE_MEMORY_HOST; 573 574 PetscFunctionBegin; 575 #if defined(PETSC_HAVE_HYPRE_DEVICE) 576 mem = hypre_CSRMatrixMemoryLocation(H); 577 PetscCall(PetscObjectTypeCompare((PetscObject)A, MATSEQAIJ, &iscpu)); 578 #endif 579 if (iscpu && mem != HYPRE_MEMORY_HOST) { 580 PetscCount nnz = hypre_CSRMatrixNumNonzeros(H); 581 PetscCall(PetscMalloc1(nnz, &a)); 582 hypre_TMemcpy(a, hypre_CSRMatrixData(H), PetscScalar, nnz, HYPRE_MEMORY_HOST, mem); 583 } else { 584 a = (PetscScalar *)hypre_CSRMatrixData(H); 585 } 586 PetscCall(MatSetValuesCOO(A, a, INSERT_VALUES)); 587 if (iscpu && mem != HYPRE_MEMORY_HOST) PetscCall(PetscFree(a)); 588 PetscFunctionReturn(PETSC_SUCCESS); 589 } 590 591 PETSC_INTERN PetscErrorCode MatConvert_AIJ_HYPRE(Mat A, MatType type, MatReuse reuse, Mat *B) 592 { 593 MPI_Comm comm = PetscObjectComm((PetscObject)A); 594 Mat M = NULL, dH = NULL, oH = NULL, dA = NULL, oA = NULL; 595 PetscBool ismpiaij, issbaij, isbaij; 596 Mat_HYPRE *hA; 597 598 PetscFunctionBegin; 599 PetscCall(PetscObjectTypeCompareAny((PetscObject)A, &issbaij, MATSEQSBAIJ, MATMPIBAIJ, "")); 600 PetscCall(PetscObjectTypeCompareAny((PetscObject)A, &isbaij, MATSEQBAIJ, MATMPIBAIJ, "")); 601 if (isbaij || issbaij) { /* handle BAIJ and SBAIJ */ 602 PetscBool ismpi; 603 MatType newtype; 604 605 PetscCall(PetscObjectTypeCompareAny((PetscObject)A, &ismpi, MATMPISBAIJ, MATMPIBAIJ, "")); 606 newtype = ismpi ? MATMPIAIJ : MATSEQAIJ; 607 if (reuse == MAT_REUSE_MATRIX) { 608 PetscCall(MatConvert(*B, newtype, MAT_INPLACE_MATRIX, B)); 609 PetscCall(MatConvert(A, newtype, MAT_REUSE_MATRIX, B)); 610 PetscCall(MatConvert(*B, MATHYPRE, MAT_INPLACE_MATRIX, B)); 611 } else if (reuse == MAT_INITIAL_MATRIX) { 612 PetscCall(MatConvert(A, newtype, MAT_INITIAL_MATRIX, B)); 613 PetscCall(MatConvert(*B, MATHYPRE, MAT_INPLACE_MATRIX, B)); 614 } else { 615 PetscCall(MatConvert(A, newtype, MAT_INPLACE_MATRIX, &A)); 616 PetscCall(MatConvert(A, MATHYPRE, MAT_INPLACE_MATRIX, &A)); 617 } 618 PetscFunctionReturn(PETSC_SUCCESS); 619 } 620 621 dA = A; 622 PetscCall(PetscObjectBaseTypeCompare((PetscObject)A, MATMPIAIJ, &ismpiaij)); 623 if (ismpiaij) PetscCall(MatMPIAIJGetSeqAIJ(A, &dA, &oA, NULL)); 624 625 if (reuse != MAT_REUSE_MATRIX) { 626 PetscCount coo_n; 627 PetscInt *coo_i, *coo_j; 628 629 PetscCall(MatCreate(comm, &M)); 630 PetscCall(MatSetType(M, MATHYPRE)); 631 PetscCall(MatSetSizes(M, A->rmap->n, A->cmap->n, A->rmap->N, A->cmap->N)); 632 PetscCall(MatSetOption(M, MAT_SORTED_FULL, PETSC_TRUE)); 633 PetscCall(MatSetOption(M, MAT_NO_OFF_PROC_ENTRIES, PETSC_TRUE)); 634 635 hA = (Mat_HYPRE *)M->data; 636 PetscCall(MatHYPRE_CreateFromMat(A, hA)); 637 PetscCall(MatHYPRE_IJMatrixCopyIJ(A, hA->ij)); 638 639 PetscCall(MatHYPRE_CreateCOOMat(M)); 640 641 dH = hA->cooMat; 642 PetscCall(PetscObjectBaseTypeCompare((PetscObject)hA->cooMat, MATMPIAIJ, &ismpiaij)); 643 if (ismpiaij) PetscCall(MatMPIAIJGetSeqAIJ(hA->cooMat, &dH, &oH, NULL)); 644 645 PetscCall(PetscObjectSetName((PetscObject)dH, "_internal_COO_mat_for_hypre")); 646 PetscCall(MatSeqAIJGetCOO_Private(dA, &coo_n, &coo_i, &coo_j)); 647 PetscCall(MatSetPreallocationCOO(dH, coo_n, coo_i, coo_j)); 648 PetscCall(PetscFree2(coo_i, coo_j)); 649 if (oH) { 650 PetscCall(PetscLayoutDestroy(&oH->cmap)); 651 PetscCall(PetscLayoutCreateFromSizes(PetscObjectComm((PetscObject)oH), oA->cmap->n, oA->cmap->n, 1, &oH->cmap)); 652 PetscCall(MatSeqAIJGetCOO_Private(oA, &coo_n, &coo_i, &coo_j)); 653 PetscCall(MatSetPreallocationCOO(oH, coo_n, coo_i, coo_j)); 654 PetscCall(PetscFree2(coo_i, coo_j)); 655 } 656 hA->cooMat->assembled = PETSC_TRUE; 657 658 M->preallocated = PETSC_TRUE; 659 PetscCall(MatAssemblyBegin(M, MAT_FINAL_ASSEMBLY)); 660 PetscCall(MatAssemblyEnd(M, MAT_FINAL_ASSEMBLY)); 661 662 PetscCall(MatHYPRE_AttachCOOMat(M)); 663 if (reuse == MAT_INITIAL_MATRIX) *B = M; 664 } else M = *B; 665 666 hA = (Mat_HYPRE *)M->data; 667 PetscCheck(hA->cooMat, PetscObjectComm((PetscObject)A), PETSC_ERR_PLIB, "HYPRE COO delegate matrix has not been created yet"); 668 669 dH = hA->cooMat; 670 PetscCall(PetscObjectBaseTypeCompare((PetscObject)hA->cooMat, MATMPIAIJ, &ismpiaij)); 671 if (ismpiaij) PetscCall(MatMPIAIJGetSeqAIJ(hA->cooMat, &dH, &oH, NULL)); 672 673 PetscScalar *a; 674 PetscCall(MatSeqAIJGetCSRAndMemType(dA, NULL, NULL, &a, NULL)); 675 PetscCall(MatSetValuesCOO(dH, a, INSERT_VALUES)); 676 if (oH) { 677 PetscCall(MatSeqAIJGetCSRAndMemType(oA, NULL, NULL, &a, NULL)); 678 PetscCall(MatSetValuesCOO(oH, a, INSERT_VALUES)); 679 } 680 681 if (reuse == MAT_INPLACE_MATRIX) PetscCall(MatHeaderReplace(A, &M)); 682 PetscFunctionReturn(PETSC_SUCCESS); 683 } 684 685 static PetscErrorCode MatConvert_HYPRE_AIJ(Mat A, MatType mtype, MatReuse reuse, Mat *B) 686 { 687 Mat M, dA = NULL, oA = NULL; 688 hypre_ParCSRMatrix *parcsr; 689 hypre_CSRMatrix *dH, *oH; 690 MPI_Comm comm; 691 PetscBool ismpiaij, isseqaij; 692 693 PetscFunctionBegin; 694 comm = PetscObjectComm((PetscObject)A); 695 if (reuse == MAT_REUSE_MATRIX) { 696 PetscCall(PetscObjectBaseTypeCompare((PetscObject)*B, MATMPIAIJ, &ismpiaij)); 697 PetscCall(PetscObjectBaseTypeCompare((PetscObject)*B, MATSEQAIJ, &isseqaij)); 698 PetscCheck(ismpiaij || isseqaij, comm, PETSC_ERR_SUP, "Only MATMPIAIJ or MATSEQAIJ base types are supported"); 699 } 700 PetscCall(MatHYPREGetParCSR(A, &parcsr)); 701 #if defined(PETSC_HAVE_HYPRE_DEVICE) 702 if (HYPRE_MEMORY_DEVICE == hypre_ParCSRMatrixMemoryLocation(parcsr)) { 703 PetscBool isaij; 704 705 PetscCall(PetscStrcmp(mtype, MATAIJ, &isaij)); 706 if (isaij) { 707 PetscMPIInt size; 708 709 PetscCallMPI(MPI_Comm_size(comm, &size)); 710 #if defined(HYPRE_USING_HIP) 711 mtype = size > 1 ? MATMPIAIJHIPSPARSE : MATSEQAIJHIPSPARSE; 712 #elif defined(HYPRE_USING_CUDA) 713 mtype = size > 1 ? MATMPIAIJCUSPARSE : MATSEQAIJCUSPARSE; 714 #else 715 mtype = size > 1 ? MATMPIAIJ : MATSEQAIJ; 716 #endif 717 } 718 } 719 #endif 720 dH = hypre_ParCSRMatrixDiag(parcsr); 721 oH = hypre_ParCSRMatrixOffd(parcsr); 722 if (reuse != MAT_REUSE_MATRIX) { 723 PetscCount coo_n; 724 PetscInt *coo_i, *coo_j; 725 726 PetscCall(MatCreate(comm, &M)); 727 PetscCall(MatSetType(M, mtype)); 728 PetscCall(MatSetSizes(M, A->rmap->n, A->cmap->n, A->rmap->N, A->cmap->N)); 729 PetscCall(MatMPIAIJSetPreallocation(M, 0, NULL, 0, NULL)); 730 731 dA = M; 732 PetscCall(PetscObjectBaseTypeCompare((PetscObject)M, MATMPIAIJ, &ismpiaij)); 733 if (ismpiaij) PetscCall(MatMPIAIJGetSeqAIJ(M, &dA, &oA, NULL)); 734 735 PetscCall(hypreCSRMatrixGetCOO_Private(dH, &coo_n, &coo_i, &coo_j)); 736 PetscCall(MatSetPreallocationCOO(dA, coo_n, coo_i, coo_j)); 737 PetscCall(PetscFree2(coo_i, coo_j)); 738 if (ismpiaij) { 739 HYPRE_Int nc = hypre_CSRMatrixNumCols(oH); 740 741 PetscCall(PetscLayoutDestroy(&oA->cmap)); 742 PetscCall(PetscLayoutCreateFromSizes(PetscObjectComm((PetscObject)oA), nc, nc, 1, &oA->cmap)); 743 PetscCall(hypreCSRMatrixGetCOO_Private(oH, &coo_n, &coo_i, &coo_j)); 744 PetscCall(MatSetPreallocationCOO(oA, coo_n, coo_i, coo_j)); 745 PetscCall(PetscFree2(coo_i, coo_j)); 746 747 /* garray */ 748 Mat_MPIAIJ *aij = (Mat_MPIAIJ *)(M->data); 749 HYPRE_BigInt *harray = hypre_ParCSRMatrixColMapOffd(parcsr); 750 PetscInt *garray; 751 752 PetscCall(PetscFree(aij->garray)); 753 PetscCall(PetscMalloc1(nc, &garray)); 754 for (HYPRE_Int i = 0; i < nc; i++) garray[i] = (PetscInt)harray[i]; 755 aij->garray = garray; 756 PetscCall(MatSetUpMultiply_MPIAIJ(M)); 757 } 758 if (reuse == MAT_INITIAL_MATRIX) *B = M; 759 } else M = *B; 760 761 dA = M; 762 PetscCall(PetscObjectBaseTypeCompare((PetscObject)M, MATMPIAIJ, &ismpiaij)); 763 if (ismpiaij) PetscCall(MatMPIAIJGetSeqAIJ(M, &dA, &oA, NULL)); 764 PetscCall(MatSetValuesCOOFromCSRMatrix_Private(dA, dH)); 765 if (oA) PetscCall(MatSetValuesCOOFromCSRMatrix_Private(oA, oH)); 766 M->assembled = PETSC_TRUE; 767 if (reuse == MAT_INPLACE_MATRIX) PetscCall(MatHeaderReplace(A, &M)); 768 PetscFunctionReturn(PETSC_SUCCESS); 769 } 770 771 static PetscErrorCode MatAIJGetParCSR_Private(Mat A, hypre_ParCSRMatrix **hA) 772 { 773 hypre_ParCSRMatrix *tA; 774 hypre_CSRMatrix *hdiag, *hoffd; 775 Mat_SeqAIJ *diag, *offd; 776 PetscInt *garray, i, noffd, dnnz, onnz, *row_starts, *col_starts; 777 MPI_Comm comm = PetscObjectComm((PetscObject)A); 778 PetscBool ismpiaij, isseqaij; 779 PetscBool sameint = (PetscBool)(sizeof(PetscInt) == sizeof(HYPRE_Int)); 780 HYPRE_Int *hdi = NULL, *hdj = NULL, *hoi = NULL, *hoj = NULL; 781 PetscInt *pdi = NULL, *pdj = NULL, *poi = NULL, *poj = NULL; 782 PetscBool iscuda, iship; 783 #if defined(PETSC_HAVE_DEVICE) && defined(PETSC_HAVE_HYPRE_DEVICE) 784 PetscBool boundtocpu = A->boundtocpu; 785 #else 786 PetscBool boundtocpu = PETSC_TRUE; 787 #endif 788 789 PetscFunctionBegin; 790 PetscCall(PetscObjectBaseTypeCompare((PetscObject)A, MATMPIAIJ, &ismpiaij)); 791 PetscCall(PetscObjectBaseTypeCompare((PetscObject)A, MATSEQAIJ, &isseqaij)); 792 PetscCheck(ismpiaij || isseqaij, comm, PETSC_ERR_SUP, "Unsupported type %s", ((PetscObject)A)->type_name); 793 PetscCall(PetscObjectTypeCompareAny((PetscObject)A, &iscuda, MATSEQAIJHIPSPARSE, MATMPIAIJCUSPARSE, "")); 794 PetscCall(PetscObjectTypeCompareAny((PetscObject)A, &iship, MATSEQAIJCUSPARSE, MATMPIAIJHIPSPARSE, "")); 795 PetscHYPREInitialize(); 796 if (ismpiaij) { 797 Mat_MPIAIJ *a = (Mat_MPIAIJ *)(A->data); 798 799 diag = (Mat_SeqAIJ *)a->A->data; 800 offd = (Mat_SeqAIJ *)a->B->data; 801 if (!boundtocpu && (iscuda || iship)) { 802 #if defined(HYPRE_USING_CUDA) && defined(PETSC_HAVE_CUDA) 803 if (iscuda) { 804 sameint = PETSC_TRUE; 805 PetscCall(MatSeqAIJCUSPARSEGetIJ(a->A, PETSC_FALSE, (const HYPRE_Int **)&hdi, (const HYPRE_Int **)&hdj)); 806 PetscCall(MatSeqAIJCUSPARSEGetIJ(a->B, PETSC_FALSE, (const HYPRE_Int **)&hoi, (const HYPRE_Int **)&hoj)); 807 } 808 #endif 809 #if defined(HYPRE_USING_HIP) && defined(PETSC_HAVE_HIP) 810 if (iship) { 811 sameint = PETSC_TRUE; 812 PetscCall(MatSeqAIJHIPSPARSEGetIJ(a->A, PETSC_FALSE, (const HYPRE_Int **)&hdi, (const HYPRE_Int **)&hdj)); 813 PetscCall(MatSeqAIJHIPSPARSEGetIJ(a->B, PETSC_FALSE, (const HYPRE_Int **)&hoi, (const HYPRE_Int **)&hoj)); 814 } 815 #endif 816 } else { 817 boundtocpu = PETSC_TRUE; 818 pdi = diag->i; 819 pdj = diag->j; 820 poi = offd->i; 821 poj = offd->j; 822 if (sameint) { 823 hdi = (HYPRE_Int *)pdi; 824 hdj = (HYPRE_Int *)pdj; 825 hoi = (HYPRE_Int *)poi; 826 hoj = (HYPRE_Int *)poj; 827 } 828 } 829 garray = a->garray; 830 noffd = a->B->cmap->N; 831 dnnz = diag->nz; 832 onnz = offd->nz; 833 } else { 834 diag = (Mat_SeqAIJ *)A->data; 835 offd = NULL; 836 if (!boundtocpu && (iscuda || iship)) { 837 #if defined(HYPRE_USING_CUDA) && defined(PETSC_HAVE_CUDA) 838 if (iscuda) { 839 sameint = PETSC_TRUE; 840 PetscCall(MatSeqAIJCUSPARSEGetIJ(A, PETSC_FALSE, (const HYPRE_Int **)&hdi, (const HYPRE_Int **)&hdj)); 841 } 842 #endif 843 #if defined(HYPRE_USING_HIP) && defined(PETSC_HAVE_HIP) 844 if (iship) { 845 sameint = PETSC_TRUE; 846 PetscCall(MatSeqAIJHIPSPARSEGetIJ(A, PETSC_FALSE, (const HYPRE_Int **)&hdi, (const HYPRE_Int **)&hdj)); 847 } 848 #endif 849 } else { 850 boundtocpu = PETSC_TRUE; 851 pdi = diag->i; 852 pdj = diag->j; 853 if (sameint) { 854 hdi = (HYPRE_Int *)pdi; 855 hdj = (HYPRE_Int *)pdj; 856 } 857 } 858 garray = NULL; 859 noffd = 0; 860 dnnz = diag->nz; 861 onnz = 0; 862 } 863 864 /* create a temporary ParCSR */ 865 if (HYPRE_AssumedPartitionCheck()) { 866 PetscMPIInt myid; 867 868 PetscCallMPI(MPI_Comm_rank(comm, &myid)); 869 row_starts = A->rmap->range + myid; 870 col_starts = A->cmap->range + myid; 871 } else { 872 row_starts = A->rmap->range; 873 col_starts = A->cmap->range; 874 } 875 tA = hypre_ParCSRMatrixCreate(comm, A->rmap->N, A->cmap->N, (HYPRE_BigInt *)row_starts, (HYPRE_BigInt *)col_starts, noffd, dnnz, onnz); 876 #if defined(hypre_ParCSRMatrixOwnsRowStarts) 877 hypre_ParCSRMatrixSetRowStartsOwner(tA, 0); 878 hypre_ParCSRMatrixSetColStartsOwner(tA, 0); 879 #endif 880 881 /* set diagonal part */ 882 hdiag = hypre_ParCSRMatrixDiag(tA); 883 if (!sameint) { /* malloc CSR pointers */ 884 PetscCall(PetscMalloc2(A->rmap->n + 1, &hdi, dnnz, &hdj)); 885 for (i = 0; i < A->rmap->n + 1; i++) hdi[i] = (HYPRE_Int)(pdi[i]); 886 for (i = 0; i < dnnz; i++) hdj[i] = (HYPRE_Int)(pdj[i]); 887 } 888 hypre_CSRMatrixI(hdiag) = hdi; 889 hypre_CSRMatrixJ(hdiag) = hdj; 890 hypre_CSRMatrixData(hdiag) = (HYPRE_Complex *)diag->a; 891 hypre_CSRMatrixNumNonzeros(hdiag) = diag->nz; 892 hypre_CSRMatrixSetRownnz(hdiag); 893 hypre_CSRMatrixSetDataOwner(hdiag, 0); 894 895 /* set off-diagonal part */ 896 hoffd = hypre_ParCSRMatrixOffd(tA); 897 if (offd) { 898 if (!sameint) { /* malloc CSR pointers */ 899 PetscCall(PetscMalloc2(A->rmap->n + 1, &hoi, onnz, &hoj)); 900 for (i = 0; i < A->rmap->n + 1; i++) hoi[i] = (HYPRE_Int)(poi[i]); 901 for (i = 0; i < onnz; i++) hoj[i] = (HYPRE_Int)(poj[i]); 902 } 903 hypre_CSRMatrixI(hoffd) = hoi; 904 hypre_CSRMatrixJ(hoffd) = hoj; 905 hypre_CSRMatrixData(hoffd) = (HYPRE_Complex *)offd->a; 906 hypre_CSRMatrixNumNonzeros(hoffd) = offd->nz; 907 hypre_CSRMatrixSetRownnz(hoffd); 908 hypre_CSRMatrixSetDataOwner(hoffd, 0); 909 } 910 #if defined(PETSC_HAVE_HYPRE_DEVICE) 911 PetscCallExternal(hypre_ParCSRMatrixInitialize_v2, tA, !boundtocpu ? HYPRE_MEMORY_DEVICE : HYPRE_MEMORY_HOST); 912 #else 913 #if PETSC_PKG_HYPRE_VERSION_LT(2, 18, 0) 914 PetscCallExternal(hypre_ParCSRMatrixInitialize, tA); 915 #else 916 PetscCallExternal(hypre_ParCSRMatrixInitialize_v2, tA, HYPRE_MEMORY_HOST); 917 #endif 918 #endif 919 hypre_TFree(hypre_ParCSRMatrixColMapOffd(tA), HYPRE_MEMORY_HOST); 920 hypre_ParCSRMatrixSetNumNonzeros(tA); 921 hypre_ParCSRMatrixColMapOffd(tA) = (HYPRE_BigInt *)garray; 922 if (!hypre_ParCSRMatrixCommPkg(tA)) PetscCallExternal(hypre_MatvecCommPkgCreate, tA); 923 *hA = tA; 924 PetscFunctionReturn(PETSC_SUCCESS); 925 } 926 927 static PetscErrorCode MatAIJRestoreParCSR_Private(Mat A, hypre_ParCSRMatrix **hA) 928 { 929 hypre_CSRMatrix *hdiag, *hoffd; 930 PetscBool ismpiaij, sameint = (PetscBool)(sizeof(PetscInt) == sizeof(HYPRE_Int)); 931 #if defined(PETSC_HAVE_HYPRE_DEVICE) 932 PetscBool iscuda = PETSC_FALSE; 933 #endif 934 935 PetscFunctionBegin; 936 PetscCall(PetscObjectBaseTypeCompare((PetscObject)A, MATMPIAIJ, &ismpiaij)); 937 #if defined(PETSC_HAVE_HYPRE_DEVICE) 938 PetscCall(PetscObjectTypeCompareAny((PetscObject)A, &iscuda, MATSEQAIJCUSPARSE, MATMPIAIJCUSPARSE, "")); 939 if (iscuda) sameint = PETSC_TRUE; 940 #endif 941 hdiag = hypre_ParCSRMatrixDiag(*hA); 942 hoffd = hypre_ParCSRMatrixOffd(*hA); 943 /* free temporary memory allocated by PETSc 944 set pointers to NULL before destroying tA */ 945 if (!sameint) { 946 HYPRE_Int *hi, *hj; 947 948 hi = hypre_CSRMatrixI(hdiag); 949 hj = hypre_CSRMatrixJ(hdiag); 950 PetscCall(PetscFree2(hi, hj)); 951 if (ismpiaij) { 952 hi = hypre_CSRMatrixI(hoffd); 953 hj = hypre_CSRMatrixJ(hoffd); 954 PetscCall(PetscFree2(hi, hj)); 955 } 956 } 957 hypre_CSRMatrixI(hdiag) = NULL; 958 hypre_CSRMatrixJ(hdiag) = NULL; 959 hypre_CSRMatrixData(hdiag) = NULL; 960 if (ismpiaij) { 961 hypre_CSRMatrixI(hoffd) = NULL; 962 hypre_CSRMatrixJ(hoffd) = NULL; 963 hypre_CSRMatrixData(hoffd) = NULL; 964 } 965 hypre_ParCSRMatrixColMapOffd(*hA) = NULL; 966 hypre_ParCSRMatrixDestroy(*hA); 967 *hA = NULL; 968 PetscFunctionReturn(PETSC_SUCCESS); 969 } 970 971 /* calls RAP from BoomerAMG: 972 the resulting ParCSR will not own the column and row starts 973 It looks like we don't need to have the diagonal entries ordered first */ 974 static PetscErrorCode MatHYPRE_ParCSR_RAP(hypre_ParCSRMatrix *hR, hypre_ParCSRMatrix *hA, hypre_ParCSRMatrix *hP, hypre_ParCSRMatrix **hRAP) 975 { 976 #if defined(hypre_ParCSRMatrixOwnsRowStarts) 977 HYPRE_Int P_owns_col_starts, R_owns_row_starts; 978 #endif 979 980 PetscFunctionBegin; 981 #if defined(hypre_ParCSRMatrixOwnsRowStarts) 982 P_owns_col_starts = hypre_ParCSRMatrixOwnsColStarts(hP); 983 R_owns_row_starts = hypre_ParCSRMatrixOwnsRowStarts(hR); 984 #endif 985 /* can be replaced by version test later */ 986 #if defined(PETSC_HAVE_HYPRE_DEVICE) 987 PetscStackPushExternal("hypre_ParCSRMatrixRAP"); 988 *hRAP = hypre_ParCSRMatrixRAP(hR, hA, hP); 989 PetscStackPop; 990 #else 991 PetscCallExternal(hypre_BoomerAMGBuildCoarseOperator, hR, hA, hP, hRAP); 992 PetscCallExternal(hypre_ParCSRMatrixSetNumNonzeros, *hRAP); 993 #endif 994 /* hypre_BoomerAMGBuildCoarseOperator steals the col_starts from P and the row_starts from R */ 995 #if defined(hypre_ParCSRMatrixOwnsRowStarts) 996 hypre_ParCSRMatrixSetRowStartsOwner(*hRAP, 0); 997 hypre_ParCSRMatrixSetColStartsOwner(*hRAP, 0); 998 if (P_owns_col_starts) hypre_ParCSRMatrixSetColStartsOwner(hP, 1); 999 if (R_owns_row_starts) hypre_ParCSRMatrixSetRowStartsOwner(hR, 1); 1000 #endif 1001 PetscFunctionReturn(PETSC_SUCCESS); 1002 } 1003 1004 static PetscErrorCode MatPtAPNumeric_AIJ_AIJ_wHYPRE(Mat A, Mat P, Mat C) 1005 { 1006 Mat B; 1007 hypre_ParCSRMatrix *hA, *hP, *hPtAP = NULL; 1008 Mat_Product *product = C->product; 1009 1010 PetscFunctionBegin; 1011 PetscCall(MatAIJGetParCSR_Private(A, &hA)); 1012 PetscCall(MatAIJGetParCSR_Private(P, &hP)); 1013 PetscCall(MatHYPRE_ParCSR_RAP(hP, hA, hP, &hPtAP)); 1014 PetscCall(MatCreateFromParCSR(hPtAP, MATAIJ, PETSC_OWN_POINTER, &B)); 1015 1016 PetscCall(MatHeaderMerge(C, &B)); 1017 C->product = product; 1018 1019 PetscCall(MatAIJRestoreParCSR_Private(A, &hA)); 1020 PetscCall(MatAIJRestoreParCSR_Private(P, &hP)); 1021 PetscFunctionReturn(PETSC_SUCCESS); 1022 } 1023 1024 PETSC_INTERN PetscErrorCode MatPtAPSymbolic_AIJ_AIJ_wHYPRE(Mat A, Mat P, PetscReal fill, Mat C) 1025 { 1026 PetscFunctionBegin; 1027 PetscCall(MatSetType(C, MATAIJ)); 1028 C->ops->ptapnumeric = MatPtAPNumeric_AIJ_AIJ_wHYPRE; 1029 C->ops->productnumeric = MatProductNumeric_PtAP; 1030 PetscFunctionReturn(PETSC_SUCCESS); 1031 } 1032 1033 static PetscErrorCode MatPtAPNumeric_AIJ_HYPRE(Mat A, Mat P, Mat C) 1034 { 1035 Mat B; 1036 Mat_HYPRE *hP; 1037 hypre_ParCSRMatrix *hA = NULL, *Pparcsr, *ptapparcsr = NULL; 1038 HYPRE_Int type; 1039 MPI_Comm comm = PetscObjectComm((PetscObject)A); 1040 PetscBool ishypre; 1041 1042 PetscFunctionBegin; 1043 PetscCall(PetscObjectTypeCompare((PetscObject)P, MATHYPRE, &ishypre)); 1044 PetscCheck(ishypre, comm, PETSC_ERR_USER, "P should be of type %s", MATHYPRE); 1045 hP = (Mat_HYPRE *)P->data; 1046 PetscCallExternal(HYPRE_IJMatrixGetObjectType, hP->ij, &type); 1047 PetscCheck(type == HYPRE_PARCSR, comm, PETSC_ERR_SUP, "Only HYPRE_PARCSR is supported"); 1048 PetscCallExternal(HYPRE_IJMatrixGetObject, hP->ij, (void **)&Pparcsr); 1049 1050 PetscCall(MatAIJGetParCSR_Private(A, &hA)); 1051 PetscCall(MatHYPRE_ParCSR_RAP(Pparcsr, hA, Pparcsr, &ptapparcsr)); 1052 PetscCall(MatAIJRestoreParCSR_Private(A, &hA)); 1053 1054 /* create temporary matrix and merge to C */ 1055 PetscCall(MatCreateFromParCSR(ptapparcsr, ((PetscObject)C)->type_name, PETSC_OWN_POINTER, &B)); 1056 PetscCall(MatHeaderMerge(C, &B)); 1057 PetscFunctionReturn(PETSC_SUCCESS); 1058 } 1059 1060 static PetscErrorCode MatPtAPNumeric_HYPRE_HYPRE(Mat A, Mat P, Mat C) 1061 { 1062 Mat B; 1063 hypre_ParCSRMatrix *Aparcsr, *Pparcsr, *ptapparcsr = NULL; 1064 Mat_HYPRE *hA, *hP; 1065 PetscBool ishypre; 1066 HYPRE_Int type; 1067 1068 PetscFunctionBegin; 1069 PetscCall(PetscObjectTypeCompare((PetscObject)P, MATHYPRE, &ishypre)); 1070 PetscCheck(ishypre, PetscObjectComm((PetscObject)P), PETSC_ERR_USER, "P should be of type %s", MATHYPRE); 1071 PetscCall(PetscObjectTypeCompare((PetscObject)A, MATHYPRE, &ishypre)); 1072 PetscCheck(ishypre, PetscObjectComm((PetscObject)A), PETSC_ERR_USER, "A should be of type %s", MATHYPRE); 1073 hA = (Mat_HYPRE *)A->data; 1074 hP = (Mat_HYPRE *)P->data; 1075 PetscCallExternal(HYPRE_IJMatrixGetObjectType, hA->ij, &type); 1076 PetscCheck(type == HYPRE_PARCSR, PetscObjectComm((PetscObject)A), PETSC_ERR_SUP, "Only HYPRE_PARCSR is supported"); 1077 PetscCallExternal(HYPRE_IJMatrixGetObjectType, hP->ij, &type); 1078 PetscCheck(type == HYPRE_PARCSR, PetscObjectComm((PetscObject)P), PETSC_ERR_SUP, "Only HYPRE_PARCSR is supported"); 1079 PetscCallExternal(HYPRE_IJMatrixGetObject, hA->ij, (void **)&Aparcsr); 1080 PetscCallExternal(HYPRE_IJMatrixGetObject, hP->ij, (void **)&Pparcsr); 1081 PetscCall(MatHYPRE_ParCSR_RAP(Pparcsr, Aparcsr, Pparcsr, &ptapparcsr)); 1082 PetscCall(MatCreateFromParCSR(ptapparcsr, MATHYPRE, PETSC_OWN_POINTER, &B)); 1083 PetscCall(MatHeaderMerge(C, &B)); 1084 PetscFunctionReturn(PETSC_SUCCESS); 1085 } 1086 1087 /* calls hypre_ParMatmul 1088 hypre_ParMatMul uses hypre_ParMatrixCreate with the communicator of hA 1089 hypre_ParMatrixCreate does not duplicate the communicator 1090 It looks like we don't need to have the diagonal entries ordered first */ 1091 static PetscErrorCode MatHYPRE_ParCSR_MatMatMult(hypre_ParCSRMatrix *hA, hypre_ParCSRMatrix *hB, hypre_ParCSRMatrix **hAB) 1092 { 1093 PetscFunctionBegin; 1094 /* can be replaced by version test later */ 1095 #if defined(PETSC_HAVE_HYPRE_DEVICE) 1096 PetscStackPushExternal("hypre_ParCSRMatMat"); 1097 *hAB = hypre_ParCSRMatMat(hA, hB); 1098 #else 1099 PetscStackPushExternal("hypre_ParMatmul"); 1100 *hAB = hypre_ParMatmul(hA, hB); 1101 #endif 1102 PetscStackPop; 1103 PetscFunctionReturn(PETSC_SUCCESS); 1104 } 1105 1106 static PetscErrorCode MatMatMultNumeric_AIJ_AIJ_wHYPRE(Mat A, Mat B, Mat C) 1107 { 1108 Mat D; 1109 hypre_ParCSRMatrix *hA, *hB, *hAB = NULL; 1110 Mat_Product *product = C->product; 1111 1112 PetscFunctionBegin; 1113 PetscCall(MatAIJGetParCSR_Private(A, &hA)); 1114 PetscCall(MatAIJGetParCSR_Private(B, &hB)); 1115 PetscCall(MatHYPRE_ParCSR_MatMatMult(hA, hB, &hAB)); 1116 PetscCall(MatCreateFromParCSR(hAB, MATAIJ, PETSC_OWN_POINTER, &D)); 1117 1118 PetscCall(MatHeaderMerge(C, &D)); 1119 C->product = product; 1120 1121 PetscCall(MatAIJRestoreParCSR_Private(A, &hA)); 1122 PetscCall(MatAIJRestoreParCSR_Private(B, &hB)); 1123 PetscFunctionReturn(PETSC_SUCCESS); 1124 } 1125 1126 PETSC_INTERN PetscErrorCode MatMatMultSymbolic_AIJ_AIJ_wHYPRE(Mat A, Mat B, PetscReal fill, Mat C) 1127 { 1128 PetscFunctionBegin; 1129 PetscCall(MatSetType(C, MATAIJ)); 1130 C->ops->matmultnumeric = MatMatMultNumeric_AIJ_AIJ_wHYPRE; 1131 C->ops->productnumeric = MatProductNumeric_AB; 1132 PetscFunctionReturn(PETSC_SUCCESS); 1133 } 1134 1135 static PetscErrorCode MatMatMultNumeric_HYPRE_HYPRE(Mat A, Mat B, Mat C) 1136 { 1137 Mat D; 1138 hypre_ParCSRMatrix *Aparcsr, *Bparcsr, *ABparcsr = NULL; 1139 Mat_HYPRE *hA, *hB; 1140 PetscBool ishypre; 1141 HYPRE_Int type; 1142 Mat_Product *product; 1143 1144 PetscFunctionBegin; 1145 PetscCall(PetscObjectTypeCompare((PetscObject)B, MATHYPRE, &ishypre)); 1146 PetscCheck(ishypre, PetscObjectComm((PetscObject)B), PETSC_ERR_USER, "B should be of type %s", MATHYPRE); 1147 PetscCall(PetscObjectTypeCompare((PetscObject)A, MATHYPRE, &ishypre)); 1148 PetscCheck(ishypre, PetscObjectComm((PetscObject)A), PETSC_ERR_USER, "A should be of type %s", MATHYPRE); 1149 hA = (Mat_HYPRE *)A->data; 1150 hB = (Mat_HYPRE *)B->data; 1151 PetscCallExternal(HYPRE_IJMatrixGetObjectType, hA->ij, &type); 1152 PetscCheck(type == HYPRE_PARCSR, PetscObjectComm((PetscObject)A), PETSC_ERR_SUP, "Only HYPRE_PARCSR is supported"); 1153 PetscCallExternal(HYPRE_IJMatrixGetObjectType, hB->ij, &type); 1154 PetscCheck(type == HYPRE_PARCSR, PetscObjectComm((PetscObject)B), PETSC_ERR_SUP, "Only HYPRE_PARCSR is supported"); 1155 PetscCallExternal(HYPRE_IJMatrixGetObject, hA->ij, (void **)&Aparcsr); 1156 PetscCallExternal(HYPRE_IJMatrixGetObject, hB->ij, (void **)&Bparcsr); 1157 PetscCall(MatHYPRE_ParCSR_MatMatMult(Aparcsr, Bparcsr, &ABparcsr)); 1158 PetscCall(MatCreateFromParCSR(ABparcsr, MATHYPRE, PETSC_OWN_POINTER, &D)); 1159 1160 /* need to use HeaderReplace because HeaderMerge messes up with the communicator */ 1161 product = C->product; /* save it from MatHeaderReplace() */ 1162 C->product = NULL; 1163 PetscCall(MatHeaderReplace(C, &D)); 1164 C->product = product; 1165 C->ops->matmultnumeric = MatMatMultNumeric_HYPRE_HYPRE; 1166 C->ops->productnumeric = MatProductNumeric_AB; 1167 PetscFunctionReturn(PETSC_SUCCESS); 1168 } 1169 1170 PETSC_INTERN PetscErrorCode MatTransposeMatMatMultNumeric_AIJ_AIJ_AIJ_wHYPRE(Mat A, Mat B, Mat C, Mat D) 1171 { 1172 Mat E; 1173 hypre_ParCSRMatrix *hA, *hB, *hC, *hABC = NULL; 1174 1175 PetscFunctionBegin; 1176 PetscCall(MatAIJGetParCSR_Private(A, &hA)); 1177 PetscCall(MatAIJGetParCSR_Private(B, &hB)); 1178 PetscCall(MatAIJGetParCSR_Private(C, &hC)); 1179 PetscCall(MatHYPRE_ParCSR_RAP(hA, hB, hC, &hABC)); 1180 PetscCall(MatCreateFromParCSR(hABC, MATAIJ, PETSC_OWN_POINTER, &E)); 1181 PetscCall(MatHeaderMerge(D, &E)); 1182 PetscCall(MatAIJRestoreParCSR_Private(A, &hA)); 1183 PetscCall(MatAIJRestoreParCSR_Private(B, &hB)); 1184 PetscCall(MatAIJRestoreParCSR_Private(C, &hC)); 1185 PetscFunctionReturn(PETSC_SUCCESS); 1186 } 1187 1188 PETSC_INTERN PetscErrorCode MatTransposeMatMatMultSymbolic_AIJ_AIJ_AIJ_wHYPRE(Mat A, Mat B, Mat C, PetscReal fill, Mat D) 1189 { 1190 PetscFunctionBegin; 1191 PetscCall(MatSetType(D, MATAIJ)); 1192 PetscFunctionReturn(PETSC_SUCCESS); 1193 } 1194 1195 static PetscErrorCode MatProductSymbolic_AB_HYPRE(Mat C) 1196 { 1197 PetscFunctionBegin; 1198 C->ops->productnumeric = MatProductNumeric_AB; 1199 PetscFunctionReturn(PETSC_SUCCESS); 1200 } 1201 1202 static PetscErrorCode MatProductSetFromOptions_HYPRE_AB(Mat C) 1203 { 1204 Mat_Product *product = C->product; 1205 PetscBool Ahypre; 1206 1207 PetscFunctionBegin; 1208 PetscCall(PetscObjectTypeCompare((PetscObject)product->A, MATHYPRE, &Ahypre)); 1209 if (Ahypre) { /* A is a Hypre matrix */ 1210 PetscCall(MatSetType(C, MATHYPRE)); 1211 C->ops->productsymbolic = MatProductSymbolic_AB_HYPRE; 1212 C->ops->matmultnumeric = MatMatMultNumeric_HYPRE_HYPRE; 1213 PetscFunctionReturn(PETSC_SUCCESS); 1214 } 1215 PetscFunctionReturn(PETSC_SUCCESS); 1216 } 1217 1218 static PetscErrorCode MatProductSymbolic_PtAP_HYPRE(Mat C) 1219 { 1220 PetscFunctionBegin; 1221 C->ops->productnumeric = MatProductNumeric_PtAP; 1222 PetscFunctionReturn(PETSC_SUCCESS); 1223 } 1224 1225 static PetscErrorCode MatProductSetFromOptions_HYPRE_PtAP(Mat C) 1226 { 1227 Mat_Product *product = C->product; 1228 PetscBool flg; 1229 PetscInt type = 0; 1230 const char *outTypes[4] = {"aij", "seqaij", "mpiaij", "hypre"}; 1231 PetscInt ntype = 4; 1232 Mat A = product->A; 1233 PetscBool Ahypre; 1234 1235 PetscFunctionBegin; 1236 PetscCall(PetscObjectTypeCompare((PetscObject)A, MATHYPRE, &Ahypre)); 1237 if (Ahypre) { /* A is a Hypre matrix */ 1238 PetscCall(MatSetType(C, MATHYPRE)); 1239 C->ops->productsymbolic = MatProductSymbolic_PtAP_HYPRE; 1240 C->ops->ptapnumeric = MatPtAPNumeric_HYPRE_HYPRE; 1241 PetscFunctionReturn(PETSC_SUCCESS); 1242 } 1243 1244 /* A is AIJ, P is Hypre, C = PtAP can be either AIJ or Hypre format */ 1245 /* Get runtime option */ 1246 if (product->api_user) { 1247 PetscOptionsBegin(PetscObjectComm((PetscObject)C), ((PetscObject)C)->prefix, "MatPtAP_HYPRE", "Mat"); 1248 PetscCall(PetscOptionsEList("-matptap_hypre_outtype", "MatPtAP outtype", "MatPtAP outtype", outTypes, ntype, outTypes[type], &type, &flg)); 1249 PetscOptionsEnd(); 1250 } else { 1251 PetscOptionsBegin(PetscObjectComm((PetscObject)C), ((PetscObject)C)->prefix, "MatProduct_PtAP_HYPRE", "Mat"); 1252 PetscCall(PetscOptionsEList("-mat_product_algorithm_hypre_outtype", "MatProduct_PtAP outtype", "MatProduct_PtAP", outTypes, ntype, outTypes[type], &type, &flg)); 1253 PetscOptionsEnd(); 1254 } 1255 1256 if (type == 0 || type == 1 || type == 2) { 1257 PetscCall(MatSetType(C, MATAIJ)); 1258 } else if (type == 3) { 1259 PetscCall(MatSetType(C, MATHYPRE)); 1260 } else SETERRQ(PetscObjectComm((PetscObject)C), PETSC_ERR_SUP, "MatPtAP outtype is not supported"); 1261 C->ops->productsymbolic = MatProductSymbolic_PtAP_HYPRE; 1262 C->ops->ptapnumeric = MatPtAPNumeric_AIJ_HYPRE; 1263 PetscFunctionReturn(PETSC_SUCCESS); 1264 } 1265 1266 static PetscErrorCode MatProductSetFromOptions_HYPRE(Mat C) 1267 { 1268 Mat_Product *product = C->product; 1269 1270 PetscFunctionBegin; 1271 switch (product->type) { 1272 case MATPRODUCT_AB: 1273 PetscCall(MatProductSetFromOptions_HYPRE_AB(C)); 1274 break; 1275 case MATPRODUCT_PtAP: 1276 PetscCall(MatProductSetFromOptions_HYPRE_PtAP(C)); 1277 break; 1278 default: 1279 break; 1280 } 1281 PetscFunctionReturn(PETSC_SUCCESS); 1282 } 1283 1284 static PetscErrorCode MatMultTranspose_HYPRE(Mat A, Vec x, Vec y) 1285 { 1286 PetscFunctionBegin; 1287 PetscCall(MatHYPRE_MultKernel_Private(A, 1.0, x, 0.0, y, PETSC_TRUE)); 1288 PetscFunctionReturn(PETSC_SUCCESS); 1289 } 1290 1291 static PetscErrorCode MatMult_HYPRE(Mat A, Vec x, Vec y) 1292 { 1293 PetscFunctionBegin; 1294 PetscCall(MatHYPRE_MultKernel_Private(A, 1.0, x, 0.0, y, PETSC_FALSE)); 1295 PetscFunctionReturn(PETSC_SUCCESS); 1296 } 1297 1298 static PetscErrorCode MatMultAdd_HYPRE(Mat A, Vec x, Vec y, Vec z) 1299 { 1300 PetscFunctionBegin; 1301 if (y != z) PetscCall(VecCopy(y, z)); 1302 PetscCall(MatHYPRE_MultKernel_Private(A, 1.0, x, 1.0, z, PETSC_FALSE)); 1303 PetscFunctionReturn(PETSC_SUCCESS); 1304 } 1305 1306 static PetscErrorCode MatMultTransposeAdd_HYPRE(Mat A, Vec x, Vec y, Vec z) 1307 { 1308 PetscFunctionBegin; 1309 if (y != z) PetscCall(VecCopy(y, z)); 1310 PetscCall(MatHYPRE_MultKernel_Private(A, 1.0, x, 1.0, z, PETSC_TRUE)); 1311 PetscFunctionReturn(PETSC_SUCCESS); 1312 } 1313 1314 /* y = a * A * x + b * y or y = a * A^t * x + b * y depending on trans */ 1315 static PetscErrorCode MatHYPRE_MultKernel_Private(Mat A, HYPRE_Complex a, Vec x, HYPRE_Complex b, Vec y, PetscBool trans) 1316 { 1317 Mat_HYPRE *hA = (Mat_HYPRE *)A->data; 1318 hypre_ParCSRMatrix *parcsr; 1319 hypre_ParVector *hx, *hy; 1320 1321 PetscFunctionBegin; 1322 if (trans) { 1323 PetscCall(VecHYPRE_IJVectorPushVecRead(hA->b, x)); 1324 if (b != 0.0) PetscCall(VecHYPRE_IJVectorPushVec(hA->x, y)); 1325 else PetscCall(VecHYPRE_IJVectorPushVecWrite(hA->x, y)); 1326 PetscCallExternal(HYPRE_IJVectorGetObject, hA->b->ij, (void **)&hx); 1327 PetscCallExternal(HYPRE_IJVectorGetObject, hA->x->ij, (void **)&hy); 1328 } else { 1329 PetscCall(VecHYPRE_IJVectorPushVecRead(hA->x, x)); 1330 if (b != 0.0) PetscCall(VecHYPRE_IJVectorPushVec(hA->b, y)); 1331 else PetscCall(VecHYPRE_IJVectorPushVecWrite(hA->b, y)); 1332 PetscCallExternal(HYPRE_IJVectorGetObject, hA->x->ij, (void **)&hx); 1333 PetscCallExternal(HYPRE_IJVectorGetObject, hA->b->ij, (void **)&hy); 1334 } 1335 PetscCallExternal(HYPRE_IJMatrixGetObject, hA->ij, (void **)&parcsr); 1336 if (trans) { 1337 PetscCallExternal(hypre_ParCSRMatrixMatvecT, a, parcsr, hx, b, hy); 1338 } else { 1339 PetscCallExternal(hypre_ParCSRMatrixMatvec, a, parcsr, hx, b, hy); 1340 } 1341 PetscCall(VecHYPRE_IJVectorPopVec(hA->x)); 1342 PetscCall(VecHYPRE_IJVectorPopVec(hA->b)); 1343 PetscFunctionReturn(PETSC_SUCCESS); 1344 } 1345 1346 static PetscErrorCode MatDestroy_HYPRE(Mat A) 1347 { 1348 Mat_HYPRE *hA = (Mat_HYPRE *)A->data; 1349 1350 PetscFunctionBegin; 1351 PetscCall(VecHYPRE_IJVectorDestroy(&hA->x)); 1352 PetscCall(VecHYPRE_IJVectorDestroy(&hA->b)); 1353 PetscCall(MatHYPRE_DestroyCOOMat(A)); /* must be called before destroying the individual CSR */ 1354 if (hA->ij) { 1355 if (!hA->inner_free) hypre_IJMatrixObject(hA->ij) = NULL; 1356 PetscCallExternal(HYPRE_IJMatrixDestroy, hA->ij); 1357 } 1358 if (hA->comm) PetscCall(PetscCommRestoreComm(PetscObjectComm((PetscObject)A), &hA->comm)); 1359 1360 PetscCall(MatStashDestroy_Private(&A->stash)); 1361 PetscCall(PetscFree(hA->array)); 1362 if (hA->rows_d) PetscStackCallExternalVoid("hypre_Free", hypre_Free(hA->rows_d, HYPRE_MEMORY_DEVICE)); 1363 1364 PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatConvert_hypre_aij_C", NULL)); 1365 PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatConvert_hypre_is_C", NULL)); 1366 PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatProductSetFromOptions_seqaij_hypre_C", NULL)); 1367 PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatProductSetFromOptions_mpiaij_hypre_C", NULL)); 1368 PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatProductSetFromOptions_seqaijhipsparse_hypre_C", NULL)); 1369 PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatProductSetFromOptions_mpiaijhipsparse_hypre_C", NULL)); 1370 PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatProductSetFromOptions_seqaijcusparse_hypre_C", NULL)); 1371 PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatProductSetFromOptions_mpiaijcusparse_hypre_C", NULL)); 1372 PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatHYPRESetPreallocation_C", NULL)); 1373 PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatHYPREGetParCSR_C", NULL)); 1374 PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatSetPreallocationCOO_C", NULL)); 1375 PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatSetValuesCOO_C", NULL)); 1376 PetscCall(PetscFree(A->data)); 1377 PetscFunctionReturn(PETSC_SUCCESS); 1378 } 1379 1380 static PetscErrorCode MatSetUp_HYPRE(Mat A) 1381 { 1382 PetscFunctionBegin; 1383 if (!A->preallocated) PetscCall(MatHYPRESetPreallocation(A, PETSC_DEFAULT, NULL, PETSC_DEFAULT, NULL)); 1384 PetscFunctionReturn(PETSC_SUCCESS); 1385 } 1386 1387 //TODO FIX hypre_CSRMatrixMatvecOutOfPlace 1388 #if defined(PETSC_HAVE_HYPRE_DEVICE) 1389 static PetscErrorCode MatBindToCPU_HYPRE(Mat A, PetscBool bind) 1390 { 1391 Mat_HYPRE *hA = (Mat_HYPRE *)A->data; 1392 HYPRE_MemoryLocation hmem = bind ? HYPRE_MEMORY_HOST : HYPRE_MEMORY_DEVICE; 1393 1394 PetscFunctionBegin; 1395 A->boundtocpu = bind; 1396 if (hA->ij && hypre_IJMatrixAssembleFlag(hA->ij) && hmem != hypre_IJMatrixMemoryLocation(hA->ij)) { 1397 hypre_ParCSRMatrix *parcsr; 1398 PetscCallExternal(HYPRE_IJMatrixGetObject, hA->ij, (void **)&parcsr); 1399 PetscCallExternal(hypre_ParCSRMatrixMigrate, parcsr, hmem); 1400 } 1401 if (hA->x) PetscCall(VecHYPRE_IJBindToCPU(hA->x, bind)); 1402 if (hA->b) PetscCall(VecHYPRE_IJBindToCPU(hA->b, bind)); 1403 PetscFunctionReturn(PETSC_SUCCESS); 1404 } 1405 #endif 1406 1407 static PetscErrorCode MatAssemblyEnd_HYPRE(Mat A, MatAssemblyType mode) 1408 { 1409 Mat_HYPRE *hA = (Mat_HYPRE *)A->data; 1410 PetscMPIInt n; 1411 PetscInt i, j, rstart, ncols, flg; 1412 PetscInt *row, *col; 1413 PetscScalar *val; 1414 1415 PetscFunctionBegin; 1416 PetscCheck(mode != MAT_FLUSH_ASSEMBLY, PetscObjectComm((PetscObject)A), PETSC_ERR_SUP, "MAT_FLUSH_ASSEMBLY currently not supported with MATHYPRE"); 1417 1418 if (!A->nooffprocentries) { 1419 while (1) { 1420 PetscCall(MatStashScatterGetMesg_Private(&A->stash, &n, &row, &col, &val, &flg)); 1421 if (!flg) break; 1422 1423 for (i = 0; i < n;) { 1424 /* Now identify the consecutive vals belonging to the same row */ 1425 for (j = i, rstart = row[j]; j < n; j++) { 1426 if (row[j] != rstart) break; 1427 } 1428 if (j < n) ncols = j - i; 1429 else ncols = n - i; 1430 /* Now assemble all these values with a single function call */ 1431 PetscCall(MatSetValues_HYPRE(A, 1, row + i, ncols, col + i, val + i, A->insertmode)); 1432 1433 i = j; 1434 } 1435 } 1436 PetscCall(MatStashScatterEnd_Private(&A->stash)); 1437 } 1438 1439 PetscCallExternal(HYPRE_IJMatrixAssemble, hA->ij); 1440 /* The assembly routine destroys the aux_matrix, we recreate it here by calling HYPRE_IJMatrixInitialize */ 1441 /* If the option MAT_SORTED_FULL is set to true, the indices and values can be passed to hypre directly, so we don't need the aux_matrix */ 1442 if (!A->sortedfull) { 1443 hypre_AuxParCSRMatrix *aux_matrix; 1444 1445 /* call destroy just to make sure we do not leak anything */ 1446 aux_matrix = (hypre_AuxParCSRMatrix *)hypre_IJMatrixTranslator(hA->ij); 1447 PetscCallExternal(hypre_AuxParCSRMatrixDestroy, aux_matrix); 1448 hypre_IJMatrixTranslator(hA->ij) = NULL; 1449 1450 /* Initialize with assembled flag -> it only recreates the aux_par_matrix */ 1451 PetscCallExternal(HYPRE_IJMatrixInitialize, hA->ij); 1452 aux_matrix = (hypre_AuxParCSRMatrix *)hypre_IJMatrixTranslator(hA->ij); 1453 if (aux_matrix) { 1454 hypre_AuxParCSRMatrixNeedAux(aux_matrix) = 1; /* see comment in MatHYPRESetPreallocation_HYPRE */ 1455 #if PETSC_PKG_HYPRE_VERSION_LT(2, 19, 0) 1456 PetscCallExternal(hypre_AuxParCSRMatrixInitialize, aux_matrix); 1457 #else 1458 PetscCallExternal(hypre_AuxParCSRMatrixInitialize_v2, aux_matrix, HYPRE_MEMORY_HOST); 1459 #endif 1460 } 1461 } 1462 { 1463 hypre_ParCSRMatrix *parcsr; 1464 1465 PetscCallExternal(HYPRE_IJMatrixGetObject, hA->ij, (void **)&parcsr); 1466 if (!hypre_ParCSRMatrixCommPkg(parcsr)) PetscCallExternal(hypre_MatvecCommPkgCreate, parcsr); 1467 } 1468 if (!hA->x) PetscCall(VecHYPRE_IJVectorCreate(A->cmap, &hA->x)); 1469 if (!hA->b) PetscCall(VecHYPRE_IJVectorCreate(A->rmap, &hA->b)); 1470 #if defined(PETSC_HAVE_HYPRE_DEVICE) 1471 PetscCall(MatBindToCPU_HYPRE(A, A->boundtocpu)); 1472 #endif 1473 PetscFunctionReturn(PETSC_SUCCESS); 1474 } 1475 1476 static PetscErrorCode MatGetArray_HYPRE(Mat A, PetscInt size, void **array) 1477 { 1478 Mat_HYPRE *hA = (Mat_HYPRE *)A->data; 1479 1480 PetscFunctionBegin; 1481 PetscCheck(hA->array_available, PETSC_COMM_SELF, PETSC_ERR_ARG_NULL, "Temporary space is in use"); 1482 1483 if (hA->array_size >= size) { 1484 *array = hA->array; 1485 } else { 1486 PetscCall(PetscFree(hA->array)); 1487 hA->array_size = size; 1488 PetscCall(PetscMalloc(hA->array_size, &hA->array)); 1489 *array = hA->array; 1490 } 1491 1492 hA->array_available = PETSC_FALSE; 1493 PetscFunctionReturn(PETSC_SUCCESS); 1494 } 1495 1496 static PetscErrorCode MatRestoreArray_HYPRE(Mat A, void **array) 1497 { 1498 Mat_HYPRE *hA = (Mat_HYPRE *)A->data; 1499 1500 PetscFunctionBegin; 1501 *array = NULL; 1502 hA->array_available = PETSC_TRUE; 1503 PetscFunctionReturn(PETSC_SUCCESS); 1504 } 1505 1506 static PetscErrorCode MatSetValues_HYPRE(Mat A, PetscInt nr, const PetscInt rows[], PetscInt nc, const PetscInt cols[], const PetscScalar v[], InsertMode ins) 1507 { 1508 Mat_HYPRE *hA = (Mat_HYPRE *)A->data; 1509 PetscScalar *vals = (PetscScalar *)v; 1510 HYPRE_Complex *sscr; 1511 PetscInt *cscr[2]; 1512 PetscInt i, nzc; 1513 PetscInt rst = A->rmap->rstart, ren = A->rmap->rend; 1514 void *array = NULL; 1515 1516 PetscFunctionBegin; 1517 PetscCall(MatGetArray_HYPRE(A, sizeof(PetscInt) * (2 * nc) + sizeof(HYPRE_Complex) * nc * nr, &array)); 1518 cscr[0] = (PetscInt *)array; 1519 cscr[1] = ((PetscInt *)array) + nc; 1520 sscr = (HYPRE_Complex *)(((PetscInt *)array) + nc * 2); 1521 for (i = 0, nzc = 0; i < nc; i++) { 1522 if (cols[i] >= 0) { 1523 cscr[0][nzc] = cols[i]; 1524 cscr[1][nzc++] = i; 1525 } 1526 } 1527 if (!nzc) { 1528 PetscCall(MatRestoreArray_HYPRE(A, &array)); 1529 PetscFunctionReturn(PETSC_SUCCESS); 1530 } 1531 1532 #if 0 //defined(PETSC_HAVE_HYPRE_DEVICE) 1533 if (HYPRE_MEMORY_HOST != hypre_IJMatrixMemoryLocation(hA->ij)) { 1534 hypre_ParCSRMatrix *parcsr; 1535 1536 PetscCallExternal(HYPRE_IJMatrixGetObject,hA->ij,(void**)&parcsr); 1537 PetscCallExternal(hypre_ParCSRMatrixMigrate,parcsr, HYPRE_MEMORY_HOST); 1538 } 1539 #endif 1540 1541 if (ins == ADD_VALUES) { 1542 for (i = 0; i < nr; i++) { 1543 if (rows[i] >= 0) { 1544 PetscInt j; 1545 HYPRE_Int hnc = (HYPRE_Int)nzc; 1546 1547 if (!nzc) continue; 1548 /* nonlocal values */ 1549 if (rows[i] < rst || rows[i] >= ren) { 1550 PetscCheck(!A->nooffprocentries, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Setting off process row %" PetscInt_FMT " even though MatSetOption(,MAT_NO_OFF_PROC_ENTRIES,PETSC_TRUE) was set", rows[i]); 1551 if (hA->donotstash) continue; 1552 } 1553 PetscCheck((PetscInt)hnc == nzc, PETSC_COMM_SELF, PETSC_ERR_SUP, "Hypre overflow! number of columns %" PetscInt_FMT " for row %" PetscInt_FMT, nzc, rows[i]); 1554 for (j = 0; j < nzc; j++) PetscCall(PetscHYPREScalarCast(vals[cscr[1][j]], &sscr[j])); 1555 PetscCallExternal(HYPRE_IJMatrixAddToValues, hA->ij, 1, &hnc, (HYPRE_BigInt *)(rows + i), (HYPRE_BigInt *)cscr[0], sscr); 1556 } 1557 vals += nc; 1558 } 1559 } else { /* INSERT_VALUES */ 1560 for (i = 0; i < nr; i++) { 1561 if (rows[i] >= 0) { 1562 PetscInt j; 1563 HYPRE_Int hnc = (HYPRE_Int)nzc; 1564 1565 if (!nzc) continue; 1566 PetscCheck((PetscInt)hnc == nzc, PETSC_COMM_SELF, PETSC_ERR_SUP, "Hypre overflow! number of columns %" PetscInt_FMT " for row %" PetscInt_FMT, nzc, rows[i]); 1567 for (j = 0; j < nzc; j++) PetscCall(PetscHYPREScalarCast(vals[cscr[1][j]], &sscr[j])); 1568 /* nonlocal values */ 1569 if (rows[i] < rst || rows[i] >= ren) { 1570 PetscCheck(!A->nooffprocentries, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Setting off process row %" PetscInt_FMT " even though MatSetOption(,MAT_NO_OFF_PROC_ENTRIES,PETSC_TRUE) was set", rows[i]); 1571 if (!hA->donotstash) PetscCall(MatStashValuesRow_Private(&A->stash, rows[i], nzc, cscr[0], (PetscScalar *)sscr, PETSC_FALSE)); 1572 } 1573 /* local values */ 1574 else 1575 PetscCallExternal(HYPRE_IJMatrixSetValues, hA->ij, 1, &hnc, (HYPRE_BigInt *)(rows + i), (HYPRE_BigInt *)cscr[0], sscr); 1576 } 1577 vals += nc; 1578 } 1579 } 1580 1581 PetscCall(MatRestoreArray_HYPRE(A, &array)); 1582 PetscFunctionReturn(PETSC_SUCCESS); 1583 } 1584 1585 static PetscErrorCode MatHYPRESetPreallocation_HYPRE(Mat A, PetscInt dnz, const PetscInt dnnz[], PetscInt onz, const PetscInt onnz[]) 1586 { 1587 Mat_HYPRE *hA = (Mat_HYPRE *)A->data; 1588 HYPRE_Int *hdnnz, *honnz; 1589 PetscInt i, rs, re, cs, ce, bs; 1590 PetscMPIInt size; 1591 1592 PetscFunctionBegin; 1593 PetscCall(PetscLayoutSetUp(A->rmap)); 1594 PetscCall(PetscLayoutSetUp(A->cmap)); 1595 rs = A->rmap->rstart; 1596 re = A->rmap->rend; 1597 cs = A->cmap->rstart; 1598 ce = A->cmap->rend; 1599 if (!hA->ij) { 1600 PetscCallExternal(HYPRE_IJMatrixCreate, hA->comm, rs, re - 1, cs, ce - 1, &hA->ij); 1601 PetscCallExternal(HYPRE_IJMatrixSetObjectType, hA->ij, HYPRE_PARCSR); 1602 } else { 1603 HYPRE_BigInt hrs, hre, hcs, hce; 1604 PetscCallExternal(HYPRE_IJMatrixGetLocalRange, hA->ij, &hrs, &hre, &hcs, &hce); 1605 PetscCheck(hre - hrs + 1 == re - rs, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Inconsistent local rows: IJMatrix [%" PetscHYPRE_BigInt_FMT ",%" PetscHYPRE_BigInt_FMT "), PETSc [%" PetscInt_FMT ",%" PetscInt_FMT ")", hrs, hre + 1, rs, re); 1606 PetscCheck(hce - hcs + 1 == ce - cs, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Inconsistent local cols: IJMatrix [%" PetscHYPRE_BigInt_FMT ",%" PetscHYPRE_BigInt_FMT "), PETSc [%" PetscInt_FMT ",%" PetscInt_FMT ")", hcs, hce + 1, cs, ce); 1607 } 1608 PetscCall(MatHYPRE_DestroyCOOMat(A)); 1609 PetscCall(MatGetBlockSize(A, &bs)); 1610 if (dnz == PETSC_DEFAULT || dnz == PETSC_DECIDE) dnz = 10 * bs; 1611 if (onz == PETSC_DEFAULT || onz == PETSC_DECIDE) onz = 10 * bs; 1612 1613 if (!dnnz) { 1614 PetscCall(PetscMalloc1(A->rmap->n, &hdnnz)); 1615 for (i = 0; i < A->rmap->n; i++) hdnnz[i] = dnz; 1616 } else { 1617 hdnnz = (HYPRE_Int *)dnnz; 1618 } 1619 PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)A), &size)); 1620 if (size > 1) { 1621 hypre_AuxParCSRMatrix *aux_matrix; 1622 if (!onnz) { 1623 PetscCall(PetscMalloc1(A->rmap->n, &honnz)); 1624 for (i = 0; i < A->rmap->n; i++) honnz[i] = onz; 1625 } else honnz = (HYPRE_Int *)onnz; 1626 /* SetDiagOffdSizes sets hypre_AuxParCSRMatrixNeedAux(aux_matrix) = 0, since it seems 1627 they assume the user will input the entire row values, properly sorted 1628 In PETSc, we don't make such an assumption and set this flag to 1, 1629 unless the option MAT_SORTED_FULL is set to true. 1630 Also, to avoid possible memory leaks, we destroy and recreate the translator 1631 This has to be done here, as HYPRE_IJMatrixInitialize will properly initialize 1632 the IJ matrix for us */ 1633 aux_matrix = (hypre_AuxParCSRMatrix *)hypre_IJMatrixTranslator(hA->ij); 1634 hypre_AuxParCSRMatrixDestroy(aux_matrix); 1635 hypre_IJMatrixTranslator(hA->ij) = NULL; 1636 PetscCallExternal(HYPRE_IJMatrixSetDiagOffdSizes, hA->ij, hdnnz, honnz); 1637 aux_matrix = (hypre_AuxParCSRMatrix *)hypre_IJMatrixTranslator(hA->ij); 1638 hypre_AuxParCSRMatrixNeedAux(aux_matrix) = !A->sortedfull; 1639 } else { 1640 honnz = NULL; 1641 PetscCallExternal(HYPRE_IJMatrixSetRowSizes, hA->ij, hdnnz); 1642 } 1643 1644 /* reset assembled flag and call the initialize method */ 1645 hypre_IJMatrixAssembleFlag(hA->ij) = 0; 1646 #if PETSC_PKG_HYPRE_VERSION_LT(2, 19, 0) 1647 PetscCallExternal(HYPRE_IJMatrixInitialize, hA->ij); 1648 #else 1649 PetscCallExternal(HYPRE_IJMatrixInitialize_v2, hA->ij, HYPRE_MEMORY_HOST); 1650 #endif 1651 if (!dnnz) PetscCall(PetscFree(hdnnz)); 1652 if (!onnz && honnz) PetscCall(PetscFree(honnz)); 1653 /* Match AIJ logic */ 1654 A->preallocated = PETSC_TRUE; 1655 A->assembled = PETSC_FALSE; 1656 PetscFunctionReturn(PETSC_SUCCESS); 1657 } 1658 1659 /*@C 1660 MatHYPRESetPreallocation - Preallocates memory for a sparse parallel matrix in HYPRE IJ format 1661 1662 Collective 1663 1664 Input Parameters: 1665 + A - the matrix 1666 . dnz - number of nonzeros per row in DIAGONAL portion of local submatrix 1667 (same value is used for all local rows) 1668 . dnnz - array containing the number of nonzeros in the various rows of the 1669 DIAGONAL portion of the local submatrix (possibly different for each row) 1670 or `NULL` (`PETSC_NULL_INTEGER` in Fortran), if `d_nz` is used to specify the nonzero structure. 1671 The size of this array is equal to the number of local rows, i.e `m`. 1672 For matrices that will be factored, you must leave room for (and set) 1673 the diagonal entry even if it is zero. 1674 . onz - number of nonzeros per row in the OFF-DIAGONAL portion of local 1675 submatrix (same value is used for all local rows). 1676 - onnz - array containing the number of nonzeros in the various rows of the 1677 OFF-DIAGONAL portion of the local submatrix (possibly different for 1678 each row) or `NULL` (`PETSC_NULL_INTEGER` in Fortran), if `o_nz` is used to specify the nonzero 1679 structure. The size of this array is equal to the number 1680 of local rows, i.e `m`. 1681 1682 Level: intermediate 1683 1684 Note: 1685 If the *nnz parameter is given then the *nz parameter is ignored; for sequential matrices, `onz` and `onnz` are ignored. 1686 1687 .seealso: [](ch_matrices), `Mat`, `MatCreate()`, `MatMPIAIJSetPreallocation()`, `MATHYPRE`, `MATAIJ` 1688 @*/ 1689 PetscErrorCode MatHYPRESetPreallocation(Mat A, PetscInt dnz, const PetscInt dnnz[], PetscInt onz, const PetscInt onnz[]) 1690 { 1691 PetscFunctionBegin; 1692 PetscValidHeaderSpecific(A, MAT_CLASSID, 1); 1693 PetscValidType(A, 1); 1694 PetscTryMethod(A, "MatHYPRESetPreallocation_C", (Mat, PetscInt, const PetscInt[], PetscInt, const PetscInt[]), (A, dnz, dnnz, onz, onnz)); 1695 PetscFunctionReturn(PETSC_SUCCESS); 1696 } 1697 1698 /*@C 1699 MatCreateFromParCSR - Creates a `Mat` from a `hypre_ParCSRMatrix` 1700 1701 Collective 1702 1703 Input Parameters: 1704 + parcsr - the pointer to the `hypre_ParCSRMatrix` 1705 . mtype - matrix type to be created. Currently `MATAIJ`, `MATIS` and `MATHYPRE` are supported. 1706 - copymode - PETSc copying options, see `PetscCopyMode` 1707 1708 Output Parameter: 1709 . A - the matrix 1710 1711 Level: intermediate 1712 1713 .seealso: [](ch_matrices), `Mat`, `MatHYPRE`, `PetscCopyMode` 1714 @*/ 1715 PETSC_EXTERN PetscErrorCode MatCreateFromParCSR(hypre_ParCSRMatrix *parcsr, MatType mtype, PetscCopyMode copymode, Mat *A) 1716 { 1717 Mat T; 1718 Mat_HYPRE *hA; 1719 MPI_Comm comm; 1720 PetscInt rstart, rend, cstart, cend, M, N; 1721 PetscBool isseqaij, isseqaijmkl, ismpiaij, isaij, ishyp, isis; 1722 1723 PetscFunctionBegin; 1724 comm = hypre_ParCSRMatrixComm(parcsr); 1725 PetscCall(PetscStrcmp(mtype, MATSEQAIJ, &isseqaij)); 1726 PetscCall(PetscStrcmp(mtype, MATSEQAIJMKL, &isseqaijmkl)); 1727 PetscCall(PetscStrcmp(mtype, MATMPIAIJ, &ismpiaij)); 1728 PetscCall(PetscStrcmp(mtype, MATAIJ, &isaij)); 1729 PetscCall(PetscStrcmp(mtype, MATHYPRE, &ishyp)); 1730 PetscCall(PetscStrcmp(mtype, MATIS, &isis)); 1731 isaij = (PetscBool)(isseqaij || isseqaijmkl || ismpiaij || isaij); 1732 /* TODO */ 1733 PetscCheck(isaij || ishyp || isis, comm, PETSC_ERR_SUP, "Unsupported MatType %s! Supported types are %s, %s, %s, %s, %s, and %s", mtype, MATAIJ, MATSEQAIJ, MATSEQAIJMKL, MATMPIAIJ, MATIS, MATHYPRE); 1734 /* access ParCSRMatrix */ 1735 rstart = hypre_ParCSRMatrixFirstRowIndex(parcsr); 1736 rend = hypre_ParCSRMatrixLastRowIndex(parcsr); 1737 cstart = hypre_ParCSRMatrixFirstColDiag(parcsr); 1738 cend = hypre_ParCSRMatrixLastColDiag(parcsr); 1739 M = hypre_ParCSRMatrixGlobalNumRows(parcsr); 1740 N = hypre_ParCSRMatrixGlobalNumCols(parcsr); 1741 1742 /* fix for empty local rows/columns */ 1743 if (rend < rstart) rend = rstart; 1744 if (cend < cstart) cend = cstart; 1745 1746 /* PETSc convention */ 1747 rend++; 1748 cend++; 1749 rend = PetscMin(rend, M); 1750 cend = PetscMin(cend, N); 1751 1752 /* create PETSc matrix with MatHYPRE */ 1753 PetscCall(MatCreate(comm, &T)); 1754 PetscCall(MatSetSizes(T, rend - rstart, cend - cstart, M, N)); 1755 PetscCall(MatSetType(T, MATHYPRE)); 1756 hA = (Mat_HYPRE *)(T->data); 1757 1758 /* create HYPRE_IJMatrix */ 1759 PetscCallExternal(HYPRE_IJMatrixCreate, hA->comm, rstart, rend - 1, cstart, cend - 1, &hA->ij); 1760 PetscCallExternal(HYPRE_IJMatrixSetObjectType, hA->ij, HYPRE_PARCSR); 1761 1762 /* create new ParCSR object if needed */ 1763 if (ishyp && copymode == PETSC_COPY_VALUES) { 1764 hypre_ParCSRMatrix *new_parcsr; 1765 #if PETSC_PKG_HYPRE_VERSION_LT(2, 18, 0) 1766 hypre_CSRMatrix *hdiag, *hoffd, *ndiag, *noffd; 1767 1768 new_parcsr = hypre_ParCSRMatrixClone(parcsr, 0); 1769 hdiag = hypre_ParCSRMatrixDiag(parcsr); 1770 hoffd = hypre_ParCSRMatrixOffd(parcsr); 1771 ndiag = hypre_ParCSRMatrixDiag(new_parcsr); 1772 noffd = hypre_ParCSRMatrixOffd(new_parcsr); 1773 PetscCall(PetscArraycpy(hypre_CSRMatrixData(ndiag), hypre_CSRMatrixData(hdiag), hypre_CSRMatrixNumNonzeros(hdiag))); 1774 PetscCall(PetscArraycpy(hypre_CSRMatrixData(noffd), hypre_CSRMatrixData(hoffd), hypre_CSRMatrixNumNonzeros(hoffd))); 1775 #else 1776 new_parcsr = hypre_ParCSRMatrixClone(parcsr, 1); 1777 #endif 1778 parcsr = new_parcsr; 1779 copymode = PETSC_OWN_POINTER; 1780 } 1781 1782 /* set ParCSR object */ 1783 hypre_IJMatrixObject(hA->ij) = parcsr; 1784 T->preallocated = PETSC_TRUE; 1785 1786 /* set assembled flag */ 1787 hypre_IJMatrixAssembleFlag(hA->ij) = 1; 1788 #if 0 1789 PetscCallExternal(HYPRE_IJMatrixInitialize,hA->ij); 1790 #endif 1791 if (ishyp) { 1792 PetscMPIInt myid = 0; 1793 1794 /* make sure we always have row_starts and col_starts available */ 1795 if (HYPRE_AssumedPartitionCheck()) PetscCallMPI(MPI_Comm_rank(comm, &myid)); 1796 #if defined(hypre_ParCSRMatrixOwnsRowStarts) 1797 if (!hypre_ParCSRMatrixOwnsColStarts(parcsr)) { 1798 PetscLayout map; 1799 1800 PetscCall(MatGetLayouts(T, NULL, &map)); 1801 PetscCall(PetscLayoutSetUp(map)); 1802 hypre_ParCSRMatrixColStarts(parcsr) = (HYPRE_BigInt *)(map->range + myid); 1803 } 1804 if (!hypre_ParCSRMatrixOwnsRowStarts(parcsr)) { 1805 PetscLayout map; 1806 1807 PetscCall(MatGetLayouts(T, &map, NULL)); 1808 PetscCall(PetscLayoutSetUp(map)); 1809 hypre_ParCSRMatrixRowStarts(parcsr) = (HYPRE_BigInt *)(map->range + myid); 1810 } 1811 #endif 1812 /* prevent from freeing the pointer */ 1813 if (copymode == PETSC_USE_POINTER) hA->inner_free = PETSC_FALSE; 1814 *A = T; 1815 PetscCall(MatSetOption(*A, MAT_SORTED_FULL, PETSC_TRUE)); 1816 PetscCall(MatAssemblyBegin(*A, MAT_FINAL_ASSEMBLY)); 1817 PetscCall(MatAssemblyEnd(*A, MAT_FINAL_ASSEMBLY)); 1818 } else if (isaij) { 1819 if (copymode != PETSC_OWN_POINTER) { 1820 /* prevent from freeing the pointer */ 1821 hA->inner_free = PETSC_FALSE; 1822 PetscCall(MatConvert_HYPRE_AIJ(T, MATAIJ, MAT_INITIAL_MATRIX, A)); 1823 PetscCall(MatDestroy(&T)); 1824 } else { /* AIJ return type with PETSC_OWN_POINTER */ 1825 PetscCall(MatConvert_HYPRE_AIJ(T, MATAIJ, MAT_INPLACE_MATRIX, &T)); 1826 *A = T; 1827 } 1828 } else if (isis) { 1829 PetscCall(MatConvert_HYPRE_IS(T, MATIS, MAT_INITIAL_MATRIX, A)); 1830 if (copymode != PETSC_OWN_POINTER) hA->inner_free = PETSC_FALSE; 1831 PetscCall(MatDestroy(&T)); 1832 } 1833 PetscFunctionReturn(PETSC_SUCCESS); 1834 } 1835 1836 static PetscErrorCode MatHYPREGetParCSR_HYPRE(Mat A, hypre_ParCSRMatrix **parcsr) 1837 { 1838 Mat_HYPRE *hA = (Mat_HYPRE *)A->data; 1839 HYPRE_Int type; 1840 1841 PetscFunctionBegin; 1842 PetscCheck(hA->ij, PetscObjectComm((PetscObject)A), PETSC_ERR_PLIB, "HYPRE_IJMatrix not present"); 1843 PetscCallExternal(HYPRE_IJMatrixGetObjectType, hA->ij, &type); 1844 PetscCheck(type == HYPRE_PARCSR, PetscObjectComm((PetscObject)A), PETSC_ERR_SUP, "HYPRE_IJMatrix is not of type HYPRE_PARCSR"); 1845 PetscCallExternal(HYPRE_IJMatrixGetObject, hA->ij, (void **)parcsr); 1846 PetscFunctionReturn(PETSC_SUCCESS); 1847 } 1848 1849 /*@C 1850 MatHYPREGetParCSR - Gets the pointer to the ParCSR matrix 1851 1852 Not Collective 1853 1854 Input Parameter: 1855 . A - the `MATHYPRE` object 1856 1857 Output Parameter: 1858 . parcsr - the pointer to the `hypre_ParCSRMatrix` 1859 1860 Level: intermediate 1861 1862 .seealso: [](ch_matrices), `Mat`, `MatHYPRE`, `PetscCopyMode` 1863 @*/ 1864 PetscErrorCode MatHYPREGetParCSR(Mat A, hypre_ParCSRMatrix **parcsr) 1865 { 1866 PetscFunctionBegin; 1867 PetscValidHeaderSpecific(A, MAT_CLASSID, 1); 1868 PetscValidType(A, 1); 1869 PetscUseMethod(A, "MatHYPREGetParCSR_C", (Mat, hypre_ParCSRMatrix **), (A, parcsr)); 1870 PetscFunctionReturn(PETSC_SUCCESS); 1871 } 1872 1873 static PetscErrorCode MatMissingDiagonal_HYPRE(Mat A, PetscBool *missing, PetscInt *dd) 1874 { 1875 hypre_ParCSRMatrix *parcsr; 1876 hypre_CSRMatrix *ha; 1877 PetscInt rst; 1878 1879 PetscFunctionBegin; 1880 PetscCheck(A->rmap->n == A->cmap->n, PETSC_COMM_SELF, PETSC_ERR_SUP, "Not implemented with non-square diagonal blocks"); 1881 PetscCall(MatGetOwnershipRange(A, &rst, NULL)); 1882 PetscCall(MatHYPREGetParCSR_HYPRE(A, &parcsr)); 1883 if (missing) *missing = PETSC_FALSE; 1884 if (dd) *dd = -1; 1885 ha = hypre_ParCSRMatrixDiag(parcsr); 1886 if (ha) { 1887 PetscInt size, i; 1888 HYPRE_Int *ii, *jj; 1889 1890 size = hypre_CSRMatrixNumRows(ha); 1891 ii = hypre_CSRMatrixI(ha); 1892 jj = hypre_CSRMatrixJ(ha); 1893 for (i = 0; i < size; i++) { 1894 PetscInt j; 1895 PetscBool found = PETSC_FALSE; 1896 1897 for (j = ii[i]; j < ii[i + 1] && !found; j++) found = (jj[j] == i) ? PETSC_TRUE : PETSC_FALSE; 1898 1899 if (!found) { 1900 PetscCall(PetscInfo(A, "Matrix is missing local diagonal entry %" PetscInt_FMT "\n", i)); 1901 if (missing) *missing = PETSC_TRUE; 1902 if (dd) *dd = i + rst; 1903 PetscFunctionReturn(PETSC_SUCCESS); 1904 } 1905 } 1906 if (!size) { 1907 PetscCall(PetscInfo(A, "Matrix has no diagonal entries therefore is missing diagonal\n")); 1908 if (missing) *missing = PETSC_TRUE; 1909 if (dd) *dd = rst; 1910 } 1911 } else { 1912 PetscCall(PetscInfo(A, "Matrix has no diagonal entries therefore is missing diagonal\n")); 1913 if (missing) *missing = PETSC_TRUE; 1914 if (dd) *dd = rst; 1915 } 1916 PetscFunctionReturn(PETSC_SUCCESS); 1917 } 1918 1919 static PetscErrorCode MatScale_HYPRE(Mat A, PetscScalar s) 1920 { 1921 hypre_ParCSRMatrix *parcsr; 1922 #if PETSC_PKG_HYPRE_VERSION_LT(2, 19, 0) 1923 hypre_CSRMatrix *ha; 1924 #endif 1925 HYPRE_Complex hs; 1926 1927 PetscFunctionBegin; 1928 PetscCall(PetscHYPREScalarCast(s, &hs)); 1929 PetscCall(MatHYPREGetParCSR_HYPRE(A, &parcsr)); 1930 #if PETSC_PKG_HYPRE_VERSION_GE(2, 19, 0) 1931 PetscCallExternal(hypre_ParCSRMatrixScale, parcsr, hs); 1932 #else /* diagonal part */ 1933 ha = hypre_ParCSRMatrixDiag(parcsr); 1934 if (ha) { 1935 PetscInt size, i; 1936 HYPRE_Int *ii; 1937 HYPRE_Complex *a; 1938 1939 size = hypre_CSRMatrixNumRows(ha); 1940 a = hypre_CSRMatrixData(ha); 1941 ii = hypre_CSRMatrixI(ha); 1942 for (i = 0; i < ii[size]; i++) a[i] *= hs; 1943 } 1944 /* off-diagonal part */ 1945 ha = hypre_ParCSRMatrixOffd(parcsr); 1946 if (ha) { 1947 PetscInt size, i; 1948 HYPRE_Int *ii; 1949 HYPRE_Complex *a; 1950 1951 size = hypre_CSRMatrixNumRows(ha); 1952 a = hypre_CSRMatrixData(ha); 1953 ii = hypre_CSRMatrixI(ha); 1954 for (i = 0; i < ii[size]; i++) a[i] *= hs; 1955 } 1956 #endif 1957 PetscFunctionReturn(PETSC_SUCCESS); 1958 } 1959 1960 static PetscErrorCode MatZeroRowsColumns_HYPRE(Mat A, PetscInt numRows, const PetscInt rows[], PetscScalar diag, Vec x, Vec b) 1961 { 1962 hypre_ParCSRMatrix *parcsr; 1963 HYPRE_Int *lrows; 1964 PetscInt rst, ren, i; 1965 1966 PetscFunctionBegin; 1967 PetscCheck(!x && !b, PetscObjectComm((PetscObject)A), PETSC_ERR_SUP, "To be implemented"); 1968 PetscCall(MatHYPREGetParCSR_HYPRE(A, &parcsr)); 1969 PetscCall(PetscMalloc1(numRows, &lrows)); 1970 PetscCall(MatGetOwnershipRange(A, &rst, &ren)); 1971 for (i = 0; i < numRows; i++) { 1972 PetscCheck(rows[i] >= rst && rows[i] < ren, PETSC_COMM_SELF, PETSC_ERR_SUP, "Non-local rows not yet supported"); 1973 lrows[i] = rows[i] - rst; 1974 } 1975 PetscCallExternal(hypre_ParCSRMatrixEliminateRowsCols, parcsr, numRows, lrows); 1976 PetscCall(PetscFree(lrows)); 1977 PetscFunctionReturn(PETSC_SUCCESS); 1978 } 1979 1980 static PetscErrorCode MatZeroEntries_HYPRE_CSRMatrix(hypre_CSRMatrix *ha) 1981 { 1982 PetscFunctionBegin; 1983 if (ha) { 1984 HYPRE_Int *ii, size; 1985 HYPRE_Complex *a; 1986 1987 size = hypre_CSRMatrixNumRows(ha); 1988 a = hypre_CSRMatrixData(ha); 1989 ii = hypre_CSRMatrixI(ha); 1990 1991 if (a) PetscCall(PetscArrayzero(a, ii[size])); 1992 } 1993 PetscFunctionReturn(PETSC_SUCCESS); 1994 } 1995 1996 static PetscErrorCode MatZeroEntries_HYPRE(Mat A) 1997 { 1998 Mat_HYPRE *hA = (Mat_HYPRE *)A->data; 1999 2000 PetscFunctionBegin; 2001 if (HYPRE_MEMORY_DEVICE == hypre_IJMatrixMemoryLocation(hA->ij)) { 2002 PetscCallExternal(HYPRE_IJMatrixSetConstantValues, hA->ij, 0.0); 2003 } else { 2004 hypre_ParCSRMatrix *parcsr; 2005 2006 PetscCall(MatHYPREGetParCSR_HYPRE(A, &parcsr)); 2007 PetscCall(MatZeroEntries_HYPRE_CSRMatrix(hypre_ParCSRMatrixDiag(parcsr))); 2008 PetscCall(MatZeroEntries_HYPRE_CSRMatrix(hypre_ParCSRMatrixOffd(parcsr))); 2009 } 2010 PetscFunctionReturn(PETSC_SUCCESS); 2011 } 2012 2013 static PetscErrorCode MatZeroRows_HYPRE_CSRMatrix(hypre_CSRMatrix *hA, PetscInt N, const PetscInt rows[], HYPRE_Complex diag) 2014 { 2015 PetscInt ii; 2016 HYPRE_Int *i, *j; 2017 HYPRE_Complex *a; 2018 2019 PetscFunctionBegin; 2020 if (!hA) PetscFunctionReturn(PETSC_SUCCESS); 2021 2022 i = hypre_CSRMatrixI(hA); 2023 j = hypre_CSRMatrixJ(hA); 2024 a = hypre_CSRMatrixData(hA); 2025 #if defined(PETSC_HAVE_HYPRE_DEVICE) 2026 if (HYPRE_MEMORY_DEVICE == hypre_CSRMatrixMemoryLocation(hA)) { 2027 #if defined(HYPRE_USING_CUDA) 2028 MatZeroRows_CUDA(N, rows, i, j, a, diag); 2029 #elif defined(HYPRE_USING_HIP) 2030 MatZeroRows_HIP(N, rows, i, j, a, diag); 2031 #elif defined(PETSC_HAVE_KOKKOS) 2032 MatZeroRows_Kokkos(N, rows, i, j, a, diag); 2033 #else 2034 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "No support for MatZeroRows on a hypre matrix in this memory location"); 2035 #endif 2036 } else 2037 #endif 2038 { 2039 for (ii = 0; ii < N; ii++) { 2040 HYPRE_Int jj, ibeg, iend, irow; 2041 2042 irow = rows[ii]; 2043 ibeg = i[irow]; 2044 iend = i[irow + 1]; 2045 for (jj = ibeg; jj < iend; jj++) 2046 if (j[jj] == irow) a[jj] = diag; 2047 else a[jj] = 0.0; 2048 } 2049 } 2050 PetscFunctionReturn(PETSC_SUCCESS); 2051 } 2052 2053 static PetscErrorCode MatZeroRows_HYPRE(Mat A, PetscInt N, const PetscInt rows[], PetscScalar diag, Vec x, Vec b) 2054 { 2055 hypre_ParCSRMatrix *parcsr; 2056 PetscInt *lrows, len, *lrows2; 2057 HYPRE_Complex hdiag; 2058 2059 PetscFunctionBegin; 2060 PetscCheck(!x && !b, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Does not support to modify the solution and the right hand size"); 2061 PetscCall(PetscHYPREScalarCast(diag, &hdiag)); 2062 /* retrieve the internal matrix */ 2063 PetscCall(MatHYPREGetParCSR_HYPRE(A, &parcsr)); 2064 /* get locally owned rows */ 2065 PetscCall(MatZeroRowsMapLocal_Private(A, N, rows, &len, &lrows)); 2066 2067 #if defined(PETSC_HAVE_HYPRE_DEVICE) 2068 if (HYPRE_MEMORY_DEVICE == hypre_CSRMatrixMemoryLocation(hypre_ParCSRMatrixDiag(parcsr))) { 2069 Mat_HYPRE *hA = (Mat_HYPRE *)A->data; 2070 PetscInt m; 2071 PetscCall(MatGetLocalSize(A, &m, NULL)); 2072 if (!hA->rows_d) { 2073 hA->rows_d = hypre_TAlloc(PetscInt, m, HYPRE_MEMORY_DEVICE); 2074 if (m) PetscCheck(hA->rows_d, PETSC_COMM_SELF, PETSC_ERR_MEM, "HYPRE_TAlloc failed"); 2075 } 2076 PetscCheck(len <= m, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Too many rows in rows[]"); 2077 PetscStackCallExternalVoid("hypre_Memcpy", hypre_Memcpy(hA->rows_d, lrows, sizeof(PetscInt) * len, HYPRE_MEMORY_DEVICE, HYPRE_MEMORY_HOST)); 2078 lrows2 = hA->rows_d; 2079 } else 2080 #endif 2081 { 2082 lrows2 = lrows; 2083 } 2084 2085 /* zero diagonal part */ 2086 PetscCall(MatZeroRows_HYPRE_CSRMatrix(hypre_ParCSRMatrixDiag(parcsr), len, lrows2, hdiag)); 2087 /* zero off-diagonal part */ 2088 PetscCall(MatZeroRows_HYPRE_CSRMatrix(hypre_ParCSRMatrixOffd(parcsr), len, lrows2, 0.0)); 2089 2090 PetscCall(PetscFree(lrows)); 2091 PetscFunctionReturn(PETSC_SUCCESS); 2092 } 2093 2094 static PetscErrorCode MatAssemblyBegin_HYPRE(Mat mat, MatAssemblyType mode) 2095 { 2096 PetscFunctionBegin; 2097 if (mat->nooffprocentries) PetscFunctionReturn(PETSC_SUCCESS); 2098 2099 PetscCall(MatStashScatterBegin_Private(mat, &mat->stash, mat->rmap->range)); 2100 PetscFunctionReturn(PETSC_SUCCESS); 2101 } 2102 2103 static PetscErrorCode MatGetRow_HYPRE(Mat A, PetscInt row, PetscInt *nz, PetscInt **idx, PetscScalar **v) 2104 { 2105 hypre_ParCSRMatrix *parcsr; 2106 HYPRE_Int hnz; 2107 2108 PetscFunctionBegin; 2109 /* retrieve the internal matrix */ 2110 PetscCall(MatHYPREGetParCSR_HYPRE(A, &parcsr)); 2111 /* call HYPRE API */ 2112 PetscCallExternal(HYPRE_ParCSRMatrixGetRow, parcsr, row, &hnz, (HYPRE_BigInt **)idx, (HYPRE_Complex **)v); 2113 if (nz) *nz = (PetscInt)hnz; 2114 PetscFunctionReturn(PETSC_SUCCESS); 2115 } 2116 2117 static PetscErrorCode MatRestoreRow_HYPRE(Mat A, PetscInt row, PetscInt *nz, PetscInt **idx, PetscScalar **v) 2118 { 2119 hypre_ParCSRMatrix *parcsr; 2120 HYPRE_Int hnz; 2121 2122 PetscFunctionBegin; 2123 /* retrieve the internal matrix */ 2124 PetscCall(MatHYPREGetParCSR_HYPRE(A, &parcsr)); 2125 /* call HYPRE API */ 2126 hnz = nz ? (HYPRE_Int)(*nz) : 0; 2127 PetscCallExternal(HYPRE_ParCSRMatrixRestoreRow, parcsr, row, &hnz, (HYPRE_BigInt **)idx, (HYPRE_Complex **)v); 2128 PetscFunctionReturn(PETSC_SUCCESS); 2129 } 2130 2131 static PetscErrorCode MatGetValues_HYPRE(Mat A, PetscInt m, const PetscInt idxm[], PetscInt n, const PetscInt idxn[], PetscScalar v[]) 2132 { 2133 Mat_HYPRE *hA = (Mat_HYPRE *)A->data; 2134 PetscInt i; 2135 2136 PetscFunctionBegin; 2137 if (!m || !n) PetscFunctionReturn(PETSC_SUCCESS); 2138 /* Ignore negative row indices 2139 * And negative column indices should be automatically ignored in hypre 2140 * */ 2141 for (i = 0; i < m; i++) { 2142 if (idxm[i] >= 0) { 2143 HYPRE_Int hn = (HYPRE_Int)n; 2144 PetscCallExternal(HYPRE_IJMatrixGetValues, hA->ij, 1, &hn, (HYPRE_BigInt *)&idxm[i], (HYPRE_BigInt *)idxn, (HYPRE_Complex *)(v + i * n)); 2145 } 2146 } 2147 PetscFunctionReturn(PETSC_SUCCESS); 2148 } 2149 2150 static PetscErrorCode MatSetOption_HYPRE(Mat A, MatOption op, PetscBool flg) 2151 { 2152 Mat_HYPRE *hA = (Mat_HYPRE *)A->data; 2153 2154 PetscFunctionBegin; 2155 switch (op) { 2156 case MAT_NO_OFF_PROC_ENTRIES: 2157 if (flg) PetscCallExternal(HYPRE_IJMatrixSetMaxOffProcElmts, hA->ij, 0); 2158 break; 2159 case MAT_IGNORE_OFF_PROC_ENTRIES: 2160 hA->donotstash = flg; 2161 break; 2162 default: 2163 break; 2164 } 2165 PetscFunctionReturn(PETSC_SUCCESS); 2166 } 2167 2168 static PetscErrorCode MatView_HYPRE(Mat A, PetscViewer view) 2169 { 2170 PetscViewerFormat format; 2171 2172 PetscFunctionBegin; 2173 PetscCall(PetscViewerGetFormat(view, &format)); 2174 if (format == PETSC_VIEWER_ASCII_FACTOR_INFO || format == PETSC_VIEWER_ASCII_INFO || format == PETSC_VIEWER_ASCII_INFO_DETAIL) PetscFunctionReturn(PETSC_SUCCESS); 2175 if (format != PETSC_VIEWER_NATIVE) { 2176 Mat B; 2177 hypre_ParCSRMatrix *parcsr; 2178 PetscErrorCode (*mview)(Mat, PetscViewer) = NULL; 2179 2180 PetscCall(MatHYPREGetParCSR_HYPRE(A, &parcsr)); 2181 PetscCall(MatCreateFromParCSR(parcsr, MATAIJ, PETSC_USE_POINTER, &B)); 2182 PetscCall(MatGetOperation(B, MATOP_VIEW, (void (**)(void)) & mview)); 2183 PetscCheck(mview, PetscObjectComm((PetscObject)A), PETSC_ERR_PLIB, "Missing view operation"); 2184 PetscCall((*mview)(B, view)); 2185 PetscCall(MatDestroy(&B)); 2186 } else { 2187 Mat_HYPRE *hA = (Mat_HYPRE *)A->data; 2188 PetscMPIInt size; 2189 PetscBool isascii; 2190 const char *filename; 2191 2192 /* HYPRE uses only text files */ 2193 PetscCall(PetscObjectTypeCompare((PetscObject)view, PETSCVIEWERASCII, &isascii)); 2194 PetscCheck(isascii, PetscObjectComm((PetscObject)view), PETSC_ERR_SUP, "PetscViewerType %s: native HYPRE format needs PETSCVIEWERASCII", ((PetscObject)view)->type_name); 2195 PetscCall(PetscViewerFileGetName(view, &filename)); 2196 PetscCallExternal(HYPRE_IJMatrixPrint, hA->ij, filename); 2197 PetscCallMPI(MPI_Comm_size(hA->comm, &size)); 2198 if (size > 1) { 2199 PetscCall(PetscViewerASCIIPrintf(view, "Matrix files: %s.%05d ... %s.%05d\n", filename, 0, filename, size - 1)); 2200 } else { 2201 PetscCall(PetscViewerASCIIPrintf(view, "Matrix file: %s.%05d\n", filename, 0)); 2202 } 2203 } 2204 PetscFunctionReturn(PETSC_SUCCESS); 2205 } 2206 2207 static PetscErrorCode MatCopy_HYPRE(Mat A, Mat B, MatStructure str) 2208 { 2209 hypre_ParCSRMatrix *acsr, *bcsr; 2210 2211 PetscFunctionBegin; 2212 if (str == SAME_NONZERO_PATTERN && A->ops->copy == B->ops->copy) { 2213 PetscCall(MatHYPREGetParCSR_HYPRE(A, &acsr)); 2214 PetscCall(MatHYPREGetParCSR_HYPRE(B, &bcsr)); 2215 PetscCallExternal(hypre_ParCSRMatrixCopy, acsr, bcsr, 1); 2216 PetscCall(MatSetOption(B, MAT_SORTED_FULL, PETSC_TRUE)); /* "perfect" preallocation, so no need for hypre_AuxParCSRMatrixNeedAux */ 2217 PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY)); 2218 PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY)); 2219 } else { 2220 PetscCall(MatCopy_Basic(A, B, str)); 2221 } 2222 PetscFunctionReturn(PETSC_SUCCESS); 2223 } 2224 2225 static PetscErrorCode MatGetDiagonal_HYPRE(Mat A, Vec d) 2226 { 2227 hypre_ParCSRMatrix *parcsr; 2228 hypre_CSRMatrix *dmat; 2229 HYPRE_Complex *a; 2230 PetscBool cong; 2231 2232 PetscFunctionBegin; 2233 PetscCall(MatHasCongruentLayouts(A, &cong)); 2234 PetscCheck(cong, PetscObjectComm((PetscObject)A), PETSC_ERR_SUP, "Only for square matrices with same local distributions of rows and columns"); 2235 PetscCall(MatHYPREGetParCSR_HYPRE(A, &parcsr)); 2236 dmat = hypre_ParCSRMatrixDiag(parcsr); 2237 if (dmat) { 2238 #if defined(PETSC_HAVE_HYPRE_DEVICE) 2239 HYPRE_MemoryLocation mem = hypre_CSRMatrixMemoryLocation(dmat); 2240 #else 2241 HYPRE_MemoryLocation mem = HYPRE_MEMORY_HOST; 2242 #endif 2243 2244 if (mem != HYPRE_MEMORY_HOST) PetscCall(VecGetArrayWriteAndMemType(d, (PetscScalar **)&a, NULL)); 2245 else PetscCall(VecGetArrayWrite(d, (PetscScalar **)&a)); 2246 hypre_CSRMatrixExtractDiagonal(dmat, a, 0); 2247 if (mem != HYPRE_MEMORY_HOST) PetscCall(VecRestoreArrayWriteAndMemType(d, (PetscScalar **)&a)); 2248 else PetscCall(VecRestoreArrayWrite(d, (PetscScalar **)&a)); 2249 } 2250 PetscFunctionReturn(PETSC_SUCCESS); 2251 } 2252 2253 #include <petscblaslapack.h> 2254 2255 static PetscErrorCode MatAXPY_HYPRE(Mat Y, PetscScalar a, Mat X, MatStructure str) 2256 { 2257 PetscFunctionBegin; 2258 #if defined(PETSC_HAVE_HYPRE_DEVICE) 2259 { 2260 Mat B; 2261 hypre_ParCSRMatrix *x, *y, *z; 2262 2263 PetscCall(MatHYPREGetParCSR(Y, &y)); 2264 PetscCall(MatHYPREGetParCSR(X, &x)); 2265 PetscCallExternal(hypre_ParCSRMatrixAdd, 1.0, y, 1.0, x, &z); 2266 PetscCall(MatCreateFromParCSR(z, MATHYPRE, PETSC_OWN_POINTER, &B)); 2267 PetscCall(MatHeaderMerge(Y, &B)); 2268 } 2269 #else 2270 if (str == SAME_NONZERO_PATTERN) { 2271 hypre_ParCSRMatrix *x, *y; 2272 hypre_CSRMatrix *xloc, *yloc; 2273 PetscInt xnnz, ynnz; 2274 HYPRE_Complex *xarr, *yarr; 2275 PetscBLASInt one = 1, bnz; 2276 2277 PetscCall(MatHYPREGetParCSR(Y, &y)); 2278 PetscCall(MatHYPREGetParCSR(X, &x)); 2279 2280 /* diagonal block */ 2281 xloc = hypre_ParCSRMatrixDiag(x); 2282 yloc = hypre_ParCSRMatrixDiag(y); 2283 xnnz = 0; 2284 ynnz = 0; 2285 xarr = NULL; 2286 yarr = NULL; 2287 if (xloc) { 2288 xarr = hypre_CSRMatrixData(xloc); 2289 xnnz = hypre_CSRMatrixNumNonzeros(xloc); 2290 } 2291 if (yloc) { 2292 yarr = hypre_CSRMatrixData(yloc); 2293 ynnz = hypre_CSRMatrixNumNonzeros(yloc); 2294 } 2295 PetscCheck(xnnz == ynnz, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Different number of nonzeros in diagonal block %" PetscInt_FMT " != %" PetscInt_FMT, xnnz, ynnz); 2296 PetscCall(PetscBLASIntCast(xnnz, &bnz)); 2297 PetscCallBLAS("BLASaxpy", BLASaxpy_(&bnz, &a, (PetscScalar *)xarr, &one, (PetscScalar *)yarr, &one)); 2298 2299 /* off-diagonal block */ 2300 xloc = hypre_ParCSRMatrixOffd(x); 2301 yloc = hypre_ParCSRMatrixOffd(y); 2302 xnnz = 0; 2303 ynnz = 0; 2304 xarr = NULL; 2305 yarr = NULL; 2306 if (xloc) { 2307 xarr = hypre_CSRMatrixData(xloc); 2308 xnnz = hypre_CSRMatrixNumNonzeros(xloc); 2309 } 2310 if (yloc) { 2311 yarr = hypre_CSRMatrixData(yloc); 2312 ynnz = hypre_CSRMatrixNumNonzeros(yloc); 2313 } 2314 PetscCheck(xnnz == ynnz, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Different number of nonzeros in off-diagonal block %" PetscInt_FMT " != %" PetscInt_FMT, xnnz, ynnz); 2315 PetscCall(PetscBLASIntCast(xnnz, &bnz)); 2316 PetscCallBLAS("BLASaxpy", BLASaxpy_(&bnz, &a, (PetscScalar *)xarr, &one, (PetscScalar *)yarr, &one)); 2317 } else if (str == SUBSET_NONZERO_PATTERN) { 2318 PetscCall(MatAXPY_Basic(Y, a, X, str)); 2319 } else { 2320 Mat B; 2321 2322 PetscCall(MatAXPY_Basic_Preallocate(Y, X, &B)); 2323 PetscCall(MatAXPY_BasicWithPreallocation(B, Y, a, X, str)); 2324 PetscCall(MatHeaderReplace(Y, &B)); 2325 } 2326 #endif 2327 PetscFunctionReturn(PETSC_SUCCESS); 2328 } 2329 2330 static PetscErrorCode MatDuplicate_HYPRE(Mat A, MatDuplicateOption op, Mat *B) 2331 { 2332 hypre_ParCSRMatrix *parcsr = NULL; 2333 PetscCopyMode cpmode; 2334 Mat_HYPRE *hA; 2335 2336 PetscFunctionBegin; 2337 PetscCall(MatHYPREGetParCSR_HYPRE(A, &parcsr)); 2338 if (op == MAT_DO_NOT_COPY_VALUES || op == MAT_SHARE_NONZERO_PATTERN) { 2339 parcsr = hypre_ParCSRMatrixClone(parcsr, 0); 2340 cpmode = PETSC_OWN_POINTER; 2341 } else { 2342 cpmode = PETSC_COPY_VALUES; 2343 } 2344 PetscCall(MatCreateFromParCSR(parcsr, MATHYPRE, cpmode, B)); 2345 hA = (Mat_HYPRE *)A->data; 2346 if (hA->cooMat) { 2347 Mat_HYPRE *hB = (Mat_HYPRE *)((*B)->data); 2348 op = (op == MAT_DO_NOT_COPY_VALUES) ? op : MAT_COPY_VALUES; 2349 /* Cannot simply increase the reference count of hA->cooMat, since B needs to share cooMat's data array */ 2350 PetscCall(MatDuplicate(hA->cooMat, op, &hB->cooMat)); 2351 PetscCall(MatHYPRE_AttachCOOMat(*B)); 2352 } 2353 PetscFunctionReturn(PETSC_SUCCESS); 2354 } 2355 2356 static PetscErrorCode MatSetPreallocationCOO_HYPRE(Mat mat, PetscCount coo_n, PetscInt coo_i[], PetscInt coo_j[]) 2357 { 2358 Mat_HYPRE *hmat = (Mat_HYPRE *)mat->data; 2359 2360 PetscFunctionBegin; 2361 /* Build an agent matrix cooMat with AIJ format 2362 It has the same sparsity pattern as mat, and also shares the data array with mat. We use cooMat to do the COO work. 2363 */ 2364 PetscCall(MatHYPRE_CreateCOOMat(mat)); 2365 PetscCall(MatSetOption(hmat->cooMat, MAT_IGNORE_OFF_PROC_ENTRIES, hmat->donotstash)); 2366 PetscCall(MatSetOption(hmat->cooMat, MAT_NO_OFF_PROC_ENTRIES, mat->nooffprocentries)); 2367 2368 /* MatSetPreallocationCOO_SeqAIJ and MatSetPreallocationCOO_MPIAIJ uses this specific 2369 name to automatically put the diagonal entries first */ 2370 PetscCall(PetscObjectSetName((PetscObject)hmat->cooMat, "_internal_COO_mat_for_hypre")); 2371 PetscCall(MatSetPreallocationCOO(hmat->cooMat, coo_n, coo_i, coo_j)); 2372 hmat->cooMat->assembled = PETSC_TRUE; 2373 2374 /* Copy the sparsity pattern from cooMat to hypre IJMatrix hmat->ij */ 2375 PetscCall(MatSetOption(mat, MAT_SORTED_FULL, PETSC_TRUE)); 2376 PetscCall(MatHYPRE_CreateFromMat(hmat->cooMat, hmat)); /* Create hmat->ij and preallocate it */ 2377 PetscCall(MatHYPRE_IJMatrixCopyIJ(hmat->cooMat, hmat->ij)); /* Copy A's (i,j) to hmat->ij */ 2378 2379 mat->preallocated = PETSC_TRUE; 2380 PetscCall(MatAssemblyBegin(mat, MAT_FINAL_ASSEMBLY)); 2381 PetscCall(MatAssemblyEnd(mat, MAT_FINAL_ASSEMBLY)); /* Migrate mat to device if it is bound to. Hypre builds its own SpMV context here */ 2382 2383 /* Attach cooMat to mat */ 2384 PetscCall(MatHYPRE_AttachCOOMat(mat)); 2385 PetscFunctionReturn(PETSC_SUCCESS); 2386 } 2387 2388 static PetscErrorCode MatSetValuesCOO_HYPRE(Mat mat, const PetscScalar v[], InsertMode imode) 2389 { 2390 Mat_HYPRE *hmat = (Mat_HYPRE *)mat->data; 2391 2392 PetscFunctionBegin; 2393 PetscCheck(hmat->cooMat, PetscObjectComm((PetscObject)mat), PETSC_ERR_PLIB, "HYPRE COO delegate matrix has not been created yet"); 2394 PetscCall(MatSetValuesCOO(hmat->cooMat, v, imode)); 2395 PetscCall(MatViewFromOptions(hmat->cooMat, (PetscObject)mat, "-cooMat_view")); 2396 PetscFunctionReturn(PETSC_SUCCESS); 2397 } 2398 2399 /*MC 2400 MATHYPRE - "hypre" - A matrix type to be used for sequential and parallel sparse matrices 2401 based on the hypre IJ interface. 2402 2403 Level: intermediate 2404 2405 .seealso: [](ch_matrices), `Mat`, `MatCreate()`, `MatHYPRESetPreallocation` 2406 M*/ 2407 2408 PETSC_EXTERN PetscErrorCode MatCreate_HYPRE(Mat B) 2409 { 2410 Mat_HYPRE *hB; 2411 2412 PetscFunctionBegin; 2413 PetscCall(PetscNew(&hB)); 2414 2415 hB->inner_free = PETSC_TRUE; 2416 hB->array_available = PETSC_TRUE; 2417 2418 B->data = (void *)hB; 2419 2420 PetscCall(PetscMemzero(B->ops, sizeof(struct _MatOps))); 2421 B->ops->mult = MatMult_HYPRE; 2422 B->ops->multtranspose = MatMultTranspose_HYPRE; 2423 B->ops->multadd = MatMultAdd_HYPRE; 2424 B->ops->multtransposeadd = MatMultTransposeAdd_HYPRE; 2425 B->ops->setup = MatSetUp_HYPRE; 2426 B->ops->destroy = MatDestroy_HYPRE; 2427 B->ops->assemblyend = MatAssemblyEnd_HYPRE; 2428 B->ops->assemblybegin = MatAssemblyBegin_HYPRE; 2429 B->ops->setvalues = MatSetValues_HYPRE; 2430 B->ops->missingdiagonal = MatMissingDiagonal_HYPRE; 2431 B->ops->scale = MatScale_HYPRE; 2432 B->ops->zerorowscolumns = MatZeroRowsColumns_HYPRE; 2433 B->ops->zeroentries = MatZeroEntries_HYPRE; 2434 B->ops->zerorows = MatZeroRows_HYPRE; 2435 B->ops->getrow = MatGetRow_HYPRE; 2436 B->ops->restorerow = MatRestoreRow_HYPRE; 2437 B->ops->getvalues = MatGetValues_HYPRE; 2438 B->ops->setoption = MatSetOption_HYPRE; 2439 B->ops->duplicate = MatDuplicate_HYPRE; 2440 B->ops->copy = MatCopy_HYPRE; 2441 B->ops->view = MatView_HYPRE; 2442 B->ops->getdiagonal = MatGetDiagonal_HYPRE; 2443 B->ops->axpy = MatAXPY_HYPRE; 2444 B->ops->productsetfromoptions = MatProductSetFromOptions_HYPRE; 2445 #if defined(PETSC_HAVE_HYPRE_DEVICE) 2446 B->ops->bindtocpu = MatBindToCPU_HYPRE; 2447 B->boundtocpu = PETSC_FALSE; 2448 #endif 2449 2450 /* build cache for off array entries formed */ 2451 PetscCall(MatStashCreate_Private(PetscObjectComm((PetscObject)B), 1, &B->stash)); 2452 2453 PetscCall(PetscCommGetComm(PetscObjectComm((PetscObject)B), &hB->comm)); 2454 PetscCall(PetscObjectChangeTypeName((PetscObject)B, MATHYPRE)); 2455 PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatConvert_hypre_aij_C", MatConvert_HYPRE_AIJ)); 2456 PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatConvert_hypre_is_C", MatConvert_HYPRE_IS)); 2457 PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatProductSetFromOptions_seqaij_hypre_C", MatProductSetFromOptions_HYPRE)); 2458 PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatProductSetFromOptions_mpiaij_hypre_C", MatProductSetFromOptions_HYPRE)); 2459 PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatHYPRESetPreallocation_C", MatHYPRESetPreallocation_HYPRE)); 2460 PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatHYPREGetParCSR_C", MatHYPREGetParCSR_HYPRE)); 2461 PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatSetPreallocationCOO_C", MatSetPreallocationCOO_HYPRE)); 2462 PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatSetValuesCOO_C", MatSetValuesCOO_HYPRE)); 2463 #if defined(PETSC_HAVE_HYPRE_DEVICE) 2464 #if defined(HYPRE_USING_HIP) 2465 PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatProductSetFromOptions_seqaijhipsparse_hypre_C", MatProductSetFromOptions_HYPRE)); 2466 PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatProductSetFromOptions_mpiaijhipsparse_hypre_C", MatProductSetFromOptions_HYPRE)); 2467 PetscCall(PetscDeviceInitialize(PETSC_DEVICE_HIP)); 2468 PetscCall(MatSetVecType(B, VECHIP)); 2469 #endif 2470 #if defined(HYPRE_USING_CUDA) 2471 PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatProductSetFromOptions_seqaijcusparse_hypre_C", MatProductSetFromOptions_HYPRE)); 2472 PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatProductSetFromOptions_mpiaijcusparse_hypre_C", MatProductSetFromOptions_HYPRE)); 2473 PetscCall(PetscDeviceInitialize(PETSC_DEVICE_CUDA)); 2474 PetscCall(MatSetVecType(B, VECCUDA)); 2475 #endif 2476 #endif 2477 PetscHYPREInitialize(); 2478 PetscFunctionReturn(PETSC_SUCCESS); 2479 } 2480