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