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 use_dm_decomposition; /* 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->use_dm_decomposition && 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 = PetscOptionsInt("-pc_gasm_total_subdomains","Total number of subdomains across communicator","PCGASMSetTotalSubdomains",osm->n,&blocks,&flg);CHKERRQ(ierr); 665 if (flg) { 666 osm->create_local = PETSC_TRUE; 667 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); 668 ierr = PCGASMSetTotalSubdomains(pc,blocks,osm->create_local);CHKERRQ(ierr); 669 osm->use_dm_decomposition = PETSC_FALSE; 670 } 671 ierr = PetscOptionsInt("-pc_gasm_overlap","Number of overlapping degrees of freedom","PCGASMSetOverlap",osm->overlap,&ovl,&flg);CHKERRQ(ierr); 672 if (flg) { 673 ierr = PCGASMSetOverlap(pc,ovl);CHKERRQ(ierr); 674 osm->use_dm_decomposition = PETSC_FALSE; 675 } 676 flg = PETSC_FALSE; 677 ierr = PetscOptionsEnum("-pc_gasm_type","Type of restriction/extension","PCGASMSetType",PCGASMTypes,(PetscEnum)osm->type,(PetscEnum*)&gasmtype,&flg);CHKERRQ(ierr); 678 if (flg) {ierr = PCGASMSetType(pc,gasmtype);CHKERRQ(ierr); } 679 ierr = PetscOptionsTail();CHKERRQ(ierr); 680 PetscFunctionReturn(0); 681 } 682 683 /*------------------------------------------------------------------------------------*/ 684 685 EXTERN_C_BEGIN 686 #undef __FUNCT__ 687 #define __FUNCT__ "PCGASMSetSubdomains_GASM" 688 PetscErrorCode PCGASMSetSubdomains_GASM(PC pc,PetscInt n,IS iis[],IS ois[]) 689 { 690 PC_GASM *osm = (PC_GASM*)pc->data; 691 PetscErrorCode ierr; 692 PetscInt i; 693 694 PetscFunctionBegin; 695 if (n < 1) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Each process must have 1 or more subdomains, n = %D",n); 696 if (pc->setupcalled && (n != osm->n || iis || ois)) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONGSTATE,"PCGASMSetSubdomains() should be called before calling PCSetUp()."); 697 698 if (!pc->setupcalled) { 699 osm->n = n; 700 osm->ois = 0; 701 osm->iis = 0; 702 if (ois) { 703 for (i=0; i<n; i++) {ierr = PetscObjectReference((PetscObject)ois[i]);CHKERRQ(ierr);} 704 } 705 if (iis) { 706 for (i=0; i<n; i++) {ierr = PetscObjectReference((PetscObject)iis[i]);CHKERRQ(ierr);} 707 } 708 ierr = PCGASMDestroySubdomains(osm->n,osm->iis,osm->ois);CHKERRQ(ierr); 709 if (ois) { 710 ierr = PetscMalloc(n*sizeof(IS),&osm->ois);CHKERRQ(ierr); 711 for (i=0; i<n; i++) osm->ois[i] = ois[i]; 712 /* Flag indicating that the user has set outer subdomains, so PCGASM should not increase their size. */ 713 osm->overlap = -1; 714 if (!iis) { 715 ierr = PetscMalloc(n*sizeof(IS),&osm->iis);CHKERRQ(ierr); 716 for (i=0; i<n; i++) { 717 for (i=0; i<n; i++) {ierr = PetscObjectReference((PetscObject)ois[i]);CHKERRQ(ierr);} 718 osm->iis[i] = ois[i]; 719 } 720 } 721 } 722 if (iis) { 723 ierr = PetscMalloc(n*sizeof(IS),&osm->iis);CHKERRQ(ierr); 724 for (i=0; i<n; i++) osm->iis[i] = iis[i]; 725 if (!ois) { 726 ierr = PetscMalloc(n*sizeof(IS),&osm->ois);CHKERRQ(ierr); 727 for (i=0; i<n; i++) { 728 for (i=0; i<n; i++) { 729 ierr = PetscObjectReference((PetscObject)iis[i]);CHKERRQ(ierr); 730 osm->ois[i] = iis[i]; 731 } 732 } 733 if (osm->overlap > 0) { 734 /* Extend the "overlapping" regions by a number of steps */ 735 ierr = MatIncreaseOverlap(pc->pmat,osm->n,osm->ois,osm->overlap);CHKERRQ(ierr); 736 } 737 } 738 } 739 } 740 PetscFunctionReturn(0); 741 } 742 EXTERN_C_END 743 744 EXTERN_C_BEGIN 745 #undef __FUNCT__ 746 #define __FUNCT__ "PCGASMSetTotalSubdomains_GASM" 747 PetscErrorCode PCGASMSetTotalSubdomains_GASM(PC pc,PetscInt N, PetscBool create_local) 748 { 749 PC_GASM *osm = (PC_GASM*)pc->data; 750 PetscErrorCode ierr; 751 PetscMPIInt rank,size; 752 PetscInt n; 753 PetscInt Nmin, Nmax; 754 755 PetscFunctionBegin; 756 if (!create_local) SETERRQ(PetscObjectComm((PetscObject)pc), PETSC_ERR_SUP, "No suppor for autocreation of nonlocal subdomains."); 757 if (N < 1) SETERRQ1(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_OUTOFRANGE,"Total number of subdomains must be > 0, N = %D",N); 758 ierr = MPI_Allreduce(&N,&Nmin,1,MPIU_INT,MPIU_MIN,PetscObjectComm((PetscObject)pc));CHKERRQ(ierr); 759 ierr = MPI_Allreduce(&N,&Nmax,1,MPIU_INT,MPIU_MAX,PetscObjectComm((PetscObject)pc));CHKERRQ(ierr); 760 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); 761 762 osm->create_local = create_local; 763 /* 764 Split the subdomains equally among all processors 765 */ 766 ierr = MPI_Comm_rank(PetscObjectComm((PetscObject)pc),&rank);CHKERRQ(ierr); 767 ierr = MPI_Comm_size(PetscObjectComm((PetscObject)pc),&size);CHKERRQ(ierr); 768 n = N/size + ((N % size) > rank); 769 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); 770 if (pc->setupcalled && n != osm->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"PCGASMSetTotalSubdomains() should be called before PCSetUp()."); 771 if (!pc->setupcalled) { 772 ierr = PCGASMDestroySubdomains(osm->n,osm->iis,osm->ois);CHKERRQ(ierr); 773 osm->N = N; 774 osm->n = n; 775 osm->nmax = N/size + ((N%size) ? 1 : 0); 776 osm->ois = 0; 777 osm->iis = 0; 778 } 779 PetscFunctionReturn(0); 780 } 781 EXTERN_C_END 782 783 EXTERN_C_BEGIN 784 #undef __FUNCT__ 785 #define __FUNCT__ "PCGASMSetOverlap_GASM" 786 PetscErrorCode PCGASMSetOverlap_GASM(PC pc,PetscInt ovl) 787 { 788 PC_GASM *osm = (PC_GASM*)pc->data; 789 790 PetscFunctionBegin; 791 if (ovl < 0) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_OUTOFRANGE,"Negative overlap value requested"); 792 if (pc->setupcalled && ovl != osm->overlap) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONGSTATE,"PCGASMSetOverlap() should be called before PCSetUp()."); 793 if (!pc->setupcalled) osm->overlap = ovl; 794 PetscFunctionReturn(0); 795 } 796 EXTERN_C_END 797 798 EXTERN_C_BEGIN 799 #undef __FUNCT__ 800 #define __FUNCT__ "PCGASMSetType_GASM" 801 PetscErrorCode PCGASMSetType_GASM(PC pc,PCGASMType type) 802 { 803 PC_GASM *osm = (PC_GASM*)pc->data; 804 805 PetscFunctionBegin; 806 osm->type = type; 807 osm->type_set = PETSC_TRUE; 808 PetscFunctionReturn(0); 809 } 810 EXTERN_C_END 811 812 EXTERN_C_BEGIN 813 #undef __FUNCT__ 814 #define __FUNCT__ "PCGASMSetSortIndices_GASM" 815 PetscErrorCode PCGASMSetSortIndices_GASM(PC pc,PetscBool doSort) 816 { 817 PC_GASM *osm = (PC_GASM*)pc->data; 818 819 PetscFunctionBegin; 820 osm->sort_indices = doSort; 821 PetscFunctionReturn(0); 822 } 823 EXTERN_C_END 824 825 EXTERN_C_BEGIN 826 #undef __FUNCT__ 827 #define __FUNCT__ "PCGASMGetSubKSP_GASM" 828 /* 829 FIX: This routine might need to be modified once multiple ranks per subdomain are allowed. 830 In particular, it would upset the global subdomain number calculation. 831 */ 832 PetscErrorCode PCGASMGetSubKSP_GASM(PC pc,PetscInt *n,PetscInt *first,KSP **ksp) 833 { 834 PC_GASM *osm = (PC_GASM*)pc->data; 835 PetscErrorCode ierr; 836 837 PetscFunctionBegin; 838 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"); 839 840 if (n) *n = osm->n; 841 if (first) { 842 ierr = MPI_Scan(&osm->n,first,1,MPIU_INT,MPI_SUM,PetscObjectComm((PetscObject)pc));CHKERRQ(ierr); 843 *first -= osm->n; 844 } 845 if (ksp) { 846 /* Assume that local solves are now different; not necessarily 847 true, though! This flag is used only for PCView_GASM() */ 848 *ksp = osm->ksp; 849 osm->same_subdomain_solvers = PETSC_FALSE; 850 } 851 PetscFunctionReturn(0); 852 } /* PCGASMGetSubKSP_GASM() */ 853 EXTERN_C_END 854 855 856 #undef __FUNCT__ 857 #define __FUNCT__ "PCGASMSetSubdomains" 858 /*@C 859 PCGASMSetSubdomains - Sets the subdomains for this processor 860 for the additive Schwarz preconditioner. 861 862 Collective on PC 863 864 Input Parameters: 865 + pc - the preconditioner context 866 . n - the number of subdomains for this processor 867 . iis - the index sets that define this processor's local inner subdomains 868 (or NULL for PETSc to determine subdomains) 869 - ois- the index sets that define this processor's local outer subdomains 870 (or NULL to use the same as iis) 871 872 Notes: 873 The IS indices use the parallel, global numbering of the vector entries. 874 Inner subdomains are those where the correction is applied. 875 Outer subdomains are those where the residual necessary to obtain the 876 corrections is obtained (see PCGASMType for the use of inner/outer subdomains). 877 Both inner and outer subdomains can extend over several processors. 878 This processor's portion of a subdomain is known as a local subdomain. 879 880 By default the GASM preconditioner uses 1 (local) subdomain per processor. 881 Use PCGASMSetTotalSubdomains() to set the total number of subdomains across 882 all processors that PCGASM will create automatically, and to specify whether 883 they should be local or not. 884 885 886 Level: advanced 887 888 .keywords: PC, GASM, set, subdomains, additive Schwarz 889 890 .seealso: PCGASMSetTotalSubdomains(), PCGASMSetOverlap(), PCGASMGetSubKSP(), 891 PCGASMCreateSubdomains2D(), PCGASMGetSubdomains() 892 @*/ 893 PetscErrorCode PCGASMSetSubdomains(PC pc,PetscInt n,IS iis[],IS ois[]) 894 { 895 PetscErrorCode ierr; 896 897 PetscFunctionBegin; 898 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 899 ierr = PetscTryMethod(pc,"PCGASMSetSubdomains_C",(PC,PetscInt,IS[],IS[]),(pc,n,iis,ois));CHKERRQ(ierr); 900 PetscFunctionReturn(0); 901 } 902 903 #undef __FUNCT__ 904 #define __FUNCT__ "PCGASMSetTotalSubdomains" 905 /*@C 906 PCGASMSetTotalSubdomains - Sets the total number of subdomains to use in the generalized additive 907 Schwarz preconditioner. The number of subdomains is cumulative across all processors in pc's 908 communicator. Either all or no processors in the PC communicator must call this routine with 909 the same N. The subdomains will be created automatically during PCSetUp(). 910 911 Collective on PC 912 913 Input Parameters: 914 + pc - the preconditioner context 915 . N - the total number of subdomains cumulative across all processors 916 - create_local - whether the subdomains to be created are to be local 917 918 Options Database Key: 919 To set the total number of subdomains and let PCGASM autocreate them, rather than specify the index sets, use the following options: 920 + -pc_gasm_total_subdomains <n> - sets the total number of subdomains to be autocreated by PCGASM 921 - -pc_gasm_subdomains_create_local <true|false> - whether autocreated subdomains should be local or not (default is true) 922 923 By default the GASM preconditioner uses 1 subdomain per processor. 924 925 926 Use PCGASMSetSubdomains() to set subdomains explicitly or to set different numbers 927 of subdomains per processor. 928 929 Level: advanced 930 931 .keywords: PC, GASM, set, total, global, subdomains, additive Schwarz 932 933 .seealso: PCGASMSetSubdomains(), PCGASMSetOverlap(), PCGASMGetSubKSP(), 934 PCGASMCreateSubdomains2D() 935 @*/ 936 PetscErrorCode PCGASMSetTotalSubdomains(PC pc,PetscInt N, PetscBool create_local) 937 { 938 PetscErrorCode ierr; 939 940 PetscFunctionBegin; 941 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 942 ierr = PetscTryMethod(pc,"PCGASMSetTotalSubdomains_C",(PC,PetscInt,PetscBool),(pc,N,create_local));CHKERRQ(ierr); 943 PetscFunctionReturn(0); 944 } 945 946 #undef __FUNCT__ 947 #define __FUNCT__ "PCGASMSetOverlap" 948 /*@ 949 PCGASMSetOverlap - Sets the overlap between a pair of subdomains for the 950 additive Schwarz preconditioner. Either all or no processors in the 951 PC communicator must call this routine. 952 953 Logically Collective on PC 954 955 Input Parameters: 956 + pc - the preconditioner context 957 - ovl - the amount of overlap between subdomains (ovl >= 0, default value = 1) 958 959 Options Database Key: 960 . -pc_gasm_overlap <overlap> - Sets overlap 961 962 Notes: 963 By default the GASM preconditioner uses 1 subdomain per processor. To use 964 multiple subdomain per perocessor, see PCGASMSetTotalSubdomains() or 965 PCGASMSetSubdomains() (and the option -pc_gasm_total_subdomains <n>). 966 967 The overlap defaults to 1, so if one desires that no additional 968 overlap be computed beyond what may have been set with a call to 969 PCGASMSetTotalSubdomains() or PCGASMSetSubdomains(), then ovl 970 must be set to be 0. In particular, if one does not explicitly set 971 the subdomains in application code, then all overlap would be computed 972 internally by PETSc, and using an overlap of 0 would result in an GASM 973 variant that is equivalent to the block Jacobi preconditioner. 974 975 Note that one can define initial index sets with any overlap via 976 PCGASMSetSubdomains(); the routine PCGASMSetOverlap() merely allows 977 PETSc to extend that overlap further, if desired. 978 979 Level: intermediate 980 981 .keywords: PC, GASM, set, overlap 982 983 .seealso: PCGASMSetTotalSubdomains(), PCGASMSetSubdomains(), PCGASMGetSubKSP(), 984 PCGASMCreateSubdomains2D(), PCGASMGetSubdomains() 985 @*/ 986 PetscErrorCode PCGASMSetOverlap(PC pc,PetscInt ovl) 987 { 988 PetscErrorCode ierr; 989 990 PetscFunctionBegin; 991 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 992 PetscValidLogicalCollectiveInt(pc,ovl,2); 993 ierr = PetscTryMethod(pc,"PCGASMSetOverlap_C",(PC,PetscInt),(pc,ovl));CHKERRQ(ierr); 994 PetscFunctionReturn(0); 995 } 996 997 #undef __FUNCT__ 998 #define __FUNCT__ "PCGASMSetType" 999 /*@ 1000 PCGASMSetType - Sets the type of restriction and interpolation used 1001 for local problems in the additive Schwarz method. 1002 1003 Logically Collective on PC 1004 1005 Input Parameters: 1006 + pc - the preconditioner context 1007 - type - variant of GASM, one of 1008 .vb 1009 PC_GASM_BASIC - full interpolation and restriction 1010 PC_GASM_RESTRICT - full restriction, local processor interpolation 1011 PC_GASM_INTERPOLATE - full interpolation, local processor restriction 1012 PC_GASM_NONE - local processor restriction and interpolation 1013 .ve 1014 1015 Options Database Key: 1016 . -pc_gasm_type [basic,restrict,interpolate,none] - Sets GASM type 1017 1018 Level: intermediate 1019 1020 .keywords: PC, GASM, set, type 1021 1022 .seealso: PCGASMSetTotalSubdomains(), PCGASMSetSubdomains(), PCGASMGetSubKSP(), 1023 PCGASMCreateSubdomains2D() 1024 @*/ 1025 PetscErrorCode PCGASMSetType(PC pc,PCGASMType type) 1026 { 1027 PetscErrorCode ierr; 1028 1029 PetscFunctionBegin; 1030 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1031 PetscValidLogicalCollectiveEnum(pc,type,2); 1032 ierr = PetscTryMethod(pc,"PCGASMSetType_C",(PC,PCGASMType),(pc,type));CHKERRQ(ierr); 1033 PetscFunctionReturn(0); 1034 } 1035 1036 #undef __FUNCT__ 1037 #define __FUNCT__ "PCGASMSetSortIndices" 1038 /*@ 1039 PCGASMSetSortIndices - Determines whether subdomain indices are sorted. 1040 1041 Logically Collective on PC 1042 1043 Input Parameters: 1044 + pc - the preconditioner context 1045 - doSort - sort the subdomain indices 1046 1047 Level: intermediate 1048 1049 .keywords: PC, GASM, set, type 1050 1051 .seealso: PCGASMSetSubdomains(), PCGASMSetTotalSubdomains(), PCGASMGetSubKSP(), 1052 PCGASMCreateSubdomains2D() 1053 @*/ 1054 PetscErrorCode PCGASMSetSortIndices(PC pc,PetscBool doSort) 1055 { 1056 PetscErrorCode ierr; 1057 1058 PetscFunctionBegin; 1059 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1060 PetscValidLogicalCollectiveBool(pc,doSort,2); 1061 ierr = PetscTryMethod(pc,"PCGASMSetSortIndices_C",(PC,PetscBool),(pc,doSort));CHKERRQ(ierr); 1062 PetscFunctionReturn(0); 1063 } 1064 1065 #undef __FUNCT__ 1066 #define __FUNCT__ "PCGASMGetSubKSP" 1067 /*@C 1068 PCGASMGetSubKSP - Gets the local KSP contexts for all blocks on 1069 this processor. 1070 1071 Collective on PC iff first_local is requested 1072 1073 Input Parameter: 1074 . pc - the preconditioner context 1075 1076 Output Parameters: 1077 + n_local - the number of blocks on this processor or NULL 1078 . first_local - the global number of the first block on this processor or NULL, 1079 all processors must request or all must pass NULL 1080 - ksp - the array of KSP contexts 1081 1082 Note: 1083 After PCGASMGetSubKSP() the array of KSPes is not to be freed 1084 1085 Currently for some matrix implementations only 1 block per processor 1086 is supported. 1087 1088 You must call KSPSetUp() before calling PCGASMGetSubKSP(). 1089 1090 Level: advanced 1091 1092 .keywords: PC, GASM, additive Schwarz, get, sub, KSP, context 1093 1094 .seealso: PCGASMSetTotalSubdomains(), PCGASMSetSubdomains(), PCGASMSetOverlap(), 1095 PCGASMCreateSubdomains2D(), 1096 @*/ 1097 PetscErrorCode PCGASMGetSubKSP(PC pc,PetscInt *n_local,PetscInt *first_local,KSP *ksp[]) 1098 { 1099 PetscErrorCode ierr; 1100 1101 PetscFunctionBegin; 1102 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1103 ierr = PetscUseMethod(pc,"PCGASMGetSubKSP_C",(PC,PetscInt*,PetscInt*,KSP **),(pc,n_local,first_local,ksp));CHKERRQ(ierr); 1104 PetscFunctionReturn(0); 1105 } 1106 1107 /* -------------------------------------------------------------------------------------*/ 1108 /*MC 1109 PCGASM - Use the (restricted) additive Schwarz method, each block is (approximately) solved with 1110 its own KSP object. 1111 1112 Options Database Keys: 1113 + -pc_gasm_total_block_count <n> - Sets total number of local subdomains (known as blocks) to be distributed among processors 1114 . -pc_gasm_view_subdomains - activates the printing of subdomain indices in PCView(), -ksp_view or -snes_view 1115 . -pc_gasm_print_subdomains - activates the printing of subdomain indices in PCSetUp() 1116 . -pc_gasm_overlap <ovl> - Sets overlap by which to (automatically) extend local subdomains 1117 - -pc_gasm_type [basic,restrict,interpolate,none] - Sets GASM type 1118 1119 IMPORTANT: If you run with, for example, 3 blocks on 1 processor or 3 blocks on 3 processors you 1120 will get a different convergence rate due to the default option of -pc_gasm_type restrict. Use 1121 -pc_gasm_type basic to use the standard GASM. 1122 1123 Notes: Each processor can have one or more blocks, but a block cannot be shared by more 1124 than one processor. Defaults to one block per processor. 1125 1126 To set options on the solvers for each block append -sub_ to all the KSP, and PC 1127 options database keys. For example, -sub_pc_type ilu -sub_pc_factor_levels 1 -sub_ksp_type preonly 1128 1129 To set the options on the solvers separate for each block call PCGASMGetSubKSP() 1130 and set the options directly on the resulting KSP object (you can access its PC 1131 with KSPGetPC()) 1132 1133 1134 Level: beginner 1135 1136 Concepts: additive Schwarz method 1137 1138 References: 1139 An additive variant of the Schwarz alternating method for the case of many subregions 1140 M Dryja, OB Widlund - Courant Institute, New York University Technical report 1141 1142 Domain Decompositions: Parallel Multilevel Methods for Elliptic Partial Differential Equations, 1143 Barry Smith, Petter Bjorstad, and William Gropp, Cambridge University Press, ISBN 0-521-49589-X. 1144 1145 .seealso: PCCreate(), PCSetType(), PCType (for list of available types), PC, 1146 PCBJACOBI, PCGASMSetUseTrueLocal(), PCGASMGetSubKSP(), PCGASMSetSubdomains(), 1147 PCGASMSetTotalSubdomains(), PCSetModifySubmatrices(), PCGASMSetOverlap(), PCGASMSetType() 1148 1149 M*/ 1150 1151 EXTERN_C_BEGIN 1152 #undef __FUNCT__ 1153 #define __FUNCT__ "PCCreate_GASM" 1154 PetscErrorCode PCCreate_GASM(PC pc) 1155 { 1156 PetscErrorCode ierr; 1157 PC_GASM *osm; 1158 1159 PetscFunctionBegin; 1160 ierr = PetscNewLog(pc,PC_GASM,&osm);CHKERRQ(ierr); 1161 1162 osm->N = PETSC_DECIDE; 1163 osm->n = PETSC_DECIDE; 1164 osm->nmax = 0; 1165 osm->overlap = 1; 1166 osm->ksp = 0; 1167 osm->gorestriction = 0; 1168 osm->girestriction = 0; 1169 osm->gx = 0; 1170 osm->gy = 0; 1171 osm->x = 0; 1172 osm->y = 0; 1173 osm->ois = 0; 1174 osm->iis = 0; 1175 osm->pmat = 0; 1176 osm->type = PC_GASM_RESTRICT; 1177 osm->same_subdomain_solvers = PETSC_TRUE; 1178 osm->sort_indices = PETSC_TRUE; 1179 osm->use_dm_decomposition = PETSC_TRUE; 1180 1181 pc->data = (void*)osm; 1182 pc->ops->apply = PCApply_GASM; 1183 pc->ops->applytranspose = PCApplyTranspose_GASM; 1184 pc->ops->setup = PCSetUp_GASM; 1185 pc->ops->reset = PCReset_GASM; 1186 pc->ops->destroy = PCDestroy_GASM; 1187 pc->ops->setfromoptions = PCSetFromOptions_GASM; 1188 pc->ops->setuponblocks = PCSetUpOnBlocks_GASM; 1189 pc->ops->view = PCView_GASM; 1190 pc->ops->applyrichardson = 0; 1191 1192 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCGASMSetSubdomains_C","PCGASMSetSubdomains_GASM", 1193 PCGASMSetSubdomains_GASM);CHKERRQ(ierr); 1194 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCGASMSetTotalSubdomains_C","PCGASMSetTotalSubdomains_GASM", 1195 PCGASMSetTotalSubdomains_GASM);CHKERRQ(ierr); 1196 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCGASMSetOverlap_C","PCGASMSetOverlap_GASM", 1197 PCGASMSetOverlap_GASM);CHKERRQ(ierr); 1198 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCGASMSetType_C","PCGASMSetType_GASM", 1199 PCGASMSetType_GASM);CHKERRQ(ierr); 1200 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCGASMSetSortIndices_C","PCGASMSetSortIndices_GASM", 1201 PCGASMSetSortIndices_GASM);CHKERRQ(ierr); 1202 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCGASMGetSubKSP_C","PCGASMGetSubKSP_GASM", 1203 PCGASMGetSubKSP_GASM);CHKERRQ(ierr); 1204 PetscFunctionReturn(0); 1205 } 1206 EXTERN_C_END 1207 1208 1209 #undef __FUNCT__ 1210 #define __FUNCT__ "PCGASMCreateLocalSubdomains" 1211 /*@C 1212 PCGASMCreateLocalSubdomains - Creates n local index sets for the overlapping 1213 Schwarz preconditioner for a any problem based on its matrix. 1214 1215 Collective 1216 1217 Input Parameters: 1218 + A - The global matrix operator 1219 . overlap - amount of overlap in outer subdomains 1220 - n - the number of local subdomains 1221 1222 Output Parameters: 1223 + iis - the array of index sets defining the local inner subdomains (on which the correction is applied) 1224 - ois - the array of index sets defining the local outer subdomains (on which the residual is computed) 1225 1226 Level: advanced 1227 1228 Note: this generates n nonoverlapping local inner subdomains on PETSC_COMM_SELF; 1229 PCGASM will generate the overlap from these if you use them in PCGASMSetSubdomains() and set a 1230 nonzero overlap with PCGASMSetOverlap() 1231 1232 In the Fortran version you must provide the array outis[] already allocated of length n. 1233 1234 .keywords: PC, GASM, additive Schwarz, create, subdomains, unstructured grid 1235 1236 .seealso: PCGASMSetSubdomains(), PCGASMDestroySubdomains() 1237 @*/ 1238 PetscErrorCode PCGASMCreateLocalSubdomains(Mat A, PetscInt overlap, PetscInt n, IS *iis[], IS *ois[]) 1239 { 1240 MatPartitioning mpart; 1241 const char *prefix; 1242 PetscErrorCode (*f)(Mat,MatReuse,Mat*); 1243 PetscMPIInt size; 1244 PetscInt i,j,rstart,rend,bs; 1245 PetscBool isbaij = PETSC_FALSE,foundpart = PETSC_FALSE; 1246 Mat Ad = NULL, adj; 1247 IS ispart,isnumb,*is; 1248 PetscErrorCode ierr; 1249 1250 PetscFunctionBegin; 1251 PetscValidHeaderSpecific(A,MAT_CLASSID,1); 1252 PetscValidPointer(iis,4); 1253 if (n < 1) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"number of local blocks must be > 0, n = %D",n); 1254 1255 /* Get prefix, row distribution, and block size */ 1256 ierr = MatGetOptionsPrefix(A,&prefix);CHKERRQ(ierr); 1257 ierr = MatGetOwnershipRange(A,&rstart,&rend);CHKERRQ(ierr); 1258 ierr = MatGetBlockSize(A,&bs);CHKERRQ(ierr); 1259 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); 1260 1261 /* Get diagonal block from matrix if possible */ 1262 ierr = MPI_Comm_size(PetscObjectComm((PetscObject)A),&size);CHKERRQ(ierr); 1263 ierr = PetscObjectQueryFunction((PetscObject)A,"MatGetDiagonalBlock_C",(void (**)(void))&f);CHKERRQ(ierr); 1264 if (f) { 1265 ierr = MatGetDiagonalBlock(A,&Ad);CHKERRQ(ierr); 1266 } else if (size == 1) { 1267 Ad = A; 1268 } 1269 if (Ad) { 1270 ierr = PetscObjectTypeCompare((PetscObject)Ad,MATSEQBAIJ,&isbaij);CHKERRQ(ierr); 1271 if (!isbaij) {ierr = PetscObjectTypeCompare((PetscObject)Ad,MATSEQSBAIJ,&isbaij);CHKERRQ(ierr);} 1272 } 1273 if (Ad && n > 1) { 1274 PetscBool match,done; 1275 /* Try to setup a good matrix partitioning if available */ 1276 ierr = MatPartitioningCreate(PETSC_COMM_SELF,&mpart);CHKERRQ(ierr); 1277 ierr = PetscObjectSetOptionsPrefix((PetscObject)mpart,prefix);CHKERRQ(ierr); 1278 ierr = MatPartitioningSetFromOptions(mpart);CHKERRQ(ierr); 1279 ierr = PetscObjectTypeCompare((PetscObject)mpart,MATPARTITIONINGCURRENT,&match);CHKERRQ(ierr); 1280 if (!match) { 1281 ierr = PetscObjectTypeCompare((PetscObject)mpart,MATPARTITIONINGSQUARE,&match);CHKERRQ(ierr); 1282 } 1283 if (!match) { /* assume a "good" partitioner is available */ 1284 PetscInt na; 1285 const PetscInt *ia,*ja; 1286 ierr = MatGetRowIJ(Ad,0,PETSC_TRUE,isbaij,&na,&ia,&ja,&done);CHKERRQ(ierr); 1287 if (done) { 1288 /* Build adjacency matrix by hand. Unfortunately a call to 1289 MatConvert(Ad,MATMPIADJ,MAT_INITIAL_MATRIX,&adj) will 1290 remove the block-aij structure and we cannot expect 1291 MatPartitioning to split vertices as we need */ 1292 PetscInt i,j,len,nnz,cnt,*iia=0,*jja=0; 1293 const PetscInt *row; 1294 nnz = 0; 1295 for (i=0; i<na; i++) { /* count number of nonzeros */ 1296 len = ia[i+1] - ia[i]; 1297 row = ja + ia[i]; 1298 for (j=0; j<len; j++) { 1299 if (row[j] == i) { /* don't count diagonal */ 1300 len--; break; 1301 } 1302 } 1303 nnz += len; 1304 } 1305 ierr = PetscMalloc((na+1)*sizeof(PetscInt),&iia);CHKERRQ(ierr); 1306 ierr = PetscMalloc((nnz)*sizeof(PetscInt),&jja);CHKERRQ(ierr); 1307 nnz = 0; 1308 iia[0] = 0; 1309 for (i=0; i<na; i++) { /* fill adjacency */ 1310 cnt = 0; 1311 len = ia[i+1] - ia[i]; 1312 row = ja + ia[i]; 1313 for (j=0; j<len; j++) { 1314 if (row[j] != i) jja[nnz+cnt++] = row[j]; /* if not diagonal */ 1315 } 1316 nnz += cnt; 1317 iia[i+1] = nnz; 1318 } 1319 /* Partitioning of the adjacency matrix */ 1320 ierr = MatCreateMPIAdj(PETSC_COMM_SELF,na,na,iia,jja,NULL,&adj);CHKERRQ(ierr); 1321 ierr = MatPartitioningSetAdjacency(mpart,adj);CHKERRQ(ierr); 1322 ierr = MatPartitioningSetNParts(mpart,n);CHKERRQ(ierr); 1323 ierr = MatPartitioningApply(mpart,&ispart);CHKERRQ(ierr); 1324 ierr = ISPartitioningToNumbering(ispart,&isnumb);CHKERRQ(ierr); 1325 ierr = MatDestroy(&adj);CHKERRQ(ierr); 1326 foundpart = PETSC_TRUE; 1327 } 1328 ierr = MatRestoreRowIJ(Ad,0,PETSC_TRUE,isbaij,&na,&ia,&ja,&done);CHKERRQ(ierr); 1329 } 1330 ierr = MatPartitioningDestroy(&mpart);CHKERRQ(ierr); 1331 } 1332 ierr = PetscMalloc(n*sizeof(IS),&is);CHKERRQ(ierr); 1333 if (!foundpart) { 1334 1335 /* Partitioning by contiguous chunks of rows */ 1336 1337 PetscInt mbs = (rend-rstart)/bs; 1338 PetscInt start = rstart; 1339 for (i=0; i<n; i++) { 1340 PetscInt count = (mbs/n + ((mbs % n) > i)) * bs; 1341 ierr = ISCreateStride(PETSC_COMM_SELF,count,start,1,&is[i]);CHKERRQ(ierr); 1342 start += count; 1343 } 1344 1345 } else { 1346 1347 /* Partitioning by adjacency of diagonal block */ 1348 1349 const PetscInt *numbering; 1350 PetscInt *count,nidx,*indices,*newidx,start=0; 1351 /* Get node count in each partition */ 1352 ierr = PetscMalloc(n*sizeof(PetscInt),&count);CHKERRQ(ierr); 1353 ierr = ISPartitioningCount(ispart,n,count);CHKERRQ(ierr); 1354 if (isbaij && bs > 1) { /* adjust for the block-aij case */ 1355 for (i=0; i<n; i++) count[i] *= bs; 1356 } 1357 /* Build indices from node numbering */ 1358 ierr = ISGetLocalSize(isnumb,&nidx);CHKERRQ(ierr); 1359 ierr = PetscMalloc(nidx*sizeof(PetscInt),&indices);CHKERRQ(ierr); 1360 for (i=0; i<nidx; i++) indices[i] = i; /* needs to be initialized */ 1361 ierr = ISGetIndices(isnumb,&numbering);CHKERRQ(ierr); 1362 ierr = PetscSortIntWithPermutation(nidx,numbering,indices);CHKERRQ(ierr); 1363 ierr = ISRestoreIndices(isnumb,&numbering);CHKERRQ(ierr); 1364 if (isbaij && bs > 1) { /* adjust for the block-aij case */ 1365 ierr = PetscMalloc(nidx*bs*sizeof(PetscInt),&newidx);CHKERRQ(ierr); 1366 for (i=0; i<nidx; i++) { 1367 for (j=0; j<bs; j++) newidx[i*bs+j] = indices[i]*bs + j; 1368 } 1369 ierr = PetscFree(indices);CHKERRQ(ierr); 1370 nidx *= bs; 1371 indices = newidx; 1372 } 1373 /* Shift to get global indices */ 1374 for (i=0; i<nidx; i++) indices[i] += rstart; 1375 1376 /* Build the index sets for each block */ 1377 for (i=0; i<n; i++) { 1378 ierr = ISCreateGeneral(PETSC_COMM_SELF,count[i],&indices[start],PETSC_COPY_VALUES,&is[i]);CHKERRQ(ierr); 1379 ierr = ISSort(is[i]);CHKERRQ(ierr); 1380 start += count[i]; 1381 } 1382 1383 ierr = PetscFree(count); 1384 ierr = PetscFree(indices); 1385 ierr = ISDestroy(&isnumb);CHKERRQ(ierr); 1386 ierr = ISDestroy(&ispart);CHKERRQ(ierr); 1387 } 1388 *iis = is; 1389 if (!ois) PetscFunctionReturn(0); 1390 /* 1391 Initially make outer subdomains the same as inner subdomains. If nonzero additional overlap 1392 has been requested, copy the inner subdomains over so they can be modified. 1393 */ 1394 ierr = PetscMalloc(n*sizeof(IS),ois);CHKERRQ(ierr); 1395 for (i=0; i<n; ++i) { 1396 if (overlap > 0) { /* With positive overlap, (*iis)[i] will be modified */ 1397 ierr = ISDuplicate((*iis)[i],(*ois)+i);CHKERRQ(ierr); 1398 ierr = ISCopy((*iis)[i],(*ois)[i]);CHKERRQ(ierr); 1399 } else { 1400 ierr = PetscObjectReference((PetscObject)(*iis)[i]);CHKERRQ(ierr); 1401 (*ois)[i] = (*iis)[i]; 1402 } 1403 } 1404 if (overlap > 0) { 1405 /* Extend the "overlapping" regions by a number of steps */ 1406 ierr = MatIncreaseOverlap(A,n,*ois,overlap);CHKERRQ(ierr); 1407 } 1408 PetscFunctionReturn(0); 1409 } 1410 1411 #undef __FUNCT__ 1412 #define __FUNCT__ "PCGASMDestroySubdomains" 1413 /*@C 1414 PCGASMDestroySubdomains - Destroys the index sets created with 1415 PCGASMCreateLocalSubdomains() or PCGASMCreateSubdomains2D. Should be 1416 called after setting subdomains with PCGASMSetSubdomains(). 1417 1418 Collective 1419 1420 Input Parameters: 1421 + n - the number of index sets 1422 . iis - the array of inner subdomains, 1423 - ois - the array of outer subdomains, can be NULL 1424 1425 Level: intermediate 1426 1427 Notes: this is merely a convenience subroutine that walks each list, 1428 destroys each IS on the list, and then frees the list. 1429 1430 .keywords: PC, GASM, additive Schwarz, create, subdomains, unstructured grid 1431 1432 .seealso: PCGASMCreateLocalSubdomains(), PCGASMSetSubdomains() 1433 @*/ 1434 PetscErrorCode PCGASMDestroySubdomains(PetscInt n, IS iis[], IS ois[]) 1435 { 1436 PetscInt i; 1437 PetscErrorCode ierr; 1438 1439 PetscFunctionBegin; 1440 if (n <= 0) PetscFunctionReturn(0); 1441 if (iis) { 1442 PetscValidPointer(iis,2); 1443 for (i=0; i<n; i++) { 1444 ierr = ISDestroy(&iis[i]);CHKERRQ(ierr); 1445 } 1446 ierr = PetscFree(iis);CHKERRQ(ierr); 1447 } 1448 if (ois) { 1449 for (i=0; i<n; i++) { 1450 ierr = ISDestroy(&ois[i]);CHKERRQ(ierr); 1451 } 1452 ierr = PetscFree(ois);CHKERRQ(ierr); 1453 } 1454 PetscFunctionReturn(0); 1455 } 1456 1457 1458 #define PCGASMLocalSubdomainBounds2D(M,N,xleft,ylow,xright,yhigh,first,last,xleft_loc,ylow_loc,xright_loc,yhigh_loc,n) \ 1459 { \ 1460 PetscInt first_row = first/M, last_row = last/M+1; \ 1461 /* \ 1462 Compute ylow_loc and yhigh_loc so that (ylow_loc,xleft) and (yhigh_loc,xright) are the corners \ 1463 of the bounding box of the intersection of the subdomain with the local ownership range (local \ 1464 subdomain). \ 1465 Also compute xleft_loc and xright_loc as the lower and upper bounds on the first and last rows \ 1466 of the intersection. \ 1467 */ \ 1468 /* ylow_loc is the grid row containing the first element of the local sumbdomain */ \ 1469 *ylow_loc = PetscMax(first_row,ylow); \ 1470 /* xleft_loc is the offset of first element of the local subdomain within its grid row (might actually be outside the local subdomain) */ \ 1471 *xleft_loc = *ylow_loc==first_row ? PetscMax(first%M,xleft) : xleft; \ 1472 /* yhigh_loc is the grid row above the last local subdomain element */ \ 1473 *yhigh_loc = PetscMin(last_row,yhigh); \ 1474 /* xright is the offset of the end of the local subdomain within its grid row (might actually be outside the local subdomain) */ \ 1475 *xright_loc = *yhigh_loc==last_row ? PetscMin(xright,last%M) : xright; \ 1476 /* Now compute the size of the local subdomain n. */ \ 1477 *n = 0; \ 1478 if (*ylow_loc < *yhigh_loc) { \ 1479 PetscInt width = xright-xleft; \ 1480 *n += width*(*yhigh_loc-*ylow_loc-1); \ 1481 *n += PetscMin(PetscMax(*xright_loc-xleft,0),width); \ 1482 *n -= PetscMin(PetscMax(*xleft_loc-xleft,0), width); \ 1483 } \ 1484 } 1485 1486 1487 1488 #undef __FUNCT__ 1489 #define __FUNCT__ "PCGASMCreateSubdomains2D" 1490 /*@ 1491 PCGASMCreateSubdomains2D - Creates the index sets for the overlapping Schwarz 1492 preconditioner for a two-dimensional problem on a regular grid. 1493 1494 Collective 1495 1496 Input Parameters: 1497 + M, N - the global number of grid points in the x and y directions 1498 . Mdomains, Ndomains - the global number of subdomains in the x and y directions 1499 . dof - degrees of freedom per node 1500 - overlap - overlap in mesh lines 1501 1502 Output Parameters: 1503 + Nsub - the number of local subdomains created 1504 . iis - array of index sets defining inner (nonoverlapping) subdomains 1505 - ois - array of index sets defining outer (overlapping, if overlap > 0) subdomains 1506 1507 1508 Level: advanced 1509 1510 .keywords: PC, GASM, additive Schwarz, create, subdomains, 2D, regular grid 1511 1512 .seealso: PCGASMSetTotalSubdomains(), PCGASMSetSubdomains(), PCGASMGetSubKSP(), 1513 PCGASMSetOverlap() 1514 @*/ 1515 PetscErrorCode PCGASMCreateSubdomains2D(PC pc, PetscInt M,PetscInt N,PetscInt Mdomains,PetscInt Ndomains,PetscInt dof,PetscInt overlap, PetscInt *nsub,IS **iis,IS **ois) 1516 { 1517 PetscErrorCode ierr; 1518 PetscMPIInt size, rank; 1519 PetscInt i, j; 1520 PetscInt maxheight, maxwidth; 1521 PetscInt xstart, xleft, xright, xleft_loc, xright_loc; 1522 PetscInt ystart, ylow, yhigh, ylow_loc, yhigh_loc; 1523 PetscInt x[2][2], y[2][2], n[2]; 1524 PetscInt first, last; 1525 PetscInt nidx, *idx; 1526 PetscInt ii,jj,s,q,d; 1527 PetscInt k,kk; 1528 PetscMPIInt color; 1529 MPI_Comm comm, subcomm; 1530 IS **xis = 0, **is = ois, **is_local = iis; 1531 1532 PetscFunctionBegin; 1533 ierr = PetscObjectGetComm((PetscObject)pc, &comm);CHKERRQ(ierr); 1534 ierr = MPI_Comm_size(comm, &size);CHKERRQ(ierr); 1535 ierr = MPI_Comm_rank(comm, &rank);CHKERRQ(ierr); 1536 ierr = MatGetOwnershipRange(pc->pmat, &first, &last);CHKERRQ(ierr); 1537 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) " 1538 "does not respect the number of degrees of freedom per grid point %D", first, last, dof); 1539 1540 /* Determine the number of domains with nonzero intersections with the local ownership range. */ 1541 s = 0; 1542 ystart = 0; 1543 for (j=0; j<Ndomains; ++j) { 1544 maxheight = N/Ndomains + ((N % Ndomains) > j); /* Maximal height of subdomain */ 1545 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); 1546 /* Vertical domain limits with an overlap. */ 1547 ylow = PetscMax(ystart - overlap,0); 1548 yhigh = PetscMin(ystart + maxheight + overlap,N); 1549 xstart = 0; 1550 for (i=0; i<Mdomains; ++i) { 1551 maxwidth = M/Mdomains + ((M % Mdomains) > i); /* Maximal width of subdomain */ 1552 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); 1553 /* Horizontal domain limits with an overlap. */ 1554 xleft = PetscMax(xstart - overlap,0); 1555 xright = PetscMin(xstart + maxwidth + overlap,M); 1556 /* 1557 Determine whether this subdomain intersects this processor's ownership range of pc->pmat. 1558 */ 1559 PCGASMLocalSubdomainBounds2D(M,N,xleft,ylow,xright,yhigh,first,last,(&xleft_loc),(&ylow_loc),(&xright_loc),(&yhigh_loc),(&nidx)); 1560 if (nidx) ++s; 1561 xstart += maxwidth; 1562 } /* for (i = 0; i < Mdomains; ++i) */ 1563 ystart += maxheight; 1564 } /* for (j = 0; j < Ndomains; ++j) */ 1565 1566 /* Now we can allocate the necessary number of ISs. */ 1567 *nsub = s; 1568 ierr = PetscMalloc((*nsub)*sizeof(IS*),is);CHKERRQ(ierr); 1569 ierr = PetscMalloc((*nsub)*sizeof(IS*),is_local);CHKERRQ(ierr); 1570 s = 0; 1571 ystart = 0; 1572 for (j=0; j<Ndomains; ++j) { 1573 maxheight = N/Ndomains + ((N % Ndomains) > j); /* Maximal height of subdomain */ 1574 if (maxheight < 2) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Too many %D subdomains in the vertical directon for mesh height %D", Ndomains, N); 1575 /* Vertical domain limits with an overlap. */ 1576 y[0][0] = PetscMax(ystart - overlap,0); 1577 y[0][1] = PetscMin(ystart + maxheight + overlap,N); 1578 /* Vertical domain limits without an overlap. */ 1579 y[1][0] = ystart; 1580 y[1][1] = PetscMin(ystart + maxheight,N); 1581 xstart = 0; 1582 for (i=0; i<Mdomains; ++i) { 1583 maxwidth = M/Mdomains + ((M % Mdomains) > i); /* Maximal width of subdomain */ 1584 if (maxwidth < 2) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Too many %D subdomains in the horizontal direction for mesh width %D", Mdomains, M); 1585 /* Horizontal domain limits with an overlap. */ 1586 x[0][0] = PetscMax(xstart - overlap,0); 1587 x[0][1] = PetscMin(xstart + maxwidth + overlap,M); 1588 /* Horizontal domain limits without an overlap. */ 1589 x[1][0] = xstart; 1590 x[1][1] = PetscMin(xstart+maxwidth,M); 1591 /* 1592 Determine whether this domain intersects this processor's ownership range of pc->pmat. 1593 Do this twice: first for the domains with overlaps, and once without. 1594 During the first pass create the subcommunicators, and use them on the second pass as well. 1595 */ 1596 for (q = 0; q < 2; ++q) { 1597 /* 1598 domain limits, (xleft, xright) and (ylow, yheigh) are adjusted 1599 according to whether the domain with an overlap or without is considered. 1600 */ 1601 xleft = x[q][0]; xright = x[q][1]; 1602 ylow = y[q][0]; yhigh = y[q][1]; 1603 PCGASMLocalSubdomainBounds2D(M,N,xleft,ylow,xright,yhigh,first,last,(&xleft_loc),(&ylow_loc),(&xright_loc),(&yhigh_loc),(&nidx)); 1604 nidx *= dof; 1605 n[q] = nidx; 1606 /* 1607 Based on the counted number of indices in the local domain *with an overlap*, 1608 construct a subcommunicator of all the processors supporting this domain. 1609 Observe that a domain with an overlap might have nontrivial local support, 1610 while the domain without an overlap might not. Hence, the decision to participate 1611 in the subcommunicator must be based on the domain with an overlap. 1612 */ 1613 if (q == 0) { 1614 if (nidx) color = 1; 1615 else color = MPI_UNDEFINED; 1616 1617 ierr = MPI_Comm_split(comm, color, rank, &subcomm);CHKERRQ(ierr); 1618 } 1619 /* 1620 Proceed only if the number of local indices *with an overlap* is nonzero. 1621 */ 1622 if (n[0]) { 1623 if (q == 0) xis = is; 1624 if (q == 1) { 1625 /* 1626 The IS for the no-overlap subdomain shares a communicator with the overlapping domain. 1627 Moreover, if the overlap is zero, the two ISs are identical. 1628 */ 1629 if (overlap == 0) { 1630 (*is_local)[s] = (*is)[s]; 1631 ierr = PetscObjectReference((PetscObject)(*is)[s]);CHKERRQ(ierr); 1632 continue; 1633 } else { 1634 xis = is_local; 1635 subcomm = ((PetscObject)(*is)[s])->comm; 1636 } 1637 } /* if (q == 1) */ 1638 idx = NULL; 1639 ierr = PetscMalloc(nidx*sizeof(PetscInt),&idx);CHKERRQ(ierr); 1640 if (nidx) { 1641 k = 0; 1642 for (jj=ylow_loc; jj<yhigh_loc; ++jj) { 1643 PetscInt x0 = (jj==ylow_loc) ? xleft_loc : xleft; 1644 PetscInt x1 = (jj==yhigh_loc-1) ? xright_loc : xright; 1645 kk = dof*(M*jj + x0); 1646 for (ii=x0; ii<x1; ++ii) { 1647 for (d = 0; d < dof; ++d) { 1648 idx[k++] = kk++; 1649 } 1650 } 1651 } 1652 } 1653 ierr = ISCreateGeneral(subcomm,nidx,idx,PETSC_OWN_POINTER,(*xis)+s);CHKERRQ(ierr); 1654 }/* if (n[0]) */ 1655 }/* for (q = 0; q < 2; ++q) */ 1656 if (n[0]) ++s; 1657 xstart += maxwidth; 1658 } /* for (i = 0; i < Mdomains; ++i) */ 1659 ystart += maxheight; 1660 } /* for (j = 0; j < Ndomains; ++j) */ 1661 PetscFunctionReturn(0); 1662 } 1663 1664 #undef __FUNCT__ 1665 #define __FUNCT__ "PCGASMGetSubdomains" 1666 /*@C 1667 PCGASMGetSubdomains - Gets the subdomains supported on this processor 1668 for the additive Schwarz preconditioner. 1669 1670 Not Collective 1671 1672 Input Parameter: 1673 . pc - the preconditioner context 1674 1675 Output Parameters: 1676 + n - the number of subdomains for this processor (default value = 1) 1677 . iis - the index sets that define the inner subdomains (without overlap) supported on this processor (can be NULL) 1678 - ois - the index sets that define the outer subdomains (with overlap) supported on this processor (can be NULL) 1679 1680 1681 Notes: 1682 The user is responsible for destroying the ISs and freeing the returned arrays. 1683 The IS numbering is in the parallel, global numbering of the vector. 1684 1685 Level: advanced 1686 1687 .keywords: PC, GASM, get, subdomains, additive Schwarz 1688 1689 .seealso: PCGASMSetTotalSubdomains(), PCGASMSetOverlap(), PCGASMGetSubKSP(), 1690 PCGASMCreateSubdomains2D(), PCGASMSetSubdomains(), PCGASMGetSubmatrices() 1691 @*/ 1692 PetscErrorCode PCGASMGetSubdomains(PC pc,PetscInt *n,IS *iis[],IS *ois[]) 1693 { 1694 PC_GASM *osm; 1695 PetscErrorCode ierr; 1696 PetscBool match; 1697 PetscInt i; 1698 1699 PetscFunctionBegin; 1700 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1701 ierr = PetscObjectTypeCompare((PetscObject)pc,PCGASM,&match);CHKERRQ(ierr); 1702 if (!match) SETERRQ2(PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_WRONG, "Incorrect object type: expected %s, got %s instead", PCGASM, ((PetscObject)pc)->type_name); 1703 osm = (PC_GASM*)pc->data; 1704 if (n) *n = osm->n; 1705 if (iis) { 1706 ierr = PetscMalloc(osm->n*sizeof(IS), iis);CHKERRQ(ierr); 1707 } 1708 if (ois) { 1709 ierr = PetscMalloc(osm->n*sizeof(IS), ois);CHKERRQ(ierr); 1710 } 1711 if (iis || ois) { 1712 for (i = 0; i < osm->n; ++i) { 1713 if (iis) (*iis)[i] = osm->iis[i]; 1714 if (ois) (*ois)[i] = osm->ois[i]; 1715 } 1716 } 1717 PetscFunctionReturn(0); 1718 } 1719 1720 #undef __FUNCT__ 1721 #define __FUNCT__ "PCGASMGetSubmatrices" 1722 /*@C 1723 PCGASMGetSubmatrices - Gets the local submatrices (for this processor 1724 only) for the additive Schwarz preconditioner. 1725 1726 Not Collective 1727 1728 Input Parameter: 1729 . pc - the preconditioner context 1730 1731 Output Parameters: 1732 + n - the number of matrices for this processor (default value = 1) 1733 - mat - the matrices 1734 1735 Notes: matrices returned by this routine have the same communicators as the index sets (IS) 1736 used to define subdomains in PCGASMSetSubdomains(), or PETSC_COMM_SELF, if the 1737 subdomains were defined using PCGASMSetTotalSubdomains(). 1738 Level: advanced 1739 1740 .keywords: PC, GASM, set, local, subdomains, additive Schwarz, block Jacobi 1741 1742 .seealso: PCGASMSetTotalSubdomain(), PCGASMSetOverlap(), PCGASMGetSubKSP(), 1743 PCGASMCreateSubdomains2D(), PCGASMSetSubdomains(), PCGASMGetSubdomains() 1744 @*/ 1745 PetscErrorCode PCGASMGetSubmatrices(PC pc,PetscInt *n,Mat *mat[]) 1746 { 1747 PC_GASM *osm; 1748 PetscErrorCode ierr; 1749 PetscBool match; 1750 1751 PetscFunctionBegin; 1752 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1753 PetscValidIntPointer(n,2); 1754 if (mat) PetscValidPointer(mat,3); 1755 if (!pc->setupcalled) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONGSTATE,"Must call after KSPSetUP() or PCSetUp()."); 1756 ierr = PetscObjectTypeCompare((PetscObject)pc,PCGASM,&match);CHKERRQ(ierr); 1757 if (!match) SETERRQ2(PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_WRONG, "Expected %s, got %s instead", PCGASM, ((PetscObject)pc)->type_name); 1758 osm = (PC_GASM*)pc->data; 1759 if (n) *n = osm->n; 1760 if (mat) *mat = osm->pmat; 1761 PetscFunctionReturn(0); 1762 } 1763