1 /* 2 This file defines an "generalized" additive Schwarz preconditioner for any Mat implementation. 3 In this version each MPI rank may intersect multiple subdomains and any subdomain may 4 intersect multiple MPI ranks. Intersections of subdomains with MPI ranks are called *local 5 subdomains*. 6 7 N - total number of distinct global subdomains (set explicitly in PCGASMSetTotalSubdomains() or implicitly PCGASMSetSubdomains() and then calculated in PCSetUp_GASM()) 8 n - actual number of local subdomains on this rank (set in PCGASMSetSubdomains() or calculated in PCGASMSetTotalSubdomains()) 9 nmax - maximum number of local subdomains per rank (calculated in PCSetUp_GASM()) 10 */ 11 #include <petsc/private/pcimpl.h> /*I "petscpc.h" I*/ 12 #include <petscdm.h> 13 14 typedef struct { 15 PetscInt N, n, nmax; 16 PetscInt overlap; /* overlap requested by user */ 17 PCGASMType type; /* use reduced interpolation, restriction or both */ 18 PetscBool type_set; /* if user set this value (so won't change it for symmetric problems) */ 19 PetscBool same_subdomain_solvers; /* flag indicating whether all local solvers are same */ 20 PetscBool sort_indices; /* flag to sort subdomain indices */ 21 PetscBool user_subdomains; /* whether the user set explicit subdomain index sets -- keep them on PCReset() */ 22 PetscBool dm_subdomains; /* whether DM is allowed to define subdomains */ 23 PetscBool hierarchicalpartitioning; 24 IS *ois; /* index sets that define the outer (conceptually, overlapping) subdomains */ 25 IS *iis; /* index sets that define the inner (conceptually, nonoverlapping) subdomains */ 26 KSP *ksp; /* linear solvers for each subdomain */ 27 Mat *pmat; /* subdomain block matrices */ 28 Vec gx, gy; /* Merged work vectors */ 29 Vec *x, *y; /* Split work vectors; storage aliases pieces of storage of the above merged vectors. */ 30 VecScatter gorestriction; /* merged restriction to disjoint union of outer subdomains */ 31 VecScatter girestriction; /* merged restriction to disjoint union of inner subdomains */ 32 VecScatter pctoouter; 33 IS permutationIS; 34 Mat permutationP; 35 Mat pcmat; 36 Vec pcx, pcy; 37 } PC_GASM; 38 39 static PetscErrorCode PCGASMComputeGlobalSubdomainNumbering_Private(PC pc, PetscInt **numbering, PetscInt **permutation) 40 { 41 PC_GASM *osm = (PC_GASM *)pc->data; 42 PetscInt i; 43 44 PetscFunctionBegin; 45 /* Determine the number of globally-distinct subdomains and compute a global numbering for them. */ 46 PetscCall(PetscMalloc2(osm->n, numbering, osm->n, permutation)); 47 PetscCall(PetscObjectsListGetGlobalNumbering(PetscObjectComm((PetscObject)pc), osm->n, (PetscObject *)osm->iis, NULL, *numbering)); 48 for (i = 0; i < osm->n; ++i) (*permutation)[i] = i; 49 PetscCall(PetscSortIntWithPermutation(osm->n, *numbering, *permutation)); 50 PetscFunctionReturn(PETSC_SUCCESS); 51 } 52 53 static PetscErrorCode PCGASMSubdomainView_Private(PC pc, PetscInt i, PetscViewer viewer) 54 { 55 PC_GASM *osm = (PC_GASM *)pc->data; 56 PetscInt j, nidx; 57 const PetscInt *idx; 58 PetscViewer sviewer; 59 char *cidx; 60 61 PetscFunctionBegin; 62 PetscCheck(i >= 0 && i < osm->n, PetscObjectComm((PetscObject)viewer), PETSC_ERR_ARG_WRONG, "Invalid subdomain %" PetscInt_FMT ": must nonnegative and less than %" PetscInt_FMT, i, osm->n); 63 /* Inner subdomains. */ 64 PetscCall(ISGetLocalSize(osm->iis[i], &nidx)); 65 /* 66 No more than 15 characters per index plus a space. 67 PetscViewerStringSPrintf requires a string of size at least 2, so use (nidx+1) instead of nidx, 68 in case nidx == 0. That will take care of the space for the trailing '\0' as well. 69 For nidx == 0, the whole string 16 '\0'. 70 */ 71 #define len 16 * (nidx + 1) + 1 72 PetscCall(PetscMalloc1(len, &cidx)); 73 PetscCall(PetscViewerStringOpen(PETSC_COMM_SELF, cidx, len, &sviewer)); 74 #undef len 75 PetscCall(ISGetIndices(osm->iis[i], &idx)); 76 for (j = 0; j < nidx; ++j) PetscCall(PetscViewerStringSPrintf(sviewer, "%" PetscInt_FMT " ", idx[j])); 77 PetscCall(ISRestoreIndices(osm->iis[i], &idx)); 78 PetscCall(PetscViewerDestroy(&sviewer)); 79 PetscCall(PetscViewerASCIIPrintf(viewer, "Inner subdomain:\n")); 80 PetscCall(PetscViewerFlush(viewer)); 81 PetscCall(PetscViewerASCIIPushSynchronized(viewer)); 82 PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, "%s", cidx)); 83 PetscCall(PetscViewerFlush(viewer)); 84 PetscCall(PetscViewerASCIIPopSynchronized(viewer)); 85 PetscCall(PetscViewerASCIIPrintf(viewer, "\n")); 86 PetscCall(PetscViewerFlush(viewer)); 87 PetscCall(PetscFree(cidx)); 88 /* Outer subdomains. */ 89 PetscCall(ISGetLocalSize(osm->ois[i], &nidx)); 90 /* 91 No more than 15 characters per index plus a space. 92 PetscViewerStringSPrintf requires a string of size at least 2, so use (nidx+1) instead of nidx, 93 in case nidx == 0. That will take care of the space for the trailing '\0' as well. 94 For nidx == 0, the whole string 16 '\0'. 95 */ 96 #define len 16 * (nidx + 1) + 1 97 PetscCall(PetscMalloc1(len, &cidx)); 98 PetscCall(PetscViewerStringOpen(PETSC_COMM_SELF, cidx, len, &sviewer)); 99 #undef len 100 PetscCall(ISGetIndices(osm->ois[i], &idx)); 101 for (j = 0; j < nidx; ++j) PetscCall(PetscViewerStringSPrintf(sviewer, "%" PetscInt_FMT " ", idx[j])); 102 PetscCall(PetscViewerDestroy(&sviewer)); 103 PetscCall(ISRestoreIndices(osm->ois[i], &idx)); 104 PetscCall(PetscViewerASCIIPrintf(viewer, "Outer subdomain:\n")); 105 PetscCall(PetscViewerFlush(viewer)); 106 PetscCall(PetscViewerASCIIPushSynchronized(viewer)); 107 PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, "%s", cidx)); 108 PetscCall(PetscViewerFlush(viewer)); 109 PetscCall(PetscViewerASCIIPopSynchronized(viewer)); 110 PetscCall(PetscViewerASCIIPrintf(viewer, "\n")); 111 PetscCall(PetscViewerFlush(viewer)); 112 PetscCall(PetscFree(cidx)); 113 PetscFunctionReturn(PETSC_SUCCESS); 114 } 115 116 static PetscErrorCode PCGASMPrintSubdomains(PC pc) 117 { 118 PC_GASM *osm = (PC_GASM *)pc->data; 119 const char *prefix; 120 char fname[PETSC_MAX_PATH_LEN + 1]; 121 PetscInt l, d, count; 122 PetscBool found; 123 PetscViewer viewer, sviewer = NULL; 124 PetscInt *numbering, *permutation; /* global numbering of locally-supported subdomains and the permutation from the local ordering */ 125 126 PetscFunctionBegin; 127 PetscCall(PCGetOptionsPrefix(pc, &prefix)); 128 PetscCall(PetscOptionsHasName(NULL, prefix, "-pc_gasm_print_subdomains", &found)); 129 if (!found) PetscFunctionReturn(PETSC_SUCCESS); 130 PetscCall(PetscOptionsGetString(NULL, prefix, "-pc_gasm_print_subdomains", fname, sizeof(fname), &found)); 131 if (!found) PetscCall(PetscStrcpy(fname, "stdout")); 132 PetscCall(PetscViewerASCIIOpen(PetscObjectComm((PetscObject)pc), fname, &viewer)); 133 /* 134 Make sure the viewer has a name. Otherwise this may cause a deadlock or other weird errors when creating a subcomm viewer: 135 the subcomm viewer will attempt to inherit the viewer's name, which, if not set, will be constructed collectively on the comm. 136 */ 137 PetscCall(PetscObjectName((PetscObject)viewer)); 138 l = 0; 139 PetscCall(PCGASMComputeGlobalSubdomainNumbering_Private(pc, &numbering, &permutation)); 140 for (count = 0; count < osm->N; ++count) { 141 /* Now let subdomains go one at a time in the global numbering order and print their subdomain/solver info. */ 142 if (l < osm->n) { 143 d = permutation[l]; /* d is the local number of the l-th smallest (in the global ordering) among the locally supported subdomains */ 144 if (numbering[d] == count) { 145 PetscCall(PetscViewerGetSubViewer(viewer, ((PetscObject)osm->ois[d])->comm, &sviewer)); 146 PetscCall(PCGASMSubdomainView_Private(pc, d, sviewer)); 147 PetscCall(PetscViewerRestoreSubViewer(viewer, ((PetscObject)osm->ois[d])->comm, &sviewer)); 148 ++l; 149 } 150 } 151 PetscCallMPI(MPI_Barrier(PetscObjectComm((PetscObject)pc))); 152 } 153 PetscCall(PetscFree2(numbering, permutation)); 154 PetscCall(PetscViewerDestroy(&viewer)); 155 PetscFunctionReturn(PETSC_SUCCESS); 156 } 157 158 static PetscErrorCode PCView_GASM(PC pc, PetscViewer viewer) 159 { 160 PC_GASM *osm = (PC_GASM *)pc->data; 161 const char *prefix; 162 PetscMPIInt rank, size; 163 PetscInt bsz; 164 PetscBool iascii, view_subdomains = PETSC_FALSE; 165 PetscViewer sviewer; 166 PetscInt count, l; 167 char overlap[256] = "user-defined overlap"; 168 char gsubdomains[256] = "unknown total number of subdomains"; 169 char msubdomains[256] = "unknown max number of local subdomains"; 170 PetscInt *numbering, *permutation; /* global numbering of locally-supported subdomains and the permutation from the local ordering */ 171 172 PetscFunctionBegin; 173 PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)pc), &size)); 174 PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)pc), &rank)); 175 176 if (osm->overlap >= 0) PetscCall(PetscSNPrintf(overlap, sizeof(overlap), "requested amount of overlap = %" PetscInt_FMT, osm->overlap)); 177 if (osm->N != PETSC_DETERMINE) PetscCall(PetscSNPrintf(gsubdomains, sizeof(gsubdomains), "total number of subdomains = %" PetscInt_FMT, osm->N)); 178 if (osm->nmax != PETSC_DETERMINE) PetscCall(PetscSNPrintf(msubdomains, sizeof(msubdomains), "max number of local subdomains = %" PetscInt_FMT, osm->nmax)); 179 180 PetscCall(PCGetOptionsPrefix(pc, &prefix)); 181 PetscCall(PetscOptionsGetBool(NULL, prefix, "-pc_gasm_view_subdomains", &view_subdomains, NULL)); 182 183 PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii)); 184 if (iascii) { 185 /* 186 Make sure the viewer has a name. Otherwise this may cause a deadlock when creating a subcomm viewer: 187 the subcomm viewer will attempt to inherit the viewer's name, which, if not set, will be constructed 188 collectively on the comm. 189 */ 190 PetscCall(PetscObjectName((PetscObject)viewer)); 191 PetscCall(PetscViewerASCIIPrintf(viewer, " Restriction/interpolation type: %s\n", PCGASMTypes[osm->type])); 192 PetscCall(PetscViewerASCIIPrintf(viewer, " %s\n", overlap)); 193 PetscCall(PetscViewerASCIIPrintf(viewer, " %s\n", gsubdomains)); 194 PetscCall(PetscViewerASCIIPrintf(viewer, " %s\n", msubdomains)); 195 PetscCall(PetscViewerASCIIPushSynchronized(viewer)); 196 PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, " [%d|%d] number of locally-supported subdomains = %" PetscInt_FMT "\n", rank, size, osm->n)); 197 PetscCall(PetscViewerFlush(viewer)); 198 PetscCall(PetscViewerASCIIPopSynchronized(viewer)); 199 /* Cannot take advantage of osm->same_subdomain_solvers without a global numbering of subdomains. */ 200 PetscCall(PetscViewerASCIIPrintf(viewer, " Subdomain solver info is as follows:\n")); 201 PetscCall(PetscViewerASCIIPushTab(viewer)); 202 PetscCall(PetscViewerASCIIPrintf(viewer, " - - - - - - - - - - - - - - - - - -\n")); 203 /* Make sure that everybody waits for the banner to be printed. */ 204 PetscCallMPI(MPI_Barrier(PetscObjectComm((PetscObject)viewer))); 205 /* Now let subdomains go one at a time in the global numbering order and print their subdomain/solver info. */ 206 PetscCall(PCGASMComputeGlobalSubdomainNumbering_Private(pc, &numbering, &permutation)); 207 l = 0; 208 for (count = 0; count < osm->N; ++count) { 209 PetscMPIInt srank, ssize; 210 if (l < osm->n) { 211 PetscInt d = permutation[l]; /* d is the local number of the l-th smallest (in the global ordering) among the locally supported subdomains */ 212 if (numbering[d] == count) { 213 PetscCallMPI(MPI_Comm_size(((PetscObject)osm->ois[d])->comm, &ssize)); 214 PetscCallMPI(MPI_Comm_rank(((PetscObject)osm->ois[d])->comm, &srank)); 215 PetscCall(PetscViewerGetSubViewer(viewer, ((PetscObject)osm->ois[d])->comm, &sviewer)); 216 PetscCall(ISGetLocalSize(osm->ois[d], &bsz)); 217 PetscCall(PetscViewerASCIISynchronizedPrintf(sviewer, " [%d|%d] (subcomm [%d|%d]) local subdomain number %" PetscInt_FMT ", local size = %" PetscInt_FMT "\n", rank, size, srank, ssize, d, bsz)); 218 PetscCall(PetscViewerFlush(sviewer)); 219 if (view_subdomains) PetscCall(PCGASMSubdomainView_Private(pc, d, sviewer)); 220 if (!pc->setupcalled) { 221 PetscCall(PetscViewerASCIIPrintf(sviewer, " Solver not set up yet: PCSetUp() not yet called\n")); 222 } else { 223 PetscCall(KSPView(osm->ksp[d], sviewer)); 224 } 225 PetscCall(PetscViewerASCIIPrintf(sviewer, " - - - - - - - - - - - - - - - - - -\n")); 226 PetscCall(PetscViewerFlush(sviewer)); 227 PetscCall(PetscViewerRestoreSubViewer(viewer, ((PetscObject)osm->ois[d])->comm, &sviewer)); 228 ++l; 229 } 230 } 231 PetscCallMPI(MPI_Barrier(PetscObjectComm((PetscObject)pc))); 232 } 233 PetscCall(PetscFree2(numbering, permutation)); 234 PetscCall(PetscViewerASCIIPopTab(viewer)); 235 PetscCall(PetscViewerFlush(viewer)); 236 /* this line is needed to match the extra PetscViewerASCIIPushSynchronized() in PetscViewerGetSubViewer() */ 237 PetscCall(PetscViewerASCIIPopSynchronized(viewer)); 238 } 239 PetscFunctionReturn(PETSC_SUCCESS); 240 } 241 242 PETSC_INTERN PetscErrorCode PCGASMCreateLocalSubdomains(Mat A, PetscInt nloc, IS *iis[]); 243 244 PetscErrorCode PCGASMSetHierarchicalPartitioning(PC pc) 245 { 246 PC_GASM *osm = (PC_GASM *)pc->data; 247 MatPartitioning part; 248 MPI_Comm comm; 249 PetscMPIInt size; 250 PetscInt nlocalsubdomains, fromrows_localsize; 251 IS partitioning, fromrows, isn; 252 Vec outervec; 253 254 PetscFunctionBegin; 255 PetscCall(PetscObjectGetComm((PetscObject)pc, &comm)); 256 PetscCallMPI(MPI_Comm_size(comm, &size)); 257 /* we do not need a hierarchical partitioning when 258 * the total number of subdomains is consistent with 259 * the number of MPI tasks. 260 * For the following cases, we do not need to use HP 261 * */ 262 if (osm->N == PETSC_DETERMINE || osm->N >= size || osm->N == 1) PetscFunctionReturn(PETSC_SUCCESS); 263 PetscCheck(size % osm->N == 0, PETSC_COMM_WORLD, PETSC_ERR_ARG_INCOMP, "have to specify the total number of subdomains %" PetscInt_FMT " to be a factor of the number of ranks %d ", osm->N, size); 264 nlocalsubdomains = size / osm->N; 265 osm->n = 1; 266 PetscCall(MatPartitioningCreate(comm, &part)); 267 PetscCall(MatPartitioningSetAdjacency(part, pc->pmat)); 268 PetscCall(MatPartitioningSetType(part, MATPARTITIONINGHIERARCH)); 269 PetscCall(MatPartitioningHierarchicalSetNcoarseparts(part, osm->N)); 270 PetscCall(MatPartitioningHierarchicalSetNfineparts(part, nlocalsubdomains)); 271 PetscCall(MatPartitioningSetFromOptions(part)); 272 /* get new rank owner number of each vertex */ 273 PetscCall(MatPartitioningApply(part, &partitioning)); 274 PetscCall(ISBuildTwoSided(partitioning, NULL, &fromrows)); 275 PetscCall(ISPartitioningToNumbering(partitioning, &isn)); 276 PetscCall(ISDestroy(&isn)); 277 PetscCall(ISGetLocalSize(fromrows, &fromrows_localsize)); 278 PetscCall(MatPartitioningDestroy(&part)); 279 PetscCall(MatCreateVecs(pc->pmat, &outervec, NULL)); 280 PetscCall(VecCreateMPI(comm, fromrows_localsize, PETSC_DETERMINE, &(osm->pcx))); 281 PetscCall(VecDuplicate(osm->pcx, &(osm->pcy))); 282 PetscCall(VecScatterCreate(osm->pcx, NULL, outervec, fromrows, &(osm->pctoouter))); 283 PetscCall(MatCreateSubMatrix(pc->pmat, fromrows, fromrows, MAT_INITIAL_MATRIX, &(osm->permutationP))); 284 PetscCall(PetscObjectReference((PetscObject)fromrows)); 285 osm->permutationIS = fromrows; 286 osm->pcmat = pc->pmat; 287 PetscCall(PetscObjectReference((PetscObject)osm->permutationP)); 288 pc->pmat = osm->permutationP; 289 PetscCall(VecDestroy(&outervec)); 290 PetscCall(ISDestroy(&fromrows)); 291 PetscCall(ISDestroy(&partitioning)); 292 osm->n = PETSC_DETERMINE; 293 PetscFunctionReturn(PETSC_SUCCESS); 294 } 295 296 static PetscErrorCode PCSetUp_GASM(PC pc) 297 { 298 PC_GASM *osm = (PC_GASM *)pc->data; 299 PetscInt i, nInnerIndices, nTotalInnerIndices; 300 PetscMPIInt rank, size; 301 MatReuse scall = MAT_REUSE_MATRIX; 302 KSP ksp; 303 PC subpc; 304 const char *prefix, *pprefix; 305 Vec x, y; 306 PetscInt oni; /* Number of indices in the i-th local outer subdomain. */ 307 const PetscInt *oidxi; /* Indices from the i-th subdomain local outer subdomain. */ 308 PetscInt on; /* Number of indices in the disjoint union of local outer subdomains. */ 309 PetscInt *oidx; /* Indices in the disjoint union of local outer subdomains. */ 310 IS gois; /* Disjoint union the global indices of outer subdomains. */ 311 IS goid; /* Identity IS of the size of the disjoint union of outer subdomains. */ 312 PetscScalar *gxarray, *gyarray; 313 PetscInt gostart; /* Start of locally-owned indices in the vectors -- osm->gx,osm->gy -- over the disjoint union of outer subdomains. */ 314 PetscInt num_subdomains = 0; 315 DM *subdomain_dm = NULL; 316 char **subdomain_names = NULL; 317 PetscInt *numbering; 318 319 PetscFunctionBegin; 320 PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)pc), &size)); 321 PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)pc), &rank)); 322 if (!pc->setupcalled) { 323 /* use a hierarchical partitioning */ 324 if (osm->hierarchicalpartitioning) PetscCall(PCGASMSetHierarchicalPartitioning(pc)); 325 if (osm->n == PETSC_DETERMINE) { 326 if (osm->N != PETSC_DETERMINE) { 327 /* No local subdomains given, but the desired number of total subdomains is known, so construct them accordingly. */ 328 PetscCall(PCGASMCreateSubdomains(pc->pmat, osm->N, &osm->n, &osm->iis)); 329 } else if (osm->dm_subdomains && pc->dm) { 330 /* try pc->dm next, if allowed */ 331 PetscInt d; 332 IS *inner_subdomain_is, *outer_subdomain_is; 333 PetscCall(DMCreateDomainDecomposition(pc->dm, &num_subdomains, &subdomain_names, &inner_subdomain_is, &outer_subdomain_is, &subdomain_dm)); 334 if (num_subdomains) PetscCall(PCGASMSetSubdomains(pc, num_subdomains, inner_subdomain_is, outer_subdomain_is)); 335 for (d = 0; d < num_subdomains; ++d) { 336 if (inner_subdomain_is) PetscCall(ISDestroy(&inner_subdomain_is[d])); 337 if (outer_subdomain_is) PetscCall(ISDestroy(&outer_subdomain_is[d])); 338 } 339 PetscCall(PetscFree(inner_subdomain_is)); 340 PetscCall(PetscFree(outer_subdomain_is)); 341 } else { 342 /* still no subdomains; use one per rank */ 343 osm->nmax = osm->n = 1; 344 PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)pc), &size)); 345 osm->N = size; 346 PetscCall(PCGASMCreateLocalSubdomains(pc->pmat, osm->n, &osm->iis)); 347 } 348 } 349 if (!osm->iis) { 350 /* 351 osm->n was set in PCGASMSetSubdomains(), but the actual subdomains have not been supplied. 352 We create the requisite number of local inner subdomains and then expand them into 353 out subdomains, if necessary. 354 */ 355 PetscCall(PCGASMCreateLocalSubdomains(pc->pmat, osm->n, &osm->iis)); 356 } 357 if (!osm->ois) { 358 /* 359 Initially make outer subdomains the same as inner subdomains. If nonzero additional overlap 360 has been requested, copy the inner subdomains over so they can be modified. 361 */ 362 PetscCall(PetscMalloc1(osm->n, &osm->ois)); 363 for (i = 0; i < osm->n; ++i) { 364 if (osm->overlap > 0 && osm->N > 1) { /* With positive overlap, osm->iis[i] will be modified */ 365 PetscCall(ISDuplicate(osm->iis[i], (osm->ois) + i)); 366 PetscCall(ISCopy(osm->iis[i], osm->ois[i])); 367 } else { 368 PetscCall(PetscObjectReference((PetscObject)((osm->iis)[i]))); 369 osm->ois[i] = osm->iis[i]; 370 } 371 } 372 if (osm->overlap > 0 && osm->N > 1) { 373 /* Extend the "overlapping" regions by a number of steps */ 374 PetscCall(MatIncreaseOverlapSplit(pc->pmat, osm->n, osm->ois, osm->overlap)); 375 } 376 } 377 378 /* Now the subdomains are defined. Determine their global and max local numbers, if necessary. */ 379 if (osm->nmax == PETSC_DETERMINE) { 380 PetscMPIInt inwork, outwork; 381 /* determine global number of subdomains and the max number of local subdomains */ 382 inwork = osm->n; 383 PetscCall(MPIU_Allreduce(&inwork, &outwork, 1, MPI_INT, MPI_MAX, PetscObjectComm((PetscObject)pc))); 384 osm->nmax = outwork; 385 } 386 if (osm->N == PETSC_DETERMINE) { 387 /* Determine the number of globally-distinct subdomains and compute a global numbering for them. */ 388 PetscCall(PetscObjectsListGetGlobalNumbering(PetscObjectComm((PetscObject)pc), osm->n, (PetscObject *)osm->ois, &osm->N, NULL)); 389 } 390 391 if (osm->sort_indices) { 392 for (i = 0; i < osm->n; i++) { 393 PetscCall(ISSort(osm->ois[i])); 394 PetscCall(ISSort(osm->iis[i])); 395 } 396 } 397 PetscCall(PCGetOptionsPrefix(pc, &prefix)); 398 PetscCall(PCGASMPrintSubdomains(pc)); 399 400 /* 401 Merge the ISs, create merged vectors and restrictions. 402 */ 403 /* Merge outer subdomain ISs and construct a restriction onto the disjoint union of local outer subdomains. */ 404 on = 0; 405 for (i = 0; i < osm->n; i++) { 406 PetscCall(ISGetLocalSize(osm->ois[i], &oni)); 407 on += oni; 408 } 409 PetscCall(PetscMalloc1(on, &oidx)); 410 on = 0; 411 /* Merge local indices together */ 412 for (i = 0; i < osm->n; i++) { 413 PetscCall(ISGetLocalSize(osm->ois[i], &oni)); 414 PetscCall(ISGetIndices(osm->ois[i], &oidxi)); 415 PetscCall(PetscArraycpy(oidx + on, oidxi, oni)); 416 PetscCall(ISRestoreIndices(osm->ois[i], &oidxi)); 417 on += oni; 418 } 419 PetscCall(ISCreateGeneral(((PetscObject)(pc))->comm, on, oidx, PETSC_OWN_POINTER, &gois)); 420 nTotalInnerIndices = 0; 421 for (i = 0; i < osm->n; i++) { 422 PetscCall(ISGetLocalSize(osm->iis[i], &nInnerIndices)); 423 nTotalInnerIndices += nInnerIndices; 424 } 425 PetscCall(VecCreateMPI(((PetscObject)(pc))->comm, nTotalInnerIndices, PETSC_DETERMINE, &x)); 426 PetscCall(VecDuplicate(x, &y)); 427 428 PetscCall(VecCreateMPI(PetscObjectComm((PetscObject)pc), on, PETSC_DECIDE, &osm->gx)); 429 PetscCall(VecDuplicate(osm->gx, &osm->gy)); 430 PetscCall(VecGetOwnershipRange(osm->gx, &gostart, NULL)); 431 PetscCall(ISCreateStride(PetscObjectComm((PetscObject)pc), on, gostart, 1, &goid)); 432 /* gois might indices not on local */ 433 PetscCall(VecScatterCreate(x, gois, osm->gx, goid, &(osm->gorestriction))); 434 PetscCall(PetscMalloc1(osm->n, &numbering)); 435 PetscCall(PetscObjectsListGetGlobalNumbering(PetscObjectComm((PetscObject)pc), osm->n, (PetscObject *)osm->ois, NULL, numbering)); 436 PetscCall(VecDestroy(&x)); 437 PetscCall(ISDestroy(&gois)); 438 439 /* Merge inner subdomain ISs and construct a restriction onto the disjoint union of local inner subdomains. */ 440 { 441 PetscInt ini; /* Number of indices the i-th a local inner subdomain. */ 442 PetscInt in; /* Number of indices in the disjoint union of local inner subdomains. */ 443 PetscInt *iidx; /* Global indices in the merged local inner subdomain. */ 444 PetscInt *ioidx; /* Global indices of the disjoint union of inner subdomains within the disjoint union of outer subdomains. */ 445 IS giis; /* IS for the disjoint union of inner subdomains. */ 446 IS giois; /* IS for the disjoint union of inner subdomains within the disjoint union of outer subdomains. */ 447 PetscScalar *array; 448 const PetscInt *indices; 449 PetscInt k; 450 on = 0; 451 for (i = 0; i < osm->n; i++) { 452 PetscCall(ISGetLocalSize(osm->ois[i], &oni)); 453 on += oni; 454 } 455 PetscCall(PetscMalloc1(on, &iidx)); 456 PetscCall(PetscMalloc1(on, &ioidx)); 457 PetscCall(VecGetArray(y, &array)); 458 /* set communicator id to determine where overlap is */ 459 in = 0; 460 for (i = 0; i < osm->n; i++) { 461 PetscCall(ISGetLocalSize(osm->iis[i], &ini)); 462 for (k = 0; k < ini; ++k) array[in + k] = numbering[i]; 463 in += ini; 464 } 465 PetscCall(VecRestoreArray(y, &array)); 466 PetscCall(VecScatterBegin(osm->gorestriction, y, osm->gy, INSERT_VALUES, SCATTER_FORWARD)); 467 PetscCall(VecScatterEnd(osm->gorestriction, y, osm->gy, INSERT_VALUES, SCATTER_FORWARD)); 468 PetscCall(VecGetOwnershipRange(osm->gy, &gostart, NULL)); 469 PetscCall(VecGetArray(osm->gy, &array)); 470 on = 0; 471 in = 0; 472 for (i = 0; i < osm->n; i++) { 473 PetscCall(ISGetLocalSize(osm->ois[i], &oni)); 474 PetscCall(ISGetIndices(osm->ois[i], &indices)); 475 for (k = 0; k < oni; k++) { 476 /* skip overlapping indices to get inner domain */ 477 if (PetscRealPart(array[on + k]) != numbering[i]) continue; 478 iidx[in] = indices[k]; 479 ioidx[in++] = gostart + on + k; 480 } 481 PetscCall(ISRestoreIndices(osm->ois[i], &indices)); 482 on += oni; 483 } 484 PetscCall(VecRestoreArray(osm->gy, &array)); 485 PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject)pc), in, iidx, PETSC_OWN_POINTER, &giis)); 486 PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject)pc), in, ioidx, PETSC_OWN_POINTER, &giois)); 487 PetscCall(VecScatterCreate(y, giis, osm->gy, giois, &osm->girestriction)); 488 PetscCall(VecDestroy(&y)); 489 PetscCall(ISDestroy(&giis)); 490 PetscCall(ISDestroy(&giois)); 491 } 492 PetscCall(ISDestroy(&goid)); 493 PetscCall(PetscFree(numbering)); 494 495 /* Create the subdomain work vectors. */ 496 PetscCall(PetscMalloc1(osm->n, &osm->x)); 497 PetscCall(PetscMalloc1(osm->n, &osm->y)); 498 PetscCall(VecGetArray(osm->gx, &gxarray)); 499 PetscCall(VecGetArray(osm->gy, &gyarray)); 500 for (i = 0, on = 0; i < osm->n; ++i, on += oni) { 501 PetscInt oNi; 502 PetscCall(ISGetLocalSize(osm->ois[i], &oni)); 503 /* on a sub communicator */ 504 PetscCall(ISGetSize(osm->ois[i], &oNi)); 505 PetscCall(VecCreateMPIWithArray(((PetscObject)(osm->ois[i]))->comm, 1, oni, oNi, gxarray + on, &osm->x[i])); 506 PetscCall(VecCreateMPIWithArray(((PetscObject)(osm->ois[i]))->comm, 1, oni, oNi, gyarray + on, &osm->y[i])); 507 } 508 PetscCall(VecRestoreArray(osm->gx, &gxarray)); 509 PetscCall(VecRestoreArray(osm->gy, &gyarray)); 510 /* Create the subdomain solvers */ 511 PetscCall(PetscMalloc1(osm->n, &osm->ksp)); 512 for (i = 0; i < osm->n; i++) { 513 char subprefix[PETSC_MAX_PATH_LEN + 1]; 514 PetscCall(KSPCreate(((PetscObject)(osm->ois[i]))->comm, &ksp)); 515 PetscCall(KSPSetErrorIfNotConverged(ksp, pc->erroriffailure)); 516 PetscCall(PetscObjectIncrementTabLevel((PetscObject)ksp, (PetscObject)pc, 1)); 517 PetscCall(KSPSetType(ksp, KSPPREONLY)); 518 PetscCall(KSPGetPC(ksp, &subpc)); /* Why do we need this here? */ 519 if (subdomain_dm) { 520 PetscCall(KSPSetDM(ksp, subdomain_dm[i])); 521 PetscCall(DMDestroy(subdomain_dm + i)); 522 } 523 PetscCall(PCGetOptionsPrefix(pc, &prefix)); 524 PetscCall(KSPSetOptionsPrefix(ksp, prefix)); 525 if (subdomain_names && subdomain_names[i]) { 526 PetscCall(PetscSNPrintf(subprefix, PETSC_MAX_PATH_LEN, "sub_%s_", subdomain_names[i])); 527 PetscCall(KSPAppendOptionsPrefix(ksp, subprefix)); 528 PetscCall(PetscFree(subdomain_names[i])); 529 } 530 PetscCall(KSPAppendOptionsPrefix(ksp, "sub_")); 531 osm->ksp[i] = ksp; 532 } 533 PetscCall(PetscFree(subdomain_dm)); 534 PetscCall(PetscFree(subdomain_names)); 535 scall = MAT_INITIAL_MATRIX; 536 } else { /* if (pc->setupcalled) */ 537 /* 538 Destroy the submatrices from the previous iteration 539 */ 540 if (pc->flag == DIFFERENT_NONZERO_PATTERN) { 541 PetscCall(MatDestroyMatrices(osm->n, &osm->pmat)); 542 scall = MAT_INITIAL_MATRIX; 543 } 544 if (osm->permutationIS) { 545 PetscCall(MatCreateSubMatrix(pc->pmat, osm->permutationIS, osm->permutationIS, scall, &osm->permutationP)); 546 PetscCall(PetscObjectReference((PetscObject)osm->permutationP)); 547 osm->pcmat = pc->pmat; 548 pc->pmat = osm->permutationP; 549 } 550 } 551 552 /* 553 Extract the submatrices. 554 */ 555 if (size > 1) { 556 PetscCall(MatCreateSubMatricesMPI(pc->pmat, osm->n, osm->ois, osm->ois, scall, &osm->pmat)); 557 } else { 558 PetscCall(MatCreateSubMatrices(pc->pmat, osm->n, osm->ois, osm->ois, scall, &osm->pmat)); 559 } 560 if (scall == MAT_INITIAL_MATRIX) { 561 PetscCall(PetscObjectGetOptionsPrefix((PetscObject)pc->pmat, &pprefix)); 562 for (i = 0; i < osm->n; i++) { PetscCall(PetscObjectSetOptionsPrefix((PetscObject)osm->pmat[i], pprefix)); } 563 } 564 565 /* Return control to the user so that the submatrices can be modified (e.g., to apply 566 different boundary conditions for the submatrices than for the global problem) */ 567 PetscCall(PCModifySubMatrices(pc, osm->n, osm->ois, osm->ois, osm->pmat, pc->modifysubmatricesP)); 568 569 /* 570 Loop over submatrices putting them into local ksps 571 */ 572 for (i = 0; i < osm->n; i++) { 573 PetscCall(KSPSetOperators(osm->ksp[i], osm->pmat[i], osm->pmat[i])); 574 PetscCall(KSPGetOptionsPrefix(osm->ksp[i], &prefix)); 575 PetscCall(MatSetOptionsPrefix(osm->pmat[i], prefix)); 576 if (!pc->setupcalled) PetscCall(KSPSetFromOptions(osm->ksp[i])); 577 } 578 if (osm->pcmat) { 579 PetscCall(MatDestroy(&pc->pmat)); 580 pc->pmat = osm->pcmat; 581 osm->pcmat = NULL; 582 } 583 PetscFunctionReturn(PETSC_SUCCESS); 584 } 585 586 static PetscErrorCode PCSetUpOnBlocks_GASM(PC pc) 587 { 588 PC_GASM *osm = (PC_GASM *)pc->data; 589 PetscInt i; 590 591 PetscFunctionBegin; 592 for (i = 0; i < osm->n; i++) PetscCall(KSPSetUp(osm->ksp[i])); 593 PetscFunctionReturn(PETSC_SUCCESS); 594 } 595 596 static PetscErrorCode PCApply_GASM(PC pc, Vec xin, Vec yout) 597 { 598 PC_GASM *osm = (PC_GASM *)pc->data; 599 PetscInt i; 600 Vec x, y; 601 ScatterMode forward = SCATTER_FORWARD, reverse = SCATTER_REVERSE; 602 603 PetscFunctionBegin; 604 if (osm->pctoouter) { 605 PetscCall(VecScatterBegin(osm->pctoouter, xin, osm->pcx, INSERT_VALUES, SCATTER_REVERSE)); 606 PetscCall(VecScatterEnd(osm->pctoouter, xin, osm->pcx, INSERT_VALUES, SCATTER_REVERSE)); 607 x = osm->pcx; 608 y = osm->pcy; 609 } else { 610 x = xin; 611 y = yout; 612 } 613 /* 614 support for limiting the restriction or interpolation only to the inner 615 subdomain values (leaving the other values 0). 616 */ 617 if (!(osm->type & PC_GASM_RESTRICT)) { 618 /* have to zero the work RHS since scatter may leave some slots empty */ 619 PetscCall(VecZeroEntries(osm->gx)); 620 PetscCall(VecScatterBegin(osm->girestriction, x, osm->gx, INSERT_VALUES, forward)); 621 } else { 622 PetscCall(VecScatterBegin(osm->gorestriction, x, osm->gx, INSERT_VALUES, forward)); 623 } 624 PetscCall(VecZeroEntries(osm->gy)); 625 if (!(osm->type & PC_GASM_RESTRICT)) { 626 PetscCall(VecScatterEnd(osm->girestriction, x, osm->gx, INSERT_VALUES, forward)); 627 } else { 628 PetscCall(VecScatterEnd(osm->gorestriction, x, osm->gx, INSERT_VALUES, forward)); 629 } 630 /* do the subdomain solves */ 631 for (i = 0; i < osm->n; ++i) { 632 PetscCall(KSPSolve(osm->ksp[i], osm->x[i], osm->y[i])); 633 PetscCall(KSPCheckSolve(osm->ksp[i], pc, osm->y[i])); 634 } 635 /* do we need to zero y? */ 636 PetscCall(VecZeroEntries(y)); 637 if (!(osm->type & PC_GASM_INTERPOLATE)) { 638 PetscCall(VecScatterBegin(osm->girestriction, osm->gy, y, ADD_VALUES, reverse)); 639 PetscCall(VecScatterEnd(osm->girestriction, osm->gy, y, ADD_VALUES, reverse)); 640 } else { 641 PetscCall(VecScatterBegin(osm->gorestriction, osm->gy, y, ADD_VALUES, reverse)); 642 PetscCall(VecScatterEnd(osm->gorestriction, osm->gy, y, ADD_VALUES, reverse)); 643 } 644 if (osm->pctoouter) { 645 PetscCall(VecScatterBegin(osm->pctoouter, y, yout, INSERT_VALUES, SCATTER_FORWARD)); 646 PetscCall(VecScatterEnd(osm->pctoouter, y, yout, INSERT_VALUES, SCATTER_FORWARD)); 647 } 648 PetscFunctionReturn(PETSC_SUCCESS); 649 } 650 651 static PetscErrorCode PCMatApply_GASM(PC pc, Mat Xin, Mat Yout) 652 { 653 PC_GASM *osm = (PC_GASM *)pc->data; 654 Mat X, Y, O = NULL, Z, W; 655 Vec x, y; 656 PetscInt i, m, M, N; 657 ScatterMode forward = SCATTER_FORWARD, reverse = SCATTER_REVERSE; 658 659 PetscFunctionBegin; 660 PetscCheck(osm->n == 1, PetscObjectComm((PetscObject)pc), PETSC_ERR_SUP, "Not yet implemented"); 661 PetscCall(MatGetSize(Xin, NULL, &N)); 662 if (osm->pctoouter) { 663 PetscCall(VecGetLocalSize(osm->pcx, &m)); 664 PetscCall(VecGetSize(osm->pcx, &M)); 665 PetscCall(MatCreateDense(PetscObjectComm((PetscObject)osm->ois[0]), m, PETSC_DECIDE, M, N, NULL, &O)); 666 for (i = 0; i < N; ++i) { 667 PetscCall(MatDenseGetColumnVecRead(Xin, i, &x)); 668 PetscCall(MatDenseGetColumnVecWrite(O, i, &y)); 669 PetscCall(VecScatterBegin(osm->pctoouter, x, y, INSERT_VALUES, SCATTER_REVERSE)); 670 PetscCall(VecScatterEnd(osm->pctoouter, x, y, INSERT_VALUES, SCATTER_REVERSE)); 671 PetscCall(MatDenseRestoreColumnVecWrite(O, i, &y)); 672 PetscCall(MatDenseRestoreColumnVecRead(Xin, i, &x)); 673 } 674 X = Y = O; 675 } else { 676 X = Xin; 677 Y = Yout; 678 } 679 /* 680 support for limiting the restriction or interpolation only to the inner 681 subdomain values (leaving the other values 0). 682 */ 683 PetscCall(VecGetLocalSize(osm->x[0], &m)); 684 PetscCall(VecGetSize(osm->x[0], &M)); 685 PetscCall(MatCreateDense(PetscObjectComm((PetscObject)osm->ois[0]), m, PETSC_DECIDE, M, N, NULL, &Z)); 686 for (i = 0; i < N; ++i) { 687 PetscCall(MatDenseGetColumnVecRead(X, i, &x)); 688 PetscCall(MatDenseGetColumnVecWrite(Z, i, &y)); 689 if (!(osm->type & PC_GASM_RESTRICT)) { 690 /* have to zero the work RHS since scatter may leave some slots empty */ 691 PetscCall(VecZeroEntries(y)); 692 PetscCall(VecScatterBegin(osm->girestriction, x, y, INSERT_VALUES, forward)); 693 PetscCall(VecScatterEnd(osm->girestriction, x, y, INSERT_VALUES, forward)); 694 } else { 695 PetscCall(VecScatterBegin(osm->gorestriction, x, y, INSERT_VALUES, forward)); 696 PetscCall(VecScatterEnd(osm->gorestriction, x, y, INSERT_VALUES, forward)); 697 } 698 PetscCall(MatDenseRestoreColumnVecWrite(Z, i, &y)); 699 PetscCall(MatDenseRestoreColumnVecRead(X, i, &x)); 700 } 701 PetscCall(MatCreateDense(PetscObjectComm((PetscObject)osm->ois[0]), m, PETSC_DECIDE, M, N, NULL, &W)); 702 PetscCall(MatSetOption(Z, MAT_NO_OFF_PROC_ENTRIES, PETSC_TRUE)); 703 PetscCall(MatAssemblyBegin(Z, MAT_FINAL_ASSEMBLY)); 704 PetscCall(MatAssemblyEnd(Z, MAT_FINAL_ASSEMBLY)); 705 /* do the subdomain solve */ 706 PetscCall(KSPMatSolve(osm->ksp[0], Z, W)); 707 PetscCall(KSPCheckSolve(osm->ksp[0], pc, NULL)); 708 PetscCall(MatDestroy(&Z)); 709 /* do we need to zero y? */ 710 PetscCall(MatZeroEntries(Y)); 711 for (i = 0; i < N; ++i) { 712 PetscCall(MatDenseGetColumnVecWrite(Y, i, &y)); 713 PetscCall(MatDenseGetColumnVecRead(W, i, &x)); 714 if (!(osm->type & PC_GASM_INTERPOLATE)) { 715 PetscCall(VecScatterBegin(osm->girestriction, x, y, ADD_VALUES, reverse)); 716 PetscCall(VecScatterEnd(osm->girestriction, x, y, ADD_VALUES, reverse)); 717 } else { 718 PetscCall(VecScatterBegin(osm->gorestriction, x, y, ADD_VALUES, reverse)); 719 PetscCall(VecScatterEnd(osm->gorestriction, x, y, ADD_VALUES, reverse)); 720 } 721 PetscCall(MatDenseRestoreColumnVecRead(W, i, &x)); 722 if (osm->pctoouter) { 723 PetscCall(MatDenseGetColumnVecWrite(Yout, i, &x)); 724 PetscCall(VecScatterBegin(osm->pctoouter, y, x, INSERT_VALUES, SCATTER_FORWARD)); 725 PetscCall(VecScatterEnd(osm->pctoouter, y, x, INSERT_VALUES, SCATTER_FORWARD)); 726 PetscCall(MatDenseRestoreColumnVecRead(Yout, i, &x)); 727 } 728 PetscCall(MatDenseRestoreColumnVecWrite(Y, i, &y)); 729 } 730 PetscCall(MatDestroy(&W)); 731 PetscCall(MatDestroy(&O)); 732 PetscFunctionReturn(PETSC_SUCCESS); 733 } 734 735 static PetscErrorCode PCApplyTranspose_GASM(PC pc, Vec xin, Vec yout) 736 { 737 PC_GASM *osm = (PC_GASM *)pc->data; 738 PetscInt i; 739 Vec x, y; 740 ScatterMode forward = SCATTER_FORWARD, reverse = SCATTER_REVERSE; 741 742 PetscFunctionBegin; 743 if (osm->pctoouter) { 744 PetscCall(VecScatterBegin(osm->pctoouter, xin, osm->pcx, INSERT_VALUES, SCATTER_REVERSE)); 745 PetscCall(VecScatterEnd(osm->pctoouter, xin, osm->pcx, INSERT_VALUES, SCATTER_REVERSE)); 746 x = osm->pcx; 747 y = osm->pcy; 748 } else { 749 x = xin; 750 y = yout; 751 } 752 /* 753 Support for limiting the restriction or interpolation to only local 754 subdomain values (leaving the other values 0). 755 756 Note: these are reversed from the PCApply_GASM() because we are applying the 757 transpose of the three terms 758 */ 759 if (!(osm->type & PC_GASM_INTERPOLATE)) { 760 /* have to zero the work RHS since scatter may leave some slots empty */ 761 PetscCall(VecZeroEntries(osm->gx)); 762 PetscCall(VecScatterBegin(osm->girestriction, x, osm->gx, INSERT_VALUES, forward)); 763 } else { 764 PetscCall(VecScatterBegin(osm->gorestriction, x, osm->gx, INSERT_VALUES, forward)); 765 } 766 PetscCall(VecZeroEntries(osm->gy)); 767 if (!(osm->type & PC_GASM_INTERPOLATE)) { 768 PetscCall(VecScatterEnd(osm->girestriction, x, osm->gx, INSERT_VALUES, forward)); 769 } else { 770 PetscCall(VecScatterEnd(osm->gorestriction, x, osm->gx, INSERT_VALUES, forward)); 771 } 772 /* do the local solves */ 773 for (i = 0; i < osm->n; ++i) { /* Note that the solves are local, so we can go to osm->n, rather than osm->nmax. */ 774 PetscCall(KSPSolveTranspose(osm->ksp[i], osm->x[i], osm->y[i])); 775 PetscCall(KSPCheckSolve(osm->ksp[i], pc, osm->y[i])); 776 } 777 PetscCall(VecZeroEntries(y)); 778 if (!(osm->type & PC_GASM_RESTRICT)) { 779 PetscCall(VecScatterBegin(osm->girestriction, osm->gy, y, ADD_VALUES, reverse)); 780 PetscCall(VecScatterEnd(osm->girestriction, osm->gy, y, ADD_VALUES, reverse)); 781 } else { 782 PetscCall(VecScatterBegin(osm->gorestriction, osm->gy, y, ADD_VALUES, reverse)); 783 PetscCall(VecScatterEnd(osm->gorestriction, osm->gy, y, ADD_VALUES, reverse)); 784 } 785 if (osm->pctoouter) { 786 PetscCall(VecScatterBegin(osm->pctoouter, y, yout, INSERT_VALUES, SCATTER_FORWARD)); 787 PetscCall(VecScatterEnd(osm->pctoouter, y, yout, INSERT_VALUES, SCATTER_FORWARD)); 788 } 789 PetscFunctionReturn(PETSC_SUCCESS); 790 } 791 792 static PetscErrorCode PCReset_GASM(PC pc) 793 { 794 PC_GASM *osm = (PC_GASM *)pc->data; 795 PetscInt i; 796 797 PetscFunctionBegin; 798 if (osm->ksp) { 799 for (i = 0; i < osm->n; i++) PetscCall(KSPReset(osm->ksp[i])); 800 } 801 if (osm->pmat) { 802 if (osm->n > 0) { 803 PetscMPIInt size; 804 PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)pc), &size)); 805 if (size > 1) { 806 /* osm->pmat is created by MatCreateSubMatricesMPI(), cannot use MatDestroySubMatrices() */ 807 PetscCall(MatDestroyMatrices(osm->n, &osm->pmat)); 808 } else { 809 PetscCall(MatDestroySubMatrices(osm->n, &osm->pmat)); 810 } 811 } 812 } 813 if (osm->x) { 814 for (i = 0; i < osm->n; i++) { 815 PetscCall(VecDestroy(&osm->x[i])); 816 PetscCall(VecDestroy(&osm->y[i])); 817 } 818 } 819 PetscCall(VecDestroy(&osm->gx)); 820 PetscCall(VecDestroy(&osm->gy)); 821 822 PetscCall(VecScatterDestroy(&osm->gorestriction)); 823 PetscCall(VecScatterDestroy(&osm->girestriction)); 824 if (!osm->user_subdomains) { 825 PetscCall(PCGASMDestroySubdomains(osm->n, &osm->ois, &osm->iis)); 826 osm->N = PETSC_DETERMINE; 827 osm->nmax = PETSC_DETERMINE; 828 } 829 if (osm->pctoouter) PetscCall(VecScatterDestroy(&(osm->pctoouter))); 830 if (osm->permutationIS) PetscCall(ISDestroy(&(osm->permutationIS))); 831 if (osm->pcx) PetscCall(VecDestroy(&(osm->pcx))); 832 if (osm->pcy) PetscCall(VecDestroy(&(osm->pcy))); 833 if (osm->permutationP) PetscCall(MatDestroy(&(osm->permutationP))); 834 if (osm->pcmat) PetscCall(MatDestroy(&osm->pcmat)); 835 PetscFunctionReturn(PETSC_SUCCESS); 836 } 837 838 static PetscErrorCode PCDestroy_GASM(PC pc) 839 { 840 PC_GASM *osm = (PC_GASM *)pc->data; 841 PetscInt i; 842 843 PetscFunctionBegin; 844 PetscCall(PCReset_GASM(pc)); 845 /* PCReset will not destroy subdomains, if user_subdomains is true. */ 846 PetscCall(PCGASMDestroySubdomains(osm->n, &osm->ois, &osm->iis)); 847 if (osm->ksp) { 848 for (i = 0; i < osm->n; i++) PetscCall(KSPDestroy(&osm->ksp[i])); 849 PetscCall(PetscFree(osm->ksp)); 850 } 851 PetscCall(PetscFree(osm->x)); 852 PetscCall(PetscFree(osm->y)); 853 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGASMSetSubdomains_C", NULL)); 854 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGASMSetOverlap_C", NULL)); 855 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGASMSetType_C", NULL)); 856 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGASMSetSortIndices_C", NULL)); 857 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGASMGetSubKSP_C", NULL)); 858 PetscCall(PetscFree(pc->data)); 859 PetscFunctionReturn(PETSC_SUCCESS); 860 } 861 862 static PetscErrorCode PCSetFromOptions_GASM(PC pc, PetscOptionItems *PetscOptionsObject) 863 { 864 PC_GASM *osm = (PC_GASM *)pc->data; 865 PetscInt blocks, ovl; 866 PetscBool flg; 867 PCGASMType gasmtype; 868 869 PetscFunctionBegin; 870 PetscOptionsHeadBegin(PetscOptionsObject, "Generalized additive Schwarz options"); 871 PetscCall(PetscOptionsBool("-pc_gasm_use_dm_subdomains", "If subdomains aren't set, use DMCreateDomainDecomposition() to define subdomains.", "PCGASMSetUseDMSubdomains", osm->dm_subdomains, &osm->dm_subdomains, &flg)); 872 PetscCall(PetscOptionsInt("-pc_gasm_total_subdomains", "Total number of subdomains across communicator", "PCGASMSetTotalSubdomains", osm->N, &blocks, &flg)); 873 if (flg) PetscCall(PCGASMSetTotalSubdomains(pc, blocks)); 874 PetscCall(PetscOptionsInt("-pc_gasm_overlap", "Number of overlapping degrees of freedom", "PCGASMSetOverlap", osm->overlap, &ovl, &flg)); 875 if (flg) { 876 PetscCall(PCGASMSetOverlap(pc, ovl)); 877 osm->dm_subdomains = PETSC_FALSE; 878 } 879 flg = PETSC_FALSE; 880 PetscCall(PetscOptionsEnum("-pc_gasm_type", "Type of restriction/extension", "PCGASMSetType", PCGASMTypes, (PetscEnum)osm->type, (PetscEnum *)&gasmtype, &flg)); 881 if (flg) PetscCall(PCGASMSetType(pc, gasmtype)); 882 PetscCall(PetscOptionsBool("-pc_gasm_use_hierachical_partitioning", "use hierarchical partitioning", NULL, osm->hierarchicalpartitioning, &osm->hierarchicalpartitioning, &flg)); 883 PetscOptionsHeadEnd(); 884 PetscFunctionReturn(PETSC_SUCCESS); 885 } 886 887 /*------------------------------------------------------------------------------------*/ 888 889 /*@ 890 PCGASMSetTotalSubdomains - sets the total number of subdomains to use across the communicator for `PCGASM` 891 892 Logically Collective 893 894 Input Parameters: 895 + pc - the preconditioner 896 - N - total number of subdomains 897 898 Level: beginner 899 900 .seealso: `PCGASM`, `PCGASMSetSubdomains()`, `PCGASMSetOverlap()` 901 `PCGASMCreateSubdomains2D()` 902 @*/ 903 PetscErrorCode PCGASMSetTotalSubdomains(PC pc, PetscInt N) 904 { 905 PC_GASM *osm = (PC_GASM *)pc->data; 906 PetscMPIInt size, rank; 907 908 PetscFunctionBegin; 909 PetscCheck(N >= 1, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Total number of subdomains must be 1 or more, got N = %" PetscInt_FMT, N); 910 PetscCheck(!pc->setupcalled, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_WRONGSTATE, "PCGASMSetTotalSubdomains() should be called before calling PCSetUp()."); 911 912 PetscCall(PCGASMDestroySubdomains(osm->n, &osm->iis, &osm->ois)); 913 osm->ois = osm->iis = NULL; 914 915 PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)pc), &size)); 916 PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)pc), &rank)); 917 osm->N = N; 918 osm->n = PETSC_DETERMINE; 919 osm->nmax = PETSC_DETERMINE; 920 osm->dm_subdomains = PETSC_FALSE; 921 PetscFunctionReturn(PETSC_SUCCESS); 922 } 923 924 static PetscErrorCode PCGASMSetSubdomains_GASM(PC pc, PetscInt n, IS iis[], IS ois[]) 925 { 926 PC_GASM *osm = (PC_GASM *)pc->data; 927 PetscInt i; 928 929 PetscFunctionBegin; 930 PetscCheck(n >= 1, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Each MPI rank must have 1 or more subdomains, got n = %" PetscInt_FMT, n); 931 PetscCheck(!pc->setupcalled, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_WRONGSTATE, "PCGASMSetSubdomains() should be called before calling PCSetUp()."); 932 933 PetscCall(PCGASMDestroySubdomains(osm->n, &osm->iis, &osm->ois)); 934 osm->iis = osm->ois = NULL; 935 osm->n = n; 936 osm->N = PETSC_DETERMINE; 937 osm->nmax = PETSC_DETERMINE; 938 if (ois) { 939 PetscCall(PetscMalloc1(n, &osm->ois)); 940 for (i = 0; i < n; i++) { 941 PetscCall(PetscObjectReference((PetscObject)ois[i])); 942 osm->ois[i] = ois[i]; 943 } 944 /* 945 Since the user set the outer subdomains, even if nontrivial overlap was requested via PCGASMSetOverlap(), 946 it will be ignored. To avoid confusion later on (e.g., when viewing the PC), the overlap size is set to -1. 947 */ 948 osm->overlap = -1; 949 /* inner subdomains must be provided */ 950 PetscCheck(iis, PETSC_COMM_SELF, PETSC_ERR_ARG_NULL, "inner indices have to be provided "); 951 } /* end if */ 952 if (iis) { 953 PetscCall(PetscMalloc1(n, &osm->iis)); 954 for (i = 0; i < n; i++) { 955 PetscCall(PetscObjectReference((PetscObject)iis[i])); 956 osm->iis[i] = iis[i]; 957 } 958 if (!ois) { 959 osm->ois = NULL; 960 /* if user does not provide outer indices, we will create the corresponding outer indices using osm->overlap =1 in PCSetUp_GASM */ 961 } 962 } 963 if (PetscDefined(USE_DEBUG)) { 964 PetscInt j, rstart, rend, *covered, lsize; 965 const PetscInt *indices; 966 /* check if the inner indices cover and only cover the local portion of the preconditioning matrix */ 967 PetscCall(MatGetOwnershipRange(pc->pmat, &rstart, &rend)); 968 PetscCall(PetscCalloc1(rend - rstart, &covered)); 969 /* check if the current MPI rank owns indices from others */ 970 for (i = 0; i < n; i++) { 971 PetscCall(ISGetIndices(osm->iis[i], &indices)); 972 PetscCall(ISGetLocalSize(osm->iis[i], &lsize)); 973 for (j = 0; j < lsize; j++) { 974 PetscCheck(indices[j] >= rstart && indices[j] < rend, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "inner subdomains can not own an index %" PetscInt_FMT " from other ranks", indices[j]); 975 PetscCheck(covered[indices[j] - rstart] != 1, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "inner subdomains can not have an overlapping index %" PetscInt_FMT " ", indices[j]); 976 covered[indices[j] - rstart] = 1; 977 } 978 PetscCall(ISRestoreIndices(osm->iis[i], &indices)); 979 } 980 /* check if we miss any indices */ 981 for (i = rstart; i < rend; i++) PetscCheck(covered[i - rstart], PETSC_COMM_SELF, PETSC_ERR_ARG_NULL, "local entity %" PetscInt_FMT " was not covered by inner subdomains", i); 982 PetscCall(PetscFree(covered)); 983 } 984 if (iis) osm->user_subdomains = PETSC_TRUE; 985 PetscFunctionReturn(PETSC_SUCCESS); 986 } 987 988 static PetscErrorCode PCGASMSetOverlap_GASM(PC pc, PetscInt ovl) 989 { 990 PC_GASM *osm = (PC_GASM *)pc->data; 991 992 PetscFunctionBegin; 993 PetscCheck(ovl >= 0, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_OUTOFRANGE, "Negative overlap value requested"); 994 PetscCheck(!pc->setupcalled || ovl == osm->overlap, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_WRONGSTATE, "PCGASMSetOverlap() should be called before PCSetUp()."); 995 if (!pc->setupcalled) osm->overlap = ovl; 996 PetscFunctionReturn(PETSC_SUCCESS); 997 } 998 999 static PetscErrorCode PCGASMSetType_GASM(PC pc, PCGASMType type) 1000 { 1001 PC_GASM *osm = (PC_GASM *)pc->data; 1002 1003 PetscFunctionBegin; 1004 osm->type = type; 1005 osm->type_set = PETSC_TRUE; 1006 PetscFunctionReturn(PETSC_SUCCESS); 1007 } 1008 1009 static PetscErrorCode PCGASMSetSortIndices_GASM(PC pc, PetscBool doSort) 1010 { 1011 PC_GASM *osm = (PC_GASM *)pc->data; 1012 1013 PetscFunctionBegin; 1014 osm->sort_indices = doSort; 1015 PetscFunctionReturn(PETSC_SUCCESS); 1016 } 1017 1018 /* 1019 FIXME: This routine might need to be modified now that multiple ranks per subdomain are allowed. 1020 In particular, it would upset the global subdomain number calculation. 1021 */ 1022 static PetscErrorCode PCGASMGetSubKSP_GASM(PC pc, PetscInt *n, PetscInt *first, KSP **ksp) 1023 { 1024 PC_GASM *osm = (PC_GASM *)pc->data; 1025 1026 PetscFunctionBegin; 1027 PetscCheck(osm->n >= 1, PetscObjectComm((PetscObject)pc), PETSC_ERR_ORDER, "Need to call PCSetUp() on PC (or KSPSetUp() on the outer KSP object) before calling here"); 1028 1029 if (n) *n = osm->n; 1030 if (first) { 1031 PetscCallMPI(MPI_Scan(&osm->n, first, 1, MPIU_INT, MPI_SUM, PetscObjectComm((PetscObject)pc))); 1032 *first -= osm->n; 1033 } 1034 if (ksp) { 1035 /* Assume that local solves are now different; not necessarily 1036 true, though! This flag is used only for PCView_GASM() */ 1037 *ksp = osm->ksp; 1038 osm->same_subdomain_solvers = PETSC_FALSE; 1039 } 1040 PetscFunctionReturn(PETSC_SUCCESS); 1041 } /* PCGASMGetSubKSP_GASM() */ 1042 1043 /*@C 1044 PCGASMSetSubdomains - Sets the subdomains for this MPI rank 1045 for the additive Schwarz preconditioner with multiple MPI ranks per subdomain, `PCGASM` 1046 1047 Collective 1048 1049 Input Parameters: 1050 + pc - the preconditioner object 1051 . n - the number of subdomains for this MPI rank 1052 . iis - the index sets that define the inner subdomains (or NULL for PETSc to determine subdomains) 1053 - ois - the index sets that define the outer subdomains (or NULL to use the same as iis, or to construct by expanding iis by the requested overlap) 1054 1055 Notes: 1056 The `IS` indices use the parallel, global numbering of the vector entries. 1057 Inner subdomains are those where the correction is applied. 1058 Outer subdomains are those where the residual necessary to obtain the 1059 corrections is obtained (see `PCGASMType` for the use of inner/outer subdomains). 1060 Both inner and outer subdomains can extend over several MPI ranks. 1061 This rank's portion of a subdomain is known as a local subdomain. 1062 1063 Inner subdomains can not overlap with each other, do not have any entities from remote ranks, 1064 and have to cover the entire local subdomain owned by the current rank. The index sets on each 1065 rank should be ordered such that the ith local subdomain is connected to the ith remote subdomain 1066 on another MPI rank. 1067 1068 By default the `PGASM` preconditioner uses 1 (local) subdomain per MPI rank. 1069 1070 Level: advanced 1071 1072 .seealso: `PCGASM`, `PCGASMSetOverlap()`, `PCGASMGetSubKSP()`, 1073 `PCGASMCreateSubdomains2D()`, `PCGASMGetSubdomains()` 1074 @*/ 1075 PetscErrorCode PCGASMSetSubdomains(PC pc, PetscInt n, IS iis[], IS ois[]) 1076 { 1077 PC_GASM *osm = (PC_GASM *)pc->data; 1078 1079 PetscFunctionBegin; 1080 PetscValidHeaderSpecific(pc, PC_CLASSID, 1); 1081 PetscTryMethod(pc, "PCGASMSetSubdomains_C", (PC, PetscInt, IS[], IS[]), (pc, n, iis, ois)); 1082 osm->dm_subdomains = PETSC_FALSE; 1083 PetscFunctionReturn(PETSC_SUCCESS); 1084 } 1085 1086 /*@ 1087 PCGASMSetOverlap - Sets the overlap between a pair of subdomains for the 1088 additive Schwarz preconditioner `PCGASM`. Either all or no MPI ranks in the 1089 pc communicator must call this routine. 1090 1091 Logically Collective 1092 1093 Input Parameters: 1094 + pc - the preconditioner context 1095 - ovl - the amount of overlap between subdomains (ovl >= 0, default value = 0) 1096 1097 Options Database Key: 1098 . -pc_gasm_overlap <overlap> - Sets overlap 1099 1100 Notes: 1101 By default the `PCGASM` preconditioner uses 1 subdomain per rank. To use 1102 multiple subdomain per perocessor or "straddling" subdomains that intersect 1103 multiple ranks use `PCGASMSetSubdomains()` (or option -pc_gasm_total_subdomains <n>). 1104 1105 The overlap defaults to 0, so if one desires that no additional 1106 overlap be computed beyond what may have been set with a call to 1107 `PCGASMSetSubdomains()`, then ovl must be set to be 0. In particular, if one does 1108 not explicitly set the subdomains in application code, then all overlap would be computed 1109 internally by PETSc, and using an overlap of 0 would result in an `PCGASM` 1110 variant that is equivalent to the block Jacobi preconditioner. 1111 1112 One can define initial index sets with any overlap via 1113 `PCGASMSetSubdomains()`; the routine `PCGASMSetOverlap()` merely allows 1114 PETSc to extend that overlap further, if desired. 1115 1116 Level: intermediate 1117 1118 .seealso: `PCGASM`, `PCGASMSetSubdomains()`, `PCGASMGetSubKSP()`, 1119 `PCGASMCreateSubdomains2D()`, `PCGASMGetSubdomains()` 1120 @*/ 1121 PetscErrorCode PCGASMSetOverlap(PC pc, PetscInt ovl) 1122 { 1123 PC_GASM *osm = (PC_GASM *)pc->data; 1124 1125 PetscFunctionBegin; 1126 PetscValidHeaderSpecific(pc, PC_CLASSID, 1); 1127 PetscValidLogicalCollectiveInt(pc, ovl, 2); 1128 PetscTryMethod(pc, "PCGASMSetOverlap_C", (PC, PetscInt), (pc, ovl)); 1129 osm->dm_subdomains = PETSC_FALSE; 1130 PetscFunctionReturn(PETSC_SUCCESS); 1131 } 1132 1133 /*@ 1134 PCGASMSetType - Sets the type of restriction and interpolation used 1135 for local problems in the `PCGASM` additive Schwarz method. 1136 1137 Logically Collective 1138 1139 Input Parameters: 1140 + pc - the preconditioner context 1141 - type - variant of `PCGASM`, one of 1142 .vb 1143 `PC_GASM_BASIC` - full interpolation and restriction 1144 `PC_GASM_RESTRICT` - full restriction, local MPI rank interpolation 1145 `PC_GASM_INTERPOLATE` - full interpolation, local MPI rank restriction 1146 `PC_GASM_NONE` - local MPI rank restriction and interpolation 1147 .ve 1148 1149 Options Database Key: 1150 . -pc_gasm_type [basic,restrict,interpolate,none] - Sets GASM type 1151 1152 Level: intermediate 1153 1154 .seealso: `PCGASM`, `PCGASMSetSubdomains()`, `PCGASMGetSubKSP()`, 1155 `PCGASMCreateSubdomains2D()`, `PCASM`, `PCASMSetType()` 1156 @*/ 1157 PetscErrorCode PCGASMSetType(PC pc, PCGASMType type) 1158 { 1159 PetscFunctionBegin; 1160 PetscValidHeaderSpecific(pc, PC_CLASSID, 1); 1161 PetscValidLogicalCollectiveEnum(pc, type, 2); 1162 PetscTryMethod(pc, "PCGASMSetType_C", (PC, PCGASMType), (pc, type)); 1163 PetscFunctionReturn(PETSC_SUCCESS); 1164 } 1165 1166 /*@ 1167 PCGASMSetSortIndices - Determines whether subdomain indices are sorted. 1168 1169 Logically Collective 1170 1171 Input Parameters: 1172 + pc - the preconditioner context 1173 - doSort - sort the subdomain indices 1174 1175 Level: intermediate 1176 1177 .seealso: `PCGASM`, `PCGASMSetSubdomains()`, `PCGASMGetSubKSP()`, 1178 `PCGASMCreateSubdomains2D()` 1179 @*/ 1180 PetscErrorCode PCGASMSetSortIndices(PC pc, PetscBool doSort) 1181 { 1182 PetscFunctionBegin; 1183 PetscValidHeaderSpecific(pc, PC_CLASSID, 1); 1184 PetscValidLogicalCollectiveBool(pc, doSort, 2); 1185 PetscTryMethod(pc, "PCGASMSetSortIndices_C", (PC, PetscBool), (pc, doSort)); 1186 PetscFunctionReturn(PETSC_SUCCESS); 1187 } 1188 1189 /*@C 1190 PCGASMGetSubKSP - Gets the local `KSP` contexts for all subdomains on 1191 this MPI rank. 1192 1193 Collective iff first_local is requested 1194 1195 Input Parameter: 1196 . pc - the preconditioner context 1197 1198 Output Parameters: 1199 + n_local - the number of blocks on this MPI rank or NULL 1200 . first_local - the global number of the first block on this rank or NULL, 1201 all ranks must request or all must pass NULL 1202 - ksp - the array of `KSP` contexts 1203 1204 Note: 1205 After `PCGASMGetSubKSP()` the array of `KSP`es is not to be freed 1206 1207 Currently for some matrix implementations only 1 block per MPI rank 1208 is supported. 1209 1210 You must call `KSPSetUp()` before calling `PCGASMGetSubKSP()`. 1211 1212 Level: advanced 1213 1214 .seealso: `PCGASM`, `PCGASMSetSubdomains()`, `PCGASMSetOverlap()`, 1215 `PCGASMCreateSubdomains2D()`, 1216 @*/ 1217 PetscErrorCode PCGASMGetSubKSP(PC pc, PetscInt *n_local, PetscInt *first_local, KSP *ksp[]) 1218 { 1219 PetscFunctionBegin; 1220 PetscValidHeaderSpecific(pc, PC_CLASSID, 1); 1221 PetscUseMethod(pc, "PCGASMGetSubKSP_C", (PC, PetscInt *, PetscInt *, KSP **), (pc, n_local, first_local, ksp)); 1222 PetscFunctionReturn(PETSC_SUCCESS); 1223 } 1224 1225 /*MC 1226 PCGASM - Use the (restricted) additive Schwarz method, each block is (approximately) solved with 1227 its own `KSP` object on a subset of MPI ranks 1228 1229 Options Database Keys: 1230 + -pc_gasm_total_subdomains <n> - Sets total number of local subdomains to be distributed among MPI ranks 1231 . -pc_gasm_view_subdomains - activates the printing of subdomain indices in `PCView()`, -ksp_view or -snes_view 1232 . -pc_gasm_print_subdomains - activates the printing of subdomain indices in `PCSetUp()` 1233 . -pc_gasm_overlap <ovl> - Sets overlap by which to (automatically) extend local subdomains 1234 - -pc_gasm_type [basic,restrict,interpolate,none] - Sets `PCGASMType` 1235 1236 Notes: 1237 To set options on the solvers for each block append -sub_ to all the `KSP`, and `PC` 1238 options database keys. For example, -sub_pc_type ilu -sub_pc_factor_levels 1 -sub_ksp_type preonly 1239 1240 To set the options on the solvers separate for each block call `PCGASMGetSubKSP()` 1241 and set the options directly on the resulting `KSP` object (you can access its `PC` 1242 with `KSPGetPC())` 1243 1244 Level: beginner 1245 1246 References: 1247 + * - M Dryja, OB Widlund, An additive variant of the Schwarz alternating method for the case of many subregions 1248 Courant Institute, New York University Technical report 1249 - * - Barry Smith, Petter Bjorstad, and William Gropp, Domain Decompositions: Parallel Multilevel Methods for Elliptic Partial Differential Equations, 1250 Cambridge University Press. 1251 1252 .seealso: `PCCreate()`, `PCSetType()`, `PCType`, `PC`, `PCASM`, `PCGASMType`, `PCGASMSetType()`, 1253 `PCBJACOBI`, `PCGASMGetSubKSP()`, `PCGASMSetSubdomains()`, 1254 `PCSetModifySubMatrices()`, `PCGASMSetOverlap()`, `PCGASMSetType()` 1255 M*/ 1256 1257 PETSC_EXTERN PetscErrorCode PCCreate_GASM(PC pc) 1258 { 1259 PC_GASM *osm; 1260 1261 PetscFunctionBegin; 1262 PetscCall(PetscNew(&osm)); 1263 1264 osm->N = PETSC_DETERMINE; 1265 osm->n = PETSC_DECIDE; 1266 osm->nmax = PETSC_DETERMINE; 1267 osm->overlap = 0; 1268 osm->ksp = NULL; 1269 osm->gorestriction = NULL; 1270 osm->girestriction = NULL; 1271 osm->pctoouter = NULL; 1272 osm->gx = NULL; 1273 osm->gy = NULL; 1274 osm->x = NULL; 1275 osm->y = NULL; 1276 osm->pcx = NULL; 1277 osm->pcy = NULL; 1278 osm->permutationIS = NULL; 1279 osm->permutationP = NULL; 1280 osm->pcmat = NULL; 1281 osm->ois = NULL; 1282 osm->iis = NULL; 1283 osm->pmat = NULL; 1284 osm->type = PC_GASM_RESTRICT; 1285 osm->same_subdomain_solvers = PETSC_TRUE; 1286 osm->sort_indices = PETSC_TRUE; 1287 osm->dm_subdomains = PETSC_FALSE; 1288 osm->hierarchicalpartitioning = PETSC_FALSE; 1289 1290 pc->data = (void *)osm; 1291 pc->ops->apply = PCApply_GASM; 1292 pc->ops->matapply = PCMatApply_GASM; 1293 pc->ops->applytranspose = PCApplyTranspose_GASM; 1294 pc->ops->setup = PCSetUp_GASM; 1295 pc->ops->reset = PCReset_GASM; 1296 pc->ops->destroy = PCDestroy_GASM; 1297 pc->ops->setfromoptions = PCSetFromOptions_GASM; 1298 pc->ops->setuponblocks = PCSetUpOnBlocks_GASM; 1299 pc->ops->view = PCView_GASM; 1300 pc->ops->applyrichardson = NULL; 1301 1302 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGASMSetSubdomains_C", PCGASMSetSubdomains_GASM)); 1303 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGASMSetOverlap_C", PCGASMSetOverlap_GASM)); 1304 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGASMSetType_C", PCGASMSetType_GASM)); 1305 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGASMSetSortIndices_C", PCGASMSetSortIndices_GASM)); 1306 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGASMGetSubKSP_C", PCGASMGetSubKSP_GASM)); 1307 PetscFunctionReturn(PETSC_SUCCESS); 1308 } 1309 1310 PetscErrorCode PCGASMCreateLocalSubdomains(Mat A, PetscInt nloc, IS *iis[]) 1311 { 1312 MatPartitioning mpart; 1313 const char *prefix; 1314 PetscInt i, j, rstart, rend, bs; 1315 PetscBool hasop, isbaij = PETSC_FALSE, foundpart = PETSC_FALSE; 1316 Mat Ad = NULL, adj; 1317 IS ispart, isnumb, *is; 1318 1319 PetscFunctionBegin; 1320 PetscCheck(nloc >= 1, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "number of local subdomains must > 0, got nloc = %" PetscInt_FMT, nloc); 1321 1322 /* Get prefix, row distribution, and block size */ 1323 PetscCall(MatGetOptionsPrefix(A, &prefix)); 1324 PetscCall(MatGetOwnershipRange(A, &rstart, &rend)); 1325 PetscCall(MatGetBlockSize(A, &bs)); 1326 PetscCheck(rstart / bs * bs == rstart && rend / bs * bs == rend, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "bad row distribution [%" PetscInt_FMT ",%" PetscInt_FMT ") for matrix block size %" PetscInt_FMT, rstart, rend, bs); 1327 1328 /* Get diagonal block from matrix if possible */ 1329 PetscCall(MatHasOperation(A, MATOP_GET_DIAGONAL_BLOCK, &hasop)); 1330 if (hasop) PetscCall(MatGetDiagonalBlock(A, &Ad)); 1331 if (Ad) { 1332 PetscCall(PetscObjectBaseTypeCompare((PetscObject)Ad, MATSEQBAIJ, &isbaij)); 1333 if (!isbaij) PetscCall(PetscObjectBaseTypeCompare((PetscObject)Ad, MATSEQSBAIJ, &isbaij)); 1334 } 1335 if (Ad && nloc > 1) { 1336 PetscBool match, done; 1337 /* Try to setup a good matrix partitioning if available */ 1338 PetscCall(MatPartitioningCreate(PETSC_COMM_SELF, &mpart)); 1339 PetscCall(PetscObjectSetOptionsPrefix((PetscObject)mpart, prefix)); 1340 PetscCall(MatPartitioningSetFromOptions(mpart)); 1341 PetscCall(PetscObjectTypeCompare((PetscObject)mpart, MATPARTITIONINGCURRENT, &match)); 1342 if (!match) PetscCall(PetscObjectTypeCompare((PetscObject)mpart, MATPARTITIONINGSQUARE, &match)); 1343 if (!match) { /* assume a "good" partitioner is available */ 1344 PetscInt na; 1345 const PetscInt *ia, *ja; 1346 PetscCall(MatGetRowIJ(Ad, 0, PETSC_TRUE, isbaij, &na, &ia, &ja, &done)); 1347 if (done) { 1348 /* Build adjacency matrix by hand. Unfortunately a call to 1349 MatConvert(Ad,MATMPIADJ,MAT_INITIAL_MATRIX,&adj) will 1350 remove the block-aij structure and we cannot expect 1351 MatPartitioning to split vertices as we need */ 1352 PetscInt i, j, len, nnz, cnt, *iia = NULL, *jja = NULL; 1353 const PetscInt *row; 1354 nnz = 0; 1355 for (i = 0; i < na; i++) { /* count number of nonzeros */ 1356 len = ia[i + 1] - ia[i]; 1357 row = ja + ia[i]; 1358 for (j = 0; j < len; j++) { 1359 if (row[j] == i) { /* don't count diagonal */ 1360 len--; 1361 break; 1362 } 1363 } 1364 nnz += len; 1365 } 1366 PetscCall(PetscMalloc1(na + 1, &iia)); 1367 PetscCall(PetscMalloc1(nnz, &jja)); 1368 nnz = 0; 1369 iia[0] = 0; 1370 for (i = 0; i < na; i++) { /* fill adjacency */ 1371 cnt = 0; 1372 len = ia[i + 1] - ia[i]; 1373 row = ja + ia[i]; 1374 for (j = 0; j < len; j++) { 1375 if (row[j] != i) jja[nnz + cnt++] = row[j]; /* if not diagonal */ 1376 } 1377 nnz += cnt; 1378 iia[i + 1] = nnz; 1379 } 1380 /* Partitioning of the adjacency matrix */ 1381 PetscCall(MatCreateMPIAdj(PETSC_COMM_SELF, na, na, iia, jja, NULL, &adj)); 1382 PetscCall(MatPartitioningSetAdjacency(mpart, adj)); 1383 PetscCall(MatPartitioningSetNParts(mpart, nloc)); 1384 PetscCall(MatPartitioningApply(mpart, &ispart)); 1385 PetscCall(ISPartitioningToNumbering(ispart, &isnumb)); 1386 PetscCall(MatDestroy(&adj)); 1387 foundpart = PETSC_TRUE; 1388 } 1389 PetscCall(MatRestoreRowIJ(Ad, 0, PETSC_TRUE, isbaij, &na, &ia, &ja, &done)); 1390 } 1391 PetscCall(MatPartitioningDestroy(&mpart)); 1392 } 1393 PetscCall(PetscMalloc1(nloc, &is)); 1394 if (!foundpart) { 1395 /* Partitioning by contiguous chunks of rows */ 1396 1397 PetscInt mbs = (rend - rstart) / bs; 1398 PetscInt start = rstart; 1399 for (i = 0; i < nloc; i++) { 1400 PetscInt count = (mbs / nloc + ((mbs % nloc) > i)) * bs; 1401 PetscCall(ISCreateStride(PETSC_COMM_SELF, count, start, 1, &is[i])); 1402 start += count; 1403 } 1404 1405 } else { 1406 /* Partitioning by adjacency of diagonal block */ 1407 1408 const PetscInt *numbering; 1409 PetscInt *count, nidx, *indices, *newidx, start = 0; 1410 /* Get node count in each partition */ 1411 PetscCall(PetscMalloc1(nloc, &count)); 1412 PetscCall(ISPartitioningCount(ispart, nloc, count)); 1413 if (isbaij && bs > 1) { /* adjust for the block-aij case */ 1414 for (i = 0; i < nloc; i++) count[i] *= bs; 1415 } 1416 /* Build indices from node numbering */ 1417 PetscCall(ISGetLocalSize(isnumb, &nidx)); 1418 PetscCall(PetscMalloc1(nidx, &indices)); 1419 for (i = 0; i < nidx; i++) indices[i] = i; /* needs to be initialized */ 1420 PetscCall(ISGetIndices(isnumb, &numbering)); 1421 PetscCall(PetscSortIntWithPermutation(nidx, numbering, indices)); 1422 PetscCall(ISRestoreIndices(isnumb, &numbering)); 1423 if (isbaij && bs > 1) { /* adjust for the block-aij case */ 1424 PetscCall(PetscMalloc1(nidx * bs, &newidx)); 1425 for (i = 0; i < nidx; i++) { 1426 for (j = 0; j < bs; j++) newidx[i * bs + j] = indices[i] * bs + j; 1427 } 1428 PetscCall(PetscFree(indices)); 1429 nidx *= bs; 1430 indices = newidx; 1431 } 1432 /* Shift to get global indices */ 1433 for (i = 0; i < nidx; i++) indices[i] += rstart; 1434 1435 /* Build the index sets for each block */ 1436 for (i = 0; i < nloc; i++) { 1437 PetscCall(ISCreateGeneral(PETSC_COMM_SELF, count[i], &indices[start], PETSC_COPY_VALUES, &is[i])); 1438 PetscCall(ISSort(is[i])); 1439 start += count[i]; 1440 } 1441 1442 PetscCall(PetscFree(count)); 1443 PetscCall(PetscFree(indices)); 1444 PetscCall(ISDestroy(&isnumb)); 1445 PetscCall(ISDestroy(&ispart)); 1446 } 1447 *iis = is; 1448 PetscFunctionReturn(PETSC_SUCCESS); 1449 } 1450 1451 PETSC_INTERN PetscErrorCode PCGASMCreateStraddlingSubdomains(Mat A, PetscInt N, PetscInt *n, IS *iis[]) 1452 { 1453 PetscFunctionBegin; 1454 PetscCall(MatSubdomainsCreateCoalesce(A, N, n, iis)); 1455 PetscFunctionReturn(PETSC_SUCCESS); 1456 } 1457 1458 /*@C 1459 PCGASMCreateSubdomains - Creates n index sets defining n nonoverlapping subdomains for the `PCGASM` additive 1460 Schwarz preconditioner for a any problem based on its matrix. 1461 1462 Collective 1463 1464 Input Parameters: 1465 + A - The global matrix operator 1466 - N - the number of global subdomains requested 1467 1468 Output Parameters: 1469 + n - the number of subdomains created on this MPI rank 1470 - iis - the array of index sets defining the local inner subdomains (on which the correction is applied) 1471 1472 Level: advanced 1473 1474 Note: 1475 When N >= A's communicator size, each subdomain is local -- contained within a single MPI rank. 1476 When N < size, the subdomains are 'straddling' (rank boundaries) and are no longer local. 1477 The resulting subdomains can be use in `PCGASMSetSubdomains`(pc,n,iss,NULL). The overlapping 1478 outer subdomains will be automatically generated from these according to the requested amount of 1479 overlap; this is currently supported only with local subdomains. 1480 1481 .seealso: `PCGASM`, `PCGASMSetSubdomains()`, `PCGASMDestroySubdomains()` 1482 @*/ 1483 PetscErrorCode PCGASMCreateSubdomains(Mat A, PetscInt N, PetscInt *n, IS *iis[]) 1484 { 1485 PetscMPIInt size; 1486 1487 PetscFunctionBegin; 1488 PetscValidHeaderSpecific(A, MAT_CLASSID, 1); 1489 PetscValidPointer(iis, 4); 1490 1491 PetscCheck(N >= 1, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Number of subdomains must be > 0, N = %" PetscInt_FMT, N); 1492 PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)A), &size)); 1493 if (N >= size) { 1494 *n = N / size + (N % size); 1495 PetscCall(PCGASMCreateLocalSubdomains(A, *n, iis)); 1496 } else { 1497 PetscCall(PCGASMCreateStraddlingSubdomains(A, N, n, iis)); 1498 } 1499 PetscFunctionReturn(PETSC_SUCCESS); 1500 } 1501 1502 /*@C 1503 PCGASMDestroySubdomains - Destroys the index sets created with 1504 `PCGASMCreateSubdomains()` or `PCGASMCreateSubdomains2D()`. Should be 1505 called after setting subdomains with `PCGASMSetSubdomains()`. 1506 1507 Collective 1508 1509 Input Parameters: 1510 + n - the number of index sets 1511 . iis - the array of inner subdomains, 1512 - ois - the array of outer subdomains, can be NULL 1513 1514 Level: intermediate 1515 1516 Note: 1517 This is a convenience subroutine that walks each list, 1518 destroys each `IS` on the list, and then frees the list. At the end the 1519 list pointers are set to NULL. 1520 1521 .seealso: `PCGASM`, `PCGASMCreateSubdomains()`, `PCGASMSetSubdomains()` 1522 @*/ 1523 PetscErrorCode PCGASMDestroySubdomains(PetscInt n, IS **iis, IS **ois) 1524 { 1525 PetscInt i; 1526 1527 PetscFunctionBegin; 1528 if (n <= 0) PetscFunctionReturn(PETSC_SUCCESS); 1529 if (ois) { 1530 PetscValidPointer(ois, 3); 1531 if (*ois) { 1532 PetscValidPointer(*ois, 3); 1533 for (i = 0; i < n; i++) PetscCall(ISDestroy(&(*ois)[i])); 1534 PetscCall(PetscFree((*ois))); 1535 } 1536 } 1537 if (iis) { 1538 PetscValidPointer(iis, 2); 1539 if (*iis) { 1540 PetscValidPointer(*iis, 2); 1541 for (i = 0; i < n; i++) PetscCall(ISDestroy(&(*iis)[i])); 1542 PetscCall(PetscFree((*iis))); 1543 } 1544 } 1545 PetscFunctionReturn(PETSC_SUCCESS); 1546 } 1547 1548 #define PCGASMLocalSubdomainBounds2D(M, N, xleft, ylow, xright, yhigh, first, last, xleft_loc, ylow_loc, xright_loc, yhigh_loc, n) \ 1549 { \ 1550 PetscInt first_row = first / M, last_row = last / M + 1; \ 1551 /* \ 1552 Compute ylow_loc and yhigh_loc so that (ylow_loc,xleft) and (yhigh_loc,xright) are the corners \ 1553 of the bounding box of the intersection of the subdomain with the local ownership range (local \ 1554 subdomain). \ 1555 Also compute xleft_loc and xright_loc as the lower and upper bounds on the first and last rows \ 1556 of the intersection. \ 1557 */ \ 1558 /* ylow_loc is the grid row containing the first element of the local sumbdomain */ \ 1559 *ylow_loc = PetscMax(first_row, ylow); \ 1560 /* xleft_loc is the offset of first element of the local subdomain within its grid row (might actually be outside the local subdomain) */ \ 1561 *xleft_loc = *ylow_loc == first_row ? PetscMax(first % M, xleft) : xleft; \ 1562 /* yhigh_loc is the grid row above the last local subdomain element */ \ 1563 *yhigh_loc = PetscMin(last_row, yhigh); \ 1564 /* xright is the offset of the end of the local subdomain within its grid row (might actually be outside the local subdomain) */ \ 1565 *xright_loc = *yhigh_loc == last_row ? PetscMin(xright, last % M) : xright; \ 1566 /* Now compute the size of the local subdomain n. */ \ 1567 *n = 0; \ 1568 if (*ylow_loc < *yhigh_loc) { \ 1569 PetscInt width = xright - xleft; \ 1570 *n += width * (*yhigh_loc - *ylow_loc - 1); \ 1571 *n += PetscMin(PetscMax(*xright_loc - xleft, 0), width); \ 1572 *n -= PetscMin(PetscMax(*xleft_loc - xleft, 0), width); \ 1573 } \ 1574 } 1575 1576 /*@ 1577 PCGASMCreateSubdomains2D - Creates the index sets for the `PCGASM` overlapping Schwarz 1578 preconditioner for a two-dimensional problem on a regular grid. 1579 1580 Collective 1581 1582 Input Parameters: 1583 + pc - the preconditioner context 1584 . M - the global number of grid points in the x direction 1585 . N - the global number of grid points in the y direction 1586 . Mdomains - the global number of subdomains in the x direction 1587 . Ndomains - the global number of subdomains in the y direction 1588 . dof - degrees of freedom per node 1589 - overlap - overlap in mesh lines 1590 1591 Output Parameters: 1592 + Nsub - the number of local subdomains created 1593 . iis - array of index sets defining inner (nonoverlapping) subdomains 1594 - ois - array of index sets defining outer (overlapping, if overlap > 0) subdomains 1595 1596 Level: advanced 1597 1598 .seealso: `PCGASM`, `PCGASMSetSubdomains()`, `PCGASMGetSubKSP()`, `PCGASMSetOverlap()` 1599 @*/ 1600 PetscErrorCode PCGASMCreateSubdomains2D(PC pc, PetscInt M, PetscInt N, PetscInt Mdomains, PetscInt Ndomains, PetscInt dof, PetscInt overlap, PetscInt *nsub, IS **iis, IS **ois) 1601 { 1602 PetscMPIInt size, rank; 1603 PetscInt i, j; 1604 PetscInt maxheight, maxwidth; 1605 PetscInt xstart, xleft, xright, xleft_loc, xright_loc; 1606 PetscInt ystart, ylow, yhigh, ylow_loc, yhigh_loc; 1607 PetscInt x[2][2], y[2][2], n[2]; 1608 PetscInt first, last; 1609 PetscInt nidx, *idx; 1610 PetscInt ii, jj, s, q, d; 1611 PetscInt k, kk; 1612 PetscMPIInt color; 1613 MPI_Comm comm, subcomm; 1614 IS **xis = NULL, **is = ois, **is_local = iis; 1615 1616 PetscFunctionBegin; 1617 PetscCall(PetscObjectGetComm((PetscObject)pc, &comm)); 1618 PetscCallMPI(MPI_Comm_size(comm, &size)); 1619 PetscCallMPI(MPI_Comm_rank(comm, &rank)); 1620 PetscCall(MatGetOwnershipRange(pc->pmat, &first, &last)); 1621 PetscCheck((first % dof) == 0 && (last % dof) == 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, 1622 "Matrix row partitioning unsuitable for domain decomposition: local row range (%" PetscInt_FMT ",%" PetscInt_FMT ") " 1623 "does not respect the number of degrees of freedom per grid point %" PetscInt_FMT, 1624 first, last, dof); 1625 1626 /* Determine the number of domains with nonzero intersections with the local ownership range. */ 1627 s = 0; 1628 ystart = 0; 1629 for (j = 0; j < Ndomains; ++j) { 1630 maxheight = N / Ndomains + ((N % Ndomains) > j); /* Maximal height of subdomain */ 1631 PetscCheck(maxheight >= 2, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Too many %" PetscInt_FMT " subdomains in the vertical direction for mesh height %" PetscInt_FMT, Ndomains, N); 1632 /* Vertical domain limits with an overlap. */ 1633 ylow = PetscMax(ystart - overlap, 0); 1634 yhigh = PetscMin(ystart + maxheight + overlap, N); 1635 xstart = 0; 1636 for (i = 0; i < Mdomains; ++i) { 1637 maxwidth = M / Mdomains + ((M % Mdomains) > i); /* Maximal width of subdomain */ 1638 PetscCheck(maxwidth >= 2, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Too many %" PetscInt_FMT " subdomains in the horizontal direction for mesh width %" PetscInt_FMT, Mdomains, M); 1639 /* Horizontal domain limits with an overlap. */ 1640 xleft = PetscMax(xstart - overlap, 0); 1641 xright = PetscMin(xstart + maxwidth + overlap, M); 1642 /* 1643 Determine whether this subdomain intersects this rank's ownership range of pc->pmat. 1644 */ 1645 PCGASMLocalSubdomainBounds2D(M, N, xleft, ylow, xright, yhigh, first, last, (&xleft_loc), (&ylow_loc), (&xright_loc), (&yhigh_loc), (&nidx)); 1646 if (nidx) ++s; 1647 xstart += maxwidth; 1648 } /* for (i = 0; i < Mdomains; ++i) */ 1649 ystart += maxheight; 1650 } /* for (j = 0; j < Ndomains; ++j) */ 1651 1652 /* Now we can allocate the necessary number of ISs. */ 1653 *nsub = s; 1654 PetscCall(PetscMalloc1(*nsub, is)); 1655 PetscCall(PetscMalloc1(*nsub, is_local)); 1656 s = 0; 1657 ystart = 0; 1658 for (j = 0; j < Ndomains; ++j) { 1659 maxheight = N / Ndomains + ((N % Ndomains) > j); /* Maximal height of subdomain */ 1660 PetscCheck(maxheight >= 2, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Too many %" PetscInt_FMT " subdomains in the vertical direction for mesh height %" PetscInt_FMT, Ndomains, N); 1661 /* Vertical domain limits with an overlap. */ 1662 y[0][0] = PetscMax(ystart - overlap, 0); 1663 y[0][1] = PetscMin(ystart + maxheight + overlap, N); 1664 /* Vertical domain limits without an overlap. */ 1665 y[1][0] = ystart; 1666 y[1][1] = PetscMin(ystart + maxheight, N); 1667 xstart = 0; 1668 for (i = 0; i < Mdomains; ++i) { 1669 maxwidth = M / Mdomains + ((M % Mdomains) > i); /* Maximal width of subdomain */ 1670 PetscCheck(maxwidth >= 2, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Too many %" PetscInt_FMT " subdomains in the horizontal direction for mesh width %" PetscInt_FMT, Mdomains, M); 1671 /* Horizontal domain limits with an overlap. */ 1672 x[0][0] = PetscMax(xstart - overlap, 0); 1673 x[0][1] = PetscMin(xstart + maxwidth + overlap, M); 1674 /* Horizontal domain limits without an overlap. */ 1675 x[1][0] = xstart; 1676 x[1][1] = PetscMin(xstart + maxwidth, M); 1677 /* 1678 Determine whether this domain intersects this rank's ownership range of pc->pmat. 1679 Do this twice: first for the domains with overlaps, and once without. 1680 During the first pass create the subcommunicators, and use them on the second pass as well. 1681 */ 1682 for (q = 0; q < 2; ++q) { 1683 PetscBool split = PETSC_FALSE; 1684 /* 1685 domain limits, (xleft, xright) and (ylow, yheigh) are adjusted 1686 according to whether the domain with an overlap or without is considered. 1687 */ 1688 xleft = x[q][0]; 1689 xright = x[q][1]; 1690 ylow = y[q][0]; 1691 yhigh = y[q][1]; 1692 PCGASMLocalSubdomainBounds2D(M, N, xleft, ylow, xright, yhigh, first, last, (&xleft_loc), (&ylow_loc), (&xright_loc), (&yhigh_loc), (&nidx)); 1693 nidx *= dof; 1694 n[q] = nidx; 1695 /* 1696 Based on the counted number of indices in the local domain *with an overlap*, 1697 construct a subcommunicator of all the MPI ranks supporting this domain. 1698 Observe that a domain with an overlap might have nontrivial local support, 1699 while the domain without an overlap might not. Hence, the decision to participate 1700 in the subcommunicator must be based on the domain with an overlap. 1701 */ 1702 if (q == 0) { 1703 if (nidx) color = 1; 1704 else color = MPI_UNDEFINED; 1705 PetscCallMPI(MPI_Comm_split(comm, color, rank, &subcomm)); 1706 split = PETSC_TRUE; 1707 } 1708 /* 1709 Proceed only if the number of local indices *with an overlap* is nonzero. 1710 */ 1711 if (n[0]) { 1712 if (q == 0) xis = is; 1713 if (q == 1) { 1714 /* 1715 The IS for the no-overlap subdomain shares a communicator with the overlapping domain. 1716 Moreover, if the overlap is zero, the two ISs are identical. 1717 */ 1718 if (overlap == 0) { 1719 (*is_local)[s] = (*is)[s]; 1720 PetscCall(PetscObjectReference((PetscObject)(*is)[s])); 1721 continue; 1722 } else { 1723 xis = is_local; 1724 subcomm = ((PetscObject)(*is)[s])->comm; 1725 } 1726 } /* if (q == 1) */ 1727 idx = NULL; 1728 PetscCall(PetscMalloc1(nidx, &idx)); 1729 if (nidx) { 1730 k = 0; 1731 for (jj = ylow_loc; jj < yhigh_loc; ++jj) { 1732 PetscInt x0 = (jj == ylow_loc) ? xleft_loc : xleft; 1733 PetscInt x1 = (jj == yhigh_loc - 1) ? xright_loc : xright; 1734 kk = dof * (M * jj + x0); 1735 for (ii = x0; ii < x1; ++ii) { 1736 for (d = 0; d < dof; ++d) idx[k++] = kk++; 1737 } 1738 } 1739 } 1740 PetscCall(ISCreateGeneral(subcomm, nidx, idx, PETSC_OWN_POINTER, (*xis) + s)); 1741 if (split) PetscCallMPI(MPI_Comm_free(&subcomm)); 1742 } /* if (n[0]) */ 1743 } /* for (q = 0; q < 2; ++q) */ 1744 if (n[0]) ++s; 1745 xstart += maxwidth; 1746 } /* for (i = 0; i < Mdomains; ++i) */ 1747 ystart += maxheight; 1748 } /* for (j = 0; j < Ndomains; ++j) */ 1749 PetscFunctionReturn(PETSC_SUCCESS); 1750 } 1751 1752 /*@C 1753 PCGASMGetSubdomains - Gets the subdomains supported on this MPI rank 1754 for the `PCGASM` additive Schwarz preconditioner. 1755 1756 Not Collective 1757 1758 Input Parameter: 1759 . pc - the preconditioner context 1760 1761 Output Parameters: 1762 + n - the number of subdomains for this MPI rank (default value = 1) 1763 . iis - the index sets that define the inner subdomains (without overlap) supported on this rank (can be NULL) 1764 - ois - the index sets that define the outer subdomains (with overlap) supported on this rank (can be NULL) 1765 1766 Note: 1767 The user is responsible for destroying the `IS`s and freeing the returned arrays. 1768 The `IS` numbering is in the parallel, global numbering of the vector. 1769 1770 Level: advanced 1771 1772 .seealso: `PCGASM`, `PCGASMSetOverlap()`, `PCGASMGetSubKSP()`, `PCGASMCreateSubdomains2D()`, 1773 `PCGASMSetSubdomains()`, `PCGASMGetSubmatrices()` 1774 @*/ 1775 PetscErrorCode PCGASMGetSubdomains(PC pc, PetscInt *n, IS *iis[], IS *ois[]) 1776 { 1777 PC_GASM *osm; 1778 PetscBool match; 1779 PetscInt i; 1780 1781 PetscFunctionBegin; 1782 PetscValidHeaderSpecific(pc, PC_CLASSID, 1); 1783 PetscCall(PetscObjectTypeCompare((PetscObject)pc, PCGASM, &match)); 1784 PetscCheck(match, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_WRONG, "Incorrect object type: expected %s, got %s instead", PCGASM, ((PetscObject)pc)->type_name); 1785 osm = (PC_GASM *)pc->data; 1786 if (n) *n = osm->n; 1787 if (iis) PetscCall(PetscMalloc1(osm->n, iis)); 1788 if (ois) PetscCall(PetscMalloc1(osm->n, ois)); 1789 if (iis || ois) { 1790 for (i = 0; i < osm->n; ++i) { 1791 if (iis) (*iis)[i] = osm->iis[i]; 1792 if (ois) (*ois)[i] = osm->ois[i]; 1793 } 1794 } 1795 PetscFunctionReturn(PETSC_SUCCESS); 1796 } 1797 1798 /*@C 1799 PCGASMGetSubmatrices - Gets the local submatrices (for this MPI rank 1800 only) for the `PCGASM` additive Schwarz preconditioner. 1801 1802 Not Collective 1803 1804 Input Parameter: 1805 . pc - the preconditioner context 1806 1807 Output Parameters: 1808 + n - the number of matrices for this MPI rank (default value = 1) 1809 - mat - the matrices 1810 1811 Note: 1812 Matrices returned by this routine have the same communicators as the index sets (IS) 1813 used to define subdomains in `PCGASMSetSubdomains()` 1814 1815 Level: advanced 1816 1817 .seealso: `PCGASM`, `PCGASMSetOverlap()`, `PCGASMGetSubKSP()`, 1818 `PCGASMCreateSubdomains2D()`, `PCGASMSetSubdomains()`, `PCGASMGetSubdomains()` 1819 @*/ 1820 PetscErrorCode PCGASMGetSubmatrices(PC pc, PetscInt *n, Mat *mat[]) 1821 { 1822 PC_GASM *osm; 1823 PetscBool match; 1824 1825 PetscFunctionBegin; 1826 PetscValidHeaderSpecific(pc, PC_CLASSID, 1); 1827 PetscValidIntPointer(n, 2); 1828 if (mat) PetscValidPointer(mat, 3); 1829 PetscCheck(pc->setupcalled, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_WRONGSTATE, "Must call after KSPSetUp() or PCSetUp()."); 1830 PetscCall(PetscObjectTypeCompare((PetscObject)pc, PCGASM, &match)); 1831 PetscCheck(match, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_WRONG, "Expected %s, got %s instead", PCGASM, ((PetscObject)pc)->type_name); 1832 osm = (PC_GASM *)pc->data; 1833 if (n) *n = osm->n; 1834 if (mat) *mat = osm->pmat; 1835 PetscFunctionReturn(PETSC_SUCCESS); 1836 } 1837 1838 /*@ 1839 PCGASMSetUseDMSubdomains - Indicates whether to use `DMCreateDomainDecomposition()` to define the subdomains, whenever possible for `PCGASM` 1840 1841 Logically Collective 1842 1843 Input Parameters: 1844 + pc - the preconditioner 1845 - flg - boolean indicating whether to use subdomains defined by the `DM` 1846 1847 Options Database Key: 1848 . -pc_gasm_dm_subdomains -pc_gasm_overlap -pc_gasm_total_subdomains 1849 1850 Level: intermediate 1851 1852 Note: 1853 `PCGASMSetSubdomains()`, `PCGASMSetTotalSubdomains()` or `PCGASMSetOverlap()` take precedence over `PCGASMSetUseDMSubdomains()`, 1854 so setting `PCGASMSetSubdomains()` with nontrivial subdomain ISs or any of `PCGASMSetTotalSubdomains()` and `PCGASMSetOverlap()` 1855 automatically turns the latter off. 1856 1857 .seealso: `PCGASM`, `PCGASMGetUseDMSubdomains()`, `PCGASMSetSubdomains()`, `PCGASMSetOverlap()` 1858 `PCGASMCreateSubdomains2D()` 1859 @*/ 1860 PetscErrorCode PCGASMSetUseDMSubdomains(PC pc, PetscBool flg) 1861 { 1862 PC_GASM *osm = (PC_GASM *)pc->data; 1863 PetscBool match; 1864 1865 PetscFunctionBegin; 1866 PetscValidHeaderSpecific(pc, PC_CLASSID, 1); 1867 PetscValidLogicalCollectiveBool(pc, flg, 2); 1868 PetscCheck(!pc->setupcalled, ((PetscObject)pc)->comm, PETSC_ERR_ARG_WRONGSTATE, "Not for a setup PC."); 1869 PetscCall(PetscObjectTypeCompare((PetscObject)pc, PCGASM, &match)); 1870 if (match) { 1871 if (!osm->user_subdomains && osm->N == PETSC_DETERMINE && osm->overlap < 0) osm->dm_subdomains = flg; 1872 } 1873 PetscFunctionReturn(PETSC_SUCCESS); 1874 } 1875 1876 /*@ 1877 PCGASMGetUseDMSubdomains - Returns flag indicating whether to use `DMCreateDomainDecomposition()` to define the subdomains, whenever possible with `PCGASM` 1878 1879 Not Collective 1880 1881 Input Parameter: 1882 . pc - the preconditioner 1883 1884 Output Parameter: 1885 . flg - boolean indicating whether to use subdomains defined by the `DM` 1886 1887 Level: intermediate 1888 1889 .seealso: `PCGASM`, `PCGASMSetUseDMSubdomains()`, `PCGASMSetOverlap()` 1890 `PCGASMCreateSubdomains2D()` 1891 @*/ 1892 PetscErrorCode PCGASMGetUseDMSubdomains(PC pc, PetscBool *flg) 1893 { 1894 PC_GASM *osm = (PC_GASM *)pc->data; 1895 PetscBool match; 1896 1897 PetscFunctionBegin; 1898 PetscValidHeaderSpecific(pc, PC_CLASSID, 1); 1899 PetscValidBoolPointer(flg, 2); 1900 PetscCall(PetscObjectTypeCompare((PetscObject)pc, PCGASM, &match)); 1901 if (match) { 1902 if (flg) *flg = osm->dm_subdomains; 1903 } 1904 PetscFunctionReturn(PETSC_SUCCESS); 1905 } 1906