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(PetscStrncpy(fname, "stdout", sizeof(fname))); 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 static 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(KSPSetNestLevel(ksp, pc->kspnestlevel)); 516 PetscCall(KSPSetErrorIfNotConverged(ksp, pc->erroriffailure)); 517 PetscCall(PetscObjectIncrementTabLevel((PetscObject)ksp, (PetscObject)pc, 1)); 518 PetscCall(KSPSetType(ksp, KSPPREONLY)); 519 PetscCall(KSPGetPC(ksp, &subpc)); /* Why do we need this here? */ 520 if (subdomain_dm) { 521 PetscCall(KSPSetDM(ksp, subdomain_dm[i])); 522 PetscCall(DMDestroy(subdomain_dm + i)); 523 } 524 PetscCall(PCGetOptionsPrefix(pc, &prefix)); 525 PetscCall(KSPSetOptionsPrefix(ksp, prefix)); 526 if (subdomain_names && subdomain_names[i]) { 527 PetscCall(PetscSNPrintf(subprefix, PETSC_MAX_PATH_LEN, "sub_%s_", subdomain_names[i])); 528 PetscCall(KSPAppendOptionsPrefix(ksp, subprefix)); 529 PetscCall(PetscFree(subdomain_names[i])); 530 } 531 PetscCall(KSPAppendOptionsPrefix(ksp, "sub_")); 532 osm->ksp[i] = ksp; 533 } 534 PetscCall(PetscFree(subdomain_dm)); 535 PetscCall(PetscFree(subdomain_names)); 536 scall = MAT_INITIAL_MATRIX; 537 } else { /* if (pc->setupcalled) */ 538 /* 539 Destroy the submatrices from the previous iteration 540 */ 541 if (pc->flag == DIFFERENT_NONZERO_PATTERN) { 542 PetscCall(MatDestroyMatrices(osm->n, &osm->pmat)); 543 scall = MAT_INITIAL_MATRIX; 544 } 545 if (osm->permutationIS) { 546 PetscCall(MatCreateSubMatrix(pc->pmat, osm->permutationIS, osm->permutationIS, scall, &osm->permutationP)); 547 PetscCall(PetscObjectReference((PetscObject)osm->permutationP)); 548 osm->pcmat = pc->pmat; 549 pc->pmat = osm->permutationP; 550 } 551 } 552 553 /* 554 Extract the submatrices. 555 */ 556 if (size > 1) { 557 PetscCall(MatCreateSubMatricesMPI(pc->pmat, osm->n, osm->ois, osm->ois, scall, &osm->pmat)); 558 } else { 559 PetscCall(MatCreateSubMatrices(pc->pmat, osm->n, osm->ois, osm->ois, scall, &osm->pmat)); 560 } 561 if (scall == MAT_INITIAL_MATRIX) { 562 PetscCall(PetscObjectGetOptionsPrefix((PetscObject)pc->pmat, &pprefix)); 563 for (i = 0; i < osm->n; i++) PetscCall(PetscObjectSetOptionsPrefix((PetscObject)osm->pmat[i], pprefix)); 564 } 565 566 /* Return control to the user so that the submatrices can be modified (e.g., to apply 567 different boundary conditions for the submatrices than for the global problem) */ 568 PetscCall(PCModifySubMatrices(pc, osm->n, osm->ois, osm->ois, osm->pmat, pc->modifysubmatricesP)); 569 570 /* 571 Loop over submatrices putting them into local ksps 572 */ 573 for (i = 0; i < osm->n; i++) { 574 PetscCall(KSPSetOperators(osm->ksp[i], osm->pmat[i], osm->pmat[i])); 575 PetscCall(KSPGetOptionsPrefix(osm->ksp[i], &prefix)); 576 PetscCall(MatSetOptionsPrefix(osm->pmat[i], prefix)); 577 if (!pc->setupcalled) PetscCall(KSPSetFromOptions(osm->ksp[i])); 578 } 579 if (osm->pcmat) { 580 PetscCall(MatDestroy(&pc->pmat)); 581 pc->pmat = osm->pcmat; 582 osm->pcmat = NULL; 583 } 584 PetscFunctionReturn(PETSC_SUCCESS); 585 } 586 587 static PetscErrorCode PCSetUpOnBlocks_GASM(PC pc) 588 { 589 PC_GASM *osm = (PC_GASM *)pc->data; 590 PetscInt i; 591 592 PetscFunctionBegin; 593 for (i = 0; i < osm->n; i++) PetscCall(KSPSetUp(osm->ksp[i])); 594 PetscFunctionReturn(PETSC_SUCCESS); 595 } 596 597 static PetscErrorCode PCApply_GASM(PC pc, Vec xin, Vec yout) 598 { 599 PC_GASM *osm = (PC_GASM *)pc->data; 600 PetscInt i; 601 Vec x, y; 602 ScatterMode forward = SCATTER_FORWARD, reverse = SCATTER_REVERSE; 603 604 PetscFunctionBegin; 605 if (osm->pctoouter) { 606 PetscCall(VecScatterBegin(osm->pctoouter, xin, osm->pcx, INSERT_VALUES, SCATTER_REVERSE)); 607 PetscCall(VecScatterEnd(osm->pctoouter, xin, osm->pcx, INSERT_VALUES, SCATTER_REVERSE)); 608 x = osm->pcx; 609 y = osm->pcy; 610 } else { 611 x = xin; 612 y = yout; 613 } 614 /* 615 support for limiting the restriction or interpolation only to the inner 616 subdomain values (leaving the other values 0). 617 */ 618 if (!(osm->type & PC_GASM_RESTRICT)) { 619 /* have to zero the work RHS since scatter may leave some slots empty */ 620 PetscCall(VecZeroEntries(osm->gx)); 621 PetscCall(VecScatterBegin(osm->girestriction, x, osm->gx, INSERT_VALUES, forward)); 622 } else { 623 PetscCall(VecScatterBegin(osm->gorestriction, x, osm->gx, INSERT_VALUES, forward)); 624 } 625 PetscCall(VecZeroEntries(osm->gy)); 626 if (!(osm->type & PC_GASM_RESTRICT)) { 627 PetscCall(VecScatterEnd(osm->girestriction, x, osm->gx, INSERT_VALUES, forward)); 628 } else { 629 PetscCall(VecScatterEnd(osm->gorestriction, x, osm->gx, INSERT_VALUES, forward)); 630 } 631 /* do the subdomain solves */ 632 for (i = 0; i < osm->n; ++i) { 633 PetscCall(KSPSolve(osm->ksp[i], osm->x[i], osm->y[i])); 634 PetscCall(KSPCheckSolve(osm->ksp[i], pc, osm->y[i])); 635 } 636 /* do we need to zero y? */ 637 PetscCall(VecZeroEntries(y)); 638 if (!(osm->type & PC_GASM_INTERPOLATE)) { 639 PetscCall(VecScatterBegin(osm->girestriction, osm->gy, y, ADD_VALUES, reverse)); 640 PetscCall(VecScatterEnd(osm->girestriction, osm->gy, y, ADD_VALUES, reverse)); 641 } else { 642 PetscCall(VecScatterBegin(osm->gorestriction, osm->gy, y, ADD_VALUES, reverse)); 643 PetscCall(VecScatterEnd(osm->gorestriction, osm->gy, y, ADD_VALUES, reverse)); 644 } 645 if (osm->pctoouter) { 646 PetscCall(VecScatterBegin(osm->pctoouter, y, yout, INSERT_VALUES, SCATTER_FORWARD)); 647 PetscCall(VecScatterEnd(osm->pctoouter, y, yout, INSERT_VALUES, SCATTER_FORWARD)); 648 } 649 PetscFunctionReturn(PETSC_SUCCESS); 650 } 651 652 static PetscErrorCode PCMatApply_GASM(PC pc, Mat Xin, Mat Yout) 653 { 654 PC_GASM *osm = (PC_GASM *)pc->data; 655 Mat X, Y, O = NULL, Z, W; 656 Vec x, y; 657 PetscInt i, m, M, N; 658 ScatterMode forward = SCATTER_FORWARD, reverse = SCATTER_REVERSE; 659 660 PetscFunctionBegin; 661 PetscCheck(osm->n == 1, PetscObjectComm((PetscObject)pc), PETSC_ERR_SUP, "Not yet implemented"); 662 PetscCall(MatGetSize(Xin, NULL, &N)); 663 if (osm->pctoouter) { 664 PetscCall(VecGetLocalSize(osm->pcx, &m)); 665 PetscCall(VecGetSize(osm->pcx, &M)); 666 PetscCall(MatCreateDense(PetscObjectComm((PetscObject)osm->ois[0]), m, PETSC_DECIDE, M, N, NULL, &O)); 667 for (i = 0; i < N; ++i) { 668 PetscCall(MatDenseGetColumnVecRead(Xin, i, &x)); 669 PetscCall(MatDenseGetColumnVecWrite(O, i, &y)); 670 PetscCall(VecScatterBegin(osm->pctoouter, x, y, INSERT_VALUES, SCATTER_REVERSE)); 671 PetscCall(VecScatterEnd(osm->pctoouter, x, y, INSERT_VALUES, SCATTER_REVERSE)); 672 PetscCall(MatDenseRestoreColumnVecWrite(O, i, &y)); 673 PetscCall(MatDenseRestoreColumnVecRead(Xin, i, &x)); 674 } 675 X = Y = O; 676 } else { 677 X = Xin; 678 Y = Yout; 679 } 680 /* 681 support for limiting the restriction or interpolation only to the inner 682 subdomain values (leaving the other values 0). 683 */ 684 PetscCall(VecGetLocalSize(osm->x[0], &m)); 685 PetscCall(VecGetSize(osm->x[0], &M)); 686 PetscCall(MatCreateDense(PetscObjectComm((PetscObject)osm->ois[0]), m, PETSC_DECIDE, M, N, NULL, &Z)); 687 for (i = 0; i < N; ++i) { 688 PetscCall(MatDenseGetColumnVecRead(X, i, &x)); 689 PetscCall(MatDenseGetColumnVecWrite(Z, i, &y)); 690 if (!(osm->type & PC_GASM_RESTRICT)) { 691 /* have to zero the work RHS since scatter may leave some slots empty */ 692 PetscCall(VecZeroEntries(y)); 693 PetscCall(VecScatterBegin(osm->girestriction, x, y, INSERT_VALUES, forward)); 694 PetscCall(VecScatterEnd(osm->girestriction, x, y, INSERT_VALUES, forward)); 695 } else { 696 PetscCall(VecScatterBegin(osm->gorestriction, x, y, INSERT_VALUES, forward)); 697 PetscCall(VecScatterEnd(osm->gorestriction, x, y, INSERT_VALUES, forward)); 698 } 699 PetscCall(MatDenseRestoreColumnVecWrite(Z, i, &y)); 700 PetscCall(MatDenseRestoreColumnVecRead(X, i, &x)); 701 } 702 PetscCall(MatCreateDense(PetscObjectComm((PetscObject)osm->ois[0]), m, PETSC_DECIDE, M, N, NULL, &W)); 703 PetscCall(MatSetOption(Z, MAT_NO_OFF_PROC_ENTRIES, PETSC_TRUE)); 704 PetscCall(MatAssemblyBegin(Z, MAT_FINAL_ASSEMBLY)); 705 PetscCall(MatAssemblyEnd(Z, MAT_FINAL_ASSEMBLY)); 706 /* do the subdomain solve */ 707 PetscCall(KSPMatSolve(osm->ksp[0], Z, W)); 708 PetscCall(KSPCheckSolve(osm->ksp[0], pc, NULL)); 709 PetscCall(MatDestroy(&Z)); 710 /* do we need to zero y? */ 711 PetscCall(MatZeroEntries(Y)); 712 for (i = 0; i < N; ++i) { 713 PetscCall(MatDenseGetColumnVecWrite(Y, i, &y)); 714 PetscCall(MatDenseGetColumnVecRead(W, i, &x)); 715 if (!(osm->type & PC_GASM_INTERPOLATE)) { 716 PetscCall(VecScatterBegin(osm->girestriction, x, y, ADD_VALUES, reverse)); 717 PetscCall(VecScatterEnd(osm->girestriction, x, y, ADD_VALUES, reverse)); 718 } else { 719 PetscCall(VecScatterBegin(osm->gorestriction, x, y, ADD_VALUES, reverse)); 720 PetscCall(VecScatterEnd(osm->gorestriction, x, y, ADD_VALUES, reverse)); 721 } 722 PetscCall(MatDenseRestoreColumnVecRead(W, i, &x)); 723 if (osm->pctoouter) { 724 PetscCall(MatDenseGetColumnVecWrite(Yout, i, &x)); 725 PetscCall(VecScatterBegin(osm->pctoouter, y, x, INSERT_VALUES, SCATTER_FORWARD)); 726 PetscCall(VecScatterEnd(osm->pctoouter, y, x, INSERT_VALUES, SCATTER_FORWARD)); 727 PetscCall(MatDenseRestoreColumnVecRead(Yout, i, &x)); 728 } 729 PetscCall(MatDenseRestoreColumnVecWrite(Y, i, &y)); 730 } 731 PetscCall(MatDestroy(&W)); 732 PetscCall(MatDestroy(&O)); 733 PetscFunctionReturn(PETSC_SUCCESS); 734 } 735 736 static PetscErrorCode PCApplyTranspose_GASM(PC pc, Vec xin, Vec yout) 737 { 738 PC_GASM *osm = (PC_GASM *)pc->data; 739 PetscInt i; 740 Vec x, y; 741 ScatterMode forward = SCATTER_FORWARD, reverse = SCATTER_REVERSE; 742 743 PetscFunctionBegin; 744 if (osm->pctoouter) { 745 PetscCall(VecScatterBegin(osm->pctoouter, xin, osm->pcx, INSERT_VALUES, SCATTER_REVERSE)); 746 PetscCall(VecScatterEnd(osm->pctoouter, xin, osm->pcx, INSERT_VALUES, SCATTER_REVERSE)); 747 x = osm->pcx; 748 y = osm->pcy; 749 } else { 750 x = xin; 751 y = yout; 752 } 753 /* 754 Support for limiting the restriction or interpolation to only local 755 subdomain values (leaving the other values 0). 756 757 Note: these are reversed from the PCApply_GASM() because we are applying the 758 transpose of the three terms 759 */ 760 if (!(osm->type & PC_GASM_INTERPOLATE)) { 761 /* have to zero the work RHS since scatter may leave some slots empty */ 762 PetscCall(VecZeroEntries(osm->gx)); 763 PetscCall(VecScatterBegin(osm->girestriction, x, osm->gx, INSERT_VALUES, forward)); 764 } else { 765 PetscCall(VecScatterBegin(osm->gorestriction, x, osm->gx, INSERT_VALUES, forward)); 766 } 767 PetscCall(VecZeroEntries(osm->gy)); 768 if (!(osm->type & PC_GASM_INTERPOLATE)) { 769 PetscCall(VecScatterEnd(osm->girestriction, x, osm->gx, INSERT_VALUES, forward)); 770 } else { 771 PetscCall(VecScatterEnd(osm->gorestriction, x, osm->gx, INSERT_VALUES, forward)); 772 } 773 /* do the local solves */ 774 for (i = 0; i < osm->n; ++i) { /* Note that the solves are local, so we can go to osm->n, rather than osm->nmax. */ 775 PetscCall(KSPSolveTranspose(osm->ksp[i], osm->x[i], osm->y[i])); 776 PetscCall(KSPCheckSolve(osm->ksp[i], pc, osm->y[i])); 777 } 778 PetscCall(VecZeroEntries(y)); 779 if (!(osm->type & PC_GASM_RESTRICT)) { 780 PetscCall(VecScatterBegin(osm->girestriction, osm->gy, y, ADD_VALUES, reverse)); 781 PetscCall(VecScatterEnd(osm->girestriction, osm->gy, y, ADD_VALUES, reverse)); 782 } else { 783 PetscCall(VecScatterBegin(osm->gorestriction, osm->gy, y, ADD_VALUES, reverse)); 784 PetscCall(VecScatterEnd(osm->gorestriction, osm->gy, y, ADD_VALUES, reverse)); 785 } 786 if (osm->pctoouter) { 787 PetscCall(VecScatterBegin(osm->pctoouter, y, yout, INSERT_VALUES, SCATTER_FORWARD)); 788 PetscCall(VecScatterEnd(osm->pctoouter, y, yout, INSERT_VALUES, SCATTER_FORWARD)); 789 } 790 PetscFunctionReturn(PETSC_SUCCESS); 791 } 792 793 static PetscErrorCode PCReset_GASM(PC pc) 794 { 795 PC_GASM *osm = (PC_GASM *)pc->data; 796 PetscInt i; 797 798 PetscFunctionBegin; 799 if (osm->ksp) { 800 for (i = 0; i < osm->n; i++) PetscCall(KSPReset(osm->ksp[i])); 801 } 802 if (osm->pmat) { 803 if (osm->n > 0) { 804 PetscMPIInt size; 805 PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)pc), &size)); 806 if (size > 1) { 807 /* osm->pmat is created by MatCreateSubMatricesMPI(), cannot use MatDestroySubMatrices() */ 808 PetscCall(MatDestroyMatrices(osm->n, &osm->pmat)); 809 } else { 810 PetscCall(MatDestroySubMatrices(osm->n, &osm->pmat)); 811 } 812 } 813 } 814 if (osm->x) { 815 for (i = 0; i < osm->n; i++) { 816 PetscCall(VecDestroy(&osm->x[i])); 817 PetscCall(VecDestroy(&osm->y[i])); 818 } 819 } 820 PetscCall(VecDestroy(&osm->gx)); 821 PetscCall(VecDestroy(&osm->gy)); 822 823 PetscCall(VecScatterDestroy(&osm->gorestriction)); 824 PetscCall(VecScatterDestroy(&osm->girestriction)); 825 if (!osm->user_subdomains) { 826 PetscCall(PCGASMDestroySubdomains(osm->n, &osm->ois, &osm->iis)); 827 osm->N = PETSC_DETERMINE; 828 osm->nmax = PETSC_DETERMINE; 829 } 830 if (osm->pctoouter) PetscCall(VecScatterDestroy(&(osm->pctoouter))); 831 if (osm->permutationIS) PetscCall(ISDestroy(&(osm->permutationIS))); 832 if (osm->pcx) PetscCall(VecDestroy(&(osm->pcx))); 833 if (osm->pcy) PetscCall(VecDestroy(&(osm->pcy))); 834 if (osm->permutationP) PetscCall(MatDestroy(&(osm->permutationP))); 835 if (osm->pcmat) PetscCall(MatDestroy(&osm->pcmat)); 836 PetscFunctionReturn(PETSC_SUCCESS); 837 } 838 839 static PetscErrorCode PCDestroy_GASM(PC pc) 840 { 841 PC_GASM *osm = (PC_GASM *)pc->data; 842 PetscInt i; 843 844 PetscFunctionBegin; 845 PetscCall(PCReset_GASM(pc)); 846 /* PCReset will not destroy subdomains, if user_subdomains is true. */ 847 PetscCall(PCGASMDestroySubdomains(osm->n, &osm->ois, &osm->iis)); 848 if (osm->ksp) { 849 for (i = 0; i < osm->n; i++) PetscCall(KSPDestroy(&osm->ksp[i])); 850 PetscCall(PetscFree(osm->ksp)); 851 } 852 PetscCall(PetscFree(osm->x)); 853 PetscCall(PetscFree(osm->y)); 854 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGASMSetSubdomains_C", NULL)); 855 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGASMSetOverlap_C", NULL)); 856 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGASMSetType_C", NULL)); 857 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGASMSetSortIndices_C", NULL)); 858 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGASMGetSubKSP_C", NULL)); 859 PetscCall(PetscFree(pc->data)); 860 PetscFunctionReturn(PETSC_SUCCESS); 861 } 862 863 static PetscErrorCode PCSetFromOptions_GASM(PC pc, PetscOptionItems *PetscOptionsObject) 864 { 865 PC_GASM *osm = (PC_GASM *)pc->data; 866 PetscInt blocks, ovl; 867 PetscBool flg; 868 PCGASMType gasmtype; 869 870 PetscFunctionBegin; 871 PetscOptionsHeadBegin(PetscOptionsObject, "Generalized additive Schwarz options"); 872 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)); 873 PetscCall(PetscOptionsInt("-pc_gasm_total_subdomains", "Total number of subdomains across communicator", "PCGASMSetTotalSubdomains", osm->N, &blocks, &flg)); 874 if (flg) PetscCall(PCGASMSetTotalSubdomains(pc, blocks)); 875 PetscCall(PetscOptionsInt("-pc_gasm_overlap", "Number of overlapping degrees of freedom", "PCGASMSetOverlap", osm->overlap, &ovl, &flg)); 876 if (flg) { 877 PetscCall(PCGASMSetOverlap(pc, ovl)); 878 osm->dm_subdomains = PETSC_FALSE; 879 } 880 flg = PETSC_FALSE; 881 PetscCall(PetscOptionsEnum("-pc_gasm_type", "Type of restriction/extension", "PCGASMSetType", PCGASMTypes, (PetscEnum)osm->type, (PetscEnum *)&gasmtype, &flg)); 882 if (flg) PetscCall(PCGASMSetType(pc, gasmtype)); 883 PetscCall(PetscOptionsBool("-pc_gasm_use_hierachical_partitioning", "use hierarchical partitioning", NULL, osm->hierarchicalpartitioning, &osm->hierarchicalpartitioning, &flg)); 884 PetscOptionsHeadEnd(); 885 PetscFunctionReturn(PETSC_SUCCESS); 886 } 887 888 /*@ 889 PCGASMSetTotalSubdomains - sets the total number of subdomains to use across the communicator for `PCGASM` 890 891 Logically Collective 892 893 Input Parameters: 894 + pc - the preconditioner 895 - N - total number of subdomains 896 897 Level: beginner 898 899 .seealso: `PCGASM`, `PCGASMSetSubdomains()`, `PCGASMSetOverlap()` 900 `PCGASMCreateSubdomains2D()` 901 @*/ 902 PetscErrorCode PCGASMSetTotalSubdomains(PC pc, PetscInt N) 903 { 904 PC_GASM *osm = (PC_GASM *)pc->data; 905 PetscMPIInt size, rank; 906 907 PetscFunctionBegin; 908 PetscCheck(N >= 1, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Total number of subdomains must be 1 or more, got N = %" PetscInt_FMT, N); 909 PetscCheck(!pc->setupcalled, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_WRONGSTATE, "PCGASMSetTotalSubdomains() should be called before calling PCSetUp()."); 910 911 PetscCall(PCGASMDestroySubdomains(osm->n, &osm->iis, &osm->ois)); 912 osm->ois = osm->iis = NULL; 913 914 PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)pc), &size)); 915 PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)pc), &rank)); 916 osm->N = N; 917 osm->n = PETSC_DETERMINE; 918 osm->nmax = PETSC_DETERMINE; 919 osm->dm_subdomains = PETSC_FALSE; 920 PetscFunctionReturn(PETSC_SUCCESS); 921 } 922 923 static PetscErrorCode PCGASMSetSubdomains_GASM(PC pc, PetscInt n, IS iis[], IS ois[]) 924 { 925 PC_GASM *osm = (PC_GASM *)pc->data; 926 PetscInt i; 927 928 PetscFunctionBegin; 929 PetscCheck(n >= 1, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Each MPI rank must have 1 or more subdomains, got n = %" PetscInt_FMT, n); 930 PetscCheck(!pc->setupcalled, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_WRONGSTATE, "PCGASMSetSubdomains() should be called before calling PCSetUp()."); 931 932 PetscCall(PCGASMDestroySubdomains(osm->n, &osm->iis, &osm->ois)); 933 osm->iis = osm->ois = NULL; 934 osm->n = n; 935 osm->N = PETSC_DETERMINE; 936 osm->nmax = PETSC_DETERMINE; 937 if (ois) { 938 PetscCall(PetscMalloc1(n, &osm->ois)); 939 for (i = 0; i < n; i++) { 940 PetscCall(PetscObjectReference((PetscObject)ois[i])); 941 osm->ois[i] = ois[i]; 942 } 943 /* 944 Since the user set the outer subdomains, even if nontrivial overlap was requested via PCGASMSetOverlap(), 945 it will be ignored. To avoid confusion later on (e.g., when viewing the PC), the overlap size is set to -1. 946 */ 947 osm->overlap = -1; 948 /* inner subdomains must be provided */ 949 PetscCheck(iis, PETSC_COMM_SELF, PETSC_ERR_ARG_NULL, "inner indices have to be provided "); 950 } /* end if */ 951 if (iis) { 952 PetscCall(PetscMalloc1(n, &osm->iis)); 953 for (i = 0; i < n; i++) { 954 PetscCall(PetscObjectReference((PetscObject)iis[i])); 955 osm->iis[i] = iis[i]; 956 } 957 if (!ois) { 958 osm->ois = NULL; 959 /* if user does not provide outer indices, we will create the corresponding outer indices using osm->overlap =1 in PCSetUp_GASM */ 960 } 961 } 962 if (PetscDefined(USE_DEBUG)) { 963 PetscInt j, rstart, rend, *covered, lsize; 964 const PetscInt *indices; 965 /* check if the inner indices cover and only cover the local portion of the preconditioning matrix */ 966 PetscCall(MatGetOwnershipRange(pc->pmat, &rstart, &rend)); 967 PetscCall(PetscCalloc1(rend - rstart, &covered)); 968 /* check if the current MPI rank owns indices from others */ 969 for (i = 0; i < n; i++) { 970 PetscCall(ISGetIndices(osm->iis[i], &indices)); 971 PetscCall(ISGetLocalSize(osm->iis[i], &lsize)); 972 for (j = 0; j < lsize; j++) { 973 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]); 974 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]); 975 covered[indices[j] - rstart] = 1; 976 } 977 PetscCall(ISRestoreIndices(osm->iis[i], &indices)); 978 } 979 /* check if we miss any indices */ 980 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); 981 PetscCall(PetscFree(covered)); 982 } 983 if (iis) osm->user_subdomains = PETSC_TRUE; 984 PetscFunctionReturn(PETSC_SUCCESS); 985 } 986 987 static PetscErrorCode PCGASMSetOverlap_GASM(PC pc, PetscInt ovl) 988 { 989 PC_GASM *osm = (PC_GASM *)pc->data; 990 991 PetscFunctionBegin; 992 PetscCheck(ovl >= 0, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_OUTOFRANGE, "Negative overlap value requested"); 993 PetscCheck(!pc->setupcalled || ovl == osm->overlap, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_WRONGSTATE, "PCGASMSetOverlap() should be called before PCSetUp()."); 994 if (!pc->setupcalled) osm->overlap = ovl; 995 PetscFunctionReturn(PETSC_SUCCESS); 996 } 997 998 static PetscErrorCode PCGASMSetType_GASM(PC pc, PCGASMType type) 999 { 1000 PC_GASM *osm = (PC_GASM *)pc->data; 1001 1002 PetscFunctionBegin; 1003 osm->type = type; 1004 osm->type_set = PETSC_TRUE; 1005 PetscFunctionReturn(PETSC_SUCCESS); 1006 } 1007 1008 static PetscErrorCode PCGASMSetSortIndices_GASM(PC pc, PetscBool doSort) 1009 { 1010 PC_GASM *osm = (PC_GASM *)pc->data; 1011 1012 PetscFunctionBegin; 1013 osm->sort_indices = doSort; 1014 PetscFunctionReturn(PETSC_SUCCESS); 1015 } 1016 1017 /* 1018 FIXME: This routine might need to be modified now that multiple ranks per subdomain are allowed. 1019 In particular, it would upset the global subdomain number calculation. 1020 */ 1021 static PetscErrorCode PCGASMGetSubKSP_GASM(PC pc, PetscInt *n, PetscInt *first, KSP **ksp) 1022 { 1023 PC_GASM *osm = (PC_GASM *)pc->data; 1024 1025 PetscFunctionBegin; 1026 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"); 1027 1028 if (n) *n = osm->n; 1029 if (first) { 1030 PetscCallMPI(MPI_Scan(&osm->n, first, 1, MPIU_INT, MPI_SUM, PetscObjectComm((PetscObject)pc))); 1031 *first -= osm->n; 1032 } 1033 if (ksp) { 1034 /* Assume that local solves are now different; not necessarily 1035 true, though! This flag is used only for PCView_GASM() */ 1036 *ksp = osm->ksp; 1037 osm->same_subdomain_solvers = PETSC_FALSE; 1038 } 1039 PetscFunctionReturn(PETSC_SUCCESS); 1040 } /* PCGASMGetSubKSP_GASM() */ 1041 1042 /*@C 1043 PCGASMSetSubdomains - Sets the subdomains for this MPI rank 1044 for the additive Schwarz preconditioner with multiple MPI ranks per subdomain, `PCGASM` 1045 1046 Collective 1047 1048 Input Parameters: 1049 + pc - the preconditioner object 1050 . n - the number of subdomains for this MPI rank 1051 . iis - the index sets that define the inner subdomains (or `NULL` for PETSc to determine subdomains) 1052 - ois - the index sets that define the outer subdomains (or `NULL` to use the same as `iis`, or to construct by expanding `iis` by 1053 the requested overlap) 1054 1055 Level: advanced 1056 1057 Notes: 1058 The `IS` indices use the parallel, global numbering of the vector entries. 1059 1060 Inner subdomains are those where the correction is applied. 1061 1062 Outer subdomains are those where the residual necessary to obtain the 1063 corrections is obtained (see `PCGASMType` for the use of inner/outer subdomains). 1064 1065 Both inner and outer subdomains can extend over several MPI ranks. 1066 This rank's portion of a subdomain is known as a local subdomain. 1067 1068 Inner subdomains can not overlap with each other, do not have any entities from remote ranks, 1069 and have to cover the entire local subdomain owned by the current rank. The index sets on each 1070 rank should be ordered such that the ith local subdomain is connected to the ith remote subdomain 1071 on another MPI rank. 1072 1073 By default the `PGASM` preconditioner uses 1 (local) subdomain per MPI rank. 1074 1075 The `iis` and `ois` arrays may be freed after this call using `PCGASMDestroySubdomains()` 1076 1077 .seealso: `PCGASM`, `PCGASMSetOverlap()`, `PCGASMGetSubKSP()`, `PCGASMDestroySubdomains()`, 1078 `PCGASMCreateSubdomains2D()`, `PCGASMGetSubdomains()` 1079 @*/ 1080 PetscErrorCode PCGASMSetSubdomains(PC pc, PetscInt n, IS iis[], IS ois[]) 1081 { 1082 PC_GASM *osm = (PC_GASM *)pc->data; 1083 1084 PetscFunctionBegin; 1085 PetscValidHeaderSpecific(pc, PC_CLASSID, 1); 1086 PetscTryMethod(pc, "PCGASMSetSubdomains_C", (PC, PetscInt, IS[], IS[]), (pc, n, iis, ois)); 1087 osm->dm_subdomains = PETSC_FALSE; 1088 PetscFunctionReturn(PETSC_SUCCESS); 1089 } 1090 1091 /*@ 1092 PCGASMSetOverlap - Sets the overlap between a pair of subdomains for the 1093 additive Schwarz preconditioner `PCGASM`. Either all or no MPI ranks in the 1094 pc communicator must call this routine. 1095 1096 Logically Collective 1097 1098 Input Parameters: 1099 + pc - the preconditioner context 1100 - ovl - the amount of overlap between subdomains (ovl >= 0, default value = 0) 1101 1102 Options Database Key: 1103 . -pc_gasm_overlap <overlap> - Sets overlap 1104 1105 Level: intermediate 1106 1107 Notes: 1108 By default the `PCGASM` preconditioner uses 1 subdomain per rank. To use 1109 multiple subdomain per perocessor or "straddling" subdomains that intersect 1110 multiple ranks use `PCGASMSetSubdomains()` (or option `-pc_gasm_total_subdomains` <n>). 1111 1112 The overlap defaults to 0, so if one desires that no additional 1113 overlap be computed beyond what may have been set with a call to 1114 `PCGASMSetSubdomains()`, then `ovl` must be set to be 0. In particular, if one does 1115 not explicitly set the subdomains in application code, then all overlap would be computed 1116 internally by PETSc, and using an overlap of 0 would result in an `PCGASM` 1117 variant that is equivalent to the block Jacobi preconditioner. 1118 1119 One can define initial index sets with any overlap via 1120 `PCGASMSetSubdomains()`; the routine `PCGASMSetOverlap()` merely allows 1121 PETSc to extend that overlap further, if desired. 1122 1123 .seealso: `PCGASM`, `PCGASMSetSubdomains()`, `PCGASMGetSubKSP()`, 1124 `PCGASMCreateSubdomains2D()`, `PCGASMGetSubdomains()` 1125 @*/ 1126 PetscErrorCode PCGASMSetOverlap(PC pc, PetscInt ovl) 1127 { 1128 PC_GASM *osm = (PC_GASM *)pc->data; 1129 1130 PetscFunctionBegin; 1131 PetscValidHeaderSpecific(pc, PC_CLASSID, 1); 1132 PetscValidLogicalCollectiveInt(pc, ovl, 2); 1133 PetscTryMethod(pc, "PCGASMSetOverlap_C", (PC, PetscInt), (pc, ovl)); 1134 osm->dm_subdomains = PETSC_FALSE; 1135 PetscFunctionReturn(PETSC_SUCCESS); 1136 } 1137 1138 /*@ 1139 PCGASMSetType - Sets the type of restriction and interpolation used 1140 for local problems in the `PCGASM` additive Schwarz method. 1141 1142 Logically Collective 1143 1144 Input Parameters: 1145 + pc - the preconditioner context 1146 - type - variant of `PCGASM`, one of 1147 .vb 1148 `PC_GASM_BASIC` - full interpolation and restriction 1149 `PC_GASM_RESTRICT` - full restriction, local MPI rank interpolation 1150 `PC_GASM_INTERPOLATE` - full interpolation, local MPI rank restriction 1151 `PC_GASM_NONE` - local MPI rank restriction and interpolation 1152 .ve 1153 1154 Options Database Key: 1155 . -pc_gasm_type [basic,restrict,interpolate,none] - Sets `PCGASM` type 1156 1157 Level: intermediate 1158 1159 .seealso: `PCGASM`, `PCGASMSetSubdomains()`, `PCGASMGetSubKSP()`, 1160 `PCGASMCreateSubdomains2D()`, `PCASM`, `PCASMSetType()` 1161 @*/ 1162 PetscErrorCode PCGASMSetType(PC pc, PCGASMType type) 1163 { 1164 PetscFunctionBegin; 1165 PetscValidHeaderSpecific(pc, PC_CLASSID, 1); 1166 PetscValidLogicalCollectiveEnum(pc, type, 2); 1167 PetscTryMethod(pc, "PCGASMSetType_C", (PC, PCGASMType), (pc, type)); 1168 PetscFunctionReturn(PETSC_SUCCESS); 1169 } 1170 1171 /*@ 1172 PCGASMSetSortIndices - Determines whether subdomain indices are sorted. 1173 1174 Logically Collective 1175 1176 Input Parameters: 1177 + pc - the preconditioner context 1178 - doSort - sort the subdomain indices 1179 1180 Level: intermediate 1181 1182 .seealso: `PCGASM`, `PCGASMSetSubdomains()`, `PCGASMGetSubKSP()`, 1183 `PCGASMCreateSubdomains2D()` 1184 @*/ 1185 PetscErrorCode PCGASMSetSortIndices(PC pc, PetscBool doSort) 1186 { 1187 PetscFunctionBegin; 1188 PetscValidHeaderSpecific(pc, PC_CLASSID, 1); 1189 PetscValidLogicalCollectiveBool(pc, doSort, 2); 1190 PetscTryMethod(pc, "PCGASMSetSortIndices_C", (PC, PetscBool), (pc, doSort)); 1191 PetscFunctionReturn(PETSC_SUCCESS); 1192 } 1193 1194 /*@C 1195 PCGASMGetSubKSP - Gets the local `KSP` contexts for all subdomains on 1196 this MPI rank. 1197 1198 Collective iff first_local is requested 1199 1200 Input Parameter: 1201 . pc - the preconditioner context 1202 1203 Output Parameters: 1204 + n_local - the number of blocks on this MPI rank or `NULL` 1205 . first_local - the global number of the first block on this rank or `NULL`, 1206 all ranks must request or all must pass `NULL` 1207 - ksp - the array of `KSP` contexts 1208 1209 Level: advanced 1210 1211 Note: 1212 After `PCGASMGetSubKSP()` the array of `KSP`es is not to be freed 1213 1214 Currently for some matrix implementations only 1 block per MPI process 1215 is supported. 1216 1217 You must call `KSPSetUp()` before calling `PCGASMGetSubKSP()`. 1218 1219 .seealso: `PCGASM`, `PCGASMSetSubdomains()`, `PCGASMSetOverlap()`, 1220 `PCGASMCreateSubdomains2D()`, 1221 @*/ 1222 PetscErrorCode PCGASMGetSubKSP(PC pc, PetscInt *n_local, PetscInt *first_local, KSP *ksp[]) 1223 { 1224 PetscFunctionBegin; 1225 PetscValidHeaderSpecific(pc, PC_CLASSID, 1); 1226 PetscUseMethod(pc, "PCGASMGetSubKSP_C", (PC, PetscInt *, PetscInt *, KSP **), (pc, n_local, first_local, ksp)); 1227 PetscFunctionReturn(PETSC_SUCCESS); 1228 } 1229 1230 /*MC 1231 PCGASM - Use the (restricted) additive Schwarz method, each block is (approximately) solved with 1232 its own `KSP` object on a subset of MPI ranks 1233 1234 Options Database Keys: 1235 + -pc_gasm_total_subdomains <n> - Sets total number of local subdomains to be distributed among MPI ranks 1236 . -pc_gasm_view_subdomains - activates the printing of subdomain indices in `PCView()`, -ksp_view or -snes_view 1237 . -pc_gasm_print_subdomains - activates the printing of subdomain indices in `PCSetUp()` 1238 . -pc_gasm_overlap <ovl> - Sets overlap by which to (automatically) extend local subdomains 1239 - -pc_gasm_type [basic,restrict,interpolate,none] - Sets `PCGASMType` 1240 1241 Level: beginner 1242 1243 Notes: 1244 To set options on the solvers for each block append `-sub_` to all the `KSP`, and `PC` 1245 options database keys. For example, `-sub_pc_type ilu -sub_pc_factor_levels 1 -sub_ksp_type preonly` 1246 1247 To set the options on the solvers separate for each block call `PCGASMGetSubKSP()` 1248 and set the options directly on the resulting `KSP` object (you can access its `PC` 1249 with `KSPGetPC()`) 1250 1251 References: 1252 + * - M Dryja, OB Widlund, An additive variant of the Schwarz alternating method for the case of many subregions 1253 Courant Institute, New York University Technical report 1254 - * - Barry Smith, Petter Bjorstad, and William Gropp, Domain Decompositions: Parallel Multilevel Methods for Elliptic Partial Differential Equations, 1255 Cambridge University Press. 1256 1257 .seealso: `PCCreate()`, `PCSetType()`, `PCType`, `PC`, `PCASM`, `PCGASMType`, `PCGASMSetType()`, 1258 `PCBJACOBI`, `PCGASMGetSubKSP()`, `PCGASMSetSubdomains()`, 1259 `PCSetModifySubMatrices()`, `PCGASMSetOverlap()`, `PCGASMSetType()` 1260 M*/ 1261 1262 PETSC_EXTERN PetscErrorCode PCCreate_GASM(PC pc) 1263 { 1264 PC_GASM *osm; 1265 1266 PetscFunctionBegin; 1267 PetscCall(PetscNew(&osm)); 1268 1269 osm->N = PETSC_DETERMINE; 1270 osm->n = PETSC_DECIDE; 1271 osm->nmax = PETSC_DETERMINE; 1272 osm->overlap = 0; 1273 osm->ksp = NULL; 1274 osm->gorestriction = NULL; 1275 osm->girestriction = NULL; 1276 osm->pctoouter = NULL; 1277 osm->gx = NULL; 1278 osm->gy = NULL; 1279 osm->x = NULL; 1280 osm->y = NULL; 1281 osm->pcx = NULL; 1282 osm->pcy = NULL; 1283 osm->permutationIS = NULL; 1284 osm->permutationP = NULL; 1285 osm->pcmat = NULL; 1286 osm->ois = NULL; 1287 osm->iis = NULL; 1288 osm->pmat = NULL; 1289 osm->type = PC_GASM_RESTRICT; 1290 osm->same_subdomain_solvers = PETSC_TRUE; 1291 osm->sort_indices = PETSC_TRUE; 1292 osm->dm_subdomains = PETSC_FALSE; 1293 osm->hierarchicalpartitioning = PETSC_FALSE; 1294 1295 pc->data = (void *)osm; 1296 pc->ops->apply = PCApply_GASM; 1297 pc->ops->matapply = PCMatApply_GASM; 1298 pc->ops->applytranspose = PCApplyTranspose_GASM; 1299 pc->ops->setup = PCSetUp_GASM; 1300 pc->ops->reset = PCReset_GASM; 1301 pc->ops->destroy = PCDestroy_GASM; 1302 pc->ops->setfromoptions = PCSetFromOptions_GASM; 1303 pc->ops->setuponblocks = PCSetUpOnBlocks_GASM; 1304 pc->ops->view = PCView_GASM; 1305 pc->ops->applyrichardson = NULL; 1306 1307 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGASMSetSubdomains_C", PCGASMSetSubdomains_GASM)); 1308 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGASMSetOverlap_C", PCGASMSetOverlap_GASM)); 1309 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGASMSetType_C", PCGASMSetType_GASM)); 1310 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGASMSetSortIndices_C", PCGASMSetSortIndices_GASM)); 1311 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGASMGetSubKSP_C", PCGASMGetSubKSP_GASM)); 1312 PetscFunctionReturn(PETSC_SUCCESS); 1313 } 1314 1315 PetscErrorCode PCGASMCreateLocalSubdomains(Mat A, PetscInt nloc, IS *iis[]) 1316 { 1317 MatPartitioning mpart; 1318 const char *prefix; 1319 PetscInt i, j, rstart, rend, bs; 1320 PetscBool hasop, isbaij = PETSC_FALSE, foundpart = PETSC_FALSE; 1321 Mat Ad = NULL, adj; 1322 IS ispart, isnumb, *is; 1323 1324 PetscFunctionBegin; 1325 PetscCheck(nloc >= 1, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "number of local subdomains must > 0, got nloc = %" PetscInt_FMT, nloc); 1326 1327 /* Get prefix, row distribution, and block size */ 1328 PetscCall(MatGetOptionsPrefix(A, &prefix)); 1329 PetscCall(MatGetOwnershipRange(A, &rstart, &rend)); 1330 PetscCall(MatGetBlockSize(A, &bs)); 1331 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); 1332 1333 /* Get diagonal block from matrix if possible */ 1334 PetscCall(MatHasOperation(A, MATOP_GET_DIAGONAL_BLOCK, &hasop)); 1335 if (hasop) PetscCall(MatGetDiagonalBlock(A, &Ad)); 1336 if (Ad) { 1337 PetscCall(PetscObjectBaseTypeCompare((PetscObject)Ad, MATSEQBAIJ, &isbaij)); 1338 if (!isbaij) PetscCall(PetscObjectBaseTypeCompare((PetscObject)Ad, MATSEQSBAIJ, &isbaij)); 1339 } 1340 if (Ad && nloc > 1) { 1341 PetscBool match, done; 1342 /* Try to setup a good matrix partitioning if available */ 1343 PetscCall(MatPartitioningCreate(PETSC_COMM_SELF, &mpart)); 1344 PetscCall(PetscObjectSetOptionsPrefix((PetscObject)mpart, prefix)); 1345 PetscCall(MatPartitioningSetFromOptions(mpart)); 1346 PetscCall(PetscObjectTypeCompare((PetscObject)mpart, MATPARTITIONINGCURRENT, &match)); 1347 if (!match) PetscCall(PetscObjectTypeCompare((PetscObject)mpart, MATPARTITIONINGSQUARE, &match)); 1348 if (!match) { /* assume a "good" partitioner is available */ 1349 PetscInt na; 1350 const PetscInt *ia, *ja; 1351 PetscCall(MatGetRowIJ(Ad, 0, PETSC_TRUE, isbaij, &na, &ia, &ja, &done)); 1352 if (done) { 1353 /* Build adjacency matrix by hand. Unfortunately a call to 1354 MatConvert(Ad,MATMPIADJ,MAT_INITIAL_MATRIX,&adj) will 1355 remove the block-aij structure and we cannot expect 1356 MatPartitioning to split vertices as we need */ 1357 PetscInt i, j, len, nnz, cnt, *iia = NULL, *jja = NULL; 1358 const PetscInt *row; 1359 nnz = 0; 1360 for (i = 0; i < na; i++) { /* count number of nonzeros */ 1361 len = ia[i + 1] - ia[i]; 1362 row = ja + ia[i]; 1363 for (j = 0; j < len; j++) { 1364 if (row[j] == i) { /* don't count diagonal */ 1365 len--; 1366 break; 1367 } 1368 } 1369 nnz += len; 1370 } 1371 PetscCall(PetscMalloc1(na + 1, &iia)); 1372 PetscCall(PetscMalloc1(nnz, &jja)); 1373 nnz = 0; 1374 iia[0] = 0; 1375 for (i = 0; i < na; i++) { /* fill adjacency */ 1376 cnt = 0; 1377 len = ia[i + 1] - ia[i]; 1378 row = ja + ia[i]; 1379 for (j = 0; j < len; j++) { 1380 if (row[j] != i) jja[nnz + cnt++] = row[j]; /* if not diagonal */ 1381 } 1382 nnz += cnt; 1383 iia[i + 1] = nnz; 1384 } 1385 /* Partitioning of the adjacency matrix */ 1386 PetscCall(MatCreateMPIAdj(PETSC_COMM_SELF, na, na, iia, jja, NULL, &adj)); 1387 PetscCall(MatPartitioningSetAdjacency(mpart, adj)); 1388 PetscCall(MatPartitioningSetNParts(mpart, nloc)); 1389 PetscCall(MatPartitioningApply(mpart, &ispart)); 1390 PetscCall(ISPartitioningToNumbering(ispart, &isnumb)); 1391 PetscCall(MatDestroy(&adj)); 1392 foundpart = PETSC_TRUE; 1393 } 1394 PetscCall(MatRestoreRowIJ(Ad, 0, PETSC_TRUE, isbaij, &na, &ia, &ja, &done)); 1395 } 1396 PetscCall(MatPartitioningDestroy(&mpart)); 1397 } 1398 PetscCall(PetscMalloc1(nloc, &is)); 1399 if (!foundpart) { 1400 /* Partitioning by contiguous chunks of rows */ 1401 1402 PetscInt mbs = (rend - rstart) / bs; 1403 PetscInt start = rstart; 1404 for (i = 0; i < nloc; i++) { 1405 PetscInt count = (mbs / nloc + ((mbs % nloc) > i)) * bs; 1406 PetscCall(ISCreateStride(PETSC_COMM_SELF, count, start, 1, &is[i])); 1407 start += count; 1408 } 1409 1410 } else { 1411 /* Partitioning by adjacency of diagonal block */ 1412 1413 const PetscInt *numbering; 1414 PetscInt *count, nidx, *indices, *newidx, start = 0; 1415 /* Get node count in each partition */ 1416 PetscCall(PetscMalloc1(nloc, &count)); 1417 PetscCall(ISPartitioningCount(ispart, nloc, count)); 1418 if (isbaij && bs > 1) { /* adjust for the block-aij case */ 1419 for (i = 0; i < nloc; i++) count[i] *= bs; 1420 } 1421 /* Build indices from node numbering */ 1422 PetscCall(ISGetLocalSize(isnumb, &nidx)); 1423 PetscCall(PetscMalloc1(nidx, &indices)); 1424 for (i = 0; i < nidx; i++) indices[i] = i; /* needs to be initialized */ 1425 PetscCall(ISGetIndices(isnumb, &numbering)); 1426 PetscCall(PetscSortIntWithPermutation(nidx, numbering, indices)); 1427 PetscCall(ISRestoreIndices(isnumb, &numbering)); 1428 if (isbaij && bs > 1) { /* adjust for the block-aij case */ 1429 PetscCall(PetscMalloc1(nidx * bs, &newidx)); 1430 for (i = 0; i < nidx; i++) { 1431 for (j = 0; j < bs; j++) newidx[i * bs + j] = indices[i] * bs + j; 1432 } 1433 PetscCall(PetscFree(indices)); 1434 nidx *= bs; 1435 indices = newidx; 1436 } 1437 /* Shift to get global indices */ 1438 for (i = 0; i < nidx; i++) indices[i] += rstart; 1439 1440 /* Build the index sets for each block */ 1441 for (i = 0; i < nloc; i++) { 1442 PetscCall(ISCreateGeneral(PETSC_COMM_SELF, count[i], &indices[start], PETSC_COPY_VALUES, &is[i])); 1443 PetscCall(ISSort(is[i])); 1444 start += count[i]; 1445 } 1446 1447 PetscCall(PetscFree(count)); 1448 PetscCall(PetscFree(indices)); 1449 PetscCall(ISDestroy(&isnumb)); 1450 PetscCall(ISDestroy(&ispart)); 1451 } 1452 *iis = is; 1453 PetscFunctionReturn(PETSC_SUCCESS); 1454 } 1455 1456 PETSC_INTERN PetscErrorCode PCGASMCreateStraddlingSubdomains(Mat A, PetscInt N, PetscInt *n, IS *iis[]) 1457 { 1458 PetscFunctionBegin; 1459 PetscCall(MatSubdomainsCreateCoalesce(A, N, n, iis)); 1460 PetscFunctionReturn(PETSC_SUCCESS); 1461 } 1462 1463 /*@C 1464 PCGASMCreateSubdomains - Creates `n` index sets defining `n` nonoverlapping subdomains on this MPI process for the `PCGASM` additive 1465 Schwarz preconditioner for a any problem based on its matrix. 1466 1467 Collective 1468 1469 Input Parameters: 1470 + A - The global matrix operator 1471 - N - the number of global subdomains requested 1472 1473 Output Parameters: 1474 + n - the number of subdomains created on this MPI rank 1475 - iis - the array of index sets defining the local inner subdomains (on which the correction is applied) 1476 1477 Level: advanced 1478 1479 Notes: 1480 When `N` >= A's communicator size, each subdomain is local -- contained within a single MPI process. 1481 When `N` < size, the subdomains are 'straddling' (rank boundaries) and are no longer local. 1482 The resulting subdomains can be use in `PCGASMSetSubdomains`(pc,n,iss,`NULL`). The overlapping 1483 outer subdomains will be automatically generated from these according to the requested amount of 1484 overlap; this is currently supported only with local subdomains. 1485 1486 Use `PCGASMDestroySubdomains()` to free the array and the list of index sets. 1487 1488 .seealso: `PCGASM`, `PCGASMSetSubdomains()`, `PCGASMDestroySubdomains()` 1489 @*/ 1490 PetscErrorCode PCGASMCreateSubdomains(Mat A, PetscInt N, PetscInt *n, IS *iis[]) 1491 { 1492 PetscMPIInt size; 1493 1494 PetscFunctionBegin; 1495 PetscValidHeaderSpecific(A, MAT_CLASSID, 1); 1496 PetscAssertPointer(iis, 4); 1497 1498 PetscCheck(N >= 1, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Number of subdomains must be > 0, N = %" PetscInt_FMT, N); 1499 PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)A), &size)); 1500 if (N >= size) { 1501 *n = N / size + (N % size); 1502 PetscCall(PCGASMCreateLocalSubdomains(A, *n, iis)); 1503 } else { 1504 PetscCall(PCGASMCreateStraddlingSubdomains(A, N, n, iis)); 1505 } 1506 PetscFunctionReturn(PETSC_SUCCESS); 1507 } 1508 1509 /*@C 1510 PCGASMDestroySubdomains - Destroys the index sets created with 1511 `PCGASMCreateSubdomains()` or `PCGASMCreateSubdomains2D()`. Should be 1512 called after setting subdomains with `PCGASMSetSubdomains()`. 1513 1514 Collective 1515 1516 Input Parameters: 1517 + n - the number of index sets 1518 . iis - the array of inner subdomains 1519 - ois - the array of outer subdomains, can be `NULL` 1520 1521 Level: intermediate 1522 1523 Note: 1524 This is a convenience subroutine that walks each list, 1525 destroys each `IS` on the list, and then frees the list. At the end the 1526 list pointers are set to `NULL`. 1527 1528 .seealso: `PCGASM`, `PCGASMCreateSubdomains()`, `PCGASMSetSubdomains()` 1529 @*/ 1530 PetscErrorCode PCGASMDestroySubdomains(PetscInt n, IS **iis, IS **ois) 1531 { 1532 PetscInt i; 1533 1534 PetscFunctionBegin; 1535 if (n <= 0) PetscFunctionReturn(PETSC_SUCCESS); 1536 if (ois) { 1537 PetscAssertPointer(ois, 3); 1538 if (*ois) { 1539 PetscAssertPointer(*ois, 3); 1540 for (i = 0; i < n; i++) PetscCall(ISDestroy(&(*ois)[i])); 1541 PetscCall(PetscFree((*ois))); 1542 } 1543 } 1544 if (iis) { 1545 PetscAssertPointer(iis, 2); 1546 if (*iis) { 1547 PetscAssertPointer(*iis, 2); 1548 for (i = 0; i < n; i++) PetscCall(ISDestroy(&(*iis)[i])); 1549 PetscCall(PetscFree((*iis))); 1550 } 1551 } 1552 PetscFunctionReturn(PETSC_SUCCESS); 1553 } 1554 1555 #define PCGASMLocalSubdomainBounds2D(M, N, xleft, ylow, xright, yhigh, first, last, xleft_loc, ylow_loc, xright_loc, yhigh_loc, n) \ 1556 do { \ 1557 PetscInt first_row = first / M, last_row = last / M + 1; \ 1558 /* \ 1559 Compute ylow_loc and yhigh_loc so that (ylow_loc,xleft) and (yhigh_loc,xright) are the corners \ 1560 of the bounding box of the intersection of the subdomain with the local ownership range (local \ 1561 subdomain). \ 1562 Also compute xleft_loc and xright_loc as the lower and upper bounds on the first and last rows \ 1563 of the intersection. \ 1564 */ \ 1565 /* ylow_loc is the grid row containing the first element of the local sumbdomain */ \ 1566 *ylow_loc = PetscMax(first_row, ylow); \ 1567 /* xleft_loc is the offset of first element of the local subdomain within its grid row (might actually be outside the local subdomain) */ \ 1568 *xleft_loc = *ylow_loc == first_row ? PetscMax(first % M, xleft) : xleft; \ 1569 /* yhigh_loc is the grid row above the last local subdomain element */ \ 1570 *yhigh_loc = PetscMin(last_row, yhigh); \ 1571 /* xright is the offset of the end of the local subdomain within its grid row (might actually be outside the local subdomain) */ \ 1572 *xright_loc = *yhigh_loc == last_row ? PetscMin(xright, last % M) : xright; \ 1573 /* Now compute the size of the local subdomain n. */ \ 1574 *n = 0; \ 1575 if (*ylow_loc < *yhigh_loc) { \ 1576 PetscInt width = xright - xleft; \ 1577 *n += width * (*yhigh_loc - *ylow_loc - 1); \ 1578 *n += PetscMin(PetscMax(*xright_loc - xleft, 0), width); \ 1579 *n -= PetscMin(PetscMax(*xleft_loc - xleft, 0), width); \ 1580 } \ 1581 } while (0) 1582 1583 /*@C 1584 PCGASMCreateSubdomains2D - Creates the index sets for the `PCGASM` overlapping Schwarz 1585 preconditioner for a two-dimensional problem on a regular grid. 1586 1587 Collective 1588 1589 Input Parameters: 1590 + pc - the preconditioner context 1591 . M - the global number of grid points in the x direction 1592 . N - the global number of grid points in the y direction 1593 . Mdomains - the global number of subdomains in the x direction 1594 . Ndomains - the global number of subdomains in the y direction 1595 . dof - degrees of freedom per node 1596 - overlap - overlap in mesh lines 1597 1598 Output Parameters: 1599 + nsub - the number of local subdomains created 1600 . iis - array of index sets defining inner (nonoverlapping) subdomains 1601 - ois - array of index sets defining outer (overlapping, if overlap > 0) subdomains 1602 1603 Level: advanced 1604 1605 Note: 1606 Use `PCGASMDestroySubdomains()` to free the index sets and the arrays 1607 1608 Fortran Notes: 1609 The `IS` must be declared as an array of length long enough to hold `Nsub` entries 1610 1611 .seealso: `PCGASM`, `PCGASMSetSubdomains()`, `PCGASMGetSubKSP()`, `PCGASMSetOverlap()`, `PCASMCreateSubdomains2D()`, 1612 `PCGASMDestroySubdomains()` 1613 @*/ 1614 PetscErrorCode PCGASMCreateSubdomains2D(PC pc, PetscInt M, PetscInt N, PetscInt Mdomains, PetscInt Ndomains, PetscInt dof, PetscInt overlap, PetscInt *nsub, IS **iis, IS **ois) 1615 { 1616 PetscMPIInt size, rank; 1617 PetscInt i, j; 1618 PetscInt maxheight, maxwidth; 1619 PetscInt xstart, xleft, xright, xleft_loc, xright_loc; 1620 PetscInt ystart, ylow, yhigh, ylow_loc, yhigh_loc; 1621 PetscInt x[2][2], y[2][2], n[2]; 1622 PetscInt first, last; 1623 PetscInt nidx, *idx; 1624 PetscInt ii, jj, s, q, d; 1625 PetscInt k, kk; 1626 PetscMPIInt color; 1627 MPI_Comm comm, subcomm; 1628 IS **xis = NULL, **is = ois, **is_local = iis; 1629 1630 PetscFunctionBegin; 1631 PetscCall(PetscObjectGetComm((PetscObject)pc, &comm)); 1632 PetscCallMPI(MPI_Comm_size(comm, &size)); 1633 PetscCallMPI(MPI_Comm_rank(comm, &rank)); 1634 PetscCall(MatGetOwnershipRange(pc->pmat, &first, &last)); 1635 PetscCheck((first % dof) == 0 && (last % dof) == 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, 1636 "Matrix row partitioning unsuitable for domain decomposition: local row range (%" PetscInt_FMT ",%" PetscInt_FMT ") " 1637 "does not respect the number of degrees of freedom per grid point %" PetscInt_FMT, 1638 first, last, dof); 1639 1640 /* Determine the number of domains with nonzero intersections with the local ownership range. */ 1641 s = 0; 1642 ystart = 0; 1643 for (j = 0; j < Ndomains; ++j) { 1644 maxheight = N / Ndomains + ((N % Ndomains) > j); /* Maximal height of subdomain */ 1645 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); 1646 /* Vertical domain limits with an overlap. */ 1647 ylow = PetscMax(ystart - overlap, 0); 1648 yhigh = PetscMin(ystart + maxheight + overlap, N); 1649 xstart = 0; 1650 for (i = 0; i < Mdomains; ++i) { 1651 maxwidth = M / Mdomains + ((M % Mdomains) > i); /* Maximal width of subdomain */ 1652 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); 1653 /* Horizontal domain limits with an overlap. */ 1654 xleft = PetscMax(xstart - overlap, 0); 1655 xright = PetscMin(xstart + maxwidth + overlap, M); 1656 /* 1657 Determine whether this subdomain intersects this rank's ownership range of pc->pmat. 1658 */ 1659 PCGASMLocalSubdomainBounds2D(M, N, xleft, ylow, xright, yhigh, first, last, (&xleft_loc), (&ylow_loc), (&xright_loc), (&yhigh_loc), (&nidx)); 1660 if (nidx) ++s; 1661 xstart += maxwidth; 1662 } /* for (i = 0; i < Mdomains; ++i) */ 1663 ystart += maxheight; 1664 } /* for (j = 0; j < Ndomains; ++j) */ 1665 1666 /* Now we can allocate the necessary number of ISs. */ 1667 *nsub = s; 1668 PetscCall(PetscMalloc1(*nsub, is)); 1669 PetscCall(PetscMalloc1(*nsub, is_local)); 1670 s = 0; 1671 ystart = 0; 1672 for (j = 0; j < Ndomains; ++j) { 1673 maxheight = N / Ndomains + ((N % Ndomains) > j); /* Maximal height of subdomain */ 1674 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); 1675 /* Vertical domain limits with an overlap. */ 1676 y[0][0] = PetscMax(ystart - overlap, 0); 1677 y[0][1] = PetscMin(ystart + maxheight + overlap, N); 1678 /* Vertical domain limits without an overlap. */ 1679 y[1][0] = ystart; 1680 y[1][1] = PetscMin(ystart + maxheight, N); 1681 xstart = 0; 1682 for (i = 0; i < Mdomains; ++i) { 1683 maxwidth = M / Mdomains + ((M % Mdomains) > i); /* Maximal width of subdomain */ 1684 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); 1685 /* Horizontal domain limits with an overlap. */ 1686 x[0][0] = PetscMax(xstart - overlap, 0); 1687 x[0][1] = PetscMin(xstart + maxwidth + overlap, M); 1688 /* Horizontal domain limits without an overlap. */ 1689 x[1][0] = xstart; 1690 x[1][1] = PetscMin(xstart + maxwidth, M); 1691 /* 1692 Determine whether this domain intersects this rank's ownership range of pc->pmat. 1693 Do this twice: first for the domains with overlaps, and once without. 1694 During the first pass create the subcommunicators, and use them on the second pass as well. 1695 */ 1696 for (q = 0; q < 2; ++q) { 1697 PetscBool split = PETSC_FALSE; 1698 /* 1699 domain limits, (xleft, xright) and (ylow, yheigh) are adjusted 1700 according to whether the domain with an overlap or without is considered. 1701 */ 1702 xleft = x[q][0]; 1703 xright = x[q][1]; 1704 ylow = y[q][0]; 1705 yhigh = y[q][1]; 1706 PCGASMLocalSubdomainBounds2D(M, N, xleft, ylow, xright, yhigh, first, last, (&xleft_loc), (&ylow_loc), (&xright_loc), (&yhigh_loc), (&nidx)); 1707 nidx *= dof; 1708 n[q] = nidx; 1709 /* 1710 Based on the counted number of indices in the local domain *with an overlap*, 1711 construct a subcommunicator of all the MPI ranks supporting this domain. 1712 Observe that a domain with an overlap might have nontrivial local support, 1713 while the domain without an overlap might not. Hence, the decision to participate 1714 in the subcommunicator must be based on the domain with an overlap. 1715 */ 1716 if (q == 0) { 1717 if (nidx) color = 1; 1718 else color = MPI_UNDEFINED; 1719 PetscCallMPI(MPI_Comm_split(comm, color, rank, &subcomm)); 1720 split = PETSC_TRUE; 1721 } 1722 /* 1723 Proceed only if the number of local indices *with an overlap* is nonzero. 1724 */ 1725 if (n[0]) { 1726 if (q == 0) xis = is; 1727 if (q == 1) { 1728 /* 1729 The IS for the no-overlap subdomain shares a communicator with the overlapping domain. 1730 Moreover, if the overlap is zero, the two ISs are identical. 1731 */ 1732 if (overlap == 0) { 1733 (*is_local)[s] = (*is)[s]; 1734 PetscCall(PetscObjectReference((PetscObject)(*is)[s])); 1735 continue; 1736 } else { 1737 xis = is_local; 1738 subcomm = ((PetscObject)(*is)[s])->comm; 1739 } 1740 } /* if (q == 1) */ 1741 idx = NULL; 1742 PetscCall(PetscMalloc1(nidx, &idx)); 1743 if (nidx) { 1744 k = 0; 1745 for (jj = ylow_loc; jj < yhigh_loc; ++jj) { 1746 PetscInt x0 = (jj == ylow_loc) ? xleft_loc : xleft; 1747 PetscInt x1 = (jj == yhigh_loc - 1) ? xright_loc : xright; 1748 kk = dof * (M * jj + x0); 1749 for (ii = x0; ii < x1; ++ii) { 1750 for (d = 0; d < dof; ++d) idx[k++] = kk++; 1751 } 1752 } 1753 } 1754 PetscCall(ISCreateGeneral(subcomm, nidx, idx, PETSC_OWN_POINTER, (*xis) + s)); 1755 if (split) PetscCallMPI(MPI_Comm_free(&subcomm)); 1756 } /* if (n[0]) */ 1757 } /* for (q = 0; q < 2; ++q) */ 1758 if (n[0]) ++s; 1759 xstart += maxwidth; 1760 } /* for (i = 0; i < Mdomains; ++i) */ 1761 ystart += maxheight; 1762 } /* for (j = 0; j < Ndomains; ++j) */ 1763 PetscFunctionReturn(PETSC_SUCCESS); 1764 } 1765 1766 /*@C 1767 PCGASMGetSubdomains - Gets the subdomains supported on this MPI rank 1768 for the `PCGASM` additive Schwarz preconditioner. 1769 1770 Not Collective 1771 1772 Input Parameter: 1773 . pc - the preconditioner context 1774 1775 Output Parameters: 1776 + n - the number of subdomains for this MPI rank (default value = 1) 1777 . iis - the index sets that define the inner subdomains (without overlap) supported on this rank (can be `NULL`) 1778 - ois - the index sets that define the outer subdomains (with overlap) supported on this rank (can be `NULL`) 1779 1780 Level: advanced 1781 1782 Notes: 1783 The user is responsible for destroying the `IS`s and freeing the returned arrays, this can be done with 1784 `PCGASMDestroySubdomains()` 1785 1786 The `IS` numbering is in the parallel, global numbering of the vector. 1787 1788 .seealso: `PCGASM`, `PCGASMSetOverlap()`, `PCGASMGetSubKSP()`, `PCGASMCreateSubdomains2D()`, 1789 `PCGASMSetSubdomains()`, `PCGASMGetSubmatrices()`, `PCGASMDestroySubdomains()` 1790 @*/ 1791 PetscErrorCode PCGASMGetSubdomains(PC pc, PetscInt *n, IS *iis[], IS *ois[]) 1792 { 1793 PC_GASM *osm; 1794 PetscBool match; 1795 PetscInt i; 1796 1797 PetscFunctionBegin; 1798 PetscValidHeaderSpecific(pc, PC_CLASSID, 1); 1799 PetscCall(PetscObjectTypeCompare((PetscObject)pc, PCGASM, &match)); 1800 PetscCheck(match, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_WRONG, "Incorrect object type: expected %s, got %s instead", PCGASM, ((PetscObject)pc)->type_name); 1801 osm = (PC_GASM *)pc->data; 1802 if (n) *n = osm->n; 1803 if (iis) PetscCall(PetscMalloc1(osm->n, iis)); 1804 if (ois) PetscCall(PetscMalloc1(osm->n, ois)); 1805 if (iis || ois) { 1806 for (i = 0; i < osm->n; ++i) { 1807 if (iis) (*iis)[i] = osm->iis[i]; 1808 if (ois) (*ois)[i] = osm->ois[i]; 1809 } 1810 } 1811 PetscFunctionReturn(PETSC_SUCCESS); 1812 } 1813 1814 /*@C 1815 PCGASMGetSubmatrices - Gets the local submatrices (for this MPI rank 1816 only) for the `PCGASM` additive Schwarz preconditioner. 1817 1818 Not Collective 1819 1820 Input Parameter: 1821 . pc - the preconditioner context 1822 1823 Output Parameters: 1824 + n - the number of matrices for this MPI rank (default value = 1) 1825 - mat - the matrices 1826 1827 Level: advanced 1828 1829 Note: 1830 Matrices returned by this routine have the same communicators as the index sets (`IS`) 1831 used to define subdomains in `PCGASMSetSubdomains()` 1832 1833 .seealso: `PCGASM`, `PCGASMSetOverlap()`, `PCGASMGetSubKSP()`, 1834 `PCGASMCreateSubdomains2D()`, `PCGASMSetSubdomains()`, `PCGASMGetSubdomains()` 1835 @*/ 1836 PetscErrorCode PCGASMGetSubmatrices(PC pc, PetscInt *n, Mat *mat[]) 1837 { 1838 PC_GASM *osm; 1839 PetscBool match; 1840 1841 PetscFunctionBegin; 1842 PetscValidHeaderSpecific(pc, PC_CLASSID, 1); 1843 PetscAssertPointer(n, 2); 1844 if (mat) PetscAssertPointer(mat, 3); 1845 PetscCheck(pc->setupcalled, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_WRONGSTATE, "Must call after KSPSetUp() or PCSetUp()."); 1846 PetscCall(PetscObjectTypeCompare((PetscObject)pc, PCGASM, &match)); 1847 PetscCheck(match, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_WRONG, "Expected %s, got %s instead", PCGASM, ((PetscObject)pc)->type_name); 1848 osm = (PC_GASM *)pc->data; 1849 if (n) *n = osm->n; 1850 if (mat) *mat = osm->pmat; 1851 PetscFunctionReturn(PETSC_SUCCESS); 1852 } 1853 1854 /*@ 1855 PCGASMSetUseDMSubdomains - Indicates whether to use `DMCreateDomainDecomposition()` to define the subdomains, whenever possible for `PCGASM` 1856 1857 Logically Collective 1858 1859 Input Parameters: 1860 + pc - the preconditioner 1861 - flg - boolean indicating whether to use subdomains defined by the `DM` 1862 1863 Options Database Key: 1864 + -pc_gasm_dm_subdomains - configure subdomains 1865 . -pc_gasm_overlap - set overlap 1866 - -pc_gasm_total_subdomains - set number of subdomains 1867 1868 Level: intermediate 1869 1870 Note: 1871 `PCGASMSetSubdomains()`, `PCGASMSetTotalSubdomains()` or `PCGASMSetOverlap()` take precedence over `PCGASMSetUseDMSubdomains()`, 1872 so setting `PCGASMSetSubdomains()` with nontrivial subdomain ISs or any of `PCGASMSetTotalSubdomains()` and `PCGASMSetOverlap()` 1873 automatically turns the latter off. 1874 1875 .seealso: `PCGASM`, `PCGASMGetUseDMSubdomains()`, `PCGASMSetSubdomains()`, `PCGASMSetOverlap()` 1876 `PCGASMCreateSubdomains2D()` 1877 @*/ 1878 PetscErrorCode PCGASMSetUseDMSubdomains(PC pc, PetscBool flg) 1879 { 1880 PC_GASM *osm = (PC_GASM *)pc->data; 1881 PetscBool match; 1882 1883 PetscFunctionBegin; 1884 PetscValidHeaderSpecific(pc, PC_CLASSID, 1); 1885 PetscValidLogicalCollectiveBool(pc, flg, 2); 1886 PetscCheck(!pc->setupcalled, ((PetscObject)pc)->comm, PETSC_ERR_ARG_WRONGSTATE, "Not for a setup PC."); 1887 PetscCall(PetscObjectTypeCompare((PetscObject)pc, PCGASM, &match)); 1888 if (match) { 1889 if (!osm->user_subdomains && osm->N == PETSC_DETERMINE && osm->overlap < 0) osm->dm_subdomains = flg; 1890 } 1891 PetscFunctionReturn(PETSC_SUCCESS); 1892 } 1893 1894 /*@ 1895 PCGASMGetUseDMSubdomains - Returns flag indicating whether to use `DMCreateDomainDecomposition()` to define the subdomains, whenever possible with `PCGASM` 1896 1897 Not Collective 1898 1899 Input Parameter: 1900 . pc - the preconditioner 1901 1902 Output Parameter: 1903 . flg - boolean indicating whether to use subdomains defined by the `DM` 1904 1905 Level: intermediate 1906 1907 .seealso: `PCGASM`, `PCGASMSetUseDMSubdomains()`, `PCGASMSetOverlap()` 1908 `PCGASMCreateSubdomains2D()` 1909 @*/ 1910 PetscErrorCode PCGASMGetUseDMSubdomains(PC pc, PetscBool *flg) 1911 { 1912 PC_GASM *osm = (PC_GASM *)pc->data; 1913 PetscBool match; 1914 1915 PetscFunctionBegin; 1916 PetscValidHeaderSpecific(pc, PC_CLASSID, 1); 1917 PetscAssertPointer(flg, 2); 1918 PetscCall(PetscObjectTypeCompare((PetscObject)pc, PCGASM, &match)); 1919 if (match) { 1920 if (flg) *flg = osm->dm_subdomains; 1921 } 1922 PetscFunctionReturn(PETSC_SUCCESS); 1923 } 1924