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 HYPRE_Int index[3], entries[9]; 30 PetscInt i, j, stencil, row; 31 HYPRE_Complex *values = (HYPRE_Complex *)y; 32 Mat_HYPREStruct *ex = (Mat_HYPREStruct *)mat->data; 33 34 PetscFunctionBegin; 35 PetscCheck(ncol <= 9, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "ncol %" PetscInt_FMT " > 9 too large", ncol); 36 for (i = 0; i < nrow; i++) { 37 for (j = 0; j < ncol; j++) { 38 stencil = icol[j] - irow[i]; 39 if (!stencil) { 40 entries[j] = 3; 41 } else if (stencil == -1) { 42 entries[j] = 2; 43 } else if (stencil == 1) { 44 entries[j] = 4; 45 } else if (stencil == -ex->gnx) { 46 entries[j] = 1; 47 } else if (stencil == ex->gnx) { 48 entries[j] = 5; 49 } else if (stencil == -ex->gnxgny) { 50 entries[j] = 0; 51 } else if (stencil == ex->gnxgny) { 52 entries[j] = 6; 53 } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Local row %" PetscInt_FMT " local column %" PetscInt_FMT " have bad stencil %" PetscInt_FMT, irow[i], icol[j], stencil); 54 } 55 row = ex->gindices[irow[i]] - ex->rstart; 56 index[0] = (HYPRE_Int)(ex->xs + (row % ex->nx)); 57 index[1] = (HYPRE_Int)(ex->ys + ((row / ex->nx) % ex->ny)); 58 index[2] = (HYPRE_Int)(ex->zs + (row / (ex->nxny))); 59 if (addv == ADD_VALUES) PetscCallExternal(HYPRE_StructMatrixAddToValues, ex->hmat, index, ncol, entries, values); 60 else PetscCallExternal(HYPRE_StructMatrixSetValues, ex->hmat, index, ncol, entries, values); 61 values += ncol; 62 } 63 PetscFunctionReturn(0); 64 } 65 66 PetscErrorCode MatZeroRowsLocal_HYPREStruct_3d(Mat mat, PetscInt nrow, const PetscInt irow[], PetscScalar d, Vec x, Vec b) { 67 HYPRE_Int index[3], entries[7] = {0, 1, 2, 3, 4, 5, 6}; 68 PetscInt row, i; 69 HYPRE_Complex values[7]; 70 Mat_HYPREStruct *ex = (Mat_HYPREStruct *)mat->data; 71 72 PetscFunctionBegin; 73 PetscCheck(!x || !b, PetscObjectComm((PetscObject)mat), PETSC_ERR_SUP, "No support"); 74 PetscCall(PetscArrayzero(values, 7)); 75 PetscCall(PetscHYPREScalarCast(d, &values[3])); 76 for (i = 0; i < nrow; i++) { 77 row = ex->gindices[irow[i]] - ex->rstart; 78 index[0] = (HYPRE_Int)(ex->xs + (row % ex->nx)); 79 index[1] = (HYPRE_Int)(ex->ys + ((row / ex->nx) % ex->ny)); 80 index[2] = (HYPRE_Int)(ex->zs + (row / (ex->nxny))); 81 PetscCallExternal(HYPRE_StructMatrixSetValues, ex->hmat, index, 7, entries, values); 82 } 83 PetscCallExternal(HYPRE_StructMatrixAssemble, ex->hmat); 84 PetscFunctionReturn(0); 85 } 86 87 PetscErrorCode MatZeroEntries_HYPREStruct_3d(Mat mat) { 88 HYPRE_Int indices[7] = {0, 1, 2, 3, 4, 5, 6}; 89 Mat_HYPREStruct *ex = (Mat_HYPREStruct *)mat->data; 90 91 PetscFunctionBegin; 92 /* hypre has no public interface to do this */ 93 PetscCallExternal(hypre_StructMatrixClearBoxValues, ex->hmat, &ex->hbox, 7, indices, 0, 1); 94 PetscCallExternal(HYPRE_StructMatrixAssemble, ex->hmat); 95 PetscFunctionReturn(0); 96 } 97 98 static PetscErrorCode MatSetUp_HYPREStruct(Mat mat) { 99 Mat_HYPREStruct *ex = (Mat_HYPREStruct *)mat->data; 100 HYPRE_Int sw[6]; 101 HYPRE_Int hlower[3], hupper[3], period[3] = {0, 0, 0}; 102 PetscInt dim, dof, psw, Nx, Ny, Nz, nx, ny, nz, ilower[3], iupper[3], ssize, i; 103 DMBoundaryType px, py, pz; 104 DMDAStencilType st; 105 ISLocalToGlobalMapping ltog; 106 DM da; 107 108 PetscFunctionBegin; 109 PetscCall(MatGetDM(mat, (DM *)&da)); 110 ex->da = da; 111 PetscCall(PetscObjectReference((PetscObject)da)); 112 113 PetscCall(DMDAGetInfo(ex->da, &dim, &Nx, &Ny, &Nz, 0, 0, 0, &dof, &psw, &px, &py, &pz, &st)); 114 PetscCall(DMDAGetCorners(ex->da, &ilower[0], &ilower[1], &ilower[2], &iupper[0], &iupper[1], &iupper[2])); 115 116 /* when HYPRE_MIXEDINT is defined, sizeof(HYPRE_Int) == 32 */ 117 iupper[0] += ilower[0] - 1; 118 iupper[1] += ilower[1] - 1; 119 iupper[2] += ilower[2] - 1; 120 hlower[0] = (HYPRE_Int)ilower[0]; 121 hlower[1] = (HYPRE_Int)ilower[1]; 122 hlower[2] = (HYPRE_Int)ilower[2]; 123 hupper[0] = (HYPRE_Int)iupper[0]; 124 hupper[1] = (HYPRE_Int)iupper[1]; 125 hupper[2] = (HYPRE_Int)iupper[2]; 126 127 /* the hypre_Box is used to zero out the matrix entries in MatZeroValues() */ 128 ex->hbox.imin[0] = hlower[0]; 129 ex->hbox.imin[1] = hlower[1]; 130 ex->hbox.imin[2] = hlower[2]; 131 ex->hbox.imax[0] = hupper[0]; 132 ex->hbox.imax[1] = hupper[1]; 133 ex->hbox.imax[2] = hupper[2]; 134 135 /* create the hypre grid object and set its information */ 136 PetscCheck(dof <= 1, PetscObjectComm((PetscObject)da), PETSC_ERR_SUP, "Currently only support for scalar problems"); 137 if (px || py || pz) { 138 if (px == DM_BOUNDARY_PERIODIC) period[0] = (HYPRE_Int)Nx; 139 if (py == DM_BOUNDARY_PERIODIC) period[1] = (HYPRE_Int)Ny; 140 if (pz == DM_BOUNDARY_PERIODIC) period[2] = (HYPRE_Int)Nz; 141 } 142 PetscCallExternal(HYPRE_StructGridCreate, ex->hcomm, dim, &ex->hgrid); 143 PetscCallExternal(HYPRE_StructGridSetExtents, ex->hgrid, hlower, hupper); 144 PetscCallExternal(HYPRE_StructGridSetPeriodic, ex->hgrid, period); 145 PetscCallExternal(HYPRE_StructGridAssemble, ex->hgrid); 146 147 sw[5] = sw[4] = sw[3] = sw[2] = sw[1] = sw[0] = psw; 148 PetscCallExternal(HYPRE_StructGridSetNumGhost, ex->hgrid, sw); 149 150 /* create the hypre stencil object and set its information */ 151 PetscCheck(sw[0] <= 1, PetscObjectComm((PetscObject)da), PETSC_ERR_SUP, "Ask us to add support for wider stencils"); 152 PetscCheck(st != DMDA_STENCIL_BOX, PetscObjectComm((PetscObject)da), PETSC_ERR_SUP, "Ask us to add support for box stencils"); 153 if (dim == 1) { 154 HYPRE_Int offsets[3][1] = {{-1}, {0}, {1}}; 155 ssize = 3; 156 PetscCallExternal(HYPRE_StructStencilCreate, dim, ssize, &ex->hstencil); 157 for (i = 0; i < ssize; i++) PetscCallExternal(HYPRE_StructStencilSetElement, ex->hstencil, i, offsets[i]); 158 } else if (dim == 2) { 159 HYPRE_Int offsets[5][2] = { 160 {0, -1}, 161 {-1, 0 }, 162 {0, 0 }, 163 {1, 0 }, 164 {0, 1 } 165 }; 166 ssize = 5; 167 PetscCallExternal(HYPRE_StructStencilCreate, dim, ssize, &ex->hstencil); 168 for (i = 0; i < ssize; i++) PetscCallExternal(HYPRE_StructStencilSetElement, ex->hstencil, i, offsets[i]); 169 } else if (dim == 3) { 170 HYPRE_Int offsets[7][3] = { 171 {0, 0, -1}, 172 {0, -1, 0 }, 173 {-1, 0, 0 }, 174 {0, 0, 0 }, 175 {1, 0, 0 }, 176 {0, 1, 0 }, 177 {0, 0, 1 } 178 }; 179 ssize = 7; 180 PetscCallExternal(HYPRE_StructStencilCreate, dim, ssize, &ex->hstencil); 181 for (i = 0; i < ssize; i++) PetscCallExternal(HYPRE_StructStencilSetElement, ex->hstencil, i, offsets[i]); 182 } 183 184 /* create the HYPRE vector for rhs and solution */ 185 PetscCallExternal(HYPRE_StructVectorCreate, ex->hcomm, ex->hgrid, &ex->hb); 186 PetscCallExternal(HYPRE_StructVectorCreate, ex->hcomm, ex->hgrid, &ex->hx); 187 PetscCallExternal(HYPRE_StructVectorInitialize, ex->hb); 188 PetscCallExternal(HYPRE_StructVectorInitialize, ex->hx); 189 PetscCallExternal(HYPRE_StructVectorAssemble, ex->hb); 190 PetscCallExternal(HYPRE_StructVectorAssemble, ex->hx); 191 192 /* create the hypre matrix object and set its information */ 193 PetscCallExternal(HYPRE_StructMatrixCreate, ex->hcomm, ex->hgrid, ex->hstencil, &ex->hmat); 194 PetscCallExternal(HYPRE_StructGridDestroy, ex->hgrid); 195 PetscCallExternal(HYPRE_StructStencilDestroy, ex->hstencil); 196 if (ex->needsinitialization) { 197 PetscCallExternal(HYPRE_StructMatrixInitialize, ex->hmat); 198 ex->needsinitialization = PETSC_FALSE; 199 } 200 201 /* set the global and local sizes of the matrix */ 202 PetscCall(DMDAGetCorners(da, 0, 0, 0, &nx, &ny, &nz)); 203 PetscCall(MatSetSizes(mat, dof * nx * ny * nz, dof * nx * ny * nz, PETSC_DECIDE, PETSC_DECIDE)); 204 PetscCall(PetscLayoutSetBlockSize(mat->rmap, 1)); 205 PetscCall(PetscLayoutSetBlockSize(mat->cmap, 1)); 206 PetscCall(PetscLayoutSetUp(mat->rmap)); 207 PetscCall(PetscLayoutSetUp(mat->cmap)); 208 mat->preallocated = PETSC_TRUE; 209 210 if (dim == 3) { 211 mat->ops->setvalueslocal = MatSetValuesLocal_HYPREStruct_3d; 212 mat->ops->zerorowslocal = MatZeroRowsLocal_HYPREStruct_3d; 213 mat->ops->zeroentries = MatZeroEntries_HYPREStruct_3d; 214 215 /* PetscCall(MatZeroEntries_HYPREStruct_3d(mat)); */ 216 } else SETERRQ(PetscObjectComm((PetscObject)da), PETSC_ERR_SUP, "Only support for 3d DMDA currently"); 217 218 /* get values that will be used repeatedly in MatSetValuesLocal() and MatZeroRowsLocal() repeatedly */ 219 PetscCall(MatGetOwnershipRange(mat, &ex->rstart, NULL)); 220 PetscCall(DMGetLocalToGlobalMapping(ex->da, <og)); 221 PetscCall(ISLocalToGlobalMappingGetIndices(ltog, (const PetscInt **)&ex->gindices)); 222 PetscCall(DMDAGetGhostCorners(ex->da, 0, 0, 0, &ex->gnx, &ex->gnxgny, 0)); 223 ex->gnxgny *= ex->gnx; 224 PetscCall(DMDAGetCorners(ex->da, &ex->xs, &ex->ys, &ex->zs, &ex->nx, &ex->ny, 0)); 225 ex->nxny = ex->nx * ex->ny; 226 PetscFunctionReturn(0); 227 } 228 229 PetscErrorCode MatMult_HYPREStruct(Mat A, Vec x, Vec y) { 230 const PetscScalar *xx; 231 PetscScalar *yy; 232 PetscInt ilower[3], iupper[3]; 233 HYPRE_Int hlower[3], hupper[3]; 234 Mat_HYPREStruct *mx = (Mat_HYPREStruct *)(A->data); 235 236 PetscFunctionBegin; 237 PetscCall(DMDAGetCorners(mx->da, &ilower[0], &ilower[1], &ilower[2], &iupper[0], &iupper[1], &iupper[2])); 238 /* when HYPRE_MIXEDINT is defined, sizeof(HYPRE_Int) == 32 */ 239 iupper[0] += ilower[0] - 1; 240 iupper[1] += ilower[1] - 1; 241 iupper[2] += ilower[2] - 1; 242 hlower[0] = (HYPRE_Int)ilower[0]; 243 hlower[1] = (HYPRE_Int)ilower[1]; 244 hlower[2] = (HYPRE_Int)ilower[2]; 245 hupper[0] = (HYPRE_Int)iupper[0]; 246 hupper[1] = (HYPRE_Int)iupper[1]; 247 hupper[2] = (HYPRE_Int)iupper[2]; 248 249 /* copy x values over to hypre */ 250 PetscCallExternal(HYPRE_StructVectorSetConstantValues, mx->hb, 0.0); 251 PetscCall(VecGetArrayRead(x, &xx)); 252 PetscCallExternal(HYPRE_StructVectorSetBoxValues, mx->hb, hlower, hupper, (HYPRE_Complex *)xx); 253 PetscCall(VecRestoreArrayRead(x, &xx)); 254 PetscCallExternal(HYPRE_StructVectorAssemble, mx->hb); 255 PetscCallExternal(HYPRE_StructMatrixMatvec, 1.0, mx->hmat, mx->hb, 0.0, mx->hx); 256 257 /* copy solution values back to PETSc */ 258 PetscCall(VecGetArray(y, &yy)); 259 PetscCallExternal(HYPRE_StructVectorGetBoxValues, mx->hx, hlower, hupper, (HYPRE_Complex *)yy); 260 PetscCall(VecRestoreArray(y, &yy)); 261 PetscFunctionReturn(0); 262 } 263 264 PetscErrorCode MatAssemblyEnd_HYPREStruct(Mat mat, MatAssemblyType mode) { 265 Mat_HYPREStruct *ex = (Mat_HYPREStruct *)mat->data; 266 267 PetscFunctionBegin; 268 PetscCallExternal(HYPRE_StructMatrixAssemble, ex->hmat); 269 /* PetscCallExternal(HYPRE_StructMatrixPrint,"dummy",ex->hmat,0); */ 270 PetscFunctionReturn(0); 271 } 272 273 PetscErrorCode MatZeroEntries_HYPREStruct(Mat mat) { 274 PetscFunctionBegin; 275 /* before the DMDA is set to the matrix the zero doesn't need to do anything */ 276 PetscFunctionReturn(0); 277 } 278 279 PetscErrorCode MatDestroy_HYPREStruct(Mat mat) { 280 Mat_HYPREStruct *ex = (Mat_HYPREStruct *)mat->data; 281 282 PetscFunctionBegin; 283 PetscCallExternal(HYPRE_StructMatrixDestroy, ex->hmat); 284 PetscCallExternal(HYPRE_StructVectorDestroy, ex->hx); 285 PetscCallExternal(HYPRE_StructVectorDestroy, ex->hb); 286 PetscCall(PetscObjectDereference((PetscObject)ex->da)); 287 PetscCallMPI(MPI_Comm_free(&(ex->hcomm))); 288 PetscCall(PetscFree(ex)); 289 PetscFunctionReturn(0); 290 } 291 292 PETSC_EXTERN PetscErrorCode MatCreate_HYPREStruct(Mat B) { 293 Mat_HYPREStruct *ex; 294 295 PetscFunctionBegin; 296 PetscCall(PetscNew(&ex)); 297 B->data = (void *)ex; 298 B->rmap->bs = 1; 299 B->assembled = PETSC_FALSE; 300 301 B->insertmode = NOT_SET_VALUES; 302 303 B->ops->assemblyend = MatAssemblyEnd_HYPREStruct; 304 B->ops->mult = MatMult_HYPREStruct; 305 B->ops->zeroentries = MatZeroEntries_HYPREStruct; 306 B->ops->destroy = MatDestroy_HYPREStruct; 307 B->ops->setup = MatSetUp_HYPREStruct; 308 309 ex->needsinitialization = PETSC_TRUE; 310 311 PetscCallMPI(MPI_Comm_dup(PetscObjectComm((PetscObject)B), &(ex->hcomm))); 312 PetscCall(PetscObjectChangeTypeName((PetscObject)B, MATHYPRESTRUCT)); 313 PetscFunctionReturn(0); 314 } 315 316 /*MC 317 MATHYPRESSTRUCT - MATHYPRESSTRUCT = "hypresstruct" - A matrix type to be used for parallel sparse matrices 318 based on the hypre HYPRE_SStructMatrix. 319 320 Level: intermediate 321 322 Notes: 323 Unlike hypre's general semi-struct object consisting of a collection of structured-grid objects and unstructured 324 grid objects, we restrict the semi-struct objects to consist of only structured-grid components. 325 326 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 327 be defined by a DMDA. 328 329 The matrix needs a DMDA associated with it by either a call to MatSetDM() or if the matrix is obtained from DMCreateMatrix() 330 331 M*/ 332 333 PetscErrorCode MatSetValuesLocal_HYPRESStruct_3d(Mat mat, PetscInt nrow, const PetscInt irow[], PetscInt ncol, const PetscInt icol[], const PetscScalar y[], InsertMode addv) { 334 HYPRE_Int index[3], *entries; 335 PetscInt i, j, stencil; 336 HYPRE_Complex *values = (HYPRE_Complex *)y; 337 Mat_HYPRESStruct *ex = (Mat_HYPRESStruct *)mat->data; 338 339 PetscInt part = 0; /* Petsc sstruct interface only allows 1 part */ 340 PetscInt ordering; 341 PetscInt grid_rank, to_grid_rank; 342 PetscInt var_type, to_var_type; 343 PetscInt to_var_entry = 0; 344 PetscInt nvars = ex->nvars; 345 PetscInt row; 346 347 PetscFunctionBegin; 348 PetscCall(PetscMalloc1(7 * nvars, &entries)); 349 ordering = ex->dofs_order; /* ordering= 0 nodal ordering 350 1 variable ordering */ 351 /* stencil entries are ordered by variables: var0_stencil0, var0_stencil1, ..., var0_stencil6, var1_stencil0, var1_stencil1, ... */ 352 if (!ordering) { 353 for (i = 0; i < nrow; i++) { 354 grid_rank = irow[i] / nvars; 355 var_type = (irow[i] % nvars); 356 357 for (j = 0; j < ncol; j++) { 358 to_grid_rank = icol[j] / nvars; 359 to_var_type = (icol[j] % nvars); 360 361 to_var_entry = to_var_entry * 7; 362 entries[j] = to_var_entry; 363 364 stencil = to_grid_rank - grid_rank; 365 if (!stencil) { 366 entries[j] += 3; 367 } else if (stencil == -1) { 368 entries[j] += 2; 369 } else if (stencil == 1) { 370 entries[j] += 4; 371 } else if (stencil == -ex->gnx) { 372 entries[j] += 1; 373 } else if (stencil == ex->gnx) { 374 entries[j] += 5; 375 } else if (stencil == -ex->gnxgny) { 376 entries[j] += 0; 377 } else if (stencil == ex->gnxgny) { 378 entries[j] += 6; 379 } 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); 380 } 381 382 row = ex->gindices[grid_rank] - ex->rstart; 383 index[0] = (HYPRE_Int)(ex->xs + (row % ex->nx)); 384 index[1] = (HYPRE_Int)(ex->ys + ((row / ex->nx) % ex->ny)); 385 index[2] = (HYPRE_Int)(ex->zs + (row / (ex->nxny))); 386 387 if (addv == ADD_VALUES) PetscCallExternal(HYPRE_SStructMatrixAddToValues, ex->ss_mat, part, index, var_type, ncol, entries, values); 388 else PetscCallExternal(HYPRE_SStructMatrixSetValues, ex->ss_mat, part, index, var_type, ncol, entries, values); 389 values += ncol; 390 } 391 } else { 392 for (i = 0; i < nrow; i++) { 393 var_type = irow[i] / (ex->gnxgnygnz); 394 grid_rank = irow[i] - var_type * (ex->gnxgnygnz); 395 396 for (j = 0; j < ncol; j++) { 397 to_var_type = icol[j] / (ex->gnxgnygnz); 398 to_grid_rank = icol[j] - to_var_type * (ex->gnxgnygnz); 399 400 to_var_entry = to_var_entry * 7; 401 entries[j] = to_var_entry; 402 403 stencil = to_grid_rank - grid_rank; 404 if (!stencil) { 405 entries[j] += 3; 406 } else if (stencil == -1) { 407 entries[j] += 2; 408 } else if (stencil == 1) { 409 entries[j] += 4; 410 } else if (stencil == -ex->gnx) { 411 entries[j] += 1; 412 } else if (stencil == ex->gnx) { 413 entries[j] += 5; 414 } else if (stencil == -ex->gnxgny) { 415 entries[j] += 0; 416 } else if (stencil == ex->gnxgny) { 417 entries[j] += 6; 418 } 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); 419 } 420 421 row = ex->gindices[grid_rank] - ex->rstart; 422 index[0] = (HYPRE_Int)(ex->xs + (row % ex->nx)); 423 index[1] = (HYPRE_Int)(ex->ys + ((row / ex->nx) % ex->ny)); 424 index[2] = (HYPRE_Int)(ex->zs + (row / (ex->nxny))); 425 426 if (addv == ADD_VALUES) PetscCallExternal(HYPRE_SStructMatrixAddToValues, ex->ss_mat, part, index, var_type, ncol, entries, values); 427 else PetscCallExternal(HYPRE_SStructMatrixSetValues, ex->ss_mat, part, index, var_type, ncol, entries, values); 428 values += ncol; 429 } 430 } 431 PetscCall(PetscFree(entries)); 432 PetscFunctionReturn(0); 433 } 434 435 PetscErrorCode MatZeroRowsLocal_HYPRESStruct_3d(Mat mat, PetscInt nrow, const PetscInt irow[], PetscScalar d, Vec x, Vec b) { 436 HYPRE_Int index[3], *entries; 437 PetscInt i; 438 HYPRE_Complex **values; 439 Mat_HYPRESStruct *ex = (Mat_HYPRESStruct *)mat->data; 440 441 PetscInt part = 0; /* Petsc sstruct interface only allows 1 part */ 442 PetscInt ordering = ex->dofs_order; 443 PetscInt grid_rank; 444 PetscInt var_type; 445 PetscInt nvars = ex->nvars; 446 PetscInt row; 447 448 PetscFunctionBegin; 449 PetscCheck(!x || !b, PetscObjectComm((PetscObject)mat), PETSC_ERR_SUP, "No support"); 450 PetscCall(PetscMalloc1(7 * nvars, &entries)); 451 452 PetscCall(PetscMalloc1(nvars, &values)); 453 PetscCall(PetscMalloc1(7 * nvars * nvars, &values[0])); 454 for (i = 1; i < nvars; i++) values[i] = values[i - 1] + nvars * 7; 455 456 for (i = 0; i < nvars; i++) { 457 PetscCall(PetscArrayzero(values[i], nvars * 7 * sizeof(HYPRE_Complex))); 458 PetscCall(PetscHYPREScalarCast(d, values[i] + 3)); 459 } 460 461 for (i = 0; i < nvars * 7; i++) entries[i] = i; 462 463 if (!ordering) { 464 for (i = 0; i < nrow; i++) { 465 grid_rank = irow[i] / nvars; 466 var_type = (irow[i] % nvars); 467 468 row = ex->gindices[grid_rank] - ex->rstart; 469 index[0] = (HYPRE_Int)(ex->xs + (row % ex->nx)); 470 index[1] = (HYPRE_Int)(ex->ys + ((row / ex->nx) % ex->ny)); 471 index[2] = (HYPRE_Int)(ex->zs + (row / (ex->nxny))); 472 PetscCallExternal(HYPRE_SStructMatrixSetValues, ex->ss_mat, part, index, var_type, 7 * nvars, entries, values[var_type]); 473 } 474 } else { 475 for (i = 0; i < nrow; i++) { 476 var_type = irow[i] / (ex->gnxgnygnz); 477 grid_rank = irow[i] - var_type * (ex->gnxgnygnz); 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 } 486 PetscCallExternal(HYPRE_SStructMatrixAssemble, ex->ss_mat); 487 PetscCall(PetscFree(values[0])); 488 PetscCall(PetscFree(values)); 489 PetscCall(PetscFree(entries)); 490 PetscFunctionReturn(0); 491 } 492 493 PetscErrorCode MatZeroEntries_HYPRESStruct_3d(Mat mat) { 494 Mat_HYPRESStruct *ex = (Mat_HYPRESStruct *)mat->data; 495 PetscInt nvars = ex->nvars; 496 PetscInt size; 497 PetscInt part = 0; /* only one part */ 498 499 PetscFunctionBegin; 500 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); 501 { 502 HYPRE_Int i, *entries, iupper[3], ilower[3]; 503 HYPRE_Complex *values; 504 505 for (i = 0; i < 3; i++) { 506 ilower[i] = ex->hbox.imin[i]; 507 iupper[i] = ex->hbox.imax[i]; 508 } 509 510 PetscCall(PetscMalloc2(nvars * 7, &entries, nvars * 7 * size, &values)); 511 for (i = 0; i < nvars * 7; i++) entries[i] = i; 512 PetscCall(PetscArrayzero(values, nvars * 7 * size)); 513 514 for (i = 0; i < nvars; i++) PetscCallExternal(HYPRE_SStructMatrixSetBoxValues, ex->ss_mat, part, ilower, iupper, i, nvars * 7, entries, values); 515 PetscCall(PetscFree2(entries, values)); 516 } 517 PetscCallExternal(HYPRE_SStructMatrixAssemble, ex->ss_mat); 518 PetscFunctionReturn(0); 519 } 520 521 static PetscErrorCode MatSetUp_HYPRESStruct(Mat mat) { 522 Mat_HYPRESStruct *ex = (Mat_HYPRESStruct *)mat->data; 523 PetscInt dim, dof, sw[3], nx, ny, nz; 524 PetscInt ilower[3], iupper[3], ssize, i; 525 DMBoundaryType px, py, pz; 526 DMDAStencilType st; 527 PetscInt nparts = 1; /* assuming only one part */ 528 PetscInt part = 0; 529 ISLocalToGlobalMapping ltog; 530 DM da; 531 532 PetscFunctionBegin; 533 PetscCall(MatGetDM(mat, (DM *)&da)); 534 ex->da = da; 535 PetscCall(PetscObjectReference((PetscObject)da)); 536 537 PetscCall(DMDAGetInfo(ex->da, &dim, 0, 0, 0, 0, 0, 0, &dof, &sw[0], &px, &py, &pz, &st)); 538 PetscCall(DMDAGetCorners(ex->da, &ilower[0], &ilower[1], &ilower[2], &iupper[0], &iupper[1], &iupper[2])); 539 iupper[0] += ilower[0] - 1; 540 iupper[1] += ilower[1] - 1; 541 iupper[2] += ilower[2] - 1; 542 /* the hypre_Box is used to zero out the matrix entries in MatZeroValues() */ 543 ex->hbox.imin[0] = ilower[0]; 544 ex->hbox.imin[1] = ilower[1]; 545 ex->hbox.imin[2] = ilower[2]; 546 ex->hbox.imax[0] = iupper[0]; 547 ex->hbox.imax[1] = iupper[1]; 548 ex->hbox.imax[2] = iupper[2]; 549 550 ex->dofs_order = 0; 551 552 /* assuming that the same number of dofs on each gridpoint. Also assume all cell-centred based */ 553 ex->nvars = dof; 554 555 /* create the hypre grid object and set its information */ 556 PetscCheck(!px && !py && !pz, PetscObjectComm((PetscObject)da), PETSC_ERR_SUP, "Ask us to add periodic support by calling HYPRE_SStructGridSetPeriodic()"); 557 PetscCallExternal(HYPRE_SStructGridCreate, ex->hcomm, dim, nparts, &ex->ss_grid); 558 PetscCallExternal(HYPRE_SStructGridSetExtents, ex->ss_grid, part, ex->hbox.imin, ex->hbox.imax); 559 { 560 HYPRE_SStructVariable *vartypes; 561 PetscCall(PetscMalloc1(ex->nvars, &vartypes)); 562 for (i = 0; i < ex->nvars; i++) vartypes[i] = HYPRE_SSTRUCT_VARIABLE_CELL; 563 PetscCallExternal(HYPRE_SStructGridSetVariables, ex->ss_grid, part, ex->nvars, vartypes); 564 PetscCall(PetscFree(vartypes)); 565 } 566 PetscCallExternal(HYPRE_SStructGridAssemble, ex->ss_grid); 567 568 sw[1] = sw[0]; 569 sw[2] = sw[1]; 570 /* PetscCallExternal(HYPRE_SStructGridSetNumGhost,ex->ss_grid,sw); */ 571 572 /* create the hypre stencil object and set its information */ 573 PetscCheck(sw[0] <= 1, PetscObjectComm((PetscObject)da), PETSC_ERR_SUP, "Ask us to add support for wider stencils"); 574 PetscCheck(st != DMDA_STENCIL_BOX, PetscObjectComm((PetscObject)da), PETSC_ERR_SUP, "Ask us to add support for box stencils"); 575 576 if (dim == 1) { 577 HYPRE_Int offsets[3][1] = {{-1}, {0}, {1}}; 578 PetscInt j, cnt; 579 580 ssize = 3 * (ex->nvars); 581 PetscCallExternal(HYPRE_SStructStencilCreate, dim, ssize, &ex->ss_stencil); 582 cnt = 0; 583 for (i = 0; i < (ex->nvars); i++) { 584 for (j = 0; j < 3; j++) { 585 PetscCallExternal(HYPRE_SStructStencilSetEntry, ex->ss_stencil, cnt, offsets[j], i); 586 cnt++; 587 } 588 } 589 } else if (dim == 2) { 590 HYPRE_Int offsets[5][2] = { 591 {0, -1}, 592 {-1, 0 }, 593 {0, 0 }, 594 {1, 0 }, 595 {0, 1 } 596 }; 597 PetscInt j, cnt; 598 599 ssize = 5 * (ex->nvars); 600 PetscCallExternal(HYPRE_SStructStencilCreate, dim, ssize, &ex->ss_stencil); 601 cnt = 0; 602 for (i = 0; i < (ex->nvars); i++) { 603 for (j = 0; j < 5; j++) { 604 PetscCallExternal(HYPRE_SStructStencilSetEntry, ex->ss_stencil, cnt, offsets[j], i); 605 cnt++; 606 } 607 } 608 } else if (dim == 3) { 609 HYPRE_Int offsets[7][3] = { 610 {0, 0, -1}, 611 {0, -1, 0 }, 612 {-1, 0, 0 }, 613 {0, 0, 0 }, 614 {1, 0, 0 }, 615 {0, 1, 0 }, 616 {0, 0, 1 } 617 }; 618 PetscInt j, cnt; 619 620 ssize = 7 * (ex->nvars); 621 PetscCallExternal(HYPRE_SStructStencilCreate, dim, ssize, &ex->ss_stencil); 622 cnt = 0; 623 for (i = 0; i < (ex->nvars); i++) { 624 for (j = 0; j < 7; j++) { 625 PetscCallExternal(HYPRE_SStructStencilSetEntry, ex->ss_stencil, cnt, offsets[j], i); 626 cnt++; 627 } 628 } 629 } 630 631 /* create the HYPRE graph */ 632 PetscCallExternal(HYPRE_SStructGraphCreate, ex->hcomm, ex->ss_grid, &(ex->ss_graph)); 633 634 /* set the stencil graph. Note that each variable has the same graph. This means that each 635 variable couples to all the other variable and with the same stencil pattern. */ 636 for (i = 0; i < (ex->nvars); i++) PetscCallExternal(HYPRE_SStructGraphSetStencil, ex->ss_graph, part, i, ex->ss_stencil); 637 PetscCallExternal(HYPRE_SStructGraphAssemble, ex->ss_graph); 638 639 /* create the HYPRE sstruct vectors for rhs and solution */ 640 PetscCallExternal(HYPRE_SStructVectorCreate, ex->hcomm, ex->ss_grid, &ex->ss_b); 641 PetscCallExternal(HYPRE_SStructVectorCreate, ex->hcomm, ex->ss_grid, &ex->ss_x); 642 PetscCallExternal(HYPRE_SStructVectorInitialize, ex->ss_b); 643 PetscCallExternal(HYPRE_SStructVectorInitialize, ex->ss_x); 644 PetscCallExternal(HYPRE_SStructVectorAssemble, ex->ss_b); 645 PetscCallExternal(HYPRE_SStructVectorAssemble, ex->ss_x); 646 647 /* create the hypre matrix object and set its information */ 648 PetscCallExternal(HYPRE_SStructMatrixCreate, ex->hcomm, ex->ss_graph, &ex->ss_mat); 649 PetscCallExternal(HYPRE_SStructGridDestroy, ex->ss_grid); 650 PetscCallExternal(HYPRE_SStructStencilDestroy, ex->ss_stencil); 651 if (ex->needsinitialization) { 652 PetscCallExternal(HYPRE_SStructMatrixInitialize, ex->ss_mat); 653 ex->needsinitialization = PETSC_FALSE; 654 } 655 656 /* set the global and local sizes of the matrix */ 657 PetscCall(DMDAGetCorners(da, 0, 0, 0, &nx, &ny, &nz)); 658 PetscCall(MatSetSizes(mat, dof * nx * ny * nz, dof * nx * ny * nz, PETSC_DECIDE, PETSC_DECIDE)); 659 PetscCall(PetscLayoutSetBlockSize(mat->rmap, dof)); 660 PetscCall(PetscLayoutSetBlockSize(mat->cmap, dof)); 661 PetscCall(PetscLayoutSetUp(mat->rmap)); 662 PetscCall(PetscLayoutSetUp(mat->cmap)); 663 mat->preallocated = PETSC_TRUE; 664 665 if (dim == 3) { 666 mat->ops->setvalueslocal = MatSetValuesLocal_HYPRESStruct_3d; 667 mat->ops->zerorowslocal = MatZeroRowsLocal_HYPRESStruct_3d; 668 mat->ops->zeroentries = MatZeroEntries_HYPRESStruct_3d; 669 670 /* PetscCall(MatZeroEntries_HYPRESStruct_3d(mat)); */ 671 } else SETERRQ(PetscObjectComm((PetscObject)da), PETSC_ERR_SUP, "Only support for 3d DMDA currently"); 672 673 /* get values that will be used repeatedly in MatSetValuesLocal() and MatZeroRowsLocal() repeatedly */ 674 PetscCall(MatGetOwnershipRange(mat, &ex->rstart, NULL)); 675 PetscCall(DMGetLocalToGlobalMapping(ex->da, <og)); 676 PetscCall(ISLocalToGlobalMappingGetIndices(ltog, (const PetscInt **)&ex->gindices)); 677 PetscCall(DMDAGetGhostCorners(ex->da, 0, 0, 0, &ex->gnx, &ex->gnxgny, &ex->gnxgnygnz)); 678 679 ex->gnxgny *= ex->gnx; 680 ex->gnxgnygnz *= ex->gnxgny; 681 682 PetscCall(DMDAGetCorners(ex->da, &ex->xs, &ex->ys, &ex->zs, &ex->nx, &ex->ny, &ex->nz)); 683 684 ex->nxny = ex->nx * ex->ny; 685 ex->nxnynz = ex->nz * ex->nxny; 686 PetscFunctionReturn(0); 687 } 688 689 PetscErrorCode MatMult_HYPRESStruct(Mat A, Vec x, Vec y) { 690 const PetscScalar *xx; 691 PetscScalar *yy; 692 HYPRE_Int hlower[3], hupper[3]; 693 PetscInt ilower[3], iupper[3]; 694 Mat_HYPRESStruct *mx = (Mat_HYPRESStruct *)(A->data); 695 PetscInt ordering = mx->dofs_order; 696 PetscInt nvars = mx->nvars; 697 PetscInt part = 0; 698 PetscInt size; 699 PetscInt i; 700 701 PetscFunctionBegin; 702 PetscCall(DMDAGetCorners(mx->da, &ilower[0], &ilower[1], &ilower[2], &iupper[0], &iupper[1], &iupper[2])); 703 704 /* when HYPRE_MIXEDINT is defined, sizeof(HYPRE_Int) == 32 */ 705 iupper[0] += ilower[0] - 1; 706 iupper[1] += ilower[1] - 1; 707 iupper[2] += ilower[2] - 1; 708 hlower[0] = (HYPRE_Int)ilower[0]; 709 hlower[1] = (HYPRE_Int)ilower[1]; 710 hlower[2] = (HYPRE_Int)ilower[2]; 711 hupper[0] = (HYPRE_Int)iupper[0]; 712 hupper[1] = (HYPRE_Int)iupper[1]; 713 hupper[2] = (HYPRE_Int)iupper[2]; 714 715 size = 1; 716 for (i = 0; i < 3; i++) size *= (iupper[i] - ilower[i] + 1); 717 718 /* copy x values over to hypre for variable ordering */ 719 if (ordering) { 720 PetscCallExternal(HYPRE_SStructVectorSetConstantValues, mx->ss_b, 0.0); 721 PetscCall(VecGetArrayRead(x, &xx)); 722 for (i = 0; i < nvars; i++) PetscCallExternal(HYPRE_SStructVectorSetBoxValues, mx->ss_b, part, hlower, hupper, i, (HYPRE_Complex *)(xx + (size * i))); 723 PetscCall(VecRestoreArrayRead(x, &xx)); 724 PetscCallExternal(HYPRE_SStructVectorAssemble, mx->ss_b); 725 PetscCallExternal(HYPRE_SStructMatrixMatvec, 1.0, mx->ss_mat, mx->ss_b, 0.0, mx->ss_x); 726 727 /* copy solution values back to PETSc */ 728 PetscCall(VecGetArray(y, &yy)); 729 for (i = 0; i < nvars; i++) PetscCallExternal(HYPRE_SStructVectorGetBoxValues, mx->ss_x, part, hlower, hupper, i, (HYPRE_Complex *)(yy + (size * i))); 730 PetscCall(VecRestoreArray(y, &yy)); 731 } else { /* nodal ordering must be mapped to variable ordering for sys_pfmg */ 732 PetscScalar *z; 733 PetscInt j, k; 734 735 PetscCall(PetscMalloc1(nvars * size, &z)); 736 PetscCallExternal(HYPRE_SStructVectorSetConstantValues, mx->ss_b, 0.0); 737 PetscCall(VecGetArrayRead(x, &xx)); 738 739 /* transform nodal to hypre's variable ordering for sys_pfmg */ 740 for (i = 0; i < size; i++) { 741 k = i * nvars; 742 for (j = 0; j < nvars; j++) z[j * size + i] = xx[k + j]; 743 } 744 for (i = 0; i < nvars; i++) PetscCallExternal(HYPRE_SStructVectorSetBoxValues, mx->ss_b, part, hlower, hupper, i, (HYPRE_Complex *)(z + (size * i))); 745 PetscCall(VecRestoreArrayRead(x, &xx)); 746 PetscCallExternal(HYPRE_SStructVectorAssemble, mx->ss_b); 747 PetscCallExternal(HYPRE_SStructMatrixMatvec, 1.0, mx->ss_mat, mx->ss_b, 0.0, mx->ss_x); 748 749 /* copy solution values back to PETSc */ 750 PetscCall(VecGetArray(y, &yy)); 751 for (i = 0; i < nvars; i++) PetscCallExternal(HYPRE_SStructVectorGetBoxValues, mx->ss_x, part, hlower, hupper, i, (HYPRE_Complex *)(z + (size * i))); 752 /* transform hypre's variable ordering for sys_pfmg to nodal ordering */ 753 for (i = 0; i < size; i++) { 754 k = i * nvars; 755 for (j = 0; j < nvars; j++) yy[k + j] = z[j * size + i]; 756 } 757 PetscCall(VecRestoreArray(y, &yy)); 758 PetscCall(PetscFree(z)); 759 } 760 PetscFunctionReturn(0); 761 } 762 763 PetscErrorCode MatAssemblyEnd_HYPRESStruct(Mat mat, MatAssemblyType mode) { 764 Mat_HYPRESStruct *ex = (Mat_HYPRESStruct *)mat->data; 765 766 PetscFunctionBegin; 767 PetscCallExternal(HYPRE_SStructMatrixAssemble, ex->ss_mat); 768 PetscFunctionReturn(0); 769 } 770 771 PetscErrorCode MatZeroEntries_HYPRESStruct(Mat mat) { 772 PetscFunctionBegin; 773 /* before the DMDA is set to the matrix the zero doesn't need to do anything */ 774 PetscFunctionReturn(0); 775 } 776 777 PetscErrorCode MatDestroy_HYPRESStruct(Mat mat) { 778 Mat_HYPRESStruct *ex = (Mat_HYPRESStruct *)mat->data; 779 ISLocalToGlobalMapping ltog; 780 781 PetscFunctionBegin; 782 PetscCall(DMGetLocalToGlobalMapping(ex->da, <og)); 783 PetscCall(ISLocalToGlobalMappingRestoreIndices(ltog, (const PetscInt **)&ex->gindices)); 784 PetscCallExternal(HYPRE_SStructGraphDestroy, ex->ss_graph); 785 PetscCallExternal(HYPRE_SStructMatrixDestroy, ex->ss_mat); 786 PetscCallExternal(HYPRE_SStructVectorDestroy, ex->ss_x); 787 PetscCallExternal(HYPRE_SStructVectorDestroy, ex->ss_b); 788 PetscCall(PetscObjectDereference((PetscObject)ex->da)); 789 PetscCallMPI(MPI_Comm_free(&(ex->hcomm))); 790 PetscCall(PetscFree(ex)); 791 PetscFunctionReturn(0); 792 } 793 794 PETSC_EXTERN PetscErrorCode MatCreate_HYPRESStruct(Mat B) { 795 Mat_HYPRESStruct *ex; 796 797 PetscFunctionBegin; 798 PetscCall(PetscNew(&ex)); 799 B->data = (void *)ex; 800 B->rmap->bs = 1; 801 B->assembled = PETSC_FALSE; 802 803 B->insertmode = NOT_SET_VALUES; 804 805 B->ops->assemblyend = MatAssemblyEnd_HYPRESStruct; 806 B->ops->mult = MatMult_HYPRESStruct; 807 B->ops->zeroentries = MatZeroEntries_HYPRESStruct; 808 B->ops->destroy = MatDestroy_HYPRESStruct; 809 B->ops->setup = MatSetUp_HYPRESStruct; 810 811 ex->needsinitialization = PETSC_TRUE; 812 813 PetscCallMPI(MPI_Comm_dup(PetscObjectComm((PetscObject)B), &(ex->hcomm))); 814 PetscCall(PetscObjectChangeTypeName((PetscObject)B, MATHYPRESSTRUCT)); 815 PetscFunctionReturn(0); 816 } 817