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