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