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