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