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