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