1 2 /* 3 Creates hypre ijmatrix from PETSc matrix 4 */ 5 #include <petscsys.h> 6 #include <petsc/private/matimpl.h> 7 #include <petscdmda.h> /*I "petscdmda.h" I*/ 8 #include <../src/dm/impls/da/hypre/mhyp.h> 9 10 /* -----------------------------------------------------------------------------------------------------------------*/ 11 12 /*MC 13 MATHYPRESTRUCT - MATHYPRESTRUCT = "hyprestruct" - A matrix type to be used for parallel sparse matrices 14 based on the hypre HYPRE_StructMatrix. 15 16 Level: intermediate 17 18 Notes: Unlike the more general support for blocks in hypre this allows only one block per process and requires the block 19 be defined by a DMDA. 20 21 The matrix needs a DMDA associated with it by either a call to MatSetDM() or if the matrix is obtained from DMCreateMatrix() 22 23 .seealso: MatCreate(), PCPFMG, MatSetDM(), DMCreateMatrix() 24 M*/ 25 26 27 PetscErrorCode MatSetValuesLocal_HYPREStruct_3d(Mat mat,PetscInt nrow,const PetscInt irow[],PetscInt ncol,const PetscInt icol[],const PetscScalar y[],InsertMode addv) 28 { 29 PetscErrorCode ierr; 30 PetscInt i,j,stencil,index[3],row,entries[7]; 31 const PetscScalar *values = y; 32 Mat_HYPREStruct *ex = (Mat_HYPREStruct*) mat->data; 33 34 PetscFunctionBegin; 35 if (PetscUnlikely(ncol > 7)) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"ncol %D > 7 too large",ncol); 36 for (i=0; i<nrow; i++) { 37 for (j=0; j<ncol; j++) { 38 stencil = icol[j] - irow[i]; 39 if (!stencil) { 40 entries[j] = 3; 41 } else if (stencil == -1) { 42 entries[j] = 2; 43 } else if (stencil == 1) { 44 entries[j] = 4; 45 } else if (stencil == -ex->gnx) { 46 entries[j] = 1; 47 } else if (stencil == ex->gnx) { 48 entries[j] = 5; 49 } else if (stencil == -ex->gnxgny) { 50 entries[j] = 0; 51 } else if (stencil == ex->gnxgny) { 52 entries[j] = 6; 53 } else SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Local row %D local column %D have bad stencil %D",irow[i],icol[j],stencil); 54 } 55 row = ex->gindices[irow[i]] - ex->rstart; 56 index[0] = ex->xs + (row % ex->nx); 57 index[1] = ex->ys + ((row/ex->nx) % ex->ny); 58 index[2] = ex->zs + (row/(ex->nxny)); 59 if (addv == ADD_VALUES) { 60 PetscStackCallStandard(HYPRE_StructMatrixAddToValues,(ex->hmat,(HYPRE_Int *)index,ncol,(HYPRE_Int *)entries,(PetscScalar*)values)); 61 } else { 62 PetscStackCallStandard(HYPRE_StructMatrixSetValues,(ex->hmat,(HYPRE_Int *)index,ncol,(HYPRE_Int *)entries,(PetscScalar*)values)); 63 } 64 values += ncol; 65 } 66 PetscFunctionReturn(0); 67 } 68 69 PetscErrorCode MatZeroRowsLocal_HYPREStruct_3d(Mat mat,PetscInt nrow,const PetscInt irow[],PetscScalar d,Vec x,Vec b) 70 { 71 PetscErrorCode ierr; 72 PetscInt i,index[3],row,entries[7] = {0,1,2,3,4,5,6}; 73 PetscScalar values[7]; 74 Mat_HYPREStruct *ex = (Mat_HYPREStruct*) mat->data; 75 76 PetscFunctionBegin; 77 if (x && b) SETERRQ(PetscObjectComm((PetscObject)mat),PETSC_ERR_SUP,"No support"); 78 ierr = PetscMemzero(values,7*sizeof(PetscScalar));CHKERRQ(ierr); 79 values[3] = d; 80 for (i=0; i<nrow; i++) { 81 row = ex->gindices[irow[i]] - ex->rstart; 82 index[0] = ex->xs + (row % ex->nx); 83 index[1] = ex->ys + ((row/ex->nx) % ex->ny); 84 index[2] = ex->zs + (row/(ex->nxny)); 85 PetscStackCallStandard(HYPRE_StructMatrixSetValues,(ex->hmat,(HYPRE_Int *)index,7,(HYPRE_Int *)entries,values)); 86 } 87 PetscStackCallStandard(HYPRE_StructMatrixAssemble,(ex->hmat)); 88 PetscFunctionReturn(0); 89 } 90 91 PetscErrorCode MatZeroEntries_HYPREStruct_3d(Mat mat) 92 { 93 PetscErrorCode ierr; 94 PetscInt indices[7] = {0,1,2,3,4,5,6}; 95 Mat_HYPREStruct *ex = (Mat_HYPREStruct*) mat->data; 96 97 PetscFunctionBegin; 98 /* hypre has no public interface to do this */ 99 PetscStackCallStandard(hypre_StructMatrixClearBoxValues,(ex->hmat,&ex->hbox,7,(HYPRE_Int *)indices,0,1)); 100 PetscStackCallStandard(HYPRE_StructMatrixAssemble,(ex->hmat)); 101 PetscFunctionReturn(0); 102 } 103 104 static PetscErrorCode MatSetUp_HYPREStruct(Mat mat) 105 { 106 PetscErrorCode ierr; 107 Mat_HYPREStruct *ex = (Mat_HYPREStruct*) mat->data; 108 PetscInt dim,dof,sw[6],nx,ny,nz,ilower[3],iupper[3],ssize,i; 109 DMBoundaryType px,py,pz; 110 DMDAStencilType st; 111 ISLocalToGlobalMapping ltog; 112 DM da; 113 114 PetscFunctionBegin; 115 ierr = MatGetDM(mat,(DM*)&da);CHKERRQ(ierr); 116 ex->da = da; 117 ierr = PetscObjectReference((PetscObject)da);CHKERRQ(ierr); 118 119 ierr = DMDAGetInfo(ex->da,&dim,0,0,0,0,0,0,&dof,&sw[0],&px,&py,&pz,&st);CHKERRQ(ierr); 120 ierr = DMDAGetCorners(ex->da,&ilower[0],&ilower[1],&ilower[2],&iupper[0],&iupper[1],&iupper[2]);CHKERRQ(ierr); 121 iupper[0] += ilower[0] - 1; 122 iupper[1] += ilower[1] - 1; 123 iupper[2] += ilower[2] - 1; 124 125 /* the hypre_Box is used to zero out the matrix entries in MatZeroValues() */ 126 ex->hbox.imin[0] = ilower[0]; 127 ex->hbox.imin[1] = ilower[1]; 128 ex->hbox.imin[2] = ilower[2]; 129 ex->hbox.imax[0] = iupper[0]; 130 ex->hbox.imax[1] = iupper[1]; 131 ex->hbox.imax[2] = iupper[2]; 132 133 /* create the hypre grid object and set its information */ 134 if (dof > 1) SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_SUP,"Currently only support for scalar problems"); 135 if (px || py || pz) SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_SUP,"Ask us to add periodic support by calling HYPRE_StructGridSetPeriodic()"); 136 PetscStackCallStandard(HYPRE_StructGridCreate,(ex->hcomm,dim,&ex->hgrid)); 137 PetscStackCallStandard(HYPRE_StructGridSetExtents,(ex->hgrid,(HYPRE_Int *)ilower,(HYPRE_Int *)iupper)); 138 PetscStackCallStandard(HYPRE_StructGridAssemble,(ex->hgrid)); 139 140 sw[5] = sw[4] = sw[3] = sw[2] = sw[1] = sw[0]; 141 PetscStackCallStandard(HYPRE_StructGridSetNumGhost,(ex->hgrid,(HYPRE_Int *)sw)); 142 143 /* create the hypre stencil object and set its information */ 144 if (sw[0] > 1) SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_SUP,"Ask us to add support for wider stencils"); 145 if (st == DMDA_STENCIL_BOX) SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_SUP,"Ask us to add support for box stencils"); 146 if (dim == 1) { 147 PetscInt offsets[3][1] = {{-1},{0},{1}}; 148 ssize = 3; 149 PetscStackCallStandard(HYPRE_StructStencilCreate,(dim,ssize,&ex->hstencil)); 150 for (i=0; i<ssize; i++) { 151 PetscStackCallStandard(HYPRE_StructStencilSetElement,(ex->hstencil,i,(HYPRE_Int *)offsets[i])); 152 } 153 } else if (dim == 2) { 154 PetscInt offsets[5][2] = {{0,-1},{-1,0},{0,0},{1,0},{0,1}}; 155 ssize = 5; 156 PetscStackCallStandard(HYPRE_StructStencilCreate,(dim,ssize,&ex->hstencil)); 157 for (i=0; i<ssize; i++) { 158 PetscStackCallStandard(HYPRE_StructStencilSetElement,(ex->hstencil,i,(HYPRE_Int *)offsets[i])); 159 } 160 } else if (dim == 3) { 161 PetscInt offsets[7][3] = {{0,0,-1},{0,-1,0},{-1,0,0},{0,0,0},{1,0,0},{0,1,0},{0,0,1}}; 162 ssize = 7; 163 PetscStackCallStandard(HYPRE_StructStencilCreate,(dim,ssize,&ex->hstencil)); 164 for (i=0; i<ssize; i++) { 165 PetscStackCallStandard(HYPRE_StructStencilSetElement,(ex->hstencil,i,(HYPRE_Int *)offsets[i])); 166 } 167 } 168 169 /* create the HYPRE vector for rhs and solution */ 170 PetscStackCallStandard(HYPRE_StructVectorCreate,(ex->hcomm,ex->hgrid,&ex->hb)); 171 PetscStackCallStandard(HYPRE_StructVectorCreate,(ex->hcomm,ex->hgrid,&ex->hx)); 172 PetscStackCallStandard(HYPRE_StructVectorInitialize,(ex->hb)); 173 PetscStackCallStandard(HYPRE_StructVectorInitialize,(ex->hx)); 174 PetscStackCallStandard(HYPRE_StructVectorAssemble,(ex->hb)); 175 PetscStackCallStandard(HYPRE_StructVectorAssemble,(ex->hx)); 176 177 /* create the hypre matrix object and set its information */ 178 PetscStackCallStandard(HYPRE_StructMatrixCreate,(ex->hcomm,ex->hgrid,ex->hstencil,&ex->hmat)); 179 PetscStackCallStandard(HYPRE_StructGridDestroy,(ex->hgrid)); 180 PetscStackCallStandard(HYPRE_StructStencilDestroy,(ex->hstencil)); 181 if (ex->needsinitialization) { 182 PetscStackCallStandard(HYPRE_StructMatrixInitialize,(ex->hmat)); 183 ex->needsinitialization = PETSC_FALSE; 184 } 185 186 /* set the global and local sizes of the matrix */ 187 ierr = DMDAGetCorners(da,0,0,0,&nx,&ny,&nz);CHKERRQ(ierr); 188 ierr = MatSetSizes(mat,dof*nx*ny*nz,dof*nx*ny*nz,PETSC_DECIDE,PETSC_DECIDE);CHKERRQ(ierr); 189 ierr = PetscLayoutSetBlockSize(mat->rmap,1);CHKERRQ(ierr); 190 ierr = PetscLayoutSetBlockSize(mat->cmap,1);CHKERRQ(ierr); 191 ierr = PetscLayoutSetUp(mat->rmap);CHKERRQ(ierr); 192 ierr = PetscLayoutSetUp(mat->cmap);CHKERRQ(ierr); 193 mat->preallocated = PETSC_TRUE; 194 195 if (dim == 3) { 196 mat->ops->setvalueslocal = MatSetValuesLocal_HYPREStruct_3d; 197 mat->ops->zerorowslocal = MatZeroRowsLocal_HYPREStruct_3d; 198 mat->ops->zeroentries = MatZeroEntries_HYPREStruct_3d; 199 200 /* ierr = MatZeroEntries_HYPREStruct_3d(mat);CHKERRQ(ierr); */ 201 } else SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_SUP,"Only support for 3d DMDA currently"); 202 203 /* get values that will be used repeatedly in MatSetValuesLocal() and MatZeroRowsLocal() repeatedly */ 204 ierr = MatGetOwnershipRange(mat,&ex->rstart,NULL);CHKERRQ(ierr); 205 ierr = DMGetLocalToGlobalMapping(ex->da,<og);CHKERRQ(ierr); 206 ierr = ISLocalToGlobalMappingGetIndices(ltog, (const PetscInt **) &ex->gindices);CHKERRQ(ierr); 207 ierr = DMDAGetGhostCorners(ex->da,0,0,0,&ex->gnx,&ex->gnxgny,0);CHKERRQ(ierr); 208 ex->gnxgny *= ex->gnx; 209 ierr = DMDAGetCorners(ex->da,&ex->xs,&ex->ys,&ex->zs,&ex->nx,&ex->ny,0);CHKERRQ(ierr); 210 ex->nxny = ex->nx*ex->ny; 211 PetscFunctionReturn(0); 212 } 213 214 PetscErrorCode MatMult_HYPREStruct(Mat A,Vec x,Vec y) 215 { 216 PetscErrorCode ierr; 217 const PetscScalar *xx; 218 PetscScalar *yy; 219 PetscInt ilower[3],iupper[3]; 220 Mat_HYPREStruct *mx = (Mat_HYPREStruct*)(A->data); 221 222 PetscFunctionBegin; 223 ierr = DMDAGetCorners(mx->da,&ilower[0],&ilower[1],&ilower[2],&iupper[0],&iupper[1],&iupper[2]);CHKERRQ(ierr); 224 225 iupper[0] += ilower[0] - 1; 226 iupper[1] += ilower[1] - 1; 227 iupper[2] += ilower[2] - 1; 228 229 /* copy x values over to hypre */ 230 PetscStackCallStandard(HYPRE_StructVectorSetConstantValues,(mx->hb,0.0)); 231 ierr = VecGetArrayRead(x,&xx);CHKERRQ(ierr); 232 PetscStackCallStandard(HYPRE_StructVectorSetBoxValues,(mx->hb,(HYPRE_Int *)ilower,(HYPRE_Int *)iupper,(PetscScalar*)xx)); 233 ierr = VecRestoreArrayRead(x,&xx);CHKERRQ(ierr); 234 PetscStackCallStandard(HYPRE_StructVectorAssemble,(mx->hb)); 235 PetscStackCallStandard(HYPRE_StructMatrixMatvec,(1.0,mx->hmat,mx->hb,0.0,mx->hx)); 236 237 /* copy solution values back to PETSc */ 238 ierr = VecGetArray(y,&yy);CHKERRQ(ierr); 239 PetscStackCallStandard(HYPRE_StructVectorGetBoxValues,(mx->hx,(HYPRE_Int *)ilower,(HYPRE_Int *)iupper,yy)); 240 ierr = VecRestoreArray(y,&yy);CHKERRQ(ierr); 241 PetscFunctionReturn(0); 242 } 243 244 PetscErrorCode MatAssemblyEnd_HYPREStruct(Mat mat,MatAssemblyType mode) 245 { 246 Mat_HYPREStruct *ex = (Mat_HYPREStruct*) mat->data; 247 PetscErrorCode ierr; 248 249 PetscFunctionBegin; 250 PetscStackCallStandard(HYPRE_StructMatrixAssemble,(ex->hmat)); 251 /* PetscStackCallStandard(HYPRE_StructMatrixPrint,("dummy",ex->hmat,0)); */ 252 PetscFunctionReturn(0); 253 } 254 255 PetscErrorCode MatZeroEntries_HYPREStruct(Mat mat) 256 { 257 PetscFunctionBegin; 258 /* before the DMDA is set to the matrix the zero doesn't need to do anything */ 259 PetscFunctionReturn(0); 260 } 261 262 PetscErrorCode MatDestroy_HYPREStruct(Mat mat) 263 { 264 Mat_HYPREStruct *ex = (Mat_HYPREStruct*) mat->data; 265 PetscErrorCode ierr; 266 267 PetscFunctionBegin; 268 PetscStackCallStandard(HYPRE_StructMatrixDestroy,(ex->hmat)); 269 PetscStackCallStandard(HYPRE_StructVectorDestroy,(ex->hx)); 270 PetscStackCallStandard(HYPRE_StructVectorDestroy,(ex->hb)); 271 ierr = PetscObjectDereference((PetscObject)ex->da);CHKERRQ(ierr); 272 ierr = MPI_Comm_free(&(ex->hcomm));CHKERRQ(ierr); 273 ierr = PetscFree(ex);CHKERRQ(ierr); 274 PetscFunctionReturn(0); 275 } 276 277 PETSC_EXTERN PetscErrorCode MatCreate_HYPREStruct(Mat B) 278 { 279 Mat_HYPREStruct *ex; 280 PetscErrorCode ierr; 281 282 PetscFunctionBegin; 283 ierr = PetscNewLog(B,&ex);CHKERRQ(ierr); 284 B->data = (void*)ex; 285 B->rmap->bs = 1; 286 B->assembled = PETSC_FALSE; 287 288 B->insertmode = NOT_SET_VALUES; 289 290 B->ops->assemblyend = MatAssemblyEnd_HYPREStruct; 291 B->ops->mult = MatMult_HYPREStruct; 292 B->ops->zeroentries = MatZeroEntries_HYPREStruct; 293 B->ops->destroy = MatDestroy_HYPREStruct; 294 B->ops->setup = MatSetUp_HYPREStruct; 295 296 ex->needsinitialization = PETSC_TRUE; 297 298 ierr = MPI_Comm_dup(PetscObjectComm((PetscObject)B),&(ex->hcomm));CHKERRQ(ierr); 299 ierr = PetscObjectChangeTypeName((PetscObject)B,MATHYPRESTRUCT);CHKERRQ(ierr); 300 PetscFunctionReturn(0); 301 } 302 303 /*MC 304 MATHYPRESSTRUCT - MATHYPRESSTRUCT = "hypresstruct" - A matrix type to be used for parallel sparse matrices 305 based on the hypre HYPRE_SStructMatrix. 306 307 308 Level: intermediate 309 310 Notes: Unlike hypre's general semi-struct object consisting of a collection of structured-grid objects and unstructured 311 grid objects, we restrict the semi-struct objects to consist of only structured-grid components. 312 313 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 314 be defined by a DMDA. 315 316 The matrix needs a DMDA associated with it by either a call to MatSetDM() or if the matrix is obtained from DMCreateMatrix() 317 318 M*/ 319 320 PetscErrorCode MatSetValuesLocal_HYPRESStruct_3d(Mat mat,PetscInt nrow,const PetscInt irow[],PetscInt ncol,const PetscInt icol[],const PetscScalar y[],InsertMode addv) 321 { 322 PetscErrorCode ierr; 323 PetscInt i,j,stencil,index[3]; 324 const PetscScalar *values = y; 325 Mat_HYPRESStruct *ex = (Mat_HYPRESStruct*) mat->data; 326 327 PetscInt part = 0; /* Petsc sstruct interface only allows 1 part */ 328 PetscInt ordering; 329 PetscInt grid_rank, to_grid_rank; 330 PetscInt var_type, to_var_type; 331 PetscInt to_var_entry = 0; 332 PetscInt nvars= ex->nvars; 333 PetscInt row,*entries; 334 335 PetscFunctionBegin; 336 ierr = PetscMalloc1(7*nvars,&entries);CHKERRQ(ierr); 337 ordering = ex->dofs_order; /* ordering= 0 nodal ordering 338 1 variable ordering */ 339 /* stencil entries are orderer by variables: var0_stencil0, var0_stencil1, ..., var0_stencil6, var1_stencil0, var1_stencil1, ... */ 340 341 if (!ordering) { 342 for (i=0; i<nrow; i++) { 343 grid_rank= irow[i]/nvars; 344 var_type = (irow[i] % nvars); 345 346 for (j=0; j<ncol; j++) { 347 to_grid_rank= icol[j]/nvars; 348 to_var_type = (icol[j] % nvars); 349 350 to_var_entry= to_var_entry*7; 351 entries[j] = to_var_entry; 352 353 stencil = to_grid_rank-grid_rank; 354 if (!stencil) { 355 entries[j] += 3; 356 } else if (stencil == -1) { 357 entries[j] += 2; 358 } else if (stencil == 1) { 359 entries[j] += 4; 360 } else if (stencil == -ex->gnx) { 361 entries[j] += 1; 362 } else if (stencil == ex->gnx) { 363 entries[j] += 5; 364 } else if (stencil == -ex->gnxgny) { 365 entries[j] += 0; 366 } else if (stencil == ex->gnxgny) { 367 entries[j] += 6; 368 } else SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Local row %D local column %D have bad stencil %D",irow[i],icol[j],stencil); 369 } 370 371 row = ex->gindices[grid_rank] - ex->rstart; 372 index[0] = ex->xs + (row % ex->nx); 373 index[1] = ex->ys + ((row/ex->nx) % ex->ny); 374 index[2] = ex->zs + (row/(ex->nxny)); 375 376 if (addv == ADD_VALUES) { 377 PetscStackCallStandard(HYPRE_SStructMatrixAddToValues,(ex->ss_mat,part,(HYPRE_Int *)index,var_type,ncol,(HYPRE_Int *)entries,(PetscScalar*)values)); 378 } else { 379 PetscStackCallStandard(HYPRE_SStructMatrixSetValues,(ex->ss_mat,part,(HYPRE_Int *)index,var_type,ncol,(HYPRE_Int *)entries,(PetscScalar*)values)); 380 } 381 values += ncol; 382 } 383 } else { 384 for (i=0; i<nrow; i++) { 385 var_type = irow[i]/(ex->gnxgnygnz); 386 grid_rank= irow[i] - var_type*(ex->gnxgnygnz); 387 388 for (j=0; j<ncol; j++) { 389 to_var_type = icol[j]/(ex->gnxgnygnz); 390 to_grid_rank= icol[j] - to_var_type*(ex->gnxgnygnz); 391 392 to_var_entry= to_var_entry*7; 393 entries[j] = to_var_entry; 394 395 stencil = to_grid_rank-grid_rank; 396 if (!stencil) { 397 entries[j] += 3; 398 } else if (stencil == -1) { 399 entries[j] += 2; 400 } else if (stencil == 1) { 401 entries[j] += 4; 402 } else if (stencil == -ex->gnx) { 403 entries[j] += 1; 404 } else if (stencil == ex->gnx) { 405 entries[j] += 5; 406 } else if (stencil == -ex->gnxgny) { 407 entries[j] += 0; 408 } else if (stencil == ex->gnxgny) { 409 entries[j] += 6; 410 } else SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Local row %D local column %D have bad stencil %D",irow[i],icol[j],stencil); 411 } 412 413 row = ex->gindices[grid_rank] - ex->rstart; 414 index[0] = ex->xs + (row % ex->nx); 415 index[1] = ex->ys + ((row/ex->nx) % ex->ny); 416 index[2] = ex->zs + (row/(ex->nxny)); 417 418 if (addv == ADD_VALUES) { 419 PetscStackCallStandard(HYPRE_SStructMatrixAddToValues,(ex->ss_mat,part,(HYPRE_Int *)index,var_type,ncol,(HYPRE_Int *)entries,(PetscScalar*)values)); 420 } else { 421 PetscStackCallStandard(HYPRE_SStructMatrixSetValues,(ex->ss_mat,part,(HYPRE_Int *)index,var_type,ncol,(HYPRE_Int *)entries,(PetscScalar*)values)); 422 } 423 values += ncol; 424 } 425 } 426 ierr = PetscFree(entries);CHKERRQ(ierr); 427 PetscFunctionReturn(0); 428 } 429 430 PetscErrorCode MatZeroRowsLocal_HYPRESStruct_3d(Mat mat,PetscInt nrow,const PetscInt irow[],PetscScalar d,Vec x,Vec b) 431 { 432 PetscErrorCode ierr; 433 PetscInt i,index[3]; 434 PetscScalar **values; 435 Mat_HYPRESStruct *ex = (Mat_HYPRESStruct*) mat->data; 436 437 PetscInt part = 0; /* Petsc sstruct interface only allows 1 part */ 438 PetscInt ordering = ex->dofs_order; 439 PetscInt grid_rank; 440 PetscInt var_type; 441 PetscInt nvars= ex->nvars; 442 PetscInt row,*entries; 443 444 PetscFunctionBegin; 445 if (x && b) SETERRQ(PetscObjectComm((PetscObject)mat),PETSC_ERR_SUP,"No support"); 446 ierr = PetscMalloc1(7*nvars,&entries);CHKERRQ(ierr); 447 448 ierr = PetscMalloc1(nvars,&values);CHKERRQ(ierr); 449 ierr = PetscMalloc1(7*nvars*nvars,&values[0]);CHKERRQ(ierr); 450 for (i=1; i<nvars; i++) { 451 values[i] = values[i-1] + nvars*7; 452 } 453 454 for (i=0; i< nvars; i++) { 455 ierr = PetscMemzero(values[i],nvars*7*sizeof(PetscScalar));CHKERRQ(ierr); 456 *(values[i]+3) = d; 457 } 458 459 for (i= 0; i< nvars*7; i++) entries[i] = i; 460 461 if (!ordering) { 462 for (i=0; i<nrow; i++) { 463 grid_rank = irow[i] / nvars; 464 var_type =(irow[i] % nvars); 465 466 row = ex->gindices[grid_rank] - ex->rstart; 467 index[0] = ex->xs + (row % ex->nx); 468 index[1] = ex->ys + ((row/ex->nx) % ex->ny); 469 index[2] = ex->zs + (row/(ex->nxny)); 470 PetscStackCallStandard(HYPRE_SStructMatrixSetValues,(ex->ss_mat,part,(HYPRE_Int *)index,var_type,7*nvars,(HYPRE_Int *)entries,values[var_type])); 471 } 472 } else { 473 for (i=0; i<nrow; i++) { 474 var_type = irow[i]/(ex->gnxgnygnz); 475 grid_rank= irow[i] - var_type*(ex->gnxgnygnz); 476 477 row = ex->gindices[grid_rank] - ex->rstart; 478 index[0] = ex->xs + (row % ex->nx); 479 index[1] = ex->ys + ((row/ex->nx) % ex->ny); 480 index[2] = ex->zs + (row/(ex->nxny)); 481 PetscStackCallStandard(HYPRE_SStructMatrixSetValues,(ex->ss_mat,part,(HYPRE_Int *)index,var_type,7*nvars,(HYPRE_Int *)entries,values[var_type])); 482 } 483 } 484 PetscStackCallStandard(HYPRE_SStructMatrixAssemble,(ex->ss_mat)); 485 ierr = PetscFree(values[0]);CHKERRQ(ierr); 486 ierr = PetscFree(values);CHKERRQ(ierr); 487 ierr = PetscFree(entries);CHKERRQ(ierr); 488 PetscFunctionReturn(0); 489 } 490 491 PetscErrorCode MatZeroEntries_HYPRESStruct_3d(Mat mat) 492 { 493 PetscErrorCode ierr; 494 Mat_HYPRESStruct *ex = (Mat_HYPRESStruct*) mat->data; 495 PetscInt nvars= ex->nvars; 496 PetscInt size; 497 PetscInt part= 0; /* only one part */ 498 499 PetscFunctionBegin; 500 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); 501 { 502 PetscInt i,*entries,iupper[3],ilower[3]; 503 PetscScalar *values; 504 505 506 for (i= 0; i< 3; i++) { 507 ilower[i]= ex->hbox.imin[i]; 508 iupper[i]= ex->hbox.imax[i]; 509 } 510 511 ierr = PetscMalloc2(nvars*7,&entries,nvars*7*size,&values);CHKERRQ(ierr); 512 for (i= 0; i< nvars*7; i++) entries[i]= i; 513 ierr = PetscMemzero(values,nvars*7*size*sizeof(PetscScalar));CHKERRQ(ierr); 514 515 for (i= 0; i< nvars; i++) { 516 PetscStackCallStandard(HYPRE_SStructMatrixSetBoxValues,(ex->ss_mat,part,(HYPRE_Int *)ilower,(HYPRE_Int *)iupper,i,nvars*7,(HYPRE_Int *)entries,values)); 517 } 518 ierr = PetscFree2(entries,values);CHKERRQ(ierr); 519 } 520 PetscStackCallStandard(HYPRE_SStructMatrixAssemble,(ex->ss_mat)); 521 PetscFunctionReturn(0); 522 } 523 524 525 static PetscErrorCode MatSetUp_HYPRESStruct(Mat mat) 526 { 527 PetscErrorCode ierr; 528 Mat_HYPRESStruct *ex = (Mat_HYPRESStruct*) mat->data; 529 PetscInt dim,dof,sw[3],nx,ny,nz; 530 PetscInt ilower[3],iupper[3],ssize,i; 531 DMBoundaryType px,py,pz; 532 DMDAStencilType st; 533 PetscInt nparts= 1; /* assuming only one part */ 534 PetscInt part = 0; 535 ISLocalToGlobalMapping ltog; 536 DM da; 537 538 PetscFunctionBegin; 539 ierr = MatGetDM(mat,(DM*)&da);CHKERRQ(ierr); 540 ex->da = da; 541 ierr = PetscObjectReference((PetscObject)da);CHKERRQ(ierr); 542 543 ierr = DMDAGetInfo(ex->da,&dim,0,0,0,0,0,0,&dof,&sw[0],&px,&py,&pz,&st);CHKERRQ(ierr); 544 ierr = DMDAGetCorners(ex->da,&ilower[0],&ilower[1],&ilower[2],&iupper[0],&iupper[1],&iupper[2]);CHKERRQ(ierr); 545 iupper[0] += ilower[0] - 1; 546 iupper[1] += ilower[1] - 1; 547 iupper[2] += ilower[2] - 1; 548 /* the hypre_Box is used to zero out the matrix entries in MatZeroValues() */ 549 ex->hbox.imin[0] = ilower[0]; 550 ex->hbox.imin[1] = ilower[1]; 551 ex->hbox.imin[2] = ilower[2]; 552 ex->hbox.imax[0] = iupper[0]; 553 ex->hbox.imax[1] = iupper[1]; 554 ex->hbox.imax[2] = iupper[2]; 555 556 ex->dofs_order = 0; 557 558 /* assuming that the same number of dofs on each gridpoint. Also assume all cell-centred based */ 559 ex->nvars= dof; 560 561 /* create the hypre grid object and set its information */ 562 if (px || py || pz) SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_SUP,"Ask us to add periodic support by calling HYPRE_SStructGridSetPeriodic()"); 563 PetscStackCallStandard(HYPRE_SStructGridCreate,(ex->hcomm,dim,nparts,&ex->ss_grid)); 564 PetscStackCallStandard(HYPRE_SStructGridSetExtents,(ex->ss_grid,part,ex->hbox.imin,ex->hbox.imax)); 565 { 566 HYPRE_SStructVariable *vartypes; 567 ierr = PetscMalloc1(ex->nvars,&vartypes);CHKERRQ(ierr); 568 for (i= 0; i< ex->nvars; i++) vartypes[i]= HYPRE_SSTRUCT_VARIABLE_CELL; 569 PetscStackCallStandard(HYPRE_SStructGridSetVariables,(ex->ss_grid, part, ex->nvars,vartypes)); 570 ierr = PetscFree(vartypes);CHKERRQ(ierr); 571 } 572 PetscStackCallStandard(HYPRE_SStructGridAssemble,(ex->ss_grid)); 573 574 sw[1] = sw[0]; 575 sw[2] = sw[1]; 576 /* PetscStackCallStandard(HYPRE_SStructGridSetNumGhost,(ex->ss_grid,sw)); */ 577 578 /* create the hypre stencil object and set its information */ 579 if (sw[0] > 1) SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_SUP,"Ask us to add support for wider stencils"); 580 if (st == DMDA_STENCIL_BOX) SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_SUP,"Ask us to add support for box stencils"); 581 582 if (dim == 1) { 583 PetscInt offsets[3][1] = {{-1},{0},{1}}; 584 PetscInt j, cnt; 585 586 ssize = 3*(ex->nvars); 587 PetscStackCallStandard(HYPRE_SStructStencilCreate,(dim,ssize,&ex->ss_stencil)); 588 cnt= 0; 589 for (i = 0; i < (ex->nvars); i++) { 590 for (j = 0; j < 3; j++) { 591 PetscStackCallStandard(HYPRE_SStructStencilSetEntry,(ex->ss_stencil, cnt, (HYPRE_Int *)offsets[j], i)); 592 cnt++; 593 } 594 } 595 } else if (dim == 2) { 596 PetscInt offsets[5][2] = {{0,-1},{-1,0},{0,0},{1,0},{0,1}}; 597 PetscInt j, cnt; 598 599 ssize = 5*(ex->nvars); 600 PetscStackCallStandard(HYPRE_SStructStencilCreate,(dim,ssize,&ex->ss_stencil)); 601 cnt= 0; 602 for (i= 0; i< (ex->nvars); i++) { 603 for (j= 0; j< 5; j++) { 604 PetscStackCallStandard(HYPRE_SStructStencilSetEntry,(ex->ss_stencil, cnt, (HYPRE_Int *)offsets[j], i)); 605 cnt++; 606 } 607 } 608 } else if (dim == 3) { 609 PetscInt offsets[7][3] = {{0,0,-1},{0,-1,0},{-1,0,0},{0,0,0},{1,0,0},{0,1,0},{0,0,1}}; 610 PetscInt j, cnt; 611 612 ssize = 7*(ex->nvars); 613 PetscStackCallStandard(HYPRE_SStructStencilCreate,(dim,ssize,&ex->ss_stencil)); 614 cnt= 0; 615 for (i= 0; i< (ex->nvars); i++) { 616 for (j= 0; j< 7; j++) { 617 PetscStackCallStandard(HYPRE_SStructStencilSetEntry,(ex->ss_stencil, cnt, (HYPRE_Int *)offsets[j], i)); 618 cnt++; 619 } 620 } 621 } 622 623 /* create the HYPRE graph */ 624 PetscStackCallStandard(HYPRE_SStructGraphCreate,(ex->hcomm, ex->ss_grid, &(ex->ss_graph))); 625 626 /* set the stencil graph. Note that each variable has the same graph. This means that each 627 variable couples to all the other variable and with the same stencil pattern. */ 628 for (i= 0; i< (ex->nvars); i++) { 629 PetscStackCallStandard(HYPRE_SStructGraphSetStencil,(ex->ss_graph,part,i,ex->ss_stencil)); 630 } 631 PetscStackCallStandard(HYPRE_SStructGraphAssemble,(ex->ss_graph)); 632 633 /* create the HYPRE sstruct vectors for rhs and solution */ 634 PetscStackCallStandard(HYPRE_SStructVectorCreate,(ex->hcomm,ex->ss_grid,&ex->ss_b)); 635 PetscStackCallStandard(HYPRE_SStructVectorCreate,(ex->hcomm,ex->ss_grid,&ex->ss_x)); 636 PetscStackCallStandard(HYPRE_SStructVectorInitialize,(ex->ss_b)); 637 PetscStackCallStandard(HYPRE_SStructVectorInitialize,(ex->ss_x)); 638 PetscStackCallStandard(HYPRE_SStructVectorAssemble,(ex->ss_b)); 639 PetscStackCallStandard(HYPRE_SStructVectorAssemble,(ex->ss_x)); 640 641 /* create the hypre matrix object and set its information */ 642 PetscStackCallStandard(HYPRE_SStructMatrixCreate,(ex->hcomm,ex->ss_graph,&ex->ss_mat)); 643 PetscStackCallStandard(HYPRE_SStructGridDestroy,(ex->ss_grid)); 644 PetscStackCallStandard(HYPRE_SStructStencilDestroy,(ex->ss_stencil)); 645 if (ex->needsinitialization) { 646 PetscStackCallStandard(HYPRE_SStructMatrixInitialize,(ex->ss_mat)); 647 ex->needsinitialization = PETSC_FALSE; 648 } 649 650 /* set the global and local sizes of the matrix */ 651 ierr = DMDAGetCorners(da,0,0,0,&nx,&ny,&nz);CHKERRQ(ierr); 652 ierr = MatSetSizes(mat,dof*nx*ny*nz,dof*nx*ny*nz,PETSC_DECIDE,PETSC_DECIDE);CHKERRQ(ierr); 653 ierr = PetscLayoutSetBlockSize(mat->rmap,1);CHKERRQ(ierr); 654 ierr = PetscLayoutSetBlockSize(mat->cmap,1);CHKERRQ(ierr); 655 ierr = PetscLayoutSetUp(mat->rmap);CHKERRQ(ierr); 656 ierr = PetscLayoutSetUp(mat->cmap);CHKERRQ(ierr); 657 658 if (dim == 3) { 659 mat->ops->setvalueslocal = MatSetValuesLocal_HYPRESStruct_3d; 660 mat->ops->zerorowslocal = MatZeroRowsLocal_HYPRESStruct_3d; 661 mat->ops->zeroentries = MatZeroEntries_HYPRESStruct_3d; 662 663 ierr = MatZeroEntries_HYPRESStruct_3d(mat);CHKERRQ(ierr); 664 } else SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_SUP,"Only support for 3d DMDA currently"); 665 666 /* get values that will be used repeatedly in MatSetValuesLocal() and MatZeroRowsLocal() repeatedly */ 667 ierr = MatGetOwnershipRange(mat,&ex->rstart,NULL);CHKERRQ(ierr); 668 ierr = DMGetLocalToGlobalMapping(ex->da,<og);CHKERRQ(ierr); 669 ierr = ISLocalToGlobalMappingGetIndices(ltog, (const PetscInt **) &ex->gindices);CHKERRQ(ierr); 670 ierr = DMDAGetGhostCorners(ex->da,0,0,0,&ex->gnx,&ex->gnxgny,&ex->gnxgnygnz);CHKERRQ(ierr); 671 672 ex->gnxgny *= ex->gnx; 673 ex->gnxgnygnz *= ex->gnxgny; 674 675 ierr = DMDAGetCorners(ex->da,&ex->xs,&ex->ys,&ex->zs,&ex->nx,&ex->ny,&ex->nz);CHKERRQ(ierr); 676 677 ex->nxny = ex->nx*ex->ny; 678 ex->nxnynz = ex->nz*ex->nxny; 679 PetscFunctionReturn(0); 680 } 681 682 PetscErrorCode MatMult_HYPRESStruct(Mat A,Vec x,Vec y) 683 { 684 PetscErrorCode ierr; 685 const PetscScalar *xx; 686 PetscScalar *yy; 687 PetscInt ilower[3],iupper[3]; 688 Mat_HYPRESStruct *mx = (Mat_HYPRESStruct*)(A->data); 689 PetscInt ordering= mx->dofs_order; 690 PetscInt nvars = mx->nvars; 691 PetscInt part = 0; 692 PetscInt size; 693 PetscInt i; 694 695 PetscFunctionBegin; 696 ierr = DMDAGetCorners(mx->da,&ilower[0],&ilower[1],&ilower[2],&iupper[0],&iupper[1],&iupper[2]);CHKERRQ(ierr); 697 iupper[0] += ilower[0] - 1; 698 iupper[1] += ilower[1] - 1; 699 iupper[2] += ilower[2] - 1; 700 701 size= 1; 702 for (i = 0; i < 3; i++) size *= (iupper[i]-ilower[i]+1); 703 704 /* copy x values over to hypre for variable ordering */ 705 if (ordering) { 706 PetscStackCallStandard(HYPRE_SStructVectorSetConstantValues,(mx->ss_b,0.0)); 707 ierr = VecGetArrayRead(x,&xx);CHKERRQ(ierr); 708 for (i= 0; i< nvars; i++) { 709 PetscStackCallStandard(HYPRE_SStructVectorSetBoxValues,(mx->ss_b,part,(HYPRE_Int *)ilower,(HYPRE_Int *)iupper,i,(PetscScalar*)xx+(size*i))); 710 } 711 ierr = VecRestoreArrayRead(x,&xx);CHKERRQ(ierr); 712 PetscStackCallStandard(HYPRE_SStructVectorAssemble,(mx->ss_b)); 713 PetscStackCallStandard(HYPRE_SStructMatrixMatvec,(1.0,mx->ss_mat,mx->ss_b,0.0,mx->ss_x)); 714 715 /* copy solution values back to PETSc */ 716 ierr = VecGetArray(y,&yy);CHKERRQ(ierr); 717 for (i= 0; i< nvars; i++) { 718 PetscStackCallStandard(HYPRE_SStructVectorGetBoxValues,(mx->ss_x,part,(HYPRE_Int *)ilower,(HYPRE_Int *)iupper,i,yy+(size*i))); 719 } 720 ierr = VecRestoreArray(y,&yy);CHKERRQ(ierr); 721 } else { /* nodal ordering must be mapped to variable ordering for sys_pfmg */ 722 PetscScalar *z; 723 PetscInt j, k; 724 725 ierr = PetscMalloc1(nvars*size,&z);CHKERRQ(ierr); 726 PetscStackCallStandard(HYPRE_SStructVectorSetConstantValues,(mx->ss_b,0.0)); 727 ierr = VecGetArrayRead(x,&xx);CHKERRQ(ierr); 728 729 /* transform nodal to hypre's variable ordering for sys_pfmg */ 730 for (i= 0; i< size; i++) { 731 k= i*nvars; 732 for (j= 0; j< nvars; j++) z[j*size+i]= xx[k+j]; 733 } 734 for (i= 0; i< nvars; i++) { 735 PetscStackCallStandard(HYPRE_SStructVectorSetBoxValues,(mx->ss_b,part,(HYPRE_Int *)ilower,(HYPRE_Int *)iupper,i,z+(size*i))); 736 } 737 ierr = VecRestoreArrayRead(x,&xx);CHKERRQ(ierr); 738 PetscStackCallStandard(HYPRE_SStructVectorAssemble,(mx->ss_b)); 739 PetscStackCallStandard(HYPRE_SStructMatrixMatvec,(1.0,mx->ss_mat,mx->ss_b,0.0,mx->ss_x)); 740 741 /* copy solution values back to PETSc */ 742 ierr = VecGetArray(y,&yy);CHKERRQ(ierr); 743 for (i= 0; i< nvars; i++) { 744 PetscStackCallStandard(HYPRE_SStructVectorGetBoxValues,(mx->ss_x,part,(HYPRE_Int *)ilower,(HYPRE_Int *)iupper,i,z+(size*i))); 745 } 746 /* transform hypre's variable ordering for sys_pfmg to nodal ordering */ 747 for (i= 0; i< size; i++) { 748 k= i*nvars; 749 for (j= 0; j< nvars; j++) yy[k+j]= z[j*size+i]; 750 } 751 ierr = VecRestoreArray(y,&yy);CHKERRQ(ierr); 752 ierr = PetscFree(z);CHKERRQ(ierr); 753 } 754 PetscFunctionReturn(0); 755 } 756 757 PetscErrorCode MatAssemblyEnd_HYPRESStruct(Mat mat,MatAssemblyType mode) 758 { 759 Mat_HYPRESStruct *ex = (Mat_HYPRESStruct*) mat->data; 760 PetscErrorCode ierr; 761 762 PetscFunctionBegin; 763 PetscStackCallStandard(HYPRE_SStructMatrixAssemble,(ex->ss_mat)); 764 PetscFunctionReturn(0); 765 } 766 767 PetscErrorCode MatZeroEntries_HYPRESStruct(Mat mat) 768 { 769 PetscFunctionBegin; 770 /* before the DMDA is set to the matrix the zero doesn't need to do anything */ 771 PetscFunctionReturn(0); 772 } 773 774 PetscErrorCode MatDestroy_HYPRESStruct(Mat mat) 775 { 776 Mat_HYPRESStruct *ex = (Mat_HYPRESStruct*) mat->data; 777 PetscErrorCode ierr; 778 779 PetscFunctionBegin; 780 PetscStackCallStandard(HYPRE_SStructGraphDestroy,(ex->ss_graph)); 781 PetscStackCallStandard(HYPRE_SStructMatrixDestroy,(ex->ss_mat)); 782 PetscStackCallStandard(HYPRE_SStructVectorDestroy,(ex->ss_x)); 783 PetscStackCallStandard(HYPRE_SStructVectorDestroy,(ex->ss_b)); 784 ierr = PetscObjectDereference((PetscObject)ex->da);CHKERRQ(ierr); 785 ierr = MPI_Comm_free(&(ex->hcomm));CHKERRQ(ierr); 786 ierr = PetscFree(ex);CHKERRQ(ierr); 787 PetscFunctionReturn(0); 788 } 789 790 PETSC_EXTERN PetscErrorCode MatCreate_HYPRESStruct(Mat B) 791 { 792 Mat_HYPRESStruct *ex; 793 PetscErrorCode ierr; 794 795 PetscFunctionBegin; 796 ierr = PetscNewLog(B,&ex);CHKERRQ(ierr); 797 B->data = (void*)ex; 798 B->rmap->bs = 1; 799 B->assembled = PETSC_FALSE; 800 801 B->insertmode = NOT_SET_VALUES; 802 803 B->ops->assemblyend = MatAssemblyEnd_HYPRESStruct; 804 B->ops->mult = MatMult_HYPRESStruct; 805 B->ops->zeroentries = MatZeroEntries_HYPRESStruct; 806 B->ops->destroy = MatDestroy_HYPRESStruct; 807 B->ops->setup = MatSetUp_HYPRESStruct; 808 809 ex->needsinitialization = PETSC_TRUE; 810 811 ierr = MPI_Comm_dup(PetscObjectComm((PetscObject)B),&(ex->hcomm));CHKERRQ(ierr); 812 ierr = PetscObjectChangeTypeName((PetscObject)B,MATHYPRESSTRUCT);CHKERRQ(ierr); 813 PetscFunctionReturn(0); 814 } 815 816 817 818