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