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