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