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