1 /*
2 This file defines an additive Schwarz preconditioner for any Mat implementation.
3
4 Note that each processor may have any number of subdomains. But in order to
5 deal easily with the VecScatter(), we treat each processor as if it has the
6 same number of subdomains.
7
8 n - total number of true subdomains on all processors
9 n_local_true - actual number of subdomains on this processor
10 n_local = maximum over all processors of n_local_true
11 */
12
13 #include <../src/ksp/pc/impls/asm/asm.h> /*I "petscpc.h" I*/
14
PCView_ASM(PC pc,PetscViewer viewer)15 static PetscErrorCode PCView_ASM(PC pc,PetscViewer viewer)
16 {
17 PC_ASM *osm = (PC_ASM*)pc->data;
18 PetscErrorCode ierr;
19 PetscMPIInt rank;
20 PetscInt i,bsz;
21 PetscBool iascii,isstring;
22 PetscViewer sviewer;
23
24 PetscFunctionBegin;
25 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr);
26 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERSTRING,&isstring);CHKERRQ(ierr);
27 if (iascii) {
28 char overlaps[256] = "user-defined overlap",blocks[256] = "total subdomain blocks not yet set";
29 if (osm->overlap >= 0) {ierr = PetscSNPrintf(overlaps,sizeof(overlaps),"amount of overlap = %D",osm->overlap);CHKERRQ(ierr);}
30 if (osm->n > 0) {ierr = PetscSNPrintf(blocks,sizeof(blocks),"total subdomain blocks = %D",osm->n);CHKERRQ(ierr);}
31 ierr = PetscViewerASCIIPrintf(viewer," %s, %s\n",blocks,overlaps);CHKERRQ(ierr);
32 ierr = PetscViewerASCIIPrintf(viewer," restriction/interpolation type - %s\n",PCASMTypes[osm->type]);CHKERRQ(ierr);
33 if (osm->dm_subdomains) {ierr = PetscViewerASCIIPrintf(viewer," Additive Schwarz: using DM to define subdomains\n");CHKERRQ(ierr);}
34 if (osm->loctype != PC_COMPOSITE_ADDITIVE) {ierr = PetscViewerASCIIPrintf(viewer," Additive Schwarz: local solve composition type - %s\n",PCCompositeTypes[osm->loctype]);CHKERRQ(ierr);}
35 ierr = MPI_Comm_rank(PetscObjectComm((PetscObject)pc),&rank);CHKERRQ(ierr);
36 if (osm->same_local_solves) {
37 if (osm->ksp) {
38 ierr = PetscViewerASCIIPrintf(viewer," Local solver is the same for all blocks, as in the following KSP and PC objects on rank 0:\n");CHKERRQ(ierr);
39 ierr = PetscViewerGetSubViewer(viewer,PETSC_COMM_SELF,&sviewer);CHKERRQ(ierr);
40 if (!rank) {
41 ierr = PetscViewerASCIIPushTab(viewer);CHKERRQ(ierr);
42 ierr = KSPView(osm->ksp[0],sviewer);CHKERRQ(ierr);
43 ierr = PetscViewerASCIIPopTab(viewer);CHKERRQ(ierr);
44 }
45 ierr = PetscViewerRestoreSubViewer(viewer,PETSC_COMM_SELF,&sviewer);CHKERRQ(ierr);
46 }
47 } else {
48 ierr = PetscViewerASCIIPushSynchronized(viewer);CHKERRQ(ierr);
49 ierr = PetscViewerASCIISynchronizedPrintf(viewer," [%d] number of local blocks = %D\n",(int)rank,osm->n_local_true);CHKERRQ(ierr);
50 ierr = PetscViewerFlush(viewer);CHKERRQ(ierr);
51 ierr = PetscViewerASCIIPrintf(viewer," Local solve info for each block is in the following KSP and PC objects:\n");CHKERRQ(ierr);
52 ierr = PetscViewerASCIIPushTab(viewer);CHKERRQ(ierr);
53 ierr = PetscViewerASCIIPrintf(viewer,"- - - - - - - - - - - - - - - - - -\n");CHKERRQ(ierr);
54 ierr = PetscViewerGetSubViewer(viewer,PETSC_COMM_SELF,&sviewer);CHKERRQ(ierr);
55 for (i=0; i<osm->n_local_true; i++) {
56 ierr = ISGetLocalSize(osm->is[i],&bsz);CHKERRQ(ierr);
57 ierr = PetscViewerASCIISynchronizedPrintf(sviewer,"[%d] local block number %D, size = %D\n",(int)rank,i,bsz);CHKERRQ(ierr);
58 ierr = KSPView(osm->ksp[i],sviewer);CHKERRQ(ierr);
59 ierr = PetscViewerASCIISynchronizedPrintf(sviewer,"- - - - - - - - - - - - - - - - - -\n");CHKERRQ(ierr);
60 }
61 ierr = PetscViewerRestoreSubViewer(viewer,PETSC_COMM_SELF,&sviewer);CHKERRQ(ierr);
62 ierr = PetscViewerASCIIPopTab(viewer);CHKERRQ(ierr);
63 ierr = PetscViewerFlush(viewer);CHKERRQ(ierr);
64 ierr = PetscViewerASCIIPopSynchronized(viewer);CHKERRQ(ierr);
65 }
66 } else if (isstring) {
67 ierr = PetscViewerStringSPrintf(viewer," blocks=%D, overlap=%D, type=%s",osm->n,osm->overlap,PCASMTypes[osm->type]);CHKERRQ(ierr);
68 ierr = PetscViewerGetSubViewer(viewer,PETSC_COMM_SELF,&sviewer);CHKERRQ(ierr);
69 if (osm->ksp) {ierr = KSPView(osm->ksp[0],sviewer);CHKERRQ(ierr);}
70 ierr = PetscViewerRestoreSubViewer(viewer,PETSC_COMM_SELF,&sviewer);CHKERRQ(ierr);
71 }
72 PetscFunctionReturn(0);
73 }
74
PCASMPrintSubdomains(PC pc)75 static PetscErrorCode PCASMPrintSubdomains(PC pc)
76 {
77 PC_ASM *osm = (PC_ASM*)pc->data;
78 const char *prefix;
79 char fname[PETSC_MAX_PATH_LEN+1];
80 PetscViewer viewer, sviewer;
81 char *s;
82 PetscInt i,j,nidx;
83 const PetscInt *idx;
84 PetscMPIInt rank, size;
85 PetscErrorCode ierr;
86
87 PetscFunctionBegin;
88 ierr = MPI_Comm_size(PetscObjectComm((PetscObject)pc), &size);CHKERRQ(ierr);
89 ierr = MPI_Comm_rank(PetscObjectComm((PetscObject)pc), &rank);CHKERRQ(ierr);
90 ierr = PCGetOptionsPrefix(pc,&prefix);CHKERRQ(ierr);
91 ierr = PetscOptionsGetString(NULL,prefix,"-pc_asm_print_subdomains",fname,sizeof(fname),NULL);CHKERRQ(ierr);
92 if (fname[0] == 0) { ierr = PetscStrcpy(fname,"stdout");CHKERRQ(ierr); };
93 ierr = PetscViewerASCIIOpen(PetscObjectComm((PetscObject)pc),fname,&viewer);CHKERRQ(ierr);
94 for (i=0; i<osm->n_local; i++) {
95 if (i < osm->n_local_true) {
96 ierr = ISGetLocalSize(osm->is[i],&nidx);CHKERRQ(ierr);
97 ierr = ISGetIndices(osm->is[i],&idx);CHKERRQ(ierr);
98 /* Print to a string viewer; no more than 15 characters per index plus 512 char for the header.*/
99 #define len 16*(nidx+1)+512
100 ierr = PetscMalloc1(len,&s);CHKERRQ(ierr);
101 ierr = PetscViewerStringOpen(PETSC_COMM_SELF, s, len, &sviewer);CHKERRQ(ierr);
102 #undef len
103 ierr = PetscViewerStringSPrintf(sviewer, "[%D:%D] Subdomain %D with overlap:\n", rank, size, i);CHKERRQ(ierr);
104 for (j=0; j<nidx; j++) {
105 ierr = PetscViewerStringSPrintf(sviewer,"%D ",idx[j]);CHKERRQ(ierr);
106 }
107 ierr = ISRestoreIndices(osm->is[i],&idx);CHKERRQ(ierr);
108 ierr = PetscViewerStringSPrintf(sviewer,"\n");CHKERRQ(ierr);
109 ierr = PetscViewerDestroy(&sviewer);CHKERRQ(ierr);
110 ierr = PetscViewerASCIIPushSynchronized(viewer);CHKERRQ(ierr);
111 ierr = PetscViewerASCIISynchronizedPrintf(viewer, s);CHKERRQ(ierr);
112 ierr = PetscViewerFlush(viewer);CHKERRQ(ierr);
113 ierr = PetscViewerASCIIPopSynchronized(viewer);CHKERRQ(ierr);
114 ierr = PetscFree(s);CHKERRQ(ierr);
115 if (osm->is_local) {
116 /* Print to a string viewer; no more than 15 characters per index plus 512 char for the header.*/
117 #define len 16*(nidx+1)+512
118 ierr = PetscMalloc1(len, &s);CHKERRQ(ierr);
119 ierr = PetscViewerStringOpen(PETSC_COMM_SELF, s, len, &sviewer);CHKERRQ(ierr);
120 #undef len
121 ierr = PetscViewerStringSPrintf(sviewer, "[%D:%D] Subdomain %D without overlap:\n", rank, size, i);CHKERRQ(ierr);
122 ierr = ISGetLocalSize(osm->is_local[i],&nidx);CHKERRQ(ierr);
123 ierr = ISGetIndices(osm->is_local[i],&idx);CHKERRQ(ierr);
124 for (j=0; j<nidx; j++) {
125 ierr = PetscViewerStringSPrintf(sviewer,"%D ",idx[j]);CHKERRQ(ierr);
126 }
127 ierr = ISRestoreIndices(osm->is_local[i],&idx);CHKERRQ(ierr);
128 ierr = PetscViewerStringSPrintf(sviewer,"\n");CHKERRQ(ierr);
129 ierr = PetscViewerDestroy(&sviewer);CHKERRQ(ierr);
130 ierr = PetscViewerASCIIPushSynchronized(viewer);CHKERRQ(ierr);
131 ierr = PetscViewerASCIISynchronizedPrintf(viewer, s);CHKERRQ(ierr);
132 ierr = PetscViewerFlush(viewer);CHKERRQ(ierr);
133 ierr = PetscViewerASCIIPopSynchronized(viewer);CHKERRQ(ierr);
134 ierr = PetscFree(s);CHKERRQ(ierr);
135 }
136 } else {
137 /* Participate in collective viewer calls. */
138 ierr = PetscViewerASCIIPushSynchronized(viewer);CHKERRQ(ierr);
139 ierr = PetscViewerFlush(viewer);CHKERRQ(ierr);
140 ierr = PetscViewerASCIIPopSynchronized(viewer);CHKERRQ(ierr);
141 /* Assume either all ranks have is_local or none do. */
142 if (osm->is_local) {
143 ierr = PetscViewerASCIIPushSynchronized(viewer);CHKERRQ(ierr);
144 ierr = PetscViewerFlush(viewer);CHKERRQ(ierr);
145 ierr = PetscViewerASCIIPopSynchronized(viewer);CHKERRQ(ierr);
146 }
147 }
148 }
149 ierr = PetscViewerFlush(viewer);CHKERRQ(ierr);
150 ierr = PetscViewerDestroy(&viewer);CHKERRQ(ierr);
151 PetscFunctionReturn(0);
152 }
153
PCSetUp_ASM(PC pc)154 static PetscErrorCode PCSetUp_ASM(PC pc)
155 {
156 PC_ASM *osm = (PC_ASM*)pc->data;
157 PetscErrorCode ierr;
158 PetscBool flg;
159 PetscInt i,m,m_local;
160 MatReuse scall = MAT_REUSE_MATRIX;
161 IS isl;
162 KSP ksp;
163 PC subpc;
164 const char *prefix,*pprefix;
165 Vec vec;
166 DM *domain_dm = NULL;
167
168 PetscFunctionBegin;
169 if (!pc->setupcalled) {
170 PetscInt m;
171
172 /* Note: if subdomains have been set either via PCASMSetTotalSubdomains() or via PCASMSetLocalSubdomains(), osm->n_local_true will not be PETSC_DECIDE */
173 if (osm->n_local_true == PETSC_DECIDE) {
174 /* no subdomains given */
175 /* try pc->dm first, if allowed */
176 if (osm->dm_subdomains && pc->dm) {
177 PetscInt num_domains, d;
178 char **domain_names;
179 IS *inner_domain_is, *outer_domain_is;
180 ierr = DMCreateDomainDecomposition(pc->dm, &num_domains, &domain_names, &inner_domain_is, &outer_domain_is, &domain_dm);CHKERRQ(ierr);
181 osm->overlap = -1; /* We do not want to increase the overlap of the IS.
182 A future improvement of this code might allow one to use
183 DM-defined subdomains and also increase the overlap,
184 but that is not currently supported */
185 if (num_domains) {
186 ierr = PCASMSetLocalSubdomains(pc, num_domains, outer_domain_is, inner_domain_is);CHKERRQ(ierr);
187 }
188 for (d = 0; d < num_domains; ++d) {
189 if (domain_names) {ierr = PetscFree(domain_names[d]);CHKERRQ(ierr);}
190 if (inner_domain_is) {ierr = ISDestroy(&inner_domain_is[d]);CHKERRQ(ierr);}
191 if (outer_domain_is) {ierr = ISDestroy(&outer_domain_is[d]);CHKERRQ(ierr);}
192 }
193 ierr = PetscFree(domain_names);CHKERRQ(ierr);
194 ierr = PetscFree(inner_domain_is);CHKERRQ(ierr);
195 ierr = PetscFree(outer_domain_is);CHKERRQ(ierr);
196 }
197 if (osm->n_local_true == PETSC_DECIDE) {
198 /* still no subdomains; use one subdomain per processor */
199 osm->n_local_true = 1;
200 }
201 }
202 { /* determine the global and max number of subdomains */
203 struct {PetscInt max,sum;} inwork,outwork;
204 PetscMPIInt size;
205
206 inwork.max = osm->n_local_true;
207 inwork.sum = osm->n_local_true;
208 ierr = MPIU_Allreduce(&inwork,&outwork,1,MPIU_2INT,MPIU_MAXSUM_OP,PetscObjectComm((PetscObject)pc));CHKERRQ(ierr);
209 osm->n_local = outwork.max;
210 osm->n = outwork.sum;
211
212 ierr = MPI_Comm_size(PetscObjectComm((PetscObject)pc),&size);CHKERRQ(ierr);
213 if (outwork.max == 1 && outwork.sum == size) {
214 /* osm->n_local_true = 1 on all processes, set this option may enable use of optimized MatCreateSubMatrices() implementation */
215 ierr = MatSetOption(pc->pmat,MAT_SUBMAT_SINGLEIS,PETSC_TRUE);CHKERRQ(ierr);
216 }
217 }
218 if (!osm->is) { /* create the index sets */
219 ierr = PCASMCreateSubdomains(pc->pmat,osm->n_local_true,&osm->is);CHKERRQ(ierr);
220 }
221 if (osm->n_local_true > 1 && !osm->is_local) {
222 ierr = PetscMalloc1(osm->n_local_true,&osm->is_local);CHKERRQ(ierr);
223 for (i=0; i<osm->n_local_true; i++) {
224 if (osm->overlap > 0) { /* With positive overlap, osm->is[i] will be modified */
225 ierr = ISDuplicate(osm->is[i],&osm->is_local[i]);CHKERRQ(ierr);
226 ierr = ISCopy(osm->is[i],osm->is_local[i]);CHKERRQ(ierr);
227 } else {
228 ierr = PetscObjectReference((PetscObject)osm->is[i]);CHKERRQ(ierr);
229 osm->is_local[i] = osm->is[i];
230 }
231 }
232 }
233 ierr = PCGetOptionsPrefix(pc,&prefix);CHKERRQ(ierr);
234 flg = PETSC_FALSE;
235 ierr = PetscOptionsGetBool(NULL,prefix,"-pc_asm_print_subdomains",&flg,NULL);CHKERRQ(ierr);
236 if (flg) { ierr = PCASMPrintSubdomains(pc);CHKERRQ(ierr); }
237
238 if (osm->overlap > 0) {
239 /* Extend the "overlapping" regions by a number of steps */
240 ierr = MatIncreaseOverlap(pc->pmat,osm->n_local_true,osm->is,osm->overlap);CHKERRQ(ierr);
241 }
242 if (osm->sort_indices) {
243 for (i=0; i<osm->n_local_true; i++) {
244 ierr = ISSort(osm->is[i]);CHKERRQ(ierr);
245 if (osm->is_local) {
246 ierr = ISSort(osm->is_local[i]);CHKERRQ(ierr);
247 }
248 }
249 }
250
251 if (!osm->ksp) {
252 /* Create the local solvers */
253 ierr = PetscMalloc1(osm->n_local_true,&osm->ksp);CHKERRQ(ierr);
254 if (domain_dm) {
255 ierr = PetscInfo(pc,"Setting up ASM subproblems using the embedded DM\n");CHKERRQ(ierr);
256 }
257 for (i=0; i<osm->n_local_true; i++) {
258 ierr = KSPCreate(PETSC_COMM_SELF,&ksp);CHKERRQ(ierr);
259 ierr = KSPSetErrorIfNotConverged(ksp,pc->erroriffailure);CHKERRQ(ierr);
260 ierr = PetscLogObjectParent((PetscObject)pc,(PetscObject)ksp);CHKERRQ(ierr);
261 ierr = PetscObjectIncrementTabLevel((PetscObject)ksp,(PetscObject)pc,1);CHKERRQ(ierr);
262 ierr = KSPSetType(ksp,KSPPREONLY);CHKERRQ(ierr);
263 ierr = KSPGetPC(ksp,&subpc);CHKERRQ(ierr);
264 ierr = PCGetOptionsPrefix(pc,&prefix);CHKERRQ(ierr);
265 ierr = KSPSetOptionsPrefix(ksp,prefix);CHKERRQ(ierr);
266 ierr = KSPAppendOptionsPrefix(ksp,"sub_");CHKERRQ(ierr);
267 if (domain_dm) {
268 ierr = KSPSetDM(ksp, domain_dm[i]);CHKERRQ(ierr);
269 ierr = KSPSetDMActive(ksp, PETSC_FALSE);CHKERRQ(ierr);
270 ierr = DMDestroy(&domain_dm[i]);CHKERRQ(ierr);
271 }
272 osm->ksp[i] = ksp;
273 }
274 if (domain_dm) {
275 ierr = PetscFree(domain_dm);CHKERRQ(ierr);
276 }
277 }
278
279 ierr = ISConcatenate(PETSC_COMM_SELF, osm->n_local_true, osm->is, &osm->lis);CHKERRQ(ierr);
280 ierr = ISSortRemoveDups(osm->lis);CHKERRQ(ierr);
281 ierr = ISGetLocalSize(osm->lis, &m);CHKERRQ(ierr);
282
283 scall = MAT_INITIAL_MATRIX;
284 } else {
285 /*
286 Destroy the blocks from the previous iteration
287 */
288 if (pc->flag == DIFFERENT_NONZERO_PATTERN) {
289 ierr = MatDestroyMatrices(osm->n_local_true,&osm->pmat);CHKERRQ(ierr);
290 scall = MAT_INITIAL_MATRIX;
291 }
292 }
293
294 /*
295 Extract out the submatrices
296 */
297 ierr = MatCreateSubMatrices(pc->pmat,osm->n_local_true,osm->is,osm->is,scall,&osm->pmat);CHKERRQ(ierr);
298 if (scall == MAT_INITIAL_MATRIX) {
299 ierr = PetscObjectGetOptionsPrefix((PetscObject)pc->pmat,&pprefix);CHKERRQ(ierr);
300 for (i=0; i<osm->n_local_true; i++) {
301 ierr = PetscLogObjectParent((PetscObject)pc,(PetscObject)osm->pmat[i]);CHKERRQ(ierr);
302 ierr = PetscObjectSetOptionsPrefix((PetscObject)osm->pmat[i],pprefix);CHKERRQ(ierr);
303 }
304 }
305
306 /* Convert the types of the submatrices (if needbe) */
307 if (osm->sub_mat_type) {
308 for (i=0; i<osm->n_local_true; i++) {
309 ierr = MatConvert(osm->pmat[i],osm->sub_mat_type,MAT_INPLACE_MATRIX,&(osm->pmat[i]));CHKERRQ(ierr);
310 }
311 }
312
313 if (!pc->setupcalled) {
314 VecType vtype;
315
316 /* Create the local work vectors (from the local matrices) and scatter contexts */
317 ierr = MatCreateVecs(pc->pmat,&vec,NULL);CHKERRQ(ierr);
318
319 if (osm->is_local && (osm->type == PC_ASM_INTERPOLATE || osm->type == PC_ASM_NONE)) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP,"Cannot use interpolate or none PCASMType if is_local was provided to PCASMSetLocalSubdomains()");
320 if (osm->is_local && osm->type == PC_ASM_RESTRICT && osm->loctype == PC_COMPOSITE_ADDITIVE) {
321 ierr = PetscMalloc1(osm->n_local_true,&osm->lprolongation);CHKERRQ(ierr);
322 }
323 ierr = PetscMalloc1(osm->n_local_true,&osm->lrestriction);CHKERRQ(ierr);
324 ierr = PetscMalloc1(osm->n_local_true,&osm->x);CHKERRQ(ierr);
325 ierr = PetscMalloc1(osm->n_local_true,&osm->y);CHKERRQ(ierr);
326
327 ierr = ISGetLocalSize(osm->lis,&m);CHKERRQ(ierr);
328 ierr = ISCreateStride(PETSC_COMM_SELF,m,0,1,&isl);CHKERRQ(ierr);
329 ierr = MatGetVecType(osm->pmat[0],&vtype);CHKERRQ(ierr);
330 ierr = VecCreate(PETSC_COMM_SELF,&osm->lx);CHKERRQ(ierr);
331 ierr = VecSetSizes(osm->lx,m,m);CHKERRQ(ierr);
332 ierr = VecSetType(osm->lx,vtype);CHKERRQ(ierr);
333 ierr = VecDuplicate(osm->lx, &osm->ly);CHKERRQ(ierr);
334 ierr = VecScatterCreate(vec,osm->lis,osm->lx,isl,&osm->restriction);CHKERRQ(ierr);
335 ierr = ISDestroy(&isl);CHKERRQ(ierr);
336
337 for (i=0; i<osm->n_local_true; ++i) {
338 ISLocalToGlobalMapping ltog;
339 IS isll;
340 const PetscInt *idx_is;
341 PetscInt *idx_lis,nout;
342
343 ierr = ISGetLocalSize(osm->is[i],&m);CHKERRQ(ierr);
344 ierr = MatCreateVecs(osm->pmat[i],&osm->x[i],NULL);CHKERRQ(ierr);
345 ierr = VecDuplicate(osm->x[i],&osm->y[i]);CHKERRQ(ierr);
346
347 /* generate a scatter from ly to y[i] picking all the overlapping is[i] entries */
348 ierr = ISLocalToGlobalMappingCreateIS(osm->lis,<og);CHKERRQ(ierr);
349 ierr = ISGetLocalSize(osm->is[i],&m);CHKERRQ(ierr);
350 ierr = ISGetIndices(osm->is[i], &idx_is);CHKERRQ(ierr);
351 ierr = PetscMalloc1(m,&idx_lis);CHKERRQ(ierr);
352 ierr = ISGlobalToLocalMappingApply(ltog,IS_GTOLM_DROP,m,idx_is,&nout,idx_lis);CHKERRQ(ierr);
353 if (nout != m) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"is not a subset of lis");
354 ierr = ISRestoreIndices(osm->is[i], &idx_is);CHKERRQ(ierr);
355 ierr = ISCreateGeneral(PETSC_COMM_SELF,m,idx_lis,PETSC_OWN_POINTER,&isll);CHKERRQ(ierr);
356 ierr = ISLocalToGlobalMappingDestroy(<og);CHKERRQ(ierr);
357 ierr = ISCreateStride(PETSC_COMM_SELF,m,0,1,&isl);CHKERRQ(ierr);
358 ierr = VecScatterCreate(osm->ly,isll,osm->y[i],isl,&osm->lrestriction[i]);CHKERRQ(ierr);
359 ierr = ISDestroy(&isll);CHKERRQ(ierr);
360 ierr = ISDestroy(&isl);CHKERRQ(ierr);
361 if (osm->lprolongation) { /* generate a scatter from y[i] to ly picking only the the non-overlapping is_local[i] entries */
362 ISLocalToGlobalMapping ltog;
363 IS isll,isll_local;
364 const PetscInt *idx_local;
365 PetscInt *idx1, *idx2, nout;
366
367 ierr = ISGetLocalSize(osm->is_local[i],&m_local);CHKERRQ(ierr);
368 ierr = ISGetIndices(osm->is_local[i], &idx_local);CHKERRQ(ierr);
369
370 ierr = ISLocalToGlobalMappingCreateIS(osm->is[i],<og);CHKERRQ(ierr);
371 ierr = PetscMalloc1(m_local,&idx1);CHKERRQ(ierr);
372 ierr = ISGlobalToLocalMappingApply(ltog,IS_GTOLM_DROP,m_local,idx_local,&nout,idx1);CHKERRQ(ierr);
373 ierr = ISLocalToGlobalMappingDestroy(<og);CHKERRQ(ierr);
374 if (nout != m_local) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"is_local not a subset of is");
375 ierr = ISCreateGeneral(PETSC_COMM_SELF,m_local,idx1,PETSC_OWN_POINTER,&isll);CHKERRQ(ierr);
376
377 ierr = ISLocalToGlobalMappingCreateIS(osm->lis,<og);CHKERRQ(ierr);
378 ierr = PetscMalloc1(m_local,&idx2);CHKERRQ(ierr);
379 ierr = ISGlobalToLocalMappingApply(ltog,IS_GTOLM_DROP,m_local,idx_local,&nout,idx2);CHKERRQ(ierr);
380 ierr = ISLocalToGlobalMappingDestroy(<og);CHKERRQ(ierr);
381 if (nout != m_local) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"is_local not a subset of lis");
382 ierr = ISCreateGeneral(PETSC_COMM_SELF,m_local,idx2,PETSC_OWN_POINTER,&isll_local);CHKERRQ(ierr);
383
384 ierr = ISRestoreIndices(osm->is_local[i], &idx_local);CHKERRQ(ierr);
385 ierr = VecScatterCreate(osm->y[i],isll,osm->ly,isll_local,&osm->lprolongation[i]);CHKERRQ(ierr);
386
387 ierr = ISDestroy(&isll);CHKERRQ(ierr);
388 ierr = ISDestroy(&isll_local);CHKERRQ(ierr);
389 }
390 }
391 ierr = VecDestroy(&vec);CHKERRQ(ierr);
392 }
393
394 if (osm->loctype == PC_COMPOSITE_MULTIPLICATIVE) {
395 IS *cis;
396 PetscInt c;
397
398 ierr = PetscMalloc1(osm->n_local_true, &cis);CHKERRQ(ierr);
399 for (c = 0; c < osm->n_local_true; ++c) cis[c] = osm->lis;
400 ierr = MatCreateSubMatrices(pc->pmat, osm->n_local_true, osm->is, cis, scall, &osm->lmats);CHKERRQ(ierr);
401 ierr = PetscFree(cis);CHKERRQ(ierr);
402 }
403
404 /* Return control to the user so that the submatrices can be modified (e.g., to apply
405 different boundary conditions for the submatrices than for the global problem) */
406 ierr = PCModifySubMatrices(pc,osm->n_local_true,osm->is,osm->is,osm->pmat,pc->modifysubmatricesP);CHKERRQ(ierr);
407
408 /*
409 Loop over subdomains putting them into local ksp
410 */
411 for (i=0; i<osm->n_local_true; i++) {
412 ierr = KSPSetOperators(osm->ksp[i],osm->pmat[i],osm->pmat[i]);CHKERRQ(ierr);
413 if (!pc->setupcalled) {
414 ierr = KSPSetFromOptions(osm->ksp[i]);CHKERRQ(ierr);
415 }
416 }
417 PetscFunctionReturn(0);
418 }
419
PCSetUpOnBlocks_ASM(PC pc)420 static PetscErrorCode PCSetUpOnBlocks_ASM(PC pc)
421 {
422 PC_ASM *osm = (PC_ASM*)pc->data;
423 PetscErrorCode ierr;
424 PetscInt i;
425 KSPConvergedReason reason;
426
427 PetscFunctionBegin;
428 for (i=0; i<osm->n_local_true; i++) {
429 ierr = KSPSetUp(osm->ksp[i]);CHKERRQ(ierr);
430 ierr = KSPGetConvergedReason(osm->ksp[i],&reason);CHKERRQ(ierr);
431 if (reason == KSP_DIVERGED_PC_FAILED) {
432 pc->failedreason = PC_SUBPC_ERROR;
433 }
434 }
435 PetscFunctionReturn(0);
436 }
437
PCApply_ASM(PC pc,Vec x,Vec y)438 static PetscErrorCode PCApply_ASM(PC pc,Vec x,Vec y)
439 {
440 PC_ASM *osm = (PC_ASM*)pc->data;
441 PetscErrorCode ierr;
442 PetscInt i,n_local_true = osm->n_local_true;
443 ScatterMode forward = SCATTER_FORWARD,reverse = SCATTER_REVERSE;
444
445 PetscFunctionBegin;
446 /*
447 support for limiting the restriction or interpolation to only local
448 subdomain values (leaving the other values 0).
449 */
450 if (!(osm->type & PC_ASM_RESTRICT)) {
451 forward = SCATTER_FORWARD_LOCAL;
452 /* have to zero the work RHS since scatter may leave some slots empty */
453 ierr = VecSet(osm->lx, 0.0);CHKERRQ(ierr);
454 }
455 if (!(osm->type & PC_ASM_INTERPOLATE)) {
456 reverse = SCATTER_REVERSE_LOCAL;
457 }
458
459 if (osm->loctype == PC_COMPOSITE_MULTIPLICATIVE || osm->loctype == PC_COMPOSITE_ADDITIVE) {
460 /* zero the global and the local solutions */
461 ierr = VecSet(y, 0.0);CHKERRQ(ierr);
462 ierr = VecSet(osm->ly, 0.0);CHKERRQ(ierr);
463
464 /* copy the global RHS to local RHS including the ghost nodes */
465 ierr = VecScatterBegin(osm->restriction, x, osm->lx, INSERT_VALUES, forward);CHKERRQ(ierr);
466 ierr = VecScatterEnd(osm->restriction, x, osm->lx, INSERT_VALUES, forward);CHKERRQ(ierr);
467
468 /* restrict local RHS to the overlapping 0-block RHS */
469 ierr = VecScatterBegin(osm->lrestriction[0], osm->lx, osm->x[0], INSERT_VALUES, forward);CHKERRQ(ierr);
470 ierr = VecScatterEnd(osm->lrestriction[0], osm->lx, osm->x[0], INSERT_VALUES, forward);CHKERRQ(ierr);
471
472 /* do the local solves */
473 for (i = 0; i < n_local_true; ++i) {
474
475 /* solve the overlapping i-block */
476 ierr = PetscLogEventBegin(PC_ApplyOnBlocks, osm->ksp[i], osm->x[i], osm->y[i],0);CHKERRQ(ierr);
477 ierr = KSPSolve(osm->ksp[i], osm->x[i], osm->y[i]);CHKERRQ(ierr);
478 ierr = KSPCheckSolve(osm->ksp[i], pc, osm->y[i]);CHKERRQ(ierr);
479 ierr = PetscLogEventEnd(PC_ApplyOnBlocks, osm->ksp[i], osm->x[i], osm->y[i], 0);CHKERRQ(ierr);
480
481 if (osm->lprolongation) { /* interpolate the non-overlapping i-block solution to the local solution (only for restrictive additive) */
482 ierr = VecScatterBegin(osm->lprolongation[i], osm->y[i], osm->ly, ADD_VALUES, forward);CHKERRQ(ierr);
483 ierr = VecScatterEnd(osm->lprolongation[i], osm->y[i], osm->ly, ADD_VALUES, forward);CHKERRQ(ierr);
484 } else { /* interpolate the overlapping i-block solution to the local solution */
485 ierr = VecScatterBegin(osm->lrestriction[i], osm->y[i], osm->ly, ADD_VALUES, reverse);CHKERRQ(ierr);
486 ierr = VecScatterEnd(osm->lrestriction[i], osm->y[i], osm->ly, ADD_VALUES, reverse);CHKERRQ(ierr);
487 }
488
489 if (i < n_local_true-1) {
490 /* restrict local RHS to the overlapping (i+1)-block RHS */
491 ierr = VecScatterBegin(osm->lrestriction[i+1], osm->lx, osm->x[i+1], INSERT_VALUES, forward);CHKERRQ(ierr);
492 ierr = VecScatterEnd(osm->lrestriction[i+1], osm->lx, osm->x[i+1], INSERT_VALUES, forward);CHKERRQ(ierr);
493
494 if (osm->loctype == PC_COMPOSITE_MULTIPLICATIVE) {
495 /* update the overlapping (i+1)-block RHS using the current local solution */
496 ierr = MatMult(osm->lmats[i+1], osm->ly, osm->y[i+1]);CHKERRQ(ierr);
497 ierr = VecAXPBY(osm->x[i+1],-1.,1., osm->y[i+1]);CHKERRQ(ierr);
498 }
499 }
500 }
501 /* add the local solution to the global solution including the ghost nodes */
502 ierr = VecScatterBegin(osm->restriction, osm->ly, y, ADD_VALUES, reverse);CHKERRQ(ierr);
503 ierr = VecScatterEnd(osm->restriction, osm->ly, y, ADD_VALUES, reverse);CHKERRQ(ierr);
504 } else SETERRQ1(PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_WRONG, "Invalid local composition type: %s", PCCompositeTypes[osm->loctype]);
505 PetscFunctionReturn(0);
506 }
507
PCMatApply_ASM(PC pc,Mat X,Mat Y)508 static PetscErrorCode PCMatApply_ASM(PC pc,Mat X,Mat Y)
509 {
510 PC_ASM *osm = (PC_ASM*)pc->data;
511 Mat Z,W;
512 Vec x;
513 PetscInt i,m,N;
514 ScatterMode forward = SCATTER_FORWARD,reverse = SCATTER_REVERSE;
515 PetscErrorCode ierr;
516
517 PetscFunctionBegin;
518 if (osm->n_local_true > 1) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP,"Not yet implemented");
519 /*
520 support for limiting the restriction or interpolation to only local
521 subdomain values (leaving the other values 0).
522 */
523 if (!(osm->type & PC_ASM_RESTRICT)) {
524 forward = SCATTER_FORWARD_LOCAL;
525 /* have to zero the work RHS since scatter may leave some slots empty */
526 ierr = VecSet(osm->lx, 0.0);CHKERRQ(ierr);
527 }
528 if (!(osm->type & PC_ASM_INTERPOLATE)) {
529 reverse = SCATTER_REVERSE_LOCAL;
530 }
531 ierr = VecGetLocalSize(osm->x[0], &m);CHKERRQ(ierr);
532 ierr = MatGetSize(X, NULL, &N);CHKERRQ(ierr);
533 ierr = MatCreateSeqDense(PETSC_COMM_SELF, m, N, NULL, &Z);CHKERRQ(ierr);
534 if (osm->loctype == PC_COMPOSITE_MULTIPLICATIVE || osm->loctype == PC_COMPOSITE_ADDITIVE) {
535 /* zero the global and the local solutions */
536 ierr = MatZeroEntries(Y);CHKERRQ(ierr);
537 ierr = VecSet(osm->ly, 0.0);CHKERRQ(ierr);
538
539 for (i = 0; i < N; ++i) {
540 ierr = MatDenseGetColumnVecRead(X, i, &x);
541 /* copy the global RHS to local RHS including the ghost nodes */
542 ierr = VecScatterBegin(osm->restriction, x, osm->lx, INSERT_VALUES, forward);CHKERRQ(ierr);
543 ierr = VecScatterEnd(osm->restriction, x, osm->lx, INSERT_VALUES, forward);CHKERRQ(ierr);
544 ierr = MatDenseRestoreColumnVecRead(X, i, &x);
545
546 ierr = MatDenseGetColumnVecWrite(Z, i, &x);
547 /* restrict local RHS to the overlapping 0-block RHS */
548 ierr = VecScatterBegin(osm->lrestriction[0], osm->lx, x, INSERT_VALUES, forward);CHKERRQ(ierr);
549 ierr = VecScatterEnd(osm->lrestriction[0], osm->lx, x, INSERT_VALUES, forward);CHKERRQ(ierr);
550 ierr = MatDenseRestoreColumnVecWrite(Z, i, &x);
551 }
552 ierr = MatCreateSeqDense(PETSC_COMM_SELF, m, N, NULL, &W);CHKERRQ(ierr);
553 /* solve the overlapping 0-block */
554 ierr = PetscLogEventBegin(PC_ApplyOnBlocks, osm->ksp[0], Z, W, 0);CHKERRQ(ierr);
555 ierr = KSPMatSolve(osm->ksp[0], Z, W);CHKERRQ(ierr);
556 ierr = KSPCheckSolve(osm->ksp[0], pc, NULL);CHKERRQ(ierr);
557 ierr = PetscLogEventEnd(PC_ApplyOnBlocks, osm->ksp[0], Z, W,0);CHKERRQ(ierr);
558 ierr = MatDestroy(&Z);CHKERRQ(ierr);
559
560 for (i = 0; i < N; ++i) {
561 ierr = VecSet(osm->ly, 0.0);CHKERRQ(ierr);
562 ierr = MatDenseGetColumnVecRead(W, i, &x);
563 if (osm->lprolongation) { /* interpolate the non-overlapping 0-block solution to the local solution (only for restrictive additive) */
564 ierr = VecScatterBegin(osm->lprolongation[0], x, osm->ly, ADD_VALUES, forward);CHKERRQ(ierr);
565 ierr = VecScatterEnd(osm->lprolongation[0], x, osm->ly, ADD_VALUES, forward);CHKERRQ(ierr);
566 } else { /* interpolate the overlapping 0-block solution to the local solution */
567 ierr = VecScatterBegin(osm->lrestriction[0], x, osm->ly, ADD_VALUES, reverse);CHKERRQ(ierr);
568 ierr = VecScatterEnd(osm->lrestriction[0], x, osm->ly, ADD_VALUES, reverse);CHKERRQ(ierr);
569 }
570 ierr = MatDenseRestoreColumnVecRead(W, i, &x);
571
572 ierr = MatDenseGetColumnVecWrite(Y, i, &x);
573 /* add the local solution to the global solution including the ghost nodes */
574 ierr = VecScatterBegin(osm->restriction, osm->ly, x, ADD_VALUES, reverse);CHKERRQ(ierr);
575 ierr = VecScatterEnd(osm->restriction, osm->ly, x, ADD_VALUES, reverse);CHKERRQ(ierr);
576 ierr = MatDenseRestoreColumnVecWrite(Y, i, &x);
577 }
578 ierr = MatDestroy(&W);CHKERRQ(ierr);
579 } else SETERRQ1(PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_WRONG, "Invalid local composition type: %s", PCCompositeTypes[osm->loctype]);
580 PetscFunctionReturn(0);
581 }
582
PCApplyTranspose_ASM(PC pc,Vec x,Vec y)583 static PetscErrorCode PCApplyTranspose_ASM(PC pc,Vec x,Vec y)
584 {
585 PC_ASM *osm = (PC_ASM*)pc->data;
586 PetscErrorCode ierr;
587 PetscInt i,n_local_true = osm->n_local_true;
588 ScatterMode forward = SCATTER_FORWARD,reverse = SCATTER_REVERSE;
589
590 PetscFunctionBegin;
591 /*
592 Support for limiting the restriction or interpolation to only local
593 subdomain values (leaving the other values 0).
594
595 Note: these are reversed from the PCApply_ASM() because we are applying the
596 transpose of the three terms
597 */
598
599 if (!(osm->type & PC_ASM_INTERPOLATE)) {
600 forward = SCATTER_FORWARD_LOCAL;
601 /* have to zero the work RHS since scatter may leave some slots empty */
602 ierr = VecSet(osm->lx, 0.0);CHKERRQ(ierr);
603 }
604 if (!(osm->type & PC_ASM_RESTRICT)) reverse = SCATTER_REVERSE_LOCAL;
605
606 /* zero the global and the local solutions */
607 ierr = VecSet(y, 0.0);CHKERRQ(ierr);
608 ierr = VecSet(osm->ly, 0.0);CHKERRQ(ierr);
609
610 /* Copy the global RHS to local RHS including the ghost nodes */
611 ierr = VecScatterBegin(osm->restriction, x, osm->lx, INSERT_VALUES, forward);CHKERRQ(ierr);
612 ierr = VecScatterEnd(osm->restriction, x, osm->lx, INSERT_VALUES, forward);CHKERRQ(ierr);
613
614 /* Restrict local RHS to the overlapping 0-block RHS */
615 ierr = VecScatterBegin(osm->lrestriction[0], osm->lx, osm->x[0], INSERT_VALUES, forward);CHKERRQ(ierr);
616 ierr = VecScatterEnd(osm->lrestriction[0], osm->lx, osm->x[0], INSERT_VALUES, forward);CHKERRQ(ierr);
617
618 /* do the local solves */
619 for (i = 0; i < n_local_true; ++i) {
620
621 /* solve the overlapping i-block */
622 ierr = PetscLogEventBegin(PC_ApplyOnBlocks,osm->ksp[i],osm->x[i],osm->y[i],0);CHKERRQ(ierr);
623 ierr = KSPSolveTranspose(osm->ksp[i], osm->x[i], osm->y[i]);CHKERRQ(ierr);
624 ierr = KSPCheckSolve(osm->ksp[i],pc,osm->y[i]);CHKERRQ(ierr);
625 ierr = PetscLogEventEnd(PC_ApplyOnBlocks,osm->ksp[i],osm->x[i],osm->y[i],0);CHKERRQ(ierr);
626
627 if (osm->lprolongation) { /* interpolate the non-overlapping i-block solution to the local solution */
628 ierr = VecScatterBegin(osm->lprolongation[i], osm->y[i], osm->ly, ADD_VALUES, forward);CHKERRQ(ierr);
629 ierr = VecScatterEnd(osm->lprolongation[i], osm->y[i], osm->ly, ADD_VALUES, forward);CHKERRQ(ierr);
630 } else { /* interpolate the overlapping i-block solution to the local solution */
631 ierr = VecScatterBegin(osm->lrestriction[i], osm->y[i], osm->ly, ADD_VALUES, reverse);CHKERRQ(ierr);
632 ierr = VecScatterEnd(osm->lrestriction[i], osm->y[i], osm->ly, ADD_VALUES, reverse);CHKERRQ(ierr);
633 }
634
635 if (i < n_local_true-1) {
636 /* Restrict local RHS to the overlapping (i+1)-block RHS */
637 ierr = VecScatterBegin(osm->lrestriction[i+1], osm->lx, osm->x[i+1], INSERT_VALUES, forward);CHKERRQ(ierr);
638 ierr = VecScatterEnd(osm->lrestriction[i+1], osm->lx, osm->x[i+1], INSERT_VALUES, forward);CHKERRQ(ierr);
639 }
640 }
641 /* Add the local solution to the global solution including the ghost nodes */
642 ierr = VecScatterBegin(osm->restriction, osm->ly, y, ADD_VALUES, reverse);CHKERRQ(ierr);
643 ierr = VecScatterEnd(osm->restriction, osm->ly, y, ADD_VALUES, reverse);CHKERRQ(ierr);
644
645 PetscFunctionReturn(0);
646
647 }
648
PCReset_ASM(PC pc)649 static PetscErrorCode PCReset_ASM(PC pc)
650 {
651 PC_ASM *osm = (PC_ASM*)pc->data;
652 PetscErrorCode ierr;
653 PetscInt i;
654
655 PetscFunctionBegin;
656 if (osm->ksp) {
657 for (i=0; i<osm->n_local_true; i++) {
658 ierr = KSPReset(osm->ksp[i]);CHKERRQ(ierr);
659 }
660 }
661 if (osm->pmat) {
662 if (osm->n_local_true > 0) {
663 ierr = MatDestroySubMatrices(osm->n_local_true,&osm->pmat);CHKERRQ(ierr);
664 }
665 }
666 if (osm->lrestriction) {
667 ierr = VecScatterDestroy(&osm->restriction);CHKERRQ(ierr);
668 for (i=0; i<osm->n_local_true; i++) {
669 ierr = VecScatterDestroy(&osm->lrestriction[i]);CHKERRQ(ierr);
670 if (osm->lprolongation) {ierr = VecScatterDestroy(&osm->lprolongation[i]);CHKERRQ(ierr);}
671 ierr = VecDestroy(&osm->x[i]);CHKERRQ(ierr);
672 ierr = VecDestroy(&osm->y[i]);CHKERRQ(ierr);
673 }
674 ierr = PetscFree(osm->lrestriction);CHKERRQ(ierr);
675 if (osm->lprolongation) {ierr = PetscFree(osm->lprolongation);CHKERRQ(ierr);}
676 ierr = PetscFree(osm->x);CHKERRQ(ierr);
677 ierr = PetscFree(osm->y);CHKERRQ(ierr);
678
679 }
680 ierr = PCASMDestroySubdomains(osm->n_local_true,osm->is,osm->is_local);CHKERRQ(ierr);
681 ierr = ISDestroy(&osm->lis);CHKERRQ(ierr);
682 ierr = VecDestroy(&osm->lx);CHKERRQ(ierr);
683 ierr = VecDestroy(&osm->ly);CHKERRQ(ierr);
684 if (osm->loctype == PC_COMPOSITE_MULTIPLICATIVE) {
685 ierr = MatDestroyMatrices(osm->n_local_true, &osm->lmats);CHKERRQ(ierr);
686 }
687
688 ierr = PetscFree(osm->sub_mat_type);CHKERRQ(ierr);
689
690 osm->is = NULL;
691 osm->is_local = NULL;
692 PetscFunctionReturn(0);
693 }
694
PCDestroy_ASM(PC pc)695 static PetscErrorCode PCDestroy_ASM(PC pc)
696 {
697 PC_ASM *osm = (PC_ASM*)pc->data;
698 PetscErrorCode ierr;
699 PetscInt i;
700
701 PetscFunctionBegin;
702 ierr = PCReset_ASM(pc);CHKERRQ(ierr);
703 if (osm->ksp) {
704 for (i=0; i<osm->n_local_true; i++) {
705 ierr = KSPDestroy(&osm->ksp[i]);CHKERRQ(ierr);
706 }
707 ierr = PetscFree(osm->ksp);CHKERRQ(ierr);
708 }
709 ierr = PetscFree(pc->data);CHKERRQ(ierr);
710 PetscFunctionReturn(0);
711 }
712
PCSetFromOptions_ASM(PetscOptionItems * PetscOptionsObject,PC pc)713 static PetscErrorCode PCSetFromOptions_ASM(PetscOptionItems *PetscOptionsObject,PC pc)
714 {
715 PC_ASM *osm = (PC_ASM*)pc->data;
716 PetscErrorCode ierr;
717 PetscInt blocks,ovl;
718 PetscBool flg;
719 PCASMType asmtype;
720 PCCompositeType loctype;
721 char sub_mat_type[256];
722
723 PetscFunctionBegin;
724 ierr = PetscOptionsHead(PetscOptionsObject,"Additive Schwarz options");CHKERRQ(ierr);
725 ierr = PetscOptionsBool("-pc_asm_dm_subdomains","Use DMCreateDomainDecomposition() to define subdomains","PCASMSetDMSubdomains",osm->dm_subdomains,&osm->dm_subdomains,&flg);CHKERRQ(ierr);
726 ierr = PetscOptionsInt("-pc_asm_blocks","Number of subdomains","PCASMSetTotalSubdomains",osm->n,&blocks,&flg);CHKERRQ(ierr);
727 if (flg) {
728 ierr = PCASMSetTotalSubdomains(pc,blocks,NULL,NULL);CHKERRQ(ierr);
729 osm->dm_subdomains = PETSC_FALSE;
730 }
731 ierr = PetscOptionsInt("-pc_asm_local_blocks","Number of local subdomains","PCASMSetLocalSubdomains",osm->n_local_true,&blocks,&flg);CHKERRQ(ierr);
732 if (flg) {
733 ierr = PCASMSetLocalSubdomains(pc,blocks,NULL,NULL);CHKERRQ(ierr);
734 osm->dm_subdomains = PETSC_FALSE;
735 }
736 ierr = PetscOptionsInt("-pc_asm_overlap","Number of grid points overlap","PCASMSetOverlap",osm->overlap,&ovl,&flg);CHKERRQ(ierr);
737 if (flg) {
738 ierr = PCASMSetOverlap(pc,ovl);CHKERRQ(ierr);
739 osm->dm_subdomains = PETSC_FALSE;
740 }
741 flg = PETSC_FALSE;
742 ierr = PetscOptionsEnum("-pc_asm_type","Type of restriction/extension","PCASMSetType",PCASMTypes,(PetscEnum)osm->type,(PetscEnum*)&asmtype,&flg);CHKERRQ(ierr);
743 if (flg) {ierr = PCASMSetType(pc,asmtype);CHKERRQ(ierr); }
744 flg = PETSC_FALSE;
745 ierr = PetscOptionsEnum("-pc_asm_local_type","Type of local solver composition","PCASMSetLocalType",PCCompositeTypes,(PetscEnum)osm->loctype,(PetscEnum*)&loctype,&flg);CHKERRQ(ierr);
746 if (flg) {ierr = PCASMSetLocalType(pc,loctype);CHKERRQ(ierr); }
747 ierr = PetscOptionsFList("-pc_asm_sub_mat_type","Subsolve Matrix Type","PCASMSetSubMatType",MatList,NULL,sub_mat_type,256,&flg);CHKERRQ(ierr);
748 if (flg) {
749 ierr = PCASMSetSubMatType(pc,sub_mat_type);CHKERRQ(ierr);
750 }
751 ierr = PetscOptionsTail();CHKERRQ(ierr);
752 PetscFunctionReturn(0);
753 }
754
755 /*------------------------------------------------------------------------------------*/
756
PCASMSetLocalSubdomains_ASM(PC pc,PetscInt n,IS is[],IS is_local[])757 static PetscErrorCode PCASMSetLocalSubdomains_ASM(PC pc,PetscInt n,IS is[],IS is_local[])
758 {
759 PC_ASM *osm = (PC_ASM*)pc->data;
760 PetscErrorCode ierr;
761 PetscInt i;
762
763 PetscFunctionBegin;
764 if (n < 1) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Each process must have 1 or more blocks, n = %D",n);
765 if (pc->setupcalled && (n != osm->n_local_true || is)) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONGSTATE,"PCASMSetLocalSubdomains() should be called before calling PCSetUp().");
766
767 if (!pc->setupcalled) {
768 if (is) {
769 for (i=0; i<n; i++) {ierr = PetscObjectReference((PetscObject)is[i]);CHKERRQ(ierr);}
770 }
771 if (is_local) {
772 for (i=0; i<n; i++) {ierr = PetscObjectReference((PetscObject)is_local[i]);CHKERRQ(ierr);}
773 }
774 ierr = PCASMDestroySubdomains(osm->n_local_true,osm->is,osm->is_local);CHKERRQ(ierr);
775
776 osm->n_local_true = n;
777 osm->is = NULL;
778 osm->is_local = NULL;
779 if (is) {
780 ierr = PetscMalloc1(n,&osm->is);CHKERRQ(ierr);
781 for (i=0; i<n; i++) osm->is[i] = is[i];
782 /* Flag indicating that the user has set overlapping subdomains so PCASM should not increase their size. */
783 osm->overlap = -1;
784 }
785 if (is_local) {
786 ierr = PetscMalloc1(n,&osm->is_local);CHKERRQ(ierr);
787 for (i=0; i<n; i++) osm->is_local[i] = is_local[i];
788 if (!is) {
789 ierr = PetscMalloc1(osm->n_local_true,&osm->is);CHKERRQ(ierr);
790 for (i=0; i<osm->n_local_true; i++) {
791 if (osm->overlap > 0) { /* With positive overlap, osm->is[i] will be modified */
792 ierr = ISDuplicate(osm->is_local[i],&osm->is[i]);CHKERRQ(ierr);
793 ierr = ISCopy(osm->is_local[i],osm->is[i]);CHKERRQ(ierr);
794 } else {
795 ierr = PetscObjectReference((PetscObject)osm->is_local[i]);CHKERRQ(ierr);
796 osm->is[i] = osm->is_local[i];
797 }
798 }
799 }
800 }
801 }
802 PetscFunctionReturn(0);
803 }
804
PCASMSetTotalSubdomains_ASM(PC pc,PetscInt N,IS * is,IS * is_local)805 static PetscErrorCode PCASMSetTotalSubdomains_ASM(PC pc,PetscInt N,IS *is,IS *is_local)
806 {
807 PC_ASM *osm = (PC_ASM*)pc->data;
808 PetscErrorCode ierr;
809 PetscMPIInt rank,size;
810 PetscInt n;
811
812 PetscFunctionBegin;
813 if (N < 1) SETERRQ1(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_OUTOFRANGE,"Number of total blocks must be > 0, N = %D",N);
814 if (is || is_local) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP,"Use PCASMSetLocalSubdomains() to set specific index sets\n\they cannot be set globally yet.");
815
816 /*
817 Split the subdomains equally among all processors
818 */
819 ierr = MPI_Comm_rank(PetscObjectComm((PetscObject)pc),&rank);CHKERRQ(ierr);
820 ierr = MPI_Comm_size(PetscObjectComm((PetscObject)pc),&size);CHKERRQ(ierr);
821 n = N/size + ((N % size) > rank);
822 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);
823 if (pc->setupcalled && n != osm->n_local_true) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"PCASMSetTotalSubdomains() should be called before PCSetUp().");
824 if (!pc->setupcalled) {
825 ierr = PCASMDestroySubdomains(osm->n_local_true,osm->is,osm->is_local);CHKERRQ(ierr);
826
827 osm->n_local_true = n;
828 osm->is = NULL;
829 osm->is_local = NULL;
830 }
831 PetscFunctionReturn(0);
832 }
833
PCASMSetOverlap_ASM(PC pc,PetscInt ovl)834 static PetscErrorCode PCASMSetOverlap_ASM(PC pc,PetscInt ovl)
835 {
836 PC_ASM *osm = (PC_ASM*)pc->data;
837
838 PetscFunctionBegin;
839 if (ovl < 0) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_OUTOFRANGE,"Negative overlap value requested");
840 if (pc->setupcalled && ovl != osm->overlap) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONGSTATE,"PCASMSetOverlap() should be called before PCSetUp().");
841 if (!pc->setupcalled) osm->overlap = ovl;
842 PetscFunctionReturn(0);
843 }
844
PCASMSetType_ASM(PC pc,PCASMType type)845 static PetscErrorCode PCASMSetType_ASM(PC pc,PCASMType type)
846 {
847 PC_ASM *osm = (PC_ASM*)pc->data;
848
849 PetscFunctionBegin;
850 osm->type = type;
851 osm->type_set = PETSC_TRUE;
852 PetscFunctionReturn(0);
853 }
854
PCASMGetType_ASM(PC pc,PCASMType * type)855 static PetscErrorCode PCASMGetType_ASM(PC pc,PCASMType *type)
856 {
857 PC_ASM *osm = (PC_ASM*)pc->data;
858
859 PetscFunctionBegin;
860 *type = osm->type;
861 PetscFunctionReturn(0);
862 }
863
PCASMSetLocalType_ASM(PC pc,PCCompositeType type)864 static PetscErrorCode PCASMSetLocalType_ASM(PC pc, PCCompositeType type)
865 {
866 PC_ASM *osm = (PC_ASM *) pc->data;
867
868 PetscFunctionBegin;
869 if (type != PC_COMPOSITE_ADDITIVE && type != PC_COMPOSITE_MULTIPLICATIVE) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP,"Only supports additive or multiplicative as the local type");
870 osm->loctype = type;
871 PetscFunctionReturn(0);
872 }
873
PCASMGetLocalType_ASM(PC pc,PCCompositeType * type)874 static PetscErrorCode PCASMGetLocalType_ASM(PC pc, PCCompositeType *type)
875 {
876 PC_ASM *osm = (PC_ASM *) pc->data;
877
878 PetscFunctionBegin;
879 *type = osm->loctype;
880 PetscFunctionReturn(0);
881 }
882
PCASMSetSortIndices_ASM(PC pc,PetscBool doSort)883 static PetscErrorCode PCASMSetSortIndices_ASM(PC pc,PetscBool doSort)
884 {
885 PC_ASM *osm = (PC_ASM*)pc->data;
886
887 PetscFunctionBegin;
888 osm->sort_indices = doSort;
889 PetscFunctionReturn(0);
890 }
891
PCASMGetSubKSP_ASM(PC pc,PetscInt * n_local,PetscInt * first_local,KSP ** ksp)892 static PetscErrorCode PCASMGetSubKSP_ASM(PC pc,PetscInt *n_local,PetscInt *first_local,KSP **ksp)
893 {
894 PC_ASM *osm = (PC_ASM*)pc->data;
895 PetscErrorCode ierr;
896
897 PetscFunctionBegin;
898 if (osm->n_local_true < 1) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ORDER,"Need to call PCSetUp() on PC (or KSPSetUp() on the outer KSP object) before calling here");
899
900 if (n_local) *n_local = osm->n_local_true;
901 if (first_local) {
902 ierr = MPI_Scan(&osm->n_local_true,first_local,1,MPIU_INT,MPI_SUM,PetscObjectComm((PetscObject)pc));CHKERRQ(ierr);
903 *first_local -= osm->n_local_true;
904 }
905 if (ksp) {
906 /* Assume that local solves are now different; not necessarily
907 true though! This flag is used only for PCView_ASM() */
908 *ksp = osm->ksp;
909 osm->same_local_solves = PETSC_FALSE;
910 }
911 PetscFunctionReturn(0);
912 }
913
PCASMGetSubMatType_ASM(PC pc,MatType * sub_mat_type)914 static PetscErrorCode PCASMGetSubMatType_ASM(PC pc,MatType *sub_mat_type)
915 {
916 PC_ASM *osm = (PC_ASM*)pc->data;
917
918 PetscFunctionBegin;
919 PetscValidHeaderSpecific(pc,PC_CLASSID,1);
920 PetscValidPointer(sub_mat_type,2);
921 *sub_mat_type = osm->sub_mat_type;
922 PetscFunctionReturn(0);
923 }
924
PCASMSetSubMatType_ASM(PC pc,MatType sub_mat_type)925 static PetscErrorCode PCASMSetSubMatType_ASM(PC pc,MatType sub_mat_type)
926 {
927 PetscErrorCode ierr;
928 PC_ASM *osm = (PC_ASM*)pc->data;
929
930 PetscFunctionBegin;
931 PetscValidHeaderSpecific(pc,PC_CLASSID,1);
932 ierr = PetscFree(osm->sub_mat_type);CHKERRQ(ierr);
933 ierr = PetscStrallocpy(sub_mat_type,(char**)&osm->sub_mat_type);CHKERRQ(ierr);
934 PetscFunctionReturn(0);
935 }
936
937 /*@C
938 PCASMSetLocalSubdomains - Sets the local subdomains (for this processor only) for the additive Schwarz preconditioner.
939
940 Collective on pc
941
942 Input Parameters:
943 + pc - the preconditioner context
944 . n - the number of subdomains for this processor (default value = 1)
945 . is - the index set that defines the subdomains for this processor
946 (or NULL for PETSc to determine subdomains)
947 - is_local - the index sets that define the local part of the subdomains for this processor, not used unless PCASMType is PC_ASM_RESTRICT
948 (or NULL to not provide these)
949
950 Options Database Key:
951 To set the total number of subdomain blocks rather than specify the
952 index sets, use the option
953 . -pc_asm_local_blocks <blks> - Sets local blocks
954
955 Notes:
956 The IS numbering is in the parallel, global numbering of the vector for both is and is_local
957
958 By default the ASM preconditioner uses 1 block per processor.
959
960 Use PCASMSetTotalSubdomains() to set the subdomains for all processors.
961
962 If is_local is provided and PCASMType is PC_ASM_RESTRICT then the solution only over the is_local region is interpolated
963 back to form the global solution (this is the standard restricted additive Schwarz method)
964
965 If the is_local is provided and PCASMType is PC_ASM_INTERPOLATE or PC_ASM_NONE then an error is generated since there is
966 no code to handle that case.
967
968 Level: advanced
969
970 .seealso: PCASMSetTotalSubdomains(), PCASMSetOverlap(), PCASMGetSubKSP(),
971 PCASMCreateSubdomains2D(), PCASMGetLocalSubdomains(), PCASMType, PCASMSetType()
972 @*/
PCASMSetLocalSubdomains(PC pc,PetscInt n,IS is[],IS is_local[])973 PetscErrorCode PCASMSetLocalSubdomains(PC pc,PetscInt n,IS is[],IS is_local[])
974 {
975 PetscErrorCode ierr;
976
977 PetscFunctionBegin;
978 PetscValidHeaderSpecific(pc,PC_CLASSID,1);
979 ierr = PetscTryMethod(pc,"PCASMSetLocalSubdomains_C",(PC,PetscInt,IS[],IS[]),(pc,n,is,is_local));CHKERRQ(ierr);
980 PetscFunctionReturn(0);
981 }
982
983 /*@C
984 PCASMSetTotalSubdomains - Sets the subdomains for all processors for the
985 additive Schwarz preconditioner. Either all or no processors in the
986 PC communicator must call this routine, with the same index sets.
987
988 Collective on pc
989
990 Input Parameters:
991 + pc - the preconditioner context
992 . N - the number of subdomains for all processors
993 . is - the index sets that define the subdomains for all processors
994 (or NULL to ask PETSc to determine the subdomains)
995 - is_local - the index sets that define the local part of the subdomains for this processor
996 (or NULL to not provide this information)
997
998 Options Database Key:
999 To set the total number of subdomain blocks rather than specify the
1000 index sets, use the option
1001 . -pc_asm_blocks <blks> - Sets total blocks
1002
1003 Notes:
1004 Currently you cannot use this to set the actual subdomains with the argument is or is_local.
1005
1006 By default the ASM preconditioner uses 1 block per processor.
1007
1008 These index sets cannot be destroyed until after completion of the
1009 linear solves for which the ASM preconditioner is being used.
1010
1011 Use PCASMSetLocalSubdomains() to set local subdomains.
1012
1013 The IS numbering is in the parallel, global numbering of the vector for both is and is_local
1014
1015 Level: advanced
1016
1017 .seealso: PCASMSetLocalSubdomains(), PCASMSetOverlap(), PCASMGetSubKSP(),
1018 PCASMCreateSubdomains2D()
1019 @*/
PCASMSetTotalSubdomains(PC pc,PetscInt N,IS is[],IS is_local[])1020 PetscErrorCode PCASMSetTotalSubdomains(PC pc,PetscInt N,IS is[],IS is_local[])
1021 {
1022 PetscErrorCode ierr;
1023
1024 PetscFunctionBegin;
1025 PetscValidHeaderSpecific(pc,PC_CLASSID,1);
1026 ierr = PetscTryMethod(pc,"PCASMSetTotalSubdomains_C",(PC,PetscInt,IS[],IS[]),(pc,N,is,is_local));CHKERRQ(ierr);
1027 PetscFunctionReturn(0);
1028 }
1029
1030 /*@
1031 PCASMSetOverlap - Sets the overlap between a pair of subdomains for the
1032 additive Schwarz preconditioner. Either all or no processors in the
1033 PC communicator must call this routine.
1034
1035 Logically Collective on pc
1036
1037 Input Parameters:
1038 + pc - the preconditioner context
1039 - ovl - the amount of overlap between subdomains (ovl >= 0, default value = 1)
1040
1041 Options Database Key:
1042 . -pc_asm_overlap <ovl> - Sets overlap
1043
1044 Notes:
1045 By default the ASM preconditioner uses 1 block per processor. To use
1046 multiple blocks per perocessor, see PCASMSetTotalSubdomains() and
1047 PCASMSetLocalSubdomains() (and the option -pc_asm_blocks <blks>).
1048
1049 The overlap defaults to 1, so if one desires that no additional
1050 overlap be computed beyond what may have been set with a call to
1051 PCASMSetTotalSubdomains() or PCASMSetLocalSubdomains(), then ovl
1052 must be set to be 0. In particular, if one does not explicitly set
1053 the subdomains an application code, then all overlap would be computed
1054 internally by PETSc, and using an overlap of 0 would result in an ASM
1055 variant that is equivalent to the block Jacobi preconditioner.
1056
1057 The default algorithm used by PETSc to increase overlap is fast, but not scalable,
1058 use the option -mat_increase_overlap_scalable when the problem and number of processes is large.
1059
1060 Note that one can define initial index sets with any overlap via
1061 PCASMSetLocalSubdomains(); the routine
1062 PCASMSetOverlap() merely allows PETSc to extend that overlap further
1063 if desired.
1064
1065 Level: intermediate
1066
1067 .seealso: PCASMSetTotalSubdomains(), PCASMSetLocalSubdomains(), PCASMGetSubKSP(),
1068 PCASMCreateSubdomains2D(), PCASMGetLocalSubdomains(), MatIncreaseOverlap()
1069 @*/
PCASMSetOverlap(PC pc,PetscInt ovl)1070 PetscErrorCode PCASMSetOverlap(PC pc,PetscInt ovl)
1071 {
1072 PetscErrorCode ierr;
1073
1074 PetscFunctionBegin;
1075 PetscValidHeaderSpecific(pc,PC_CLASSID,1);
1076 PetscValidLogicalCollectiveInt(pc,ovl,2);
1077 ierr = PetscTryMethod(pc,"PCASMSetOverlap_C",(PC,PetscInt),(pc,ovl));CHKERRQ(ierr);
1078 PetscFunctionReturn(0);
1079 }
1080
1081 /*@
1082 PCASMSetType - Sets the type of restriction and interpolation used
1083 for local problems in the additive Schwarz method.
1084
1085 Logically Collective on pc
1086
1087 Input Parameters:
1088 + pc - the preconditioner context
1089 - type - variant of ASM, one of
1090 .vb
1091 PC_ASM_BASIC - full interpolation and restriction
1092 PC_ASM_RESTRICT - full restriction, local processor interpolation
1093 PC_ASM_INTERPOLATE - full interpolation, local processor restriction
1094 PC_ASM_NONE - local processor restriction and interpolation
1095 .ve
1096
1097 Options Database Key:
1098 . -pc_asm_type [basic,restrict,interpolate,none] - Sets ASM type
1099
1100 Notes:
1101 if the is_local arguments are passed to PCASMSetLocalSubdomains() then they are used when PC_ASM_RESTRICT has been selected
1102 to limit the local processor interpolation
1103
1104 Level: intermediate
1105
1106 .seealso: PCASMSetTotalSubdomains(), PCASMSetTotalSubdomains(), PCASMGetSubKSP(),
1107 PCASMCreateSubdomains2D(), PCASMType, PCASMSetLocalType(), PCASMGetLocalType()
1108 @*/
PCASMSetType(PC pc,PCASMType type)1109 PetscErrorCode PCASMSetType(PC pc,PCASMType type)
1110 {
1111 PetscErrorCode ierr;
1112
1113 PetscFunctionBegin;
1114 PetscValidHeaderSpecific(pc,PC_CLASSID,1);
1115 PetscValidLogicalCollectiveEnum(pc,type,2);
1116 ierr = PetscTryMethod(pc,"PCASMSetType_C",(PC,PCASMType),(pc,type));CHKERRQ(ierr);
1117 PetscFunctionReturn(0);
1118 }
1119
1120 /*@
1121 PCASMGetType - Gets the type of restriction and interpolation used
1122 for local problems in the additive Schwarz method.
1123
1124 Logically Collective on pc
1125
1126 Input Parameter:
1127 . pc - the preconditioner context
1128
1129 Output Parameter:
1130 . type - variant of ASM, one of
1131
1132 .vb
1133 PC_ASM_BASIC - full interpolation and restriction
1134 PC_ASM_RESTRICT - full restriction, local processor interpolation
1135 PC_ASM_INTERPOLATE - full interpolation, local processor restriction
1136 PC_ASM_NONE - local processor restriction and interpolation
1137 .ve
1138
1139 Options Database Key:
1140 . -pc_asm_type [basic,restrict,interpolate,none] - Sets ASM type
1141
1142 Level: intermediate
1143
1144 .seealso: PCASMSetTotalSubdomains(), PCASMSetTotalSubdomains(), PCASMGetSubKSP(),
1145 PCASMCreateSubdomains2D(), PCASMType, PCASMSetType(), PCASMSetLocalType(), PCASMGetLocalType()
1146 @*/
PCASMGetType(PC pc,PCASMType * type)1147 PetscErrorCode PCASMGetType(PC pc,PCASMType *type)
1148 {
1149 PetscErrorCode ierr;
1150
1151 PetscFunctionBegin;
1152 PetscValidHeaderSpecific(pc,PC_CLASSID,1);
1153 ierr = PetscUseMethod(pc,"PCASMGetType_C",(PC,PCASMType*),(pc,type));CHKERRQ(ierr);
1154 PetscFunctionReturn(0);
1155 }
1156
1157 /*@
1158 PCASMSetLocalType - Sets the type of composition used for local problems in the additive Schwarz method.
1159
1160 Logically Collective on pc
1161
1162 Input Parameters:
1163 + pc - the preconditioner context
1164 - type - type of composition, one of
1165 .vb
1166 PC_COMPOSITE_ADDITIVE - local additive combination
1167 PC_COMPOSITE_MULTIPLICATIVE - local multiplicative combination
1168 .ve
1169
1170 Options Database Key:
1171 . -pc_asm_local_type [additive,multiplicative] - Sets local solver composition type
1172
1173 Level: intermediate
1174
1175 .seealso: PCASMSetType(), PCASMGetType(), PCASMGetLocalType(), PCASM, PCASMType, PCASMSetType(), PCASMGetType(), PCCompositeType
1176 @*/
PCASMSetLocalType(PC pc,PCCompositeType type)1177 PetscErrorCode PCASMSetLocalType(PC pc, PCCompositeType type)
1178 {
1179 PetscErrorCode ierr;
1180
1181 PetscFunctionBegin;
1182 PetscValidHeaderSpecific(pc, PC_CLASSID, 1);
1183 PetscValidLogicalCollectiveEnum(pc, type, 2);
1184 ierr = PetscTryMethod(pc, "PCASMSetLocalType_C", (PC, PCCompositeType), (pc, type));CHKERRQ(ierr);
1185 PetscFunctionReturn(0);
1186 }
1187
1188 /*@
1189 PCASMGetLocalType - Gets the type of composition used for local problems in the additive Schwarz method.
1190
1191 Logically Collective on pc
1192
1193 Input Parameter:
1194 . pc - the preconditioner context
1195
1196 Output Parameter:
1197 . type - type of composition, one of
1198 .vb
1199 PC_COMPOSITE_ADDITIVE - local additive combination
1200 PC_COMPOSITE_MULTIPLICATIVE - local multiplicative combination
1201 .ve
1202
1203 Options Database Key:
1204 . -pc_asm_local_type [additive,multiplicative] - Sets local solver composition type
1205
1206 Level: intermediate
1207
1208 .seealso: PCASMSetType(), PCASMGetType(), PCASMSetLocalType(), PCASMCreate(), PCASMType, PCASMSetType(), PCASMGetType(), PCCompositeType
1209 @*/
PCASMGetLocalType(PC pc,PCCompositeType * type)1210 PetscErrorCode PCASMGetLocalType(PC pc, PCCompositeType *type)
1211 {
1212 PetscErrorCode ierr;
1213
1214 PetscFunctionBegin;
1215 PetscValidHeaderSpecific(pc, PC_CLASSID, 1);
1216 PetscValidPointer(type, 2);
1217 ierr = PetscUseMethod(pc, "PCASMGetLocalType_C", (PC, PCCompositeType *), (pc, type));CHKERRQ(ierr);
1218 PetscFunctionReturn(0);
1219 }
1220
1221 /*@
1222 PCASMSetSortIndices - Determines whether subdomain indices are sorted.
1223
1224 Logically Collective on pc
1225
1226 Input Parameters:
1227 + pc - the preconditioner context
1228 - doSort - sort the subdomain indices
1229
1230 Level: intermediate
1231
1232 .seealso: PCASMSetLocalSubdomains(), PCASMSetTotalSubdomains(), PCASMGetSubKSP(),
1233 PCASMCreateSubdomains2D()
1234 @*/
PCASMSetSortIndices(PC pc,PetscBool doSort)1235 PetscErrorCode PCASMSetSortIndices(PC pc,PetscBool doSort)
1236 {
1237 PetscErrorCode ierr;
1238
1239 PetscFunctionBegin;
1240 PetscValidHeaderSpecific(pc,PC_CLASSID,1);
1241 PetscValidLogicalCollectiveBool(pc,doSort,2);
1242 ierr = PetscTryMethod(pc,"PCASMSetSortIndices_C",(PC,PetscBool),(pc,doSort));CHKERRQ(ierr);
1243 PetscFunctionReturn(0);
1244 }
1245
1246 /*@C
1247 PCASMGetSubKSP - Gets the local KSP contexts for all blocks on
1248 this processor.
1249
1250 Collective on pc iff first_local is requested
1251
1252 Input Parameter:
1253 . pc - the preconditioner context
1254
1255 Output Parameters:
1256 + n_local - the number of blocks on this processor or NULL
1257 . first_local - the global number of the first block on this processor or NULL,
1258 all processors must request or all must pass NULL
1259 - ksp - the array of KSP contexts
1260
1261 Note:
1262 After PCASMGetSubKSP() the array of KSPes is not to be freed.
1263
1264 You must call KSPSetUp() before calling PCASMGetSubKSP().
1265
1266 Fortran note:
1267 The output argument 'ksp' must be an array of sufficient length or PETSC_NULL_KSP. The latter can be used to learn the necessary length.
1268
1269 Level: advanced
1270
1271 .seealso: PCASMSetTotalSubdomains(), PCASMSetTotalSubdomains(), PCASMSetOverlap(),
1272 PCASMCreateSubdomains2D(),
1273 @*/
PCASMGetSubKSP(PC pc,PetscInt * n_local,PetscInt * first_local,KSP * ksp[])1274 PetscErrorCode PCASMGetSubKSP(PC pc,PetscInt *n_local,PetscInt *first_local,KSP *ksp[])
1275 {
1276 PetscErrorCode ierr;
1277
1278 PetscFunctionBegin;
1279 PetscValidHeaderSpecific(pc,PC_CLASSID,1);
1280 ierr = PetscUseMethod(pc,"PCASMGetSubKSP_C",(PC,PetscInt*,PetscInt*,KSP **),(pc,n_local,first_local,ksp));CHKERRQ(ierr);
1281 PetscFunctionReturn(0);
1282 }
1283
1284 /* -------------------------------------------------------------------------------------*/
1285 /*MC
1286 PCASM - Use the (restricted) additive Schwarz method, each block is (approximately) solved with
1287 its own KSP object.
1288
1289 Options Database Keys:
1290 + -pc_asm_blocks <blks> - Sets total blocks
1291 . -pc_asm_overlap <ovl> - Sets overlap
1292 . -pc_asm_type [basic,restrict,interpolate,none] - Sets ASM type, default is restrict
1293 - -pc_asm_local_type [additive, multiplicative] - Sets ASM type, default is additive
1294
1295 IMPORTANT: If you run with, for example, 3 blocks on 1 processor or 3 blocks on 3 processors you
1296 will get a different convergence rate due to the default option of -pc_asm_type restrict. Use
1297 -pc_asm_type basic to use the standard ASM.
1298
1299 Notes:
1300 Each processor can have one or more blocks, but a block cannot be shared by more
1301 than one processor. Use PCGASM for subdomains shared by multiple processes. Defaults to one block per processor.
1302
1303 To set options on the solvers for each block append -sub_ to all the KSP, and PC
1304 options database keys. For example, -sub_pc_type ilu -sub_pc_factor_levels 1 -sub_ksp_type preonly
1305
1306 To set the options on the solvers separate for each block call PCASMGetSubKSP()
1307 and set the options directly on the resulting KSP object (you can access its PC
1308 with KSPGetPC())
1309
1310 Level: beginner
1311
1312 References:
1313 + 1. - M Dryja, OB Widlund, An additive variant of the Schwarz alternating method for the case of many subregions
1314 Courant Institute, New York University Technical report
1315 - 2. - Barry Smith, Petter Bjorstad, and William Gropp, Domain Decompositions: Parallel Multilevel Methods for Elliptic Partial Differential Equations,
1316 Cambridge University Press.
1317
1318 .seealso: PCCreate(), PCSetType(), PCType (for list of available types), PC,
1319 PCBJACOBI, PCASMGetSubKSP(), PCASMSetLocalSubdomains(), PCASMType, PCASMGetType(), PCASMSetLocalType(), PCASMGetLocalType()
1320 PCASMSetTotalSubdomains(), PCSetModifySubMatrices(), PCASMSetOverlap(), PCASMSetType(), PCCompositeType
1321
1322 M*/
1323
PCCreate_ASM(PC pc)1324 PETSC_EXTERN PetscErrorCode PCCreate_ASM(PC pc)
1325 {
1326 PetscErrorCode ierr;
1327 PC_ASM *osm;
1328
1329 PetscFunctionBegin;
1330 ierr = PetscNewLog(pc,&osm);CHKERRQ(ierr);
1331
1332 osm->n = PETSC_DECIDE;
1333 osm->n_local = 0;
1334 osm->n_local_true = PETSC_DECIDE;
1335 osm->overlap = 1;
1336 osm->ksp = NULL;
1337 osm->restriction = NULL;
1338 osm->lprolongation = NULL;
1339 osm->lrestriction = NULL;
1340 osm->x = NULL;
1341 osm->y = NULL;
1342 osm->is = NULL;
1343 osm->is_local = NULL;
1344 osm->mat = NULL;
1345 osm->pmat = NULL;
1346 osm->type = PC_ASM_RESTRICT;
1347 osm->loctype = PC_COMPOSITE_ADDITIVE;
1348 osm->same_local_solves = PETSC_TRUE;
1349 osm->sort_indices = PETSC_TRUE;
1350 osm->dm_subdomains = PETSC_FALSE;
1351 osm->sub_mat_type = NULL;
1352
1353 pc->data = (void*)osm;
1354 pc->ops->apply = PCApply_ASM;
1355 pc->ops->matapply = PCMatApply_ASM;
1356 pc->ops->applytranspose = PCApplyTranspose_ASM;
1357 pc->ops->setup = PCSetUp_ASM;
1358 pc->ops->reset = PCReset_ASM;
1359 pc->ops->destroy = PCDestroy_ASM;
1360 pc->ops->setfromoptions = PCSetFromOptions_ASM;
1361 pc->ops->setuponblocks = PCSetUpOnBlocks_ASM;
1362 pc->ops->view = PCView_ASM;
1363 pc->ops->applyrichardson = NULL;
1364
1365 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCASMSetLocalSubdomains_C",PCASMSetLocalSubdomains_ASM);CHKERRQ(ierr);
1366 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCASMSetTotalSubdomains_C",PCASMSetTotalSubdomains_ASM);CHKERRQ(ierr);
1367 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCASMSetOverlap_C",PCASMSetOverlap_ASM);CHKERRQ(ierr);
1368 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCASMSetType_C",PCASMSetType_ASM);CHKERRQ(ierr);
1369 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCASMGetType_C",PCASMGetType_ASM);CHKERRQ(ierr);
1370 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCASMSetLocalType_C",PCASMSetLocalType_ASM);CHKERRQ(ierr);
1371 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCASMGetLocalType_C",PCASMGetLocalType_ASM);CHKERRQ(ierr);
1372 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCASMSetSortIndices_C",PCASMSetSortIndices_ASM);CHKERRQ(ierr);
1373 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCASMGetSubKSP_C",PCASMGetSubKSP_ASM);CHKERRQ(ierr);
1374 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCASMGetSubMatType_C",PCASMGetSubMatType_ASM);CHKERRQ(ierr);
1375 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCASMSetSubMatType_C",PCASMSetSubMatType_ASM);CHKERRQ(ierr);
1376 PetscFunctionReturn(0);
1377 }
1378
1379 /*@C
1380 PCASMCreateSubdomains - Creates the index sets for the overlapping Schwarz
1381 preconditioner for a any problem on a general grid.
1382
1383 Collective
1384
1385 Input Parameters:
1386 + A - The global matrix operator
1387 - n - the number of local blocks
1388
1389 Output Parameters:
1390 . outis - the array of index sets defining the subdomains
1391
1392 Level: advanced
1393
1394 Note: this generates nonoverlapping subdomains; the PCASM will generate the overlap
1395 from these if you use PCASMSetLocalSubdomains()
1396
1397 In the Fortran version you must provide the array outis[] already allocated of length n.
1398
1399 .seealso: PCASMSetLocalSubdomains(), PCASMDestroySubdomains()
1400 @*/
PCASMCreateSubdomains(Mat A,PetscInt n,IS * outis[])1401 PetscErrorCode PCASMCreateSubdomains(Mat A, PetscInt n, IS* outis[])
1402 {
1403 MatPartitioning mpart;
1404 const char *prefix;
1405 PetscInt i,j,rstart,rend,bs;
1406 PetscBool hasop, isbaij = PETSC_FALSE,foundpart = PETSC_FALSE;
1407 Mat Ad = NULL, adj;
1408 IS ispart,isnumb,*is;
1409 PetscErrorCode ierr;
1410
1411 PetscFunctionBegin;
1412 PetscValidHeaderSpecific(A,MAT_CLASSID,1);
1413 PetscValidPointer(outis,3);
1414 if (n < 1) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"number of local blocks must be > 0, n = %D",n);
1415
1416 /* Get prefix, row distribution, and block size */
1417 ierr = MatGetOptionsPrefix(A,&prefix);CHKERRQ(ierr);
1418 ierr = MatGetOwnershipRange(A,&rstart,&rend);CHKERRQ(ierr);
1419 ierr = MatGetBlockSize(A,&bs);CHKERRQ(ierr);
1420 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);
1421
1422 /* Get diagonal block from matrix if possible */
1423 ierr = MatHasOperation(A,MATOP_GET_DIAGONAL_BLOCK,&hasop);CHKERRQ(ierr);
1424 if (hasop) {
1425 ierr = MatGetDiagonalBlock(A,&Ad);CHKERRQ(ierr);
1426 }
1427 if (Ad) {
1428 ierr = PetscObjectBaseTypeCompare((PetscObject)Ad,MATSEQBAIJ,&isbaij);CHKERRQ(ierr);
1429 if (!isbaij) {ierr = PetscObjectBaseTypeCompare((PetscObject)Ad,MATSEQSBAIJ,&isbaij);CHKERRQ(ierr);}
1430 }
1431 if (Ad && n > 1) {
1432 PetscBool match,done;
1433 /* Try to setup a good matrix partitioning if available */
1434 ierr = MatPartitioningCreate(PETSC_COMM_SELF,&mpart);CHKERRQ(ierr);
1435 ierr = PetscObjectSetOptionsPrefix((PetscObject)mpart,prefix);CHKERRQ(ierr);
1436 ierr = MatPartitioningSetFromOptions(mpart);CHKERRQ(ierr);
1437 ierr = PetscObjectTypeCompare((PetscObject)mpart,MATPARTITIONINGCURRENT,&match);CHKERRQ(ierr);
1438 if (!match) {
1439 ierr = PetscObjectTypeCompare((PetscObject)mpart,MATPARTITIONINGSQUARE,&match);CHKERRQ(ierr);
1440 }
1441 if (!match) { /* assume a "good" partitioner is available */
1442 PetscInt na;
1443 const PetscInt *ia,*ja;
1444 ierr = MatGetRowIJ(Ad,0,PETSC_TRUE,isbaij,&na,&ia,&ja,&done);CHKERRQ(ierr);
1445 if (done) {
1446 /* Build adjacency matrix by hand. Unfortunately a call to
1447 MatConvert(Ad,MATMPIADJ,MAT_INITIAL_MATRIX,&adj) will
1448 remove the block-aij structure and we cannot expect
1449 MatPartitioning to split vertices as we need */
1450 PetscInt i,j,len,nnz,cnt,*iia=NULL,*jja=NULL;
1451 const PetscInt *row;
1452 nnz = 0;
1453 for (i=0; i<na; i++) { /* count number of nonzeros */
1454 len = ia[i+1] - ia[i];
1455 row = ja + ia[i];
1456 for (j=0; j<len; j++) {
1457 if (row[j] == i) { /* don't count diagonal */
1458 len--; break;
1459 }
1460 }
1461 nnz += len;
1462 }
1463 ierr = PetscMalloc1(na+1,&iia);CHKERRQ(ierr);
1464 ierr = PetscMalloc1(nnz,&jja);CHKERRQ(ierr);
1465 nnz = 0;
1466 iia[0] = 0;
1467 for (i=0; i<na; i++) { /* fill adjacency */
1468 cnt = 0;
1469 len = ia[i+1] - ia[i];
1470 row = ja + ia[i];
1471 for (j=0; j<len; j++) {
1472 if (row[j] != i) { /* if not diagonal */
1473 jja[nnz+cnt++] = row[j];
1474 }
1475 }
1476 nnz += cnt;
1477 iia[i+1] = nnz;
1478 }
1479 /* Partitioning of the adjacency matrix */
1480 ierr = MatCreateMPIAdj(PETSC_COMM_SELF,na,na,iia,jja,NULL,&adj);CHKERRQ(ierr);
1481 ierr = MatPartitioningSetAdjacency(mpart,adj);CHKERRQ(ierr);
1482 ierr = MatPartitioningSetNParts(mpart,n);CHKERRQ(ierr);
1483 ierr = MatPartitioningApply(mpart,&ispart);CHKERRQ(ierr);
1484 ierr = ISPartitioningToNumbering(ispart,&isnumb);CHKERRQ(ierr);
1485 ierr = MatDestroy(&adj);CHKERRQ(ierr);
1486 foundpart = PETSC_TRUE;
1487 }
1488 ierr = MatRestoreRowIJ(Ad,0,PETSC_TRUE,isbaij,&na,&ia,&ja,&done);CHKERRQ(ierr);
1489 }
1490 ierr = MatPartitioningDestroy(&mpart);CHKERRQ(ierr);
1491 }
1492
1493 ierr = PetscMalloc1(n,&is);CHKERRQ(ierr);
1494 *outis = is;
1495
1496 if (!foundpart) {
1497
1498 /* Partitioning by contiguous chunks of rows */
1499
1500 PetscInt mbs = (rend-rstart)/bs;
1501 PetscInt start = rstart;
1502 for (i=0; i<n; i++) {
1503 PetscInt count = (mbs/n + ((mbs % n) > i)) * bs;
1504 ierr = ISCreateStride(PETSC_COMM_SELF,count,start,1,&is[i]);CHKERRQ(ierr);
1505 start += count;
1506 }
1507
1508 } else {
1509
1510 /* Partitioning by adjacency of diagonal block */
1511
1512 const PetscInt *numbering;
1513 PetscInt *count,nidx,*indices,*newidx,start=0;
1514 /* Get node count in each partition */
1515 ierr = PetscMalloc1(n,&count);CHKERRQ(ierr);
1516 ierr = ISPartitioningCount(ispart,n,count);CHKERRQ(ierr);
1517 if (isbaij && bs > 1) { /* adjust for the block-aij case */
1518 for (i=0; i<n; i++) count[i] *= bs;
1519 }
1520 /* Build indices from node numbering */
1521 ierr = ISGetLocalSize(isnumb,&nidx);CHKERRQ(ierr);
1522 ierr = PetscMalloc1(nidx,&indices);CHKERRQ(ierr);
1523 for (i=0; i<nidx; i++) indices[i] = i; /* needs to be initialized */
1524 ierr = ISGetIndices(isnumb,&numbering);CHKERRQ(ierr);
1525 ierr = PetscSortIntWithPermutation(nidx,numbering,indices);CHKERRQ(ierr);
1526 ierr = ISRestoreIndices(isnumb,&numbering);CHKERRQ(ierr);
1527 if (isbaij && bs > 1) { /* adjust for the block-aij case */
1528 ierr = PetscMalloc1(nidx*bs,&newidx);CHKERRQ(ierr);
1529 for (i=0; i<nidx; i++) {
1530 for (j=0; j<bs; j++) newidx[i*bs+j] = indices[i]*bs + j;
1531 }
1532 ierr = PetscFree(indices);CHKERRQ(ierr);
1533 nidx *= bs;
1534 indices = newidx;
1535 }
1536 /* Shift to get global indices */
1537 for (i=0; i<nidx; i++) indices[i] += rstart;
1538
1539 /* Build the index sets for each block */
1540 for (i=0; i<n; i++) {
1541 ierr = ISCreateGeneral(PETSC_COMM_SELF,count[i],&indices[start],PETSC_COPY_VALUES,&is[i]);CHKERRQ(ierr);
1542 ierr = ISSort(is[i]);CHKERRQ(ierr);
1543 start += count[i];
1544 }
1545
1546 ierr = PetscFree(count);CHKERRQ(ierr);
1547 ierr = PetscFree(indices);CHKERRQ(ierr);
1548 ierr = ISDestroy(&isnumb);CHKERRQ(ierr);
1549 ierr = ISDestroy(&ispart);CHKERRQ(ierr);
1550
1551 }
1552 PetscFunctionReturn(0);
1553 }
1554
1555 /*@C
1556 PCASMDestroySubdomains - Destroys the index sets created with
1557 PCASMCreateSubdomains(). Should be called after setting subdomains
1558 with PCASMSetLocalSubdomains().
1559
1560 Collective
1561
1562 Input Parameters:
1563 + n - the number of index sets
1564 . is - the array of index sets
1565 - is_local - the array of local index sets, can be NULL
1566
1567 Level: advanced
1568
1569 .seealso: PCASMCreateSubdomains(), PCASMSetLocalSubdomains()
1570 @*/
PCASMDestroySubdomains(PetscInt n,IS is[],IS is_local[])1571 PetscErrorCode PCASMDestroySubdomains(PetscInt n, IS is[], IS is_local[])
1572 {
1573 PetscInt i;
1574 PetscErrorCode ierr;
1575
1576 PetscFunctionBegin;
1577 if (n <= 0) PetscFunctionReturn(0);
1578 if (is) {
1579 PetscValidPointer(is,2);
1580 for (i=0; i<n; i++) { ierr = ISDestroy(&is[i]);CHKERRQ(ierr); }
1581 ierr = PetscFree(is);CHKERRQ(ierr);
1582 }
1583 if (is_local) {
1584 PetscValidPointer(is_local,3);
1585 for (i=0; i<n; i++) { ierr = ISDestroy(&is_local[i]);CHKERRQ(ierr); }
1586 ierr = PetscFree(is_local);CHKERRQ(ierr);
1587 }
1588 PetscFunctionReturn(0);
1589 }
1590
1591 /*@
1592 PCASMCreateSubdomains2D - Creates the index sets for the overlapping Schwarz
1593 preconditioner for a two-dimensional problem on a regular grid.
1594
1595 Not Collective
1596
1597 Input Parameters:
1598 + m, n - the number of mesh points in the x and y directions
1599 . M, N - the number of subdomains in the x and y directions
1600 . dof - degrees of freedom per node
1601 - overlap - overlap in mesh lines
1602
1603 Output Parameters:
1604 + Nsub - the number of subdomains created
1605 . is - array of index sets defining overlapping (if overlap > 0) subdomains
1606 - is_local - array of index sets defining non-overlapping subdomains
1607
1608 Note:
1609 Presently PCAMSCreateSubdomains2d() is valid only for sequential
1610 preconditioners. More general related routines are
1611 PCASMSetTotalSubdomains() and PCASMSetLocalSubdomains().
1612
1613 Level: advanced
1614
1615 .seealso: PCASMSetTotalSubdomains(), PCASMSetLocalSubdomains(), PCASMGetSubKSP(),
1616 PCASMSetOverlap()
1617 @*/
PCASMCreateSubdomains2D(PetscInt m,PetscInt n,PetscInt M,PetscInt N,PetscInt dof,PetscInt overlap,PetscInt * Nsub,IS ** is,IS ** is_local)1618 PetscErrorCode PCASMCreateSubdomains2D(PetscInt m,PetscInt n,PetscInt M,PetscInt N,PetscInt dof,PetscInt overlap,PetscInt *Nsub,IS **is,IS **is_local)
1619 {
1620 PetscInt i,j,height,width,ystart,xstart,yleft,yright,xleft,xright,loc_outer;
1621 PetscErrorCode ierr;
1622 PetscInt nidx,*idx,loc,ii,jj,count;
1623
1624 PetscFunctionBegin;
1625 if (dof != 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP," ");
1626
1627 *Nsub = N*M;
1628 ierr = PetscMalloc1(*Nsub,is);CHKERRQ(ierr);
1629 ierr = PetscMalloc1(*Nsub,is_local);CHKERRQ(ierr);
1630 ystart = 0;
1631 loc_outer = 0;
1632 for (i=0; i<N; i++) {
1633 height = n/N + ((n % N) > i); /* height of subdomain */
1634 if (height < 2) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Too many N subdomains for mesh dimension n");
1635 yleft = ystart - overlap; if (yleft < 0) yleft = 0;
1636 yright = ystart + height + overlap; if (yright > n) yright = n;
1637 xstart = 0;
1638 for (j=0; j<M; j++) {
1639 width = m/M + ((m % M) > j); /* width of subdomain */
1640 if (width < 2) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Too many M subdomains for mesh dimension m");
1641 xleft = xstart - overlap; if (xleft < 0) xleft = 0;
1642 xright = xstart + width + overlap; if (xright > m) xright = m;
1643 nidx = (xright - xleft)*(yright - yleft);
1644 ierr = PetscMalloc1(nidx,&idx);CHKERRQ(ierr);
1645 loc = 0;
1646 for (ii=yleft; ii<yright; ii++) {
1647 count = m*ii + xleft;
1648 for (jj=xleft; jj<xright; jj++) idx[loc++] = count++;
1649 }
1650 ierr = ISCreateGeneral(PETSC_COMM_SELF,nidx,idx,PETSC_COPY_VALUES,(*is)+loc_outer);CHKERRQ(ierr);
1651 if (overlap == 0) {
1652 ierr = PetscObjectReference((PetscObject)(*is)[loc_outer]);CHKERRQ(ierr);
1653
1654 (*is_local)[loc_outer] = (*is)[loc_outer];
1655 } else {
1656 for (loc=0,ii=ystart; ii<ystart+height; ii++) {
1657 for (jj=xstart; jj<xstart+width; jj++) {
1658 idx[loc++] = m*ii + jj;
1659 }
1660 }
1661 ierr = ISCreateGeneral(PETSC_COMM_SELF,loc,idx,PETSC_COPY_VALUES,*is_local+loc_outer);CHKERRQ(ierr);
1662 }
1663 ierr = PetscFree(idx);CHKERRQ(ierr);
1664 xstart += width;
1665 loc_outer++;
1666 }
1667 ystart += height;
1668 }
1669 for (i=0; i<*Nsub; i++) { ierr = ISSort((*is)[i]);CHKERRQ(ierr); }
1670 PetscFunctionReturn(0);
1671 }
1672
1673 /*@C
1674 PCASMGetLocalSubdomains - Gets the local subdomains (for this processor
1675 only) for the additive Schwarz preconditioner.
1676
1677 Not Collective
1678
1679 Input Parameter:
1680 . pc - the preconditioner context
1681
1682 Output Parameters:
1683 + n - the number of subdomains for this processor (default value = 1)
1684 . is - the index sets that define the subdomains for this processor
1685 - is_local - the index sets that define the local part of the subdomains for this processor (can be NULL)
1686
1687
1688 Notes:
1689 The IS numbering is in the parallel, global numbering of the vector.
1690
1691 Level: advanced
1692
1693 .seealso: PCASMSetTotalSubdomains(), PCASMSetOverlap(), PCASMGetSubKSP(),
1694 PCASMCreateSubdomains2D(), PCASMSetLocalSubdomains(), PCASMGetLocalSubmatrices()
1695 @*/
PCASMGetLocalSubdomains(PC pc,PetscInt * n,IS * is[],IS * is_local[])1696 PetscErrorCode PCASMGetLocalSubdomains(PC pc,PetscInt *n,IS *is[],IS *is_local[])
1697 {
1698 PC_ASM *osm = (PC_ASM*)pc->data;
1699 PetscErrorCode ierr;
1700 PetscBool match;
1701
1702 PetscFunctionBegin;
1703 PetscValidHeaderSpecific(pc,PC_CLASSID,1);
1704 PetscValidIntPointer(n,2);
1705 if (is) PetscValidPointer(is,3);
1706 ierr = PetscObjectTypeCompare((PetscObject)pc,PCASM,&match);CHKERRQ(ierr);
1707 if (!match) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONG,"PC is not a PCASM");
1708 if (n) *n = osm->n_local_true;
1709 if (is) *is = osm->is;
1710 if (is_local) *is_local = osm->is_local;
1711 PetscFunctionReturn(0);
1712 }
1713
1714 /*@C
1715 PCASMGetLocalSubmatrices - Gets the local submatrices (for this processor
1716 only) for the additive Schwarz preconditioner.
1717
1718 Not Collective
1719
1720 Input Parameter:
1721 . pc - the preconditioner context
1722
1723 Output Parameters:
1724 + n - the number of matrices for this processor (default value = 1)
1725 - mat - the matrices
1726
1727 Level: advanced
1728
1729 Notes:
1730 Call after PCSetUp() (or KSPSetUp()) but before PCApply() and before PCSetUpOnBlocks())
1731
1732 Usually one would use PCSetModifySubMatrices() to change the submatrices in building the preconditioner.
1733
1734 .seealso: PCASMSetTotalSubdomains(), PCASMSetOverlap(), PCASMGetSubKSP(),
1735 PCASMCreateSubdomains2D(), PCASMSetLocalSubdomains(), PCASMGetLocalSubdomains(), PCSetModifySubMatrices()
1736 @*/
PCASMGetLocalSubmatrices(PC pc,PetscInt * n,Mat * mat[])1737 PetscErrorCode PCASMGetLocalSubmatrices(PC pc,PetscInt *n,Mat *mat[])
1738 {
1739 PC_ASM *osm;
1740 PetscErrorCode ierr;
1741 PetscBool match;
1742
1743 PetscFunctionBegin;
1744 PetscValidHeaderSpecific(pc,PC_CLASSID,1);
1745 PetscValidIntPointer(n,2);
1746 if (mat) PetscValidPointer(mat,3);
1747 if (!pc->setupcalled) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONGSTATE,"Must call after KSPSetUp() or PCSetUp().");
1748 ierr = PetscObjectTypeCompare((PetscObject)pc,PCASM,&match);CHKERRQ(ierr);
1749 if (!match) {
1750 if (n) *n = 0;
1751 if (mat) *mat = NULL;
1752 } else {
1753 osm = (PC_ASM*)pc->data;
1754 if (n) *n = osm->n_local_true;
1755 if (mat) *mat = osm->pmat;
1756 }
1757 PetscFunctionReturn(0);
1758 }
1759
1760 /*@
1761 PCASMSetDMSubdomains - Indicates whether to use DMCreateDomainDecomposition() to define the subdomains, whenever possible.
1762
1763 Logically Collective
1764
1765 Input Parameter:
1766 + pc - the preconditioner
1767 - flg - boolean indicating whether to use subdomains defined by the DM
1768
1769 Options Database Key:
1770 . -pc_asm_dm_subdomains
1771
1772 Level: intermediate
1773
1774 Notes:
1775 PCASMSetTotalSubdomains() and PCASMSetOverlap() take precedence over PCASMSetDMSubdomains(),
1776 so setting either of the first two effectively turns the latter off.
1777
1778 .seealso: PCASMGetDMSubdomains(), PCASMSetTotalSubdomains(), PCASMSetOverlap()
1779 PCASMCreateSubdomains2D(), PCASMSetLocalSubdomains(), PCASMGetLocalSubdomains()
1780 @*/
PCASMSetDMSubdomains(PC pc,PetscBool flg)1781 PetscErrorCode PCASMSetDMSubdomains(PC pc,PetscBool flg)
1782 {
1783 PC_ASM *osm = (PC_ASM*)pc->data;
1784 PetscErrorCode ierr;
1785 PetscBool match;
1786
1787 PetscFunctionBegin;
1788 PetscValidHeaderSpecific(pc,PC_CLASSID,1);
1789 PetscValidLogicalCollectiveBool(pc,flg,2);
1790 if (pc->setupcalled) SETERRQ(((PetscObject)pc)->comm,PETSC_ERR_ARG_WRONGSTATE,"Not for a setup PC.");
1791 ierr = PetscObjectTypeCompare((PetscObject)pc,PCASM,&match);CHKERRQ(ierr);
1792 if (match) {
1793 osm->dm_subdomains = flg;
1794 }
1795 PetscFunctionReturn(0);
1796 }
1797
1798 /*@
1799 PCASMGetDMSubdomains - Returns flag indicating whether to use DMCreateDomainDecomposition() to define the subdomains, whenever possible.
1800 Not Collective
1801
1802 Input Parameter:
1803 . pc - the preconditioner
1804
1805 Output Parameter:
1806 . flg - boolean indicating whether to use subdomains defined by the DM
1807
1808 Level: intermediate
1809
1810 .seealso: PCASMSetDMSubdomains(), PCASMSetTotalSubdomains(), PCASMSetOverlap()
1811 PCASMCreateSubdomains2D(), PCASMSetLocalSubdomains(), PCASMGetLocalSubdomains()
1812 @*/
PCASMGetDMSubdomains(PC pc,PetscBool * flg)1813 PetscErrorCode PCASMGetDMSubdomains(PC pc,PetscBool* flg)
1814 {
1815 PC_ASM *osm = (PC_ASM*)pc->data;
1816 PetscErrorCode ierr;
1817 PetscBool match;
1818
1819 PetscFunctionBegin;
1820 PetscValidHeaderSpecific(pc,PC_CLASSID,1);
1821 PetscValidBoolPointer(flg,2);
1822 ierr = PetscObjectTypeCompare((PetscObject)pc,PCASM,&match);CHKERRQ(ierr);
1823 if (match) *flg = osm->dm_subdomains;
1824 else *flg = PETSC_FALSE;
1825 PetscFunctionReturn(0);
1826 }
1827
1828 /*@
1829 PCASMGetSubMatType - Gets the matrix type used for ASM subsolves, as a string.
1830
1831 Not Collective
1832
1833 Input Parameter:
1834 . pc - the PC
1835
1836 Output Parameter:
1837 . -pc_asm_sub_mat_type - name of matrix type
1838
1839 Level: advanced
1840
1841 .seealso: PCASMSetSubMatType(), PCASM, PCSetType(), VecSetType(), MatType, Mat
1842 @*/
PCASMGetSubMatType(PC pc,MatType * sub_mat_type)1843 PetscErrorCode PCASMGetSubMatType(PC pc,MatType *sub_mat_type)
1844 {
1845 PetscErrorCode ierr;
1846
1847 PetscFunctionBegin;
1848 PetscValidHeaderSpecific(pc,PC_CLASSID,1);
1849 ierr = PetscTryMethod(pc,"PCASMGetSubMatType_C",(PC,MatType*),(pc,sub_mat_type));CHKERRQ(ierr);
1850 PetscFunctionReturn(0);
1851 }
1852
1853 /*@
1854 PCASMSetSubMatType - Set the type of matrix used for ASM subsolves
1855
1856 Collective on Mat
1857
1858 Input Parameters:
1859 + pc - the PC object
1860 - sub_mat_type - matrix type
1861
1862 Options Database Key:
1863 . -pc_asm_sub_mat_type <sub_mat_type> - Sets the matrix type used for subsolves, for example, seqaijviennacl. If you specify a base name like aijviennacl, the corresponding sequential type is assumed.
1864
1865 Notes:
1866 See "${PETSC_DIR}/include/petscmat.h" for available types
1867
1868 Level: advanced
1869
1870 .seealso: PCASMGetSubMatType(), PCASM, PCSetType(), VecSetType(), MatType, Mat
1871 @*/
PCASMSetSubMatType(PC pc,MatType sub_mat_type)1872 PetscErrorCode PCASMSetSubMatType(PC pc,MatType sub_mat_type)
1873 {
1874 PetscErrorCode ierr;
1875
1876 PetscFunctionBegin;
1877 PetscValidHeaderSpecific(pc,PC_CLASSID,1);
1878 ierr = PetscTryMethod(pc,"PCASMSetSubMatType_C",(PC,MatType),(pc,sub_mat_type));CHKERRQ(ierr);
1879 PetscFunctionReturn(0);
1880 }
1881