1 #define PETSCMAT_DLL 2 3 /* 4 Creates hypre ijmatrix from PETSc matrix 5 */ 6 #include "petscsys.h" 7 #include "private/matimpl.h" /*I "petscmat.h" I*/ 8 #include "petscdm.h" /*I "petscdm.h" I*/ 9 #include "../src/dm/impls/da/hypre/mhyp.h" 10 11 #undef __FUNCT__ 12 #define __FUNCT__ "MatHYPRE_IJMatrixPreallocate" 13 PetscErrorCode MatHYPRE_IJMatrixPreallocate(Mat A_d, Mat A_o,HYPRE_IJMatrix ij) 14 { 15 PetscErrorCode ierr; 16 PetscInt i; 17 PetscInt n_d,*ia_d,n_o,*ia_o; 18 PetscBool done_d=PETSC_FALSE,done_o=PETSC_FALSE; 19 PetscInt *nnz_d=PETSC_NULL,*nnz_o=PETSC_NULL; 20 21 PetscFunctionBegin; 22 if (A_d) { /* determine number of nonzero entries in local diagonal part */ 23 ierr = MatGetRowIJ(A_d,0,PETSC_FALSE,PETSC_FALSE,&n_d,&ia_d,PETSC_NULL,&done_d);CHKERRQ(ierr); 24 if (done_d) { 25 ierr = PetscMalloc(n_d*sizeof(PetscInt),&nnz_d);CHKERRQ(ierr); 26 for (i=0; i<n_d; i++) { 27 nnz_d[i] = ia_d[i+1] - ia_d[i]; 28 } 29 } 30 ierr = MatRestoreRowIJ(A_d,0,PETSC_FALSE,PETSC_FALSE,&n_d,&ia_d,PETSC_NULL,&done_d);CHKERRQ(ierr); 31 } 32 if (A_o) { /* determine number of nonzero entries in local off-diagonal part */ 33 ierr = MatGetRowIJ(A_o,0,PETSC_FALSE,PETSC_FALSE,&n_o,&ia_o,PETSC_NULL,&done_o);CHKERRQ(ierr); 34 if (done_o) { 35 ierr = PetscMalloc(n_o*sizeof(PetscInt),&nnz_o);CHKERRQ(ierr); 36 for (i=0; i<n_o; i++) { 37 nnz_o[i] = ia_o[i+1] - ia_o[i]; 38 } 39 } 40 ierr = MatRestoreRowIJ(A_o,0,PETSC_FALSE,PETSC_FALSE,&n_o,&ia_o,PETSC_NULL,&done_o);CHKERRQ(ierr); 41 } 42 if (done_d) { /* set number of nonzeros in HYPRE IJ matrix */ 43 if (!done_o) { /* only diagonal part */ 44 ierr = PetscMalloc(n_d*sizeof(PetscInt),&nnz_o);CHKERRQ(ierr); 45 for (i=0; i<n_d; i++) { 46 nnz_o[i] = 0; 47 } 48 } 49 PetscStackCallHypre(0,HYPRE_IJMatrixSetDiagOffdSizes,(ij,nnz_d,nnz_o)); 50 ierr = PetscFree(nnz_d);CHKERRQ(ierr); 51 ierr = PetscFree(nnz_o);CHKERRQ(ierr); 52 } 53 PetscFunctionReturn(0); 54 } 55 56 #undef __FUNCT__ 57 #define __FUNCT__ "MatHYPRE_IJMatrixCreate" 58 PetscErrorCode MatHYPRE_IJMatrixCreate(Mat A,HYPRE_IJMatrix *ij) 59 { 60 PetscErrorCode ierr; 61 int rstart,rend,cstart,cend; 62 63 PetscFunctionBegin; 64 PetscValidHeaderSpecific(A,MAT_CLASSID,1); 65 PetscValidType(A,1); 66 PetscValidPointer(ij,2); 67 ierr = MatPreallocated(A);CHKERRQ(ierr); 68 rstart = A->rmap->rstart; 69 rend = A->rmap->rend; 70 cstart = A->cmap->rstart; 71 cend = A->cmap->rend; 72 PetscStackCallHypre(0,HYPRE_IJMatrixCreate,(((PetscObject)A)->comm,rstart,rend-1,cstart,cend-1,ij)); 73 PetscStackCallHypre(0,HYPRE_IJMatrixSetObjectType,(*ij,HYPRE_PARCSR)); 74 { 75 PetscBool same; 76 Mat A_d,A_o; 77 PetscInt *colmap; 78 ierr = PetscTypeCompare((PetscObject)A,MATMPIAIJ,&same);CHKERRQ(ierr); 79 if (same) { 80 ierr = MatMPIAIJGetSeqAIJ(A,&A_d,&A_o,&colmap);CHKERRQ(ierr); 81 ierr = MatHYPRE_IJMatrixPreallocate(A_d,A_o,*ij);CHKERRQ(ierr); 82 PetscFunctionReturn(0); 83 } 84 ierr = PetscTypeCompare((PetscObject)A,MATMPIBAIJ,&same);CHKERRQ(ierr); 85 if (same) { 86 ierr = MatMPIBAIJGetSeqBAIJ(A,&A_d,&A_o,&colmap);CHKERRQ(ierr); 87 ierr = MatHYPRE_IJMatrixPreallocate(A_d,A_o,*ij);CHKERRQ(ierr); 88 PetscFunctionReturn(0); 89 } 90 ierr = PetscTypeCompare((PetscObject)A,MATSEQAIJ,&same);CHKERRQ(ierr); 91 if (same) { 92 ierr = MatHYPRE_IJMatrixPreallocate(A,PETSC_NULL,*ij);CHKERRQ(ierr); 93 PetscFunctionReturn(0); 94 } 95 ierr = PetscTypeCompare((PetscObject)A,MATSEQBAIJ,&same);CHKERRQ(ierr); 96 if (same) { 97 ierr = MatHYPRE_IJMatrixPreallocate(A,PETSC_NULL,*ij);CHKERRQ(ierr); 98 PetscFunctionReturn(0); 99 } 100 } 101 PetscFunctionReturn(0); 102 } 103 104 extern PetscErrorCode MatHYPRE_IJMatrixFastCopy_MPIAIJ(Mat,HYPRE_IJMatrix); 105 extern PetscErrorCode MatHYPRE_IJMatrixFastCopy_SeqAIJ(Mat,HYPRE_IJMatrix); 106 /* 107 Copies the data over (column indices, numerical values) to hypre matrix 108 */ 109 110 #undef __FUNCT__ 111 #define __FUNCT__ "MatHYPRE_IJMatrixCopy" 112 PetscErrorCode MatHYPRE_IJMatrixCopy(Mat A,HYPRE_IJMatrix ij) 113 { 114 PetscErrorCode ierr; 115 PetscInt i,rstart,rend,ncols; 116 const PetscScalar *values; 117 const PetscInt *cols; 118 PetscBool flg; 119 120 PetscFunctionBegin; 121 ierr = PetscTypeCompare((PetscObject)A,MATMPIAIJ,&flg);CHKERRQ(ierr); 122 if (flg) { 123 ierr = MatHYPRE_IJMatrixFastCopy_MPIAIJ(A,ij);CHKERRQ(ierr); 124 PetscFunctionReturn(0); 125 } 126 ierr = PetscTypeCompare((PetscObject)A,MATSEQAIJ,&flg);CHKERRQ(ierr); 127 if (flg) { 128 ierr = MatHYPRE_IJMatrixFastCopy_SeqAIJ(A,ij);CHKERRQ(ierr); 129 PetscFunctionReturn(0); 130 } 131 132 ierr = PetscLogEventBegin(MAT_Convert,A,0,0,0);CHKERRQ(ierr); 133 PetscStackCallHypre(0,HYPRE_IJMatrixInitialize,(ij)); 134 ierr = MatGetOwnershipRange(A,&rstart,&rend);CHKERRQ(ierr); 135 for (i=rstart; i<rend; i++) { 136 ierr = MatGetRow(A,i,&ncols,&cols,&values);CHKERRQ(ierr); 137 PetscStackCallHypre(0,HYPRE_IJMatrixSetValues,(ij,1,&ncols,&i,cols,values)); 138 ierr = MatRestoreRow(A,i,&ncols,&cols,&values);CHKERRQ(ierr); 139 } 140 PetscStackCallHypre(0,HYPRE_IJMatrixAssemble,(ij)); 141 ierr = PetscLogEventEnd(MAT_Convert,A,0,0,0);CHKERRQ(ierr); 142 PetscFunctionReturn(0); 143 } 144 145 /* 146 This copies the CSR format directly from the PETSc data structure to the hypre 147 data structure without calls to MatGetRow() or hypre's set values. 148 149 */ 150 #include "_hypre_IJ_mv.h" 151 #include "HYPRE_IJ_mv.h" 152 #include "../src/mat/impls/aij/mpi/mpiaij.h" 153 154 #undef __FUNCT__ 155 #define __FUNCT__ "MatHYPRE_IJMatrixFastCopy_SeqAIJ" 156 PetscErrorCode MatHYPRE_IJMatrixFastCopy_SeqAIJ(Mat A,HYPRE_IJMatrix ij) 157 { 158 PetscErrorCode ierr; 159 Mat_SeqAIJ *pdiag = (Mat_SeqAIJ*)A->data;; 160 161 hypre_ParCSRMatrix *par_matrix; 162 hypre_AuxParCSRMatrix *aux_matrix; 163 hypre_CSRMatrix *hdiag,*hoffd; 164 165 PetscFunctionBegin; 166 PetscValidHeaderSpecific(A,MAT_CLASSID,1); 167 PetscValidType(A,1); 168 PetscValidPointer(ij,2); 169 170 ierr = PetscLogEventBegin(MAT_Convert,A,0,0,0);CHKERRQ(ierr); 171 PetscStackCallHypre(0,HYPRE_IJMatrixInitialize,(ij)); 172 par_matrix = (hypre_ParCSRMatrix*)hypre_IJMatrixObject(ij); 173 aux_matrix = (hypre_AuxParCSRMatrix*)hypre_IJMatrixTranslator(ij); 174 hdiag = hypre_ParCSRMatrixDiag(par_matrix); 175 hoffd = hypre_ParCSRMatrixOffd(par_matrix); 176 177 /* 178 this is the Hack part where we monkey directly with the hypre datastructures 179 */ 180 ierr = PetscMemcpy(hdiag->i,pdiag->i,(A->rmap->n + 1)*sizeof(PetscInt)); 181 ierr = PetscMemcpy(hdiag->j,pdiag->j,pdiag->nz*sizeof(PetscInt)); 182 ierr = PetscMemcpy(hdiag->data,pdiag->a,pdiag->nz*sizeof(PetscScalar)); 183 184 hypre_AuxParCSRMatrixNeedAux(aux_matrix) = 0; 185 PetscStackCallHypre(0,HYPRE_IJMatrixAssemble,(ij)); 186 ierr = PetscLogEventEnd(MAT_Convert,A,0,0,0);CHKERRQ(ierr); 187 PetscFunctionReturn(0); 188 } 189 190 #undef __FUNCT__ 191 #define __FUNCT__ "MatHYPRE_IJMatrixFastCopy_MPIAIJ" 192 PetscErrorCode MatHYPRE_IJMatrixFastCopy_MPIAIJ(Mat A,HYPRE_IJMatrix ij) 193 { 194 PetscErrorCode ierr; 195 Mat_MPIAIJ *pA = (Mat_MPIAIJ*)A->data; 196 Mat_SeqAIJ *pdiag,*poffd; 197 PetscInt i,*garray = pA->garray,*jj,cstart,*pjj; 198 199 hypre_ParCSRMatrix *par_matrix; 200 hypre_AuxParCSRMatrix *aux_matrix; 201 hypre_CSRMatrix *hdiag,*hoffd; 202 203 PetscFunctionBegin; 204 PetscValidHeaderSpecific(A,MAT_CLASSID,1); 205 PetscValidType(A,1); 206 PetscValidPointer(ij,2); 207 pdiag = (Mat_SeqAIJ*) pA->A->data; 208 poffd = (Mat_SeqAIJ*) pA->B->data; 209 /* cstart is only valid for square MPIAIJ layed out in the usual way */ 210 ierr = MatGetOwnershipRange(A,&cstart,PETSC_NULL);CHKERRQ(ierr); 211 212 ierr = PetscLogEventBegin(MAT_Convert,A,0,0,0);CHKERRQ(ierr); 213 PetscStackCallHypre(0,HYPRE_IJMatrixInitialize,(ij)); 214 par_matrix = (hypre_ParCSRMatrix*)hypre_IJMatrixObject(ij); 215 aux_matrix = (hypre_AuxParCSRMatrix*)hypre_IJMatrixTranslator(ij); 216 hdiag = hypre_ParCSRMatrixDiag(par_matrix); 217 hoffd = hypre_ParCSRMatrixOffd(par_matrix); 218 219 /* 220 this is the Hack part where we monkey directly with the hypre datastructures 221 */ 222 ierr = PetscMemcpy(hdiag->i,pdiag->i,(pA->A->rmap->n + 1)*sizeof(PetscInt)); 223 /* need to shift the diag column indices (hdiag->j) back to global numbering since hypre is expecting this */ 224 jj = hdiag->j; 225 pjj = pdiag->j; 226 for (i=0; i<pdiag->nz; i++) { 227 jj[i] = cstart + pjj[i]; 228 } 229 ierr = PetscMemcpy(hdiag->data,pdiag->a,pdiag->nz*sizeof(PetscScalar)); 230 ierr = PetscMemcpy(hoffd->i,poffd->i,(pA->A->rmap->n + 1)*sizeof(PetscInt)); 231 /* need to move the offd column indices (hoffd->j) back to global numbering since hypre is expecting this 232 If we hacked a hypre a bit more we might be able to avoid this step */ 233 jj = hoffd->j; 234 pjj = poffd->j; 235 for (i=0; i<poffd->nz; i++) { 236 jj[i] = garray[pjj[i]]; 237 } 238 ierr = PetscMemcpy(hoffd->data,poffd->a,poffd->nz*sizeof(PetscScalar)); 239 240 hypre_AuxParCSRMatrixNeedAux(aux_matrix) = 0; 241 PetscStackCallHypre(0,HYPRE_IJMatrixAssemble,(ij)); 242 ierr = PetscLogEventEnd(MAT_Convert,A,0,0,0);CHKERRQ(ierr); 243 PetscFunctionReturn(0); 244 } 245 246 /* 247 Does NOT copy the data over, instead uses DIRECTLY the pointers from the PETSc MPIAIJ format 248 249 This is UNFINISHED and does NOT work! The problem is that hypre puts the diagonal entry first 250 which will corrupt the PETSc data structure if we did this. Need a work around to this problem. 251 */ 252 #include "_hypre_IJ_mv.h" 253 #include "HYPRE_IJ_mv.h" 254 255 #undef __FUNCT__ 256 #define __FUNCT__ "MatHYPRE_IJMatrixLink" 257 PetscErrorCode MatHYPRE_IJMatrixLink(Mat A,HYPRE_IJMatrix *ij) 258 { 259 PetscErrorCode ierr; 260 int rstart,rend,cstart,cend; 261 PetscBool flg; 262 hypre_ParCSRMatrix *par_matrix; 263 hypre_AuxParCSRMatrix *aux_matrix; 264 265 PetscFunctionBegin; 266 PetscValidHeaderSpecific(A,MAT_CLASSID,1); 267 PetscValidType(A,1); 268 PetscValidPointer(ij,2); 269 ierr = PetscTypeCompare((PetscObject)A,MATMPIAIJ,&flg);CHKERRQ(ierr); 270 if (!flg) SETERRQ(((PetscObject)A)->comm,PETSC_ERR_SUP,"Can only use with PETSc MPIAIJ matrices"); 271 ierr = MatPreallocated(A);CHKERRQ(ierr); 272 273 ierr = PetscLogEventBegin(MAT_Convert,A,0,0,0);CHKERRQ(ierr); 274 rstart = A->rmap->rstart; 275 rend = A->rmap->rend; 276 cstart = A->cmap->rstart; 277 cend = A->cmap->rend; 278 PetscStackCallHypre(0,HYPRE_IJMatrixCreate,(((PetscObject)A)->comm,rstart,rend-1,cstart,cend-1,ij)); 279 PetscStackCallHypre(0,HYPRE_IJMatrixSetObjectType,(*ij,HYPRE_PARCSR)); 280 281 PetscStackCallHypre(0,HYPRE_IJMatrixInitialize,(*ij)); 282 par_matrix = (hypre_ParCSRMatrix*)hypre_IJMatrixObject(*ij); 283 aux_matrix = (hypre_AuxParCSRMatrix*)hypre_IJMatrixTranslator(*ij); 284 285 hypre_AuxParCSRMatrixNeedAux(aux_matrix) = 0; 286 287 /* this is the Hack part where we monkey directly with the hypre datastructures */ 288 289 PetscStackCallHypre(0,HYPRE_IJMatrixAssemble,(*ij)); 290 ierr = PetscLogEventEnd(MAT_Convert,A,0,0,0);CHKERRQ(ierr); 291 PetscFunctionReturn(0); 292 } 293 294 /* -----------------------------------------------------------------------------------------------------------------*/ 295 296 /*MC 297 MATHYPRESTRUCT - MATHYPRESTRUCT = "hyprestruct" - A matrix type to be used for parallel sparse matrices 298 based on the hypre HYPRE_StructMatrix. 299 300 Level: intermediate 301 302 Notes: Unlike the more general support for blocks in hypre this allows only one block per process and requires the block 303 be defined by a DMDA. 304 305 The matrix needs a DMDA associated with it by either a call to MatSetDM() or if the matrix is obtained from DMGetMatrix() 306 307 .seealso: MatCreate(), PCPFMG, MatSetDA(), DMGetMatrix() 308 M*/ 309 310 311 #undef __FUNCT__ 312 #define __FUNCT__ "MatSetValuesLocal_HYPREStruct_3d" 313 PetscErrorCode MatSetValuesLocal_HYPREStruct_3d(Mat mat,PetscInt nrow,const PetscInt irow[],PetscInt ncol,const PetscInt icol[],const PetscScalar y[],InsertMode addv) 314 { 315 PetscErrorCode ierr; 316 PetscInt i,j,stencil,index[3],row,entries[7]; 317 const PetscScalar *values = y; 318 Mat_HYPREStruct *ex = (Mat_HYPREStruct*) mat->data; 319 320 PetscFunctionBegin; 321 for (i=0; i<nrow; i++) { 322 for (j=0; j<ncol; j++) { 323 stencil = icol[j] - irow[i]; 324 if (!stencil) { 325 entries[j] = 3; 326 } else if (stencil == -1) { 327 entries[j] = 2; 328 } else if (stencil == 1) { 329 entries[j] = 4; 330 } else if (stencil == -ex->gnx) { 331 entries[j] = 1; 332 } else if (stencil == ex->gnx) { 333 entries[j] = 5; 334 } else if (stencil == -ex->gnxgny) { 335 entries[j] = 0; 336 } else if (stencil == ex->gnxgny) { 337 entries[j] = 6; 338 } else SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Local row %D local column %D have bad stencil %D",irow[i],icol[j],stencil); 339 } 340 row = ex->gindices[irow[i]] - ex->rstart; 341 index[0] = ex->xs + (row % ex->nx); 342 index[1] = ex->ys + ((row/ex->nx) % ex->ny); 343 index[2] = ex->zs + (row/(ex->nxny)); 344 if (addv == ADD_VALUES) { 345 PetscStackCallHypre(0,HYPRE_StructMatrixAddToValues,(ex->hmat,index,ncol,entries,(PetscScalar*)values)); 346 } else { 347 PetscStackCallHypre(0,HYPRE_StructMatrixSetValues,(ex->hmat,index,ncol,entries,(PetscScalar*)values)); 348 } 349 values += ncol; 350 } 351 PetscFunctionReturn(0); 352 } 353 354 #undef __FUNCT__ 355 #define __FUNCT__ "MatZeroRowsLocal_HYPREStruct_3d" 356 PetscErrorCode MatZeroRowsLocal_HYPREStruct_3d(Mat mat,PetscInt nrow,const PetscInt irow[],PetscScalar d,Vec x,Vec b) 357 { 358 PetscErrorCode ierr; 359 PetscInt i,index[3],row,entries[7] = {0,1,2,3,4,5,6}; 360 PetscScalar values[7]; 361 Mat_HYPREStruct *ex = (Mat_HYPREStruct*) mat->data; 362 363 PetscFunctionBegin; 364 if (x && b) SETERRQ(((PetscObject)mat)->comm,PETSC_ERR_SUP,"No support"); 365 ierr = PetscMemzero(values,7*sizeof(PetscScalar));CHKERRQ(ierr); 366 values[3] = d; 367 for (i=0; i<nrow; i++) { 368 row = ex->gindices[irow[i]] - ex->rstart; 369 index[0] = ex->xs + (row % ex->nx); 370 index[1] = ex->ys + ((row/ex->nx) % ex->ny); 371 index[2] = ex->zs + (row/(ex->nxny)); 372 PetscStackCallHypre(0,HYPRE_StructMatrixSetValues,(ex->hmat,index,7,entries,values)); 373 } 374 PetscStackCallHypre(0,HYPRE_StructMatrixAssemble,(ex->hmat)); 375 PetscFunctionReturn(0); 376 } 377 378 #undef __FUNCT__ 379 #define __FUNCT__ "MatZeroEntries_HYPREStruct_3d" 380 PetscErrorCode MatZeroEntries_HYPREStruct_3d(Mat mat) 381 { 382 PetscErrorCode ierr; 383 PetscInt indices[7] = {0,1,2,3,4,5,6}; 384 Mat_HYPREStruct *ex = (Mat_HYPREStruct*) mat->data; 385 386 PetscFunctionBegin; 387 /* hypre has no public interface to do this */ 388 PetscStackCallHypre(0,hypre_StructMatrixClearBoxValues,(ex->hmat,&ex->hbox,7,indices,0,1)); 389 PetscStackCallHypre(0,HYPRE_StructMatrixAssemble,(ex->hmat)); 390 PetscFunctionReturn(0); 391 } 392 393 #undef __FUNCT__ 394 #define __FUNCT__ "MatSetDA_HYPREStruct" 395 PetscErrorCode MatSetDA_HYPREStruct(Mat mat,DM da) 396 { 397 PetscErrorCode ierr; 398 Mat_HYPREStruct *ex = (Mat_HYPREStruct*) mat->data; 399 PetscInt dim,dof,sw[3],nx,ny,nz; 400 int ilower[3],iupper[3],ssize,i; 401 DMDABoundaryType p; 402 DMDAStencilType st; 403 404 PetscFunctionBegin; 405 ex->da = da; 406 ierr = PetscObjectReference((PetscObject)da);CHKERRQ(ierr); 407 408 ierr = DMDAGetInfo(ex->da,&dim,0,0,0,0,0,0,&dof,&sw[0],&p,&st);CHKERRQ(ierr); 409 ierr = DMDAGetCorners(ex->da,&ilower[0],&ilower[1],&ilower[2],&iupper[0],&iupper[1],&iupper[2]);CHKERRQ(ierr); 410 iupper[0] += ilower[0] - 1; 411 iupper[1] += ilower[1] - 1; 412 iupper[2] += ilower[2] - 1; 413 414 /* the hypre_Box is used to zero out the matrix entries in MatZeroValues() */ 415 ex->hbox.imin[0] = ilower[0]; 416 ex->hbox.imin[1] = ilower[1]; 417 ex->hbox.imin[2] = ilower[2]; 418 ex->hbox.imax[0] = iupper[0]; 419 ex->hbox.imax[1] = iupper[1]; 420 ex->hbox.imax[2] = iupper[2]; 421 422 /* create the hypre grid object and set its information */ 423 if (dof > 1) SETERRQ(((PetscObject)da)->comm,PETSC_ERR_SUP,"Currently only support for scalar problems"); 424 if (p) SETERRQ(((PetscObject)da)->comm,PETSC_ERR_SUP,"Ask us to add periodic support by calling HYPRE_StructGridSetPeriodic()"); 425 PetscStackCallHypre(0,HYPRE_StructGridCreate,(ex->hcomm,dim,&ex->hgrid)); 426 PetscStackCallHypre(0,HYPRE_StructGridSetExtents,(ex->hgrid,ilower,iupper)); 427 PetscStackCallHypre(0,HYPRE_StructGridAssemble,(ex->hgrid)); 428 429 sw[1] = sw[0]; 430 sw[2] = sw[1]; 431 PetscStackCallHypre(0,HYPRE_StructGridSetNumGhost,(ex->hgrid,sw)); 432 433 /* create the hypre stencil object and set its information */ 434 if (sw[0] > 1) SETERRQ(((PetscObject)da)->comm,PETSC_ERR_SUP,"Ask us to add support for wider stencils"); 435 if (st == DMDA_STENCIL_BOX) SETERRQ(((PetscObject)da)->comm,PETSC_ERR_SUP,"Ask us to add support for box stencils"); 436 if (dim == 1) { 437 int offsets[3][1] = {{-1},{0},{1}}; 438 ssize = 3; 439 PetscStackCallHypre(0,HYPRE_StructStencilCreate,(dim,ssize,&ex->hstencil)); 440 for (i=0; i<ssize; i++) { 441 PetscStackCallHypre(0,HYPRE_StructStencilSetElement,(ex->hstencil,i,offsets[i])); 442 } 443 } else if (dim == 2) { 444 int offsets[5][2] = {{0,-1},{-1,0},{0,0},{1,0},{0,1}}; 445 ssize = 5; 446 PetscStackCallHypre(0,HYPRE_StructStencilCreate,(dim,ssize,&ex->hstencil)); 447 for (i=0; i<ssize; i++) { 448 PetscStackCallHypre(0,HYPRE_StructStencilSetElement,(ex->hstencil,i,offsets[i])); 449 } 450 } else if (dim == 3) { 451 int offsets[7][3] = {{0,0,-1},{0,-1,0},{-1,0,0},{0,0,0},{1,0,0},{0,1,0},{0,0,1}}; 452 ssize = 7; 453 PetscStackCallHypre(0,HYPRE_StructStencilCreate,(dim,ssize,&ex->hstencil)); 454 for (i=0; i<ssize; i++) { 455 PetscStackCallHypre(0,HYPRE_StructStencilSetElement,(ex->hstencil,i,offsets[i])); 456 } 457 } 458 459 /* create the HYPRE vector for rhs and solution */ 460 PetscStackCallHypre(0,HYPRE_StructVectorCreate,(ex->hcomm,ex->hgrid,&ex->hb)); 461 PetscStackCallHypre(0,HYPRE_StructVectorCreate,(ex->hcomm,ex->hgrid,&ex->hx)); 462 PetscStackCallHypre(0,HYPRE_StructVectorInitialize,(ex->hb)); 463 PetscStackCallHypre(0,HYPRE_StructVectorInitialize,(ex->hx)); 464 PetscStackCallHypre(0,HYPRE_StructVectorAssemble,(ex->hb)); 465 PetscStackCallHypre(0,HYPRE_StructVectorAssemble,(ex->hx)); 466 467 /* create the hypre matrix object and set its information */ 468 PetscStackCallHypre(0,HYPRE_StructMatrixCreate,(ex->hcomm,ex->hgrid,ex->hstencil,&ex->hmat)); 469 PetscStackCallHypre(0,HYPRE_StructGridDestroy,(ex->hgrid)); 470 PetscStackCallHypre(0,HYPRE_StructStencilDestroy,(ex->hstencil)); 471 if (ex->needsinitialization) { 472 PetscStackCallHypre(0,HYPRE_StructMatrixInitialize,(ex->hmat)); 473 ex->needsinitialization = PETSC_FALSE; 474 } 475 476 /* set the global and local sizes of the matrix */ 477 ierr = DMDAGetCorners(da,0,0,0,&nx,&ny,&nz);CHKERRQ(ierr); 478 ierr = MatSetSizes(mat,dof*nx*ny*nz,dof*nx*ny*nz,PETSC_DECIDE,PETSC_DECIDE);CHKERRQ(ierr); 479 ierr = PetscLayoutSetBlockSize(mat->rmap,1);CHKERRQ(ierr); 480 ierr = PetscLayoutSetBlockSize(mat->cmap,1);CHKERRQ(ierr); 481 ierr = PetscLayoutSetUp(mat->rmap);CHKERRQ(ierr); 482 ierr = PetscLayoutSetUp(mat->cmap);CHKERRQ(ierr); 483 484 if (dim == 3) { 485 mat->ops->setvalueslocal = MatSetValuesLocal_HYPREStruct_3d; 486 mat->ops->zerorowslocal = MatZeroRowsLocal_HYPREStruct_3d; 487 mat->ops->zeroentries = MatZeroEntries_HYPREStruct_3d; 488 ierr = MatZeroEntries_HYPREStruct_3d(mat);CHKERRQ(ierr); 489 } else SETERRQ(((PetscObject)da)->comm,PETSC_ERR_SUP,"Only support for 3d DMDA currently"); 490 491 /* get values that will be used repeatedly in MatSetValuesLocal() and MatZeroRowsLocal() repeatedly */ 492 ierr = MatGetOwnershipRange(mat,&ex->rstart,PETSC_NULL);CHKERRQ(ierr); 493 ierr = DMDAGetGlobalIndices(ex->da,PETSC_NULL,&ex->gindices);CHKERRQ(ierr); 494 ierr = DMDAGetGhostCorners(ex->da,0,0,0,&ex->gnx,&ex->gnxgny,0);CHKERRQ(ierr); 495 ex->gnxgny *= ex->gnx; 496 ierr = DMDAGetCorners(ex->da,&ex->xs,&ex->ys,&ex->zs,&ex->nx,&ex->ny,0);CHKERRQ(ierr); 497 ex->nxny = ex->nx*ex->ny; 498 PetscFunctionReturn(0); 499 } 500 501 #undef __FUNCT__ 502 #define __FUNCT__ "MatMult_HYPREStruct" 503 PetscErrorCode MatMult_HYPREStruct(Mat A,Vec x,Vec y) 504 { 505 PetscErrorCode ierr; 506 PetscScalar *xx,*yy; 507 int ilower[3],iupper[3]; 508 Mat_HYPREStruct *mx = (Mat_HYPREStruct *)(A->data); 509 510 PetscFunctionBegin; 511 ierr = DMDAGetCorners(mx->da,&ilower[0],&ilower[1],&ilower[2],&iupper[0],&iupper[1],&iupper[2]);CHKERRQ(ierr); 512 iupper[0] += ilower[0] - 1; 513 iupper[1] += ilower[1] - 1; 514 iupper[2] += ilower[2] - 1; 515 516 /* copy x values over to hypre */ 517 PetscStackCallHypre(0,HYPRE_StructVectorSetConstantValues,(mx->hb,0.0)); 518 ierr = VecGetArray(x,&xx);CHKERRQ(ierr); 519 PetscStackCallHypre(0,HYPRE_StructVectorSetBoxValues,(mx->hb,ilower,iupper,xx)); 520 ierr = VecRestoreArray(x,&xx);CHKERRQ(ierr); 521 PetscStackCallHypre(0,HYPRE_StructVectorAssemble,(mx->hb)); 522 PetscStackCallHypre(0,HYPRE_StructMatrixMatvec,(1.0,mx->hmat,mx->hb,0.0,mx->hx)); 523 524 /* copy solution values back to PETSc */ 525 ierr = VecGetArray(y,&yy);CHKERRQ(ierr); 526 PetscStackCallHypre(0,HYPRE_StructVectorGetBoxValues,(mx->hx,ilower,iupper,yy)); 527 ierr = VecRestoreArray(y,&yy);CHKERRQ(ierr); 528 PetscFunctionReturn(0); 529 } 530 531 #undef __FUNCT__ 532 #define __FUNCT__ "MatAssemblyEnd_HYPREStruct" 533 PetscErrorCode MatAssemblyEnd_HYPREStruct(Mat mat,MatAssemblyType mode) 534 { 535 Mat_HYPREStruct *ex = (Mat_HYPREStruct*) mat->data; 536 PetscErrorCode ierr; 537 538 PetscFunctionBegin; 539 PetscStackCallHypre(0,HYPRE_StructMatrixAssemble,(ex->hmat)); 540 /* PetscStackCallHypre(0,HYPRE_StructMatrixPrint,("dummy",ex->hmat,0)); */ 541 PetscFunctionReturn(0); 542 } 543 544 #undef __FUNCT__ 545 #define __FUNCT__ "MatZeroEntries_HYPREStruct" 546 PetscErrorCode MatZeroEntries_HYPREStruct(Mat mat) 547 { 548 PetscFunctionBegin; 549 /* before the DMDA is set to the matrix the zero doesn't need to do anything */ 550 PetscFunctionReturn(0); 551 } 552 553 554 #undef __FUNCT__ 555 #define __FUNCT__ "MatDestroy_HYPREStruct" 556 PetscErrorCode MatDestroy_HYPREStruct(Mat mat) 557 { 558 Mat_HYPREStruct *ex = (Mat_HYPREStruct*) mat->data; 559 PetscErrorCode ierr; 560 561 PetscFunctionBegin; 562 PetscStackCallHypre(0,HYPRE_StructMatrixDestroy,(ex->hmat)); 563 PetscStackCallHypre(0,HYPRE_StructVectorDestroy,(ex->hx)); 564 PetscStackCallHypre(0,HYPRE_StructVectorDestroy,(ex->hb)); 565 PetscFunctionReturn(0); 566 } 567 568 569 EXTERN_C_BEGIN 570 #undef __FUNCT__ 571 #define __FUNCT__ "MatCreate_HYPREStruct" 572 PetscErrorCode MatCreate_HYPREStruct(Mat B) 573 { 574 Mat_HYPREStruct *ex; 575 PetscErrorCode ierr; 576 577 PetscFunctionBegin; 578 ierr = PetscNewLog(B,Mat_HYPREStruct,&ex);CHKERRQ(ierr); 579 B->data = (void*)ex; 580 B->rmap->bs = 1; 581 B->assembled = PETSC_FALSE; 582 583 B->insertmode = NOT_SET_VALUES; 584 585 B->ops->assemblyend = MatAssemblyEnd_HYPREStruct; 586 B->ops->mult = MatMult_HYPREStruct; 587 B->ops->zeroentries = MatZeroEntries_HYPREStruct; 588 B->ops->destroy = MatDestroy_HYPREStruct; 589 590 ex->needsinitialization = PETSC_TRUE; 591 592 ierr = MPI_Comm_dup(((PetscObject)B)->comm,&(ex->hcomm));CHKERRQ(ierr); 593 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatSetDA_C","MatSetDA_HYPREStruct",MatSetDA_HYPREStruct);CHKERRQ(ierr); 594 ierr = PetscObjectChangeTypeName((PetscObject)B,MATHYPRESTRUCT);CHKERRQ(ierr); 595 PetscFunctionReturn(0); 596 } 597 EXTERN_C_END 598 599 600 /*MC 601 MATHYPRESSTRUCT - MATHYPRESSTRUCT = "hypresstruct" - A matrix type to be used for parallel sparse matrices 602 based on the hypre HYPRE_SStructMatrix. 603 604 605 Level: intermediate 606 607 Notes: Unlike hypre's general semi-struct object consisting of a collection of structured-grid objects and unstructured 608 grid objects, we will restrict the semi-struct objects to consist of only structured-grid components. 609 610 Unlike the more general support for parts and blocks in hypre this allows only one part, and one block per process and requires the block 611 be defined by a DMDA. 612 613 The matrix needs a DMDA associated with it by either a call to MatSetDA() or if the matrix is obtained from DMGetMatrix() 614 615 M*/ 616 617 #undef __FUNCT__ 618 #define __FUNCT__ "MatSetValuesLocal_HYPRESStruct_3d" 619 PetscErrorCode MatSetValuesLocal_HYPRESStruct_3d(Mat mat,PetscInt nrow,const PetscInt irow[],PetscInt ncol,const PetscInt icol[],const PetscScalar y[],InsertMode addv) 620 { 621 PetscErrorCode ierr; 622 PetscInt i,j,stencil,index[3]; 623 const PetscScalar *values = y; 624 Mat_HYPRESStruct *ex = (Mat_HYPRESStruct*) mat->data; 625 626 int part= 0; /* Petsc sstruct interface only allows 1 part */ 627 int ordering; 628 int grid_rank, to_grid_rank; 629 int var_type, to_var_type; 630 int to_var_entry = 0; 631 632 int nvars= ex->nvars; 633 PetscInt row,*entries; 634 635 PetscFunctionBegin; 636 ierr = PetscMalloc(7*nvars*sizeof(PetscInt),&entries);CHKERRQ(ierr); 637 ordering= ex-> dofs_order; /* ordering= 0 nodal ordering 638 1 variable ordering */ 639 /* stencil entries are orderer by variables: var0_stencil0, var0_stencil1, ..., var0_stencil6, var1_stencil0, var1_stencil1, ... */ 640 641 if (!ordering) { 642 for (i=0; i<nrow; i++) { 643 grid_rank= irow[i]/nvars; 644 var_type = (irow[i] % nvars); 645 646 for (j=0; j<ncol; j++) { 647 to_grid_rank= icol[j]/nvars; 648 to_var_type = (icol[j] % nvars); 649 650 to_var_entry= to_var_entry*7; 651 entries[j]= to_var_entry; 652 653 stencil = to_grid_rank-grid_rank; 654 if (!stencil) { 655 entries[j] += 3; 656 } else if (stencil == -1) { 657 entries[j] += 2; 658 } else if (stencil == 1) { 659 entries[j] += 4; 660 } else if (stencil == -ex->gnx) { 661 entries[j] += 1; 662 } else if (stencil == ex->gnx) { 663 entries[j] += 5; 664 } else if (stencil == -ex->gnxgny) { 665 entries[j] += 0; 666 } else if (stencil == ex->gnxgny) { 667 entries[j] += 6; 668 } else SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Local row %D local column %D have bad stencil %D",irow[i],icol[j],stencil); 669 } 670 671 row = ex->gindices[grid_rank] - ex->rstart; 672 index[0] = ex->xs + (row % ex->nx); 673 index[1] = ex->ys + ((row/ex->nx) % ex->ny); 674 index[2] = ex->zs + (row/(ex->nxny)); 675 676 if (addv == ADD_VALUES) { 677 PetscStackCallHypre(0,HYPRE_SStructMatrixAddToValues,(ex->ss_mat,part,index,var_type,ncol,entries,(PetscScalar*)values)); 678 } else { 679 PetscStackCallHypre(0,HYPRE_SStructMatrixSetValues,(ex->ss_mat,part,index,var_type,ncol,entries,(PetscScalar*)values)); 680 } 681 values += ncol; 682 } 683 } else { 684 for (i=0; i<nrow; i++) { 685 var_type = irow[i]/(ex->gnxgnygnz); 686 grid_rank= irow[i] - var_type*(ex->gnxgnygnz); 687 688 for (j=0; j<ncol; j++) { 689 to_var_type = icol[j]/(ex->gnxgnygnz); 690 to_grid_rank= icol[j] - to_var_type*(ex->gnxgnygnz); 691 692 to_var_entry= to_var_entry*7; 693 entries[j]= to_var_entry; 694 695 stencil = to_grid_rank-grid_rank; 696 if (!stencil) { 697 entries[j] += 3; 698 } else if (stencil == -1) { 699 entries[j] += 2; 700 } else if (stencil == 1) { 701 entries[j] += 4; 702 } else if (stencil == -ex->gnx) { 703 entries[j] += 1; 704 } else if (stencil == ex->gnx) { 705 entries[j] += 5; 706 } else if (stencil == -ex->gnxgny) { 707 entries[j] += 0; 708 } else if (stencil == ex->gnxgny) { 709 entries[j] += 6; 710 } else SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Local row %D local column %D have bad stencil %D",irow[i],icol[j],stencil); 711 } 712 713 row = ex->gindices[grid_rank] - ex->rstart; 714 index[0] = ex->xs + (row % ex->nx); 715 index[1] = ex->ys + ((row/ex->nx) % ex->ny); 716 index[2] = ex->zs + (row/(ex->nxny)); 717 718 if (addv == ADD_VALUES) { 719 PetscStackCallHypre(0,HYPRE_SStructMatrixAddToValues,(ex->ss_mat,part,index,var_type,ncol,entries,(PetscScalar*)values)); 720 } else { 721 PetscStackCallHypre(0,HYPRE_SStructMatrixSetValues,(ex->ss_mat,part,index,var_type,ncol,entries,(PetscScalar*)values)); 722 } 723 values += ncol; 724 } 725 } 726 ierr = PetscFree(entries);CHKERRQ(ierr); 727 PetscFunctionReturn(0); 728 } 729 730 #undef __FUNCT__ 731 #define __FUNCT__ "MatZeroRowsLocal_HYPRESStruct_3d" 732 PetscErrorCode MatZeroRowsLocal_HYPRESStruct_3d(Mat mat,PetscInt nrow,const PetscInt irow[],PetscScalar d,Vec x,Vec b) 733 { 734 PetscErrorCode ierr; 735 PetscInt i,index[3]; 736 PetscScalar **values; 737 Mat_HYPRESStruct *ex = (Mat_HYPRESStruct*) mat->data; 738 739 int part= 0; /* Petsc sstruct interface only allows 1 part */ 740 int ordering= ex->dofs_order; 741 int grid_rank; 742 int var_type; 743 int nvars= ex->nvars; 744 PetscInt row,*entries; 745 746 PetscFunctionBegin; 747 if (x && b) SETERRQ(((PetscObject)mat)->comm,PETSC_ERR_SUP,"No support"); 748 ierr = PetscMalloc(7*nvars*sizeof(PetscInt),&entries);CHKERRQ(ierr); 749 750 ierr = PetscMalloc(nvars*sizeof(PetscScalar *),&values);CHKERRQ(ierr); 751 ierr = PetscMalloc(7*nvars*nvars*sizeof(PetscScalar),&values[0]);CHKERRQ(ierr); 752 for (i=1; i<nvars; i++) { 753 values[i] = values[i-1] + nvars*7; 754 } 755 756 for (i=0; i< nvars; i++) { 757 ierr = PetscMemzero(values[i],nvars*7*sizeof(PetscScalar));CHKERRQ(ierr); 758 *(values[i]+3)= d; 759 } 760 761 for (i= 0; i< nvars*7; i++) { 762 entries[i]= i; 763 } 764 765 if (!ordering) { 766 for (i=0; i<nrow; i++) { 767 grid_rank= irow[i]/nvars; 768 var_type = (irow[i] % nvars); 769 770 row = ex->gindices[grid_rank] - ex->rstart; 771 index[0] = ex->xs + (row % ex->nx); 772 index[1] = ex->ys + ((row/ex->nx) % ex->ny); 773 index[2] = ex->zs + (row/(ex->nxny)); 774 PetscStackCallHypre(0,HYPRE_SStructMatrixSetValues,(ex->ss_mat,part,index,var_type,7*nvars,entries,values[var_type])); 775 } 776 } else { 777 for (i=0; i<nrow; i++) { 778 var_type = irow[i]/(ex->gnxgnygnz); 779 grid_rank= irow[i] - var_type*(ex->gnxgnygnz); 780 781 row = ex->gindices[grid_rank] - ex->rstart; 782 index[0] = ex->xs + (row % ex->nx); 783 index[1] = ex->ys + ((row/ex->nx) % ex->ny); 784 index[2] = ex->zs + (row/(ex->nxny)); 785 PetscStackCallHypre(0,HYPRE_SStructMatrixSetValues,(ex->ss_mat,part,index,var_type,7*nvars,entries,values[var_type])); 786 } 787 } 788 PetscStackCallHypre(0,HYPRE_SStructMatrixAssemble,(ex->ss_mat)); 789 ierr = PetscFree(values[0]);CHKERRQ(ierr); 790 ierr = PetscFree(values);CHKERRQ(ierr); 791 ierr = PetscFree(entries);CHKERRQ(ierr); 792 PetscFunctionReturn(0); 793 } 794 795 #undef __FUNCT__ 796 #define __FUNCT__ "MatZeroEntries_HYPRESStruct_3d" 797 PetscErrorCode MatZeroEntries_HYPRESStruct_3d(Mat mat) 798 { 799 PetscErrorCode ierr; 800 Mat_HYPRESStruct *ex = (Mat_HYPRESStruct*) mat->data; 801 int nvars= ex->nvars; 802 int size; 803 int part= 0; /* only one part */ 804 805 PetscFunctionBegin; 806 size = ((ex->hbox.imax[0])-(ex->hbox.imin[0])+1)*((ex->hbox.imax[1])-(ex->hbox.imin[1])+1)*((ex->hbox.imax[2])-(ex->hbox.imin[2])+1); 807 { 808 PetscInt i,*entries; 809 PetscScalar *values; 810 int iupper[3], ilower[3]; 811 812 for (i= 0; i< 3; i++) { 813 ilower[i]= ex->hbox.imin[i]; 814 iupper[i]= ex->hbox.imax[i]; 815 } 816 817 ierr = PetscMalloc2(nvars*7,PetscInt,&entries,nvars*7*size,PetscScalar,&values);CHKERRQ(ierr); 818 for (i= 0; i< nvars*7; i++) { 819 entries[i]= i; 820 } 821 ierr = PetscMemzero(values,nvars*7*size*sizeof(PetscScalar));CHKERRQ(ierr); 822 823 for (i= 0; i< nvars; i++) { 824 PetscStackCallHypre(0,HYPRE_SStructMatrixSetBoxValues,(ex->ss_mat,part,ilower,iupper,i,nvars*7,entries,values)); 825 } 826 ierr = PetscFree2(entries,values);CHKERRQ(ierr); 827 } 828 PetscStackCallHypre(0,HYPRE_SStructMatrixAssemble,(ex->ss_mat)); 829 PetscFunctionReturn(0); 830 } 831 832 833 #undef __FUNCT__ 834 #define __FUNCT__ "MatSetDA_HYPRESStruct" 835 PetscErrorCode MatSetDA_HYPRESStruct(Mat mat,DM da) 836 { 837 PetscErrorCode ierr; 838 Mat_HYPRESStruct *ex = (Mat_HYPRESStruct*) mat->data; 839 PetscInt dim,dof,sw[3],nx,ny,nz; 840 int ilower[3],iupper[3],ssize,i; 841 DMDABoundaryType p; 842 DMDAStencilType st; 843 int nparts= 1; /* assuming only one part */ 844 int part = 0; 845 846 PetscFunctionBegin; 847 ex->da = da; 848 ierr = PetscObjectReference((PetscObject)da);CHKERRQ(ierr); 849 850 ierr = DMDAGetInfo(ex->da,&dim,0,0,0,0,0,0,&dof,&sw[0],&p,&st);CHKERRQ(ierr); 851 ierr = DMDAGetCorners(ex->da,&ilower[0],&ilower[1],&ilower[2],&iupper[0],&iupper[1],&iupper[2]);CHKERRQ(ierr); 852 iupper[0] += ilower[0] - 1; 853 iupper[1] += ilower[1] - 1; 854 iupper[2] += ilower[2] - 1; 855 /* the hypre_Box is used to zero out the matrix entries in MatZeroValues() */ 856 ex->hbox.imin[0] = ilower[0]; 857 ex->hbox.imin[1] = ilower[1]; 858 ex->hbox.imin[2] = ilower[2]; 859 ex->hbox.imax[0] = iupper[0]; 860 ex->hbox.imax[1] = iupper[1]; 861 ex->hbox.imax[2] = iupper[2]; 862 863 ex->dofs_order = 0; 864 865 /* assuming that the same number of dofs on each gridpoint. Also assume all cell-centred based */ 866 ex->nvars= dof; 867 868 /* create the hypre grid object and set its information */ 869 if (p) SETERRQ(((PetscObject)da)->comm,PETSC_ERR_SUP,"Ask us to add periodic support by calling HYPRE_SStructGridSetPeriodic()"); 870 PetscStackCallHypre(0,HYPRE_SStructGridCreate,(ex->hcomm,dim,nparts,&ex->ss_grid)); 871 872 PetscStackCallHypre(0,HYPRE_SStructGridSetExtents,(ex->ss_grid,part,ex->hbox.imin,ex->hbox.imax)); 873 874 { 875 HYPRE_SStructVariable *vartypes; 876 ierr = PetscMalloc(ex->nvars*sizeof(HYPRE_SStructVariable),&vartypes);CHKERRQ(ierr); 877 for (i= 0; i< ex->nvars; i++) { 878 vartypes[i]= HYPRE_SSTRUCT_VARIABLE_CELL; 879 } 880 PetscStackCallHypre(0,HYPRE_SStructGridSetVariables,(ex->ss_grid, part, ex->nvars,vartypes)); 881 ierr = PetscFree(vartypes);CHKERRQ(ierr); 882 } 883 PetscStackCallHypre(0,HYPRE_SStructGridAssemble,(ex->ss_grid)); 884 885 sw[1] = sw[0]; 886 sw[2] = sw[1]; 887 /* PetscStackCallHypre(0,HYPRE_SStructGridSetNumGhost,(ex->ss_grid,sw)); */ 888 889 /* create the hypre stencil object and set its information */ 890 if (sw[0] > 1) SETERRQ(((PetscObject)da)->comm,PETSC_ERR_SUP,"Ask us to add support for wider stencils"); 891 if (st == DMDA_STENCIL_BOX) SETERRQ(((PetscObject)da)->comm,PETSC_ERR_SUP,"Ask us to add support for box stencils"); 892 893 if (dim == 1) { 894 int offsets[3][1] = {{-1},{0},{1}}; 895 int j, cnt; 896 897 ssize = 3*(ex->nvars); 898 PetscStackCallHypre(0,HYPRE_SStructStencilCreate,(dim,ssize,&ex->ss_stencil)); 899 cnt= 0; 900 for (i= 0; i< (ex->nvars); i++) { 901 for (j= 0; j< 3; j++) { 902 PetscStackCallHypre(0,HYPRE_SStructStencilSetEntry,(ex->ss_stencil, cnt, offsets[j], i)); 903 cnt++; 904 } 905 } 906 } else if (dim == 2) { 907 int offsets[5][2] = {{0,-1},{-1,0},{0,0},{1,0},{0,1}}; 908 int j, cnt; 909 910 ssize = 5*(ex->nvars); 911 PetscStackCallHypre(0,HYPRE_SStructStencilCreate,(dim,ssize,&ex->ss_stencil)); 912 cnt= 0; 913 for (i= 0; i< (ex->nvars); i++) { 914 for (j= 0; j< 5; j++) { 915 PetscStackCallHypre(0,HYPRE_SStructStencilSetEntry,(ex->ss_stencil, cnt, offsets[j], i)); 916 cnt++; 917 } 918 } 919 } else if (dim == 3) { 920 int offsets[7][3] = {{0,0,-1},{0,-1,0},{-1,0,0},{0,0,0},{1,0,0},{0,1,0},{0,0,1}}; 921 int j, cnt; 922 923 ssize = 7*(ex->nvars); 924 PetscStackCallHypre(0,HYPRE_SStructStencilCreate,(dim,ssize,&ex->ss_stencil)); 925 cnt= 0; 926 for (i= 0; i< (ex->nvars); i++) { 927 for (j= 0; j< 7; j++) { 928 PetscStackCallHypre(0,HYPRE_SStructStencilSetEntry,(ex->ss_stencil, cnt, offsets[j], i)); 929 cnt++; 930 } 931 } 932 } 933 934 /* create the HYPRE graph */ 935 PetscStackCallHypre(0,HYPRE_SStructGraphCreate,(ex->hcomm, ex->ss_grid, &(ex->ss_graph))); 936 937 /* set the stencil graph. Note that each variable has the same graph. This means that each 938 variable couples to all the other variable and with the same stencil pattern. */ 939 for (i= 0; i< (ex->nvars); i++) { 940 PetscStackCallHypre(0,HYPRE_SStructGraphSetStencil,(ex->ss_graph,part,i,ex->ss_stencil)); 941 } 942 PetscStackCallHypre(0,HYPRE_SStructGraphAssemble,(ex->ss_graph)); 943 944 /* create the HYPRE sstruct vectors for rhs and solution */ 945 PetscStackCallHypre(0,HYPRE_SStructVectorCreate,(ex->hcomm,ex->ss_grid,&ex->ss_b)); 946 PetscStackCallHypre(0,HYPRE_SStructVectorCreate,(ex->hcomm,ex->ss_grid,&ex->ss_x)); 947 PetscStackCallHypre(0,HYPRE_SStructVectorInitialize,(ex->ss_b)); 948 PetscStackCallHypre(0,HYPRE_SStructVectorInitialize,(ex->ss_x)); 949 PetscStackCallHypre(0,HYPRE_SStructVectorAssemble,(ex->ss_b)); 950 PetscStackCallHypre(0,HYPRE_SStructVectorAssemble,(ex->ss_x)); 951 952 /* create the hypre matrix object and set its information */ 953 PetscStackCallHypre(0,HYPRE_SStructMatrixCreate,(ex->hcomm,ex->ss_graph,&ex->ss_mat)); 954 PetscStackCallHypre(0,HYPRE_SStructGridDestroy,(ex->ss_grid)); 955 PetscStackCallHypre(0,HYPRE_SStructStencilDestroy,(ex->ss_stencil)); 956 if (ex->needsinitialization) { 957 PetscStackCallHypre(0,HYPRE_SStructMatrixInitialize,(ex->ss_mat)); 958 ex->needsinitialization = PETSC_FALSE; 959 } 960 961 /* set the global and local sizes of the matrix */ 962 ierr = DMDAGetCorners(da,0,0,0,&nx,&ny,&nz);CHKERRQ(ierr); 963 ierr = MatSetSizes(mat,dof*nx*ny*nz,dof*nx*ny*nz,PETSC_DECIDE,PETSC_DECIDE);CHKERRQ(ierr); 964 ierr = PetscLayoutSetBlockSize(mat->rmap,1);CHKERRQ(ierr); 965 ierr = PetscLayoutSetBlockSize(mat->cmap,1);CHKERRQ(ierr); 966 ierr = PetscLayoutSetUp(mat->rmap);CHKERRQ(ierr); 967 ierr = PetscLayoutSetUp(mat->cmap);CHKERRQ(ierr); 968 969 if (dim == 3) { 970 mat->ops->setvalueslocal = MatSetValuesLocal_HYPRESStruct_3d; 971 mat->ops->zerorowslocal = MatZeroRowsLocal_HYPRESStruct_3d; 972 mat->ops->zeroentries = MatZeroEntries_HYPRESStruct_3d; 973 ierr = MatZeroEntries_HYPRESStruct_3d(mat);CHKERRQ(ierr); 974 } else SETERRQ(((PetscObject)da)->comm,PETSC_ERR_SUP,"Only support for 3d DMDA currently"); 975 976 /* get values that will be used repeatedly in MatSetValuesLocal() and MatZeroRowsLocal() repeatedly */ 977 ierr = MatGetOwnershipRange(mat,&ex->rstart,PETSC_NULL);CHKERRQ(ierr); 978 ierr = DMDAGetGlobalIndices(ex->da,PETSC_NULL,&ex->gindices);CHKERRQ(ierr); 979 ierr = DMDAGetGhostCorners(ex->da,0,0,0,&ex->gnx,&ex->gnxgny,&ex->gnxgnygnz);CHKERRQ(ierr); 980 ex->gnxgny *= ex->gnx; 981 ex->gnxgnygnz *= ex->gnxgny; 982 ierr = DMDAGetCorners(ex->da,&ex->xs,&ex->ys,&ex->zs,&ex->nx,&ex->ny,&ex->nz);CHKERRQ(ierr); 983 ex->nxny = ex->nx*ex->ny; 984 ex->nxnynz = ex->nz*ex->nxny; 985 PetscFunctionReturn(0); 986 } 987 988 #undef __FUNCT__ 989 #define __FUNCT__ "MatMult_HYPRESStruct" 990 PetscErrorCode MatMult_HYPRESStruct(Mat A,Vec x,Vec y) 991 { 992 PetscErrorCode ierr; 993 PetscScalar *xx,*yy; 994 int ilower[3],iupper[3]; 995 Mat_HYPRESStruct *mx = (Mat_HYPRESStruct *)(A->data); 996 int ordering= mx->dofs_order; 997 int nvars= mx->nvars; 998 int part= 0; 999 int size; 1000 int i; 1001 1002 PetscFunctionBegin; 1003 ierr = DMDAGetCorners(mx->da,&ilower[0],&ilower[1],&ilower[2],&iupper[0],&iupper[1],&iupper[2]);CHKERRQ(ierr); 1004 iupper[0] += ilower[0] - 1; 1005 iupper[1] += ilower[1] - 1; 1006 iupper[2] += ilower[2] - 1; 1007 1008 size= 1; 1009 for (i= 0; i< 3; i++) { 1010 size*= (iupper[i]-ilower[i]+1); 1011 } 1012 1013 /* copy x values over to hypre for variable ordering */ 1014 if (ordering) { 1015 PetscStackCallHypre(0,HYPRE_SStructVectorSetConstantValues,(mx->ss_b,0.0)); 1016 ierr = VecGetArray(x,&xx);CHKERRQ(ierr); 1017 for (i= 0; i< nvars; i++) { 1018 PetscStackCallHypre(0,HYPRE_SStructVectorSetBoxValues,(mx->ss_b,part,ilower,iupper,i,xx+(size*i))); 1019 } 1020 ierr = VecRestoreArray(x,&xx);CHKERRQ(ierr); 1021 PetscStackCallHypre(0,HYPRE_SStructVectorAssemble,(mx->ss_b)); 1022 PetscStackCallHypre(0,HYPRE_SStructMatrixMatvec,(1.0,mx->ss_mat,mx->ss_b,0.0,mx->ss_x)); 1023 1024 /* copy solution values back to PETSc */ 1025 ierr = VecGetArray(y,&yy);CHKERRQ(ierr); 1026 for (i= 0; i< nvars; i++) { 1027 PetscStackCallHypre(0,HYPRE_SStructVectorGetBoxValues,(mx->ss_x,part,ilower,iupper,i,yy+(size*i))); 1028 } 1029 ierr = VecRestoreArray(y,&yy);CHKERRQ(ierr); 1030 } else { /* nodal ordering must be mapped to variable ordering for sys_pfmg */ 1031 PetscScalar *z; 1032 int j, k; 1033 1034 ierr = PetscMalloc(nvars*size*sizeof(PetscScalar),&z);CHKERRQ(ierr); 1035 PetscStackCallHypre(0,HYPRE_SStructVectorSetConstantValues,(mx->ss_b,0.0)); 1036 ierr = VecGetArray(x,&xx);CHKERRQ(ierr); 1037 1038 /* transform nodal to hypre's variable ordering for sys_pfmg */ 1039 for (i= 0; i< size; i++) { 1040 k= i*nvars; 1041 for (j= 0; j< nvars; j++) { 1042 z[j*size+i]= xx[k+j]; 1043 } 1044 } 1045 for (i= 0; i< nvars; i++) { 1046 PetscStackCallHypre(0,HYPRE_SStructVectorSetBoxValues,(mx->ss_b,part,ilower,iupper,i,z+(size*i))); 1047 } 1048 ierr = VecRestoreArray(x,&xx);CHKERRQ(ierr); 1049 PetscStackCallHypre(0,HYPRE_SStructVectorAssemble,(mx->ss_b)); 1050 PetscStackCallHypre(0,HYPRE_SStructMatrixMatvec,(1.0,mx->ss_mat,mx->ss_b,0.0,mx->ss_x)); 1051 1052 /* copy solution values back to PETSc */ 1053 ierr = VecGetArray(y,&yy);CHKERRQ(ierr); 1054 for (i= 0; i< nvars; i++) { 1055 PetscStackCallHypre(0,HYPRE_SStructVectorGetBoxValues,(mx->ss_x,part,ilower,iupper,i,z+(size*i))); 1056 } 1057 /* transform hypre's variable ordering for sys_pfmg to nodal ordering */ 1058 for (i= 0; i< size; i++) { 1059 k= i*nvars; 1060 for (j= 0; j< nvars; j++) { 1061 yy[k+j]= z[j*size+i]; 1062 } 1063 } 1064 ierr = VecRestoreArray(y,&yy);CHKERRQ(ierr); 1065 ierr = PetscFree(z);CHKERRQ(ierr); 1066 } 1067 PetscFunctionReturn(0); 1068 } 1069 1070 #undef __FUNCT__ 1071 #define __FUNCT__ "MatAssemblyEnd_HYPRESStruct" 1072 PetscErrorCode MatAssemblyEnd_HYPRESStruct(Mat mat,MatAssemblyType mode) 1073 { 1074 Mat_HYPRESStruct *ex = (Mat_HYPRESStruct*) mat->data; 1075 PetscErrorCode ierr; 1076 1077 PetscFunctionBegin; 1078 PetscStackCallHypre(0,HYPRE_SStructMatrixAssemble,(ex->ss_mat)); 1079 PetscFunctionReturn(0); 1080 } 1081 1082 #undef __FUNCT__ 1083 #define __FUNCT__ "MatZeroEntries_HYPRESStruct" 1084 PetscErrorCode MatZeroEntries_HYPRESStruct(Mat mat) 1085 { 1086 PetscFunctionBegin; 1087 /* before the DMDA is set to the matrix the zero doesn't need to do anything */ 1088 PetscFunctionReturn(0); 1089 } 1090 1091 1092 #undef __FUNCT__ 1093 #define __FUNCT__ "MatDestroy_HYPRESStruct" 1094 PetscErrorCode MatDestroy_HYPRESStruct(Mat mat) 1095 { 1096 Mat_HYPRESStruct *ex = (Mat_HYPRESStruct*) mat->data; 1097 PetscErrorCode ierr; 1098 1099 PetscFunctionBegin; 1100 PetscStackCallHypre(0,HYPRE_SStructGraphDestroy,(ex->ss_graph)); 1101 PetscStackCallHypre(0,HYPRE_SStructMatrixDestroy,(ex->ss_mat)); 1102 PetscStackCallHypre(0,HYPRE_SStructVectorDestroy,(ex->ss_x)); 1103 PetscStackCallHypre(0,HYPRE_SStructVectorDestroy,(ex->ss_b)); 1104 PetscFunctionReturn(0); 1105 } 1106 1107 EXTERN_C_BEGIN 1108 #undef __FUNCT__ 1109 #define __FUNCT__ "MatCreate_HYPRESStruct" 1110 PetscErrorCode MatCreate_HYPRESStruct(Mat B) 1111 { 1112 Mat_HYPRESStruct *ex; 1113 PetscErrorCode ierr; 1114 1115 PetscFunctionBegin; 1116 ierr = PetscNewLog(B,Mat_HYPRESStruct,&ex);CHKERRQ(ierr); 1117 B->data = (void*)ex; 1118 B->rmap->bs = 1; 1119 B->assembled = PETSC_FALSE; 1120 1121 B->insertmode = NOT_SET_VALUES; 1122 1123 B->ops->assemblyend = MatAssemblyEnd_HYPRESStruct; 1124 B->ops->mult = MatMult_HYPRESStruct; 1125 B->ops->zeroentries = MatZeroEntries_HYPRESStruct; 1126 B->ops->destroy = MatDestroy_HYPRESStruct; 1127 1128 ex->needsinitialization = PETSC_TRUE; 1129 1130 ierr = MPI_Comm_dup(((PetscObject)B)->comm,&(ex->hcomm));CHKERRQ(ierr); 1131 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatSetDA_C","MatSetDA_HYPRESStruct",MatSetDA_HYPRESStruct);CHKERRQ(ierr); 1132 ierr = PetscObjectChangeTypeName((PetscObject)B,MATHYPRESSTRUCT);CHKERRQ(ierr); 1133 PetscFunctionReturn(0); 1134 } 1135 EXTERN_C_END 1136 1137 1138 1139