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