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