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