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