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 PetscStackCallStandard(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 PetscStackCallStandard(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 PetscStackCallStandard(HYPRE_IJMatrixCreate,hA->comm,rstart,rend-1,cstart,cend-1,&hA->ij); 97 PetscStackCallStandard(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 PetscStackCallStandard(HYPRE_IJMatrixInitialize,ij); 138 #else 139 PetscStackCallStandard(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 PetscStackCallStandard(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 PetscStackCallStandard(HYPRE_IJMatrixGetObjectType,ij,&type); 182 PetscCheck(type == HYPRE_PARCSR,PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"Only HYPRE_PARCSR is supported"); 183 PetscStackCallStandard(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 PetscStackCallStandard(HYPRE_IJMatrixGetObjectType,ij,&type); 226 PetscCheck(type == HYPRE_PARCSR,PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"Only HYPRE_PARCSR is supported"); 227 PetscStackCallStandard(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 PetscStackCallStandard(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 PetscStackCallStandard(HYPRE_IJMatrixGetObjectType,mhA->ij,&type); 296 PetscCheck(type == HYPRE_PARCSR,comm,PETSC_ERR_SUP,"Only HYPRE_PARCSR is supported"); 297 PetscStackCallStandard(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 PetscStackCallStandard(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 PetscStackCallStandard(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 PetscStackCallStandard(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 PetscStackCallStandard(hypre_ParCSRMatrixInitialize_v2,tA,iscuda ? HYPRE_MEMORY_DEVICE : HYPRE_MEMORY_HOST); 773 #else 774 #if PETSC_PKG_HYPRE_VERSION_LT(2,18,0) 775 PetscStackCallStandard(hypre_ParCSRMatrixInitialize,tA); 776 #else 777 PetscStackCallStandard(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)) PetscStackCallStandard(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 PetscStackPush("hypre_ParCSRMatrixRAP"); 849 *hRAP = hypre_ParCSRMatrixRAP(hR,hA,hP); 850 PetscStackPop; 851 #else 852 PetscStackCallStandard(hypre_BoomerAMGBuildCoarseOperator,hR,hA,hP,hRAP); 853 PetscStackCallStandard(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 PetscStackCallStandard(HYPRE_IJMatrixGetObjectType,hP->ij,&type); 908 PetscCheck(type == HYPRE_PARCSR,comm,PETSC_ERR_SUP,"Only HYPRE_PARCSR is supported"); 909 PetscStackCallStandard(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 PetscStackCallStandard(HYPRE_IJMatrixGetObjectType,hA->ij,&type); 937 PetscCheck(type == HYPRE_PARCSR,PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"Only HYPRE_PARCSR is supported"); 938 PetscStackCallStandard(HYPRE_IJMatrixGetObjectType,hP->ij,&type); 939 PetscCheck(type == HYPRE_PARCSR,PetscObjectComm((PetscObject)P),PETSC_ERR_SUP,"Only HYPRE_PARCSR is supported"); 940 PetscStackCallStandard(HYPRE_IJMatrixGetObject,hA->ij,(void**)&Aparcsr); 941 PetscStackCallStandard(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 PetscStackPush("hypre_ParCSRMatMat"); 958 *hAB = hypre_ParCSRMatMat(hA,hB); 959 #else 960 PetscStackPush("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 PetscStackCallStandard(HYPRE_IJMatrixGetObjectType,hA->ij,&type); 1013 PetscCheck(type == HYPRE_PARCSR,PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"Only HYPRE_PARCSR is supported"); 1014 PetscStackCallStandard(HYPRE_IJMatrixGetObjectType,hB->ij,&type); 1015 PetscCheck(type == HYPRE_PARCSR,PetscObjectComm((PetscObject)B),PETSC_ERR_SUP,"Only HYPRE_PARCSR is supported"); 1016 PetscStackCallStandard(HYPRE_IJMatrixGetObject,hA->ij,(void**)&Aparcsr); 1017 PetscStackCallStandard(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 PetscStackCallStandard(HYPRE_IJVectorGetObject,hA->b->ij,(void**)&hx); 1195 PetscStackCallStandard(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 PetscStackCallStandard(HYPRE_IJVectorGetObject,hA->x->ij,(void**)&hx); 1201 PetscStackCallStandard(HYPRE_IJVectorGetObject,hA->b->ij,(void**)&hy); 1202 } 1203 PetscStackCallStandard(HYPRE_IJMatrixGetObject,hA->ij,(void**)&parcsr); 1204 if (trans) { 1205 PetscStackCallStandard(hypre_ParCSRMatrixMatvecT,a,parcsr,hx,b,hy); 1206 } else { 1207 PetscStackCallStandard(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 PetscStackCallStandard(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 PetscStackCall("hypre_TFree",hypre_TFree(hA->diagJ,hA->memType)); 1233 PetscStackCall("hypre_TFree",hypre_TFree(hA->offdJ,hA->memType)); 1234 PetscStackCall("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 PetscStackCallStandard(HYPRE_IJMatrixGetObject,hA->ij,(void**)&parcsr); 1268 PetscStackCallStandard(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 PetscStackCallStandard(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 PetscStackCallStandard(hypre_AuxParCSRMatrixDestroy,aux_matrix); 1317 hypre_IJMatrixTranslator(hA->ij) = NULL; 1318 1319 /* Initialize with assembled flag -> it only recreates the aux_par_matrix */ 1320 PetscStackCallStandard(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 PetscStackCallStandard(hypre_AuxParCSRMatrixInitialize,aux_matrix); 1326 #else 1327 PetscStackCallStandard(hypre_AuxParCSRMatrixInitialize_v2,aux_matrix,HYPRE_MEMORY_HOST); 1328 #endif 1329 } 1330 } 1331 { 1332 hypre_ParCSRMatrix *parcsr; 1333 1334 PetscStackCallStandard(HYPRE_IJMatrixGetObject,hA->ij,(void**)&parcsr); 1335 if (!hypre_ParCSRMatrixCommPkg(parcsr)) PetscStackCallStandard(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 PetscStackCallStandard(HYPRE_IJMatrixGetObject,hA->ij,(void**)&parcsr); 1405 PetscStackCallStandard(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 PetscStackCallStandard(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 PetscStackCallStandard(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 PetscStackCallStandard(HYPRE_IJMatrixCreate,hA->comm,rs,re-1,cs,ce-1,&hA->ij); 1461 PetscStackCallStandard(HYPRE_IJMatrixSetObjectType,hA->ij,HYPRE_PARCSR); 1462 } else { 1463 HYPRE_BigInt hrs,hre,hcs,hce; 1464 PetscStackCallStandard(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 PetscStackCallStandard(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 PetscStackCallStandard(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 PetscStackCallStandard(HYPRE_IJMatrixInitialize,hA->ij); 1507 #else 1508 PetscStackCallStandard(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 PetscStackCallStandard(HYPRE_IJMatrixCreate,hA->comm,rstart,rend-1,cstart,cend-1,&hA->ij); 1623 PetscStackCallStandard(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 PetscStackCallStandard(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 PetscStackCallStandard(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 PetscStackCallStandard(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 PetscStackCallStandard(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 if (rows[i] < rst || rows[i] >= ren) 1840 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Non-local rows not yet supported"); 1841 lrows[i] = rows[i] - rst; 1842 } 1843 PetscStackCallStandard(hypre_ParCSRMatrixEliminateRowsCols,parcsr,numRows,lrows); 1844 PetscCall(PetscFree(lrows)); 1845 PetscFunctionReturn(0); 1846 } 1847 1848 static PetscErrorCode MatZeroEntries_HYPRE_CSRMatrix(hypre_CSRMatrix *ha) 1849 { 1850 PetscFunctionBegin; 1851 if (ha) { 1852 HYPRE_Int *ii, size; 1853 HYPRE_Complex *a; 1854 1855 size = hypre_CSRMatrixNumRows(ha); 1856 a = hypre_CSRMatrixData(ha); 1857 ii = hypre_CSRMatrixI(ha); 1858 1859 if (a) PetscCall(PetscArrayzero(a,ii[size])); 1860 } 1861 PetscFunctionReturn(0); 1862 } 1863 1864 PetscErrorCode MatZeroEntries_HYPRE(Mat A) 1865 { 1866 Mat_HYPRE *hA = (Mat_HYPRE*)A->data; 1867 1868 PetscFunctionBegin; 1869 if (HYPRE_MEMORY_DEVICE == hypre_IJMatrixMemoryLocation(hA->ij)) { 1870 PetscStackCallStandard(HYPRE_IJMatrixSetConstantValues,hA->ij,0.0); 1871 } else { 1872 hypre_ParCSRMatrix *parcsr; 1873 1874 PetscCall(MatHYPREGetParCSR_HYPRE(A,&parcsr)); 1875 PetscCall(MatZeroEntries_HYPRE_CSRMatrix(hypre_ParCSRMatrixDiag(parcsr))); 1876 PetscCall(MatZeroEntries_HYPRE_CSRMatrix(hypre_ParCSRMatrixOffd(parcsr))); 1877 } 1878 PetscFunctionReturn(0); 1879 } 1880 1881 static PetscErrorCode MatZeroRows_HYPRE_CSRMatrix(hypre_CSRMatrix *hA,PetscInt N,const PetscInt rows[],HYPRE_Complex diag) 1882 { 1883 PetscInt ii; 1884 HYPRE_Int *i, *j; 1885 HYPRE_Complex *a; 1886 1887 PetscFunctionBegin; 1888 if (!hA) PetscFunctionReturn(0); 1889 1890 i = hypre_CSRMatrixI(hA); 1891 j = hypre_CSRMatrixJ(hA); 1892 a = hypre_CSRMatrixData(hA); 1893 1894 for (ii = 0; ii < N; ii++) { 1895 HYPRE_Int jj, ibeg, iend, irow; 1896 1897 irow = rows[ii]; 1898 ibeg = i[irow]; 1899 iend = i[irow+1]; 1900 for (jj = ibeg; jj < iend; jj++) 1901 if (j[jj] == irow) a[jj] = diag; 1902 else a[jj] = 0.0; 1903 } 1904 PetscFunctionReturn(0); 1905 } 1906 1907 static PetscErrorCode MatZeroRows_HYPRE(Mat A,PetscInt N,const PetscInt rows[],PetscScalar diag,Vec x,Vec b) 1908 { 1909 hypre_ParCSRMatrix *parcsr; 1910 PetscInt *lrows,len; 1911 HYPRE_Complex hdiag; 1912 1913 PetscFunctionBegin; 1914 PetscCheck(!x && !b,PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Does not support to modify the solution and the right hand size"); 1915 PetscCall(PetscHYPREScalarCast(diag,&hdiag)); 1916 /* retrieve the internal matrix */ 1917 PetscCall(MatHYPREGetParCSR_HYPRE(A,&parcsr)); 1918 /* get locally owned rows */ 1919 PetscCall(MatZeroRowsMapLocal_Private(A,N,rows,&len,&lrows)); 1920 /* zero diagonal part */ 1921 PetscCall(MatZeroRows_HYPRE_CSRMatrix(hypre_ParCSRMatrixDiag(parcsr),len,lrows,hdiag)); 1922 /* zero off-diagonal part */ 1923 PetscCall(MatZeroRows_HYPRE_CSRMatrix(hypre_ParCSRMatrixOffd(parcsr),len,lrows,0.0)); 1924 1925 PetscCall(PetscFree(lrows)); 1926 PetscFunctionReturn(0); 1927 } 1928 1929 static PetscErrorCode MatAssemblyBegin_HYPRE(Mat mat,MatAssemblyType mode) 1930 { 1931 PetscFunctionBegin; 1932 if (mat->nooffprocentries) PetscFunctionReturn(0); 1933 1934 PetscCall(MatStashScatterBegin_Private(mat,&mat->stash,mat->rmap->range)); 1935 PetscFunctionReturn(0); 1936 } 1937 1938 static PetscErrorCode MatGetRow_HYPRE(Mat A,PetscInt row,PetscInt *nz,PetscInt **idx,PetscScalar **v) 1939 { 1940 hypre_ParCSRMatrix *parcsr; 1941 HYPRE_Int hnz; 1942 1943 PetscFunctionBegin; 1944 /* retrieve the internal matrix */ 1945 PetscCall(MatHYPREGetParCSR_HYPRE(A,&parcsr)); 1946 /* call HYPRE API */ 1947 PetscStackCallStandard(HYPRE_ParCSRMatrixGetRow,parcsr,row,&hnz,(HYPRE_BigInt**)idx,(HYPRE_Complex**)v); 1948 if (nz) *nz = (PetscInt)hnz; 1949 PetscFunctionReturn(0); 1950 } 1951 1952 static PetscErrorCode MatRestoreRow_HYPRE(Mat A,PetscInt row,PetscInt *nz,PetscInt **idx,PetscScalar **v) 1953 { 1954 hypre_ParCSRMatrix *parcsr; 1955 HYPRE_Int hnz; 1956 1957 PetscFunctionBegin; 1958 /* retrieve the internal matrix */ 1959 PetscCall(MatHYPREGetParCSR_HYPRE(A,&parcsr)); 1960 /* call HYPRE API */ 1961 hnz = nz ? (HYPRE_Int)(*nz) : 0; 1962 PetscStackCallStandard(HYPRE_ParCSRMatrixRestoreRow,parcsr,row,&hnz,(HYPRE_BigInt**)idx,(HYPRE_Complex**)v); 1963 PetscFunctionReturn(0); 1964 } 1965 1966 static PetscErrorCode MatGetValues_HYPRE(Mat A,PetscInt m,const PetscInt idxm[],PetscInt n,const PetscInt idxn[],PetscScalar v[]) 1967 { 1968 Mat_HYPRE *hA = (Mat_HYPRE*)A->data; 1969 PetscInt i; 1970 1971 PetscFunctionBegin; 1972 if (!m || !n) PetscFunctionReturn(0); 1973 /* Ignore negative row indices 1974 * And negative column indices should be automatically ignored in hypre 1975 * */ 1976 for (i=0; i<m; i++) { 1977 if (idxm[i] >= 0) { 1978 HYPRE_Int hn = (HYPRE_Int)n; 1979 PetscStackCallStandard(HYPRE_IJMatrixGetValues,hA->ij,1,&hn,(HYPRE_BigInt*)&idxm[i],(HYPRE_BigInt*)idxn,(HYPRE_Complex*)(v + i*n)); 1980 } 1981 } 1982 PetscFunctionReturn(0); 1983 } 1984 1985 static PetscErrorCode MatSetOption_HYPRE(Mat A,MatOption op,PetscBool flg) 1986 { 1987 Mat_HYPRE *hA = (Mat_HYPRE*)A->data; 1988 1989 PetscFunctionBegin; 1990 switch (op) { 1991 case MAT_NO_OFF_PROC_ENTRIES: 1992 if (flg) { 1993 PetscStackCallStandard(HYPRE_IJMatrixSetMaxOffProcElmts,hA->ij,0); 1994 } 1995 break; 1996 case MAT_SORTED_FULL: 1997 hA->sorted_full = flg; 1998 break; 1999 default: 2000 break; 2001 } 2002 PetscFunctionReturn(0); 2003 } 2004 2005 static PetscErrorCode MatView_HYPRE(Mat A, PetscViewer view) 2006 { 2007 PetscViewerFormat format; 2008 2009 PetscFunctionBegin; 2010 PetscCall(PetscViewerGetFormat(view,&format)); 2011 if (format == PETSC_VIEWER_ASCII_FACTOR_INFO || format == PETSC_VIEWER_ASCII_INFO || format == PETSC_VIEWER_ASCII_INFO_DETAIL) PetscFunctionReturn(0); 2012 if (format != PETSC_VIEWER_NATIVE) { 2013 Mat B; 2014 hypre_ParCSRMatrix *parcsr; 2015 PetscErrorCode (*mview)(Mat,PetscViewer) = NULL; 2016 2017 PetscCall(MatHYPREGetParCSR_HYPRE(A,&parcsr)); 2018 PetscCall(MatCreateFromParCSR(parcsr,MATAIJ,PETSC_USE_POINTER,&B)); 2019 PetscCall(MatGetOperation(B,MATOP_VIEW,(void(**)(void))&mview)); 2020 PetscCheck(mview,PetscObjectComm((PetscObject)A),PETSC_ERR_PLIB,"Missing view operation"); 2021 PetscCall((*mview)(B,view)); 2022 PetscCall(MatDestroy(&B)); 2023 } else { 2024 Mat_HYPRE *hA = (Mat_HYPRE*)A->data; 2025 PetscMPIInt size; 2026 PetscBool isascii; 2027 const char *filename; 2028 2029 /* HYPRE uses only text files */ 2030 PetscCall(PetscObjectTypeCompare((PetscObject)view,PETSCVIEWERASCII,&isascii)); 2031 PetscCheck(isascii,PetscObjectComm((PetscObject)view),PETSC_ERR_SUP,"PetscViewerType %s: native HYPRE format needs PETSCVIEWERASCII",((PetscObject)view)->type_name); 2032 PetscCall(PetscViewerFileGetName(view,&filename)); 2033 PetscStackCallStandard(HYPRE_IJMatrixPrint,hA->ij,filename); 2034 PetscCallMPI(MPI_Comm_size(hA->comm,&size)); 2035 if (size > 1) { 2036 PetscCall(PetscViewerASCIIPrintf(view,"Matrix files: %s.%05d ... %s.%05d\n",filename,0,filename,size-1)); 2037 } else { 2038 PetscCall(PetscViewerASCIIPrintf(view,"Matrix file: %s.%05d\n",filename,0)); 2039 } 2040 } 2041 PetscFunctionReturn(0); 2042 } 2043 2044 static PetscErrorCode MatDuplicate_HYPRE(Mat A,MatDuplicateOption op, Mat *B) 2045 { 2046 hypre_ParCSRMatrix *parcsr = NULL; 2047 PetscCopyMode cpmode; 2048 2049 PetscFunctionBegin; 2050 PetscCall(MatHYPREGetParCSR_HYPRE(A,&parcsr)); 2051 if (op == MAT_DO_NOT_COPY_VALUES || op == MAT_SHARE_NONZERO_PATTERN) { 2052 parcsr = hypre_ParCSRMatrixClone(parcsr,0); 2053 cpmode = PETSC_OWN_POINTER; 2054 } else { 2055 cpmode = PETSC_COPY_VALUES; 2056 } 2057 PetscCall(MatCreateFromParCSR(parcsr,MATHYPRE,cpmode,B)); 2058 PetscFunctionReturn(0); 2059 } 2060 2061 static PetscErrorCode MatCopy_HYPRE(Mat A, Mat B, MatStructure str) 2062 { 2063 hypre_ParCSRMatrix *acsr,*bcsr; 2064 2065 PetscFunctionBegin; 2066 if (str == SAME_NONZERO_PATTERN && A->ops->copy == B->ops->copy) { 2067 PetscCall(MatHYPREGetParCSR_HYPRE(A,&acsr)); 2068 PetscCall(MatHYPREGetParCSR_HYPRE(B,&bcsr)); 2069 PetscStackCallStandard(hypre_ParCSRMatrixCopy,acsr,bcsr,1); 2070 PetscCall(MatSetOption(B,MAT_SORTED_FULL,PETSC_TRUE)); /* "perfect" preallocation, so no need for hypre_AuxParCSRMatrixNeedAux */ 2071 PetscCall(MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY)); 2072 PetscCall(MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY)); 2073 } else { 2074 PetscCall(MatCopy_Basic(A,B,str)); 2075 } 2076 PetscFunctionReturn(0); 2077 } 2078 2079 static PetscErrorCode MatGetDiagonal_HYPRE(Mat A, Vec d) 2080 { 2081 hypre_ParCSRMatrix *parcsr; 2082 hypre_CSRMatrix *dmat; 2083 HYPRE_Complex *a; 2084 HYPRE_Complex *data = NULL; 2085 HYPRE_Int *diag = NULL; 2086 PetscInt i; 2087 PetscBool cong; 2088 2089 PetscFunctionBegin; 2090 PetscCall(MatHasCongruentLayouts(A,&cong)); 2091 PetscCheck(cong,PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"Only for square matrices with same local distributions of rows and columns"); 2092 if (PetscDefined(USE_DEBUG)) { 2093 PetscBool miss; 2094 PetscCall(MatMissingDiagonal(A,&miss,NULL)); 2095 PetscCheck(!miss || !A->rmap->n,PETSC_COMM_SELF,PETSC_ERR_SUP,"Not implemented when diagonal entries are missing"); 2096 } 2097 PetscCall(MatHYPREGetParCSR_HYPRE(A,&parcsr)); 2098 dmat = hypre_ParCSRMatrixDiag(parcsr); 2099 if (dmat) { 2100 /* this cast fixes the clang error: implicit conversion from 'HYPRE_Complex' (aka '_Complex double') to 'double' is not permitted in C++ */ 2101 PetscCall(VecGetArray(d,(PetscScalar**)&a)); 2102 diag = hypre_CSRMatrixI(dmat); 2103 data = hypre_CSRMatrixData(dmat); 2104 for (i=0;i<A->rmap->n;i++) a[i] = data[diag[i]]; 2105 PetscCall(VecRestoreArray(d,(PetscScalar**)&a)); 2106 } 2107 PetscFunctionReturn(0); 2108 } 2109 2110 #include <petscblaslapack.h> 2111 2112 static PetscErrorCode MatAXPY_HYPRE(Mat Y,PetscScalar a,Mat X,MatStructure str) 2113 { 2114 PetscFunctionBegin; 2115 #if defined(PETSC_HAVE_HYPRE_DEVICE) 2116 { 2117 Mat B; 2118 hypre_ParCSRMatrix *x,*y,*z; 2119 2120 PetscCall(MatHYPREGetParCSR(Y,&y)); 2121 PetscCall(MatHYPREGetParCSR(X,&x)); 2122 PetscStackCallStandard(hypre_ParCSRMatrixAdd,1.0,y,1.0,x,&z); 2123 PetscCall(MatCreateFromParCSR(z,MATHYPRE,PETSC_OWN_POINTER,&B)); 2124 PetscCall(MatHeaderMerge(Y,&B)); 2125 } 2126 #else 2127 if (str == SAME_NONZERO_PATTERN) { 2128 hypre_ParCSRMatrix *x,*y; 2129 hypre_CSRMatrix *xloc,*yloc; 2130 PetscInt xnnz,ynnz; 2131 HYPRE_Complex *xarr,*yarr; 2132 PetscBLASInt one=1,bnz; 2133 2134 PetscCall(MatHYPREGetParCSR(Y,&y)); 2135 PetscCall(MatHYPREGetParCSR(X,&x)); 2136 2137 /* diagonal block */ 2138 xloc = hypre_ParCSRMatrixDiag(x); 2139 yloc = hypre_ParCSRMatrixDiag(y); 2140 xnnz = 0; 2141 ynnz = 0; 2142 xarr = NULL; 2143 yarr = NULL; 2144 if (xloc) { 2145 xarr = hypre_CSRMatrixData(xloc); 2146 xnnz = hypre_CSRMatrixNumNonzeros(xloc); 2147 } 2148 if (yloc) { 2149 yarr = hypre_CSRMatrixData(yloc); 2150 ynnz = hypre_CSRMatrixNumNonzeros(yloc); 2151 } 2152 PetscCheck(xnnz == ynnz,PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Different number of nonzeros in diagonal block %" PetscInt_FMT " != %" PetscInt_FMT,xnnz,ynnz); 2153 PetscCall(PetscBLASIntCast(xnnz,&bnz)); 2154 PetscStackCallBLAS("BLASaxpy",BLASaxpy_(&bnz,&a,(PetscScalar*)xarr,&one,(PetscScalar*)yarr,&one)); 2155 2156 /* off-diagonal block */ 2157 xloc = hypre_ParCSRMatrixOffd(x); 2158 yloc = hypre_ParCSRMatrixOffd(y); 2159 xnnz = 0; 2160 ynnz = 0; 2161 xarr = NULL; 2162 yarr = NULL; 2163 if (xloc) { 2164 xarr = hypre_CSRMatrixData(xloc); 2165 xnnz = hypre_CSRMatrixNumNonzeros(xloc); 2166 } 2167 if (yloc) { 2168 yarr = hypre_CSRMatrixData(yloc); 2169 ynnz = hypre_CSRMatrixNumNonzeros(yloc); 2170 } 2171 PetscCheck(xnnz == ynnz,PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Different number of nonzeros in off-diagonal block %" PetscInt_FMT " != %" PetscInt_FMT,xnnz,ynnz); 2172 PetscCall(PetscBLASIntCast(xnnz,&bnz)); 2173 PetscStackCallBLAS("BLASaxpy",BLASaxpy_(&bnz,&a,(PetscScalar*)xarr,&one,(PetscScalar*)yarr,&one)); 2174 } else if (str == SUBSET_NONZERO_PATTERN) { 2175 PetscCall(MatAXPY_Basic(Y,a,X,str)); 2176 } else { 2177 Mat B; 2178 2179 PetscCall(MatAXPY_Basic_Preallocate(Y,X,&B)); 2180 PetscCall(MatAXPY_BasicWithPreallocation(B,Y,a,X,str)); 2181 PetscCall(MatHeaderReplace(Y,&B)); 2182 } 2183 #endif 2184 PetscFunctionReturn(0); 2185 } 2186 2187 static PetscErrorCode MatSetPreallocationCOO_HYPRE(Mat mat, PetscCount coo_n, PetscInt coo_i[], PetscInt coo_j[]) 2188 { 2189 MPI_Comm comm; 2190 PetscMPIInt size; 2191 PetscLayout rmap,cmap; 2192 Mat_HYPRE *hmat; 2193 hypre_ParCSRMatrix *parCSR; 2194 hypre_CSRMatrix *diag,*offd; 2195 Mat A,B,cooMat; 2196 PetscScalar *Aa,*Ba; 2197 HYPRE_MemoryLocation hypreMemtype = HYPRE_MEMORY_HOST; 2198 PetscMemType petscMemtype; 2199 MatType matType = MATAIJ; /* default type of cooMat */ 2200 2201 PetscFunctionBegin; 2202 /* Build an agent matrix cooMat whose type is either MATAIJ or MATAIJKOKKOS. 2203 It has the same sparsity pattern as mat, and also shares the data array with mat. We use cooMat to do the COO work. 2204 */ 2205 PetscCall(PetscObjectGetComm((PetscObject)mat,&comm)); 2206 PetscCallMPI(MPI_Comm_size(comm,&size)); 2207 PetscCall(PetscLayoutSetUp(mat->rmap)); 2208 PetscCall(PetscLayoutSetUp(mat->cmap)); 2209 PetscCall(MatGetLayouts(mat,&rmap,&cmap)); 2210 2211 /* I do not know how hypre_ParCSRMatrix stores diagonal elements for non-square matrices, so I just give up now */ 2212 PetscCheck(rmap->N == cmap->N,comm,PETSC_ERR_SUP,"MATHYPRE COO cannot handle non-square matrices"); 2213 2214 #if defined(PETSC_HAVE_DEVICE) 2215 if (!mat->boundtocpu) { /* mat will be on device, so will cooMat */ 2216 #if defined(PETSC_HAVE_KOKKOS) 2217 matType = MATAIJKOKKOS; 2218 #else 2219 SETERRQ(comm,PETSC_ERR_SUP,"To support MATHYPRE COO assembly on device, we need Kokkos, e.g., --download-kokkos --download-kokkos-kernels"); 2220 #endif 2221 } 2222 #endif 2223 2224 /* Do COO preallocation through cooMat */ 2225 hmat = (Mat_HYPRE*)mat->data; 2226 PetscCall(MatDestroy(&hmat->cooMat)); 2227 PetscCall(MatCreate(comm,&cooMat)); 2228 PetscCall(MatSetType(cooMat,matType)); 2229 PetscCall(MatSetLayouts(cooMat,rmap,cmap)); 2230 PetscCall(MatSetPreallocationCOO(cooMat,coo_n,coo_i,coo_j)); 2231 2232 /* Copy the sparsity pattern from cooMat to hypre IJMatrix hmat->ij */ 2233 PetscCall(MatSetOption(mat,MAT_SORTED_FULL,PETSC_TRUE)); 2234 PetscCall(MatSetOption(mat,MAT_NO_OFF_PROC_ENTRIES,PETSC_TRUE)); 2235 PetscCall(MatHYPRE_CreateFromMat(cooMat,hmat)); /* Create hmat->ij and preallocate it */ 2236 PetscCall(MatHYPRE_IJMatrixCopy(cooMat,hmat->ij)); /* Copy A's (a,i,j) to hmat->ij. To reuse code. Copying 'a' is not really needed */ 2237 2238 mat->preallocated = PETSC_TRUE; 2239 PetscCall(MatAssemblyBegin(mat,MAT_FINAL_ASSEMBLY)); 2240 PetscCall(MatAssemblyEnd(mat,MAT_FINAL_ASSEMBLY)); /* Migrate mat to device if it is bound to. Hypre builds its own SpMV context here */ 2241 2242 /* Alias cooMat's data array to IJMatrix's */ 2243 PetscStackCallStandard(HYPRE_IJMatrixGetObject,hmat->ij,(void**)&parCSR); 2244 diag = hypre_ParCSRMatrixDiag(parCSR); 2245 offd = hypre_ParCSRMatrixOffd(parCSR); 2246 2247 hypreMemtype = hypre_CSRMatrixMemoryLocation(diag); 2248 A = (size == 1)? cooMat : ((Mat_MPIAIJ*)cooMat->data)->A; 2249 PetscCall(MatSeqAIJGetCSRAndMemType(A,NULL,NULL,&Aa,&petscMemtype)); 2250 PetscAssert((PetscMemTypeHost(petscMemtype) && hypreMemtype == HYPRE_MEMORY_HOST) || 2251 (PetscMemTypeDevice(petscMemtype) && hypreMemtype == HYPRE_MEMORY_DEVICE), 2252 comm,PETSC_ERR_PLIB,"PETSc and hypre's memory types mismatch"); 2253 2254 hmat->diagJ = hypre_CSRMatrixJ(diag); 2255 PetscStackCall("hypre_TFree",hypre_TFree(hypre_CSRMatrixData(diag),hypreMemtype)); 2256 hypre_CSRMatrixData(diag) = (HYPRE_Complex*)Aa; 2257 hypre_CSRMatrixOwnsData(diag) = 0; /* Take ownership of (j,a) away from hypre. As a result, we need to free them on our own */ 2258 2259 /* Copy diagonal pointers of A to device to facilitate MatSeqAIJMoveDiagonalValuesFront_SeqAIJKokkos */ 2260 if (hypreMemtype == HYPRE_MEMORY_DEVICE) { 2261 PetscStackCall("hypre_TAlloc",hmat->diag = hypre_TAlloc(PetscInt,rmap->n,hypreMemtype)); 2262 PetscCall(MatMarkDiagonal_SeqAIJ(A)); /* We need updated diagonal positions */ 2263 PetscStackCall("hypre_TMemcpy",hypre_TMemcpy(hmat->diag,((Mat_SeqAIJ*)A->data)->diag,PetscInt,rmap->n,hypreMemtype,HYPRE_MEMORY_HOST)); 2264 } 2265 2266 if (size > 1) { 2267 B = ((Mat_MPIAIJ*)cooMat->data)->B; 2268 PetscCall(MatSeqAIJGetCSRAndMemType(B,NULL,NULL,&Ba,&petscMemtype)); 2269 hmat->offdJ = hypre_CSRMatrixJ(offd); 2270 PetscStackCall("hypre_TFree",hypre_TFree(hypre_CSRMatrixData(offd),hypreMemtype)); 2271 hypre_CSRMatrixData(offd) = (HYPRE_Complex*)Ba; 2272 hypre_CSRMatrixOwnsData(offd) = 0; 2273 } 2274 2275 /* Record cooMat for use in MatSetValuesCOO_HYPRE */ 2276 hmat->cooMat = cooMat; 2277 hmat->memType = hypreMemtype; 2278 PetscFunctionReturn(0); 2279 } 2280 2281 static PetscErrorCode MatSetValuesCOO_HYPRE(Mat mat, const PetscScalar v[], InsertMode imode) 2282 { 2283 Mat_HYPRE *hmat = (Mat_HYPRE*)mat->data; 2284 PetscMPIInt size; 2285 Mat A; 2286 2287 PetscFunctionBegin; 2288 PetscCheck(hmat->cooMat,hmat->comm,PETSC_ERR_PLIB,"HYPRE COO delegate matrix has not been created yet"); 2289 PetscCallMPI(MPI_Comm_size(hmat->comm,&size)); 2290 PetscCall(MatSetValuesCOO(hmat->cooMat,v,imode)); 2291 2292 /* Move diagonal elements of the diagonal block to the front of their row, as needed by ParCSRMatrix. So damn hacky */ 2293 A = (size == 1) ? hmat->cooMat : ((Mat_MPIAIJ*)hmat->cooMat->data)->A; 2294 if (hmat->memType == HYPRE_MEMORY_HOST) { 2295 Mat_SeqAIJ *aij = (Mat_SeqAIJ*)A->data; 2296 PetscInt i,m,*Ai = aij->i,*Adiag = aij->diag; 2297 PetscScalar *Aa = aij->a,tmp; 2298 2299 PetscCall(MatGetSize(A,&m,NULL)); 2300 for (i=0; i<m; i++) { 2301 if (Adiag[i] >= Ai[i] && Adiag[i] < Ai[i+1]) { /* Digonal element of this row exists in a[] and j[] */ 2302 tmp = Aa[Ai[i]]; 2303 Aa[Ai[i]] = Aa[Adiag[i]]; 2304 Aa[Adiag[i]] = tmp; 2305 } 2306 } 2307 } else { 2308 #if defined(PETSC_HAVE_KOKKOS_KERNELS) 2309 PetscCall(MatSeqAIJMoveDiagonalValuesFront_SeqAIJKokkos(A,hmat->diag)); 2310 #endif 2311 } 2312 PetscFunctionReturn(0); 2313 } 2314 2315 /*MC 2316 MATHYPRE - MATHYPRE = "hypre" - A matrix type to be used for sequential and parallel sparse matrices 2317 based on the hypre IJ interface. 2318 2319 Level: intermediate 2320 2321 .seealso: `MatCreate()` 2322 M*/ 2323 2324 PETSC_EXTERN PetscErrorCode MatCreate_HYPRE(Mat B) 2325 { 2326 Mat_HYPRE *hB; 2327 2328 PetscFunctionBegin; 2329 PetscCall(PetscNewLog(B,&hB)); 2330 2331 hB->inner_free = PETSC_TRUE; 2332 hB->available = PETSC_TRUE; 2333 hB->sorted_full = PETSC_FALSE; /* no assumption whether column indices are sorted or not */ 2334 hB->size = 0; 2335 hB->array = NULL; 2336 2337 B->data = (void*)hB; 2338 B->assembled = PETSC_FALSE; 2339 2340 PetscCall(PetscMemzero(B->ops,sizeof(struct _MatOps))); 2341 B->ops->mult = MatMult_HYPRE; 2342 B->ops->multtranspose = MatMultTranspose_HYPRE; 2343 B->ops->multadd = MatMultAdd_HYPRE; 2344 B->ops->multtransposeadd = MatMultTransposeAdd_HYPRE; 2345 B->ops->setup = MatSetUp_HYPRE; 2346 B->ops->destroy = MatDestroy_HYPRE; 2347 B->ops->assemblyend = MatAssemblyEnd_HYPRE; 2348 B->ops->assemblybegin = MatAssemblyBegin_HYPRE; 2349 B->ops->setvalues = MatSetValues_HYPRE; 2350 B->ops->missingdiagonal = MatMissingDiagonal_HYPRE; 2351 B->ops->scale = MatScale_HYPRE; 2352 B->ops->zerorowscolumns = MatZeroRowsColumns_HYPRE; 2353 B->ops->zeroentries = MatZeroEntries_HYPRE; 2354 B->ops->zerorows = MatZeroRows_HYPRE; 2355 B->ops->getrow = MatGetRow_HYPRE; 2356 B->ops->restorerow = MatRestoreRow_HYPRE; 2357 B->ops->getvalues = MatGetValues_HYPRE; 2358 B->ops->setoption = MatSetOption_HYPRE; 2359 B->ops->duplicate = MatDuplicate_HYPRE; 2360 B->ops->copy = MatCopy_HYPRE; 2361 B->ops->view = MatView_HYPRE; 2362 B->ops->getdiagonal = MatGetDiagonal_HYPRE; 2363 B->ops->axpy = MatAXPY_HYPRE; 2364 B->ops->productsetfromoptions = MatProductSetFromOptions_HYPRE; 2365 #if defined(PETSC_HAVE_HYPRE_DEVICE) 2366 B->ops->bindtocpu = MatBindToCPU_HYPRE; 2367 B->boundtocpu = PETSC_FALSE; 2368 #endif 2369 2370 /* build cache for off array entries formed */ 2371 PetscCall(MatStashCreate_Private(PetscObjectComm((PetscObject)B),1,&B->stash)); 2372 2373 PetscCall(PetscCommGetComm(PetscObjectComm((PetscObject)B),&hB->comm)); 2374 PetscCall(PetscObjectChangeTypeName((PetscObject)B,MATHYPRE)); 2375 PetscCall(PetscObjectComposeFunction((PetscObject)B,"MatConvert_hypre_aij_C",MatConvert_HYPRE_AIJ)); 2376 PetscCall(PetscObjectComposeFunction((PetscObject)B,"MatConvert_hypre_is_C",MatConvert_HYPRE_IS)); 2377 PetscCall(PetscObjectComposeFunction((PetscObject)B,"MatProductSetFromOptions_seqaij_hypre_C",MatProductSetFromOptions_HYPRE)); 2378 PetscCall(PetscObjectComposeFunction((PetscObject)B,"MatProductSetFromOptions_mpiaij_hypre_C",MatProductSetFromOptions_HYPRE)); 2379 PetscCall(PetscObjectComposeFunction((PetscObject)B,"MatHYPRESetPreallocation_C",MatHYPRESetPreallocation_HYPRE)); 2380 PetscCall(PetscObjectComposeFunction((PetscObject)B,"MatHYPREGetParCSR_C",MatHYPREGetParCSR_HYPRE)); 2381 PetscCall(PetscObjectComposeFunction((PetscObject)B,"MatSetPreallocationCOO_C",MatSetPreallocationCOO_HYPRE)); 2382 PetscCall(PetscObjectComposeFunction((PetscObject)B,"MatSetValuesCOO_C",MatSetValuesCOO_HYPRE)); 2383 #if defined(PETSC_HAVE_HYPRE_DEVICE) 2384 #if defined(HYPRE_USING_HIP) 2385 PetscCall(PetscDeviceInitialize(PETSC_DEVICE_HIP)); 2386 PetscCall(MatSetVecType(B,VECHIP)); 2387 #endif 2388 #if defined(HYPRE_USING_CUDA) 2389 PetscCall(PetscDeviceInitialize(PETSC_DEVICE_CUDA)); 2390 PetscCall(MatSetVecType(B,VECCUDA)); 2391 #endif 2392 #endif 2393 PetscFunctionReturn(0); 2394 } 2395 2396 static PetscErrorCode hypre_array_destroy(void *ptr) 2397 { 2398 PetscFunctionBegin; 2399 hypre_TFree(ptr,HYPRE_MEMORY_HOST); 2400 PetscFunctionReturn(0); 2401 } 2402