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