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