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