1 2 /* 3 Creates hypre ijmatrix from PETSc matrix 4 */ 5 #include <petscsys.h> 6 #include <petsc/private/petschypre.h> 7 #include <petsc/private/matimpl.h> 8 #include <petscdmda.h> /*I "petscdmda.h" I*/ 9 #include <../src/dm/impls/da/hypre/mhyp.h> 10 11 /* -----------------------------------------------------------------------------------------------------------------*/ 12 13 /*MC 14 MATHYPRESTRUCT - MATHYPRESTRUCT = "hyprestruct" - A matrix type to be used for parallel sparse matrices 15 based on the hypre HYPRE_StructMatrix. 16 17 Level: intermediate 18 19 Notes: 20 Unlike the more general support for blocks in hypre this allows only one block per process and requires the block 21 be defined by a `DMDA`. 22 23 The matrix needs a `DMDA` associated with it by either a call to `MatSetDM()` or if the matrix is obtained from `DMCreateMatrix()` 24 25 .seealso: `MatCreate()`, `PCPFMG`, `MatSetDM()`, `DMCreateMatrix()` 26 M*/ 27 28 PetscErrorCode MatSetValuesLocal_HYPREStruct_3d(Mat mat, PetscInt nrow, const PetscInt irow[], PetscInt ncol, const PetscInt icol[], const PetscScalar y[], InsertMode addv) 29 { 30 HYPRE_Int index[3], entries[9]; 31 PetscInt i, j, stencil, row; 32 HYPRE_Complex *values = (HYPRE_Complex *)y; 33 Mat_HYPREStruct *ex = (Mat_HYPREStruct *)mat->data; 34 35 PetscFunctionBegin; 36 PetscCheck(ncol <= 9, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "ncol %" PetscInt_FMT " > 9 too large", ncol); 37 for (i = 0; i < nrow; i++) { 38 for (j = 0; j < ncol; j++) { 39 stencil = icol[j] - irow[i]; 40 if (!stencil) { 41 entries[j] = 3; 42 } else if (stencil == -1) { 43 entries[j] = 2; 44 } else if (stencil == 1) { 45 entries[j] = 4; 46 } else if (stencil == -ex->gnx) { 47 entries[j] = 1; 48 } else if (stencil == ex->gnx) { 49 entries[j] = 5; 50 } else if (stencil == -ex->gnxgny) { 51 entries[j] = 0; 52 } else if (stencil == ex->gnxgny) { 53 entries[j] = 6; 54 } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Local row %" PetscInt_FMT " local column %" PetscInt_FMT " have bad stencil %" PetscInt_FMT, irow[i], icol[j], stencil); 55 } 56 row = ex->gindices[irow[i]] - ex->rstart; 57 index[0] = (HYPRE_Int)(ex->xs + (row % ex->nx)); 58 index[1] = (HYPRE_Int)(ex->ys + ((row / ex->nx) % ex->ny)); 59 index[2] = (HYPRE_Int)(ex->zs + (row / (ex->nxny))); 60 if (addv == ADD_VALUES) PetscCallExternal(HYPRE_StructMatrixAddToValues, ex->hmat, index, ncol, entries, values); 61 else PetscCallExternal(HYPRE_StructMatrixSetValues, ex->hmat, index, ncol, entries, values); 62 values += ncol; 63 } 64 PetscFunctionReturn(PETSC_SUCCESS); 65 } 66 67 PetscErrorCode MatZeroRowsLocal_HYPREStruct_3d(Mat mat, PetscInt nrow, const PetscInt irow[], PetscScalar d, Vec x, Vec b) 68 { 69 HYPRE_Int index[3], entries[7] = {0, 1, 2, 3, 4, 5, 6}; 70 PetscInt row, i; 71 HYPRE_Complex values[7]; 72 Mat_HYPREStruct *ex = (Mat_HYPREStruct *)mat->data; 73 74 PetscFunctionBegin; 75 PetscCheck(!x || !b, PetscObjectComm((PetscObject)mat), PETSC_ERR_SUP, "No support"); 76 PetscCall(PetscArrayzero(values, 7)); 77 PetscCall(PetscHYPREScalarCast(d, &values[3])); 78 for (i = 0; i < nrow; i++) { 79 row = ex->gindices[irow[i]] - ex->rstart; 80 index[0] = (HYPRE_Int)(ex->xs + (row % ex->nx)); 81 index[1] = (HYPRE_Int)(ex->ys + ((row / ex->nx) % ex->ny)); 82 index[2] = (HYPRE_Int)(ex->zs + (row / (ex->nxny))); 83 PetscCallExternal(HYPRE_StructMatrixSetValues, ex->hmat, index, 7, entries, values); 84 } 85 PetscCallExternal(HYPRE_StructMatrixAssemble, ex->hmat); 86 PetscFunctionReturn(PETSC_SUCCESS); 87 } 88 89 PetscErrorCode MatZeroEntries_HYPREStruct_3d(Mat mat) 90 { 91 HYPRE_Int indices[7] = {0, 1, 2, 3, 4, 5, 6}; 92 Mat_HYPREStruct *ex = (Mat_HYPREStruct *)mat->data; 93 94 PetscFunctionBegin; 95 /* hypre has no public interface to do this */ 96 PetscCallExternal(hypre_StructMatrixClearBoxValues, ex->hmat, &ex->hbox, 7, indices, 0, 1); 97 PetscCallExternal(HYPRE_StructMatrixAssemble, ex->hmat); 98 PetscFunctionReturn(PETSC_SUCCESS); 99 } 100 101 static PetscErrorCode MatSetUp_HYPREStruct(Mat mat) 102 { 103 Mat_HYPREStruct *ex = (Mat_HYPREStruct *)mat->data; 104 HYPRE_Int sw[6]; 105 HYPRE_Int hlower[3], hupper[3], period[3] = {0, 0, 0}; 106 PetscInt dim, dof, psw, Nx, Ny, Nz, nx, ny, nz, ilower[3], iupper[3], ssize, i; 107 DMBoundaryType px, py, pz; 108 DMDAStencilType st; 109 ISLocalToGlobalMapping ltog; 110 DM da; 111 112 PetscFunctionBegin; 113 PetscCall(MatGetDM(mat, (DM *)&da)); 114 ex->da = da; 115 PetscCall(PetscObjectReference((PetscObject)da)); 116 117 PetscCall(DMDAGetInfo(ex->da, &dim, &Nx, &Ny, &Nz, 0, 0, 0, &dof, &psw, &px, &py, &pz, &st)); 118 PetscCall(DMDAGetCorners(ex->da, &ilower[0], &ilower[1], &ilower[2], &iupper[0], &iupper[1], &iupper[2])); 119 120 /* when HYPRE_MIXEDINT is defined, sizeof(HYPRE_Int) == 32 */ 121 iupper[0] += ilower[0] - 1; 122 iupper[1] += ilower[1] - 1; 123 iupper[2] += ilower[2] - 1; 124 hlower[0] = (HYPRE_Int)ilower[0]; 125 hlower[1] = (HYPRE_Int)ilower[1]; 126 hlower[2] = (HYPRE_Int)ilower[2]; 127 hupper[0] = (HYPRE_Int)iupper[0]; 128 hupper[1] = (HYPRE_Int)iupper[1]; 129 hupper[2] = (HYPRE_Int)iupper[2]; 130 131 /* the hypre_Box is used to zero out the matrix entries in MatZeroValues() */ 132 ex->hbox.imin[0] = hlower[0]; 133 ex->hbox.imin[1] = hlower[1]; 134 ex->hbox.imin[2] = hlower[2]; 135 ex->hbox.imax[0] = hupper[0]; 136 ex->hbox.imax[1] = hupper[1]; 137 ex->hbox.imax[2] = hupper[2]; 138 139 /* create the hypre grid object and set its information */ 140 PetscCheck(dof <= 1, PetscObjectComm((PetscObject)da), PETSC_ERR_SUP, "Currently only support for scalar problems"); 141 if (px || py || pz) { 142 if (px == DM_BOUNDARY_PERIODIC) period[0] = (HYPRE_Int)Nx; 143 if (py == DM_BOUNDARY_PERIODIC) period[1] = (HYPRE_Int)Ny; 144 if (pz == DM_BOUNDARY_PERIODIC) period[2] = (HYPRE_Int)Nz; 145 } 146 PetscCallExternal(HYPRE_StructGridCreate, ex->hcomm, dim, &ex->hgrid); 147 PetscCallExternal(HYPRE_StructGridSetExtents, ex->hgrid, hlower, hupper); 148 PetscCallExternal(HYPRE_StructGridSetPeriodic, ex->hgrid, period); 149 PetscCallExternal(HYPRE_StructGridAssemble, ex->hgrid); 150 151 sw[5] = sw[4] = sw[3] = sw[2] = sw[1] = sw[0] = psw; 152 PetscCallExternal(HYPRE_StructGridSetNumGhost, ex->hgrid, sw); 153 154 /* create the hypre stencil object and set its information */ 155 PetscCheck(sw[0] <= 1, PetscObjectComm((PetscObject)da), PETSC_ERR_SUP, "Ask us to add support for wider stencils"); 156 PetscCheck(st != DMDA_STENCIL_BOX, PetscObjectComm((PetscObject)da), PETSC_ERR_SUP, "Ask us to add support for box stencils"); 157 if (dim == 1) { 158 HYPRE_Int offsets[3][1] = {{-1}, {0}, {1}}; 159 ssize = 3; 160 PetscCallExternal(HYPRE_StructStencilCreate, dim, ssize, &ex->hstencil); 161 for (i = 0; i < ssize; i++) PetscCallExternal(HYPRE_StructStencilSetElement, ex->hstencil, i, offsets[i]); 162 } else if (dim == 2) { 163 HYPRE_Int offsets[5][2] = { 164 {0, -1}, 165 {-1, 0 }, 166 {0, 0 }, 167 {1, 0 }, 168 {0, 1 } 169 }; 170 ssize = 5; 171 PetscCallExternal(HYPRE_StructStencilCreate, dim, ssize, &ex->hstencil); 172 for (i = 0; i < ssize; i++) PetscCallExternal(HYPRE_StructStencilSetElement, ex->hstencil, i, offsets[i]); 173 } else if (dim == 3) { 174 HYPRE_Int offsets[7][3] = { 175 {0, 0, -1}, 176 {0, -1, 0 }, 177 {-1, 0, 0 }, 178 {0, 0, 0 }, 179 {1, 0, 0 }, 180 {0, 1, 0 }, 181 {0, 0, 1 } 182 }; 183 ssize = 7; 184 PetscCallExternal(HYPRE_StructStencilCreate, dim, ssize, &ex->hstencil); 185 for (i = 0; i < ssize; i++) PetscCallExternal(HYPRE_StructStencilSetElement, ex->hstencil, i, offsets[i]); 186 } 187 188 /* create the HYPRE vector for rhs and solution */ 189 PetscCallExternal(HYPRE_StructVectorCreate, ex->hcomm, ex->hgrid, &ex->hb); 190 PetscCallExternal(HYPRE_StructVectorCreate, ex->hcomm, ex->hgrid, &ex->hx); 191 PetscCallExternal(HYPRE_StructVectorInitialize, ex->hb); 192 PetscCallExternal(HYPRE_StructVectorInitialize, ex->hx); 193 PetscCallExternal(HYPRE_StructVectorAssemble, ex->hb); 194 PetscCallExternal(HYPRE_StructVectorAssemble, ex->hx); 195 196 /* create the hypre matrix object and set its information */ 197 PetscCallExternal(HYPRE_StructMatrixCreate, ex->hcomm, ex->hgrid, ex->hstencil, &ex->hmat); 198 PetscCallExternal(HYPRE_StructGridDestroy, ex->hgrid); 199 PetscCallExternal(HYPRE_StructStencilDestroy, ex->hstencil); 200 if (ex->needsinitialization) { 201 PetscCallExternal(HYPRE_StructMatrixInitialize, ex->hmat); 202 ex->needsinitialization = PETSC_FALSE; 203 } 204 205 /* set the global and local sizes of the matrix */ 206 PetscCall(DMDAGetCorners(da, 0, 0, 0, &nx, &ny, &nz)); 207 PetscCall(MatSetSizes(mat, dof * nx * ny * nz, dof * nx * ny * nz, PETSC_DECIDE, PETSC_DECIDE)); 208 PetscCall(PetscLayoutSetBlockSize(mat->rmap, 1)); 209 PetscCall(PetscLayoutSetBlockSize(mat->cmap, 1)); 210 PetscCall(PetscLayoutSetUp(mat->rmap)); 211 PetscCall(PetscLayoutSetUp(mat->cmap)); 212 mat->preallocated = PETSC_TRUE; 213 214 if (dim == 3) { 215 mat->ops->setvalueslocal = MatSetValuesLocal_HYPREStruct_3d; 216 mat->ops->zerorowslocal = MatZeroRowsLocal_HYPREStruct_3d; 217 mat->ops->zeroentries = MatZeroEntries_HYPREStruct_3d; 218 219 /* PetscCall(MatZeroEntries_HYPREStruct_3d(mat)); */ 220 } else SETERRQ(PetscObjectComm((PetscObject)da), PETSC_ERR_SUP, "Only support for 3d DMDA currently"); 221 222 /* get values that will be used repeatedly in MatSetValuesLocal() and MatZeroRowsLocal() repeatedly */ 223 PetscCall(MatGetOwnershipRange(mat, &ex->rstart, NULL)); 224 PetscCall(DMGetLocalToGlobalMapping(ex->da, <og)); 225 PetscCall(ISLocalToGlobalMappingGetIndices(ltog, (const PetscInt **)&ex->gindices)); 226 PetscCall(DMDAGetGhostCorners(ex->da, 0, 0, 0, &ex->gnx, &ex->gnxgny, 0)); 227 ex->gnxgny *= ex->gnx; 228 PetscCall(DMDAGetCorners(ex->da, &ex->xs, &ex->ys, &ex->zs, &ex->nx, &ex->ny, 0)); 229 ex->nxny = ex->nx * ex->ny; 230 PetscFunctionReturn(PETSC_SUCCESS); 231 } 232 233 PetscErrorCode MatMult_HYPREStruct(Mat A, Vec x, Vec y) 234 { 235 const PetscScalar *xx; 236 PetscScalar *yy; 237 PetscInt ilower[3], iupper[3]; 238 HYPRE_Int hlower[3], hupper[3]; 239 Mat_HYPREStruct *mx = (Mat_HYPREStruct *)(A->data); 240 241 PetscFunctionBegin; 242 PetscCall(DMDAGetCorners(mx->da, &ilower[0], &ilower[1], &ilower[2], &iupper[0], &iupper[1], &iupper[2])); 243 /* when HYPRE_MIXEDINT is defined, sizeof(HYPRE_Int) == 32 */ 244 iupper[0] += ilower[0] - 1; 245 iupper[1] += ilower[1] - 1; 246 iupper[2] += ilower[2] - 1; 247 hlower[0] = (HYPRE_Int)ilower[0]; 248 hlower[1] = (HYPRE_Int)ilower[1]; 249 hlower[2] = (HYPRE_Int)ilower[2]; 250 hupper[0] = (HYPRE_Int)iupper[0]; 251 hupper[1] = (HYPRE_Int)iupper[1]; 252 hupper[2] = (HYPRE_Int)iupper[2]; 253 254 /* copy x values over to hypre */ 255 PetscCallExternal(HYPRE_StructVectorSetConstantValues, mx->hb, 0.0); 256 PetscCall(VecGetArrayRead(x, &xx)); 257 PetscCallExternal(HYPRE_StructVectorSetBoxValues, mx->hb, hlower, hupper, (HYPRE_Complex *)xx); 258 PetscCall(VecRestoreArrayRead(x, &xx)); 259 PetscCallExternal(HYPRE_StructVectorAssemble, mx->hb); 260 PetscCallExternal(HYPRE_StructMatrixMatvec, 1.0, mx->hmat, mx->hb, 0.0, mx->hx); 261 262 /* copy solution values back to PETSc */ 263 PetscCall(VecGetArray(y, &yy)); 264 PetscCallExternal(HYPRE_StructVectorGetBoxValues, mx->hx, hlower, hupper, (HYPRE_Complex *)yy); 265 PetscCall(VecRestoreArray(y, &yy)); 266 PetscFunctionReturn(PETSC_SUCCESS); 267 } 268 269 PetscErrorCode MatAssemblyEnd_HYPREStruct(Mat mat, MatAssemblyType mode) 270 { 271 Mat_HYPREStruct *ex = (Mat_HYPREStruct *)mat->data; 272 273 PetscFunctionBegin; 274 PetscCallExternal(HYPRE_StructMatrixAssemble, ex->hmat); 275 /* PetscCallExternal(HYPRE_StructMatrixPrint,"dummy",ex->hmat,0); */ 276 PetscFunctionReturn(PETSC_SUCCESS); 277 } 278 279 PetscErrorCode MatZeroEntries_HYPREStruct(Mat mat) 280 { 281 PetscFunctionBegin; 282 /* before the DMDA is set to the matrix the zero doesn't need to do anything */ 283 PetscFunctionReturn(PETSC_SUCCESS); 284 } 285 286 PetscErrorCode MatDestroy_HYPREStruct(Mat mat) 287 { 288 Mat_HYPREStruct *ex = (Mat_HYPREStruct *)mat->data; 289 290 PetscFunctionBegin; 291 PetscCallExternal(HYPRE_StructMatrixDestroy, ex->hmat); 292 PetscCallExternal(HYPRE_StructVectorDestroy, ex->hx); 293 PetscCallExternal(HYPRE_StructVectorDestroy, ex->hb); 294 PetscCall(PetscObjectDereference((PetscObject)ex->da)); 295 PetscCallMPI(MPI_Comm_free(&(ex->hcomm))); 296 PetscCall(PetscFree(ex)); 297 PetscFunctionReturn(PETSC_SUCCESS); 298 } 299 300 PETSC_EXTERN PetscErrorCode MatCreate_HYPREStruct(Mat B) 301 { 302 Mat_HYPREStruct *ex; 303 304 PetscFunctionBegin; 305 PetscCall(PetscNew(&ex)); 306 B->data = (void *)ex; 307 B->rmap->bs = 1; 308 B->assembled = PETSC_FALSE; 309 310 B->insertmode = NOT_SET_VALUES; 311 312 B->ops->assemblyend = MatAssemblyEnd_HYPREStruct; 313 B->ops->mult = MatMult_HYPREStruct; 314 B->ops->zeroentries = MatZeroEntries_HYPREStruct; 315 B->ops->destroy = MatDestroy_HYPREStruct; 316 B->ops->setup = MatSetUp_HYPREStruct; 317 318 ex->needsinitialization = PETSC_TRUE; 319 320 PetscCallMPI(MPI_Comm_dup(PetscObjectComm((PetscObject)B), &(ex->hcomm))); 321 PetscCall(PetscObjectChangeTypeName((PetscObject)B, MATHYPRESTRUCT)); 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 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 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 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 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 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 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 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 PetscFunctionReturn(PETSC_SUCCESS); 835 } 836