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 /* ---------------------------------------------------- */ 1053 static PetscErrorCode MatProductSymbolic_AB_HYPRE(Mat C) 1054 { 1055 PetscFunctionBegin; 1056 C->ops->productnumeric = MatProductNumeric_AB; 1057 PetscFunctionReturn(PETSC_SUCCESS); 1058 } 1059 1060 static PetscErrorCode MatProductSetFromOptions_HYPRE_AB(Mat C) 1061 { 1062 Mat_Product *product = C->product; 1063 PetscBool Ahypre; 1064 1065 PetscFunctionBegin; 1066 PetscCall(PetscObjectTypeCompare((PetscObject)product->A, MATHYPRE, &Ahypre)); 1067 if (Ahypre) { /* A is a Hypre matrix */ 1068 PetscCall(MatSetType(C, MATHYPRE)); 1069 C->ops->productsymbolic = MatProductSymbolic_AB_HYPRE; 1070 C->ops->matmultnumeric = MatMatMultNumeric_HYPRE_HYPRE; 1071 PetscFunctionReturn(PETSC_SUCCESS); 1072 } 1073 PetscFunctionReturn(PETSC_SUCCESS); 1074 } 1075 1076 static PetscErrorCode MatProductSymbolic_PtAP_HYPRE(Mat C) 1077 { 1078 PetscFunctionBegin; 1079 C->ops->productnumeric = MatProductNumeric_PtAP; 1080 PetscFunctionReturn(PETSC_SUCCESS); 1081 } 1082 1083 static PetscErrorCode MatProductSetFromOptions_HYPRE_PtAP(Mat C) 1084 { 1085 Mat_Product *product = C->product; 1086 PetscBool flg; 1087 PetscInt type = 0; 1088 const char *outTypes[4] = {"aij", "seqaij", "mpiaij", "hypre"}; 1089 PetscInt ntype = 4; 1090 Mat A = product->A; 1091 PetscBool Ahypre; 1092 1093 PetscFunctionBegin; 1094 PetscCall(PetscObjectTypeCompare((PetscObject)A, MATHYPRE, &Ahypre)); 1095 if (Ahypre) { /* A is a Hypre matrix */ 1096 PetscCall(MatSetType(C, MATHYPRE)); 1097 C->ops->productsymbolic = MatProductSymbolic_PtAP_HYPRE; 1098 C->ops->ptapnumeric = MatPtAPNumeric_HYPRE_HYPRE; 1099 PetscFunctionReturn(PETSC_SUCCESS); 1100 } 1101 1102 /* A is AIJ, P is Hypre, C = PtAP can be either AIJ or Hypre format */ 1103 /* Get runtime option */ 1104 if (product->api_user) { 1105 PetscOptionsBegin(PetscObjectComm((PetscObject)C), ((PetscObject)C)->prefix, "MatPtAP_HYPRE", "Mat"); 1106 PetscCall(PetscOptionsEList("-matptap_hypre_outtype", "MatPtAP outtype", "MatPtAP outtype", outTypes, ntype, outTypes[type], &type, &flg)); 1107 PetscOptionsEnd(); 1108 } else { 1109 PetscOptionsBegin(PetscObjectComm((PetscObject)C), ((PetscObject)C)->prefix, "MatProduct_PtAP_HYPRE", "Mat"); 1110 PetscCall(PetscOptionsEList("-mat_product_algorithm_hypre_outtype", "MatProduct_PtAP outtype", "MatProduct_PtAP", outTypes, ntype, outTypes[type], &type, &flg)); 1111 PetscOptionsEnd(); 1112 } 1113 1114 if (type == 0 || type == 1 || type == 2) { 1115 PetscCall(MatSetType(C, MATAIJ)); 1116 } else if (type == 3) { 1117 PetscCall(MatSetType(C, MATHYPRE)); 1118 } else SETERRQ(PetscObjectComm((PetscObject)C), PETSC_ERR_SUP, "MatPtAP outtype is not supported"); 1119 C->ops->productsymbolic = MatProductSymbolic_PtAP_HYPRE; 1120 C->ops->ptapnumeric = MatPtAPNumeric_AIJ_HYPRE; 1121 PetscFunctionReturn(PETSC_SUCCESS); 1122 } 1123 1124 static PetscErrorCode MatProductSetFromOptions_HYPRE(Mat C) 1125 { 1126 Mat_Product *product = C->product; 1127 1128 PetscFunctionBegin; 1129 switch (product->type) { 1130 case MATPRODUCT_AB: 1131 PetscCall(MatProductSetFromOptions_HYPRE_AB(C)); 1132 break; 1133 case MATPRODUCT_PtAP: 1134 PetscCall(MatProductSetFromOptions_HYPRE_PtAP(C)); 1135 break; 1136 default: 1137 break; 1138 } 1139 PetscFunctionReturn(PETSC_SUCCESS); 1140 } 1141 1142 /* -------------------------------------------------------- */ 1143 1144 static PetscErrorCode MatMultTranspose_HYPRE(Mat A, Vec x, Vec y) 1145 { 1146 PetscFunctionBegin; 1147 PetscCall(MatHYPRE_MultKernel_Private(A, 1.0, x, 0.0, y, PETSC_TRUE)); 1148 PetscFunctionReturn(PETSC_SUCCESS); 1149 } 1150 1151 static PetscErrorCode MatMult_HYPRE(Mat A, Vec x, Vec y) 1152 { 1153 PetscFunctionBegin; 1154 PetscCall(MatHYPRE_MultKernel_Private(A, 1.0, x, 0.0, y, PETSC_FALSE)); 1155 PetscFunctionReturn(PETSC_SUCCESS); 1156 } 1157 1158 static PetscErrorCode MatMultAdd_HYPRE(Mat A, Vec x, Vec y, Vec z) 1159 { 1160 PetscFunctionBegin; 1161 if (y != z) PetscCall(VecCopy(y, z)); 1162 PetscCall(MatHYPRE_MultKernel_Private(A, 1.0, x, 1.0, z, PETSC_FALSE)); 1163 PetscFunctionReturn(PETSC_SUCCESS); 1164 } 1165 1166 static PetscErrorCode MatMultTransposeAdd_HYPRE(Mat A, Vec x, Vec y, Vec z) 1167 { 1168 PetscFunctionBegin; 1169 if (y != z) PetscCall(VecCopy(y, z)); 1170 PetscCall(MatHYPRE_MultKernel_Private(A, 1.0, x, 1.0, z, PETSC_TRUE)); 1171 PetscFunctionReturn(PETSC_SUCCESS); 1172 } 1173 1174 /* y = a * A * x + b * y or y = a * A^t * x + b * y depending on trans */ 1175 static PetscErrorCode MatHYPRE_MultKernel_Private(Mat A, HYPRE_Complex a, Vec x, HYPRE_Complex b, Vec y, PetscBool trans) 1176 { 1177 Mat_HYPRE *hA = (Mat_HYPRE *)A->data; 1178 hypre_ParCSRMatrix *parcsr; 1179 hypre_ParVector *hx, *hy; 1180 1181 PetscFunctionBegin; 1182 if (trans) { 1183 PetscCall(VecHYPRE_IJVectorPushVecRead(hA->b, x)); 1184 if (b != 0.0) PetscCall(VecHYPRE_IJVectorPushVec(hA->x, y)); 1185 else PetscCall(VecHYPRE_IJVectorPushVecWrite(hA->x, y)); 1186 PetscCallExternal(HYPRE_IJVectorGetObject, hA->b->ij, (void **)&hx); 1187 PetscCallExternal(HYPRE_IJVectorGetObject, hA->x->ij, (void **)&hy); 1188 } else { 1189 PetscCall(VecHYPRE_IJVectorPushVecRead(hA->x, x)); 1190 if (b != 0.0) PetscCall(VecHYPRE_IJVectorPushVec(hA->b, y)); 1191 else PetscCall(VecHYPRE_IJVectorPushVecWrite(hA->b, y)); 1192 PetscCallExternal(HYPRE_IJVectorGetObject, hA->x->ij, (void **)&hx); 1193 PetscCallExternal(HYPRE_IJVectorGetObject, hA->b->ij, (void **)&hy); 1194 } 1195 PetscCallExternal(HYPRE_IJMatrixGetObject, hA->ij, (void **)&parcsr); 1196 if (trans) { 1197 PetscCallExternal(hypre_ParCSRMatrixMatvecT, a, parcsr, hx, b, hy); 1198 } else { 1199 PetscCallExternal(hypre_ParCSRMatrixMatvec, a, parcsr, hx, b, hy); 1200 } 1201 PetscCall(VecHYPRE_IJVectorPopVec(hA->x)); 1202 PetscCall(VecHYPRE_IJVectorPopVec(hA->b)); 1203 PetscFunctionReturn(PETSC_SUCCESS); 1204 } 1205 1206 static PetscErrorCode MatDestroy_HYPRE(Mat A) 1207 { 1208 Mat_HYPRE *hA = (Mat_HYPRE *)A->data; 1209 1210 PetscFunctionBegin; 1211 PetscCall(VecHYPRE_IJVectorDestroy(&hA->x)); 1212 PetscCall(VecHYPRE_IJVectorDestroy(&hA->b)); 1213 if (hA->ij) { 1214 if (!hA->inner_free) hypre_IJMatrixObject(hA->ij) = NULL; 1215 PetscCallExternal(HYPRE_IJMatrixDestroy, hA->ij); 1216 } 1217 if (hA->comm) PetscCall(PetscCommRestoreComm(PetscObjectComm((PetscObject)A), &hA->comm)); 1218 1219 PetscCall(MatStashDestroy_Private(&A->stash)); 1220 PetscCall(PetscFree(hA->array)); 1221 1222 if (hA->cooMat) { 1223 PetscCall(MatDestroy(&hA->cooMat)); 1224 PetscStackCallExternalVoid("hypre_TFree", hypre_TFree(hA->diagJ, hA->memType)); 1225 PetscStackCallExternalVoid("hypre_TFree", hypre_TFree(hA->offdJ, hA->memType)); 1226 PetscStackCallExternalVoid("hypre_TFree", hypre_TFree(hA->diag, hA->memType)); 1227 } 1228 1229 PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatConvert_hypre_aij_C", NULL)); 1230 PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatConvert_hypre_is_C", NULL)); 1231 PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatProductSetFromOptions_seqaij_hypre_C", NULL)); 1232 PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatProductSetFromOptions_mpiaij_hypre_C", NULL)); 1233 PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatHYPRESetPreallocation_C", NULL)); 1234 PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatHYPREGetParCSR_C", NULL)); 1235 PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatSetPreallocationCOO_C", NULL)); 1236 PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatSetValuesCOO_C", NULL)); 1237 PetscCall(PetscFree(A->data)); 1238 PetscFunctionReturn(PETSC_SUCCESS); 1239 } 1240 1241 static PetscErrorCode MatSetUp_HYPRE(Mat A) 1242 { 1243 PetscFunctionBegin; 1244 PetscCall(MatHYPRESetPreallocation(A, PETSC_DEFAULT, NULL, PETSC_DEFAULT, NULL)); 1245 PetscFunctionReturn(PETSC_SUCCESS); 1246 } 1247 1248 //TODO FIX hypre_CSRMatrixMatvecOutOfPlace 1249 #if defined(PETSC_HAVE_HYPRE_DEVICE) 1250 static PetscErrorCode MatBindToCPU_HYPRE(Mat A, PetscBool bind) 1251 { 1252 Mat_HYPRE *hA = (Mat_HYPRE *)A->data; 1253 HYPRE_MemoryLocation hmem = bind ? HYPRE_MEMORY_HOST : HYPRE_MEMORY_DEVICE; 1254 1255 PetscFunctionBegin; 1256 A->boundtocpu = bind; 1257 if (hA->ij && hypre_IJMatrixAssembleFlag(hA->ij) && hmem != hypre_IJMatrixMemoryLocation(hA->ij)) { 1258 hypre_ParCSRMatrix *parcsr; 1259 PetscCallExternal(HYPRE_IJMatrixGetObject, hA->ij, (void **)&parcsr); 1260 PetscCallExternal(hypre_ParCSRMatrixMigrate, parcsr, hmem); 1261 } 1262 if (hA->x) PetscCall(VecHYPRE_IJBindToCPU(hA->x, bind)); 1263 if (hA->b) PetscCall(VecHYPRE_IJBindToCPU(hA->b, bind)); 1264 PetscFunctionReturn(PETSC_SUCCESS); 1265 } 1266 #endif 1267 1268 static PetscErrorCode MatAssemblyEnd_HYPRE(Mat A, MatAssemblyType mode) 1269 { 1270 Mat_HYPRE *hA = (Mat_HYPRE *)A->data; 1271 PetscMPIInt n; 1272 PetscInt i, j, rstart, ncols, flg; 1273 PetscInt *row, *col; 1274 PetscScalar *val; 1275 1276 PetscFunctionBegin; 1277 PetscCheck(mode != MAT_FLUSH_ASSEMBLY, PetscObjectComm((PetscObject)A), PETSC_ERR_SUP, "MAT_FLUSH_ASSEMBLY currently not supported with MATHYPRE"); 1278 1279 if (!A->nooffprocentries) { 1280 while (1) { 1281 PetscCall(MatStashScatterGetMesg_Private(&A->stash, &n, &row, &col, &val, &flg)); 1282 if (!flg) break; 1283 1284 for (i = 0; i < n;) { 1285 /* Now identify the consecutive vals belonging to the same row */ 1286 for (j = i, rstart = row[j]; j < n; j++) { 1287 if (row[j] != rstart) break; 1288 } 1289 if (j < n) ncols = j - i; 1290 else ncols = n - i; 1291 /* Now assemble all these values with a single function call */ 1292 PetscCall(MatSetValues_HYPRE(A, 1, row + i, ncols, col + i, val + i, A->insertmode)); 1293 1294 i = j; 1295 } 1296 } 1297 PetscCall(MatStashScatterEnd_Private(&A->stash)); 1298 } 1299 1300 PetscCallExternal(HYPRE_IJMatrixAssemble, hA->ij); 1301 /* The assembly routine destroys the aux_matrix, we recreate it here by calling HYPRE_IJMatrixInitialize */ 1302 /* 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 */ 1303 if (!hA->sorted_full) { 1304 hypre_AuxParCSRMatrix *aux_matrix; 1305 1306 /* call destroy just to make sure we do not leak anything */ 1307 aux_matrix = (hypre_AuxParCSRMatrix *)hypre_IJMatrixTranslator(hA->ij); 1308 PetscCallExternal(hypre_AuxParCSRMatrixDestroy, aux_matrix); 1309 hypre_IJMatrixTranslator(hA->ij) = NULL; 1310 1311 /* Initialize with assembled flag -> it only recreates the aux_par_matrix */ 1312 PetscCallExternal(HYPRE_IJMatrixInitialize, hA->ij); 1313 aux_matrix = (hypre_AuxParCSRMatrix *)hypre_IJMatrixTranslator(hA->ij); 1314 if (aux_matrix) { 1315 hypre_AuxParCSRMatrixNeedAux(aux_matrix) = 1; /* see comment in MatHYPRESetPreallocation_HYPRE */ 1316 #if PETSC_PKG_HYPRE_VERSION_LT(2, 19, 0) 1317 PetscCallExternal(hypre_AuxParCSRMatrixInitialize, aux_matrix); 1318 #else 1319 PetscCallExternal(hypre_AuxParCSRMatrixInitialize_v2, aux_matrix, HYPRE_MEMORY_HOST); 1320 #endif 1321 } 1322 } 1323 { 1324 hypre_ParCSRMatrix *parcsr; 1325 1326 PetscCallExternal(HYPRE_IJMatrixGetObject, hA->ij, (void **)&parcsr); 1327 if (!hypre_ParCSRMatrixCommPkg(parcsr)) PetscCallExternal(hypre_MatvecCommPkgCreate, parcsr); 1328 } 1329 if (!hA->x) PetscCall(VecHYPRE_IJVectorCreate(A->cmap, &hA->x)); 1330 if (!hA->b) PetscCall(VecHYPRE_IJVectorCreate(A->rmap, &hA->b)); 1331 #if defined(PETSC_HAVE_HYPRE_DEVICE) 1332 PetscCall(MatBindToCPU_HYPRE(A, A->boundtocpu)); 1333 #endif 1334 PetscFunctionReturn(PETSC_SUCCESS); 1335 } 1336 1337 static PetscErrorCode MatGetArray_HYPRE(Mat A, PetscInt size, void **array) 1338 { 1339 Mat_HYPRE *hA = (Mat_HYPRE *)A->data; 1340 1341 PetscFunctionBegin; 1342 PetscCheck(hA->available, PETSC_COMM_SELF, PETSC_ERR_ARG_NULL, "Temporary space is in use"); 1343 1344 if (hA->size >= size) { 1345 *array = hA->array; 1346 } else { 1347 PetscCall(PetscFree(hA->array)); 1348 hA->size = size; 1349 PetscCall(PetscMalloc(hA->size, &hA->array)); 1350 *array = hA->array; 1351 } 1352 1353 hA->available = PETSC_FALSE; 1354 PetscFunctionReturn(PETSC_SUCCESS); 1355 } 1356 1357 static PetscErrorCode MatRestoreArray_HYPRE(Mat A, void **array) 1358 { 1359 Mat_HYPRE *hA = (Mat_HYPRE *)A->data; 1360 1361 PetscFunctionBegin; 1362 *array = NULL; 1363 hA->available = PETSC_TRUE; 1364 PetscFunctionReturn(PETSC_SUCCESS); 1365 } 1366 1367 static PetscErrorCode MatSetValues_HYPRE(Mat A, PetscInt nr, const PetscInt rows[], PetscInt nc, const PetscInt cols[], const PetscScalar v[], InsertMode ins) 1368 { 1369 Mat_HYPRE *hA = (Mat_HYPRE *)A->data; 1370 PetscScalar *vals = (PetscScalar *)v; 1371 HYPRE_Complex *sscr; 1372 PetscInt *cscr[2]; 1373 PetscInt i, nzc; 1374 void *array = NULL; 1375 1376 PetscFunctionBegin; 1377 PetscCall(MatGetArray_HYPRE(A, sizeof(PetscInt) * (2 * nc) + sizeof(HYPRE_Complex) * nc * nr, &array)); 1378 cscr[0] = (PetscInt *)array; 1379 cscr[1] = ((PetscInt *)array) + nc; 1380 sscr = (HYPRE_Complex *)(((PetscInt *)array) + nc * 2); 1381 for (i = 0, nzc = 0; i < nc; i++) { 1382 if (cols[i] >= 0) { 1383 cscr[0][nzc] = cols[i]; 1384 cscr[1][nzc++] = i; 1385 } 1386 } 1387 if (!nzc) { 1388 PetscCall(MatRestoreArray_HYPRE(A, &array)); 1389 PetscFunctionReturn(PETSC_SUCCESS); 1390 } 1391 1392 #if 0 //defined(PETSC_HAVE_HYPRE_DEVICE) 1393 if (HYPRE_MEMORY_HOST != hypre_IJMatrixMemoryLocation(hA->ij)) { 1394 hypre_ParCSRMatrix *parcsr; 1395 1396 PetscCallExternal(HYPRE_IJMatrixGetObject,hA->ij,(void**)&parcsr); 1397 PetscCallExternal(hypre_ParCSRMatrixMigrate,parcsr, HYPRE_MEMORY_HOST); 1398 } 1399 #endif 1400 1401 if (ins == ADD_VALUES) { 1402 for (i = 0; i < nr; i++) { 1403 if (rows[i] >= 0) { 1404 PetscInt j; 1405 HYPRE_Int hnc = (HYPRE_Int)nzc; 1406 1407 PetscCheck((PetscInt)hnc == nzc, PETSC_COMM_SELF, PETSC_ERR_SUP, "Hypre overflow! number of columns %" PetscInt_FMT " for row %" PetscInt_FMT, nzc, rows[i]); 1408 for (j = 0; j < nzc; j++) PetscCall(PetscHYPREScalarCast(vals[cscr[1][j]], &sscr[j])); 1409 PetscCallExternal(HYPRE_IJMatrixAddToValues, hA->ij, 1, &hnc, (HYPRE_BigInt *)(rows + i), (HYPRE_BigInt *)cscr[0], sscr); 1410 } 1411 vals += nc; 1412 } 1413 } else { /* INSERT_VALUES */ 1414 PetscInt rst, ren; 1415 1416 PetscCall(MatGetOwnershipRange(A, &rst, &ren)); 1417 for (i = 0; i < nr; i++) { 1418 if (rows[i] >= 0) { 1419 PetscInt j; 1420 HYPRE_Int hnc = (HYPRE_Int)nzc; 1421 1422 PetscCheck((PetscInt)hnc == nzc, PETSC_COMM_SELF, PETSC_ERR_SUP, "Hypre overflow! number of columns %" PetscInt_FMT " for row %" PetscInt_FMT, nzc, rows[i]); 1423 for (j = 0; j < nzc; j++) PetscCall(PetscHYPREScalarCast(vals[cscr[1][j]], &sscr[j])); 1424 /* nonlocal values */ 1425 if (rows[i] < rst || rows[i] >= ren) PetscCall(MatStashValuesRow_Private(&A->stash, rows[i], nzc, cscr[0], (PetscScalar *)sscr, PETSC_FALSE)); 1426 /* local values */ 1427 else PetscCallExternal(HYPRE_IJMatrixSetValues, hA->ij, 1, &hnc, (HYPRE_BigInt *)(rows + i), (HYPRE_BigInt *)cscr[0], sscr); 1428 } 1429 vals += nc; 1430 } 1431 } 1432 1433 PetscCall(MatRestoreArray_HYPRE(A, &array)); 1434 PetscFunctionReturn(PETSC_SUCCESS); 1435 } 1436 1437 static PetscErrorCode MatHYPRESetPreallocation_HYPRE(Mat A, PetscInt dnz, const PetscInt dnnz[], PetscInt onz, const PetscInt onnz[]) 1438 { 1439 Mat_HYPRE *hA = (Mat_HYPRE *)A->data; 1440 HYPRE_Int *hdnnz, *honnz; 1441 PetscInt i, rs, re, cs, ce, bs; 1442 PetscMPIInt size; 1443 1444 PetscFunctionBegin; 1445 PetscCall(PetscLayoutSetUp(A->rmap)); 1446 PetscCall(PetscLayoutSetUp(A->cmap)); 1447 rs = A->rmap->rstart; 1448 re = A->rmap->rend; 1449 cs = A->cmap->rstart; 1450 ce = A->cmap->rend; 1451 if (!hA->ij) { 1452 PetscCallExternal(HYPRE_IJMatrixCreate, hA->comm, rs, re - 1, cs, ce - 1, &hA->ij); 1453 PetscCallExternal(HYPRE_IJMatrixSetObjectType, hA->ij, HYPRE_PARCSR); 1454 } else { 1455 HYPRE_BigInt hrs, hre, hcs, hce; 1456 PetscCallExternal(HYPRE_IJMatrixGetLocalRange, hA->ij, &hrs, &hre, &hcs, &hce); 1457 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); 1458 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); 1459 } 1460 PetscCall(MatGetBlockSize(A, &bs)); 1461 if (dnz == PETSC_DEFAULT || dnz == PETSC_DECIDE) dnz = 10 * bs; 1462 if (onz == PETSC_DEFAULT || onz == PETSC_DECIDE) onz = 10 * bs; 1463 1464 if (!dnnz) { 1465 PetscCall(PetscMalloc1(A->rmap->n, &hdnnz)); 1466 for (i = 0; i < A->rmap->n; i++) hdnnz[i] = dnz; 1467 } else { 1468 hdnnz = (HYPRE_Int *)dnnz; 1469 } 1470 PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)A), &size)); 1471 if (size > 1) { 1472 hypre_AuxParCSRMatrix *aux_matrix; 1473 if (!onnz) { 1474 PetscCall(PetscMalloc1(A->rmap->n, &honnz)); 1475 for (i = 0; i < A->rmap->n; i++) honnz[i] = onz; 1476 } else honnz = (HYPRE_Int *)onnz; 1477 /* SetDiagOffdSizes sets hypre_AuxParCSRMatrixNeedAux(aux_matrix) = 0, since it seems 1478 they assume the user will input the entire row values, properly sorted 1479 In PETSc, we don't make such an assumption and set this flag to 1, 1480 unless the option MAT_SORTED_FULL is set to true. 1481 Also, to avoid possible memory leaks, we destroy and recreate the translator 1482 This has to be done here, as HYPRE_IJMatrixInitialize will properly initialize 1483 the IJ matrix for us */ 1484 aux_matrix = (hypre_AuxParCSRMatrix *)hypre_IJMatrixTranslator(hA->ij); 1485 hypre_AuxParCSRMatrixDestroy(aux_matrix); 1486 hypre_IJMatrixTranslator(hA->ij) = NULL; 1487 PetscCallExternal(HYPRE_IJMatrixSetDiagOffdSizes, hA->ij, hdnnz, honnz); 1488 aux_matrix = (hypre_AuxParCSRMatrix *)hypre_IJMatrixTranslator(hA->ij); 1489 hypre_AuxParCSRMatrixNeedAux(aux_matrix) = !hA->sorted_full; 1490 } else { 1491 honnz = NULL; 1492 PetscCallExternal(HYPRE_IJMatrixSetRowSizes, hA->ij, hdnnz); 1493 } 1494 1495 /* reset assembled flag and call the initialize method */ 1496 hypre_IJMatrixAssembleFlag(hA->ij) = 0; 1497 #if PETSC_PKG_HYPRE_VERSION_LT(2, 19, 0) 1498 PetscCallExternal(HYPRE_IJMatrixInitialize, hA->ij); 1499 #else 1500 PetscCallExternal(HYPRE_IJMatrixInitialize_v2, hA->ij, HYPRE_MEMORY_HOST); 1501 #endif 1502 if (!dnnz) PetscCall(PetscFree(hdnnz)); 1503 if (!onnz && honnz) PetscCall(PetscFree(honnz)); 1504 /* Match AIJ logic */ 1505 A->preallocated = PETSC_TRUE; 1506 A->assembled = PETSC_FALSE; 1507 PetscFunctionReturn(PETSC_SUCCESS); 1508 } 1509 1510 /*@C 1511 MatHYPRESetPreallocation - Preallocates memory for a sparse parallel matrix in HYPRE IJ format 1512 1513 Collective 1514 1515 Input Parameters: 1516 + A - the matrix 1517 . dnz - number of nonzeros per row in DIAGONAL portion of local submatrix 1518 (same value is used for all local rows) 1519 . dnnz - array containing the number of nonzeros in the various rows of the 1520 DIAGONAL portion of the local submatrix (possibly different for each row) 1521 or NULL (`PETSC_NULL_INTEGER` in Fortran), if d_nz is used to specify the nonzero structure. 1522 The size of this array is equal to the number of local rows, i.e 'm'. 1523 For matrices that will be factored, you must leave room for (and set) 1524 the diagonal entry even if it is zero. 1525 . onz - number of nonzeros per row in the OFF-DIAGONAL portion of local 1526 submatrix (same value is used for all local rows). 1527 - onnz - array containing the number of nonzeros in the various rows of the 1528 OFF-DIAGONAL portion of the local submatrix (possibly different for 1529 each row) or NULL (`PETSC_NULL_INTEGER` in Fortran), if o_nz is used to specify the nonzero 1530 structure. The size of this array is equal to the number 1531 of local rows, i.e 'm'. 1532 1533 Note: 1534 If the *nnz parameter is given then the *nz parameter is ignored; for sequential matrices, onz and onnz are ignored. 1535 1536 Level: intermediate 1537 1538 .seealso: `MatCreate()`, `MatMPIAIJSetPreallocation()`, `MATHYPRE` 1539 @*/ 1540 PetscErrorCode MatHYPRESetPreallocation(Mat A, PetscInt dnz, const PetscInt dnnz[], PetscInt onz, const PetscInt onnz[]) 1541 { 1542 PetscFunctionBegin; 1543 PetscValidHeaderSpecific(A, MAT_CLASSID, 1); 1544 PetscValidType(A, 1); 1545 PetscTryMethod(A, "MatHYPRESetPreallocation_C", (Mat, PetscInt, const PetscInt[], PetscInt, const PetscInt[]), (A, dnz, dnnz, onz, onnz)); 1546 PetscFunctionReturn(PETSC_SUCCESS); 1547 } 1548 1549 /* 1550 MatCreateFromParCSR - Creates a matrix from a hypre_ParCSRMatrix 1551 1552 Collective 1553 1554 Input Parameters: 1555 + parcsr - the pointer to the hypre_ParCSRMatrix 1556 . mtype - matrix type to be created. Currently MATAIJ, MATIS and MATHYPRE are supported. 1557 - copymode - PETSc copying options 1558 1559 Output Parameter: 1560 . A - the matrix 1561 1562 Level: intermediate 1563 1564 .seealso: `MatHYPRE`, `PetscCopyMode` 1565 */ 1566 PETSC_EXTERN PetscErrorCode MatCreateFromParCSR(hypre_ParCSRMatrix *parcsr, MatType mtype, PetscCopyMode copymode, Mat *A) 1567 { 1568 Mat T; 1569 Mat_HYPRE *hA; 1570 MPI_Comm comm; 1571 PetscInt rstart, rend, cstart, cend, M, N; 1572 PetscBool isseqaij, isseqaijmkl, ismpiaij, isaij, ishyp, isis; 1573 1574 PetscFunctionBegin; 1575 comm = hypre_ParCSRMatrixComm(parcsr); 1576 PetscCall(PetscStrcmp(mtype, MATSEQAIJ, &isseqaij)); 1577 PetscCall(PetscStrcmp(mtype, MATSEQAIJMKL, &isseqaijmkl)); 1578 PetscCall(PetscStrcmp(mtype, MATMPIAIJ, &ismpiaij)); 1579 PetscCall(PetscStrcmp(mtype, MATAIJ, &isaij)); 1580 PetscCall(PetscStrcmp(mtype, MATHYPRE, &ishyp)); 1581 PetscCall(PetscStrcmp(mtype, MATIS, &isis)); 1582 isaij = (PetscBool)(isseqaij || isseqaijmkl || ismpiaij || isaij); 1583 /* TODO */ 1584 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); 1585 /* access ParCSRMatrix */ 1586 rstart = hypre_ParCSRMatrixFirstRowIndex(parcsr); 1587 rend = hypre_ParCSRMatrixLastRowIndex(parcsr); 1588 cstart = hypre_ParCSRMatrixFirstColDiag(parcsr); 1589 cend = hypre_ParCSRMatrixLastColDiag(parcsr); 1590 M = hypre_ParCSRMatrixGlobalNumRows(parcsr); 1591 N = hypre_ParCSRMatrixGlobalNumCols(parcsr); 1592 1593 /* fix for empty local rows/columns */ 1594 if (rend < rstart) rend = rstart; 1595 if (cend < cstart) cend = cstart; 1596 1597 /* PETSc convention */ 1598 rend++; 1599 cend++; 1600 rend = PetscMin(rend, M); 1601 cend = PetscMin(cend, N); 1602 1603 /* create PETSc matrix with MatHYPRE */ 1604 PetscCall(MatCreate(comm, &T)); 1605 PetscCall(MatSetSizes(T, rend - rstart, cend - cstart, M, N)); 1606 PetscCall(MatSetType(T, MATHYPRE)); 1607 hA = (Mat_HYPRE *)(T->data); 1608 1609 /* create HYPRE_IJMatrix */ 1610 PetscCallExternal(HYPRE_IJMatrixCreate, hA->comm, rstart, rend - 1, cstart, cend - 1, &hA->ij); 1611 PetscCallExternal(HYPRE_IJMatrixSetObjectType, hA->ij, HYPRE_PARCSR); 1612 1613 // TODO DEV 1614 /* create new ParCSR object if needed */ 1615 if (ishyp && copymode == PETSC_COPY_VALUES) { 1616 hypre_ParCSRMatrix *new_parcsr; 1617 #if PETSC_PKG_HYPRE_VERSION_LT(2, 18, 0) 1618 hypre_CSRMatrix *hdiag, *hoffd, *ndiag, *noffd; 1619 1620 new_parcsr = hypre_ParCSRMatrixClone(parcsr, 0); 1621 hdiag = hypre_ParCSRMatrixDiag(parcsr); 1622 hoffd = hypre_ParCSRMatrixOffd(parcsr); 1623 ndiag = hypre_ParCSRMatrixDiag(new_parcsr); 1624 noffd = hypre_ParCSRMatrixOffd(new_parcsr); 1625 PetscCall(PetscArraycpy(hypre_CSRMatrixData(ndiag), hypre_CSRMatrixData(hdiag), hypre_CSRMatrixNumNonzeros(hdiag))); 1626 PetscCall(PetscArraycpy(hypre_CSRMatrixData(noffd), hypre_CSRMatrixData(hoffd), hypre_CSRMatrixNumNonzeros(hoffd))); 1627 #else 1628 new_parcsr = hypre_ParCSRMatrixClone(parcsr, 1); 1629 #endif 1630 parcsr = new_parcsr; 1631 copymode = PETSC_OWN_POINTER; 1632 } 1633 1634 /* set ParCSR object */ 1635 hypre_IJMatrixObject(hA->ij) = parcsr; 1636 T->preallocated = PETSC_TRUE; 1637 1638 /* set assembled flag */ 1639 hypre_IJMatrixAssembleFlag(hA->ij) = 1; 1640 #if 0 1641 PetscCallExternal(HYPRE_IJMatrixInitialize,hA->ij); 1642 #endif 1643 if (ishyp) { 1644 PetscMPIInt myid = 0; 1645 1646 /* make sure we always have row_starts and col_starts available */ 1647 if (HYPRE_AssumedPartitionCheck()) PetscCallMPI(MPI_Comm_rank(comm, &myid)); 1648 #if defined(hypre_ParCSRMatrixOwnsRowStarts) 1649 if (!hypre_ParCSRMatrixOwnsColStarts(parcsr)) { 1650 PetscLayout map; 1651 1652 PetscCall(MatGetLayouts(T, NULL, &map)); 1653 PetscCall(PetscLayoutSetUp(map)); 1654 hypre_ParCSRMatrixColStarts(parcsr) = (HYPRE_BigInt *)(map->range + myid); 1655 } 1656 if (!hypre_ParCSRMatrixOwnsRowStarts(parcsr)) { 1657 PetscLayout map; 1658 1659 PetscCall(MatGetLayouts(T, &map, NULL)); 1660 PetscCall(PetscLayoutSetUp(map)); 1661 hypre_ParCSRMatrixRowStarts(parcsr) = (HYPRE_BigInt *)(map->range + myid); 1662 } 1663 #endif 1664 /* prevent from freeing the pointer */ 1665 if (copymode == PETSC_USE_POINTER) hA->inner_free = PETSC_FALSE; 1666 *A = T; 1667 PetscCall(MatSetOption(*A, MAT_SORTED_FULL, PETSC_TRUE)); 1668 PetscCall(MatAssemblyBegin(*A, MAT_FINAL_ASSEMBLY)); 1669 PetscCall(MatAssemblyEnd(*A, MAT_FINAL_ASSEMBLY)); 1670 } else if (isaij) { 1671 if (copymode != PETSC_OWN_POINTER) { 1672 /* prevent from freeing the pointer */ 1673 hA->inner_free = PETSC_FALSE; 1674 PetscCall(MatConvert_HYPRE_AIJ(T, MATAIJ, MAT_INITIAL_MATRIX, A)); 1675 PetscCall(MatDestroy(&T)); 1676 } else { /* AIJ return type with PETSC_OWN_POINTER */ 1677 PetscCall(MatConvert_HYPRE_AIJ(T, MATAIJ, MAT_INPLACE_MATRIX, &T)); 1678 *A = T; 1679 } 1680 } else if (isis) { 1681 PetscCall(MatConvert_HYPRE_IS(T, MATIS, MAT_INITIAL_MATRIX, A)); 1682 if (copymode != PETSC_OWN_POINTER) hA->inner_free = PETSC_FALSE; 1683 PetscCall(MatDestroy(&T)); 1684 } 1685 PetscFunctionReturn(PETSC_SUCCESS); 1686 } 1687 1688 static PetscErrorCode MatHYPREGetParCSR_HYPRE(Mat A, hypre_ParCSRMatrix **parcsr) 1689 { 1690 Mat_HYPRE *hA = (Mat_HYPRE *)A->data; 1691 HYPRE_Int type; 1692 1693 PetscFunctionBegin; 1694 PetscCheck(hA->ij, PetscObjectComm((PetscObject)A), PETSC_ERR_PLIB, "HYPRE_IJMatrix not present"); 1695 PetscCallExternal(HYPRE_IJMatrixGetObjectType, hA->ij, &type); 1696 PetscCheck(type == HYPRE_PARCSR, PetscObjectComm((PetscObject)A), PETSC_ERR_SUP, "HYPRE_IJMatrix is not of type HYPRE_PARCSR"); 1697 PetscCallExternal(HYPRE_IJMatrixGetObject, hA->ij, (void **)parcsr); 1698 PetscFunctionReturn(PETSC_SUCCESS); 1699 } 1700 1701 /* 1702 MatHYPREGetParCSR - Gets the pointer to the ParCSR matrix 1703 1704 Not collective 1705 1706 Input Parameters: 1707 + A - the MATHYPRE object 1708 1709 Output Parameter: 1710 . parcsr - the pointer to the hypre_ParCSRMatrix 1711 1712 Level: intermediate 1713 1714 .seealso: `MatHYPRE`, `PetscCopyMode` 1715 */ 1716 PetscErrorCode MatHYPREGetParCSR(Mat A, hypre_ParCSRMatrix **parcsr) 1717 { 1718 PetscFunctionBegin; 1719 PetscValidHeaderSpecific(A, MAT_CLASSID, 1); 1720 PetscValidType(A, 1); 1721 PetscUseMethod(A, "MatHYPREGetParCSR_C", (Mat, hypre_ParCSRMatrix **), (A, parcsr)); 1722 PetscFunctionReturn(PETSC_SUCCESS); 1723 } 1724 1725 static PetscErrorCode MatMissingDiagonal_HYPRE(Mat A, PetscBool *missing, PetscInt *dd) 1726 { 1727 hypre_ParCSRMatrix *parcsr; 1728 hypre_CSRMatrix *ha; 1729 PetscInt rst; 1730 1731 PetscFunctionBegin; 1732 PetscCheck(A->rmap->n == A->cmap->n, PETSC_COMM_SELF, PETSC_ERR_SUP, "Not implemented with non-square diagonal blocks"); 1733 PetscCall(MatGetOwnershipRange(A, &rst, NULL)); 1734 PetscCall(MatHYPREGetParCSR_HYPRE(A, &parcsr)); 1735 if (missing) *missing = PETSC_FALSE; 1736 if (dd) *dd = -1; 1737 ha = hypre_ParCSRMatrixDiag(parcsr); 1738 if (ha) { 1739 PetscInt size, i; 1740 HYPRE_Int *ii, *jj; 1741 1742 size = hypre_CSRMatrixNumRows(ha); 1743 ii = hypre_CSRMatrixI(ha); 1744 jj = hypre_CSRMatrixJ(ha); 1745 for (i = 0; i < size; i++) { 1746 PetscInt j; 1747 PetscBool found = PETSC_FALSE; 1748 1749 for (j = ii[i]; j < ii[i + 1] && !found; j++) found = (jj[j] == i) ? PETSC_TRUE : PETSC_FALSE; 1750 1751 if (!found) { 1752 PetscCall(PetscInfo(A, "Matrix is missing local diagonal entry %" PetscInt_FMT "\n", i)); 1753 if (missing) *missing = PETSC_TRUE; 1754 if (dd) *dd = i + rst; 1755 PetscFunctionReturn(PETSC_SUCCESS); 1756 } 1757 } 1758 if (!size) { 1759 PetscCall(PetscInfo(A, "Matrix has no diagonal entries therefore is missing diagonal\n")); 1760 if (missing) *missing = PETSC_TRUE; 1761 if (dd) *dd = rst; 1762 } 1763 } else { 1764 PetscCall(PetscInfo(A, "Matrix has no diagonal entries therefore is missing diagonal\n")); 1765 if (missing) *missing = PETSC_TRUE; 1766 if (dd) *dd = rst; 1767 } 1768 PetscFunctionReturn(PETSC_SUCCESS); 1769 } 1770 1771 static PetscErrorCode MatScale_HYPRE(Mat A, PetscScalar s) 1772 { 1773 hypre_ParCSRMatrix *parcsr; 1774 #if PETSC_PKG_HYPRE_VERSION_LT(2, 19, 0) 1775 hypre_CSRMatrix *ha; 1776 #endif 1777 HYPRE_Complex hs; 1778 1779 PetscFunctionBegin; 1780 PetscCall(PetscHYPREScalarCast(s, &hs)); 1781 PetscCall(MatHYPREGetParCSR_HYPRE(A, &parcsr)); 1782 #if PETSC_PKG_HYPRE_VERSION_GE(2, 19, 0) 1783 PetscCallExternal(hypre_ParCSRMatrixScale, parcsr, hs); 1784 #else /* diagonal part */ 1785 ha = hypre_ParCSRMatrixDiag(parcsr); 1786 if (ha) { 1787 PetscInt size, i; 1788 HYPRE_Int *ii; 1789 HYPRE_Complex *a; 1790 1791 size = hypre_CSRMatrixNumRows(ha); 1792 a = hypre_CSRMatrixData(ha); 1793 ii = hypre_CSRMatrixI(ha); 1794 for (i = 0; i < ii[size]; i++) a[i] *= hs; 1795 } 1796 /* offdiagonal part */ 1797 ha = hypre_ParCSRMatrixOffd(parcsr); 1798 if (ha) { 1799 PetscInt size, i; 1800 HYPRE_Int *ii; 1801 HYPRE_Complex *a; 1802 1803 size = hypre_CSRMatrixNumRows(ha); 1804 a = hypre_CSRMatrixData(ha); 1805 ii = hypre_CSRMatrixI(ha); 1806 for (i = 0; i < ii[size]; i++) a[i] *= hs; 1807 } 1808 #endif 1809 PetscFunctionReturn(PETSC_SUCCESS); 1810 } 1811 1812 static PetscErrorCode MatZeroRowsColumns_HYPRE(Mat A, PetscInt numRows, const PetscInt rows[], PetscScalar diag, Vec x, Vec b) 1813 { 1814 hypre_ParCSRMatrix *parcsr; 1815 HYPRE_Int *lrows; 1816 PetscInt rst, ren, i; 1817 1818 PetscFunctionBegin; 1819 PetscCheck(!x && !b, PetscObjectComm((PetscObject)A), PETSC_ERR_SUP, "To be implemented"); 1820 PetscCall(MatHYPREGetParCSR_HYPRE(A, &parcsr)); 1821 PetscCall(PetscMalloc1(numRows, &lrows)); 1822 PetscCall(MatGetOwnershipRange(A, &rst, &ren)); 1823 for (i = 0; i < numRows; i++) { 1824 PetscCheck(rows[i] >= rst && rows[i] < ren, PETSC_COMM_SELF, PETSC_ERR_SUP, "Non-local rows not yet supported"); 1825 lrows[i] = rows[i] - rst; 1826 } 1827 PetscCallExternal(hypre_ParCSRMatrixEliminateRowsCols, parcsr, numRows, lrows); 1828 PetscCall(PetscFree(lrows)); 1829 PetscFunctionReturn(PETSC_SUCCESS); 1830 } 1831 1832 static PetscErrorCode MatZeroEntries_HYPRE_CSRMatrix(hypre_CSRMatrix *ha) 1833 { 1834 PetscFunctionBegin; 1835 if (ha) { 1836 HYPRE_Int *ii, size; 1837 HYPRE_Complex *a; 1838 1839 size = hypre_CSRMatrixNumRows(ha); 1840 a = hypre_CSRMatrixData(ha); 1841 ii = hypre_CSRMatrixI(ha); 1842 1843 if (a) PetscCall(PetscArrayzero(a, ii[size])); 1844 } 1845 PetscFunctionReturn(PETSC_SUCCESS); 1846 } 1847 1848 PetscErrorCode MatZeroEntries_HYPRE(Mat A) 1849 { 1850 Mat_HYPRE *hA = (Mat_HYPRE *)A->data; 1851 1852 PetscFunctionBegin; 1853 if (HYPRE_MEMORY_DEVICE == hypre_IJMatrixMemoryLocation(hA->ij)) { 1854 PetscCallExternal(HYPRE_IJMatrixSetConstantValues, hA->ij, 0.0); 1855 } else { 1856 hypre_ParCSRMatrix *parcsr; 1857 1858 PetscCall(MatHYPREGetParCSR_HYPRE(A, &parcsr)); 1859 PetscCall(MatZeroEntries_HYPRE_CSRMatrix(hypre_ParCSRMatrixDiag(parcsr))); 1860 PetscCall(MatZeroEntries_HYPRE_CSRMatrix(hypre_ParCSRMatrixOffd(parcsr))); 1861 } 1862 PetscFunctionReturn(PETSC_SUCCESS); 1863 } 1864 1865 static PetscErrorCode MatZeroRows_HYPRE_CSRMatrix(hypre_CSRMatrix *hA, PetscInt N, const PetscInt rows[], HYPRE_Complex diag) 1866 { 1867 PetscInt ii; 1868 HYPRE_Int *i, *j; 1869 HYPRE_Complex *a; 1870 1871 PetscFunctionBegin; 1872 if (!hA) PetscFunctionReturn(PETSC_SUCCESS); 1873 1874 i = hypre_CSRMatrixI(hA); 1875 j = hypre_CSRMatrixJ(hA); 1876 a = hypre_CSRMatrixData(hA); 1877 1878 for (ii = 0; ii < N; ii++) { 1879 HYPRE_Int jj, ibeg, iend, irow; 1880 1881 irow = rows[ii]; 1882 ibeg = i[irow]; 1883 iend = i[irow + 1]; 1884 for (jj = ibeg; jj < iend; jj++) 1885 if (j[jj] == irow) a[jj] = diag; 1886 else a[jj] = 0.0; 1887 } 1888 PetscFunctionReturn(PETSC_SUCCESS); 1889 } 1890 1891 static PetscErrorCode MatZeroRows_HYPRE(Mat A, PetscInt N, const PetscInt rows[], PetscScalar diag, Vec x, Vec b) 1892 { 1893 hypre_ParCSRMatrix *parcsr; 1894 PetscInt *lrows, len; 1895 HYPRE_Complex hdiag; 1896 1897 PetscFunctionBegin; 1898 PetscCheck(!x && !b, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Does not support to modify the solution and the right hand size"); 1899 PetscCall(PetscHYPREScalarCast(diag, &hdiag)); 1900 /* retrieve the internal matrix */ 1901 PetscCall(MatHYPREGetParCSR_HYPRE(A, &parcsr)); 1902 /* get locally owned rows */ 1903 PetscCall(MatZeroRowsMapLocal_Private(A, N, rows, &len, &lrows)); 1904 /* zero diagonal part */ 1905 PetscCall(MatZeroRows_HYPRE_CSRMatrix(hypre_ParCSRMatrixDiag(parcsr), len, lrows, hdiag)); 1906 /* zero off-diagonal part */ 1907 PetscCall(MatZeroRows_HYPRE_CSRMatrix(hypre_ParCSRMatrixOffd(parcsr), len, lrows, 0.0)); 1908 1909 PetscCall(PetscFree(lrows)); 1910 PetscFunctionReturn(PETSC_SUCCESS); 1911 } 1912 1913 static PetscErrorCode MatAssemblyBegin_HYPRE(Mat mat, MatAssemblyType mode) 1914 { 1915 PetscFunctionBegin; 1916 if (mat->nooffprocentries) PetscFunctionReturn(PETSC_SUCCESS); 1917 1918 PetscCall(MatStashScatterBegin_Private(mat, &mat->stash, mat->rmap->range)); 1919 PetscFunctionReturn(PETSC_SUCCESS); 1920 } 1921 1922 static PetscErrorCode MatGetRow_HYPRE(Mat A, PetscInt row, PetscInt *nz, PetscInt **idx, PetscScalar **v) 1923 { 1924 hypre_ParCSRMatrix *parcsr; 1925 HYPRE_Int hnz; 1926 1927 PetscFunctionBegin; 1928 /* retrieve the internal matrix */ 1929 PetscCall(MatHYPREGetParCSR_HYPRE(A, &parcsr)); 1930 /* call HYPRE API */ 1931 PetscCallExternal(HYPRE_ParCSRMatrixGetRow, parcsr, row, &hnz, (HYPRE_BigInt **)idx, (HYPRE_Complex **)v); 1932 if (nz) *nz = (PetscInt)hnz; 1933 PetscFunctionReturn(PETSC_SUCCESS); 1934 } 1935 1936 static PetscErrorCode MatRestoreRow_HYPRE(Mat A, PetscInt row, PetscInt *nz, PetscInt **idx, PetscScalar **v) 1937 { 1938 hypre_ParCSRMatrix *parcsr; 1939 HYPRE_Int hnz; 1940 1941 PetscFunctionBegin; 1942 /* retrieve the internal matrix */ 1943 PetscCall(MatHYPREGetParCSR_HYPRE(A, &parcsr)); 1944 /* call HYPRE API */ 1945 hnz = nz ? (HYPRE_Int)(*nz) : 0; 1946 PetscCallExternal(HYPRE_ParCSRMatrixRestoreRow, parcsr, row, &hnz, (HYPRE_BigInt **)idx, (HYPRE_Complex **)v); 1947 PetscFunctionReturn(PETSC_SUCCESS); 1948 } 1949 1950 static PetscErrorCode MatGetValues_HYPRE(Mat A, PetscInt m, const PetscInt idxm[], PetscInt n, const PetscInt idxn[], PetscScalar v[]) 1951 { 1952 Mat_HYPRE *hA = (Mat_HYPRE *)A->data; 1953 PetscInt i; 1954 1955 PetscFunctionBegin; 1956 if (!m || !n) PetscFunctionReturn(PETSC_SUCCESS); 1957 /* Ignore negative row indices 1958 * And negative column indices should be automatically ignored in hypre 1959 * */ 1960 for (i = 0; i < m; i++) { 1961 if (idxm[i] >= 0) { 1962 HYPRE_Int hn = (HYPRE_Int)n; 1963 PetscCallExternal(HYPRE_IJMatrixGetValues, hA->ij, 1, &hn, (HYPRE_BigInt *)&idxm[i], (HYPRE_BigInt *)idxn, (HYPRE_Complex *)(v + i * n)); 1964 } 1965 } 1966 PetscFunctionReturn(PETSC_SUCCESS); 1967 } 1968 1969 static PetscErrorCode MatSetOption_HYPRE(Mat A, MatOption op, PetscBool flg) 1970 { 1971 Mat_HYPRE *hA = (Mat_HYPRE *)A->data; 1972 1973 PetscFunctionBegin; 1974 switch (op) { 1975 case MAT_NO_OFF_PROC_ENTRIES: 1976 if (flg) PetscCallExternal(HYPRE_IJMatrixSetMaxOffProcElmts, hA->ij, 0); 1977 break; 1978 case MAT_SORTED_FULL: 1979 hA->sorted_full = flg; 1980 break; 1981 default: 1982 break; 1983 } 1984 PetscFunctionReturn(PETSC_SUCCESS); 1985 } 1986 1987 static PetscErrorCode MatView_HYPRE(Mat A, PetscViewer view) 1988 { 1989 PetscViewerFormat format; 1990 1991 PetscFunctionBegin; 1992 PetscCall(PetscViewerGetFormat(view, &format)); 1993 if (format == PETSC_VIEWER_ASCII_FACTOR_INFO || format == PETSC_VIEWER_ASCII_INFO || format == PETSC_VIEWER_ASCII_INFO_DETAIL) PetscFunctionReturn(PETSC_SUCCESS); 1994 if (format != PETSC_VIEWER_NATIVE) { 1995 Mat B; 1996 hypre_ParCSRMatrix *parcsr; 1997 PetscErrorCode (*mview)(Mat, PetscViewer) = NULL; 1998 1999 PetscCall(MatHYPREGetParCSR_HYPRE(A, &parcsr)); 2000 PetscCall(MatCreateFromParCSR(parcsr, MATAIJ, PETSC_USE_POINTER, &B)); 2001 PetscCall(MatGetOperation(B, MATOP_VIEW, (void (**)(void)) & mview)); 2002 PetscCheck(mview, PetscObjectComm((PetscObject)A), PETSC_ERR_PLIB, "Missing view operation"); 2003 PetscCall((*mview)(B, view)); 2004 PetscCall(MatDestroy(&B)); 2005 } else { 2006 Mat_HYPRE *hA = (Mat_HYPRE *)A->data; 2007 PetscMPIInt size; 2008 PetscBool isascii; 2009 const char *filename; 2010 2011 /* HYPRE uses only text files */ 2012 PetscCall(PetscObjectTypeCompare((PetscObject)view, PETSCVIEWERASCII, &isascii)); 2013 PetscCheck(isascii, PetscObjectComm((PetscObject)view), PETSC_ERR_SUP, "PetscViewerType %s: native HYPRE format needs PETSCVIEWERASCII", ((PetscObject)view)->type_name); 2014 PetscCall(PetscViewerFileGetName(view, &filename)); 2015 PetscCallExternal(HYPRE_IJMatrixPrint, hA->ij, filename); 2016 PetscCallMPI(MPI_Comm_size(hA->comm, &size)); 2017 if (size > 1) { 2018 PetscCall(PetscViewerASCIIPrintf(view, "Matrix files: %s.%05d ... %s.%05d\n", filename, 0, filename, size - 1)); 2019 } else { 2020 PetscCall(PetscViewerASCIIPrintf(view, "Matrix file: %s.%05d\n", filename, 0)); 2021 } 2022 } 2023 PetscFunctionReturn(PETSC_SUCCESS); 2024 } 2025 2026 static PetscErrorCode MatDuplicate_HYPRE(Mat A, MatDuplicateOption op, Mat *B) 2027 { 2028 hypre_ParCSRMatrix *parcsr = NULL; 2029 PetscCopyMode cpmode; 2030 2031 PetscFunctionBegin; 2032 PetscCall(MatHYPREGetParCSR_HYPRE(A, &parcsr)); 2033 if (op == MAT_DO_NOT_COPY_VALUES || op == MAT_SHARE_NONZERO_PATTERN) { 2034 parcsr = hypre_ParCSRMatrixClone(parcsr, 0); 2035 cpmode = PETSC_OWN_POINTER; 2036 } else { 2037 cpmode = PETSC_COPY_VALUES; 2038 } 2039 PetscCall(MatCreateFromParCSR(parcsr, MATHYPRE, cpmode, B)); 2040 PetscFunctionReturn(PETSC_SUCCESS); 2041 } 2042 2043 static PetscErrorCode MatCopy_HYPRE(Mat A, Mat B, MatStructure str) 2044 { 2045 hypre_ParCSRMatrix *acsr, *bcsr; 2046 2047 PetscFunctionBegin; 2048 if (str == SAME_NONZERO_PATTERN && A->ops->copy == B->ops->copy) { 2049 PetscCall(MatHYPREGetParCSR_HYPRE(A, &acsr)); 2050 PetscCall(MatHYPREGetParCSR_HYPRE(B, &bcsr)); 2051 PetscCallExternal(hypre_ParCSRMatrixCopy, acsr, bcsr, 1); 2052 PetscCall(MatSetOption(B, MAT_SORTED_FULL, PETSC_TRUE)); /* "perfect" preallocation, so no need for hypre_AuxParCSRMatrixNeedAux */ 2053 PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY)); 2054 PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY)); 2055 } else { 2056 PetscCall(MatCopy_Basic(A, B, str)); 2057 } 2058 PetscFunctionReturn(PETSC_SUCCESS); 2059 } 2060 2061 static PetscErrorCode MatGetDiagonal_HYPRE(Mat A, Vec d) 2062 { 2063 hypre_ParCSRMatrix *parcsr; 2064 hypre_CSRMatrix *dmat; 2065 HYPRE_Complex *a; 2066 HYPRE_Complex *data = NULL; 2067 HYPRE_Int *diag = NULL; 2068 PetscInt i; 2069 PetscBool cong; 2070 2071 PetscFunctionBegin; 2072 PetscCall(MatHasCongruentLayouts(A, &cong)); 2073 PetscCheck(cong, PetscObjectComm((PetscObject)A), PETSC_ERR_SUP, "Only for square matrices with same local distributions of rows and columns"); 2074 if (PetscDefined(USE_DEBUG)) { 2075 PetscBool miss; 2076 PetscCall(MatMissingDiagonal(A, &miss, NULL)); 2077 PetscCheck(!miss || !A->rmap->n, PETSC_COMM_SELF, PETSC_ERR_SUP, "Not implemented when diagonal entries are missing"); 2078 } 2079 PetscCall(MatHYPREGetParCSR_HYPRE(A, &parcsr)); 2080 dmat = hypre_ParCSRMatrixDiag(parcsr); 2081 if (dmat) { 2082 /* this cast fixes the clang error: implicit conversion from 'HYPRE_Complex' (aka '_Complex double') to 'double' is not permitted in C++ */ 2083 PetscCall(VecGetArray(d, (PetscScalar **)&a)); 2084 diag = hypre_CSRMatrixI(dmat); 2085 data = hypre_CSRMatrixData(dmat); 2086 for (i = 0; i < A->rmap->n; i++) a[i] = data[diag[i]]; 2087 PetscCall(VecRestoreArray(d, (PetscScalar **)&a)); 2088 } 2089 PetscFunctionReturn(PETSC_SUCCESS); 2090 } 2091 2092 #include <petscblaslapack.h> 2093 2094 static PetscErrorCode MatAXPY_HYPRE(Mat Y, PetscScalar a, Mat X, MatStructure str) 2095 { 2096 PetscFunctionBegin; 2097 #if defined(PETSC_HAVE_HYPRE_DEVICE) 2098 { 2099 Mat B; 2100 hypre_ParCSRMatrix *x, *y, *z; 2101 2102 PetscCall(MatHYPREGetParCSR(Y, &y)); 2103 PetscCall(MatHYPREGetParCSR(X, &x)); 2104 PetscCallExternal(hypre_ParCSRMatrixAdd, 1.0, y, 1.0, x, &z); 2105 PetscCall(MatCreateFromParCSR(z, MATHYPRE, PETSC_OWN_POINTER, &B)); 2106 PetscCall(MatHeaderMerge(Y, &B)); 2107 } 2108 #else 2109 if (str == SAME_NONZERO_PATTERN) { 2110 hypre_ParCSRMatrix *x, *y; 2111 hypre_CSRMatrix *xloc, *yloc; 2112 PetscInt xnnz, ynnz; 2113 HYPRE_Complex *xarr, *yarr; 2114 PetscBLASInt one = 1, bnz; 2115 2116 PetscCall(MatHYPREGetParCSR(Y, &y)); 2117 PetscCall(MatHYPREGetParCSR(X, &x)); 2118 2119 /* diagonal block */ 2120 xloc = hypre_ParCSRMatrixDiag(x); 2121 yloc = hypre_ParCSRMatrixDiag(y); 2122 xnnz = 0; 2123 ynnz = 0; 2124 xarr = NULL; 2125 yarr = NULL; 2126 if (xloc) { 2127 xarr = hypre_CSRMatrixData(xloc); 2128 xnnz = hypre_CSRMatrixNumNonzeros(xloc); 2129 } 2130 if (yloc) { 2131 yarr = hypre_CSRMatrixData(yloc); 2132 ynnz = hypre_CSRMatrixNumNonzeros(yloc); 2133 } 2134 PetscCheck(xnnz == ynnz, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Different number of nonzeros in diagonal block %" PetscInt_FMT " != %" PetscInt_FMT, xnnz, ynnz); 2135 PetscCall(PetscBLASIntCast(xnnz, &bnz)); 2136 PetscCallBLAS("BLASaxpy", BLASaxpy_(&bnz, &a, (PetscScalar *)xarr, &one, (PetscScalar *)yarr, &one)); 2137 2138 /* off-diagonal block */ 2139 xloc = hypre_ParCSRMatrixOffd(x); 2140 yloc = hypre_ParCSRMatrixOffd(y); 2141 xnnz = 0; 2142 ynnz = 0; 2143 xarr = NULL; 2144 yarr = NULL; 2145 if (xloc) { 2146 xarr = hypre_CSRMatrixData(xloc); 2147 xnnz = hypre_CSRMatrixNumNonzeros(xloc); 2148 } 2149 if (yloc) { 2150 yarr = hypre_CSRMatrixData(yloc); 2151 ynnz = hypre_CSRMatrixNumNonzeros(yloc); 2152 } 2153 PetscCheck(xnnz == ynnz, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Different number of nonzeros in off-diagonal block %" PetscInt_FMT " != %" PetscInt_FMT, xnnz, ynnz); 2154 PetscCall(PetscBLASIntCast(xnnz, &bnz)); 2155 PetscCallBLAS("BLASaxpy", BLASaxpy_(&bnz, &a, (PetscScalar *)xarr, &one, (PetscScalar *)yarr, &one)); 2156 } else if (str == SUBSET_NONZERO_PATTERN) { 2157 PetscCall(MatAXPY_Basic(Y, a, X, str)); 2158 } else { 2159 Mat B; 2160 2161 PetscCall(MatAXPY_Basic_Preallocate(Y, X, &B)); 2162 PetscCall(MatAXPY_BasicWithPreallocation(B, Y, a, X, str)); 2163 PetscCall(MatHeaderReplace(Y, &B)); 2164 } 2165 #endif 2166 PetscFunctionReturn(PETSC_SUCCESS); 2167 } 2168 2169 static PetscErrorCode MatSetPreallocationCOO_HYPRE(Mat mat, PetscCount coo_n, PetscInt coo_i[], PetscInt coo_j[]) 2170 { 2171 MPI_Comm comm; 2172 PetscMPIInt size; 2173 PetscLayout rmap, cmap; 2174 Mat_HYPRE *hmat; 2175 hypre_ParCSRMatrix *parCSR; 2176 hypre_CSRMatrix *diag, *offd; 2177 Mat A, B, cooMat; 2178 PetscScalar *Aa, *Ba; 2179 HYPRE_MemoryLocation hypreMemtype = HYPRE_MEMORY_HOST; 2180 PetscMemType petscMemtype; 2181 MatType matType = MATAIJ; /* default type of cooMat */ 2182 2183 PetscFunctionBegin; 2184 /* Build an agent matrix cooMat whose type is either MATAIJ or MATAIJKOKKOS. 2185 It has the same sparsity pattern as mat, and also shares the data array with mat. We use cooMat to do the COO work. 2186 */ 2187 PetscCall(PetscObjectGetComm((PetscObject)mat, &comm)); 2188 PetscCallMPI(MPI_Comm_size(comm, &size)); 2189 PetscCall(PetscLayoutSetUp(mat->rmap)); 2190 PetscCall(PetscLayoutSetUp(mat->cmap)); 2191 PetscCall(MatGetLayouts(mat, &rmap, &cmap)); 2192 2193 /* I do not know how hypre_ParCSRMatrix stores diagonal elements for non-square matrices, so I just give up now */ 2194 PetscCheck(rmap->N == cmap->N, comm, PETSC_ERR_SUP, "MATHYPRE COO cannot handle non-square matrices"); 2195 2196 #if defined(PETSC_HAVE_DEVICE) 2197 if (!mat->boundtocpu) { /* mat will be on device, so will cooMat */ 2198 #if defined(PETSC_HAVE_KOKKOS) 2199 matType = MATAIJKOKKOS; 2200 #else 2201 SETERRQ(comm, PETSC_ERR_SUP, "To support MATHYPRE COO assembly on device, we need Kokkos, e.g., --download-kokkos --download-kokkos-kernels"); 2202 #endif 2203 } 2204 #endif 2205 2206 /* Do COO preallocation through cooMat */ 2207 hmat = (Mat_HYPRE *)mat->data; 2208 PetscCall(MatDestroy(&hmat->cooMat)); 2209 PetscCall(MatCreate(comm, &cooMat)); 2210 PetscCall(MatSetType(cooMat, matType)); 2211 PetscCall(MatSetLayouts(cooMat, rmap, cmap)); 2212 PetscCall(MatSetPreallocationCOO(cooMat, coo_n, coo_i, coo_j)); 2213 2214 /* Copy the sparsity pattern from cooMat to hypre IJMatrix hmat->ij */ 2215 PetscCall(MatSetOption(mat, MAT_SORTED_FULL, PETSC_TRUE)); 2216 PetscCall(MatSetOption(mat, MAT_NO_OFF_PROC_ENTRIES, PETSC_TRUE)); 2217 PetscCall(MatHYPRE_CreateFromMat(cooMat, hmat)); /* Create hmat->ij and preallocate it */ 2218 PetscCall(MatHYPRE_IJMatrixCopy(cooMat, hmat->ij)); /* Copy A's (a,i,j) to hmat->ij. To reuse code. Copying 'a' is not really needed */ 2219 2220 mat->preallocated = PETSC_TRUE; 2221 PetscCall(MatAssemblyBegin(mat, MAT_FINAL_ASSEMBLY)); 2222 PetscCall(MatAssemblyEnd(mat, MAT_FINAL_ASSEMBLY)); /* Migrate mat to device if it is bound to. Hypre builds its own SpMV context here */ 2223 2224 /* Alias cooMat's data array to IJMatrix's */ 2225 PetscCallExternal(HYPRE_IJMatrixGetObject, hmat->ij, (void **)&parCSR); 2226 diag = hypre_ParCSRMatrixDiag(parCSR); 2227 offd = hypre_ParCSRMatrixOffd(parCSR); 2228 2229 hypreMemtype = hypre_CSRMatrixMemoryLocation(diag); 2230 A = (size == 1) ? cooMat : ((Mat_MPIAIJ *)cooMat->data)->A; 2231 PetscCall(MatSeqAIJGetCSRAndMemType(A, NULL, NULL, &Aa, &petscMemtype)); 2232 PetscAssert((PetscMemTypeHost(petscMemtype) && hypreMemtype == HYPRE_MEMORY_HOST) || (PetscMemTypeDevice(petscMemtype) && hypreMemtype == HYPRE_MEMORY_DEVICE), comm, PETSC_ERR_PLIB, "PETSc and hypre's memory types mismatch"); 2233 2234 hmat->diagJ = hypre_CSRMatrixJ(diag); 2235 PetscStackCallExternalVoid("hypre_TFree", hypre_TFree(hypre_CSRMatrixData(diag), hypreMemtype)); 2236 hypre_CSRMatrixData(diag) = (HYPRE_Complex *)Aa; 2237 hypre_CSRMatrixOwnsData(diag) = 0; /* Take ownership of (j,a) away from hypre. As a result, we need to free them on our own */ 2238 2239 /* Copy diagonal pointers of A to device to facilitate MatSeqAIJMoveDiagonalValuesFront_SeqAIJKokkos */ 2240 if (hypreMemtype == HYPRE_MEMORY_DEVICE) { 2241 PetscStackCallExternalVoid("hypre_TAlloc", hmat->diag = hypre_TAlloc(PetscInt, rmap->n, hypreMemtype)); 2242 PetscCall(MatMarkDiagonal_SeqAIJ(A)); /* We need updated diagonal positions */ 2243 PetscStackCallExternalVoid("hypre_TMemcpy", hypre_TMemcpy(hmat->diag, ((Mat_SeqAIJ *)A->data)->diag, PetscInt, rmap->n, hypreMemtype, HYPRE_MEMORY_HOST)); 2244 } 2245 2246 if (size > 1) { 2247 B = ((Mat_MPIAIJ *)cooMat->data)->B; 2248 PetscCall(MatSeqAIJGetCSRAndMemType(B, NULL, NULL, &Ba, &petscMemtype)); 2249 hmat->offdJ = hypre_CSRMatrixJ(offd); 2250 PetscStackCallExternalVoid("hypre_TFree", hypre_TFree(hypre_CSRMatrixData(offd), hypreMemtype)); 2251 hypre_CSRMatrixData(offd) = (HYPRE_Complex *)Ba; 2252 hypre_CSRMatrixOwnsData(offd) = 0; 2253 } 2254 2255 /* Record cooMat for use in MatSetValuesCOO_HYPRE */ 2256 hmat->cooMat = cooMat; 2257 hmat->memType = hypreMemtype; 2258 PetscFunctionReturn(PETSC_SUCCESS); 2259 } 2260 2261 static PetscErrorCode MatSetValuesCOO_HYPRE(Mat mat, const PetscScalar v[], InsertMode imode) 2262 { 2263 Mat_HYPRE *hmat = (Mat_HYPRE *)mat->data; 2264 PetscMPIInt size; 2265 Mat A; 2266 2267 PetscFunctionBegin; 2268 PetscCheck(hmat->cooMat, hmat->comm, PETSC_ERR_PLIB, "HYPRE COO delegate matrix has not been created yet"); 2269 PetscCallMPI(MPI_Comm_size(hmat->comm, &size)); 2270 PetscCall(MatSetValuesCOO(hmat->cooMat, v, imode)); 2271 2272 /* Move diagonal elements of the diagonal block to the front of their row, as needed by ParCSRMatrix. So damn hacky */ 2273 A = (size == 1) ? hmat->cooMat : ((Mat_MPIAIJ *)hmat->cooMat->data)->A; 2274 if (hmat->memType == HYPRE_MEMORY_HOST) { 2275 Mat_SeqAIJ *aij = (Mat_SeqAIJ *)A->data; 2276 PetscInt i, m, *Ai = aij->i, *Adiag = aij->diag; 2277 PetscScalar *Aa = aij->a, tmp; 2278 2279 PetscCall(MatGetSize(A, &m, NULL)); 2280 for (i = 0; i < m; i++) { 2281 if (Adiag[i] >= Ai[i] && Adiag[i] < Ai[i + 1]) { /* Digonal element of this row exists in a[] and j[] */ 2282 tmp = Aa[Ai[i]]; 2283 Aa[Ai[i]] = Aa[Adiag[i]]; 2284 Aa[Adiag[i]] = tmp; 2285 } 2286 } 2287 } else { 2288 #if defined(PETSC_HAVE_KOKKOS_KERNELS) 2289 PetscCall(MatSeqAIJMoveDiagonalValuesFront_SeqAIJKokkos(A, hmat->diag)); 2290 #endif 2291 } 2292 PetscFunctionReturn(PETSC_SUCCESS); 2293 } 2294 2295 /*MC 2296 MATHYPRE - MATHYPRE = "hypre" - A matrix type to be used for sequential and parallel sparse matrices 2297 based on the hypre IJ interface. 2298 2299 Level: intermediate 2300 2301 .seealso: `MatCreate()`, `MatHYPRESetPreallocation` 2302 M*/ 2303 2304 PETSC_EXTERN PetscErrorCode MatCreate_HYPRE(Mat B) 2305 { 2306 Mat_HYPRE *hB; 2307 2308 PetscFunctionBegin; 2309 PetscCall(PetscNew(&hB)); 2310 2311 hB->inner_free = PETSC_TRUE; 2312 hB->available = PETSC_TRUE; 2313 hB->sorted_full = PETSC_FALSE; /* no assumption whether column indices are sorted or not */ 2314 hB->size = 0; 2315 hB->array = NULL; 2316 2317 B->data = (void *)hB; 2318 B->assembled = PETSC_FALSE; 2319 2320 PetscCall(PetscMemzero(B->ops, sizeof(struct _MatOps))); 2321 B->ops->mult = MatMult_HYPRE; 2322 B->ops->multtranspose = MatMultTranspose_HYPRE; 2323 B->ops->multadd = MatMultAdd_HYPRE; 2324 B->ops->multtransposeadd = MatMultTransposeAdd_HYPRE; 2325 B->ops->setup = MatSetUp_HYPRE; 2326 B->ops->destroy = MatDestroy_HYPRE; 2327 B->ops->assemblyend = MatAssemblyEnd_HYPRE; 2328 B->ops->assemblybegin = MatAssemblyBegin_HYPRE; 2329 B->ops->setvalues = MatSetValues_HYPRE; 2330 B->ops->missingdiagonal = MatMissingDiagonal_HYPRE; 2331 B->ops->scale = MatScale_HYPRE; 2332 B->ops->zerorowscolumns = MatZeroRowsColumns_HYPRE; 2333 B->ops->zeroentries = MatZeroEntries_HYPRE; 2334 B->ops->zerorows = MatZeroRows_HYPRE; 2335 B->ops->getrow = MatGetRow_HYPRE; 2336 B->ops->restorerow = MatRestoreRow_HYPRE; 2337 B->ops->getvalues = MatGetValues_HYPRE; 2338 B->ops->setoption = MatSetOption_HYPRE; 2339 B->ops->duplicate = MatDuplicate_HYPRE; 2340 B->ops->copy = MatCopy_HYPRE; 2341 B->ops->view = MatView_HYPRE; 2342 B->ops->getdiagonal = MatGetDiagonal_HYPRE; 2343 B->ops->axpy = MatAXPY_HYPRE; 2344 B->ops->productsetfromoptions = MatProductSetFromOptions_HYPRE; 2345 #if defined(PETSC_HAVE_HYPRE_DEVICE) 2346 B->ops->bindtocpu = MatBindToCPU_HYPRE; 2347 B->boundtocpu = PETSC_FALSE; 2348 #endif 2349 2350 /* build cache for off array entries formed */ 2351 PetscCall(MatStashCreate_Private(PetscObjectComm((PetscObject)B), 1, &B->stash)); 2352 2353 PetscCall(PetscCommGetComm(PetscObjectComm((PetscObject)B), &hB->comm)); 2354 PetscCall(PetscObjectChangeTypeName((PetscObject)B, MATHYPRE)); 2355 PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatConvert_hypre_aij_C", MatConvert_HYPRE_AIJ)); 2356 PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatConvert_hypre_is_C", MatConvert_HYPRE_IS)); 2357 PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatProductSetFromOptions_seqaij_hypre_C", MatProductSetFromOptions_HYPRE)); 2358 PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatProductSetFromOptions_mpiaij_hypre_C", MatProductSetFromOptions_HYPRE)); 2359 PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatHYPRESetPreallocation_C", MatHYPRESetPreallocation_HYPRE)); 2360 PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatHYPREGetParCSR_C", MatHYPREGetParCSR_HYPRE)); 2361 PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatSetPreallocationCOO_C", MatSetPreallocationCOO_HYPRE)); 2362 PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatSetValuesCOO_C", MatSetValuesCOO_HYPRE)); 2363 #if defined(PETSC_HAVE_HYPRE_DEVICE) 2364 #if defined(HYPRE_USING_HIP) 2365 PetscCall(PetscDeviceInitialize(PETSC_DEVICE_HIP)); 2366 PetscCall(MatSetVecType(B, VECHIP)); 2367 #endif 2368 #if defined(HYPRE_USING_CUDA) 2369 PetscCall(PetscDeviceInitialize(PETSC_DEVICE_CUDA)); 2370 PetscCall(MatSetVecType(B, VECCUDA)); 2371 #endif 2372 #endif 2373 PetscFunctionReturn(PETSC_SUCCESS); 2374 } 2375 2376 static PetscErrorCode hypre_array_destroy(void *ptr) 2377 { 2378 PetscFunctionBegin; 2379 hypre_TFree(ptr, HYPRE_MEMORY_HOST); 2380 PetscFunctionReturn(PETSC_SUCCESS); 2381 } 2382