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