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