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