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