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