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