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