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