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