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