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