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