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