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