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