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