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