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 CHKERRQ(PetscMalloc2(osm->n,numbering,osm->n,permutation)); 47 CHKERRQ(PetscObjectsListGetGlobalNumbering(PetscObjectComm((PetscObject)pc),osm->n,(PetscObject*)osm->iis,NULL,*numbering)); 48 for (i = 0; i < osm->n; ++i) (*permutation)[i] = i; 49 CHKERRQ(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 CHKERRQ(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 CHKERRQ(PetscMalloc1(len, &cidx)); 73 CHKERRQ(PetscViewerStringOpen(PETSC_COMM_SELF, cidx, len, &sviewer)); 74 #undef len 75 CHKERRQ(ISGetIndices(osm->iis[i], &idx)); 76 for (j = 0; j < nidx; ++j) { 77 CHKERRQ(PetscViewerStringSPrintf(sviewer, "%D ", idx[j])); 78 } 79 CHKERRQ(ISRestoreIndices(osm->iis[i],&idx)); 80 CHKERRQ(PetscViewerDestroy(&sviewer)); 81 CHKERRQ(PetscViewerASCIIPrintf(viewer, "Inner subdomain:\n")); 82 CHKERRQ(PetscViewerFlush(viewer)); 83 CHKERRQ(PetscViewerASCIIPushSynchronized(viewer)); 84 CHKERRQ(PetscViewerASCIISynchronizedPrintf(viewer, "%s", cidx)); 85 CHKERRQ(PetscViewerFlush(viewer)); 86 CHKERRQ(PetscViewerASCIIPopSynchronized(viewer)); 87 CHKERRQ(PetscViewerASCIIPrintf(viewer, "\n")); 88 CHKERRQ(PetscViewerFlush(viewer)); 89 CHKERRQ(PetscFree(cidx)); 90 /* Outer subdomains. */ 91 CHKERRQ(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 CHKERRQ(PetscMalloc1(len, &cidx)); 100 CHKERRQ(PetscViewerStringOpen(PETSC_COMM_SELF, cidx, len, &sviewer)); 101 #undef len 102 CHKERRQ(ISGetIndices(osm->ois[i], &idx)); 103 for (j = 0; j < nidx; ++j) { 104 CHKERRQ(PetscViewerStringSPrintf(sviewer,"%D ", idx[j])); 105 } 106 CHKERRQ(PetscViewerDestroy(&sviewer)); 107 CHKERRQ(ISRestoreIndices(osm->ois[i],&idx)); 108 CHKERRQ(PetscViewerASCIIPrintf(viewer, "Outer subdomain:\n")); 109 CHKERRQ(PetscViewerFlush(viewer)); 110 CHKERRQ(PetscViewerASCIIPushSynchronized(viewer)); 111 CHKERRQ(PetscViewerASCIISynchronizedPrintf(viewer, "%s", cidx)); 112 CHKERRQ(PetscViewerFlush(viewer)); 113 CHKERRQ(PetscViewerASCIIPopSynchronized(viewer)); 114 CHKERRQ(PetscViewerASCIIPrintf(viewer, "\n")); 115 CHKERRQ(PetscViewerFlush(viewer)); 116 CHKERRQ(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 CHKERRQ(PCGetOptionsPrefix(pc,&prefix)); 132 CHKERRQ(PetscOptionsHasName(NULL,prefix,"-pc_gasm_print_subdomains",&found)); 133 if (!found) PetscFunctionReturn(0); 134 CHKERRQ(PetscOptionsGetString(NULL,prefix,"-pc_gasm_print_subdomains",fname,sizeof(fname),&found)); 135 if (!found) CHKERRQ(PetscStrcpy(fname,"stdout")); 136 CHKERRQ(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 CHKERRQ(PetscObjectName((PetscObject)viewer)); 142 l = 0; 143 CHKERRQ(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 CHKERRQ(PetscViewerGetSubViewer(viewer,((PetscObject)osm->ois[d])->comm, &sviewer)); 150 CHKERRQ(PCGASMSubdomainView_Private(pc,d,sviewer)); 151 CHKERRQ(PetscViewerRestoreSubViewer(viewer,((PetscObject)osm->ois[d])->comm, &sviewer)); 152 ++l; 153 } 154 } 155 CHKERRMPI(MPI_Barrier(PetscObjectComm((PetscObject)pc))); 156 } 157 CHKERRQ(PetscFree2(numbering,permutation)); 158 CHKERRQ(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 CHKERRMPI(MPI_Comm_size(PetscObjectComm((PetscObject)pc), &size)); 178 CHKERRMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)pc), &rank)); 179 180 if (osm->overlap >= 0) { 181 CHKERRQ(PetscSNPrintf(overlap,sizeof(overlap),"requested amount of overlap = %D",osm->overlap)); 182 } 183 if (osm->N != PETSC_DETERMINE) { 184 CHKERRQ(PetscSNPrintf(gsubdomains, sizeof(gsubdomains), "total number of subdomains = %D",osm->N)); 185 } 186 if (osm->nmax != PETSC_DETERMINE) { 187 CHKERRQ(PetscSNPrintf(msubdomains,sizeof(msubdomains),"max number of local subdomains = %D",osm->nmax)); 188 } 189 190 CHKERRQ(PCGetOptionsPrefix(pc,&prefix)); 191 CHKERRQ(PetscOptionsGetBool(NULL,prefix,"-pc_gasm_view_subdomains",&view_subdomains,NULL)); 192 193 CHKERRQ(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 CHKERRQ(PetscObjectName((PetscObject)viewer)); 201 CHKERRQ(PetscViewerASCIIPrintf(viewer," Restriction/interpolation type: %s\n",PCGASMTypes[osm->type])); 202 CHKERRQ(PetscViewerASCIIPrintf(viewer," %s\n",overlap)); 203 CHKERRQ(PetscViewerASCIIPrintf(viewer," %s\n",gsubdomains)); 204 CHKERRQ(PetscViewerASCIIPrintf(viewer," %s\n",msubdomains)); 205 CHKERRQ(PetscViewerASCIIPushSynchronized(viewer)); 206 CHKERRQ(PetscViewerASCIISynchronizedPrintf(viewer," [%d|%d] number of locally-supported subdomains = %D\n",rank,size,osm->n)); 207 CHKERRQ(PetscViewerFlush(viewer)); 208 CHKERRQ(PetscViewerASCIIPopSynchronized(viewer)); 209 /* Cannot take advantage of osm->same_subdomain_solvers without a global numbering of subdomains. */ 210 CHKERRQ(PetscViewerASCIIPrintf(viewer," Subdomain solver info is as follows:\n")); 211 CHKERRQ(PetscViewerASCIIPushTab(viewer)); 212 CHKERRQ(PetscViewerASCIIPrintf(viewer," - - - - - - - - - - - - - - - - - -\n")); 213 /* Make sure that everybody waits for the banner to be printed. */ 214 CHKERRMPI(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 CHKERRQ(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 CHKERRMPI(MPI_Comm_size(((PetscObject)osm->ois[d])->comm, &ssize)); 224 CHKERRMPI(MPI_Comm_rank(((PetscObject)osm->ois[d])->comm, &srank)); 225 CHKERRQ(PetscViewerGetSubViewer(viewer,((PetscObject)osm->ois[d])->comm, &sviewer)); 226 CHKERRQ(ISGetLocalSize(osm->ois[d],&bsz)); 227 CHKERRQ(PetscViewerASCIISynchronizedPrintf(sviewer," [%d|%d] (subcomm [%d|%d]) local subdomain number %D, local size = %D\n",rank,size,srank,ssize,d,bsz)); 228 CHKERRQ(PetscViewerFlush(sviewer)); 229 if (view_subdomains) { 230 CHKERRQ(PCGASMSubdomainView_Private(pc,d,sviewer)); 231 } 232 if (!pc->setupcalled) { 233 CHKERRQ(PetscViewerASCIIPrintf(sviewer, " Solver not set up yet: PCSetUp() not yet called\n")); 234 } else { 235 CHKERRQ(KSPView(osm->ksp[d],sviewer)); 236 } 237 CHKERRQ(PetscViewerASCIIPrintf(sviewer," - - - - - - - - - - - - - - - - - -\n")); 238 CHKERRQ(PetscViewerFlush(sviewer)); 239 CHKERRQ(PetscViewerRestoreSubViewer(viewer,((PetscObject)osm->ois[d])->comm, &sviewer)); 240 ++l; 241 } 242 } 243 CHKERRMPI(MPI_Barrier(PetscObjectComm((PetscObject)pc))); 244 } 245 CHKERRQ(PetscFree2(numbering,permutation)); 246 CHKERRQ(PetscViewerASCIIPopTab(viewer)); 247 CHKERRQ(PetscViewerFlush(viewer)); 248 /* this line is needed to match the extra PetscViewerASCIIPushSynchronized() in PetscViewerGetSubViewer() */ 249 CHKERRQ(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 CHKERRQ(PetscObjectGetComm((PetscObject)pc,&comm)); 268 CHKERRMPI(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 PetscCheckFalse(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 CHKERRQ(MatPartitioningCreate(comm,&part)); 279 CHKERRQ(MatPartitioningSetAdjacency(part,pc->pmat)); 280 CHKERRQ(MatPartitioningSetType(part,MATPARTITIONINGHIERARCH)); 281 CHKERRQ(MatPartitioningHierarchicalSetNcoarseparts(part,osm->N)); 282 CHKERRQ(MatPartitioningHierarchicalSetNfineparts(part,nlocalsubdomains)); 283 CHKERRQ(MatPartitioningSetFromOptions(part)); 284 /* get new processor owner number of each vertex */ 285 CHKERRQ(MatPartitioningApply(part,&partitioning)); 286 CHKERRQ(ISBuildTwoSided(partitioning,NULL,&fromrows)); 287 CHKERRQ(ISPartitioningToNumbering(partitioning,&isn)); 288 CHKERRQ(ISDestroy(&isn)); 289 CHKERRQ(ISGetLocalSize(fromrows,&fromrows_localsize)); 290 CHKERRQ(MatPartitioningDestroy(&part)); 291 CHKERRQ(MatCreateVecs(pc->pmat,&outervec,NULL)); 292 CHKERRQ(VecCreateMPI(comm,fromrows_localsize,PETSC_DETERMINE,&(osm->pcx))); 293 CHKERRQ(VecDuplicate(osm->pcx,&(osm->pcy))); 294 CHKERRQ(VecScatterCreate(osm->pcx,NULL,outervec,fromrows,&(osm->pctoouter))); 295 CHKERRQ(MatCreateSubMatrix(pc->pmat,fromrows,fromrows,MAT_INITIAL_MATRIX,&(osm->permutationP))); 296 CHKERRQ(PetscObjectReference((PetscObject)fromrows)); 297 osm->permutationIS = fromrows; 298 osm->pcmat = pc->pmat; 299 CHKERRQ(PetscObjectReference((PetscObject)osm->permutationP)); 300 pc->pmat = osm->permutationP; 301 CHKERRQ(VecDestroy(&outervec)); 302 CHKERRQ(ISDestroy(&fromrows)); 303 CHKERRQ(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 CHKERRMPI(MPI_Comm_size(PetscObjectComm((PetscObject)pc),&size)); 333 CHKERRMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)pc),&rank)); 334 if (!pc->setupcalled) { 335 /* use a hierarchical partitioning */ 336 if (osm->hierarchicalpartitioning) { 337 CHKERRQ(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 CHKERRQ(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 CHKERRQ(DMCreateDomainDecomposition(pc->dm, &num_subdomains, &subdomain_names, &inner_subdomain_is, &outer_subdomain_is, &subdomain_dm)); 348 if (num_subdomains) { 349 CHKERRQ(PCGASMSetSubdomains(pc, num_subdomains, inner_subdomain_is, outer_subdomain_is)); 350 } 351 for (d = 0; d < num_subdomains; ++d) { 352 if (inner_subdomain_is) CHKERRQ(ISDestroy(&inner_subdomain_is[d])); 353 if (outer_subdomain_is) CHKERRQ(ISDestroy(&outer_subdomain_is[d])); 354 } 355 CHKERRQ(PetscFree(inner_subdomain_is)); 356 CHKERRQ(PetscFree(outer_subdomain_is)); 357 } else { 358 /* still no subdomains; use one per processor */ 359 osm->nmax = osm->n = 1; 360 CHKERRMPI(MPI_Comm_size(PetscObjectComm((PetscObject)pc),&size)); 361 osm->N = size; 362 CHKERRQ(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 CHKERRQ(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 CHKERRQ(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 CHKERRQ(ISDuplicate(osm->iis[i],(osm->ois)+i)); 382 CHKERRQ(ISCopy(osm->iis[i],osm->ois[i])); 383 } else { 384 CHKERRQ(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 CHKERRQ(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 CHKERRMPI(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 CHKERRQ(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 CHKERRQ(ISSort(osm->ois[i])); 410 CHKERRQ(ISSort(osm->iis[i])); 411 } 412 } 413 CHKERRQ(PCGetOptionsPrefix(pc,&prefix)); 414 CHKERRQ(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 CHKERRQ(ISGetLocalSize(osm->ois[i],&oni)); 423 on += oni; 424 } 425 CHKERRQ(PetscMalloc1(on, &oidx)); 426 on = 0; 427 /* Merge local indices together */ 428 for (i=0; i<osm->n; i++) { 429 CHKERRQ(ISGetLocalSize(osm->ois[i],&oni)); 430 CHKERRQ(ISGetIndices(osm->ois[i],&oidxi)); 431 CHKERRQ(PetscArraycpy(oidx+on,oidxi,oni)); 432 CHKERRQ(ISRestoreIndices(osm->ois[i],&oidxi)); 433 on += oni; 434 } 435 CHKERRQ(ISCreateGeneral(((PetscObject)(pc))->comm,on,oidx,PETSC_OWN_POINTER,&gois)); 436 nTotalInnerIndices = 0; 437 for (i=0; i<osm->n; i++) { 438 CHKERRQ(ISGetLocalSize(osm->iis[i],&nInnerIndices)); 439 nTotalInnerIndices += nInnerIndices; 440 } 441 CHKERRQ(VecCreateMPI(((PetscObject)(pc))->comm,nTotalInnerIndices,PETSC_DETERMINE,&x)); 442 CHKERRQ(VecDuplicate(x,&y)); 443 444 CHKERRQ(VecCreateMPI(PetscObjectComm((PetscObject)pc),on,PETSC_DECIDE,&osm->gx)); 445 CHKERRQ(VecDuplicate(osm->gx,&osm->gy)); 446 CHKERRQ(VecGetOwnershipRange(osm->gx, &gostart, NULL)); 447 CHKERRQ(ISCreateStride(PetscObjectComm((PetscObject)pc),on,gostart,1, &goid)); 448 /* gois might indices not on local */ 449 CHKERRQ(VecScatterCreate(x,gois,osm->gx,goid, &(osm->gorestriction))); 450 CHKERRQ(PetscMalloc1(osm->n,&numbering)); 451 CHKERRQ(PetscObjectsListGetGlobalNumbering(PetscObjectComm((PetscObject)pc),osm->n,(PetscObject*)osm->ois,NULL,numbering)); 452 CHKERRQ(VecDestroy(&x)); 453 CHKERRQ(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 CHKERRQ(ISGetLocalSize(osm->ois[i],&oni)); 469 on += oni; 470 } 471 CHKERRQ(PetscMalloc1(on, &iidx)); 472 CHKERRQ(PetscMalloc1(on, &ioidx)); 473 CHKERRQ(VecGetArray(y,&array)); 474 /* set communicator id to determine where overlap is */ 475 in = 0; 476 for (i=0; i<osm->n; i++) { 477 CHKERRQ(ISGetLocalSize(osm->iis[i],&ini)); 478 for (k = 0; k < ini; ++k) { 479 array[in+k] = numbering[i]; 480 } 481 in += ini; 482 } 483 CHKERRQ(VecRestoreArray(y,&array)); 484 CHKERRQ(VecScatterBegin(osm->gorestriction,y,osm->gy,INSERT_VALUES,SCATTER_FORWARD)); 485 CHKERRQ(VecScatterEnd(osm->gorestriction,y,osm->gy,INSERT_VALUES,SCATTER_FORWARD)); 486 CHKERRQ(VecGetOwnershipRange(osm->gy,&gostart, NULL)); 487 CHKERRQ(VecGetArray(osm->gy,&array)); 488 on = 0; 489 in = 0; 490 for (i=0; i<osm->n; i++) { 491 CHKERRQ(ISGetLocalSize(osm->ois[i],&oni)); 492 CHKERRQ(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 CHKERRQ(ISRestoreIndices(osm->ois[i], &indices)); 500 on += oni; 501 } 502 CHKERRQ(VecRestoreArray(osm->gy,&array)); 503 CHKERRQ(ISCreateGeneral(PetscObjectComm((PetscObject)pc),in,iidx,PETSC_OWN_POINTER,&giis)); 504 CHKERRQ(ISCreateGeneral(PetscObjectComm((PetscObject)pc),in,ioidx,PETSC_OWN_POINTER,&giois)); 505 CHKERRQ(VecScatterCreate(y,giis,osm->gy,giois,&osm->girestriction)); 506 CHKERRQ(VecDestroy(&y)); 507 CHKERRQ(ISDestroy(&giis)); 508 CHKERRQ(ISDestroy(&giois)); 509 } 510 CHKERRQ(ISDestroy(&goid)); 511 CHKERRQ(PetscFree(numbering)); 512 513 /* Create the subdomain work vectors. */ 514 CHKERRQ(PetscMalloc1(osm->n,&osm->x)); 515 CHKERRQ(PetscMalloc1(osm->n,&osm->y)); 516 CHKERRQ(VecGetArray(osm->gx, &gxarray)); 517 CHKERRQ(VecGetArray(osm->gy, &gyarray)); 518 for (i=0, on=0; i<osm->n; ++i, on += oni) { 519 PetscInt oNi; 520 CHKERRQ(ISGetLocalSize(osm->ois[i],&oni)); 521 /* on a sub communicator */ 522 CHKERRQ(ISGetSize(osm->ois[i],&oNi)); 523 CHKERRQ(VecCreateMPIWithArray(((PetscObject)(osm->ois[i]))->comm,1,oni,oNi,gxarray+on,&osm->x[i])); 524 CHKERRQ(VecCreateMPIWithArray(((PetscObject)(osm->ois[i]))->comm,1,oni,oNi,gyarray+on,&osm->y[i])); 525 } 526 CHKERRQ(VecRestoreArray(osm->gx, &gxarray)); 527 CHKERRQ(VecRestoreArray(osm->gy, &gyarray)); 528 /* Create the subdomain solvers */ 529 CHKERRQ(PetscMalloc1(osm->n,&osm->ksp)); 530 for (i=0; i<osm->n; i++) { 531 char subprefix[PETSC_MAX_PATH_LEN+1]; 532 CHKERRQ(KSPCreate(((PetscObject)(osm->ois[i]))->comm,&ksp)); 533 CHKERRQ(KSPSetErrorIfNotConverged(ksp,pc->erroriffailure)); 534 CHKERRQ(PetscLogObjectParent((PetscObject)pc,(PetscObject)ksp)); 535 CHKERRQ(PetscObjectIncrementTabLevel((PetscObject)ksp,(PetscObject)pc,1)); 536 CHKERRQ(KSPSetType(ksp,KSPPREONLY)); 537 CHKERRQ(KSPGetPC(ksp,&subpc)); /* Why do we need this here? */ 538 if (subdomain_dm) { 539 CHKERRQ(KSPSetDM(ksp,subdomain_dm[i])); 540 CHKERRQ(DMDestroy(subdomain_dm+i)); 541 } 542 CHKERRQ(PCGetOptionsPrefix(pc,&prefix)); 543 CHKERRQ(KSPSetOptionsPrefix(ksp,prefix)); 544 if (subdomain_names && subdomain_names[i]) { 545 CHKERRQ(PetscSNPrintf(subprefix,PETSC_MAX_PATH_LEN,"sub_%s_",subdomain_names[i])); 546 CHKERRQ(KSPAppendOptionsPrefix(ksp,subprefix)); 547 CHKERRQ(PetscFree(subdomain_names[i])); 548 } 549 CHKERRQ(KSPAppendOptionsPrefix(ksp,"sub_")); 550 osm->ksp[i] = ksp; 551 } 552 CHKERRQ(PetscFree(subdomain_dm)); 553 CHKERRQ(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 CHKERRQ(MatDestroyMatrices(osm->n,&osm->pmat)); 561 scall = MAT_INITIAL_MATRIX; 562 } 563 if (osm->permutationIS) { 564 CHKERRQ(MatCreateSubMatrix(pc->pmat,osm->permutationIS,osm->permutationIS,scall,&osm->permutationP)); 565 CHKERRQ(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 CHKERRQ(MatCreateSubMatricesMPI(pc->pmat,osm->n,osm->ois,osm->ois,scall,&osm->pmat)); 576 } else { 577 CHKERRQ(MatCreateSubMatrices(pc->pmat,osm->n,osm->ois,osm->ois,scall,&osm->pmat)); 578 } 579 if (scall == MAT_INITIAL_MATRIX) { 580 CHKERRQ(PetscObjectGetOptionsPrefix((PetscObject)pc->pmat,&pprefix)); 581 for (i=0; i<osm->n; i++) { 582 CHKERRQ(PetscLogObjectParent((PetscObject)pc,(PetscObject)osm->pmat[i])); 583 CHKERRQ(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 CHKERRQ(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 CHKERRQ(KSPSetOperators(osm->ksp[i],osm->pmat[i],osm->pmat[i])); 596 CHKERRQ(KSPGetOptionsPrefix(osm->ksp[i],&prefix)); 597 CHKERRQ(MatSetOptionsPrefix(osm->pmat[i],prefix)); 598 if (!pc->setupcalled) { 599 CHKERRQ(KSPSetFromOptions(osm->ksp[i])); 600 } 601 } 602 if (osm->pcmat) { 603 CHKERRQ(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 CHKERRQ(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 CHKERRQ(VecScatterBegin(osm->pctoouter,xin,osm->pcx,INSERT_VALUES,SCATTER_REVERSE)); 632 CHKERRQ(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 CHKERRQ(VecZeroEntries(osm->gx)); 646 CHKERRQ(VecScatterBegin(osm->girestriction,x,osm->gx,INSERT_VALUES,forward)); 647 } else { 648 CHKERRQ(VecScatterBegin(osm->gorestriction,x,osm->gx,INSERT_VALUES,forward)); 649 } 650 CHKERRQ(VecZeroEntries(osm->gy)); 651 if (!(osm->type & PC_GASM_RESTRICT)) { 652 CHKERRQ(VecScatterEnd(osm->girestriction,x,osm->gx,INSERT_VALUES,forward)); 653 } else { 654 CHKERRQ(VecScatterEnd(osm->gorestriction,x,osm->gx,INSERT_VALUES,forward)); 655 } 656 /* do the subdomain solves */ 657 for (i=0; i<osm->n; ++i) { 658 CHKERRQ(KSPSolve(osm->ksp[i],osm->x[i],osm->y[i])); 659 CHKERRQ(KSPCheckSolve(osm->ksp[i],pc,osm->y[i])); 660 } 661 /* do we need to zero y? */ 662 CHKERRQ(VecZeroEntries(y)); 663 if (!(osm->type & PC_GASM_INTERPOLATE)) { 664 CHKERRQ(VecScatterBegin(osm->girestriction,osm->gy,y,ADD_VALUES,reverse)); 665 CHKERRQ(VecScatterEnd(osm->girestriction,osm->gy,y,ADD_VALUES,reverse)); 666 } else { 667 CHKERRQ(VecScatterBegin(osm->gorestriction,osm->gy,y,ADD_VALUES,reverse)); 668 CHKERRQ(VecScatterEnd(osm->gorestriction,osm->gy,y,ADD_VALUES,reverse)); 669 } 670 if (osm->pctoouter) { 671 CHKERRQ(VecScatterBegin(osm->pctoouter,y,yout,INSERT_VALUES,SCATTER_FORWARD)); 672 CHKERRQ(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 CHKERRQ(MatGetSize(Xin,NULL,&N)); 688 if (osm->pctoouter) { 689 CHKERRQ(VecGetLocalSize(osm->pcx,&m)); 690 CHKERRQ(VecGetSize(osm->pcx,&M)); 691 CHKERRQ(MatCreateDense(PetscObjectComm((PetscObject)osm->ois[0]),m,PETSC_DECIDE,M,N,NULL,&O)); 692 for (i = 0; i < N; ++i) { 693 CHKERRQ(MatDenseGetColumnVecRead(Xin,i,&x)); 694 CHKERRQ(MatDenseGetColumnVecWrite(O,i,&y)); 695 CHKERRQ(VecScatterBegin(osm->pctoouter,x,y,INSERT_VALUES,SCATTER_REVERSE)); 696 CHKERRQ(VecScatterEnd(osm->pctoouter,x,y,INSERT_VALUES,SCATTER_REVERSE)); 697 CHKERRQ(MatDenseRestoreColumnVecWrite(O,i,&y)); 698 CHKERRQ(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 CHKERRQ(VecGetLocalSize(osm->x[0],&m)); 710 CHKERRQ(VecGetSize(osm->x[0],&M)); 711 CHKERRQ(MatCreateDense(PetscObjectComm((PetscObject)osm->ois[0]),m,PETSC_DECIDE,M,N,NULL,&Z)); 712 for (i = 0; i < N; ++i) { 713 CHKERRQ(MatDenseGetColumnVecRead(X,i,&x)); 714 CHKERRQ(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 CHKERRQ(VecZeroEntries(y)); 718 CHKERRQ(VecScatterBegin(osm->girestriction,x,y,INSERT_VALUES,forward)); 719 CHKERRQ(VecScatterEnd(osm->girestriction,x,y,INSERT_VALUES,forward)); 720 } else { 721 CHKERRQ(VecScatterBegin(osm->gorestriction,x,y,INSERT_VALUES,forward)); 722 CHKERRQ(VecScatterEnd(osm->gorestriction,x,y,INSERT_VALUES,forward)); 723 } 724 CHKERRQ(MatDenseRestoreColumnVecWrite(Z,i,&y)); 725 CHKERRQ(MatDenseRestoreColumnVecRead(X,i,&x)); 726 } 727 CHKERRQ(MatCreateDense(PetscObjectComm((PetscObject)osm->ois[0]),m,PETSC_DECIDE,M,N,NULL,&W)); 728 CHKERRQ(MatSetOption(Z,MAT_NO_OFF_PROC_ENTRIES,PETSC_TRUE)); 729 CHKERRQ(MatAssemblyBegin(Z,MAT_FINAL_ASSEMBLY)); 730 CHKERRQ(MatAssemblyEnd(Z,MAT_FINAL_ASSEMBLY)); 731 /* do the subdomain solve */ 732 CHKERRQ(KSPMatSolve(osm->ksp[0],Z,W)); 733 CHKERRQ(KSPCheckSolve(osm->ksp[0],pc,NULL)); 734 CHKERRQ(MatDestroy(&Z)); 735 /* do we need to zero y? */ 736 CHKERRQ(MatZeroEntries(Y)); 737 for (i = 0; i < N; ++i) { 738 CHKERRQ(MatDenseGetColumnVecWrite(Y,i,&y)); 739 CHKERRQ(MatDenseGetColumnVecRead(W,i,&x)); 740 if (!(osm->type & PC_GASM_INTERPOLATE)) { 741 CHKERRQ(VecScatterBegin(osm->girestriction,x,y,ADD_VALUES,reverse)); 742 CHKERRQ(VecScatterEnd(osm->girestriction,x,y,ADD_VALUES,reverse)); 743 } else { 744 CHKERRQ(VecScatterBegin(osm->gorestriction,x,y,ADD_VALUES,reverse)); 745 CHKERRQ(VecScatterEnd(osm->gorestriction,x,y,ADD_VALUES,reverse)); 746 } 747 CHKERRQ(MatDenseRestoreColumnVecRead(W,i,&x)); 748 if (osm->pctoouter) { 749 CHKERRQ(MatDenseGetColumnVecWrite(Yout,i,&x)); 750 CHKERRQ(VecScatterBegin(osm->pctoouter,y,x,INSERT_VALUES,SCATTER_FORWARD)); 751 CHKERRQ(VecScatterEnd(osm->pctoouter,y,x,INSERT_VALUES,SCATTER_FORWARD)); 752 CHKERRQ(MatDenseRestoreColumnVecRead(Yout,i,&x)); 753 } 754 CHKERRQ(MatDenseRestoreColumnVecWrite(Y,i,&y)); 755 } 756 CHKERRQ(MatDestroy(&W)); 757 CHKERRQ(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 CHKERRQ(VecScatterBegin(osm->pctoouter,xin,osm->pcx,INSERT_VALUES,SCATTER_REVERSE)); 771 CHKERRQ(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 CHKERRQ(VecZeroEntries(osm->gx)); 788 CHKERRQ(VecScatterBegin(osm->girestriction,x,osm->gx,INSERT_VALUES,forward)); 789 } else { 790 CHKERRQ(VecScatterBegin(osm->gorestriction,x,osm->gx,INSERT_VALUES,forward)); 791 } 792 CHKERRQ(VecZeroEntries(osm->gy)); 793 if (!(osm->type & PC_GASM_INTERPOLATE)) { 794 CHKERRQ(VecScatterEnd(osm->girestriction,x,osm->gx,INSERT_VALUES,forward)); 795 } else { 796 CHKERRQ(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 CHKERRQ(KSPSolveTranspose(osm->ksp[i],osm->x[i],osm->y[i])); 801 CHKERRQ(KSPCheckSolve(osm->ksp[i],pc,osm->y[i])); 802 } 803 CHKERRQ(VecZeroEntries(y)); 804 if (!(osm->type & PC_GASM_RESTRICT)) { 805 CHKERRQ(VecScatterBegin(osm->girestriction,osm->gy,y,ADD_VALUES,reverse)); 806 CHKERRQ(VecScatterEnd(osm->girestriction,osm->gy,y,ADD_VALUES,reverse)); 807 } else { 808 CHKERRQ(VecScatterBegin(osm->gorestriction,osm->gy,y,ADD_VALUES,reverse)); 809 CHKERRQ(VecScatterEnd(osm->gorestriction,osm->gy,y,ADD_VALUES,reverse)); 810 } 811 if (osm->pctoouter) { 812 CHKERRQ(VecScatterBegin(osm->pctoouter,y,yout,INSERT_VALUES,SCATTER_FORWARD)); 813 CHKERRQ(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 CHKERRQ(KSPReset(osm->ksp[i])); 827 } 828 } 829 if (osm->pmat) { 830 if (osm->n > 0) { 831 PetscMPIInt size; 832 CHKERRMPI(MPI_Comm_size(PetscObjectComm((PetscObject)pc),&size)); 833 if (size > 1) { 834 /* osm->pmat is created by MatCreateSubMatricesMPI(), cannot use MatDestroySubMatrices() */ 835 CHKERRQ(MatDestroyMatrices(osm->n,&osm->pmat)); 836 } else { 837 CHKERRQ(MatDestroySubMatrices(osm->n,&osm->pmat)); 838 } 839 } 840 } 841 if (osm->x) { 842 for (i=0; i<osm->n; i++) { 843 CHKERRQ(VecDestroy(&osm->x[i])); 844 CHKERRQ(VecDestroy(&osm->y[i])); 845 } 846 } 847 CHKERRQ(VecDestroy(&osm->gx)); 848 CHKERRQ(VecDestroy(&osm->gy)); 849 850 CHKERRQ(VecScatterDestroy(&osm->gorestriction)); 851 CHKERRQ(VecScatterDestroy(&osm->girestriction)); 852 if (!osm->user_subdomains) { 853 CHKERRQ(PCGASMDestroySubdomains(osm->n,&osm->ois,&osm->iis)); 854 osm->N = PETSC_DETERMINE; 855 osm->nmax = PETSC_DETERMINE; 856 } 857 if (osm->pctoouter) { 858 CHKERRQ(VecScatterDestroy(&(osm->pctoouter))); 859 } 860 if (osm->permutationIS) { 861 CHKERRQ(ISDestroy(&(osm->permutationIS))); 862 } 863 if (osm->pcx) { 864 CHKERRQ(VecDestroy(&(osm->pcx))); 865 } 866 if (osm->pcy) { 867 CHKERRQ(VecDestroy(&(osm->pcy))); 868 } 869 if (osm->permutationP) { 870 CHKERRQ(MatDestroy(&(osm->permutationP))); 871 } 872 if (osm->pcmat) { 873 CHKERRQ(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 CHKERRQ(PCReset_GASM(pc)); 885 /* PCReset will not destroy subdomains, if user_subdomains is true. */ 886 CHKERRQ(PCGASMDestroySubdomains(osm->n,&osm->ois,&osm->iis)); 887 if (osm->ksp) { 888 for (i=0; i<osm->n; i++) { 889 CHKERRQ(KSPDestroy(&osm->ksp[i])); 890 } 891 CHKERRQ(PetscFree(osm->ksp)); 892 } 893 CHKERRQ(PetscFree(osm->x)); 894 CHKERRQ(PetscFree(osm->y)); 895 CHKERRQ(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 CHKERRQ(PetscOptionsHead(PetscOptionsObject,"Generalized additive Schwarz options")); 908 CHKERRQ(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 CHKERRQ(PetscOptionsInt("-pc_gasm_total_subdomains","Total number of subdomains across communicator","PCGASMSetTotalSubdomains",osm->N,&blocks,&flg)); 910 if (flg) { 911 CHKERRQ(PCGASMSetTotalSubdomains(pc,blocks)); 912 } 913 CHKERRQ(PetscOptionsInt("-pc_gasm_overlap","Number of overlapping degrees of freedom","PCGASMSetOverlap",osm->overlap,&ovl,&flg)); 914 if (flg) { 915 CHKERRQ(PCGASMSetOverlap(pc,ovl)); 916 osm->dm_subdomains = PETSC_FALSE; 917 } 918 flg = PETSC_FALSE; 919 CHKERRQ(PetscOptionsEnum("-pc_gasm_type","Type of restriction/extension","PCGASMSetType",PCGASMTypes,(PetscEnum)osm->type,(PetscEnum*)&gasmtype,&flg)); 920 if (flg) CHKERRQ(PCGASMSetType(pc,gasmtype)); 921 CHKERRQ(PetscOptionsBool("-pc_gasm_use_hierachical_partitioning","use hierarchical partitioning",NULL,osm->hierarchicalpartitioning,&osm->hierarchicalpartitioning,&flg)); 922 CHKERRQ(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 CHKERRQ(PCGASMDestroySubdomains(osm->n,&osm->iis,&osm->ois)); 952 osm->ois = osm->iis = NULL; 953 954 CHKERRMPI(MPI_Comm_size(PetscObjectComm((PetscObject)pc),&size)); 955 CHKERRMPI(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 CHKERRQ(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 CHKERRQ(PetscMalloc1(n,&osm->ois)); 979 for (i=0; i<n; i++) { 980 CHKERRQ(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 CHKERRQ(PetscMalloc1(n,&osm->iis)); 993 for (i=0; i<n; i++) { 994 CHKERRQ(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 CHKERRQ(MatGetOwnershipRange(pc->pmat,&rstart,&rend)); 1007 CHKERRQ(PetscCalloc1(rend-rstart,&covered)); 1008 /* check if the current processor owns indices from others */ 1009 for (i=0; i<n; i++) { 1010 CHKERRQ(ISGetIndices(osm->iis[i],&indices)); 1011 CHKERRQ(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 CHKERRQ(ISRestoreIndices(osm->iis[i],&indices)); 1018 } 1019 /* check if we miss any indices */ 1020 for (i=rstart; i<rend; i++) { 1021 PetscCheckFalse(!covered[i-rstart],PETSC_COMM_SELF,PETSC_ERR_ARG_NULL,"local entity %d was not covered by inner subdomains",i); 1022 } 1023 CHKERRQ(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 CHKERRMPI(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 CHKERRQ(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 CHKERRQ(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 CHKERRQ(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 CHKERRQ(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 CHKERRQ(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 CHKERRQ(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 CHKERRQ(PetscObjectComposeFunction((PetscObject)pc,"PCGASMSetSubdomains_C",PCGASMSetSubdomains_GASM)); 1352 CHKERRQ(PetscObjectComposeFunction((PetscObject)pc,"PCGASMSetOverlap_C",PCGASMSetOverlap_GASM)); 1353 CHKERRQ(PetscObjectComposeFunction((PetscObject)pc,"PCGASMSetType_C",PCGASMSetType_GASM)); 1354 CHKERRQ(PetscObjectComposeFunction((PetscObject)pc,"PCGASMSetSortIndices_C",PCGASMSetSortIndices_GASM)); 1355 CHKERRQ(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 CHKERRQ(MatGetOptionsPrefix(A,&prefix)); 1373 CHKERRQ(MatGetOwnershipRange(A,&rstart,&rend)); 1374 CHKERRQ(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 CHKERRQ(MatHasOperation(A,MATOP_GET_DIAGONAL_BLOCK,&hasop)); 1379 if (hasop) { 1380 CHKERRQ(MatGetDiagonalBlock(A,&Ad)); 1381 } 1382 if (Ad) { 1383 CHKERRQ(PetscObjectBaseTypeCompare((PetscObject)Ad,MATSEQBAIJ,&isbaij)); 1384 if (!isbaij) CHKERRQ(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 CHKERRQ(MatPartitioningCreate(PETSC_COMM_SELF,&mpart)); 1390 CHKERRQ(PetscObjectSetOptionsPrefix((PetscObject)mpart,prefix)); 1391 CHKERRQ(MatPartitioningSetFromOptions(mpart)); 1392 CHKERRQ(PetscObjectTypeCompare((PetscObject)mpart,MATPARTITIONINGCURRENT,&match)); 1393 if (!match) { 1394 CHKERRQ(PetscObjectTypeCompare((PetscObject)mpart,MATPARTITIONINGSQUARE,&match)); 1395 } 1396 if (!match) { /* assume a "good" partitioner is available */ 1397 PetscInt na; 1398 const PetscInt *ia,*ja; 1399 CHKERRQ(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 CHKERRQ(PetscMalloc1(na+1,&iia)); 1419 CHKERRQ(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 CHKERRQ(MatCreateMPIAdj(PETSC_COMM_SELF,na,na,iia,jja,NULL,&adj)); 1434 CHKERRQ(MatPartitioningSetAdjacency(mpart,adj)); 1435 CHKERRQ(MatPartitioningSetNParts(mpart,nloc)); 1436 CHKERRQ(MatPartitioningApply(mpart,&ispart)); 1437 CHKERRQ(ISPartitioningToNumbering(ispart,&isnumb)); 1438 CHKERRQ(MatDestroy(&adj)); 1439 foundpart = PETSC_TRUE; 1440 } 1441 CHKERRQ(MatRestoreRowIJ(Ad,0,PETSC_TRUE,isbaij,&na,&ia,&ja,&done)); 1442 } 1443 CHKERRQ(MatPartitioningDestroy(&mpart)); 1444 } 1445 CHKERRQ(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 CHKERRQ(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 CHKERRQ(PetscMalloc1(nloc,&count)); 1466 CHKERRQ(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 CHKERRQ(ISGetLocalSize(isnumb,&nidx)); 1472 CHKERRQ(PetscMalloc1(nidx,&indices)); 1473 for (i=0; i<nidx; i++) indices[i] = i; /* needs to be initialized */ 1474 CHKERRQ(ISGetIndices(isnumb,&numbering)); 1475 CHKERRQ(PetscSortIntWithPermutation(nidx,numbering,indices)); 1476 CHKERRQ(ISRestoreIndices(isnumb,&numbering)); 1477 if (isbaij && bs > 1) { /* adjust for the block-aij case */ 1478 CHKERRQ(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 CHKERRQ(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 CHKERRQ(ISCreateGeneral(PETSC_COMM_SELF,count[i],&indices[start],PETSC_COPY_VALUES,&is[i])); 1492 CHKERRQ(ISSort(is[i])); 1493 start += count[i]; 1494 } 1495 1496 CHKERRQ(PetscFree(count)); 1497 CHKERRQ(PetscFree(indices)); 1498 CHKERRQ(ISDestroy(&isnumb)); 1499 CHKERRQ(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 CHKERRQ(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 CHKERRMPI(MPI_Comm_size(PetscObjectComm((PetscObject)A),&size)); 1546 if (N >= size) { 1547 *n = N/size + (N%size); 1548 CHKERRQ(PCGASMCreateLocalSubdomains(A,*n,iis)); 1549 } else { 1550 CHKERRQ(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 CHKERRQ(ISDestroy(&(*ois)[i])); 1588 } 1589 CHKERRQ(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 CHKERRQ(ISDestroy(&(*iis)[i])); 1598 } 1599 CHKERRQ(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 CHKERRQ(PetscObjectGetComm((PetscObject)pc, &comm)); 1675 CHKERRMPI(MPI_Comm_size(comm, &size)); 1676 CHKERRMPI(MPI_Comm_rank(comm, &rank)); 1677 CHKERRQ(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 CHKERRQ(PetscMalloc1(*nsub,is)); 1710 CHKERRQ(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 CHKERRMPI(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 CHKERRQ(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 CHKERRQ(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 CHKERRQ(ISCreateGeneral(subcomm,nidx,idx,PETSC_OWN_POINTER,(*xis)+s)); 1796 if (split) { 1797 CHKERRMPI(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 CHKERRQ(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 CHKERRQ(PetscMalloc1(osm->n, iis)); 1846 } 1847 if (ois) { 1848 CHKERRQ(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 CHKERRQ(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 CHKERRQ(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 CHKERRQ(PetscObjectTypeCompare((PetscObject)pc,PCGASM,&match)); 1961 if (match) { 1962 if (flg) *flg = osm->dm_subdomains; 1963 } 1964 PetscFunctionReturn(0); 1965 } 1966