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