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