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