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