1 #define PETSCKSP_DLL 2 3 /* 4 This file defines an "generalized" additive Schwarz preconditioner for any Mat implementation. 5 In this version each processor may have any number of subdomains and a subdomain may live on multiple 6 processors. 7 8 N - total number of subdomains on all processors 9 n - actual number of subdomains on this processor 10 nmax - maximum number of subdomains per processor 11 */ 12 #include "private/pcimpl.h" /*I "petscpc.h" I*/ 13 14 typedef struct { 15 PetscInt N,n,nmax; 16 PetscInt overlap; /* overlap requested by user */ 17 KSP *ksp; /* linear solvers for each block */ 18 Vec gx,gy; /* Merged work vectors */ 19 Vec *x,*y; /* Split work vectors; storage aliases pieces of storage of the above merged vectors. */ 20 IS gis, gis_local; /* merged ISs */ 21 VecScatter grestriction; /* merged restriction */ 22 VecScatter gprolongation; /* merged prolongation */ 23 IS *is; /* index set that defines each overlapping subdomain */ 24 IS *is_local; /* index set that defines each local subdomain (same as subdomain with the overlap removed); may be NULL */ 25 Mat *pmat; /* subdomain block matrices */ 26 PCGASMType type; /* use reduced interpolation, restriction or both */ 27 PetscBool type_set; /* if user set this value (so won't change it for symmetric problems) */ 28 PetscBool same_local_solves; /* flag indicating whether all local solvers are same */ 29 PetscBool sort_indices; /* flag to sort subdomain indices */ 30 } PC_GASM; 31 32 #undef __FUNCT__ 33 #define __FUNCT__ "PCGASMPrintSubdomains" 34 static PetscErrorCode PCGASMPrintSubdomains(PC pc) 35 { 36 PC_GASM *osm = (PC_GASM*)pc->data; 37 const char *prefix; 38 char fname[PETSC_MAX_PATH_LEN+1]; 39 PetscViewer viewer, sviewer; 40 PetscInt i,j,nidx; 41 const PetscInt *idx; 42 char *cidx; 43 PetscBool found; 44 PetscErrorCode ierr; 45 46 PetscFunctionBegin; 47 ierr = PCGetOptionsPrefix(pc,&prefix);CHKERRQ(ierr); 48 ierr = PetscOptionsGetString(prefix,"-pc_gasm_print_subdomains",fname,PETSC_MAX_PATH_LEN,&found);CHKERRQ(ierr); 49 if (!found) { ierr = PetscStrcpy(fname,"stdout");CHKERRQ(ierr); }; 50 for (i=0;i<osm->n;++i) { 51 ierr = PetscViewerASCIIOpen(((PetscObject)((osm->is)[i]))->comm,fname,&viewer);CHKERRQ(ierr); 52 ierr = ISGetLocalSize(osm->is[i], &nidx); CHKERRQ(ierr); 53 ierr = ISGetIndices(osm->is[i], &idx); CHKERRQ(ierr); 54 ierr = PetscMalloc(sizeof(char)*(16*nidx+1), &cidx); CHKERRQ(ierr); /* no more than 15 characters per index plus a space*/ 55 ierr = PetscViewerStringOpen(((PetscObject)(osm->is[i]))->comm, cidx, 16*nidx, &sviewer); CHKERRQ(ierr); 56 for(j = 0; j < nidx; ++j) { 57 ierr = PetscViewerStringSPrintf(sviewer, "%D ", idx[j]); CHKERRQ(ierr); 58 } 59 ierr = PetscViewerDestroy(sviewer); CHKERRQ(ierr); 60 ierr = ISRestoreIndices(osm->is[i],&idx); CHKERRQ(ierr); 61 ierr = PetscViewerASCIIPrintf(viewer, "Subdomain with overlap\n"); CHKERRQ(ierr); 62 ierr = PetscViewerASCIISynchronizedPrintf(viewer, "%s", cidx); CHKERRQ(ierr); 63 ierr = PetscViewerFlush(viewer); CHKERRQ(ierr); 64 ierr = PetscViewerASCIIPrintf(viewer, "\n"); CHKERRQ(ierr); 65 ierr = PetscFree(cidx); CHKERRQ(ierr); 66 ierr = PetscViewerDestroy(viewer); CHKERRQ(ierr); 67 if (osm->is_local) { 68 ierr = PetscViewerASCIIOpen(((PetscObject)((osm->is)[i]))->comm,fname,&viewer);CHKERRQ(ierr); 69 ierr = ISGetLocalSize(osm->is_local[i], &nidx); CHKERRQ(ierr); 70 ierr = ISGetIndices(osm->is_local[i], &idx); CHKERRQ(ierr); 71 ierr = PetscMalloc(sizeof(char)*(16*nidx+1), &cidx); CHKERRQ(ierr); /* no more than 15 characters per index plus a space*/ 72 ierr = PetscViewerStringOpen(((PetscObject)(osm->is[i]))->comm, cidx, 16*nidx, &sviewer); CHKERRQ(ierr); 73 for(j = 0; j < nidx; ++j) { 74 ierr = PetscViewerStringSPrintf(sviewer, "%D ", idx[j]); CHKERRQ(ierr); 75 } 76 ierr = PetscViewerDestroy(sviewer); CHKERRQ(ierr); 77 ierr = ISRestoreIndices(osm->is_local[i],&idx); CHKERRQ(ierr); 78 ierr = PetscViewerASCIIPrintf(viewer, "Subdomain without overlap\n"); CHKERRQ(ierr); 79 ierr = PetscViewerASCIISynchronizedPrintf(viewer, "%s", cidx); CHKERRQ(ierr); 80 ierr = PetscViewerFlush(viewer); CHKERRQ(ierr); 81 ierr = PetscViewerASCIIPrintf(viewer, "\n"); CHKERRQ(ierr); 82 ierr = PetscFree(cidx); CHKERRQ(ierr); 83 ierr = PetscViewerDestroy(viewer); CHKERRQ(ierr); 84 } 85 } 86 PetscFunctionReturn(0); 87 } 88 89 90 #undef __FUNCT__ 91 #define __FUNCT__ "PCView_GASM" 92 static PetscErrorCode PCView_GASM(PC pc,PetscViewer viewer) 93 { 94 PC_GASM *osm = (PC_GASM*)pc->data; 95 const char *prefix; 96 PetscErrorCode ierr; 97 PetscMPIInt rank, size; 98 PetscInt i,bsz; 99 PetscBool iascii,isstring, print_subdomains; 100 PetscViewer sviewer; 101 102 103 PetscFunctionBegin; 104 ierr = MPI_Comm_size(((PetscObject)pc)->comm, &size); CHKERRQ(ierr); 105 ierr = MPI_Comm_rank(((PetscObject)pc)->comm, &rank); CHKERRQ(ierr); 106 ierr = PetscTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr); 107 ierr = PetscTypeCompare((PetscObject)viewer,PETSCVIEWERSTRING,&isstring);CHKERRQ(ierr); 108 ierr = PCGetOptionsPrefix(pc,&prefix);CHKERRQ(ierr); 109 ierr = PetscOptionsGetBool(prefix,"-pc_gasm_print_subdomains",&print_subdomains,PETSC_NULL);CHKERRQ(ierr); 110 if (iascii) { 111 char overlaps[256] = "user-defined overlap",subdomains[256] = "total subdomains set"; 112 if (osm->overlap >= 0) {ierr = PetscSNPrintf(overlaps,sizeof overlaps,"amount of overlap = %D",osm->overlap);CHKERRQ(ierr);} 113 if (osm->nmax > 0) {ierr = PetscSNPrintf(subdomains,sizeof subdomains,"max number of local subdomains = %D",osm->nmax);CHKERRQ(ierr);} 114 ierr = PetscViewerASCIISynchronizedPrintf(viewer," [%d:%d] number of locally subdomains = %D\n",(int)rank,(int)size,osm->n);CHKERRQ(ierr); 115 ierr = PetscViewerFlush(viewer); CHKERRQ(ierr); 116 ierr = PetscViewerASCIIPrintf(viewer," Generalized additive Schwarz: %s, %s\n",subdomains,overlaps);CHKERRQ(ierr); 117 ierr = PetscViewerASCIIPrintf(viewer," Generalized additive Schwarz: restriction/interpolation type - %s\n",PCGASMTypes[osm->type]);CHKERRQ(ierr); 118 ierr = MPI_Comm_rank(((PetscObject)pc)->comm,&rank);CHKERRQ(ierr); 119 if (osm->same_local_solves) { 120 if (osm->ksp) { 121 ierr = PetscViewerASCIIPrintf(viewer," Local solve is same for all subdomains, in the following KSP and PC objects:\n");CHKERRQ(ierr); 122 ierr = PetscViewerGetSingleton(viewer,&sviewer);CHKERRQ(ierr); 123 if (!rank) { 124 ierr = PetscViewerASCIIPushTab(viewer);CHKERRQ(ierr); 125 ierr = KSPView(osm->ksp[0],sviewer);CHKERRQ(ierr); 126 ierr = PetscViewerASCIIPopTab(viewer);CHKERRQ(ierr); 127 } 128 ierr = PetscViewerRestoreSingleton(viewer,&sviewer);CHKERRQ(ierr); 129 } 130 } else { 131 ierr = PetscViewerFlush(viewer);CHKERRQ(ierr); 132 ierr = PetscViewerASCIIPrintf(viewer," Local solve info for each subdomain is in the following KSP and PC objects:\n");CHKERRQ(ierr); 133 ierr = PetscViewerASCIIPushTab(viewer);CHKERRQ(ierr); 134 ierr = PetscViewerASCIIPrintf(viewer,"- - - - - - - - - - - - - - - - - -\n");CHKERRQ(ierr); 135 for (i=0; i<osm->nmax; i++) { 136 ierr = PetscViewerGetSingleton(viewer,&sviewer);CHKERRQ(ierr); 137 if (i < osm->n) { 138 ierr = ISGetLocalSize(osm->is[i],&bsz);CHKERRQ(ierr); 139 ierr = PetscViewerASCIISynchronizedPrintf(sviewer,"[%d:%d] local subdomain number %D, size = %D\n",(int)rank,(int)size,i,bsz);CHKERRQ(ierr); 140 ierr = KSPView(osm->ksp[i],sviewer);CHKERRQ(ierr); 141 ierr = PetscViewerASCIISynchronizedPrintf(sviewer,"- - - - - - - - - - - - - - - - - -\n");CHKERRQ(ierr); 142 } 143 ierr = PetscViewerRestoreSingleton(viewer,&sviewer);CHKERRQ(ierr); 144 } 145 ierr = PetscViewerASCIIPopTab(viewer);CHKERRQ(ierr); 146 ierr = PetscViewerFlush(viewer);CHKERRQ(ierr); 147 } 148 } else if (isstring) { 149 ierr = PetscViewerStringSPrintf(viewer," subdomains=%D, overlap=%D, type=%s",osm->n,osm->overlap,PCGASMTypes[osm->type]);CHKERRQ(ierr); 150 ierr = PetscViewerGetSingleton(viewer,&sviewer);CHKERRQ(ierr); 151 if (osm->ksp) {ierr = KSPView(osm->ksp[0],sviewer);CHKERRQ(ierr);} 152 ierr = PetscViewerRestoreSingleton(viewer,&sviewer);CHKERRQ(ierr); 153 } else { 154 SETERRQ1(((PetscObject)pc)->comm,PETSC_ERR_SUP,"Viewer type %s not supported for PCGASM",((PetscObject)viewer)->type_name); 155 } 156 ierr = PetscViewerFlush(viewer); CHKERRQ(ierr); 157 if(print_subdomains) { 158 ierr = PCGASMPrintSubdomains(pc); CHKERRQ(ierr); 159 } 160 PetscFunctionReturn(0); 161 } 162 163 164 165 166 167 #undef __FUNCT__ 168 #define __FUNCT__ "PCSetUp_GASM" 169 static PetscErrorCode PCSetUp_GASM(PC pc) 170 { 171 PC_GASM *osm = (PC_GASM*)pc->data; 172 PetscErrorCode ierr; 173 PetscBool symset,flg; 174 PetscInt i; 175 PetscMPIInt rank, size; 176 MatReuse scall = MAT_REUSE_MATRIX; 177 KSP ksp; 178 PC subpc; 179 const char *prefix,*pprefix; 180 PetscInt dn; /* Number of indices in a single subdomain assigned to this processor. */ 181 const PetscInt *didx; /* Indices from a single subdomain assigned to this processor. */ 182 PetscInt ddn; /* Number of indices in all subdomains assigned to this processor. */ 183 PetscInt *ddidx; /* Indices of all subdomains assigned to this processor. */ 184 IS gid; /* Identity IS of the size of all subdomains assigned to this processor. */ 185 Vec x,y; 186 PetscScalar *gxarray, *gyarray; 187 PetscInt gfirst, glast; 188 189 PetscFunctionBegin; 190 ierr = MPI_Comm_size(((PetscObject)pc)->comm,&size);CHKERRQ(ierr); 191 ierr = MPI_Comm_rank(((PetscObject)pc)->comm,&rank);CHKERRQ(ierr); 192 if (!pc->setupcalled) { 193 194 if (!osm->type_set) { 195 ierr = MatIsSymmetricKnown(pc->pmat,&symset,&flg);CHKERRQ(ierr); 196 if (symset && flg) { osm->type = PC_GASM_BASIC; } 197 } 198 199 if (osm->N == PETSC_DECIDE && osm->n < 1) { 200 /* no subdomains given, use one per processor */ 201 osm->nmax = osm->n = 1; 202 ierr = MPI_Comm_size(((PetscObject)pc)->comm,&size);CHKERRQ(ierr); 203 osm->N = size; 204 } else if (osm->N == PETSC_DECIDE) { 205 PetscInt inwork[2], outwork[2]; 206 /* determine global number of subdomains and the max number of local subdomains */ 207 inwork[0] = inwork[1] = osm->n; 208 ierr = MPI_Allreduce(inwork,outwork,1,MPIU_2INT,PetscMaxSum_Op,((PetscObject)pc)->comm);CHKERRQ(ierr); 209 osm->nmax = outwork[0]; 210 osm->N = outwork[1]; 211 } 212 if (!osm->is){ /* create the index sets */ 213 ierr = PCGASMCreateSubdomains(pc->pmat,osm->n,&osm->is);CHKERRQ(ierr); 214 } 215 if (!osm->is_local) { 216 /* 217 This indicates that osm->is should define a nonoverlapping decomposition 218 (there is no way to really guarantee that if subdomains are set by the user through PCGASMSetLocalSubdomains, 219 but the assumption is that either the user does the right thing, or subdomains in ossm->is have been created 220 via PCGASMCreateSubdomains, which guarantees a nonoverlapping decomposition). 221 Therefore, osm->is will be used to define osm->is_local. 222 If a nonzero overlap has been requested by the user, then osm->is will be expanded and will overlap, 223 so osm->is_local should obtain a copy of osm->is while they are still (presumably) nonoverlapping. 224 Otherwise (no overlap has been requested), osm->is_local are simply aliases for osm->is. 225 */ 226 ierr = PetscMalloc(osm->n*sizeof(IS),&osm->is_local);CHKERRQ(ierr); 227 for (i=0; i<osm->n; i++) { 228 if (osm->overlap > 0) { /* With positive overlap, osm->is[i] will be modified */ 229 ierr = ISDuplicate(osm->is[i],&osm->is_local[i]);CHKERRQ(ierr); 230 ierr = ISCopy(osm->is[i],osm->is_local[i]);CHKERRQ(ierr); 231 } else { 232 ierr = PetscObjectReference((PetscObject)osm->is[i]);CHKERRQ(ierr); 233 osm->is_local[i] = osm->is[i]; 234 } 235 } 236 } 237 /* Beyond this point osm->is_local is not null. */ 238 if (osm->overlap > 0) { 239 /* Extend the "overlapping" regions by a number of steps */ 240 ierr = MatIncreaseOverlap(pc->pmat,osm->n,osm->is,osm->overlap);CHKERRQ(ierr); 241 } 242 ierr = PCGetOptionsPrefix(pc,&prefix);CHKERRQ(ierr); 243 flg = PETSC_FALSE; 244 ierr = PetscOptionsGetBool(prefix,"-pc_gasm_print_subdomains",&flg,PETSC_NULL);CHKERRQ(ierr); 245 if (flg) { ierr = PCGASMPrintSubdomains(pc);CHKERRQ(ierr); } 246 247 if (osm->sort_indices) { 248 for (i=0; i<osm->n; i++) { 249 ierr = ISSort(osm->is[i]);CHKERRQ(ierr); 250 ierr = ISSort(osm->is_local[i]);CHKERRQ(ierr); 251 } 252 } 253 /* Merge the ISs, create merged vectors and scatter contexts. */ 254 /* Restriction ISs. */ 255 ddn = 0; 256 for (i=0; i<osm->n; i++) { 257 ierr = ISGetLocalSize(osm->is[i],&dn); CHKERRQ(ierr); 258 ddn += dn; 259 } 260 ierr = PetscMalloc(ddn*sizeof(PetscInt), &ddidx); CHKERRQ(ierr); 261 ddn = 0; 262 for (i=0; i<osm->n; i++) { 263 ierr = ISGetLocalSize(osm->is[i],&dn); CHKERRQ(ierr); 264 ierr = ISGetIndices(osm->is[i],&didx); CHKERRQ(ierr); 265 ierr = PetscMemcpy(ddidx+ddn, didx, sizeof(PetscInt)*dn); CHKERRQ(ierr); 266 ierr = ISRestoreIndices(osm->is[i], &didx); CHKERRQ(ierr); 267 ddn += dn; 268 } 269 ierr = ISCreateGeneral(((PetscObject)(pc))->comm, ddn, ddidx, PETSC_OWN_POINTER, &osm->gis); CHKERRQ(ierr); 270 ierr = MatGetVecs(pc->pmat,&x,&y);CHKERRQ(ierr); 271 ierr = VecCreateMPI(((PetscObject)pc)->comm, ddn, PETSC_DECIDE, &osm->gx); CHKERRQ(ierr); 272 ierr = VecDuplicate(osm->gx,&osm->gy); CHKERRQ(ierr); 273 ierr = VecGetOwnershipRange(osm->gx, &gfirst, &glast); CHKERRQ(ierr); 274 ierr = ISCreateStride(((PetscObject)pc)->comm,ddn,gfirst,1, &gid); CHKERRQ(ierr); 275 ierr = VecScatterCreate(x,osm->gis,osm->gx,gid, &(osm->grestriction)); CHKERRQ(ierr); 276 ierr = ISDestroy(gid); CHKERRQ(ierr); 277 /* Prolongation ISs */ 278 { PetscInt dn_local; /* Number of indices in the local part of a single domain assigned to this processor. */ 279 const PetscInt *didx_local; /* Global indices from the local part of a single domain assigned to this processor. */ 280 PetscInt ddn_local; /* Number of indices in the local part of the disjoint union all domains assigned to this processor. */ 281 PetscInt *ddidx_local; /* Global indices of the local part of the disjoint union of all domains assigned to this processor. */ 282 /**/ 283 ISLocalToGlobalMapping ltog; /* Map from global to local indices on the disjoint union of subdomains: "local" ind's run from 0 to ddn-1. */ 284 PetscInt *ddidx_llocal; /* Mapped local indices of the disjoint union of local parts of subdomains. */ 285 PetscInt ddn_llocal; /* Number of indices in ddidx_llocal; must equal ddn_local, or else gis_local is not a sub-IS of gis. */ 286 IS gis_llocal; /* IS with ddidx_llocal indices. */ 287 PetscInt j; 288 ddn_local = 0; 289 for (i=0; i<osm->n; i++) { 290 ierr = ISGetLocalSize(osm->is_local[i],&dn_local); CHKERRQ(ierr); 291 ddn_local += dn_local; 292 } 293 ierr = PetscMalloc(ddn_local*sizeof(PetscInt), &ddidx_local); CHKERRQ(ierr); 294 ddn_local = 0; 295 for (i=0; i<osm->n; i++) { 296 ierr = ISGetLocalSize(osm->is_local[i],&dn_local); CHKERRQ(ierr); 297 ierr = ISGetIndices(osm->is_local[i],&didx_local); CHKERRQ(ierr); 298 ierr = PetscMemcpy(ddidx_local+ddn_local, didx_local, sizeof(PetscInt)*dn_local); CHKERRQ(ierr); 299 ierr = ISRestoreIndices(osm->is_local[i], &didx_local); CHKERRQ(ierr); 300 ddn_local += dn_local; 301 } 302 ierr = PetscMalloc(sizeof(PetscInt)*ddn_local, &ddidx_llocal); CHKERRQ(ierr); 303 ierr = ISCreateGeneral(((PetscObject)pc)->comm, ddn_local, ddidx_local, PETSC_OWN_POINTER, &(osm->gis_local)); CHKERRQ(ierr); 304 ierr = ISLocalToGlobalMappingCreateIS(osm->gis,<og);CHKERRQ(ierr); 305 ierr = ISGlobalToLocalMappingApply(ltog,IS_GTOLM_DROP,ddn_local,ddidx_local,&ddn_llocal,ddidx_llocal);CHKERRQ(ierr); 306 ierr = ISLocalToGlobalMappingDestroy(ltog);CHKERRQ(ierr); 307 if (ddn_llocal != ddn_local) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_PLIB,"gis_local contains %D indices outside of gis", ddn_llocal - ddn_local); 308 /* Now convert these localized indices into the global indices into the merged output vector. */ 309 ierr = VecGetOwnershipRange(osm->gy, &gfirst, &glast); CHKERRQ(ierr); 310 for(j=0; j < ddn_llocal; ++j) { 311 ddidx_llocal[j] += gfirst; 312 } 313 ierr = ISCreateGeneral(PETSC_COMM_SELF,ddn_llocal,ddidx_llocal,PETSC_OWN_POINTER,&gis_llocal); CHKERRQ(ierr); 314 ierr = VecScatterCreate(y,osm->gis_local,osm->gy,gis_llocal,&osm->gprolongation); CHKERRQ(ierr); 315 ierr = ISDestroy(gis_llocal);CHKERRQ(ierr); 316 } 317 /* Create the subdomain work vectors. */ 318 ierr = PetscMalloc(osm->n*sizeof(Vec),&osm->x);CHKERRQ(ierr); 319 ierr = PetscMalloc(osm->n*sizeof(Vec),&osm->y);CHKERRQ(ierr); 320 ierr = VecGetOwnershipRange(osm->gx, &gfirst, &glast);CHKERRQ(ierr); 321 ierr = VecGetArray(osm->gx, &gxarray); CHKERRQ(ierr); 322 ierr = VecGetArray(osm->gy, &gyarray); CHKERRQ(ierr); 323 for (i=0, ddn=0; i<osm->n; ++i, ddn += dn) { 324 ierr = ISGetLocalSize(osm->is[i],&dn);CHKERRQ(ierr); 325 ierr = VecCreateSeqWithArray(PETSC_COMM_SELF,dn,gxarray+ddn,&osm->x[i]);CHKERRQ(ierr); 326 ierr = VecCreateSeqWithArray(PETSC_COMM_SELF,dn,gyarray+ddn,&osm->y[i]);CHKERRQ(ierr); 327 } 328 ierr = VecRestoreArray(osm->gx, &gxarray); CHKERRQ(ierr); 329 ierr = VecRestoreArray(osm->gy, &gyarray); CHKERRQ(ierr); 330 ierr = VecDestroy(x); CHKERRQ(ierr); 331 ierr = VecDestroy(y); CHKERRQ(ierr); 332 /* Create the local solvers */ 333 ierr = PetscMalloc(osm->n*sizeof(KSP *),&osm->ksp);CHKERRQ(ierr); 334 for (i=0; i<osm->n; i++) { /* KSPs are local */ 335 ierr = KSPCreate(PETSC_COMM_SELF,&ksp);CHKERRQ(ierr); 336 ierr = PetscLogObjectParent(pc,ksp);CHKERRQ(ierr); 337 ierr = PetscObjectIncrementTabLevel((PetscObject)ksp,(PetscObject)pc,1);CHKERRQ(ierr); 338 ierr = KSPSetType(ksp,KSPPREONLY);CHKERRQ(ierr); 339 ierr = KSPGetPC(ksp,&subpc);CHKERRQ(ierr); 340 ierr = PCGetOptionsPrefix(pc,&prefix);CHKERRQ(ierr); 341 ierr = KSPSetOptionsPrefix(ksp,prefix);CHKERRQ(ierr); 342 ierr = KSPAppendOptionsPrefix(ksp,"sub_");CHKERRQ(ierr); 343 osm->ksp[i] = ksp; 344 } 345 scall = MAT_INITIAL_MATRIX; 346 347 } else { 348 /* 349 Destroy the blocks from the previous iteration 350 */ 351 if (pc->flag == DIFFERENT_NONZERO_PATTERN) { 352 ierr = MatDestroyMatrices(osm->n,&osm->pmat);CHKERRQ(ierr); 353 scall = MAT_INITIAL_MATRIX; 354 } 355 } 356 357 /* 358 Extract out the submatrices. 359 */ 360 if(size > 1) { 361 ierr = MatGetSubMatricesParallel(pc->pmat,osm->n,osm->is, osm->is,scall,&osm->pmat);CHKERRQ(ierr); 362 } 363 else { 364 ierr = MatGetSubMatrices(pc->pmat,osm->n,osm->is, osm->is,scall,&osm->pmat);CHKERRQ(ierr); 365 } 366 if (scall == MAT_INITIAL_MATRIX) { 367 ierr = PetscObjectGetOptionsPrefix((PetscObject)pc->pmat,&pprefix);CHKERRQ(ierr); 368 for (i=0; i<osm->n; i++) { 369 ierr = PetscLogObjectParent(pc,osm->pmat[i]);CHKERRQ(ierr); 370 ierr = PetscObjectSetOptionsPrefix((PetscObject)osm->pmat[i],pprefix);CHKERRQ(ierr); 371 } 372 } 373 374 /* Return control to the user so that the submatrices can be modified (e.g., to apply 375 different boundary conditions for the submatrices than for the global problem) */ 376 ierr = PCModifySubMatrices(pc,osm->n,osm->is,osm->is,osm->pmat,pc->modifysubmatricesP);CHKERRQ(ierr); 377 378 /* 379 Loop over submatrices putting them into local ksp 380 */ 381 for (i=0; i<osm->n; i++) { 382 ierr = KSPSetOperators(osm->ksp[i],osm->pmat[i],osm->pmat[i],pc->flag);CHKERRQ(ierr); 383 if (!pc->setupcalled) { 384 ierr = KSPSetFromOptions(osm->ksp[i]);CHKERRQ(ierr); 385 } 386 } 387 388 PetscFunctionReturn(0); 389 } 390 391 #undef __FUNCT__ 392 #define __FUNCT__ "PCSetUpOnBlocks_GASM" 393 static PetscErrorCode PCSetUpOnBlocks_GASM(PC pc) 394 { 395 PC_GASM *osm = (PC_GASM*)pc->data; 396 PetscErrorCode ierr; 397 PetscInt i; 398 399 PetscFunctionBegin; 400 for (i=0; i<osm->n; i++) { 401 ierr = KSPSetUp(osm->ksp[i]);CHKERRQ(ierr); 402 } 403 PetscFunctionReturn(0); 404 } 405 406 #undef __FUNCT__ 407 #define __FUNCT__ "PCApply_GASM" 408 static PetscErrorCode PCApply_GASM(PC pc,Vec x,Vec y) 409 { 410 PC_GASM *osm = (PC_GASM*)pc->data; 411 PetscErrorCode ierr; 412 PetscInt i; 413 ScatterMode forward = SCATTER_FORWARD,reverse = SCATTER_REVERSE; 414 415 PetscFunctionBegin; 416 /* 417 Support for limiting the restriction or interpolation to only local 418 subdomain values (leaving the other values 0). 419 */ 420 if (!(osm->type & PC_GASM_RESTRICT)) { 421 forward = SCATTER_FORWARD_LOCAL; 422 /* have to zero the work RHS since scatter may leave some slots empty */ 423 ierr = VecZeroEntries(osm->gx); CHKERRQ(ierr); 424 } 425 if (!(osm->type & PC_GASM_INTERPOLATE)) { 426 reverse = SCATTER_REVERSE_LOCAL; 427 } 428 429 ierr = VecScatterBegin(osm->grestriction,x,osm->gx,INSERT_VALUES,forward);CHKERRQ(ierr); 430 ierr = VecZeroEntries(y);CHKERRQ(ierr); 431 ierr = VecScatterEnd(osm->grestriction,x,osm->gx,INSERT_VALUES,forward);CHKERRQ(ierr); 432 /* do the local solves */ 433 for (i=0; i<osm->n; ++i) { /* Note that the solves are local, so we can go to osm->n, rather than osm->nmax. */ 434 ierr = KSPSolve(osm->ksp[i],osm->x[i],osm->y[i]);CHKERRQ(ierr); 435 } 436 ierr = VecScatterBegin(osm->gprolongation,osm->gy,y,ADD_VALUES,reverse);CHKERRQ(ierr); 437 ierr = VecScatterEnd(osm->gprolongation,osm->gy,y,ADD_VALUES,reverse);CHKERRQ(ierr); 438 PetscFunctionReturn(0); 439 } 440 441 #undef __FUNCT__ 442 #define __FUNCT__ "PCApplyTranspose_GASM" 443 static PetscErrorCode PCApplyTranspose_GASM(PC pc,Vec x,Vec y) 444 { 445 PC_GASM *osm = (PC_GASM*)pc->data; 446 PetscErrorCode ierr; 447 PetscInt i; 448 ScatterMode forward = SCATTER_FORWARD,reverse = SCATTER_REVERSE; 449 450 PetscFunctionBegin; 451 /* 452 Support for limiting the restriction or interpolation to only local 453 subdomain values (leaving the other values 0). 454 455 Note: these are reversed from the PCApply_GASM() because we are applying the 456 transpose of the three terms 457 */ 458 if (!(osm->type & PC_GASM_INTERPOLATE)) { 459 forward = SCATTER_FORWARD_LOCAL; 460 /* have to zero the work RHS since scatter may leave some slots empty */ 461 ierr = VecZeroEntries(osm->gx);CHKERRQ(ierr); 462 } 463 if (!(osm->type & PC_GASM_RESTRICT)) { 464 reverse = SCATTER_REVERSE_LOCAL; 465 } 466 467 ierr = VecScatterBegin(osm->grestriction,x,osm->gx,INSERT_VALUES,forward);CHKERRQ(ierr); 468 ierr = VecZeroEntries(y);CHKERRQ(ierr); 469 ierr = VecScatterEnd(osm->grestriction,x,osm->gx,INSERT_VALUES,forward);CHKERRQ(ierr); 470 /* do the local solves */ 471 for (i=0; i<osm->n; ++i) { /* Note that the solves are local, so we can go to osm->n, rather than osm->nmax. */ 472 ierr = KSPSolveTranspose(osm->ksp[i],osm->x[i],osm->y[i]);CHKERRQ(ierr); 473 } 474 ierr = VecScatterBegin(osm->gprolongation,osm->gy,y,ADD_VALUES,reverse);CHKERRQ(ierr); 475 ierr = VecScatterEnd(osm->gprolongation,osm->gy,y,ADD_VALUES,reverse);CHKERRQ(ierr); 476 PetscFunctionReturn(0); 477 } 478 479 #undef __FUNCT__ 480 #define __FUNCT__ "PCDestroy_GASM" 481 static PetscErrorCode PCDestroy_GASM(PC pc) 482 { 483 PC_GASM *osm = (PC_GASM*)pc->data; 484 PetscErrorCode ierr; 485 PetscInt i; 486 487 PetscFunctionBegin; 488 if (osm->ksp) { 489 for (i=0; i<osm->n; i++) { 490 ierr = KSPDestroy(osm->ksp[i]);CHKERRQ(ierr); 491 } 492 ierr = PetscFree(osm->ksp);CHKERRQ(ierr); 493 } 494 if (osm->pmat) { 495 if (osm->n > 0) { 496 ierr = MatDestroyMatrices(osm->n,&osm->pmat);CHKERRQ(ierr); 497 } 498 } 499 for (i=0; i<osm->n; i++) { 500 ierr = VecDestroy(osm->x[i]);CHKERRQ(ierr); 501 ierr = VecDestroy(osm->y[i]);CHKERRQ(ierr); 502 } 503 ierr = PetscFree(osm->x);CHKERRQ(ierr); 504 ierr = PetscFree(osm->y);CHKERRQ(ierr); 505 ierr = VecDestroy(osm->gx); CHKERRQ(ierr); 506 ierr = VecDestroy(osm->gy); CHKERRQ(ierr); 507 508 ierr = VecScatterDestroy(osm->grestriction); CHKERRQ(ierr); 509 ierr = VecScatterDestroy(osm->gprolongation); CHKERRQ(ierr); 510 511 if (osm->is) { 512 ierr = PCGASMDestroySubdomains(osm->n,osm->is,osm->is_local);CHKERRQ(ierr); 513 } 514 ierr = PetscFree(osm);CHKERRQ(ierr); 515 PetscFunctionReturn(0); 516 } 517 518 #undef __FUNCT__ 519 #define __FUNCT__ "PCSetFromOptions_GASM" 520 static PetscErrorCode PCSetFromOptions_GASM(PC pc) { 521 PC_GASM *osm = (PC_GASM*)pc->data; 522 PetscErrorCode ierr; 523 PetscInt blocks,ovl; 524 PetscBool symset,flg; 525 PCGASMType gasmtype; 526 527 PetscFunctionBegin; 528 /* set the type to symmetric if matrix is symmetric */ 529 if (!osm->type_set && pc->pmat) { 530 ierr = MatIsSymmetricKnown(pc->pmat,&symset,&flg);CHKERRQ(ierr); 531 if (symset && flg) { osm->type = PC_GASM_BASIC; } 532 } 533 ierr = PetscOptionsHead("Generalized additive Schwarz options");CHKERRQ(ierr); 534 ierr = PetscOptionsInt("-pc_gasm_blocks","Number of subdomains","PCGASMSetTotalSubdomains",osm->n,&blocks,&flg);CHKERRQ(ierr); 535 if (flg) {ierr = PCGASMSetTotalSubdomains(pc,blocks);CHKERRQ(ierr); } 536 ierr = PetscOptionsInt("-pc_gasm_overlap","Number of grid points overlap","PCGASMSetOverlap",osm->overlap,&ovl,&flg);CHKERRQ(ierr); 537 if (flg) {ierr = PCGASMSetOverlap(pc,ovl);CHKERRQ(ierr); } 538 flg = PETSC_FALSE; 539 ierr = PetscOptionsEnum("-pc_gasm_type","Type of restriction/extension","PCGASMSetType",PCGASMTypes,(PetscEnum)osm->type,(PetscEnum*)&gasmtype,&flg);CHKERRQ(ierr); 540 if (flg) {ierr = PCGASMSetType(pc,gasmtype);CHKERRQ(ierr); } 541 ierr = PetscOptionsTail();CHKERRQ(ierr); 542 PetscFunctionReturn(0); 543 } 544 545 /*------------------------------------------------------------------------------------*/ 546 547 EXTERN_C_BEGIN 548 #undef __FUNCT__ 549 #define __FUNCT__ "PCGASMSetLocalSubdomains_GASM" 550 PetscErrorCode PETSCKSP_DLLEXPORT PCGASMSetLocalSubdomains_GASM(PC pc,PetscInt n,IS is[],IS is_local[]) 551 { 552 PC_GASM *osm = (PC_GASM*)pc->data; 553 PetscErrorCode ierr; 554 PetscInt i; 555 556 PetscFunctionBegin; 557 if (n < 1) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Each process must have 1 or more blocks, n = %D",n); 558 if (pc->setupcalled && (n != osm->n || is)) SETERRQ(((PetscObject)pc)->comm,PETSC_ERR_ARG_WRONGSTATE,"PCGASMSetLocalSubdomains() should be called before calling PCSetUp()."); 559 560 if (!pc->setupcalled) { 561 osm->n = n; 562 osm->is = 0; 563 osm->is_local = 0; 564 if (is) { 565 for (i=0; i<n; i++) {ierr = PetscObjectReference((PetscObject)is[i]);CHKERRQ(ierr);} 566 } 567 if (is_local) { 568 for (i=0; i<n; i++) {ierr = PetscObjectReference((PetscObject)is_local[i]);CHKERRQ(ierr);} 569 } 570 if (osm->is) { 571 ierr = PCGASMDestroySubdomains(osm->n,osm->is,osm->is_local);CHKERRQ(ierr); 572 } 573 if (is) { 574 ierr = PetscMalloc(n*sizeof(IS),&osm->is);CHKERRQ(ierr); 575 for (i=0; i<n; i++) { osm->is[i] = is[i]; } 576 /* Flag indicating that the user has set overlapping subdomains so PCGASM should not increase their size. */ 577 osm->overlap = -1; 578 } 579 if (is_local) { 580 ierr = PetscMalloc(n*sizeof(IS),&osm->is_local);CHKERRQ(ierr); 581 for (i=0; i<n; i++) { osm->is_local[i] = is_local[i]; } 582 } 583 } 584 PetscFunctionReturn(0); 585 } 586 EXTERN_C_END 587 588 EXTERN_C_BEGIN 589 #undef __FUNCT__ 590 #define __FUNCT__ "PCGASMSetTotalSubdomains_GASM" 591 PetscErrorCode PETSCKSP_DLLEXPORT PCGASMSetTotalSubdomains_GASM(PC pc,PetscInt N) { 592 PC_GASM *osm = (PC_GASM*)pc->data; 593 PetscErrorCode ierr; 594 PetscMPIInt rank,size; 595 PetscInt n; 596 597 PetscFunctionBegin; 598 if (N < 1) SETERRQ1(((PetscObject)pc)->comm,PETSC_ERR_ARG_OUTOFRANGE,"Number of total blocks must be > 0, N = %D",N); 599 600 /* 601 Split the subdomains equally among all processors 602 */ 603 ierr = MPI_Comm_rank(((PetscObject)pc)->comm,&rank);CHKERRQ(ierr); 604 ierr = MPI_Comm_size(((PetscObject)pc)->comm,&size);CHKERRQ(ierr); 605 n = N/size + ((N % size) > rank); 606 if (!n) SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Process %d must have at least one block: total processors %d total blocks %D",(int)rank,(int)size,N); 607 if (pc->setupcalled && n != osm->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"PCGASMSetTotalSubdomains() should be called before PCSetUp()."); 608 if (!pc->setupcalled) { 609 if (osm->is) { 610 ierr = PCGASMDestroySubdomains(osm->n,osm->is,osm->is_local);CHKERRQ(ierr); 611 } 612 osm->N = N; 613 osm->n = n; 614 osm->is = 0; 615 osm->is_local = 0; 616 } 617 PetscFunctionReturn(0); 618 }/* PCGASMSetTotalSubdomains_GASM() */ 619 EXTERN_C_END 620 621 EXTERN_C_BEGIN 622 #undef __FUNCT__ 623 #define __FUNCT__ "PCGASMSetOverlap_GASM" 624 PetscErrorCode PETSCKSP_DLLEXPORT PCGASMSetOverlap_GASM(PC pc,PetscInt ovl) 625 { 626 PC_GASM *osm = (PC_GASM*)pc->data; 627 628 PetscFunctionBegin; 629 if (ovl < 0) SETERRQ(((PetscObject)pc)->comm,PETSC_ERR_ARG_OUTOFRANGE,"Negative overlap value requested"); 630 if (pc->setupcalled && ovl != osm->overlap) SETERRQ(((PetscObject)pc)->comm,PETSC_ERR_ARG_WRONGSTATE,"PCGASMSetOverlap() should be called before PCSetUp()."); 631 if (!pc->setupcalled) { 632 osm->overlap = ovl; 633 } 634 PetscFunctionReturn(0); 635 } 636 EXTERN_C_END 637 638 EXTERN_C_BEGIN 639 #undef __FUNCT__ 640 #define __FUNCT__ "PCGASMSetType_GASM" 641 PetscErrorCode PETSCKSP_DLLEXPORT PCGASMSetType_GASM(PC pc,PCGASMType type) 642 { 643 PC_GASM *osm = (PC_GASM*)pc->data; 644 645 PetscFunctionBegin; 646 osm->type = type; 647 osm->type_set = PETSC_TRUE; 648 PetscFunctionReturn(0); 649 } 650 EXTERN_C_END 651 652 EXTERN_C_BEGIN 653 #undef __FUNCT__ 654 #define __FUNCT__ "PCGASMSetSortIndices_GASM" 655 PetscErrorCode PETSCKSP_DLLEXPORT PCGASMSetSortIndices_GASM(PC pc,PetscBool doSort) 656 { 657 PC_GASM *osm = (PC_GASM*)pc->data; 658 659 PetscFunctionBegin; 660 osm->sort_indices = doSort; 661 PetscFunctionReturn(0); 662 } 663 EXTERN_C_END 664 665 EXTERN_C_BEGIN 666 #undef __FUNCT__ 667 #define __FUNCT__ "PCGASMGetSubKSP_GASM" 668 /* 669 FIX: This routine might need to be modified once multiple ranks per subdomain are allowed. 670 In particular, it would upset the global subdomain number calculation. 671 */ 672 PetscErrorCode PETSCKSP_DLLEXPORT PCGASMGetSubKSP_GASM(PC pc,PetscInt *n,PetscInt *first,KSP **ksp) 673 { 674 PC_GASM *osm = (PC_GASM*)pc->data; 675 PetscErrorCode ierr; 676 677 PetscFunctionBegin; 678 if (osm->n < 1) SETERRQ(((PetscObject)pc)->comm,PETSC_ERR_ORDER,"Need to call PCSetUP() on PC (or KSPSetUp() on the outer KSP object) before calling here"); 679 680 if (n) { 681 *n = osm->n; 682 } 683 if (first) { 684 ierr = MPI_Scan(&osm->n,first,1,MPIU_INT,MPI_SUM,((PetscObject)pc)->comm);CHKERRQ(ierr); 685 *first -= osm->n; 686 } 687 if (ksp) { 688 /* Assume that local solves are now different; not necessarily 689 true though! This flag is used only for PCView_GASM() */ 690 *ksp = osm->ksp; 691 osm->same_local_solves = PETSC_FALSE; 692 } 693 PetscFunctionReturn(0); 694 }/* PCGASMGetSubKSP_GASM() */ 695 EXTERN_C_END 696 697 698 #undef __FUNCT__ 699 #define __FUNCT__ "PCGASMSetLocalSubdomains" 700 /*@C 701 PCGASMSetLocalSubdomains - Sets the local subdomains (for this processor 702 only) for the additive Schwarz preconditioner. 703 704 Collective on PC 705 706 Input Parameters: 707 + pc - the preconditioner context 708 . n - the number of subdomains for this processor (default value = 1) 709 . is - the index set that defines the subdomains for this processor 710 (or PETSC_NULL for PETSc to determine subdomains) 711 - is_local - the index sets that define the local part of the subdomains for this processor 712 (or PETSC_NULL to use the default of 1 subdomain per process) 713 714 Notes: 715 The IS numbering is in the parallel, global numbering of the vector. 716 717 By default the GASM preconditioner uses 1 block per processor. 718 719 Use PCGASMSetTotalSubdomains() to set the subdomains for all processors. 720 721 Level: advanced 722 723 .keywords: PC, GASM, set, local, subdomains, additive Schwarz 724 725 .seealso: PCGASMSetTotalSubdomains(), PCGASMSetOverlap(), PCGASMGetSubKSP(), 726 PCGASMCreateSubdomains2D(), PCGASMGetLocalSubdomains() 727 @*/ 728 PetscErrorCode PETSCKSP_DLLEXPORT PCGASMSetLocalSubdomains(PC pc,PetscInt n,IS is[],IS is_local[]) 729 { 730 PetscErrorCode ierr; 731 732 PetscFunctionBegin; 733 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 734 ierr = PetscTryMethod(pc,"PCGASMSetLocalSubdomains_C",(PC,PetscInt,IS[],IS[]),(pc,n,is,is_local));CHKERRQ(ierr); 735 PetscFunctionReturn(0); 736 } 737 738 #undef __FUNCT__ 739 #define __FUNCT__ "PCGASMSetTotalSubdomains" 740 /*@C 741 PCGASMSetTotalSubdomains - Sets the subdomains for all processor for the 742 additive Schwarz preconditioner. Either all or no processors in the 743 PC communicator must call this routine, with the same index sets. 744 745 Collective on PC 746 747 Input Parameters: 748 + pc - the preconditioner context 749 . n - the number of subdomains for all processors 750 . is - the index sets that define the subdomains for all processor 751 (or PETSC_NULL for PETSc to determine subdomains) 752 - is_local - the index sets that define the local part of the subdomains for this processor 753 (or PETSC_NULL to use the default of 1 subdomain per process) 754 755 Options Database Key: 756 To set the total number of subdomain blocks rather than specify the 757 index sets, use the option 758 . -pc_gasm_blocks <blks> - Sets total blocks 759 760 Notes: 761 Currently you cannot use this to set the actual subdomains with the argument is. 762 763 By default the GASM preconditioner uses 1 block per processor. 764 765 These index sets cannot be destroyed until after completion of the 766 linear solves for which the GASM preconditioner is being used. 767 768 Use PCGASMSetLocalSubdomains() to set local subdomains. 769 770 Level: advanced 771 772 .keywords: PC, GASM, set, total, global, subdomains, additive Schwarz 773 774 .seealso: PCGASMSetLocalSubdomains(), PCGASMSetOverlap(), PCGASMGetSubKSP(), 775 PCGASMCreateSubdomains2D() 776 @*/ 777 PetscErrorCode PETSCKSP_DLLEXPORT PCGASMSetTotalSubdomains(PC pc,PetscInt N) 778 { 779 PetscErrorCode ierr; 780 781 PetscFunctionBegin; 782 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 783 ierr = PetscTryMethod(pc,"PCGASMSetTotalSubdomains_C",(PC,PetscInt),(pc,N));CHKERRQ(ierr); 784 PetscFunctionReturn(0); 785 } 786 787 #undef __FUNCT__ 788 #define __FUNCT__ "PCGASMSetOverlap" 789 /*@ 790 PCGASMSetOverlap - Sets the overlap between a pair of subdomains for the 791 additive Schwarz preconditioner. Either all or no processors in the 792 PC communicator must call this routine. 793 794 Logically Collective on PC 795 796 Input Parameters: 797 + pc - the preconditioner context 798 - ovl - the amount of overlap between subdomains (ovl >= 0, default value = 1) 799 800 Options Database Key: 801 . -pc_gasm_overlap <ovl> - Sets overlap 802 803 Notes: 804 By default the GASM preconditioner uses 1 block per processor. To use 805 multiple blocks per perocessor, see PCGASMSetTotalSubdomains() and 806 PCGASMSetLocalSubdomains() (and the option -pc_gasm_blocks <blks>). 807 808 The overlap defaults to 1, so if one desires that no additional 809 overlap be computed beyond what may have been set with a call to 810 PCGASMSetTotalSubdomains() or PCGASMSetLocalSubdomains(), then ovl 811 must be set to be 0. In particular, if one does not explicitly set 812 the subdomains an application code, then all overlap would be computed 813 internally by PETSc, and using an overlap of 0 would result in an GASM 814 variant that is equivalent to the block Jacobi preconditioner. 815 816 Note that one can define initial index sets with any overlap via 817 PCGASMSetTotalSubdomains() or PCGASMSetLocalSubdomains(); the routine 818 PCGASMSetOverlap() merely allows PETSc to extend that overlap further 819 if desired. 820 821 Level: intermediate 822 823 .keywords: PC, GASM, set, overlap 824 825 .seealso: PCGASMSetTotalSubdomains(), PCGASMSetLocalSubdomains(), PCGASMGetSubKSP(), 826 PCGASMCreateSubdomains2D(), PCGASMGetLocalSubdomains() 827 @*/ 828 PetscErrorCode PETSCKSP_DLLEXPORT PCGASMSetOverlap(PC pc,PetscInt ovl) 829 { 830 PetscErrorCode ierr; 831 832 PetscFunctionBegin; 833 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 834 PetscValidLogicalCollectiveInt(pc,ovl,2); 835 ierr = PetscTryMethod(pc,"PCGASMSetOverlap_C",(PC,PetscInt),(pc,ovl));CHKERRQ(ierr); 836 PetscFunctionReturn(0); 837 } 838 839 #undef __FUNCT__ 840 #define __FUNCT__ "PCGASMSetType" 841 /*@ 842 PCGASMSetType - Sets the type of restriction and interpolation used 843 for local problems in the additive Schwarz method. 844 845 Logically Collective on PC 846 847 Input Parameters: 848 + pc - the preconditioner context 849 - type - variant of GASM, one of 850 .vb 851 PC_GASM_BASIC - full interpolation and restriction 852 PC_GASM_RESTRICT - full restriction, local processor interpolation 853 PC_GASM_INTERPOLATE - full interpolation, local processor restriction 854 PC_GASM_NONE - local processor restriction and interpolation 855 .ve 856 857 Options Database Key: 858 . -pc_gasm_type [basic,restrict,interpolate,none] - Sets GASM type 859 860 Level: intermediate 861 862 .keywords: PC, GASM, set, type 863 864 .seealso: PCGASMSetTotalSubdomains(), PCGASMSetTotalSubdomains(), PCGASMGetSubKSP(), 865 PCGASMCreateSubdomains2D() 866 @*/ 867 PetscErrorCode PETSCKSP_DLLEXPORT PCGASMSetType(PC pc,PCGASMType type) 868 { 869 PetscErrorCode ierr; 870 871 PetscFunctionBegin; 872 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 873 PetscValidLogicalCollectiveEnum(pc,type,2); 874 ierr = PetscTryMethod(pc,"PCGASMSetType_C",(PC,PCGASMType),(pc,type));CHKERRQ(ierr); 875 PetscFunctionReturn(0); 876 } 877 878 #undef __FUNCT__ 879 #define __FUNCT__ "PCGASMSetSortIndices" 880 /*@ 881 PCGASMSetSortIndices - Determines whether subdomain indices are sorted. 882 883 Logically Collective on PC 884 885 Input Parameters: 886 + pc - the preconditioner context 887 - doSort - sort the subdomain indices 888 889 Level: intermediate 890 891 .keywords: PC, GASM, set, type 892 893 .seealso: PCGASMSetLocalSubdomains(), PCGASMSetTotalSubdomains(), PCGASMGetSubKSP(), 894 PCGASMCreateSubdomains2D() 895 @*/ 896 PetscErrorCode PETSCKSP_DLLEXPORT PCGASMSetSortIndices(PC pc,PetscBool doSort) 897 { 898 PetscErrorCode ierr; 899 900 PetscFunctionBegin; 901 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 902 PetscValidLogicalCollectiveBool(pc,doSort,2); 903 ierr = PetscTryMethod(pc,"PCGASMSetSortIndices_C",(PC,PetscBool),(pc,doSort));CHKERRQ(ierr); 904 PetscFunctionReturn(0); 905 } 906 907 #undef __FUNCT__ 908 #define __FUNCT__ "PCGASMGetSubKSP" 909 /*@C 910 PCGASMGetSubKSP - Gets the local KSP contexts for all blocks on 911 this processor. 912 913 Collective on PC iff first_local is requested 914 915 Input Parameter: 916 . pc - the preconditioner context 917 918 Output Parameters: 919 + n_local - the number of blocks on this processor or PETSC_NULL 920 . first_local - the global number of the first block on this processor or PETSC_NULL, 921 all processors must request or all must pass PETSC_NULL 922 - ksp - the array of KSP contexts 923 924 Note: 925 After PCGASMGetSubKSP() the array of KSPes is not to be freed 926 927 Currently for some matrix implementations only 1 block per processor 928 is supported. 929 930 You must call KSPSetUp() before calling PCGASMGetSubKSP(). 931 932 Level: advanced 933 934 .keywords: PC, GASM, additive Schwarz, get, sub, KSP, context 935 936 .seealso: PCGASMSetTotalSubdomains(), PCGASMSetTotalSubdomains(), PCGASMSetOverlap(), 937 PCGASMCreateSubdomains2D(), 938 @*/ 939 PetscErrorCode PETSCKSP_DLLEXPORT PCGASMGetSubKSP(PC pc,PetscInt *n_local,PetscInt *first_local,KSP *ksp[]) 940 { 941 PetscErrorCode ierr; 942 943 PetscFunctionBegin; 944 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 945 ierr = PetscUseMethod(pc,"PCGASMGetSubKSP_C",(PC,PetscInt*,PetscInt*,KSP **),(pc,n_local,first_local,ksp));CHKERRQ(ierr); 946 PetscFunctionReturn(0); 947 } 948 949 /* -------------------------------------------------------------------------------------*/ 950 /*MC 951 PCGASM - Use the (restricted) additive Schwarz method, each block is (approximately) solved with 952 its own KSP object. 953 954 Options Database Keys: 955 + -pc_gasm_truelocal - Activates PCGGASMSetUseTrueLocal() 956 . -pc_gasm_blocks <blks> - Sets total blocks 957 . -pc_gasm_overlap <ovl> - Sets overlap 958 - -pc_gasm_type [basic,restrict,interpolate,none] - Sets GASM type 959 960 IMPORTANT: If you run with, for example, 3 blocks on 1 processor or 3 blocks on 3 processors you 961 will get a different convergence rate due to the default option of -pc_gasm_type restrict. Use 962 -pc_gasm_type basic to use the standard GASM. 963 964 Notes: Each processor can have one or more blocks, but a block cannot be shared by more 965 than one processor. Defaults to one block per processor. 966 967 To set options on the solvers for each block append -sub_ to all the KSP, and PC 968 options database keys. For example, -sub_pc_type ilu -sub_pc_factor_levels 1 -sub_ksp_type preonly 969 970 To set the options on the solvers separate for each block call PCGASMGetSubKSP() 971 and set the options directly on the resulting KSP object (you can access its PC 972 with KSPGetPC()) 973 974 975 Level: beginner 976 977 Concepts: additive Schwarz method 978 979 References: 980 An additive variant of the Schwarz alternating method for the case of many subregions 981 M Dryja, OB Widlund - Courant Institute, New York University Technical report 982 983 Domain Decompositions: Parallel Multilevel Methods for Elliptic Partial Differential Equations, 984 Barry Smith, Petter Bjorstad, and William Gropp, Cambridge University Press, ISBN 0-521-49589-X. 985 986 .seealso: PCCreate(), PCSetType(), PCType (for list of available types), PC, 987 PCBJACOBI, PCGASMSetUseTrueLocal(), PCGASMGetSubKSP(), PCGASMSetLocalSubdomains(), 988 PCGASMSetTotalSubdomains(), PCSetModifySubmatrices(), PCGASMSetOverlap(), PCGASMSetType() 989 990 M*/ 991 992 EXTERN_C_BEGIN 993 #undef __FUNCT__ 994 #define __FUNCT__ "PCCreate_GASM" 995 PetscErrorCode PETSCKSP_DLLEXPORT PCCreate_GASM(PC pc) 996 { 997 PetscErrorCode ierr; 998 PC_GASM *osm; 999 1000 PetscFunctionBegin; 1001 ierr = PetscNewLog(pc,PC_GASM,&osm);CHKERRQ(ierr); 1002 osm->N = PETSC_DECIDE; 1003 osm->n = 0; 1004 osm->nmax = 0; 1005 osm->overlap = 1; 1006 osm->ksp = 0; 1007 osm->grestriction = 0; 1008 osm->gprolongation = 0; 1009 osm->gx = 0; 1010 osm->gy = 0; 1011 osm->x = 0; 1012 osm->y = 0; 1013 osm->is = 0; 1014 osm->is_local = 0; 1015 osm->pmat = 0; 1016 osm->type = PC_GASM_RESTRICT; 1017 osm->same_local_solves = PETSC_TRUE; 1018 osm->sort_indices = PETSC_TRUE; 1019 1020 pc->data = (void*)osm; 1021 pc->ops->apply = PCApply_GASM; 1022 pc->ops->applytranspose = PCApplyTranspose_GASM; 1023 pc->ops->setup = PCSetUp_GASM; 1024 pc->ops->destroy = PCDestroy_GASM; 1025 pc->ops->setfromoptions = PCSetFromOptions_GASM; 1026 pc->ops->setuponblocks = PCSetUpOnBlocks_GASM; 1027 pc->ops->view = PCView_GASM; 1028 pc->ops->applyrichardson = 0; 1029 1030 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCGASMSetLocalSubdomains_C","PCGASMSetLocalSubdomains_GASM", 1031 PCGASMSetLocalSubdomains_GASM);CHKERRQ(ierr); 1032 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCGASMSetTotalSubdomains_C","PCGASMSetTotalSubdomains_GASM", 1033 PCGASMSetTotalSubdomains_GASM);CHKERRQ(ierr); 1034 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCGASMSetOverlap_C","PCGASMSetOverlap_GASM", 1035 PCGASMSetOverlap_GASM);CHKERRQ(ierr); 1036 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCGASMSetType_C","PCGASMSetType_GASM", 1037 PCGASMSetType_GASM);CHKERRQ(ierr); 1038 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCGASMSetSortIndices_C","PCGASMSetSortIndices_GASM", 1039 PCGASMSetSortIndices_GASM);CHKERRQ(ierr); 1040 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCGASMGetSubKSP_C","PCGASMGetSubKSP_GASM", 1041 PCGASMGetSubKSP_GASM);CHKERRQ(ierr); 1042 PetscFunctionReturn(0); 1043 } 1044 EXTERN_C_END 1045 1046 1047 #undef __FUNCT__ 1048 #define __FUNCT__ "PCGASMCreateSubdomains" 1049 /*@C 1050 PCGASMCreateSubdomains - Creates the index sets for the overlapping Schwarz 1051 preconditioner for a any problem on a general grid. 1052 1053 Collective 1054 1055 Input Parameters: 1056 + A - The global matrix operator 1057 - n - the number of local blocks 1058 1059 Output Parameters: 1060 . outis - the array of index sets defining the subdomains 1061 1062 Level: advanced 1063 1064 Note: this generates nonoverlapping subdomains; the PCGASM will generate the overlap 1065 from these if you use PCGASMSetLocalSubdomains() 1066 1067 In the Fortran version you must provide the array outis[] already allocated of length n. 1068 1069 .keywords: PC, GASM, additive Schwarz, create, subdomains, unstructured grid 1070 1071 .seealso: PCGASMSetLocalSubdomains(), PCGASMDestroySubdomains() 1072 @*/ 1073 PetscErrorCode PETSCKSP_DLLEXPORT PCGASMCreateSubdomains(Mat A, PetscInt n, IS* outis[]) 1074 { 1075 MatPartitioning mpart; 1076 const char *prefix; 1077 PetscErrorCode (*f)(Mat,PetscBool *,MatReuse,Mat*); 1078 PetscMPIInt size; 1079 PetscInt i,j,rstart,rend,bs; 1080 PetscBool iscopy = PETSC_FALSE,isbaij = PETSC_FALSE,foundpart = PETSC_FALSE; 1081 Mat Ad = PETSC_NULL, adj; 1082 IS ispart,isnumb,*is; 1083 PetscErrorCode ierr; 1084 1085 PetscFunctionBegin; 1086 PetscValidHeaderSpecific(A,MAT_CLASSID,1); 1087 PetscValidPointer(outis,3); 1088 if (n < 1) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"number of local blocks must be > 0, n = %D",n); 1089 1090 /* Get prefix, row distribution, and block size */ 1091 ierr = MatGetOptionsPrefix(A,&prefix);CHKERRQ(ierr); 1092 ierr = MatGetOwnershipRange(A,&rstart,&rend);CHKERRQ(ierr); 1093 ierr = MatGetBlockSize(A,&bs);CHKERRQ(ierr); 1094 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); 1095 1096 /* Get diagonal block from matrix if possible */ 1097 ierr = MPI_Comm_size(((PetscObject)A)->comm,&size);CHKERRQ(ierr); 1098 ierr = PetscObjectQueryFunction((PetscObject)A,"MatGetDiagonalBlock_C",(void (**)(void))&f);CHKERRQ(ierr); 1099 if (f) { 1100 ierr = (*f)(A,&iscopy,MAT_INITIAL_MATRIX,&Ad);CHKERRQ(ierr); 1101 } else if (size == 1) { 1102 iscopy = PETSC_FALSE; Ad = A; 1103 } else { 1104 iscopy = PETSC_FALSE; Ad = PETSC_NULL; 1105 } 1106 if (Ad) { 1107 ierr = PetscTypeCompare((PetscObject)Ad,MATSEQBAIJ,&isbaij);CHKERRQ(ierr); 1108 if (!isbaij) {ierr = PetscTypeCompare((PetscObject)Ad,MATSEQSBAIJ,&isbaij);CHKERRQ(ierr);} 1109 } 1110 if (Ad && n > 1) { 1111 PetscBool match,done; 1112 /* Try to setup a good matrix partitioning if available */ 1113 ierr = MatPartitioningCreate(PETSC_COMM_SELF,&mpart);CHKERRQ(ierr); 1114 ierr = PetscObjectSetOptionsPrefix((PetscObject)mpart,prefix);CHKERRQ(ierr); 1115 ierr = MatPartitioningSetFromOptions(mpart);CHKERRQ(ierr); 1116 ierr = PetscTypeCompare((PetscObject)mpart,MATPARTITIONINGCURRENT,&match);CHKERRQ(ierr); 1117 if (!match) { 1118 ierr = PetscTypeCompare((PetscObject)mpart,MATPARTITIONINGSQUARE,&match);CHKERRQ(ierr); 1119 } 1120 if (!match) { /* assume a "good" partitioner is available */ 1121 PetscInt na,*ia,*ja; 1122 ierr = MatGetRowIJ(Ad,0,PETSC_TRUE,isbaij,&na,&ia,&ja,&done);CHKERRQ(ierr); 1123 if (done) { 1124 /* Build adjacency matrix by hand. Unfortunately a call to 1125 MatConvert(Ad,MATMPIADJ,MAT_INITIAL_MATRIX,&adj) will 1126 remove the block-aij structure and we cannot expect 1127 MatPartitioning to split vertices as we need */ 1128 PetscInt i,j,*row,len,nnz,cnt,*iia=0,*jja=0; 1129 nnz = 0; 1130 for (i=0; i<na; i++) { /* count number of nonzeros */ 1131 len = ia[i+1] - ia[i]; 1132 row = ja + ia[i]; 1133 for (j=0; j<len; j++) { 1134 if (row[j] == i) { /* don't count diagonal */ 1135 len--; break; 1136 } 1137 } 1138 nnz += len; 1139 } 1140 ierr = PetscMalloc((na+1)*sizeof(PetscInt),&iia);CHKERRQ(ierr); 1141 ierr = PetscMalloc((nnz)*sizeof(PetscInt),&jja);CHKERRQ(ierr); 1142 nnz = 0; 1143 iia[0] = 0; 1144 for (i=0; i<na; i++) { /* fill adjacency */ 1145 cnt = 0; 1146 len = ia[i+1] - ia[i]; 1147 row = ja + ia[i]; 1148 for (j=0; j<len; j++) { 1149 if (row[j] != i) { /* if not diagonal */ 1150 jja[nnz+cnt++] = row[j]; 1151 } 1152 } 1153 nnz += cnt; 1154 iia[i+1] = nnz; 1155 } 1156 /* Partitioning of the adjacency matrix */ 1157 ierr = MatCreateMPIAdj(PETSC_COMM_SELF,na,na,iia,jja,PETSC_NULL,&adj);CHKERRQ(ierr); 1158 ierr = MatPartitioningSetAdjacency(mpart,adj);CHKERRQ(ierr); 1159 ierr = MatPartitioningSetNParts(mpart,n);CHKERRQ(ierr); 1160 ierr = MatPartitioningApply(mpart,&ispart);CHKERRQ(ierr); 1161 ierr = ISPartitioningToNumbering(ispart,&isnumb);CHKERRQ(ierr); 1162 ierr = MatDestroy(adj);CHKERRQ(ierr); 1163 foundpart = PETSC_TRUE; 1164 } 1165 ierr = MatRestoreRowIJ(Ad,0,PETSC_TRUE,isbaij,&na,&ia,&ja,&done);CHKERRQ(ierr); 1166 } 1167 ierr = MatPartitioningDestroy(mpart);CHKERRQ(ierr); 1168 } 1169 if (iscopy) {ierr = MatDestroy(Ad);CHKERRQ(ierr);} 1170 1171 ierr = PetscMalloc(n*sizeof(IS),&is);CHKERRQ(ierr); 1172 *outis = is; 1173 1174 if (!foundpart) { 1175 1176 /* Partitioning by contiguous chunks of rows */ 1177 1178 PetscInt mbs = (rend-rstart)/bs; 1179 PetscInt start = rstart; 1180 for (i=0; i<n; i++) { 1181 PetscInt count = (mbs/n + ((mbs % n) > i)) * bs; 1182 ierr = ISCreateStride(PETSC_COMM_SELF,count,start,1,&is[i]);CHKERRQ(ierr); 1183 start += count; 1184 } 1185 1186 } else { 1187 1188 /* Partitioning by adjacency of diagonal block */ 1189 1190 const PetscInt *numbering; 1191 PetscInt *count,nidx,*indices,*newidx,start=0; 1192 /* Get node count in each partition */ 1193 ierr = PetscMalloc(n*sizeof(PetscInt),&count);CHKERRQ(ierr); 1194 ierr = ISPartitioningCount(ispart,n,count);CHKERRQ(ierr); 1195 if (isbaij && bs > 1) { /* adjust for the block-aij case */ 1196 for (i=0; i<n; i++) count[i] *= bs; 1197 } 1198 /* Build indices from node numbering */ 1199 ierr = ISGetLocalSize(isnumb,&nidx);CHKERRQ(ierr); 1200 ierr = PetscMalloc(nidx*sizeof(PetscInt),&indices);CHKERRQ(ierr); 1201 for (i=0; i<nidx; i++) indices[i] = i; /* needs to be initialized */ 1202 ierr = ISGetIndices(isnumb,&numbering);CHKERRQ(ierr); 1203 ierr = PetscSortIntWithPermutation(nidx,numbering,indices);CHKERRQ(ierr); 1204 ierr = ISRestoreIndices(isnumb,&numbering);CHKERRQ(ierr); 1205 if (isbaij && bs > 1) { /* adjust for the block-aij case */ 1206 ierr = PetscMalloc(nidx*bs*sizeof(PetscInt),&newidx);CHKERRQ(ierr); 1207 for (i=0; i<nidx; i++) 1208 for (j=0; j<bs; j++) 1209 newidx[i*bs+j] = indices[i]*bs + j; 1210 ierr = PetscFree(indices);CHKERRQ(ierr); 1211 nidx *= bs; 1212 indices = newidx; 1213 } 1214 /* Shift to get global indices */ 1215 for (i=0; i<nidx; i++) indices[i] += rstart; 1216 1217 /* Build the index sets for each block */ 1218 for (i=0; i<n; i++) { 1219 ierr = ISCreateGeneral(PETSC_COMM_SELF,count[i],&indices[start],PETSC_COPY_VALUES,&is[i]);CHKERRQ(ierr); 1220 ierr = ISSort(is[i]);CHKERRQ(ierr); 1221 start += count[i]; 1222 } 1223 1224 ierr = PetscFree(count); 1225 ierr = PetscFree(indices); 1226 ierr = ISDestroy(isnumb);CHKERRQ(ierr); 1227 ierr = ISDestroy(ispart);CHKERRQ(ierr); 1228 1229 } 1230 1231 PetscFunctionReturn(0); 1232 } 1233 1234 #undef __FUNCT__ 1235 #define __FUNCT__ "PCGASMDestroySubdomains" 1236 /*@C 1237 PCGASMDestroySubdomains - Destroys the index sets created with 1238 PCGASMCreateSubdomains(). Should be called after setting subdomains 1239 with PCGASMSetLocalSubdomains(). 1240 1241 Collective 1242 1243 Input Parameters: 1244 + n - the number of index sets 1245 . is - the array of index sets 1246 - is_local - the array of local index sets, can be PETSC_NULL 1247 1248 Level: advanced 1249 1250 .keywords: PC, GASM, additive Schwarz, create, subdomains, unstructured grid 1251 1252 .seealso: PCGASMCreateSubdomains(), PCGASMSetLocalSubdomains() 1253 @*/ 1254 PetscErrorCode PETSCKSP_DLLEXPORT PCGASMDestroySubdomains(PetscInt n, IS is[], IS is_local[]) 1255 { 1256 PetscInt i; 1257 PetscErrorCode ierr; 1258 PetscFunctionBegin; 1259 if (n <= 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"n must be > 0: n = %D",n); 1260 PetscValidPointer(is,2); 1261 for (i=0; i<n; i++) { ierr = ISDestroy(is[i]);CHKERRQ(ierr); } 1262 ierr = PetscFree(is);CHKERRQ(ierr); 1263 if (is_local) { 1264 PetscValidPointer(is_local,3); 1265 for (i=0; i<n; i++) { ierr = ISDestroy(is_local[i]);CHKERRQ(ierr); } 1266 ierr = PetscFree(is_local);CHKERRQ(ierr); 1267 } 1268 PetscFunctionReturn(0); 1269 } 1270 1271 1272 #define PCGASMLocalSubdomainBounds2D(M,N,xleft,ylow,xright,yhigh,first,last,xleft_loc,ylow_loc,xright_loc,yhigh_loc,n) \ 1273 { \ 1274 PetscInt first_row = first/M, last_row = last/M+1; \ 1275 /* \ 1276 Compute ylow_loc and yhigh_loc so that (ylow_loc,xleft) and (yhigh_loc,xright) are the corners \ 1277 of the bounding box of the intersection of the subdomain with the local ownership range (local \ 1278 subdomain). \ 1279 Also compute xleft_loc and xright_loc as the lower and upper bounds on the first and last rows \ 1280 of the intersection. \ 1281 */ \ 1282 /* ylow_loc is the grid row containing the first element of the local sumbdomain */ \ 1283 ylow_loc = PetscMax(first_row,ylow); \ 1284 /* xleft_loc is the offset of first element of the local subdomain within its grid row (might actually be outside the local subdomain) */ \ 1285 xleft_loc = ylow_loc==first_row?PetscMax(first%M,xleft):xleft; \ 1286 /* yhigh_loc is the grid row above the last local subdomain element */ \ 1287 yhigh_loc = PetscMin(last_row,yhigh); \ 1288 /* xright is the offset of the end of the local subdomain within its grid row (might actually be outside the local subdomain) */ \ 1289 xright_loc = yhigh_loc==last_row?PetscMin(xright,last%M):xright; \ 1290 /* Now compute the size of the local subdomain n. */ \ 1291 n = 0; \ 1292 if(ylow_loc < yhigh_loc) { \ 1293 PetscInt width = xright-xleft; \ 1294 n += width*(yhigh_loc-ylow_loc-1); \ 1295 n += PetscMin(PetscMax(xright_loc-xleft,0),width); \ 1296 n -= PetscMin(PetscMax(xleft_loc-xleft,0), width); \ 1297 }\ 1298 } 1299 1300 #undef __FUNCT__ 1301 #define __FUNCT__ "PCGASMCreateSubdomains2D" 1302 /*@ 1303 PCGASMCreateSubdomains2D - Creates the index sets for the overlapping Schwarz 1304 preconditioner for a two-dimensional problem on a regular grid. 1305 1306 Collective 1307 1308 Input Parameters: 1309 + M, N - the global number of mesh points in the x and y directions 1310 . Mdomains, Ndomains - the global number of subdomains in the x and y directions 1311 . dof - degrees of freedom per node 1312 - overlap - overlap in mesh lines 1313 1314 Output Parameters: 1315 + Nsub - the number of local subdomains created 1316 . is - array of index sets defining overlapping (if overlap > 0) subdomains 1317 - is_local - array of index sets defining non-overlapping subdomains 1318 1319 Note: 1320 Presently PCAMSCreateSubdomains2d() is valid only for sequential 1321 preconditioners. More general related routines are 1322 PCGASMSetTotalSubdomains() and PCGASMSetLocalSubdomains(). 1323 1324 Level: advanced 1325 1326 .keywords: PC, GASM, additive Schwarz, create, subdomains, 2D, regular grid 1327 1328 .seealso: PCGASMSetTotalSubdomains(), PCGASMSetLocalSubdomains(), PCGASMGetSubKSP(), 1329 PCGASMSetOverlap() 1330 @*/ 1331 PetscErrorCode PETSCKSP_DLLEXPORT PCGASMCreateSubdomains2D(PC pc, PetscInt M,PetscInt N,PetscInt Mdomains,PetscInt Ndomains,PetscInt dof,PetscInt overlap, PetscInt *nsub,IS **is,IS **is_local) 1332 { 1333 PetscErrorCode ierr; 1334 PetscMPIInt size, rank; 1335 PetscInt i, j; 1336 PetscInt maxheight, maxwidth; 1337 PetscInt xstart, xleft, xright, xleft_loc, xright_loc; 1338 PetscInt ystart, ylow, yhigh, ylow_loc, yhigh_loc; 1339 PetscInt x[2][2], y[2][2]; 1340 PetscInt first, last; 1341 PetscInt nidx, *idx; 1342 PetscInt ii,jj,s,q,d; 1343 PetscInt k,kk; 1344 PetscMPIInt color; 1345 MPI_Comm comm, subcomm; 1346 1347 IS **iis; 1348 PetscFunctionBegin; 1349 ierr = PetscObjectGetComm((PetscObject)pc, &comm); CHKERRQ(ierr); 1350 ierr = MPI_Comm_size(comm, &size); CHKERRQ(ierr); 1351 ierr = MPI_Comm_rank(comm, &rank); CHKERRQ(ierr); 1352 ierr = MatGetOwnershipRange(pc->pmat, &first, &last); CHKERRQ(ierr); 1353 if(first%dof || last%dof) { 1354 SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Matrix row partitioning unsuitable for domain decomposition: local row range (%D,%D) " 1355 "does not respect the number of degrees of freedom per grid point %D", first, last, dof); 1356 } 1357 /* Determine the number of domains with nonzero intersections with the local ownership range. */ 1358 s = 0; 1359 ystart = 0; 1360 for (j=0; j<Ndomains; ++j) { 1361 maxheight = N/Ndomains + ((N % Ndomains) > j); /* Maximal height of subdomain */ 1362 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); 1363 /* Vertical domain limits without an overlap. */ 1364 ylow = ystart; 1365 yhigh = PetscMin(ystart + maxheight,N); 1366 xstart = 0; 1367 for (i=0; i<Mdomains; ++i) { 1368 maxwidth = M/Mdomains + ((M % Mdomains) > i); /* Maximal width of subdomain */ 1369 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); 1370 /* Horizontal domain limits without an overlap. */ 1371 xleft = xstart; 1372 xright = PetscMin(xstart+maxwidth,M); 1373 /* 1374 Determine whether this subdomain intersects this processor's ownership range of pc->pmat. 1375 */ 1376 PCGASMLocalSubdomainBounds2D(M,N,xleft,ylow,xright,yhigh,first,last,xleft_loc,ylow_loc,xright_loc,yhigh_loc,nidx); 1377 if(nidx) { 1378 ++s; 1379 } 1380 xstart += maxwidth; 1381 }/* for(i = 0; i < Mdomains; ++i) */ 1382 ystart += maxheight; 1383 }/* for(j = 0; j < Ndomains; ++j) */ 1384 /* Now we can allocate the necessary number of ISs. */ 1385 *nsub = s; 1386 ierr = PetscMalloc((*nsub)*sizeof(IS*),is); CHKERRQ(ierr); 1387 ierr = PetscMalloc((*nsub)*sizeof(IS*),is_local); CHKERRQ(ierr); 1388 s = 0; 1389 ystart = 0; 1390 for (j=0; j<Ndomains; ++j) { 1391 maxheight = N/Ndomains + ((N % Ndomains) > j); /* Maximal height of subdomain */ 1392 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); 1393 /* Vertical domain limits without an overlap. */ 1394 y[0][0] = ystart; 1395 y[0][1] = PetscMin(ystart + maxheight,N); 1396 /* Vertical domain limits with an overlap. */ 1397 y[1][1] = PetscMin(ystart + maxheight + overlap,N); 1398 y[1][0] = PetscMax(ystart - overlap,0); xstart = 0; 1399 for (i=0; i<Mdomains; ++i) { 1400 maxwidth = M/Mdomains + ((M % Mdomains) > i); /* Maximal width of subdomain */ 1401 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); 1402 /* Horizontal domain limits without an overlap. */ 1403 x[0][0] = xstart; 1404 x[0][1] = PetscMin(xstart+maxwidth,M); 1405 /* Horizontal domain limits with an overlap. */ 1406 x[1][0] = PetscMax(xstart - overlap,0); 1407 x[1][1] = PetscMin(xstart + maxwidth + overlap,M); 1408 /* 1409 Determine whether this domain intersects this processor's ownership range of pc->pmat. 1410 Do this twice: first for the domains with overlaps, and once without. 1411 During the first pass create the subcommunicators, and use them on the second pass as well. 1412 */ 1413 for(q = 0; q < 2; ++q) { 1414 /* 1415 domain limits, (xleft, xright) and (ylow, yheigh) are adjusted 1416 according to whether the domain with an overlap or without is considered. 1417 */ 1418 xleft = x[q][0]; xright = x[q][1]; 1419 ylow = y[q][0]; yhigh = y[q][1]; 1420 PCGASMLocalSubdomainBounds2D(M,N,xleft,ylow,xright,yhigh,first,last,xleft_loc,ylow_loc,xright_loc,yhigh_loc,nidx); 1421 nidx *= dof; 1422 /* 1423 Based on the counted number of indices in the local domain (without an overlap), 1424 construct a subcommunicator of all the processors supporting this domain. 1425 */ 1426 if(q == 0) { 1427 if(nidx) { 1428 color = 1; 1429 } 1430 else { 1431 color = MPI_UNDEFINED; 1432 } 1433 ierr = MPI_Comm_split(comm, color, rank, &subcomm); CHKERRQ(ierr); 1434 } 1435 /* 1436 Proceed only if the number of local indices is nonzero. 1437 */ 1438 if(nidx) { 1439 if(q == 0) { 1440 iis = is_local; 1441 } 1442 if(q == 1) { 1443 /* 1444 The IS for the no-overlap subdomain shares a communicator with the nonoverlapping domain. 1445 Moreover, if the overlap is zero, the two ISs are identical. 1446 */ 1447 if(overlap == 0) { 1448 (*is)[s] = (*is_local)[s]; 1449 ierr = PetscObjectReference((PetscObject)(*is)[s]);CHKERRQ(ierr); 1450 continue; 1451 } 1452 else { 1453 iis = is; 1454 subcomm = ((PetscObject)(*is_local)[s])->comm; 1455 } 1456 }/* if(q == 1) */ 1457 ierr = PetscMalloc(nidx*sizeof(PetscInt),&idx);CHKERRQ(ierr); 1458 k = 0; 1459 for (jj=ylow_loc; jj<yhigh_loc; ++jj) { 1460 PetscInt x0 = (jj==ylow_loc)?xleft_loc:xleft; 1461 PetscInt x1 = (jj==yhigh_loc-1)?xright_loc:xright; 1462 kk = dof*(M*jj + x0); 1463 for (ii=x0; ii<x1; ++ii) { 1464 for(d = 0; d < dof; ++d) { 1465 idx[k++] = kk++; 1466 } 1467 } 1468 } 1469 ierr = ISCreateGeneral(subcomm,nidx,idx,PETSC_OWN_POINTER,(*iis)+s);CHKERRQ(ierr); 1470 }/* if(nidx) */ 1471 }/* for(q = 0; q < 2; ++q) */ 1472 if(nidx) { 1473 ++s; 1474 } 1475 xstart += maxwidth; 1476 }/* for(i = 0; i < Mdomains; ++i) */ 1477 ystart += maxheight; 1478 }/* for(j = 0; j < Ndomains; ++j) */ 1479 PetscFunctionReturn(0); 1480 } 1481 1482 #undef __FUNCT__ 1483 #define __FUNCT__ "PCGASMGetLocalSubdomains" 1484 /*@C 1485 PCGASMGetLocalSubdomains - Gets the local subdomains (for this processor 1486 only) for the additive Schwarz preconditioner. 1487 1488 Not Collective 1489 1490 Input Parameter: 1491 . pc - the preconditioner context 1492 1493 Output Parameters: 1494 + n - the number of subdomains for this processor (default value = 1) 1495 . is - the index sets that define the subdomains for this processor 1496 - is_local - the index sets that define the local part of the subdomains for this processor (can be PETSC_NULL) 1497 1498 1499 Notes: 1500 The IS numbering is in the parallel, global numbering of the vector. 1501 1502 Level: advanced 1503 1504 .keywords: PC, GASM, set, local, subdomains, additive Schwarz 1505 1506 .seealso: PCGASMSetTotalSubdomains(), PCGASMSetOverlap(), PCGASMGetSubKSP(), 1507 PCGASMCreateSubdomains2D(), PCGASMSetLocalSubdomains(), PCGASMGetLocalSubmatrices() 1508 @*/ 1509 PetscErrorCode PETSCKSP_DLLEXPORT PCGASMGetLocalSubdomains(PC pc,PetscInt *n,IS *is[],IS *is_local[]) 1510 { 1511 PC_GASM *osm; 1512 PetscErrorCode ierr; 1513 PetscBool match; 1514 1515 PetscFunctionBegin; 1516 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1517 PetscValidIntPointer(n,2); 1518 if (is) PetscValidPointer(is,3); 1519 ierr = PetscTypeCompare((PetscObject)pc,PCGASM,&match);CHKERRQ(ierr); 1520 if (!match) { 1521 if (n) *n = 0; 1522 if (is) *is = PETSC_NULL; 1523 } else { 1524 osm = (PC_GASM*)pc->data; 1525 if (n) *n = osm->n; 1526 if (is) *is = osm->is; 1527 if (is_local) *is_local = osm->is_local; 1528 } 1529 PetscFunctionReturn(0); 1530 } 1531 1532 #undef __FUNCT__ 1533 #define __FUNCT__ "PCGASMGetLocalSubmatrices" 1534 /*@C 1535 PCGASMGetLocalSubmatrices - Gets the local submatrices (for this processor 1536 only) for the additive Schwarz preconditioner. 1537 1538 Not Collective 1539 1540 Input Parameter: 1541 . pc - the preconditioner context 1542 1543 Output Parameters: 1544 + n - the number of matrices for this processor (default value = 1) 1545 - mat - the matrices 1546 1547 1548 Level: advanced 1549 1550 .keywords: PC, GASM, set, local, subdomains, additive Schwarz, block Jacobi 1551 1552 .seealso: PCGASMSetTotalSubdomains(), PCGASMSetOverlap(), PCGASMGetSubKSP(), 1553 PCGASMCreateSubdomains2D(), PCGASMSetLocalSubdomains(), PCGASMGetLocalSubdomains() 1554 @*/ 1555 PetscErrorCode PETSCKSP_DLLEXPORT PCGASMGetLocalSubmatrices(PC pc,PetscInt *n,Mat *mat[]) 1556 { 1557 PC_GASM *osm; 1558 PetscErrorCode ierr; 1559 PetscBool match; 1560 1561 PetscFunctionBegin; 1562 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1563 PetscValidIntPointer(n,2); 1564 if (mat) PetscValidPointer(mat,3); 1565 if (!pc->setupcalled) SETERRQ(((PetscObject)pc)->comm,PETSC_ERR_ARG_WRONGSTATE,"Must call after KSPSetUP() or PCSetUp()."); 1566 ierr = PetscTypeCompare((PetscObject)pc,PCGASM,&match);CHKERRQ(ierr); 1567 if (!match) { 1568 if (n) *n = 0; 1569 if (mat) *mat = PETSC_NULL; 1570 } else { 1571 osm = (PC_GASM*)pc->data; 1572 if (n) *n = osm->n; 1573 if (mat) *mat = osm->pmat; 1574 } 1575 PetscFunctionReturn(0); 1576 } 1577