1 /* 2 This file defines an "generalized" additive Schwarz preconditioner for any Mat implementation. 3 In this version each processor may intersect multiple subdomains and any subdomain may 4 intersect multiple processors. Intersections of subdomains with processors 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 processor (set in PCGASMSetSubdomains() or calculated in PCGASMSetTotalSubdomains()) 9 nmax - maximum number of local subdomains per processor (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 processors %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 processor 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 processor */ 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 882 communicator. 883 Logically collective on pc 884 885 Input Parameters: 886 + pc - the preconditioner 887 - N - total number of subdomains 888 889 Level: beginner 890 891 .seealso: `PCGASMSetSubdomains()`, `PCGASMSetOverlap()` 892 `PCGASMCreateSubdomains2D()` 893 @*/ 894 PetscErrorCode PCGASMSetTotalSubdomains(PC pc, PetscInt N) { 895 PC_GASM *osm = (PC_GASM *)pc->data; 896 PetscMPIInt size, rank; 897 898 PetscFunctionBegin; 899 PetscCheck(N >= 1, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Total number of subdomains must be 1 or more, got N = %" PetscInt_FMT, N); 900 PetscCheck(!pc->setupcalled, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_WRONGSTATE, "PCGASMSetTotalSubdomains() should be called before calling PCSetUp()."); 901 902 PetscCall(PCGASMDestroySubdomains(osm->n, &osm->iis, &osm->ois)); 903 osm->ois = osm->iis = NULL; 904 905 PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)pc), &size)); 906 PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)pc), &rank)); 907 osm->N = N; 908 osm->n = PETSC_DETERMINE; 909 osm->nmax = PETSC_DETERMINE; 910 osm->dm_subdomains = PETSC_FALSE; 911 PetscFunctionReturn(0); 912 } 913 914 static PetscErrorCode PCGASMSetSubdomains_GASM(PC pc, PetscInt n, IS iis[], IS ois[]) { 915 PC_GASM *osm = (PC_GASM *)pc->data; 916 PetscInt i; 917 918 PetscFunctionBegin; 919 PetscCheck(n >= 1, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Each process must have 1 or more subdomains, got n = %" PetscInt_FMT, n); 920 PetscCheck(!pc->setupcalled, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_WRONGSTATE, "PCGASMSetSubdomains() should be called before calling PCSetUp()."); 921 922 PetscCall(PCGASMDestroySubdomains(osm->n, &osm->iis, &osm->ois)); 923 osm->iis = osm->ois = NULL; 924 osm->n = n; 925 osm->N = PETSC_DETERMINE; 926 osm->nmax = PETSC_DETERMINE; 927 if (ois) { 928 PetscCall(PetscMalloc1(n, &osm->ois)); 929 for (i = 0; i < n; i++) { 930 PetscCall(PetscObjectReference((PetscObject)ois[i])); 931 osm->ois[i] = ois[i]; 932 } 933 /* 934 Since the user set the outer subdomains, even if nontrivial overlap was requested via PCGASMSetOverlap(), 935 it will be ignored. To avoid confusion later on (e.g., when viewing the PC), the overlap size is set to -1. 936 */ 937 osm->overlap = -1; 938 /* inner subdomains must be provided */ 939 PetscCheck(iis, PETSC_COMM_SELF, PETSC_ERR_ARG_NULL, "inner indices have to be provided "); 940 } /* end if */ 941 if (iis) { 942 PetscCall(PetscMalloc1(n, &osm->iis)); 943 for (i = 0; i < n; i++) { 944 PetscCall(PetscObjectReference((PetscObject)iis[i])); 945 osm->iis[i] = iis[i]; 946 } 947 if (!ois) { 948 osm->ois = NULL; 949 /* if user does not provide outer indices, we will create the corresponding outer indices using osm->overlap =1 in PCSetUp_GASM */ 950 } 951 } 952 if (PetscDefined(USE_DEBUG)) { 953 PetscInt j, rstart, rend, *covered, lsize; 954 const PetscInt *indices; 955 /* check if the inner indices cover and only cover the local portion of the preconditioning matrix */ 956 PetscCall(MatGetOwnershipRange(pc->pmat, &rstart, &rend)); 957 PetscCall(PetscCalloc1(rend - rstart, &covered)); 958 /* check if the current processor owns indices from others */ 959 for (i = 0; i < n; i++) { 960 PetscCall(ISGetIndices(osm->iis[i], &indices)); 961 PetscCall(ISGetLocalSize(osm->iis[i], &lsize)); 962 for (j = 0; j < lsize; j++) { 963 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 processors", indices[j]); 964 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]); 965 covered[indices[j] - rstart] = 1; 966 } 967 PetscCall(ISRestoreIndices(osm->iis[i], &indices)); 968 } 969 /* check if we miss any indices */ 970 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); } 971 PetscCall(PetscFree(covered)); 972 } 973 if (iis) osm->user_subdomains = PETSC_TRUE; 974 PetscFunctionReturn(0); 975 } 976 977 static PetscErrorCode PCGASMSetOverlap_GASM(PC pc, PetscInt ovl) { 978 PC_GASM *osm = (PC_GASM *)pc->data; 979 980 PetscFunctionBegin; 981 PetscCheck(ovl >= 0, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_OUTOFRANGE, "Negative overlap value requested"); 982 PetscCheck(!pc->setupcalled || ovl == osm->overlap, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_WRONGSTATE, "PCGASMSetOverlap() should be called before PCSetUp()."); 983 if (!pc->setupcalled) osm->overlap = ovl; 984 PetscFunctionReturn(0); 985 } 986 987 static PetscErrorCode PCGASMSetType_GASM(PC pc, PCGASMType type) { 988 PC_GASM *osm = (PC_GASM *)pc->data; 989 990 PetscFunctionBegin; 991 osm->type = type; 992 osm->type_set = PETSC_TRUE; 993 PetscFunctionReturn(0); 994 } 995 996 static PetscErrorCode PCGASMSetSortIndices_GASM(PC pc, PetscBool doSort) { 997 PC_GASM *osm = (PC_GASM *)pc->data; 998 999 PetscFunctionBegin; 1000 osm->sort_indices = doSort; 1001 PetscFunctionReturn(0); 1002 } 1003 1004 /* 1005 FIXME: This routine might need to be modified now that multiple ranks per subdomain are allowed. 1006 In particular, it would upset the global subdomain number calculation. 1007 */ 1008 static PetscErrorCode PCGASMGetSubKSP_GASM(PC pc, PetscInt *n, PetscInt *first, KSP **ksp) { 1009 PC_GASM *osm = (PC_GASM *)pc->data; 1010 1011 PetscFunctionBegin; 1012 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"); 1013 1014 if (n) *n = osm->n; 1015 if (first) { 1016 PetscCallMPI(MPI_Scan(&osm->n, first, 1, MPIU_INT, MPI_SUM, PetscObjectComm((PetscObject)pc))); 1017 *first -= osm->n; 1018 } 1019 if (ksp) { 1020 /* Assume that local solves are now different; not necessarily 1021 true, though! This flag is used only for PCView_GASM() */ 1022 *ksp = osm->ksp; 1023 osm->same_subdomain_solvers = PETSC_FALSE; 1024 } 1025 PetscFunctionReturn(0); 1026 } /* PCGASMGetSubKSP_GASM() */ 1027 1028 /*@C 1029 PCGASMSetSubdomains - Sets the subdomains for this processor 1030 for the additive Schwarz preconditioner. 1031 1032 Collective on pc 1033 1034 Input Parameters: 1035 + pc - the preconditioner object 1036 . n - the number of subdomains for this processor 1037 . iis - the index sets that define the inner subdomains (or NULL for PETSc to determine subdomains) 1038 - 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) 1039 1040 Notes: 1041 The IS indices use the parallel, global numbering of the vector entries. 1042 Inner subdomains are those where the correction is applied. 1043 Outer subdomains are those where the residual necessary to obtain the 1044 corrections is obtained (see PCGASMType for the use of inner/outer subdomains). 1045 Both inner and outer subdomains can extend over several processors. 1046 This processor's portion of a subdomain is known as a local subdomain. 1047 1048 Inner subdomains can not overlap with each other, do not have any entities from remote processors, 1049 and have to cover the entire local subdomain owned by the current processor. The index sets on each 1050 process should be ordered such that the ith local subdomain is connected to the ith remote subdomain 1051 on another MPI process. 1052 1053 By default the GASM preconditioner uses 1 (local) subdomain per processor. 1054 1055 Level: advanced 1056 1057 .seealso: `PCGASMSetOverlap()`, `PCGASMGetSubKSP()`, 1058 `PCGASMCreateSubdomains2D()`, `PCGASMGetSubdomains()` 1059 @*/ 1060 PetscErrorCode PCGASMSetSubdomains(PC pc, PetscInt n, IS iis[], IS ois[]) { 1061 PC_GASM *osm = (PC_GASM *)pc->data; 1062 1063 PetscFunctionBegin; 1064 PetscValidHeaderSpecific(pc, PC_CLASSID, 1); 1065 PetscTryMethod(pc, "PCGASMSetSubdomains_C", (PC, PetscInt, IS[], IS[]), (pc, n, iis, ois)); 1066 osm->dm_subdomains = PETSC_FALSE; 1067 PetscFunctionReturn(0); 1068 } 1069 1070 /*@ 1071 PCGASMSetOverlap - Sets the overlap between a pair of subdomains for the 1072 additive Schwarz preconditioner. Either all or no processors in the 1073 pc communicator must call this routine. 1074 1075 Logically Collective on pc 1076 1077 Input Parameters: 1078 + pc - the preconditioner context 1079 - ovl - the amount of overlap between subdomains (ovl >= 0, default value = 0) 1080 1081 Options Database Key: 1082 . -pc_gasm_overlap <overlap> - Sets overlap 1083 1084 Notes: 1085 By default the GASM preconditioner uses 1 subdomain per processor. To use 1086 multiple subdomain per perocessor or "straddling" subdomains that intersect 1087 multiple processors use PCGASMSetSubdomains() (or option -pc_gasm_total_subdomains <n>). 1088 1089 The overlap defaults to 0, so if one desires that no additional 1090 overlap be computed beyond what may have been set with a call to 1091 PCGASMSetSubdomains(), then ovl must be set to be 0. In particular, if one does 1092 not explicitly set the subdomains in application code, then all overlap would be computed 1093 internally by PETSc, and using an overlap of 0 would result in an GASM 1094 variant that is equivalent to the block Jacobi preconditioner. 1095 1096 Note that one can define initial index sets with any overlap via 1097 PCGASMSetSubdomains(); the routine PCGASMSetOverlap() merely allows 1098 PETSc to extend that overlap further, if desired. 1099 1100 Level: intermediate 1101 1102 .seealso: `PCGASMSetSubdomains()`, `PCGASMGetSubKSP()`, 1103 `PCGASMCreateSubdomains2D()`, `PCGASMGetSubdomains()` 1104 @*/ 1105 PetscErrorCode PCGASMSetOverlap(PC pc, PetscInt ovl) { 1106 PC_GASM *osm = (PC_GASM *)pc->data; 1107 1108 PetscFunctionBegin; 1109 PetscValidHeaderSpecific(pc, PC_CLASSID, 1); 1110 PetscValidLogicalCollectiveInt(pc, ovl, 2); 1111 PetscTryMethod(pc, "PCGASMSetOverlap_C", (PC, PetscInt), (pc, ovl)); 1112 osm->dm_subdomains = PETSC_FALSE; 1113 PetscFunctionReturn(0); 1114 } 1115 1116 /*@ 1117 PCGASMSetType - Sets the type of restriction and interpolation used 1118 for local problems in the additive Schwarz method. 1119 1120 Logically Collective on PC 1121 1122 Input Parameters: 1123 + pc - the preconditioner context 1124 - type - variant of GASM, one of 1125 .vb 1126 PC_GASM_BASIC - full interpolation and restriction 1127 PC_GASM_RESTRICT - full restriction, local processor interpolation 1128 PC_GASM_INTERPOLATE - full interpolation, local processor restriction 1129 PC_GASM_NONE - local processor restriction and interpolation 1130 .ve 1131 1132 Options Database Key: 1133 . -pc_gasm_type [basic,restrict,interpolate,none] - Sets GASM type 1134 1135 Level: intermediate 1136 1137 .seealso: `PCGASMSetSubdomains()`, `PCGASMGetSubKSP()`, 1138 `PCGASMCreateSubdomains2D()` 1139 @*/ 1140 PetscErrorCode PCGASMSetType(PC pc, PCGASMType type) { 1141 PetscFunctionBegin; 1142 PetscValidHeaderSpecific(pc, PC_CLASSID, 1); 1143 PetscValidLogicalCollectiveEnum(pc, type, 2); 1144 PetscTryMethod(pc, "PCGASMSetType_C", (PC, PCGASMType), (pc, type)); 1145 PetscFunctionReturn(0); 1146 } 1147 1148 /*@ 1149 PCGASMSetSortIndices - Determines whether subdomain indices are sorted. 1150 1151 Logically Collective on PC 1152 1153 Input Parameters: 1154 + pc - the preconditioner context 1155 - doSort - sort the subdomain indices 1156 1157 Level: intermediate 1158 1159 .seealso: `PCGASMSetSubdomains()`, `PCGASMGetSubKSP()`, 1160 `PCGASMCreateSubdomains2D()` 1161 @*/ 1162 PetscErrorCode PCGASMSetSortIndices(PC pc, PetscBool doSort) { 1163 PetscFunctionBegin; 1164 PetscValidHeaderSpecific(pc, PC_CLASSID, 1); 1165 PetscValidLogicalCollectiveBool(pc, doSort, 2); 1166 PetscTryMethod(pc, "PCGASMSetSortIndices_C", (PC, PetscBool), (pc, doSort)); 1167 PetscFunctionReturn(0); 1168 } 1169 1170 /*@C 1171 PCGASMGetSubKSP - Gets the local KSP contexts for all blocks on 1172 this processor. 1173 1174 Collective on PC iff first_local is requested 1175 1176 Input Parameter: 1177 . pc - the preconditioner context 1178 1179 Output Parameters: 1180 + n_local - the number of blocks on this processor or NULL 1181 . first_local - the global number of the first block on this processor or NULL, 1182 all processors must request or all must pass NULL 1183 - ksp - the array of KSP contexts 1184 1185 Note: 1186 After PCGASMGetSubKSP() the array of KSPes is not to be freed 1187 1188 Currently for some matrix implementations only 1 block per processor 1189 is supported. 1190 1191 You must call KSPSetUp() before calling PCGASMGetSubKSP(). 1192 1193 Level: advanced 1194 1195 .seealso: `PCGASMSetSubdomains()`, `PCGASMSetOverlap()`, 1196 `PCGASMCreateSubdomains2D()`, 1197 @*/ 1198 PetscErrorCode PCGASMGetSubKSP(PC pc, PetscInt *n_local, PetscInt *first_local, KSP *ksp[]) { 1199 PetscFunctionBegin; 1200 PetscValidHeaderSpecific(pc, PC_CLASSID, 1); 1201 PetscUseMethod(pc, "PCGASMGetSubKSP_C", (PC, PetscInt *, PetscInt *, KSP **), (pc, n_local, first_local, ksp)); 1202 PetscFunctionReturn(0); 1203 } 1204 1205 /* -------------------------------------------------------------------------------------*/ 1206 /*MC 1207 PCGASM - Use the (restricted) additive Schwarz method, each block is (approximately) solved with 1208 its own KSP object. 1209 1210 Options Database Keys: 1211 + -pc_gasm_total_subdomains <n> - Sets total number of local subdomains to be distributed among processors 1212 . -pc_gasm_view_subdomains - activates the printing of subdomain indices in PCView(), -ksp_view or -snes_view 1213 . -pc_gasm_print_subdomains - activates the printing of subdomain indices in PCSetUp() 1214 . -pc_gasm_overlap <ovl> - Sets overlap by which to (automatically) extend local subdomains 1215 - -pc_gasm_type [basic,restrict,interpolate,none] - Sets GASM type 1216 1217 IMPORTANT: If you run with, for example, 3 blocks on 1 processor or 3 blocks on 3 processors you 1218 will get a different convergence rate due to the default option of -pc_gasm_type restrict. Use 1219 -pc_gasm_type basic to use the standard GASM. 1220 1221 Notes: 1222 Blocks can be shared by multiple processes. 1223 1224 To set options on the solvers for each block append -sub_ to all the KSP, and PC 1225 options database keys. For example, -sub_pc_type ilu -sub_pc_factor_levels 1 -sub_ksp_type preonly 1226 1227 To set the options on the solvers separate for each block call PCGASMGetSubKSP() 1228 and set the options directly on the resulting KSP object (you can access its PC 1229 with KSPGetPC()) 1230 1231 Level: beginner 1232 1233 References: 1234 + * - M Dryja, OB Widlund, An additive variant of the Schwarz alternating method for the case of many subregions 1235 Courant Institute, New York University Technical report 1236 - * - Barry Smith, Petter Bjorstad, and William Gropp, Domain Decompositions: Parallel Multilevel Methods for Elliptic Partial Differential Equations, 1237 Cambridge University Press. 1238 1239 .seealso: `PCCreate()`, `PCSetType()`, `PCType`, `PC`, 1240 `PCBJACOBI`, `PCGASMGetSubKSP()`, `PCGASMSetSubdomains()`, 1241 `PCSetModifySubMatrices()`, `PCGASMSetOverlap()`, `PCGASMSetType()` 1242 1243 M*/ 1244 1245 PETSC_EXTERN PetscErrorCode PCCreate_GASM(PC pc) { 1246 PC_GASM *osm; 1247 1248 PetscFunctionBegin; 1249 PetscCall(PetscNewLog(pc, &osm)); 1250 1251 osm->N = PETSC_DETERMINE; 1252 osm->n = PETSC_DECIDE; 1253 osm->nmax = PETSC_DETERMINE; 1254 osm->overlap = 0; 1255 osm->ksp = NULL; 1256 osm->gorestriction = NULL; 1257 osm->girestriction = NULL; 1258 osm->pctoouter = NULL; 1259 osm->gx = NULL; 1260 osm->gy = NULL; 1261 osm->x = NULL; 1262 osm->y = NULL; 1263 osm->pcx = NULL; 1264 osm->pcy = NULL; 1265 osm->permutationIS = NULL; 1266 osm->permutationP = NULL; 1267 osm->pcmat = NULL; 1268 osm->ois = NULL; 1269 osm->iis = NULL; 1270 osm->pmat = NULL; 1271 osm->type = PC_GASM_RESTRICT; 1272 osm->same_subdomain_solvers = PETSC_TRUE; 1273 osm->sort_indices = PETSC_TRUE; 1274 osm->dm_subdomains = PETSC_FALSE; 1275 osm->hierarchicalpartitioning = PETSC_FALSE; 1276 1277 pc->data = (void *)osm; 1278 pc->ops->apply = PCApply_GASM; 1279 pc->ops->matapply = PCMatApply_GASM; 1280 pc->ops->applytranspose = PCApplyTranspose_GASM; 1281 pc->ops->setup = PCSetUp_GASM; 1282 pc->ops->reset = PCReset_GASM; 1283 pc->ops->destroy = PCDestroy_GASM; 1284 pc->ops->setfromoptions = PCSetFromOptions_GASM; 1285 pc->ops->setuponblocks = PCSetUpOnBlocks_GASM; 1286 pc->ops->view = PCView_GASM; 1287 pc->ops->applyrichardson = NULL; 1288 1289 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGASMSetSubdomains_C", PCGASMSetSubdomains_GASM)); 1290 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGASMSetOverlap_C", PCGASMSetOverlap_GASM)); 1291 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGASMSetType_C", PCGASMSetType_GASM)); 1292 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGASMSetSortIndices_C", PCGASMSetSortIndices_GASM)); 1293 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGASMGetSubKSP_C", PCGASMGetSubKSP_GASM)); 1294 PetscFunctionReturn(0); 1295 } 1296 1297 PetscErrorCode PCGASMCreateLocalSubdomains(Mat A, PetscInt nloc, IS *iis[]) { 1298 MatPartitioning mpart; 1299 const char *prefix; 1300 PetscInt i, j, rstart, rend, bs; 1301 PetscBool hasop, isbaij = PETSC_FALSE, foundpart = PETSC_FALSE; 1302 Mat Ad = NULL, adj; 1303 IS ispart, isnumb, *is; 1304 1305 PetscFunctionBegin; 1306 PetscCheck(nloc >= 1, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "number of local subdomains must > 0, got nloc = %" PetscInt_FMT, nloc); 1307 1308 /* Get prefix, row distribution, and block size */ 1309 PetscCall(MatGetOptionsPrefix(A, &prefix)); 1310 PetscCall(MatGetOwnershipRange(A, &rstart, &rend)); 1311 PetscCall(MatGetBlockSize(A, &bs)); 1312 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); 1313 1314 /* Get diagonal block from matrix if possible */ 1315 PetscCall(MatHasOperation(A, MATOP_GET_DIAGONAL_BLOCK, &hasop)); 1316 if (hasop) { PetscCall(MatGetDiagonalBlock(A, &Ad)); } 1317 if (Ad) { 1318 PetscCall(PetscObjectBaseTypeCompare((PetscObject)Ad, MATSEQBAIJ, &isbaij)); 1319 if (!isbaij) PetscCall(PetscObjectBaseTypeCompare((PetscObject)Ad, MATSEQSBAIJ, &isbaij)); 1320 } 1321 if (Ad && nloc > 1) { 1322 PetscBool match, done; 1323 /* Try to setup a good matrix partitioning if available */ 1324 PetscCall(MatPartitioningCreate(PETSC_COMM_SELF, &mpart)); 1325 PetscCall(PetscObjectSetOptionsPrefix((PetscObject)mpart, prefix)); 1326 PetscCall(MatPartitioningSetFromOptions(mpart)); 1327 PetscCall(PetscObjectTypeCompare((PetscObject)mpart, MATPARTITIONINGCURRENT, &match)); 1328 if (!match) { PetscCall(PetscObjectTypeCompare((PetscObject)mpart, MATPARTITIONINGSQUARE, &match)); } 1329 if (!match) { /* assume a "good" partitioner is available */ 1330 PetscInt na; 1331 const PetscInt *ia, *ja; 1332 PetscCall(MatGetRowIJ(Ad, 0, PETSC_TRUE, isbaij, &na, &ia, &ja, &done)); 1333 if (done) { 1334 /* Build adjacency matrix by hand. Unfortunately a call to 1335 MatConvert(Ad,MATMPIADJ,MAT_INITIAL_MATRIX,&adj) will 1336 remove the block-aij structure and we cannot expect 1337 MatPartitioning to split vertices as we need */ 1338 PetscInt i, j, len, nnz, cnt, *iia = NULL, *jja = NULL; 1339 const PetscInt *row; 1340 nnz = 0; 1341 for (i = 0; i < na; i++) { /* count number of nonzeros */ 1342 len = ia[i + 1] - ia[i]; 1343 row = ja + ia[i]; 1344 for (j = 0; j < len; j++) { 1345 if (row[j] == i) { /* don't count diagonal */ 1346 len--; 1347 break; 1348 } 1349 } 1350 nnz += len; 1351 } 1352 PetscCall(PetscMalloc1(na + 1, &iia)); 1353 PetscCall(PetscMalloc1(nnz, &jja)); 1354 nnz = 0; 1355 iia[0] = 0; 1356 for (i = 0; i < na; i++) { /* fill adjacency */ 1357 cnt = 0; 1358 len = ia[i + 1] - ia[i]; 1359 row = ja + ia[i]; 1360 for (j = 0; j < len; j++) { 1361 if (row[j] != i) jja[nnz + cnt++] = row[j]; /* if not diagonal */ 1362 } 1363 nnz += cnt; 1364 iia[i + 1] = nnz; 1365 } 1366 /* Partitioning of the adjacency matrix */ 1367 PetscCall(MatCreateMPIAdj(PETSC_COMM_SELF, na, na, iia, jja, NULL, &adj)); 1368 PetscCall(MatPartitioningSetAdjacency(mpart, adj)); 1369 PetscCall(MatPartitioningSetNParts(mpart, nloc)); 1370 PetscCall(MatPartitioningApply(mpart, &ispart)); 1371 PetscCall(ISPartitioningToNumbering(ispart, &isnumb)); 1372 PetscCall(MatDestroy(&adj)); 1373 foundpart = PETSC_TRUE; 1374 } 1375 PetscCall(MatRestoreRowIJ(Ad, 0, PETSC_TRUE, isbaij, &na, &ia, &ja, &done)); 1376 } 1377 PetscCall(MatPartitioningDestroy(&mpart)); 1378 } 1379 PetscCall(PetscMalloc1(nloc, &is)); 1380 if (!foundpart) { 1381 /* Partitioning by contiguous chunks of rows */ 1382 1383 PetscInt mbs = (rend - rstart) / bs; 1384 PetscInt start = rstart; 1385 for (i = 0; i < nloc; i++) { 1386 PetscInt count = (mbs / nloc + ((mbs % nloc) > i)) * bs; 1387 PetscCall(ISCreateStride(PETSC_COMM_SELF, count, start, 1, &is[i])); 1388 start += count; 1389 } 1390 1391 } else { 1392 /* Partitioning by adjacency of diagonal block */ 1393 1394 const PetscInt *numbering; 1395 PetscInt *count, nidx, *indices, *newidx, start = 0; 1396 /* Get node count in each partition */ 1397 PetscCall(PetscMalloc1(nloc, &count)); 1398 PetscCall(ISPartitioningCount(ispart, nloc, count)); 1399 if (isbaij && bs > 1) { /* adjust for the block-aij case */ 1400 for (i = 0; i < nloc; i++) count[i] *= bs; 1401 } 1402 /* Build indices from node numbering */ 1403 PetscCall(ISGetLocalSize(isnumb, &nidx)); 1404 PetscCall(PetscMalloc1(nidx, &indices)); 1405 for (i = 0; i < nidx; i++) indices[i] = i; /* needs to be initialized */ 1406 PetscCall(ISGetIndices(isnumb, &numbering)); 1407 PetscCall(PetscSortIntWithPermutation(nidx, numbering, indices)); 1408 PetscCall(ISRestoreIndices(isnumb, &numbering)); 1409 if (isbaij && bs > 1) { /* adjust for the block-aij case */ 1410 PetscCall(PetscMalloc1(nidx * bs, &newidx)); 1411 for (i = 0; i < nidx; i++) { 1412 for (j = 0; j < bs; j++) newidx[i * bs + j] = indices[i] * bs + j; 1413 } 1414 PetscCall(PetscFree(indices)); 1415 nidx *= bs; 1416 indices = newidx; 1417 } 1418 /* Shift to get global indices */ 1419 for (i = 0; i < nidx; i++) indices[i] += rstart; 1420 1421 /* Build the index sets for each block */ 1422 for (i = 0; i < nloc; i++) { 1423 PetscCall(ISCreateGeneral(PETSC_COMM_SELF, count[i], &indices[start], PETSC_COPY_VALUES, &is[i])); 1424 PetscCall(ISSort(is[i])); 1425 start += count[i]; 1426 } 1427 1428 PetscCall(PetscFree(count)); 1429 PetscCall(PetscFree(indices)); 1430 PetscCall(ISDestroy(&isnumb)); 1431 PetscCall(ISDestroy(&ispart)); 1432 } 1433 *iis = is; 1434 PetscFunctionReturn(0); 1435 } 1436 1437 PETSC_INTERN PetscErrorCode PCGASMCreateStraddlingSubdomains(Mat A, PetscInt N, PetscInt *n, IS *iis[]) { 1438 PetscFunctionBegin; 1439 PetscCall(MatSubdomainsCreateCoalesce(A, N, n, iis)); 1440 PetscFunctionReturn(0); 1441 } 1442 1443 /*@C 1444 PCGASMCreateSubdomains - Creates n index sets defining n nonoverlapping subdomains for the additive 1445 Schwarz preconditioner for a any problem based on its matrix. 1446 1447 Collective 1448 1449 Input Parameters: 1450 + A - The global matrix operator 1451 - N - the number of global subdomains requested 1452 1453 Output Parameters: 1454 + n - the number of subdomains created on this processor 1455 - iis - the array of index sets defining the local inner subdomains (on which the correction is applied) 1456 1457 Level: advanced 1458 1459 Note: When N >= A's communicator size, each subdomain is local -- contained within a single processor. 1460 When N < size, the subdomains are 'straddling' (processor boundaries) and are no longer local. 1461 The resulting subdomains can be use in PCGASMSetSubdomains(pc,n,iss,NULL). The overlapping 1462 outer subdomains will be automatically generated from these according to the requested amount of 1463 overlap; this is currently supported only with local subdomains. 1464 1465 .seealso: `PCGASMSetSubdomains()`, `PCGASMDestroySubdomains()` 1466 @*/ 1467 PetscErrorCode PCGASMCreateSubdomains(Mat A, PetscInt N, PetscInt *n, IS *iis[]) { 1468 PetscMPIInt size; 1469 1470 PetscFunctionBegin; 1471 PetscValidHeaderSpecific(A, MAT_CLASSID, 1); 1472 PetscValidPointer(iis, 4); 1473 1474 PetscCheck(N >= 1, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Number of subdomains must be > 0, N = %" PetscInt_FMT, N); 1475 PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)A), &size)); 1476 if (N >= size) { 1477 *n = N / size + (N % size); 1478 PetscCall(PCGASMCreateLocalSubdomains(A, *n, iis)); 1479 } else { 1480 PetscCall(PCGASMCreateStraddlingSubdomains(A, N, n, iis)); 1481 } 1482 PetscFunctionReturn(0); 1483 } 1484 1485 /*@C 1486 PCGASMDestroySubdomains - Destroys the index sets created with 1487 PCGASMCreateSubdomains() or PCGASMCreateSubdomains2D. Should be 1488 called after setting subdomains with PCGASMSetSubdomains(). 1489 1490 Collective 1491 1492 Input Parameters: 1493 + n - the number of index sets 1494 . iis - the array of inner subdomains, 1495 - ois - the array of outer subdomains, can be NULL 1496 1497 Level: intermediate 1498 1499 Notes: 1500 this is merely a convenience subroutine that walks each list, 1501 destroys each IS on the list, and then frees the list. At the end the 1502 list pointers are set to NULL. 1503 1504 .seealso: `PCGASMCreateSubdomains()`, `PCGASMSetSubdomains()` 1505 @*/ 1506 PetscErrorCode PCGASMDestroySubdomains(PetscInt n, IS **iis, IS **ois) { 1507 PetscInt i; 1508 1509 PetscFunctionBegin; 1510 if (n <= 0) PetscFunctionReturn(0); 1511 if (ois) { 1512 PetscValidPointer(ois, 3); 1513 if (*ois) { 1514 PetscValidPointer(*ois, 3); 1515 for (i = 0; i < n; i++) { PetscCall(ISDestroy(&(*ois)[i])); } 1516 PetscCall(PetscFree((*ois))); 1517 } 1518 } 1519 if (iis) { 1520 PetscValidPointer(iis, 2); 1521 if (*iis) { 1522 PetscValidPointer(*iis, 2); 1523 for (i = 0; i < n; i++) { PetscCall(ISDestroy(&(*iis)[i])); } 1524 PetscCall(PetscFree((*iis))); 1525 } 1526 } 1527 PetscFunctionReturn(0); 1528 } 1529 1530 #define PCGASMLocalSubdomainBounds2D(M, N, xleft, ylow, xright, yhigh, first, last, xleft_loc, ylow_loc, xright_loc, yhigh_loc, n) \ 1531 { \ 1532 PetscInt first_row = first / M, last_row = last / M + 1; \ 1533 /* \ 1534 Compute ylow_loc and yhigh_loc so that (ylow_loc,xleft) and (yhigh_loc,xright) are the corners \ 1535 of the bounding box of the intersection of the subdomain with the local ownership range (local \ 1536 subdomain). \ 1537 Also compute xleft_loc and xright_loc as the lower and upper bounds on the first and last rows \ 1538 of the intersection. \ 1539 */ \ 1540 /* ylow_loc is the grid row containing the first element of the local sumbdomain */ \ 1541 *ylow_loc = PetscMax(first_row, ylow); \ 1542 /* xleft_loc is the offset of first element of the local subdomain within its grid row (might actually be outside the local subdomain) */ \ 1543 *xleft_loc = *ylow_loc == first_row ? PetscMax(first % M, xleft) : xleft; \ 1544 /* yhigh_loc is the grid row above the last local subdomain element */ \ 1545 *yhigh_loc = PetscMin(last_row, yhigh); \ 1546 /* xright is the offset of the end of the local subdomain within its grid row (might actually be outside the local subdomain) */ \ 1547 *xright_loc = *yhigh_loc == last_row ? PetscMin(xright, last % M) : xright; \ 1548 /* Now compute the size of the local subdomain n. */ \ 1549 *n = 0; \ 1550 if (*ylow_loc < *yhigh_loc) { \ 1551 PetscInt width = xright - xleft; \ 1552 *n += width * (*yhigh_loc - *ylow_loc - 1); \ 1553 *n += PetscMin(PetscMax(*xright_loc - xleft, 0), width); \ 1554 *n -= PetscMin(PetscMax(*xleft_loc - xleft, 0), width); \ 1555 } \ 1556 } 1557 1558 /*@ 1559 PCGASMCreateSubdomains2D - Creates the index sets for the overlapping Schwarz 1560 preconditioner for a two-dimensional problem on a regular grid. 1561 1562 Collective 1563 1564 Input Parameters: 1565 + pc - the preconditioner context 1566 . M - the global number of grid points in the x direction 1567 . N - the global number of grid points in the y direction 1568 . Mdomains - the global number of subdomains in the x direction 1569 . Ndomains - the global number of subdomains in the y direction 1570 . dof - degrees of freedom per node 1571 - overlap - overlap in mesh lines 1572 1573 Output Parameters: 1574 + Nsub - the number of local subdomains created 1575 . iis - array of index sets defining inner (nonoverlapping) subdomains 1576 - ois - array of index sets defining outer (overlapping, if overlap > 0) subdomains 1577 1578 Level: advanced 1579 1580 .seealso: `PCGASMSetSubdomains()`, `PCGASMGetSubKSP()`, `PCGASMSetOverlap()` 1581 @*/ 1582 PetscErrorCode PCGASMCreateSubdomains2D(PC pc, PetscInt M, PetscInt N, PetscInt Mdomains, PetscInt Ndomains, PetscInt dof, PetscInt overlap, PetscInt *nsub, IS **iis, IS **ois) { 1583 PetscMPIInt size, rank; 1584 PetscInt i, j; 1585 PetscInt maxheight, maxwidth; 1586 PetscInt xstart, xleft, xright, xleft_loc, xright_loc; 1587 PetscInt ystart, ylow, yhigh, ylow_loc, yhigh_loc; 1588 PetscInt x[2][2], y[2][2], n[2]; 1589 PetscInt first, last; 1590 PetscInt nidx, *idx; 1591 PetscInt ii, jj, s, q, d; 1592 PetscInt k, kk; 1593 PetscMPIInt color; 1594 MPI_Comm comm, subcomm; 1595 IS **xis = NULL, **is = ois, **is_local = iis; 1596 1597 PetscFunctionBegin; 1598 PetscCall(PetscObjectGetComm((PetscObject)pc, &comm)); 1599 PetscCallMPI(MPI_Comm_size(comm, &size)); 1600 PetscCallMPI(MPI_Comm_rank(comm, &rank)); 1601 PetscCall(MatGetOwnershipRange(pc->pmat, &first, &last)); 1602 PetscCheck((first % dof) == 0 && (last % dof) == 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, 1603 "Matrix row partitioning unsuitable for domain decomposition: local row range (%" PetscInt_FMT ",%" PetscInt_FMT ") " 1604 "does not respect the number of degrees of freedom per grid point %" PetscInt_FMT, 1605 first, last, dof); 1606 1607 /* Determine the number of domains with nonzero intersections with the local ownership range. */ 1608 s = 0; 1609 ystart = 0; 1610 for (j = 0; j < Ndomains; ++j) { 1611 maxheight = N / Ndomains + ((N % Ndomains) > j); /* Maximal height of subdomain */ 1612 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); 1613 /* Vertical domain limits with an overlap. */ 1614 ylow = PetscMax(ystart - overlap, 0); 1615 yhigh = PetscMin(ystart + maxheight + overlap, N); 1616 xstart = 0; 1617 for (i = 0; i < Mdomains; ++i) { 1618 maxwidth = M / Mdomains + ((M % Mdomains) > i); /* Maximal width of subdomain */ 1619 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); 1620 /* Horizontal domain limits with an overlap. */ 1621 xleft = PetscMax(xstart - overlap, 0); 1622 xright = PetscMin(xstart + maxwidth + overlap, M); 1623 /* 1624 Determine whether this subdomain intersects this processor's ownership range of pc->pmat. 1625 */ 1626 PCGASMLocalSubdomainBounds2D(M, N, xleft, ylow, xright, yhigh, first, last, (&xleft_loc), (&ylow_loc), (&xright_loc), (&yhigh_loc), (&nidx)); 1627 if (nidx) ++s; 1628 xstart += maxwidth; 1629 } /* for (i = 0; i < Mdomains; ++i) */ 1630 ystart += maxheight; 1631 } /* for (j = 0; j < Ndomains; ++j) */ 1632 1633 /* Now we can allocate the necessary number of ISs. */ 1634 *nsub = s; 1635 PetscCall(PetscMalloc1(*nsub, is)); 1636 PetscCall(PetscMalloc1(*nsub, is_local)); 1637 s = 0; 1638 ystart = 0; 1639 for (j = 0; j < Ndomains; ++j) { 1640 maxheight = N / Ndomains + ((N % Ndomains) > j); /* Maximal height of subdomain */ 1641 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); 1642 /* Vertical domain limits with an overlap. */ 1643 y[0][0] = PetscMax(ystart - overlap, 0); 1644 y[0][1] = PetscMin(ystart + maxheight + overlap, N); 1645 /* Vertical domain limits without an overlap. */ 1646 y[1][0] = ystart; 1647 y[1][1] = PetscMin(ystart + maxheight, N); 1648 xstart = 0; 1649 for (i = 0; i < Mdomains; ++i) { 1650 maxwidth = M / Mdomains + ((M % Mdomains) > i); /* Maximal width of subdomain */ 1651 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); 1652 /* Horizontal domain limits with an overlap. */ 1653 x[0][0] = PetscMax(xstart - overlap, 0); 1654 x[0][1] = PetscMin(xstart + maxwidth + overlap, M); 1655 /* Horizontal domain limits without an overlap. */ 1656 x[1][0] = xstart; 1657 x[1][1] = PetscMin(xstart + maxwidth, M); 1658 /* 1659 Determine whether this domain intersects this processor's ownership range of pc->pmat. 1660 Do this twice: first for the domains with overlaps, and once without. 1661 During the first pass create the subcommunicators, and use them on the second pass as well. 1662 */ 1663 for (q = 0; q < 2; ++q) { 1664 PetscBool split = PETSC_FALSE; 1665 /* 1666 domain limits, (xleft, xright) and (ylow, yheigh) are adjusted 1667 according to whether the domain with an overlap or without is considered. 1668 */ 1669 xleft = x[q][0]; 1670 xright = x[q][1]; 1671 ylow = y[q][0]; 1672 yhigh = y[q][1]; 1673 PCGASMLocalSubdomainBounds2D(M, N, xleft, ylow, xright, yhigh, first, last, (&xleft_loc), (&ylow_loc), (&xright_loc), (&yhigh_loc), (&nidx)); 1674 nidx *= dof; 1675 n[q] = nidx; 1676 /* 1677 Based on the counted number of indices in the local domain *with an overlap*, 1678 construct a subcommunicator of all the processors supporting this domain. 1679 Observe that a domain with an overlap might have nontrivial local support, 1680 while the domain without an overlap might not. Hence, the decision to participate 1681 in the subcommunicator must be based on the domain with an overlap. 1682 */ 1683 if (q == 0) { 1684 if (nidx) color = 1; 1685 else color = MPI_UNDEFINED; 1686 PetscCallMPI(MPI_Comm_split(comm, color, rank, &subcomm)); 1687 split = PETSC_TRUE; 1688 } 1689 /* 1690 Proceed only if the number of local indices *with an overlap* is nonzero. 1691 */ 1692 if (n[0]) { 1693 if (q == 0) xis = is; 1694 if (q == 1) { 1695 /* 1696 The IS for the no-overlap subdomain shares a communicator with the overlapping domain. 1697 Moreover, if the overlap is zero, the two ISs are identical. 1698 */ 1699 if (overlap == 0) { 1700 (*is_local)[s] = (*is)[s]; 1701 PetscCall(PetscObjectReference((PetscObject)(*is)[s])); 1702 continue; 1703 } else { 1704 xis = is_local; 1705 subcomm = ((PetscObject)(*is)[s])->comm; 1706 } 1707 } /* if (q == 1) */ 1708 idx = NULL; 1709 PetscCall(PetscMalloc1(nidx, &idx)); 1710 if (nidx) { 1711 k = 0; 1712 for (jj = ylow_loc; jj < yhigh_loc; ++jj) { 1713 PetscInt x0 = (jj == ylow_loc) ? xleft_loc : xleft; 1714 PetscInt x1 = (jj == yhigh_loc - 1) ? xright_loc : xright; 1715 kk = dof * (M * jj + x0); 1716 for (ii = x0; ii < x1; ++ii) { 1717 for (d = 0; d < dof; ++d) { idx[k++] = kk++; } 1718 } 1719 } 1720 } 1721 PetscCall(ISCreateGeneral(subcomm, nidx, idx, PETSC_OWN_POINTER, (*xis) + s)); 1722 if (split) { PetscCallMPI(MPI_Comm_free(&subcomm)); } 1723 } /* if (n[0]) */ 1724 } /* for (q = 0; q < 2; ++q) */ 1725 if (n[0]) ++s; 1726 xstart += maxwidth; 1727 } /* for (i = 0; i < Mdomains; ++i) */ 1728 ystart += maxheight; 1729 } /* for (j = 0; j < Ndomains; ++j) */ 1730 PetscFunctionReturn(0); 1731 } 1732 1733 /*@C 1734 PCGASMGetSubdomains - Gets the subdomains supported on this processor 1735 for the additive Schwarz preconditioner. 1736 1737 Not Collective 1738 1739 Input Parameter: 1740 . pc - the preconditioner context 1741 1742 Output Parameters: 1743 + n - the number of subdomains for this processor (default value = 1) 1744 . iis - the index sets that define the inner subdomains (without overlap) supported on this processor (can be NULL) 1745 - ois - the index sets that define the outer subdomains (with overlap) supported on this processor (can be NULL) 1746 1747 Notes: 1748 The user is responsible for destroying the ISs and freeing the returned arrays. 1749 The IS numbering is in the parallel, global numbering of the vector. 1750 1751 Level: advanced 1752 1753 .seealso: `PCGASMSetOverlap()`, `PCGASMGetSubKSP()`, `PCGASMCreateSubdomains2D()`, 1754 `PCGASMSetSubdomains()`, `PCGASMGetSubmatrices()` 1755 @*/ 1756 PetscErrorCode PCGASMGetSubdomains(PC pc, PetscInt *n, IS *iis[], IS *ois[]) { 1757 PC_GASM *osm; 1758 PetscBool match; 1759 PetscInt i; 1760 1761 PetscFunctionBegin; 1762 PetscValidHeaderSpecific(pc, PC_CLASSID, 1); 1763 PetscCall(PetscObjectTypeCompare((PetscObject)pc, PCGASM, &match)); 1764 PetscCheck(match, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_WRONG, "Incorrect object type: expected %s, got %s instead", PCGASM, ((PetscObject)pc)->type_name); 1765 osm = (PC_GASM *)pc->data; 1766 if (n) *n = osm->n; 1767 if (iis) PetscCall(PetscMalloc1(osm->n, iis)); 1768 if (ois) PetscCall(PetscMalloc1(osm->n, ois)); 1769 if (iis || ois) { 1770 for (i = 0; i < osm->n; ++i) { 1771 if (iis) (*iis)[i] = osm->iis[i]; 1772 if (ois) (*ois)[i] = osm->ois[i]; 1773 } 1774 } 1775 PetscFunctionReturn(0); 1776 } 1777 1778 /*@C 1779 PCGASMGetSubmatrices - Gets the local submatrices (for this processor 1780 only) for the additive Schwarz preconditioner. 1781 1782 Not Collective 1783 1784 Input Parameter: 1785 . pc - the preconditioner context 1786 1787 Output Parameters: 1788 + n - the number of matrices for this processor (default value = 1) 1789 - mat - the matrices 1790 1791 Notes: 1792 matrices returned by this routine have the same communicators as the index sets (IS) 1793 used to define subdomains in PCGASMSetSubdomains() 1794 Level: advanced 1795 1796 .seealso: `PCGASMSetOverlap()`, `PCGASMGetSubKSP()`, 1797 `PCGASMCreateSubdomains2D()`, `PCGASMSetSubdomains()`, `PCGASMGetSubdomains()` 1798 @*/ 1799 PetscErrorCode PCGASMGetSubmatrices(PC pc, PetscInt *n, Mat *mat[]) { 1800 PC_GASM *osm; 1801 PetscBool match; 1802 1803 PetscFunctionBegin; 1804 PetscValidHeaderSpecific(pc, PC_CLASSID, 1); 1805 PetscValidIntPointer(n, 2); 1806 if (mat) PetscValidPointer(mat, 3); 1807 PetscCheck(pc->setupcalled, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_WRONGSTATE, "Must call after KSPSetUp() or PCSetUp()."); 1808 PetscCall(PetscObjectTypeCompare((PetscObject)pc, PCGASM, &match)); 1809 PetscCheck(match, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_WRONG, "Expected %s, got %s instead", PCGASM, ((PetscObject)pc)->type_name); 1810 osm = (PC_GASM *)pc->data; 1811 if (n) *n = osm->n; 1812 if (mat) *mat = osm->pmat; 1813 PetscFunctionReturn(0); 1814 } 1815 1816 /*@ 1817 PCGASMSetUseDMSubdomains - Indicates whether to use DMCreateDomainDecomposition() to define the subdomains, whenever possible. 1818 Logically Collective 1819 1820 Input Parameters: 1821 + pc - the preconditioner 1822 - flg - boolean indicating whether to use subdomains defined by the DM 1823 1824 Options Database Key: 1825 . -pc_gasm_dm_subdomains -pc_gasm_overlap -pc_gasm_total_subdomains 1826 1827 Level: intermediate 1828 1829 Notes: 1830 PCGASMSetSubdomains(), PCGASMSetTotalSubdomains() or PCGASMSetOverlap() take precedence over PCGASMSetUseDMSubdomains(), 1831 so setting PCGASMSetSubdomains() with nontrivial subdomain ISs or any of PCGASMSetTotalSubdomains() and PCGASMSetOverlap() 1832 automatically turns the latter off. 1833 1834 .seealso: `PCGASMGetUseDMSubdomains()`, `PCGASMSetSubdomains()`, `PCGASMSetOverlap()` 1835 `PCGASMCreateSubdomains2D()` 1836 @*/ 1837 PetscErrorCode PCGASMSetUseDMSubdomains(PC pc, PetscBool flg) { 1838 PC_GASM *osm = (PC_GASM *)pc->data; 1839 PetscBool match; 1840 1841 PetscFunctionBegin; 1842 PetscValidHeaderSpecific(pc, PC_CLASSID, 1); 1843 PetscValidLogicalCollectiveBool(pc, flg, 2); 1844 PetscCheck(!pc->setupcalled, ((PetscObject)pc)->comm, PETSC_ERR_ARG_WRONGSTATE, "Not for a setup PC."); 1845 PetscCall(PetscObjectTypeCompare((PetscObject)pc, PCGASM, &match)); 1846 if (match) { 1847 if (!osm->user_subdomains && osm->N == PETSC_DETERMINE && osm->overlap < 0) { osm->dm_subdomains = flg; } 1848 } 1849 PetscFunctionReturn(0); 1850 } 1851 1852 /*@ 1853 PCGASMGetUseDMSubdomains - Returns flag indicating whether to use DMCreateDomainDecomposition() to define the subdomains, whenever possible. 1854 Not Collective 1855 1856 Input Parameter: 1857 . pc - the preconditioner 1858 1859 Output Parameter: 1860 . flg - boolean indicating whether to use subdomains defined by the DM 1861 1862 Level: intermediate 1863 1864 .seealso: `PCGASMSetUseDMSubdomains()`, `PCGASMSetOverlap()` 1865 `PCGASMCreateSubdomains2D()` 1866 @*/ 1867 PetscErrorCode PCGASMGetUseDMSubdomains(PC pc, PetscBool *flg) { 1868 PC_GASM *osm = (PC_GASM *)pc->data; 1869 PetscBool match; 1870 1871 PetscFunctionBegin; 1872 PetscValidHeaderSpecific(pc, PC_CLASSID, 1); 1873 PetscValidBoolPointer(flg, 2); 1874 PetscCall(PetscObjectTypeCompare((PetscObject)pc, PCGASM, &match)); 1875 if (match) { 1876 if (flg) *flg = osm->dm_subdomains; 1877 } 1878 PetscFunctionReturn(0); 1879 } 1880