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