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