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 + 1. - 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 - 2. - Barry Smith, Petter Bjorstad, and William Gropp, Domain Decompositions: Parallel Multilevel Methods for Elliptic Partial Differential Equations, 1297 Cambridge University Press. 1298 1299 .seealso: PCCreate(), PCSetType(), PCType (for list of available types), PC, 1300 PCBJACOBI, PCGASMGetSubKSP(), PCGASMSetSubdomains(), 1301 PCSetModifySubmatrices(), PCGASMSetOverlap(), PCGASMSetType() 1302 1303 M*/ 1304 1305 #undef __FUNCT__ 1306 #define __FUNCT__ "PCCreate_GASM" 1307 PETSC_EXTERN PetscErrorCode PCCreate_GASM(PC pc) 1308 { 1309 PetscErrorCode ierr; 1310 PC_GASM *osm; 1311 1312 PetscFunctionBegin; 1313 ierr = PetscNewLog(pc,&osm);CHKERRQ(ierr); 1314 1315 osm->N = PETSC_DETERMINE; 1316 osm->n = PETSC_DECIDE; 1317 osm->nmax = PETSC_DETERMINE; 1318 osm->overlap = 0; 1319 osm->ksp = 0; 1320 osm->gorestriction = 0; 1321 osm->girestriction = 0; 1322 osm->pctoouter = 0; 1323 osm->gx = 0; 1324 osm->gy = 0; 1325 osm->x = 0; 1326 osm->y = 0; 1327 osm->pcx = 0; 1328 osm->pcy = 0; 1329 osm->permutationIS = 0; 1330 osm->permutationP = 0; 1331 osm->pcmat = 0; 1332 osm->ois = 0; 1333 osm->iis = 0; 1334 osm->pmat = 0; 1335 osm->type = PC_GASM_RESTRICT; 1336 osm->same_subdomain_solvers = PETSC_TRUE; 1337 osm->sort_indices = PETSC_TRUE; 1338 osm->dm_subdomains = PETSC_FALSE; 1339 osm->hierarchicalpartitioning = PETSC_FALSE; 1340 1341 pc->data = (void*)osm; 1342 pc->ops->apply = PCApply_GASM; 1343 pc->ops->applytranspose = PCApplyTranspose_GASM; 1344 pc->ops->setup = PCSetUp_GASM; 1345 pc->ops->reset = PCReset_GASM; 1346 pc->ops->destroy = PCDestroy_GASM; 1347 pc->ops->setfromoptions = PCSetFromOptions_GASM; 1348 pc->ops->setuponblocks = PCSetUpOnBlocks_GASM; 1349 pc->ops->view = PCView_GASM; 1350 pc->ops->applyrichardson = 0; 1351 1352 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCGASMSetSubdomains_C",PCGASMSetSubdomains_GASM);CHKERRQ(ierr); 1353 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCGASMSetOverlap_C",PCGASMSetOverlap_GASM);CHKERRQ(ierr); 1354 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCGASMSetType_C",PCGASMSetType_GASM);CHKERRQ(ierr); 1355 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCGASMSetSortIndices_C",PCGASMSetSortIndices_GASM);CHKERRQ(ierr); 1356 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCGASMGetSubKSP_C",PCGASMGetSubKSP_GASM);CHKERRQ(ierr); 1357 PetscFunctionReturn(0); 1358 } 1359 1360 1361 #undef __FUNCT__ 1362 #define __FUNCT__ "PCGASMCreateLocalSubdomains" 1363 PetscErrorCode PCGASMCreateLocalSubdomains(Mat A, PetscInt nloc, IS *iis[]) 1364 { 1365 MatPartitioning mpart; 1366 const char *prefix; 1367 PetscErrorCode (*f)(Mat,MatReuse,Mat*); 1368 PetscMPIInt size; 1369 PetscInt i,j,rstart,rend,bs; 1370 PetscBool isbaij = PETSC_FALSE,foundpart = PETSC_FALSE; 1371 Mat Ad = NULL, adj; 1372 IS ispart,isnumb,*is; 1373 PetscErrorCode ierr; 1374 1375 PetscFunctionBegin; 1376 if (nloc < 1) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"number of local subdomains must > 0, got nloc = %D",nloc); 1377 1378 /* Get prefix, row distribution, and block size */ 1379 ierr = MatGetOptionsPrefix(A,&prefix);CHKERRQ(ierr); 1380 ierr = MatGetOwnershipRange(A,&rstart,&rend);CHKERRQ(ierr); 1381 ierr = MatGetBlockSize(A,&bs);CHKERRQ(ierr); 1382 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); 1383 1384 /* Get diagonal block from matrix if possible */ 1385 ierr = MPI_Comm_size(PetscObjectComm((PetscObject)A),&size);CHKERRQ(ierr); 1386 ierr = PetscObjectQueryFunction((PetscObject)A,"MatGetDiagonalBlock_C",&f);CHKERRQ(ierr); 1387 if (f) { 1388 ierr = MatGetDiagonalBlock(A,&Ad);CHKERRQ(ierr); 1389 } else if (size == 1) { 1390 Ad = A; 1391 } 1392 if (Ad) { 1393 ierr = PetscObjectTypeCompare((PetscObject)Ad,MATSEQBAIJ,&isbaij);CHKERRQ(ierr); 1394 if (!isbaij) {ierr = PetscObjectTypeCompare((PetscObject)Ad,MATSEQSBAIJ,&isbaij);CHKERRQ(ierr);} 1395 } 1396 if (Ad && nloc > 1) { 1397 PetscBool match,done; 1398 /* Try to setup a good matrix partitioning if available */ 1399 ierr = MatPartitioningCreate(PETSC_COMM_SELF,&mpart);CHKERRQ(ierr); 1400 ierr = PetscObjectSetOptionsPrefix((PetscObject)mpart,prefix);CHKERRQ(ierr); 1401 ierr = MatPartitioningSetFromOptions(mpart);CHKERRQ(ierr); 1402 ierr = PetscObjectTypeCompare((PetscObject)mpart,MATPARTITIONINGCURRENT,&match);CHKERRQ(ierr); 1403 if (!match) { 1404 ierr = PetscObjectTypeCompare((PetscObject)mpart,MATPARTITIONINGSQUARE,&match);CHKERRQ(ierr); 1405 } 1406 if (!match) { /* assume a "good" partitioner is available */ 1407 PetscInt na; 1408 const PetscInt *ia,*ja; 1409 ierr = MatGetRowIJ(Ad,0,PETSC_TRUE,isbaij,&na,&ia,&ja,&done);CHKERRQ(ierr); 1410 if (done) { 1411 /* Build adjacency matrix by hand. Unfortunately a call to 1412 MatConvert(Ad,MATMPIADJ,MAT_INITIAL_MATRIX,&adj) will 1413 remove the block-aij structure and we cannot expect 1414 MatPartitioning to split vertices as we need */ 1415 PetscInt i,j,len,nnz,cnt,*iia=0,*jja=0; 1416 const PetscInt *row; 1417 nnz = 0; 1418 for (i=0; i<na; i++) { /* count number of nonzeros */ 1419 len = ia[i+1] - ia[i]; 1420 row = ja + ia[i]; 1421 for (j=0; j<len; j++) { 1422 if (row[j] == i) { /* don't count diagonal */ 1423 len--; break; 1424 } 1425 } 1426 nnz += len; 1427 } 1428 ierr = PetscMalloc1(na+1,&iia);CHKERRQ(ierr); 1429 ierr = PetscMalloc1(nnz,&jja);CHKERRQ(ierr); 1430 nnz = 0; 1431 iia[0] = 0; 1432 for (i=0; i<na; i++) { /* fill adjacency */ 1433 cnt = 0; 1434 len = ia[i+1] - ia[i]; 1435 row = ja + ia[i]; 1436 for (j=0; j<len; j++) { 1437 if (row[j] != i) jja[nnz+cnt++] = row[j]; /* if not diagonal */ 1438 } 1439 nnz += cnt; 1440 iia[i+1] = nnz; 1441 } 1442 /* Partitioning of the adjacency matrix */ 1443 ierr = MatCreateMPIAdj(PETSC_COMM_SELF,na,na,iia,jja,NULL,&adj);CHKERRQ(ierr); 1444 ierr = MatPartitioningSetAdjacency(mpart,adj);CHKERRQ(ierr); 1445 ierr = MatPartitioningSetNParts(mpart,nloc);CHKERRQ(ierr); 1446 ierr = MatPartitioningApply(mpart,&ispart);CHKERRQ(ierr); 1447 ierr = ISPartitioningToNumbering(ispart,&isnumb);CHKERRQ(ierr); 1448 ierr = MatDestroy(&adj);CHKERRQ(ierr); 1449 foundpart = PETSC_TRUE; 1450 } 1451 ierr = MatRestoreRowIJ(Ad,0,PETSC_TRUE,isbaij,&na,&ia,&ja,&done);CHKERRQ(ierr); 1452 } 1453 ierr = MatPartitioningDestroy(&mpart);CHKERRQ(ierr); 1454 } 1455 ierr = PetscMalloc1(nloc,&is);CHKERRQ(ierr); 1456 if (!foundpart) { 1457 1458 /* Partitioning by contiguous chunks of rows */ 1459 1460 PetscInt mbs = (rend-rstart)/bs; 1461 PetscInt start = rstart; 1462 for (i=0; i<nloc; i++) { 1463 PetscInt count = (mbs/nloc + ((mbs % nloc) > i)) * bs; 1464 ierr = ISCreateStride(PETSC_COMM_SELF,count,start,1,&is[i]);CHKERRQ(ierr); 1465 start += count; 1466 } 1467 1468 } else { 1469 1470 /* Partitioning by adjacency of diagonal block */ 1471 1472 const PetscInt *numbering; 1473 PetscInt *count,nidx,*indices,*newidx,start=0; 1474 /* Get node count in each partition */ 1475 ierr = PetscMalloc1(nloc,&count);CHKERRQ(ierr); 1476 ierr = ISPartitioningCount(ispart,nloc,count);CHKERRQ(ierr); 1477 if (isbaij && bs > 1) { /* adjust for the block-aij case */ 1478 for (i=0; i<nloc; i++) count[i] *= bs; 1479 } 1480 /* Build indices from node numbering */ 1481 ierr = ISGetLocalSize(isnumb,&nidx);CHKERRQ(ierr); 1482 ierr = PetscMalloc1(nidx,&indices);CHKERRQ(ierr); 1483 for (i=0; i<nidx; i++) indices[i] = i; /* needs to be initialized */ 1484 ierr = ISGetIndices(isnumb,&numbering);CHKERRQ(ierr); 1485 ierr = PetscSortIntWithPermutation(nidx,numbering,indices);CHKERRQ(ierr); 1486 ierr = ISRestoreIndices(isnumb,&numbering);CHKERRQ(ierr); 1487 if (isbaij && bs > 1) { /* adjust for the block-aij case */ 1488 ierr = PetscMalloc1(nidx*bs,&newidx);CHKERRQ(ierr); 1489 for (i=0; i<nidx; i++) { 1490 for (j=0; j<bs; j++) newidx[i*bs+j] = indices[i]*bs + j; 1491 } 1492 ierr = PetscFree(indices);CHKERRQ(ierr); 1493 nidx *= bs; 1494 indices = newidx; 1495 } 1496 /* Shift to get global indices */ 1497 for (i=0; i<nidx; i++) indices[i] += rstart; 1498 1499 /* Build the index sets for each block */ 1500 for (i=0; i<nloc; i++) { 1501 ierr = ISCreateGeneral(PETSC_COMM_SELF,count[i],&indices[start],PETSC_COPY_VALUES,&is[i]);CHKERRQ(ierr); 1502 ierr = ISSort(is[i]);CHKERRQ(ierr); 1503 start += count[i]; 1504 } 1505 1506 ierr = PetscFree(count);CHKERRQ(ierr); 1507 ierr = PetscFree(indices);CHKERRQ(ierr); 1508 ierr = ISDestroy(&isnumb);CHKERRQ(ierr); 1509 ierr = ISDestroy(&ispart);CHKERRQ(ierr); 1510 } 1511 *iis = is; 1512 PetscFunctionReturn(0); 1513 } 1514 1515 #undef __FUNCT__ 1516 #define __FUNCT__ "PCGASMCreateStraddlingSubdomains" 1517 PETSC_INTERN PetscErrorCode PCGASMCreateStraddlingSubdomains(Mat A,PetscInt N,PetscInt *n,IS *iis[]) 1518 { 1519 PetscErrorCode ierr; 1520 1521 PetscFunctionBegin; 1522 ierr = MatSubdomainsCreateCoalesce(A,N,n,iis);CHKERRQ(ierr); 1523 PetscFunctionReturn(0); 1524 } 1525 1526 1527 1528 #undef __FUNCT__ 1529 #define __FUNCT__ "PCGASMCreateSubdomains" 1530 /*@C 1531 PCGASMCreateSubdomains - Creates n index sets defining n nonoverlapping subdomains for the additive 1532 Schwarz preconditioner for a any problem based on its matrix. 1533 1534 Collective 1535 1536 Input Parameters: 1537 + A - The global matrix operator 1538 - N - the number of global subdomains requested 1539 1540 Output Parameters: 1541 + n - the number of subdomains created on this processor 1542 - iis - the array of index sets defining the local inner subdomains (on which the correction is applied) 1543 1544 Level: advanced 1545 1546 Note: When N >= A's communicator size, each subdomain is local -- contained within a single processor. 1547 When N < size, the subdomains are 'straddling' (processor boundaries) and are no longer local. 1548 The resulting subdomains can be use in PCGASMSetSubdomains(pc,n,iss,NULL). The overlapping 1549 outer subdomains will be automatically generated from these according to the requested amount of 1550 overlap; this is currently supported only with local subdomains. 1551 1552 1553 .keywords: PC, GASM, additive Schwarz, create, subdomains, unstructured grid 1554 1555 .seealso: PCGASMSetSubdomains(), PCGASMDestroySubdomains() 1556 @*/ 1557 PetscErrorCode PCGASMCreateSubdomains(Mat A,PetscInt N,PetscInt *n,IS *iis[]) 1558 { 1559 PetscMPIInt size; 1560 PetscErrorCode ierr; 1561 1562 PetscFunctionBegin; 1563 PetscValidHeaderSpecific(A,MAT_CLASSID,1); 1564 PetscValidPointer(iis,4); 1565 1566 if (N < 1) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Number of subdomains must be > 0, N = %D",N); 1567 ierr = MPI_Comm_size(PetscObjectComm((PetscObject)A),&size);CHKERRQ(ierr); 1568 if (N >= size) { 1569 *n = N/size + (N%size); 1570 ierr = PCGASMCreateLocalSubdomains(A,*n,iis);CHKERRQ(ierr); 1571 } else { 1572 ierr = PCGASMCreateStraddlingSubdomains(A,N,n,iis);CHKERRQ(ierr); 1573 } 1574 PetscFunctionReturn(0); 1575 } 1576 1577 #undef __FUNCT__ 1578 #define __FUNCT__ "PCGASMDestroySubdomains" 1579 /*@C 1580 PCGASMDestroySubdomains - Destroys the index sets created with 1581 PCGASMCreateSubdomains() or PCGASMCreateSubdomains2D. Should be 1582 called after setting subdomains with PCGASMSetSubdomains(). 1583 1584 Collective 1585 1586 Input Parameters: 1587 + n - the number of index sets 1588 . iis - the array of inner subdomains, 1589 - ois - the array of outer subdomains, can be NULL 1590 1591 Level: intermediate 1592 1593 Notes: this is merely a convenience subroutine that walks each list, 1594 destroys each IS on the list, and then frees the list. At the end the 1595 list pointers are set to NULL. 1596 1597 .keywords: PC, GASM, additive Schwarz, create, subdomains, unstructured grid 1598 1599 .seealso: PCGASMCreateSubdomains(), PCGASMSetSubdomains() 1600 @*/ 1601 PetscErrorCode PCGASMDestroySubdomains(PetscInt n,IS **iis,IS **ois) 1602 { 1603 PetscInt i; 1604 PetscErrorCode ierr; 1605 1606 PetscFunctionBegin; 1607 if (n <= 0) PetscFunctionReturn(0); 1608 if (ois) { 1609 PetscValidPointer(ois,3); 1610 if (*ois) { 1611 PetscValidPointer(*ois,3); 1612 for (i=0; i<n; i++) { 1613 ierr = ISDestroy(&(*ois)[i]);CHKERRQ(ierr); 1614 } 1615 ierr = PetscFree((*ois));CHKERRQ(ierr); 1616 } 1617 } 1618 if (iis) { 1619 PetscValidPointer(iis,2); 1620 if (*iis) { 1621 PetscValidPointer(*iis,2); 1622 for (i=0; i<n; i++) { 1623 ierr = ISDestroy(&(*iis)[i]);CHKERRQ(ierr); 1624 } 1625 ierr = PetscFree((*iis));CHKERRQ(ierr); 1626 } 1627 } 1628 PetscFunctionReturn(0); 1629 } 1630 1631 1632 #define PCGASMLocalSubdomainBounds2D(M,N,xleft,ylow,xright,yhigh,first,last,xleft_loc,ylow_loc,xright_loc,yhigh_loc,n) \ 1633 { \ 1634 PetscInt first_row = first/M, last_row = last/M+1; \ 1635 /* \ 1636 Compute ylow_loc and yhigh_loc so that (ylow_loc,xleft) and (yhigh_loc,xright) are the corners \ 1637 of the bounding box of the intersection of the subdomain with the local ownership range (local \ 1638 subdomain). \ 1639 Also compute xleft_loc and xright_loc as the lower and upper bounds on the first and last rows \ 1640 of the intersection. \ 1641 */ \ 1642 /* ylow_loc is the grid row containing the first element of the local sumbdomain */ \ 1643 *ylow_loc = PetscMax(first_row,ylow); \ 1644 /* xleft_loc is the offset of first element of the local subdomain within its grid row (might actually be outside the local subdomain) */ \ 1645 *xleft_loc = *ylow_loc==first_row ? PetscMax(first%M,xleft) : xleft; \ 1646 /* yhigh_loc is the grid row above the last local subdomain element */ \ 1647 *yhigh_loc = PetscMin(last_row,yhigh); \ 1648 /* xright is the offset of the end of the local subdomain within its grid row (might actually be outside the local subdomain) */ \ 1649 *xright_loc = *yhigh_loc==last_row ? PetscMin(xright,last%M) : xright; \ 1650 /* Now compute the size of the local subdomain n. */ \ 1651 *n = 0; \ 1652 if (*ylow_loc < *yhigh_loc) { \ 1653 PetscInt width = xright-xleft; \ 1654 *n += width*(*yhigh_loc-*ylow_loc-1); \ 1655 *n += PetscMin(PetscMax(*xright_loc-xleft,0),width); \ 1656 *n -= PetscMin(PetscMax(*xleft_loc-xleft,0), width); \ 1657 } \ 1658 } 1659 1660 1661 1662 #undef __FUNCT__ 1663 #define __FUNCT__ "PCGASMCreateSubdomains2D" 1664 /*@ 1665 PCGASMCreateSubdomains2D - Creates the index sets for the overlapping Schwarz 1666 preconditioner for a two-dimensional problem on a regular grid. 1667 1668 Collective 1669 1670 Input Parameters: 1671 + M, N - the global number of grid points in the x and y directions 1672 . Mdomains, Ndomains - the global number of subdomains in the x and y directions 1673 . dof - degrees of freedom per node 1674 - overlap - overlap in mesh lines 1675 1676 Output Parameters: 1677 + Nsub - the number of local subdomains created 1678 . iis - array of index sets defining inner (nonoverlapping) subdomains 1679 - ois - array of index sets defining outer (overlapping, if overlap > 0) subdomains 1680 1681 1682 Level: advanced 1683 1684 .keywords: PC, GASM, additive Schwarz, create, subdomains, 2D, regular grid 1685 1686 .seealso: PCGASMSetSubdomains(), PCGASMGetSubKSP(), PCGASMSetOverlap() 1687 @*/ 1688 PetscErrorCode PCGASMCreateSubdomains2D(PC pc,PetscInt M,PetscInt N,PetscInt Mdomains,PetscInt Ndomains,PetscInt dof,PetscInt overlap,PetscInt *nsub,IS **iis,IS **ois) 1689 { 1690 PetscErrorCode ierr; 1691 PetscMPIInt size, rank; 1692 PetscInt i, j; 1693 PetscInt maxheight, maxwidth; 1694 PetscInt xstart, xleft, xright, xleft_loc, xright_loc; 1695 PetscInt ystart, ylow, yhigh, ylow_loc, yhigh_loc; 1696 PetscInt x[2][2], y[2][2], n[2]; 1697 PetscInt first, last; 1698 PetscInt nidx, *idx; 1699 PetscInt ii,jj,s,q,d; 1700 PetscInt k,kk; 1701 PetscMPIInt color; 1702 MPI_Comm comm, subcomm; 1703 IS **xis = 0, **is = ois, **is_local = iis; 1704 1705 PetscFunctionBegin; 1706 ierr = PetscObjectGetComm((PetscObject)pc, &comm);CHKERRQ(ierr); 1707 ierr = MPI_Comm_size(comm, &size);CHKERRQ(ierr); 1708 ierr = MPI_Comm_rank(comm, &rank);CHKERRQ(ierr); 1709 ierr = MatGetOwnershipRange(pc->pmat, &first, &last);CHKERRQ(ierr); 1710 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) " 1711 "does not respect the number of degrees of freedom per grid point %D", first, last, dof); 1712 1713 /* Determine the number of domains with nonzero intersections with the local ownership range. */ 1714 s = 0; 1715 ystart = 0; 1716 for (j=0; j<Ndomains; ++j) { 1717 maxheight = N/Ndomains + ((N % Ndomains) > j); /* Maximal height of subdomain */ 1718 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); 1719 /* Vertical domain limits with an overlap. */ 1720 ylow = PetscMax(ystart - overlap,0); 1721 yhigh = PetscMin(ystart + maxheight + overlap,N); 1722 xstart = 0; 1723 for (i=0; i<Mdomains; ++i) { 1724 maxwidth = M/Mdomains + ((M % Mdomains) > i); /* Maximal width of subdomain */ 1725 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); 1726 /* Horizontal domain limits with an overlap. */ 1727 xleft = PetscMax(xstart - overlap,0); 1728 xright = PetscMin(xstart + maxwidth + overlap,M); 1729 /* 1730 Determine whether this subdomain intersects this processor's ownership range of pc->pmat. 1731 */ 1732 PCGASMLocalSubdomainBounds2D(M,N,xleft,ylow,xright,yhigh,first,last,(&xleft_loc),(&ylow_loc),(&xright_loc),(&yhigh_loc),(&nidx)); 1733 if (nidx) ++s; 1734 xstart += maxwidth; 1735 } /* for (i = 0; i < Mdomains; ++i) */ 1736 ystart += maxheight; 1737 } /* for (j = 0; j < Ndomains; ++j) */ 1738 1739 /* Now we can allocate the necessary number of ISs. */ 1740 *nsub = s; 1741 ierr = PetscMalloc1(*nsub,is);CHKERRQ(ierr); 1742 ierr = PetscMalloc1(*nsub,is_local);CHKERRQ(ierr); 1743 s = 0; 1744 ystart = 0; 1745 for (j=0; j<Ndomains; ++j) { 1746 maxheight = N/Ndomains + ((N % Ndomains) > j); /* Maximal height of subdomain */ 1747 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); 1748 /* Vertical domain limits with an overlap. */ 1749 y[0][0] = PetscMax(ystart - overlap,0); 1750 y[0][1] = PetscMin(ystart + maxheight + overlap,N); 1751 /* Vertical domain limits without an overlap. */ 1752 y[1][0] = ystart; 1753 y[1][1] = PetscMin(ystart + maxheight,N); 1754 xstart = 0; 1755 for (i=0; i<Mdomains; ++i) { 1756 maxwidth = M/Mdomains + ((M % Mdomains) > i); /* Maximal width of subdomain */ 1757 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); 1758 /* Horizontal domain limits with an overlap. */ 1759 x[0][0] = PetscMax(xstart - overlap,0); 1760 x[0][1] = PetscMin(xstart + maxwidth + overlap,M); 1761 /* Horizontal domain limits without an overlap. */ 1762 x[1][0] = xstart; 1763 x[1][1] = PetscMin(xstart+maxwidth,M); 1764 /* 1765 Determine whether this domain intersects this processor's ownership range of pc->pmat. 1766 Do this twice: first for the domains with overlaps, and once without. 1767 During the first pass create the subcommunicators, and use them on the second pass as well. 1768 */ 1769 for (q = 0; q < 2; ++q) { 1770 PetscBool split = PETSC_FALSE; 1771 /* 1772 domain limits, (xleft, xright) and (ylow, yheigh) are adjusted 1773 according to whether the domain with an overlap or without is considered. 1774 */ 1775 xleft = x[q][0]; xright = x[q][1]; 1776 ylow = y[q][0]; yhigh = y[q][1]; 1777 PCGASMLocalSubdomainBounds2D(M,N,xleft,ylow,xright,yhigh,first,last,(&xleft_loc),(&ylow_loc),(&xright_loc),(&yhigh_loc),(&nidx)); 1778 nidx *= dof; 1779 n[q] = nidx; 1780 /* 1781 Based on the counted number of indices in the local domain *with an overlap*, 1782 construct a subcommunicator of all the processors supporting this domain. 1783 Observe that a domain with an overlap might have nontrivial local support, 1784 while the domain without an overlap might not. Hence, the decision to participate 1785 in the subcommunicator must be based on the domain with an overlap. 1786 */ 1787 if (q == 0) { 1788 if (nidx) color = 1; 1789 else color = MPI_UNDEFINED; 1790 ierr = MPI_Comm_split(comm, color, rank, &subcomm);CHKERRQ(ierr); 1791 split = PETSC_TRUE; 1792 } 1793 /* 1794 Proceed only if the number of local indices *with an overlap* is nonzero. 1795 */ 1796 if (n[0]) { 1797 if (q == 0) xis = is; 1798 if (q == 1) { 1799 /* 1800 The IS for the no-overlap subdomain shares a communicator with the overlapping domain. 1801 Moreover, if the overlap is zero, the two ISs are identical. 1802 */ 1803 if (overlap == 0) { 1804 (*is_local)[s] = (*is)[s]; 1805 ierr = PetscObjectReference((PetscObject)(*is)[s]);CHKERRQ(ierr); 1806 continue; 1807 } else { 1808 xis = is_local; 1809 subcomm = ((PetscObject)(*is)[s])->comm; 1810 } 1811 } /* if (q == 1) */ 1812 idx = NULL; 1813 ierr = PetscMalloc1(nidx,&idx);CHKERRQ(ierr); 1814 if (nidx) { 1815 k = 0; 1816 for (jj=ylow_loc; jj<yhigh_loc; ++jj) { 1817 PetscInt x0 = (jj==ylow_loc) ? xleft_loc : xleft; 1818 PetscInt x1 = (jj==yhigh_loc-1) ? xright_loc : xright; 1819 kk = dof*(M*jj + x0); 1820 for (ii=x0; ii<x1; ++ii) { 1821 for (d = 0; d < dof; ++d) { 1822 idx[k++] = kk++; 1823 } 1824 } 1825 } 1826 } 1827 ierr = ISCreateGeneral(subcomm,nidx,idx,PETSC_OWN_POINTER,(*xis)+s);CHKERRQ(ierr); 1828 if (split) { 1829 ierr = MPI_Comm_free(&subcomm);CHKERRQ(ierr); 1830 } 1831 }/* if (n[0]) */ 1832 }/* for (q = 0; q < 2; ++q) */ 1833 if (n[0]) ++s; 1834 xstart += maxwidth; 1835 } /* for (i = 0; i < Mdomains; ++i) */ 1836 ystart += maxheight; 1837 } /* for (j = 0; j < Ndomains; ++j) */ 1838 PetscFunctionReturn(0); 1839 } 1840 1841 #undef __FUNCT__ 1842 #define __FUNCT__ "PCGASMGetSubdomains" 1843 /*@C 1844 PCGASMGetSubdomains - Gets the subdomains supported on this processor 1845 for the additive Schwarz preconditioner. 1846 1847 Not Collective 1848 1849 Input Parameter: 1850 . pc - the preconditioner context 1851 1852 Output Parameters: 1853 + n - the number of subdomains for this processor (default value = 1) 1854 . iis - the index sets that define the inner subdomains (without overlap) supported on this processor (can be NULL) 1855 - ois - the index sets that define the outer subdomains (with overlap) supported on this processor (can be NULL) 1856 1857 1858 Notes: 1859 The user is responsible for destroying the ISs and freeing the returned arrays. 1860 The IS numbering is in the parallel, global numbering of the vector. 1861 1862 Level: advanced 1863 1864 .keywords: PC, GASM, get, subdomains, additive Schwarz 1865 1866 .seealso: PCGASMSetOverlap(), PCGASMGetSubKSP(), PCGASMCreateSubdomains2D(), 1867 PCGASMSetSubdomains(), PCGASMGetSubmatrices() 1868 @*/ 1869 PetscErrorCode PCGASMGetSubdomains(PC pc,PetscInt *n,IS *iis[],IS *ois[]) 1870 { 1871 PC_GASM *osm; 1872 PetscErrorCode ierr; 1873 PetscBool match; 1874 PetscInt i; 1875 1876 PetscFunctionBegin; 1877 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1878 ierr = PetscObjectTypeCompare((PetscObject)pc,PCGASM,&match);CHKERRQ(ierr); 1879 if (!match) SETERRQ2(PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_WRONG, "Incorrect object type: expected %s, got %s instead", PCGASM, ((PetscObject)pc)->type_name); 1880 osm = (PC_GASM*)pc->data; 1881 if (n) *n = osm->n; 1882 if (iis) { 1883 ierr = PetscMalloc1(osm->n, iis);CHKERRQ(ierr); 1884 } 1885 if (ois) { 1886 ierr = PetscMalloc1(osm->n, ois);CHKERRQ(ierr); 1887 } 1888 if (iis || ois) { 1889 for (i = 0; i < osm->n; ++i) { 1890 if (iis) (*iis)[i] = osm->iis[i]; 1891 if (ois) (*ois)[i] = osm->ois[i]; 1892 } 1893 } 1894 PetscFunctionReturn(0); 1895 } 1896 1897 #undef __FUNCT__ 1898 #define __FUNCT__ "PCGASMGetSubmatrices" 1899 /*@C 1900 PCGASMGetSubmatrices - Gets the local submatrices (for this processor 1901 only) for the additive Schwarz preconditioner. 1902 1903 Not Collective 1904 1905 Input Parameter: 1906 . pc - the preconditioner context 1907 1908 Output Parameters: 1909 + n - the number of matrices for this processor (default value = 1) 1910 - mat - the matrices 1911 1912 Notes: matrices returned by this routine have the same communicators as the index sets (IS) 1913 used to define subdomains in PCGASMSetSubdomains() 1914 Level: advanced 1915 1916 .keywords: PC, GASM, set, local, subdomains, additive Schwarz, block Jacobi 1917 1918 .seealso: PCGASMSetOverlap(), PCGASMGetSubKSP(), 1919 PCGASMCreateSubdomains2D(), PCGASMSetSubdomains(), PCGASMGetSubdomains() 1920 @*/ 1921 PetscErrorCode PCGASMGetSubmatrices(PC pc,PetscInt *n,Mat *mat[]) 1922 { 1923 PC_GASM *osm; 1924 PetscErrorCode ierr; 1925 PetscBool match; 1926 1927 PetscFunctionBegin; 1928 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1929 PetscValidIntPointer(n,2); 1930 if (mat) PetscValidPointer(mat,3); 1931 if (!pc->setupcalled) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONGSTATE,"Must call after KSPSetUP() or PCSetUp()."); 1932 ierr = PetscObjectTypeCompare((PetscObject)pc,PCGASM,&match);CHKERRQ(ierr); 1933 if (!match) SETERRQ2(PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_WRONG, "Expected %s, got %s instead", PCGASM, ((PetscObject)pc)->type_name); 1934 osm = (PC_GASM*)pc->data; 1935 if (n) *n = osm->n; 1936 if (mat) *mat = osm->pmat; 1937 PetscFunctionReturn(0); 1938 } 1939 1940 #undef __FUNCT__ 1941 #define __FUNCT__ "PCGASMSetUseDMSubdomains" 1942 /*@ 1943 PCGASMSetUseDMSubdomains - Indicates whether to use DMCreateDomainDecomposition() to define the subdomains, whenever possible. 1944 Logically Collective 1945 1946 Input Parameter: 1947 + pc - the preconditioner 1948 - flg - boolean indicating whether to use subdomains defined by the DM 1949 1950 Options Database Key: 1951 . -pc_gasm_dm_subdomains -pc_gasm_overlap -pc_gasm_total_subdomains 1952 1953 Level: intermediate 1954 1955 Notes: 1956 PCGASMSetSubdomains(), PCGASMSetTotalSubdomains() or PCGASMSetOverlap() take precedence over PCGASMSetUseDMSubdomains(), 1957 so setting PCGASMSetSubdomains() with nontrivial subdomain ISs or any of PCGASMSetTotalSubdomains() and PCGASMSetOverlap() 1958 automatically turns the latter off. 1959 1960 .keywords: PC, ASM, DM, set, subdomains, additive Schwarz 1961 1962 .seealso: PCGASMGetUseDMSubdomains(), PCGASMSetSubdomains(), PCGASMSetOverlap() 1963 PCGASMCreateSubdomains2D() 1964 @*/ 1965 PetscErrorCode PCGASMSetUseDMSubdomains(PC pc,PetscBool flg) 1966 { 1967 PC_GASM *osm = (PC_GASM*)pc->data; 1968 PetscErrorCode ierr; 1969 PetscBool match; 1970 1971 PetscFunctionBegin; 1972 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1973 PetscValidLogicalCollectiveBool(pc,flg,2); 1974 if (pc->setupcalled) SETERRQ(((PetscObject)pc)->comm,PETSC_ERR_ARG_WRONGSTATE,"Not for a setup PC."); 1975 ierr = PetscObjectTypeCompare((PetscObject)pc,PCGASM,&match);CHKERRQ(ierr); 1976 if (match) { 1977 if (!osm->user_subdomains && osm->N == PETSC_DETERMINE && osm->overlap < 0) { 1978 osm->dm_subdomains = flg; 1979 } 1980 } 1981 PetscFunctionReturn(0); 1982 } 1983 1984 #undef __FUNCT__ 1985 #define __FUNCT__ "PCGASMGetUseDMSubdomains" 1986 /*@ 1987 PCGASMGetUseDMSubdomains - Returns flag indicating whether to use DMCreateDomainDecomposition() to define the subdomains, whenever possible. 1988 Not Collective 1989 1990 Input Parameter: 1991 . pc - the preconditioner 1992 1993 Output Parameter: 1994 . flg - boolean indicating whether to use subdomains defined by the DM 1995 1996 Level: intermediate 1997 1998 .keywords: PC, ASM, DM, set, subdomains, additive Schwarz 1999 2000 .seealso: PCGASMSetUseDMSubdomains(), PCGASMSetOverlap() 2001 PCGASMCreateSubdomains2D() 2002 @*/ 2003 PetscErrorCode PCGASMGetUseDMSubdomains(PC pc,PetscBool* flg) 2004 { 2005 PC_GASM *osm = (PC_GASM*)pc->data; 2006 PetscErrorCode ierr; 2007 PetscBool match; 2008 2009 PetscFunctionBegin; 2010 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 2011 PetscValidPointer(flg,2); 2012 ierr = PetscObjectTypeCompare((PetscObject)pc,PCGASM,&match);CHKERRQ(ierr); 2013 if (match) { 2014 if (flg) *flg = osm->dm_subdomains; 2015 } 2016 PetscFunctionReturn(0); 2017 } 2018