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