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