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