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