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