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