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