1 /*
2 This file contains some simple default routines.
3 These routines should be SHORT, since they will be included in every
4 executable image that uses the iterative routines (note that, through
5 the registry system, we provide a way to load only the truly necessary
6 files)
7 */
8 #include <petsc/private/kspimpl.h> /*I "petscksp.h" I*/
9 #include <petscdmshell.h>
10
11 /*@
12 KSPGetResidualNorm - Gets the last (approximate preconditioned)
13 residual norm that has been computed.
14
15 Not Collective
16
17 Input Parameters:
18 . ksp - the iterative context
19
20 Output Parameters:
21 . rnorm - residual norm
22
23 Level: intermediate
24
25 .seealso: KSPBuildResidual()
26 @*/
KSPGetResidualNorm(KSP ksp,PetscReal * rnorm)27 PetscErrorCode KSPGetResidualNorm(KSP ksp,PetscReal *rnorm)
28 {
29 PetscFunctionBegin;
30 PetscValidHeaderSpecific(ksp,KSP_CLASSID,1);
31 PetscValidRealPointer(rnorm,2);
32 *rnorm = ksp->rnorm;
33 PetscFunctionReturn(0);
34 }
35
36 /*@
37 KSPGetIterationNumber - Gets the current iteration number; if the
38 KSPSolve() is complete, returns the number of iterations
39 used.
40
41 Not Collective
42
43 Input Parameters:
44 . ksp - the iterative context
45
46 Output Parameters:
47 . its - number of iterations
48
49 Level: intermediate
50
51 Notes:
52 During the ith iteration this returns i-1
53 .seealso: KSPBuildResidual(), KSPGetResidualNorm(), KSPGetTotalIterations()
54 @*/
KSPGetIterationNumber(KSP ksp,PetscInt * its)55 PetscErrorCode KSPGetIterationNumber(KSP ksp,PetscInt *its)
56 {
57 PetscFunctionBegin;
58 PetscValidHeaderSpecific(ksp,KSP_CLASSID,1);
59 PetscValidIntPointer(its,2);
60 *its = ksp->its;
61 PetscFunctionReturn(0);
62 }
63
64 /*@
65 KSPGetTotalIterations - Gets the total number of iterations this KSP object has performed since was created, counted over all linear solves
66
67 Not Collective
68
69 Input Parameters:
70 . ksp - the iterative context
71
72 Output Parameters:
73 . its - total number of iterations
74
75 Level: intermediate
76
77 Notes:
78 Use KSPGetIterationNumber() to get the count for the most recent solve only
79 If this is called within a linear solve (such as in a KSPMonitor routine) then it does not include iterations within that current solve
80
81 .seealso: KSPBuildResidual(), KSPGetResidualNorm(), KSPGetIterationNumber()
82 @*/
KSPGetTotalIterations(KSP ksp,PetscInt * its)83 PetscErrorCode KSPGetTotalIterations(KSP ksp,PetscInt *its)
84 {
85 PetscFunctionBegin;
86 PetscValidHeaderSpecific(ksp,KSP_CLASSID,1);
87 PetscValidIntPointer(its,2);
88 *its = ksp->totalits;
89 PetscFunctionReturn(0);
90 }
91
92 /*@C
93 KSPMonitorSingularValue - Prints the two norm of the true residual and
94 estimation of the extreme singular values of the preconditioned problem
95 at each iteration.
96
97 Logically Collective on ksp
98
99 Input Parameters:
100 + ksp - the iterative context
101 . n - the iteration
102 - rnorm - the two norm of the residual
103
104 Options Database Key:
105 . -ksp_monitor_singular_value - Activates KSPMonitorSingularValue()
106
107 Notes:
108 The CG solver uses the Lanczos technique for eigenvalue computation,
109 while GMRES uses the Arnoldi technique; other iterative methods do
110 not currently compute singular values.
111
112 Level: intermediate
113
114 .seealso: KSPComputeExtremeSingularValues()
115 @*/
KSPMonitorSingularValue(KSP ksp,PetscInt n,PetscReal rnorm,PetscViewerAndFormat * dummy)116 PetscErrorCode KSPMonitorSingularValue(KSP ksp,PetscInt n,PetscReal rnorm,PetscViewerAndFormat *dummy)
117 {
118 PetscReal emin,emax,c;
119 PetscErrorCode ierr;
120 PetscViewer viewer = dummy->viewer;
121
122 PetscFunctionBegin;
123 PetscValidHeaderSpecific(ksp,KSP_CLASSID,1);
124 PetscValidHeaderSpecific(viewer,PETSC_VIEWER_CLASSID,4);
125 ierr = PetscViewerPushFormat(viewer,dummy->format);CHKERRQ(ierr);
126 ierr = PetscViewerASCIIAddTab(viewer,((PetscObject)ksp)->tablevel);CHKERRQ(ierr);
127 if (!ksp->calc_sings) {
128 ierr = PetscViewerASCIIPrintf(viewer,"%3D KSP Residual norm %14.12e \n",n,(double)rnorm);CHKERRQ(ierr);
129 } else {
130 ierr = KSPComputeExtremeSingularValues(ksp,&emax,&emin);CHKERRQ(ierr);
131 c = emax/emin;
132 ierr = PetscViewerASCIIPrintf(viewer,"%3D KSP Residual norm %14.12e %% max %14.12e min %14.12e max/min %14.12e\n",n,(double)rnorm,(double)emax,(double)emin,(double)c);CHKERRQ(ierr);
133 }
134 ierr = PetscViewerASCIISubtractTab(viewer,((PetscObject)ksp)->tablevel);CHKERRQ(ierr);
135 ierr = PetscViewerPopFormat(viewer);CHKERRQ(ierr);
136 PetscFunctionReturn(0);
137 }
138
139 /*@C
140 KSPMonitorSolution - Monitors progress of the KSP solvers by calling
141 VecView() for the approximate solution at each iteration.
142
143 Collective on ksp
144
145 Input Parameters:
146 + ksp - the KSP context
147 . its - iteration number
148 . fgnorm - 2-norm of residual (or gradient)
149 - dummy - a viewer
150
151 Level: intermediate
152
153 Notes:
154 For some Krylov methods such as GMRES constructing the solution at
155 each iteration is expensive, hence using this will slow the code.
156
157 .seealso: KSPMonitorSet(), KSPMonitorDefault(), VecView()
158 @*/
KSPMonitorSolution(KSP ksp,PetscInt its,PetscReal fgnorm,PetscViewerAndFormat * dummy)159 PetscErrorCode KSPMonitorSolution(KSP ksp,PetscInt its,PetscReal fgnorm,PetscViewerAndFormat *dummy)
160 {
161 PetscErrorCode ierr;
162 Vec x;
163 PetscViewer viewer = dummy->viewer;
164
165 PetscFunctionBegin;
166 PetscValidHeaderSpecific(viewer,PETSC_VIEWER_CLASSID,2);
167 ierr = KSPBuildSolution(ksp,NULL,&x);CHKERRQ(ierr);
168 ierr = PetscViewerPushFormat(viewer,dummy->format);CHKERRQ(ierr);
169 ierr = VecView(x,viewer);CHKERRQ(ierr);
170 ierr = PetscViewerPopFormat(viewer);CHKERRQ(ierr);
171 PetscFunctionReturn(0);
172 }
173
174 /*@C
175 KSPMonitorDefault - Print the residual norm at each iteration of an
176 iterative solver.
177
178 Collective on ksp
179
180 Input Parameters:
181 + ksp - iterative context
182 . n - iteration number
183 . rnorm - 2-norm (preconditioned) residual value (may be estimated).
184 - dummy - an ASCII PetscViewer
185
186 Level: intermediate
187
188 .seealso: KSPMonitorSet(), KSPMonitorTrueResidualNorm(), KSPMonitorLGResidualNormCreate()
189 @*/
KSPMonitorDefault(KSP ksp,PetscInt n,PetscReal rnorm,PetscViewerAndFormat * dummy)190 PetscErrorCode KSPMonitorDefault(KSP ksp,PetscInt n,PetscReal rnorm,PetscViewerAndFormat *dummy)
191 {
192 PetscErrorCode ierr;
193 PetscViewer viewer = dummy->viewer;
194
195 PetscFunctionBegin;
196 PetscValidHeaderSpecific(viewer,PETSC_VIEWER_CLASSID,4);
197 ierr = PetscViewerPushFormat(viewer,dummy->format);CHKERRQ(ierr);
198 ierr = PetscViewerASCIIAddTab(viewer,((PetscObject)ksp)->tablevel);CHKERRQ(ierr);
199 if (n == 0 && ((PetscObject)ksp)->prefix) {
200 ierr = PetscViewerASCIIPrintf(viewer," Residual norms for %s solve.\n",((PetscObject)ksp)->prefix);CHKERRQ(ierr);
201 }
202 ierr = PetscViewerASCIIPrintf(viewer,"%3D KSP Residual norm %14.12e \n",n,(double)rnorm);CHKERRQ(ierr);
203 ierr = PetscViewerASCIISubtractTab(viewer,((PetscObject)ksp)->tablevel);CHKERRQ(ierr);
204 ierr = PetscViewerPopFormat(viewer);CHKERRQ(ierr);
205 PetscFunctionReturn(0);
206 }
207
208 /*@C
209 KSPMonitorTrueResidualNorm - Prints the true residual norm as well as the preconditioned
210 residual norm at each iteration of an iterative solver.
211
212 Collective on ksp
213
214 Input Parameters:
215 + ksp - iterative context
216 . n - iteration number
217 . rnorm - 2-norm (preconditioned) residual value (may be estimated).
218 - dummy - an ASCII PetscViewer
219
220 Options Database Key:
221 . -ksp_monitor_true_residual - Activates KSPMonitorTrueResidualNorm()
222
223 Notes:
224 When using right preconditioning, these values are equivalent.
225
226 Level: intermediate
227
228 .seealso: KSPMonitorSet(), KSPMonitorDefault(), KSPMonitorLGResidualNormCreate(),KSPMonitorTrueResidualMaxNorm()
229 @*/
KSPMonitorTrueResidualNorm(KSP ksp,PetscInt n,PetscReal rnorm,PetscViewerAndFormat * dummy)230 PetscErrorCode KSPMonitorTrueResidualNorm(KSP ksp,PetscInt n,PetscReal rnorm,PetscViewerAndFormat *dummy)
231 {
232 PetscErrorCode ierr;
233 Vec resid;
234 PetscReal truenorm,bnorm;
235 PetscViewer viewer = dummy->viewer;
236 char normtype[256];
237
238 PetscFunctionBegin;
239 PetscValidHeaderSpecific(viewer,PETSC_VIEWER_CLASSID,4);
240 ierr = PetscViewerPushFormat(viewer,dummy->format);CHKERRQ(ierr);
241 ierr = PetscViewerASCIIAddTab(viewer,((PetscObject)ksp)->tablevel);CHKERRQ(ierr);
242 if (n == 0 && ((PetscObject)ksp)->prefix) {
243 ierr = PetscViewerASCIIPrintf(viewer," Residual norms for %s solve.\n",((PetscObject)ksp)->prefix);CHKERRQ(ierr);
244 }
245 ierr = KSPBuildResidual(ksp,NULL,NULL,&resid);CHKERRQ(ierr);
246 ierr = VecNorm(resid,NORM_2,&truenorm);CHKERRQ(ierr);
247 ierr = VecDestroy(&resid);CHKERRQ(ierr);
248 ierr = VecNorm(ksp->vec_rhs,NORM_2,&bnorm);CHKERRQ(ierr);
249 ierr = PetscStrncpy(normtype,KSPNormTypes[ksp->normtype],sizeof(normtype));CHKERRQ(ierr);
250 ierr = PetscStrtolower(normtype);CHKERRQ(ierr);
251 ierr = PetscViewerASCIIPrintf(viewer,"%3D KSP %s resid norm %14.12e true resid norm %14.12e ||r(i)||/||b|| %14.12e\n",n,normtype,(double)rnorm,(double)truenorm,(double)(truenorm/bnorm));CHKERRQ(ierr);
252 ierr = PetscViewerASCIISubtractTab(viewer,((PetscObject)ksp)->tablevel);CHKERRQ(ierr);
253 ierr = PetscViewerPopFormat(viewer);CHKERRQ(ierr);
254 PetscFunctionReturn(0);
255 }
256
257 /*@C
258 KSPMonitorTrueResidualMaxNorm - Prints the true residual max norm each iteration of an iterative solver.
259
260 Collective on ksp
261
262 Input Parameters:
263 + ksp - iterative context
264 . n - iteration number
265 . rnorm - norm (preconditioned) residual value (may be estimated).
266 - dummy - an ASCII viewer
267
268 Options Database Key:
269 . -ksp_monitor_max - Activates KSPMonitorTrueResidualMaxNorm()
270
271 Notes:
272 This could be implemented (better) with a flag in ksp.
273
274 Level: intermediate
275
276 .seealso: KSPMonitorSet(), KSPMonitorDefault(), KSPMonitorLGResidualNormCreate(),KSPMonitorTrueResidualNorm()
277 @*/
KSPMonitorTrueResidualMaxNorm(KSP ksp,PetscInt n,PetscReal rnorm,PetscViewerAndFormat * dummy)278 PetscErrorCode KSPMonitorTrueResidualMaxNorm(KSP ksp,PetscInt n,PetscReal rnorm,PetscViewerAndFormat *dummy)
279 {
280 PetscErrorCode ierr;
281 Vec resid;
282 PetscReal truenorm,bnorm;
283 PetscViewer viewer = dummy->viewer;
284 char normtype[256];
285
286 PetscFunctionBegin;
287 PetscValidHeaderSpecific(viewer,PETSC_VIEWER_CLASSID,4);
288 ierr = PetscViewerPushFormat(viewer,dummy->format);CHKERRQ(ierr);
289 ierr = PetscViewerASCIIAddTab(viewer,((PetscObject)ksp)->tablevel);CHKERRQ(ierr);
290 if (n == 0 && ((PetscObject)ksp)->prefix) {
291 ierr = PetscViewerASCIIPrintf(viewer," Residual norms (max) for %s solve.\n",((PetscObject)ksp)->prefix);CHKERRQ(ierr);
292 }
293 ierr = KSPBuildResidual(ksp,NULL,NULL,&resid);CHKERRQ(ierr);
294 ierr = VecNorm(resid,NORM_INFINITY,&truenorm);CHKERRQ(ierr);
295 ierr = VecDestroy(&resid);CHKERRQ(ierr);
296 ierr = VecNorm(ksp->vec_rhs,NORM_INFINITY,&bnorm);CHKERRQ(ierr);
297 ierr = PetscStrncpy(normtype,KSPNormTypes[ksp->normtype],sizeof(normtype));CHKERRQ(ierr);
298 ierr = PetscStrtolower(normtype);CHKERRQ(ierr);
299 ierr = PetscViewerASCIIPrintf(viewer,"%3D KSP true resid max norm %14.12e ||r(i)||/||b|| %14.12e\n",n,(double)truenorm,(double)(truenorm/bnorm));CHKERRQ(ierr);
300 ierr = PetscViewerASCIISubtractTab(viewer,((PetscObject)ksp)->tablevel);CHKERRQ(ierr);
301 ierr = PetscViewerPopFormat(viewer);CHKERRQ(ierr);
302 PetscFunctionReturn(0);
303 }
304
KSPMonitorRange_Private(KSP ksp,PetscInt it,PetscReal * per)305 PetscErrorCode KSPMonitorRange_Private(KSP ksp,PetscInt it,PetscReal *per)
306 {
307 PetscErrorCode ierr;
308 Vec resid;
309 PetscReal rmax,pwork;
310 PetscInt i,n,N;
311 const PetscScalar *r;
312
313 PetscFunctionBegin;
314 ierr = KSPBuildResidual(ksp,NULL,NULL,&resid);CHKERRQ(ierr);
315 ierr = VecNorm(resid,NORM_INFINITY,&rmax);CHKERRQ(ierr);
316 ierr = VecGetLocalSize(resid,&n);CHKERRQ(ierr);
317 ierr = VecGetSize(resid,&N);CHKERRQ(ierr);
318 ierr = VecGetArrayRead(resid,&r);CHKERRQ(ierr);
319 pwork = 0.0;
320 for (i=0; i<n; i++) pwork += (PetscAbsScalar(r[i]) > .20*rmax);
321 ierr = VecRestoreArrayRead(resid,&r);CHKERRQ(ierr);
322 ierr = VecDestroy(&resid);CHKERRQ(ierr);
323 ierr = MPIU_Allreduce(&pwork,per,1,MPIU_REAL,MPIU_SUM,PetscObjectComm((PetscObject)ksp));CHKERRQ(ierr);
324 *per = *per/N;
325 PetscFunctionReturn(0);
326 }
327
328 /*@C
329 KSPMonitorRange - Prints the percentage of residual elements that are more then 10 percent of the maximum value.
330
331 Collective on ksp
332
333 Input Parameters:
334 + ksp - iterative context
335 . it - iteration number
336 . rnorm - 2-norm (preconditioned) residual value (may be estimated).
337 - dummy - an ASCII viewer
338
339 Options Database Key:
340 . -ksp_monitor_range - Activates KSPMonitorRange()
341
342 Level: intermediate
343
344 .seealso: KSPMonitorSet(), KSPMonitorDefault(), KSPMonitorLGResidualNormCreate()
345 @*/
KSPMonitorRange(KSP ksp,PetscInt it,PetscReal rnorm,PetscViewerAndFormat * dummy)346 PetscErrorCode KSPMonitorRange(KSP ksp,PetscInt it,PetscReal rnorm,PetscViewerAndFormat *dummy)
347 {
348 PetscErrorCode ierr;
349 PetscReal perc,rel;
350 PetscViewer viewer = dummy->viewer;
351 /* should be in a MonitorRangeContext */
352 static PetscReal prev;
353
354 PetscFunctionBegin;
355 PetscValidHeaderSpecific(viewer,PETSC_VIEWER_CLASSID,4);
356 ierr = PetscViewerPushFormat(viewer,dummy->format);CHKERRQ(ierr);
357 ierr = PetscViewerASCIIAddTab(viewer,((PetscObject)ksp)->tablevel);CHKERRQ(ierr);
358 if (!it) prev = rnorm;
359 if (it == 0 && ((PetscObject)ksp)->prefix) {
360 ierr = PetscViewerASCIIPrintf(viewer," Residual norms for %s solve.\n",((PetscObject)ksp)->prefix);CHKERRQ(ierr);
361 }
362 ierr = KSPMonitorRange_Private(ksp,it,&perc);CHKERRQ(ierr);
363
364 rel = (prev - rnorm)/prev;
365 prev = rnorm;
366 ierr = PetscViewerASCIIPrintf(viewer,"%3D KSP preconditioned resid norm %14.12e Percent values above 20 percent of maximum %5.2f relative decrease %5.2e ratio %5.2e \n",it,(double)rnorm,(double)(100.0*perc),(double)rel,(double)(rel/perc));CHKERRQ(ierr);
367 ierr = PetscViewerASCIISubtractTab(viewer,((PetscObject)ksp)->tablevel);CHKERRQ(ierr);
368 ierr = PetscViewerPopFormat(viewer);CHKERRQ(ierr);
369 PetscFunctionReturn(0);
370 }
371
372 /*@C
373 KSPMonitorDynamicTolerance - Recompute the inner tolerance in every
374 outer iteration in an adaptive way.
375
376 Collective on ksp
377
378 Input Parameters:
379 + ksp - iterative context
380 . n - iteration number (not used)
381 . fnorm - the current residual norm
382 - dummy - some context as a C struct. fields:
383 coef: a scaling coefficient. default 1.0. can be passed through
384 -sub_ksp_dynamic_tolerance_param
385 bnrm: norm of the right-hand side. store it to avoid repeated calculation
386
387 Notes:
388 This may be useful for a flexibly preconditioner Krylov method to
389 control the accuracy of the inner solves needed to gaurantee the
390 convergence of the outer iterations.
391
392 Level: advanced
393
394 .seealso: KSPMonitorDynamicToleranceDestroy()
395 @*/
KSPMonitorDynamicTolerance(KSP ksp,PetscInt its,PetscReal fnorm,void * dummy)396 PetscErrorCode KSPMonitorDynamicTolerance(KSP ksp,PetscInt its,PetscReal fnorm,void *dummy)
397 {
398 PetscErrorCode ierr;
399 PC pc;
400 PetscReal outer_rtol, outer_abstol, outer_dtol, inner_rtol;
401 PetscInt outer_maxits,nksp,first,i;
402 KSPDynTolCtx *scale = (KSPDynTolCtx*)dummy;
403 KSP *subksp = NULL;
404 KSP kspinner;
405 PetscBool flg;
406
407 PetscFunctionBegin;
408 ierr = KSPGetPC(ksp, &pc);CHKERRQ(ierr);
409
410 /* compute inner_rtol */
411 if (scale->bnrm < 0.0) {
412 Vec b;
413 ierr = KSPGetRhs(ksp, &b);CHKERRQ(ierr);
414 ierr = VecNorm(b, NORM_2, &(scale->bnrm));CHKERRQ(ierr);
415 }
416 ierr = KSPGetTolerances(ksp, &outer_rtol, &outer_abstol, &outer_dtol, &outer_maxits);CHKERRQ(ierr);
417 inner_rtol = PetscMin(scale->coef * scale->bnrm * outer_rtol / fnorm, 0.999);
418 /*ierr = PetscPrintf(PETSC_COMM_WORLD, " Inner rtol = %g\n", (double)inner_rtol);CHKERRQ(ierr);*/
419
420 /* if pc is ksp */
421 ierr = PetscObjectTypeCompare((PetscObject)pc,PCKSP,&flg);CHKERRQ(ierr);
422 if (flg) {
423 ierr = PCKSPGetKSP(pc, &kspinner);CHKERRQ(ierr);
424 ierr = KSPSetTolerances(kspinner, inner_rtol, outer_abstol, outer_dtol, outer_maxits);CHKERRQ(ierr);
425 PetscFunctionReturn(0);
426 }
427
428 /* if pc is bjacobi */
429 ierr = PetscObjectTypeCompare((PetscObject)pc,PCBJACOBI,&flg);CHKERRQ(ierr);
430 if (flg) {
431 ierr = PCBJacobiGetSubKSP(pc, &nksp, &first, &subksp);CHKERRQ(ierr);
432 if (subksp) {
433 for (i=0; i<nksp; i++) {
434 ierr = KSPSetTolerances(subksp[i], inner_rtol, outer_abstol, outer_dtol, outer_maxits);CHKERRQ(ierr);
435 }
436 PetscFunctionReturn(0);
437 }
438 }
439
440 /* if pc is deflation*/
441 ierr = PetscObjectTypeCompare((PetscObject)pc,PCDEFLATION,&flg);CHKERRQ(ierr);
442 if (flg) {
443 ierr = PCDeflationGetCoarseKSP(pc,&kspinner);CHKERRQ(ierr);
444 ierr = KSPSetTolerances(kspinner,inner_rtol,outer_abstol,outer_dtol,PETSC_DEFAULT);CHKERRQ(ierr);
445 PetscFunctionReturn(0);
446 }
447
448 /* todo: dynamic tolerance may apply to other types of pc too */
449 PetscFunctionReturn(0);
450 }
451
452 /*
453 Destroy the dummy context used in KSPMonitorDynamicTolerance()
454 */
KSPMonitorDynamicToleranceDestroy(void ** dummy)455 PetscErrorCode KSPMonitorDynamicToleranceDestroy(void **dummy)
456 {
457 PetscErrorCode ierr;
458
459 PetscFunctionBegin;
460 ierr = PetscFree(*dummy);CHKERRQ(ierr);
461 PetscFunctionReturn(0);
462 }
463
464 /*
465 Default (short) KSP Monitor, same as KSPMonitorDefault() except
466 it prints fewer digits of the residual as the residual gets smaller.
467 This is because the later digits are meaningless and are often
468 different on different machines; by using this routine different
469 machines will usually generate the same output.
470
471 Deprecated: Intentionally has no manual page
472 */
KSPMonitorDefaultShort(KSP ksp,PetscInt its,PetscReal fnorm,PetscViewerAndFormat * dummy)473 PetscErrorCode KSPMonitorDefaultShort(KSP ksp,PetscInt its,PetscReal fnorm,PetscViewerAndFormat *dummy)
474 {
475 PetscErrorCode ierr;
476 PetscViewer viewer = dummy->viewer;
477
478 PetscFunctionBegin;
479 PetscValidHeaderSpecific(viewer,PETSC_VIEWER_CLASSID,4);
480 ierr = PetscViewerPushFormat(viewer,dummy->format);CHKERRQ(ierr);
481 ierr = PetscViewerASCIIAddTab(viewer,((PetscObject)ksp)->tablevel);CHKERRQ(ierr);
482 if (its == 0 && ((PetscObject)ksp)->prefix) {
483 ierr = PetscViewerASCIIPrintf(viewer," Residual norms for %s solve.\n",((PetscObject)ksp)->prefix);CHKERRQ(ierr);
484 }
485
486 if (fnorm > 1.e-9) {
487 ierr = PetscViewerASCIIPrintf(viewer,"%3D KSP Residual norm %g \n",its,(double)fnorm);CHKERRQ(ierr);
488 } else if (fnorm > 1.e-11) {
489 ierr = PetscViewerASCIIPrintf(viewer,"%3D KSP Residual norm %5.3e \n",its,(double)fnorm);CHKERRQ(ierr);
490 } else {
491 ierr = PetscViewerASCIIPrintf(viewer,"%3D KSP Residual norm < 1.e-11\n",its);CHKERRQ(ierr);
492 }
493 ierr = PetscViewerASCIISubtractTab(viewer,((PetscObject)ksp)->tablevel);CHKERRQ(ierr);
494 ierr = PetscViewerPopFormat(viewer);CHKERRQ(ierr);
495 PetscFunctionReturn(0);
496 }
497
498 /*@C
499 KSPConvergedSkip - Convergence test that do not return as converged
500 until the maximum number of iterations is reached.
501
502 Collective on ksp
503
504 Input Parameters:
505 + ksp - iterative context
506 . n - iteration number
507 . rnorm - 2-norm residual value (may be estimated)
508 - dummy - unused convergence context
509
510 Returns:
511 . reason - KSP_CONVERGED_ITERATING, KSP_CONVERGED_ITS
512
513 Notes:
514 This should be used as the convergence test with the option
515 KSPSetNormType(ksp,KSP_NORM_NONE), since norms of the residual are
516 not computed. Convergence is then declared after the maximum number
517 of iterations have been reached. Useful when one is using CG or
518 BiCGStab as a smoother.
519
520 Level: advanced
521
522 .seealso: KSPSetConvergenceTest(), KSPSetTolerances(), KSPSetNormType()
523 @*/
KSPConvergedSkip(KSP ksp,PetscInt n,PetscReal rnorm,KSPConvergedReason * reason,void * dummy)524 PetscErrorCode KSPConvergedSkip(KSP ksp,PetscInt n,PetscReal rnorm,KSPConvergedReason *reason,void *dummy)
525 {
526 PetscFunctionBegin;
527 PetscValidHeaderSpecific(ksp,KSP_CLASSID,1);
528 PetscValidPointer(reason,4);
529 *reason = KSP_CONVERGED_ITERATING;
530 if (n >= ksp->max_it) *reason = KSP_CONVERGED_ITS;
531 PetscFunctionReturn(0);
532 }
533
534
535 /*@C
536 KSPConvergedDefaultCreate - Creates and initializes the space used by the KSPConvergedDefault() function context
537
538 Note Collective
539
540 Output Parameter:
541 . ctx - convergence context
542
543 Level: intermediate
544
545 .seealso: KSPConvergedDefault(), KSPConvergedDefaultDestroy(), KSPSetConvergenceTest(), KSPSetTolerances(),
546 KSPConvergedSkip(), KSPConvergedReason, KSPGetConvergedReason(), KSPConvergedDefaultSetUIRNorm(), KSPConvergedDefaultSetUMIRNorm(), KSPConvergedDefaultSetConvergedMaxits()
547 @*/
KSPConvergedDefaultCreate(void ** ctx)548 PetscErrorCode KSPConvergedDefaultCreate(void **ctx)
549 {
550 PetscErrorCode ierr;
551 KSPConvergedDefaultCtx *cctx;
552
553 PetscFunctionBegin;
554 ierr = PetscNew(&cctx);CHKERRQ(ierr);
555 *ctx = cctx;
556 PetscFunctionReturn(0);
557 }
558
559 /*@
560 KSPConvergedDefaultSetUIRNorm - makes the default convergence test use || B*(b - A*(initial guess))||
561 instead of || B*b ||. In the case of right preconditioner or if KSPSetNormType(ksp,KSP_NORM_UNPRECONDIITONED)
562 is used there is no B in the above formula. UIRNorm is short for Use Initial Residual Norm.
563
564 Collective on ksp
565
566 Input Parameters:
567 . ksp - iterative context
568
569 Options Database:
570 . -ksp_converged_use_initial_residual_norm
571
572 Notes:
573 Use KSPSetTolerances() to alter the defaults for rtol, abstol, dtol.
574
575 The precise values of reason are macros such as KSP_CONVERGED_RTOL, which
576 are defined in petscksp.h.
577
578 If the convergence test is not KSPConvergedDefault() then this is ignored.
579
580 If right preconditioning is being used then B does not appear in the above formula.
581
582
583 Level: intermediate
584
585 .seealso: KSPSetConvergenceTest(), KSPSetTolerances(), KSPConvergedSkip(), KSPConvergedReason, KSPGetConvergedReason(), KSPConvergedDefaultSetUMIRNorm(), KSPConvergedDefaultSetConvergedMaxits()
586 @*/
KSPConvergedDefaultSetUIRNorm(KSP ksp)587 PetscErrorCode KSPConvergedDefaultSetUIRNorm(KSP ksp)
588 {
589 KSPConvergedDefaultCtx *ctx = (KSPConvergedDefaultCtx*) ksp->cnvP;
590
591 PetscFunctionBegin;
592 PetscValidHeaderSpecific(ksp,KSP_CLASSID,1);
593 if (ksp->converged != KSPConvergedDefault) PetscFunctionReturn(0);
594 if (ctx->mininitialrtol) SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_ARG_WRONGSTATE,"Cannot use KSPConvergedDefaultSetUIRNorm() and KSPConvergedDefaultSetUMIRNorm() together");
595 ctx->initialrtol = PETSC_TRUE;
596 PetscFunctionReturn(0);
597 }
598
599 /*@
600 KSPConvergedDefaultSetUMIRNorm - makes the default convergence test use min(|| B*(b - A*(initial guess))||,|| B*b ||)
601 In the case of right preconditioner or if KSPSetNormType(ksp,KSP_NORM_UNPRECONDIITONED)
602 is used there is no B in the above formula. UMIRNorm is short for Use Minimum Initial Residual Norm.
603
604 Collective on ksp
605
606 Input Parameters:
607 . ksp - iterative context
608
609 Options Database:
610 . -ksp_converged_use_min_initial_residual_norm
611
612 Use KSPSetTolerances() to alter the defaults for rtol, abstol, dtol.
613
614 The precise values of reason are macros such as KSP_CONVERGED_RTOL, which
615 are defined in petscksp.h.
616
617 Level: intermediate
618
619 .seealso: KSPSetConvergenceTest(), KSPSetTolerances(), KSPConvergedSkip(), KSPConvergedReason, KSPGetConvergedReason(), KSPConvergedDefaultSetUIRNorm(), KSPConvergedDefaultSetConvergedMaxits()
620 @*/
KSPConvergedDefaultSetUMIRNorm(KSP ksp)621 PetscErrorCode KSPConvergedDefaultSetUMIRNorm(KSP ksp)
622 {
623 KSPConvergedDefaultCtx *ctx = (KSPConvergedDefaultCtx*) ksp->cnvP;
624
625 PetscFunctionBegin;
626 PetscValidHeaderSpecific(ksp,KSP_CLASSID,1);
627 if (ksp->converged != KSPConvergedDefault) PetscFunctionReturn(0);
628 if (ctx->initialrtol) SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_ARG_WRONGSTATE,"Cannot use KSPConvergedDefaultSetUIRNorm() and KSPConvergedDefaultSetUMIRNorm() together");
629 ctx->mininitialrtol = PETSC_TRUE;
630 PetscFunctionReturn(0);
631 }
632
633 /*@
634 KSPConvergedDefaultSetConvergedMaxits - allows the default convergence test to declare convergence and return KSP_CONVERGED_ITS if the maximum number of iterations is reached
635
636 Collective on ksp
637
638 Input Parameters:
639 + ksp - iterative context
640 - flg - boolean flag
641
642 Options Database:
643 . -ksp_converged_maxits
644
645 Use KSPSetTolerances() to alter the defaults for rtol, abstol, dtol.
646
647 The precise values of reason are macros such as KSP_CONVERGED_RTOL, which
648 are defined in petscksp.h.
649
650 Level: intermediate
651
652 .seealso: KSPSetConvergenceTest(), KSPSetTolerances(), KSPConvergedSkip(), KSPConvergedReason, KSPGetConvergedReason(), KSPConvergedDefaultSetUMIRNorm(), KSPConvergedDefaultSetUIRNorm()
653 @*/
KSPConvergedDefaultSetConvergedMaxits(KSP ksp,PetscBool flg)654 PetscErrorCode KSPConvergedDefaultSetConvergedMaxits(KSP ksp, PetscBool flg)
655 {
656 KSPConvergedDefaultCtx *ctx = (KSPConvergedDefaultCtx*) ksp->cnvP;
657
658 PetscFunctionBegin;
659 PetscValidHeaderSpecific(ksp,KSP_CLASSID,1);
660 PetscValidLogicalCollectiveBool(ksp,flg,2);
661 if (ksp->converged != KSPConvergedDefault) PetscFunctionReturn(0);
662 ctx->convmaxits = flg;
663 PetscFunctionReturn(0);
664 }
665
666 /*@C
667 KSPConvergedDefault - Determines convergence of the linear iterative solvers by default
668
669 Collective on ksp
670
671 Input Parameters:
672 + ksp - iterative context
673 . n - iteration number
674 . rnorm - residual norm (may be estimated, depending on the method may be the preconditioned residual norm)
675 - ctx - convergence context which must be created by KSPConvergedDefaultCreate()
676
677 Output Parameter:
678 + positive - if the iteration has converged
679 . negative - if the iteration has diverged
680 - KSP_CONVERGED_ITERATING - otherwise.
681
682 Notes:
683 KSPConvergedDefault() reaches convergence when rnorm < MAX (rtol * rnorm_0, abstol);
684 Divergence is detected if rnorm > dtol * rnorm_0, or when failures are detected throughout the iteration.
685 By default, reaching the maximum number of iterations is considered divergence (i.e. KSP_DIVERGED_ITS).
686 In order to have PETSc declaring convergence in such a case (i.e. KSP_CONVERGED_ITS), users can use KSPConvergedDefaultSetConvergedMaxits()
687
688 where:
689 + rtol = relative tolerance,
690 . abstol = absolute tolerance.
691 . dtol = divergence tolerance,
692 - rnorm_0 is the two norm of the right hand side (or the preconditioned norm, depending on what was set with
693 KSPSetNormType(). When initial guess is non-zero you
694 can call KSPConvergedDefaultSetUIRNorm() to use the norm of (b - A*(initial guess))
695 as the starting point for relative norm convergence testing, that is as rnorm_0
696
697 Use KSPSetTolerances() to alter the defaults for rtol, abstol, dtol.
698
699 Use KSPSetNormType() (or -ksp_norm_type <none,preconditioned,unpreconditioned,natural>) to change the norm used for computing rnorm
700
701 The precise values of reason are macros such as KSP_CONVERGED_RTOL, which are defined in petscksp.h.
702
703 This routine is used by KSP by default so the user generally never needs call it directly.
704
705 Use KSPSetConvergenceTest() to provide your own test instead of using this one.
706
707 Level: intermediate
708
709 .seealso: KSPSetConvergenceTest(), KSPSetTolerances(), KSPConvergedSkip(), KSPConvergedReason, KSPGetConvergedReason(),
710 KSPConvergedDefaultSetUIRNorm(), KSPConvergedDefaultSetUMIRNorm(), KSPConvergedDefaultSetConvergedMaxits(), KSPConvergedDefaultCreate(), KSPConvergedDefaultDestroy()
711 @*/
KSPConvergedDefault(KSP ksp,PetscInt n,PetscReal rnorm,KSPConvergedReason * reason,void * ctx)712 PetscErrorCode KSPConvergedDefault(KSP ksp,PetscInt n,PetscReal rnorm,KSPConvergedReason *reason,void *ctx)
713 {
714 PetscErrorCode ierr;
715 KSPConvergedDefaultCtx *cctx = (KSPConvergedDefaultCtx*) ctx;
716 KSPNormType normtype;
717
718 PetscFunctionBegin;
719 PetscValidHeaderSpecific(ksp,KSP_CLASSID,1);
720 PetscValidLogicalCollectiveInt(ksp,n,2);
721 PetscValidPointer(reason,4);
722 if (!cctx) SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_ARG_NULL,"Convergence context must have been created with KSPConvergedDefaultCreate()");
723 *reason = KSP_CONVERGED_ITERATING;
724
725 if (cctx->convmaxits && n >= ksp->max_it) {
726 *reason = KSP_CONVERGED_ITS;
727 ierr = PetscInfo1(ksp,"Linear solver has converged. Maximum number of iterations reached %D\n",n);CHKERRQ(ierr);
728 PetscFunctionReturn(0);
729 }
730 ierr = KSPGetNormType(ksp,&normtype);CHKERRQ(ierr);
731 if (normtype == KSP_NORM_NONE) PetscFunctionReturn(0);
732
733 if (!n) {
734 /* if user gives initial guess need to compute norm of b */
735 if (!ksp->guess_zero && !cctx->initialrtol) {
736 PetscReal snorm = 0.0;
737 if (ksp->normtype == KSP_NORM_UNPRECONDITIONED || ksp->pc_side == PC_RIGHT) {
738 ierr = PetscInfo(ksp,"user has provided nonzero initial guess, computing 2-norm of RHS\n");CHKERRQ(ierr);
739 ierr = VecNorm(ksp->vec_rhs,NORM_2,&snorm);CHKERRQ(ierr); /* <- b'*b */
740 } else {
741 Vec z;
742 /* Should avoid allocating the z vector each time but cannot stash it in cctx because if KSPReset() is called the vector size might change */
743 ierr = VecDuplicate(ksp->vec_rhs,&z);CHKERRQ(ierr);
744 ierr = KSP_PCApply(ksp,ksp->vec_rhs,z);CHKERRQ(ierr);
745 if (ksp->normtype == KSP_NORM_PRECONDITIONED) {
746 ierr = PetscInfo(ksp,"user has provided nonzero initial guess, computing 2-norm of preconditioned RHS\n");CHKERRQ(ierr);
747 ierr = VecNorm(z,NORM_2,&snorm);CHKERRQ(ierr); /* dp <- b'*B'*B*b */
748 } else if (ksp->normtype == KSP_NORM_NATURAL) {
749 PetscScalar norm;
750 ierr = PetscInfo(ksp,"user has provided nonzero initial guess, computing natural norm of RHS\n");CHKERRQ(ierr);
751 ierr = VecDot(ksp->vec_rhs,z,&norm);CHKERRQ(ierr);
752 snorm = PetscSqrtReal(PetscAbsScalar(norm)); /* dp <- b'*B*b */
753 }
754 ierr = VecDestroy(&z);CHKERRQ(ierr);
755 }
756 /* handle special case of zero RHS and nonzero guess */
757 if (!snorm) {
758 ierr = PetscInfo(ksp,"Special case, user has provided nonzero initial guess and zero RHS\n");CHKERRQ(ierr);
759 snorm = rnorm;
760 }
761 if (cctx->mininitialrtol) ksp->rnorm0 = PetscMin(snorm,rnorm);
762 else ksp->rnorm0 = snorm;
763 } else {
764 ksp->rnorm0 = rnorm;
765 }
766 ksp->ttol = PetscMax(ksp->rtol*ksp->rnorm0,ksp->abstol);
767 }
768
769 if (n <= ksp->chknorm) PetscFunctionReturn(0);
770
771 if (PetscIsInfOrNanReal(rnorm)) {
772 PCFailedReason pcreason;
773 PetscInt sendbuf,recvbuf;
774 ierr = PCGetFailedReasonRank(ksp->pc,&pcreason);CHKERRQ(ierr);
775 sendbuf = (PetscInt)pcreason;
776 ierr = MPI_Allreduce(&sendbuf,&recvbuf,1,MPIU_INT,MPIU_MAX,PetscObjectComm((PetscObject)ksp));CHKERRQ(ierr);
777 if (recvbuf) {
778 *reason = KSP_DIVERGED_PC_FAILED;
779 ierr = PCSetFailedReason(ksp->pc,(PCFailedReason)recvbuf);CHKERRQ(ierr);
780 ierr = PetscInfo(ksp,"Linear solver pcsetup fails, declaring divergence \n");CHKERRQ(ierr);
781 } else {
782 *reason = KSP_DIVERGED_NANORINF;
783 ierr = PetscInfo(ksp,"Linear solver has created a not a number (NaN) as the residual norm, declaring divergence \n");CHKERRQ(ierr);
784 }
785 } else if (rnorm <= ksp->ttol) {
786 if (rnorm < ksp->abstol) {
787 ierr = PetscInfo3(ksp,"Linear solver has converged. Residual norm %14.12e is less than absolute tolerance %14.12e at iteration %D\n",(double)rnorm,(double)ksp->abstol,n);CHKERRQ(ierr);
788 *reason = KSP_CONVERGED_ATOL;
789 } else {
790 if (cctx->initialrtol) {
791 ierr = PetscInfo4(ksp,"Linear solver has converged. Residual norm %14.12e is less than relative tolerance %14.12e times initial residual norm %14.12e at iteration %D\n",(double)rnorm,(double)ksp->rtol,(double)ksp->rnorm0,n);CHKERRQ(ierr);
792 } else {
793 ierr = PetscInfo4(ksp,"Linear solver has converged. Residual norm %14.12e is less than relative tolerance %14.12e times initial right hand side norm %14.12e at iteration %D\n",(double)rnorm,(double)ksp->rtol,(double)ksp->rnorm0,n);CHKERRQ(ierr);
794 }
795 *reason = KSP_CONVERGED_RTOL;
796 }
797 } else if (rnorm >= ksp->divtol*ksp->rnorm0) {
798 ierr = PetscInfo3(ksp,"Linear solver is diverging. Initial right hand size norm %14.12e, current residual norm %14.12e at iteration %D\n",(double)ksp->rnorm0,(double)rnorm,n);CHKERRQ(ierr);
799 *reason = KSP_DIVERGED_DTOL;
800 }
801 PetscFunctionReturn(0);
802 }
803
804 /*@C
805 KSPConvergedDefaultDestroy - Frees the space used by the KSPConvergedDefault() function context
806
807 Not Collective
808
809 Input Parameters:
810 . ctx - convergence context
811
812 Level: intermediate
813
814 .seealso: KSPConvergedDefault(), KSPConvergedDefaultCreate(), KSPSetConvergenceTest(), KSPSetTolerances(), KSPConvergedSkip(),
815 KSPConvergedReason, KSPGetConvergedReason(), KSPConvergedDefaultSetUIRNorm(), KSPConvergedDefaultSetUMIRNorm()
816 @*/
KSPConvergedDefaultDestroy(void * ctx)817 PetscErrorCode KSPConvergedDefaultDestroy(void *ctx)
818 {
819 PetscErrorCode ierr;
820 KSPConvergedDefaultCtx *cctx = (KSPConvergedDefaultCtx*) ctx;
821
822 PetscFunctionBegin;
823 ierr = VecDestroy(&cctx->work);CHKERRQ(ierr);
824 ierr = PetscFree(ctx);CHKERRQ(ierr);
825 PetscFunctionReturn(0);
826 }
827
828 /*
829 KSPBuildSolutionDefault - Default code to create/move the solution.
830
831 Collective on ksp
832
833 Input Parameters:
834 + ksp - iterative context
835 - v - pointer to the user's vector
836
837 Output Parameter:
838 . V - pointer to a vector containing the solution
839
840 Level: advanced
841
842 Developers Note: This is PETSC_EXTERN because it may be used by user written plugin KSP implementations
843
844 .seealso: KSPGetSolution(), KSPBuildResidualDefault()
845 */
KSPBuildSolutionDefault(KSP ksp,Vec v,Vec * V)846 PetscErrorCode KSPBuildSolutionDefault(KSP ksp,Vec v,Vec *V)
847 {
848 PetscErrorCode ierr;
849
850 PetscFunctionBegin;
851 if (ksp->pc_side == PC_RIGHT) {
852 if (ksp->pc) {
853 if (v) {
854 ierr = KSP_PCApply(ksp,ksp->vec_sol,v);CHKERRQ(ierr); *V = v;
855 } else SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_SUP,"Not working with right preconditioner");
856 } else {
857 if (v) {
858 ierr = VecCopy(ksp->vec_sol,v);CHKERRQ(ierr); *V = v;
859 } else *V = ksp->vec_sol;
860 }
861 } else if (ksp->pc_side == PC_SYMMETRIC) {
862 if (ksp->pc) {
863 if (ksp->transpose_solve) SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_SUP,"Not working with symmetric preconditioner and transpose solve");
864 if (v) {
865 ierr = PCApplySymmetricRight(ksp->pc,ksp->vec_sol,v);CHKERRQ(ierr);
866 *V = v;
867 } else SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_SUP,"Not working with symmetric preconditioner");
868 } else {
869 if (v) {
870 ierr = VecCopy(ksp->vec_sol,v);CHKERRQ(ierr); *V = v;
871 } else *V = ksp->vec_sol;
872 }
873 } else {
874 if (v) {
875 ierr = VecCopy(ksp->vec_sol,v);CHKERRQ(ierr); *V = v;
876 } else *V = ksp->vec_sol;
877 }
878 PetscFunctionReturn(0);
879 }
880
881 /*
882 KSPBuildResidualDefault - Default code to compute the residual.
883
884 Collecive on ksp
885
886 Input Parameters:
887 . ksp - iterative context
888 . t - pointer to temporary vector
889 . v - pointer to user vector
890
891 Output Parameter:
892 . V - pointer to a vector containing the residual
893
894 Level: advanced
895
896 Developers Note: This is PETSC_EXTERN because it may be used by user written plugin KSP implementations
897
898 .seealso: KSPBuildSolutionDefault()
899 */
KSPBuildResidualDefault(KSP ksp,Vec t,Vec v,Vec * V)900 PetscErrorCode KSPBuildResidualDefault(KSP ksp,Vec t,Vec v,Vec *V)
901 {
902 PetscErrorCode ierr;
903 Mat Amat,Pmat;
904
905 PetscFunctionBegin;
906 if (!ksp->pc) {ierr = KSPGetPC(ksp,&ksp->pc);CHKERRQ(ierr);}
907 ierr = PCGetOperators(ksp->pc,&Amat,&Pmat);CHKERRQ(ierr);
908 ierr = KSPBuildSolution(ksp,t,NULL);CHKERRQ(ierr);
909 ierr = KSP_MatMult(ksp,Amat,t,v);CHKERRQ(ierr);
910 ierr = VecAYPX(v,-1.0,ksp->vec_rhs);CHKERRQ(ierr);
911 *V = v;
912 PetscFunctionReturn(0);
913 }
914
915 /*@C
916 KSPCreateVecs - Gets a number of work vectors.
917
918 Collective on ksp
919
920 Input Parameters:
921 + ksp - iterative context
922 . rightn - number of right work vectors
923 - leftn - number of left work vectors to allocate
924
925 Output Parameter:
926 + right - the array of vectors created
927 - left - the array of left vectors
928
929 Note: The right vector has as many elements as the matrix has columns. The left
930 vector has as many elements as the matrix has rows.
931
932 The vectors are new vectors that are not owned by the KSP, they should be destroyed with calls to VecDestroyVecs() when no longer needed.
933
934 Developers Note: First tries to duplicate the rhs and solution vectors of the KSP, if they do not exist tries to get them from the matrix, if
935 that does not exist tries to get them from the DM (if it is provided).
936
937 Level: advanced
938
939 .seealso: MatCreateVecs(), VecDestroyVecs()
940
941 @*/
KSPCreateVecs(KSP ksp,PetscInt rightn,Vec ** right,PetscInt leftn,Vec ** left)942 PetscErrorCode KSPCreateVecs(KSP ksp,PetscInt rightn, Vec **right,PetscInt leftn,Vec **left)
943 {
944 PetscErrorCode ierr;
945 Vec vecr = NULL,vecl = NULL;
946 PetscBool matset,pmatset;
947 Mat mat = NULL;
948
949 PetscFunctionBegin;
950 if (rightn) {
951 if (!right) SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_ARG_INCOMP,"You asked for right vectors but did not pass a pointer to hold them");
952 if (ksp->vec_sol) vecr = ksp->vec_sol;
953 else {
954 if (ksp->pc) {
955 ierr = PCGetOperatorsSet(ksp->pc,&matset,&pmatset);CHKERRQ(ierr);
956 /* check for mat before pmat because for KSPLSQR pmat may be a different size than mat since pmat maybe mat'*mat */
957 if (matset) {
958 ierr = PCGetOperators(ksp->pc,&mat,NULL);CHKERRQ(ierr);
959 ierr = MatCreateVecs(mat,&vecr,NULL);CHKERRQ(ierr);
960 } else if (pmatset) {
961 ierr = PCGetOperators(ksp->pc,NULL,&mat);CHKERRQ(ierr);
962 ierr = MatCreateVecs(mat,&vecr,NULL);CHKERRQ(ierr);
963 }
964 }
965 if (!vecr) {
966 if (ksp->dm) {
967 ierr = DMGetGlobalVector(ksp->dm,&vecr);CHKERRQ(ierr);
968 } else SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_ARG_WRONGSTATE,"You requested a vector from a KSP that cannot provide one");
969 }
970 }
971 ierr = VecDuplicateVecs(vecr,rightn,right);CHKERRQ(ierr);
972 if (!ksp->vec_sol) {
973 if (mat) {
974 ierr = VecDestroy(&vecr);CHKERRQ(ierr);
975 } else if (ksp->dm) {
976 ierr = DMRestoreGlobalVector(ksp->dm,&vecr);CHKERRQ(ierr);
977 }
978 }
979 }
980 if (leftn) {
981 if (!left) SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_ARG_INCOMP,"You asked for left vectors but did not pass a pointer to hold them");
982 if (ksp->vec_rhs) vecl = ksp->vec_rhs;
983 else {
984 if (ksp->pc) {
985 ierr = PCGetOperatorsSet(ksp->pc,&matset,&pmatset);CHKERRQ(ierr);
986 /* check for mat before pmat because for KSPLSQR pmat may be a different size than mat since pmat maybe mat'*mat */
987 if (matset) {
988 ierr = PCGetOperators(ksp->pc,&mat,NULL);CHKERRQ(ierr);
989 ierr = MatCreateVecs(mat,NULL,&vecl);CHKERRQ(ierr);
990 } else if (pmatset) {
991 ierr = PCGetOperators(ksp->pc,NULL,&mat);CHKERRQ(ierr);
992 ierr = MatCreateVecs(mat,NULL,&vecl);CHKERRQ(ierr);
993 }
994 }
995 if (!vecl) {
996 if (ksp->dm) {
997 ierr = DMGetGlobalVector(ksp->dm,&vecl);CHKERRQ(ierr);
998 } else SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_ARG_WRONGSTATE,"You requested a vector from a KSP that cannot provide one");
999 }
1000 }
1001 ierr = VecDuplicateVecs(vecl,leftn,left);CHKERRQ(ierr);
1002 if (!ksp->vec_rhs) {
1003 if (mat) {
1004 ierr = VecDestroy(&vecl);CHKERRQ(ierr);
1005 } else if (ksp->dm) {
1006 ierr = DMRestoreGlobalVector(ksp->dm,&vecl);CHKERRQ(ierr);
1007 }
1008 }
1009 }
1010 PetscFunctionReturn(0);
1011 }
1012
1013 /*@C
1014 KSPSetWorkVecs - Sets a number of work vectors into a KSP object
1015
1016 Collective on ksp
1017
1018 Input Parameters:
1019 + ksp - iterative context
1020 - nw - number of work vectors to allocate
1021
1022 Level: developer
1023
1024 Developers Note: This is PETSC_EXTERN because it may be used by user written plugin KSP implementations
1025
1026 @*/
KSPSetWorkVecs(KSP ksp,PetscInt nw)1027 PetscErrorCode KSPSetWorkVecs(KSP ksp,PetscInt nw)
1028 {
1029 PetscErrorCode ierr;
1030
1031 PetscFunctionBegin;
1032 ierr = VecDestroyVecs(ksp->nwork,&ksp->work);CHKERRQ(ierr);
1033 ksp->nwork = nw;
1034 ierr = KSPCreateVecs(ksp,nw,&ksp->work,0,NULL);CHKERRQ(ierr);
1035 ierr = PetscLogObjectParents(ksp,nw,ksp->work);CHKERRQ(ierr);
1036 PetscFunctionReturn(0);
1037 }
1038
1039 /*
1040 KSPDestroyDefault - Destroys a iterative context variable for methods with
1041 no separate context. Preferred calling sequence KSPDestroy().
1042
1043 Input Parameter:
1044 . ksp - the iterative context
1045
1046 Developers Note: This is PETSC_EXTERN because it may be used by user written plugin KSP implementations
1047
1048 */
KSPDestroyDefault(KSP ksp)1049 PetscErrorCode KSPDestroyDefault(KSP ksp)
1050 {
1051 PetscErrorCode ierr;
1052
1053 PetscFunctionBegin;
1054 PetscValidHeaderSpecific(ksp,KSP_CLASSID,1);
1055 ierr = PetscFree(ksp->data);CHKERRQ(ierr);
1056 PetscFunctionReturn(0);
1057 }
1058
1059 /*@
1060 KSPGetConvergedReason - Gets the reason the KSP iteration was stopped.
1061
1062 Not Collective
1063
1064 Input Parameter:
1065 . ksp - the KSP context
1066
1067 Output Parameter:
1068 . reason - negative value indicates diverged, positive value converged, see KSPConvergedReason
1069
1070 Possible values for reason: See also manual page for each reason
1071 $ KSP_CONVERGED_RTOL (residual 2-norm decreased by a factor of rtol, from 2-norm of right hand side)
1072 $ KSP_CONVERGED_ATOL (residual 2-norm less than abstol)
1073 $ KSP_CONVERGED_ITS (used by the preonly preconditioner that always uses ONE iteration, or when the KSPConvergedSkip() convergence test routine is set.
1074 $ KSP_CONVERGED_CG_NEG_CURVE (see note below)
1075 $ KSP_CONVERGED_CG_CONSTRAINED (see note below)
1076 $ KSP_CONVERGED_STEP_LENGTH (see note below)
1077 $ KSP_CONVERGED_ITERATING (returned if the solver is not yet finished)
1078 $ KSP_DIVERGED_ITS (required more than its to reach convergence)
1079 $ KSP_DIVERGED_DTOL (residual norm increased by a factor of divtol)
1080 $ KSP_DIVERGED_NANORINF (residual norm became Not-a-number or Inf likely due to 0/0)
1081 $ KSP_DIVERGED_BREAKDOWN (generic breakdown in method)
1082 $ KSP_DIVERGED_BREAKDOWN_BICG (Initial residual is orthogonal to preconditioned initial residual. Try a different preconditioner, or a different initial Level.)
1083
1084 Options Database:
1085 . -ksp_converged_reason - prints the reason to standard out
1086
1087 Notes:
1088 If this routine is called before or doing the KSPSolve() the value of KSP_CONVERGED_ITERATING is returned
1089
1090 The values KSP_CONVERGED_CG_NEG_CURVE, KSP_CONVERGED_CG_CONSTRAINED, and KSP_CONVERGED_STEP_LENGTH are returned only by the special KSPNASH, KSPSTCG, and KSPGLTR
1091 solvers which are used by the SNESNEWTONTR (trust region) solver.
1092
1093 Level: intermediate
1094
1095 .seealso: KSPSetConvergenceTest(), KSPConvergedDefault(), KSPSetTolerances(), KSPConvergedReason,
1096 KSPConvergedReasonView()
1097 @*/
KSPGetConvergedReason(KSP ksp,KSPConvergedReason * reason)1098 PetscErrorCode KSPGetConvergedReason(KSP ksp,KSPConvergedReason *reason)
1099 {
1100 PetscFunctionBegin;
1101 PetscValidHeaderSpecific(ksp,KSP_CLASSID,1);
1102 PetscValidPointer(reason,2);
1103 *reason = ksp->reason;
1104 PetscFunctionReturn(0);
1105 }
1106
1107 #include <petsc/private/dmimpl.h>
1108 /*@
1109 KSPSetDM - Sets the DM that may be used by some preconditioners
1110
1111 Logically Collective on ksp
1112
1113 Input Parameters:
1114 + ksp - the preconditioner context
1115 - dm - the dm, cannot be NULL
1116
1117 Notes:
1118 If this is used then the KSP will attempt to use the DM to create the matrix and use the routine set with
1119 DMKSPSetComputeOperators(). Use KSPSetDMActive(ksp,PETSC_FALSE) to instead use the matrix you've provided with
1120 KSPSetOperators().
1121
1122 A DM can only be used for solving one problem at a time because information about the problem is stored on the DM,
1123 even when not using interfaces like DMKSPSetComputeOperators(). Use DMClone() to get a distinct DM when solving
1124 different problems using the same function space.
1125
1126 Level: intermediate
1127
1128 .seealso: KSPGetDM(), KSPSetDMActive(), KSPSetComputeOperators(), KSPSetComputeRHS(), KSPSetComputeInitialGuess(), DMKSPSetComputeOperators(), DMKSPSetComputeRHS(), DMKSPSetComputeInitialGuess()
1129 @*/
KSPSetDM(KSP ksp,DM dm)1130 PetscErrorCode KSPSetDM(KSP ksp,DM dm)
1131 {
1132 PetscErrorCode ierr;
1133 PC pc;
1134
1135 PetscFunctionBegin;
1136 PetscValidHeaderSpecific(ksp,KSP_CLASSID,1);
1137 PetscValidHeaderSpecific(dm,DM_CLASSID,2);
1138 ierr = PetscObjectReference((PetscObject)dm);CHKERRQ(ierr);
1139 if (ksp->dm) { /* Move the DMSNES context over to the new DM unless the new DM already has one */
1140 if (ksp->dm->dmksp && !dm->dmksp) {
1141 DMKSP kdm;
1142 ierr = DMCopyDMKSP(ksp->dm,dm);CHKERRQ(ierr);
1143 ierr = DMGetDMKSP(ksp->dm,&kdm);CHKERRQ(ierr);
1144 if (kdm->originaldm == ksp->dm) kdm->originaldm = dm; /* Grant write privileges to the replacement DM */
1145 }
1146 ierr = DMDestroy(&ksp->dm);CHKERRQ(ierr);
1147 }
1148 ksp->dm = dm;
1149 ksp->dmAuto = PETSC_FALSE;
1150 ierr = KSPGetPC(ksp,&pc);CHKERRQ(ierr);
1151 ierr = PCSetDM(pc,dm);CHKERRQ(ierr);
1152 ksp->dmActive = PETSC_TRUE;
1153 PetscFunctionReturn(0);
1154 }
1155
1156 /*@
1157 KSPSetDMActive - Indicates the DM should be used to generate the linear system matrix and right hand side
1158
1159 Logically Collective on ksp
1160
1161 Input Parameters:
1162 + ksp - the preconditioner context
1163 - flg - use the DM
1164
1165 Notes:
1166 By default KSPSetDM() sets the DM as active, call KSPSetDMActive(ksp,PETSC_FALSE); after KSPSetDM(ksp,dm) to not have the KSP object use the DM to generate the matrices.
1167
1168 Level: intermediate
1169
1170 .seealso: KSPGetDM(), KSPSetDM(), SNESSetDM(), KSPSetComputeOperators(), KSPSetComputeRHS(), KSPSetComputeInitialGuess()
1171 @*/
KSPSetDMActive(KSP ksp,PetscBool flg)1172 PetscErrorCode KSPSetDMActive(KSP ksp,PetscBool flg)
1173 {
1174 PetscFunctionBegin;
1175 PetscValidHeaderSpecific(ksp,KSP_CLASSID,1);
1176 PetscValidLogicalCollectiveBool(ksp,flg,2);
1177 ksp->dmActive = flg;
1178 PetscFunctionReturn(0);
1179 }
1180
1181 /*@
1182 KSPGetDM - Gets the DM that may be used by some preconditioners
1183
1184 Not Collective
1185
1186 Input Parameter:
1187 . ksp - the preconditioner context
1188
1189 Output Parameter:
1190 . dm - the dm
1191
1192 Level: intermediate
1193
1194
1195 .seealso: KSPSetDM(), KSPSetDMActive()
1196 @*/
KSPGetDM(KSP ksp,DM * dm)1197 PetscErrorCode KSPGetDM(KSP ksp,DM *dm)
1198 {
1199 PetscErrorCode ierr;
1200
1201 PetscFunctionBegin;
1202 PetscValidHeaderSpecific(ksp,KSP_CLASSID,1);
1203 if (!ksp->dm) {
1204 ierr = DMShellCreate(PetscObjectComm((PetscObject)ksp),&ksp->dm);CHKERRQ(ierr);
1205 ksp->dmAuto = PETSC_TRUE;
1206 }
1207 *dm = ksp->dm;
1208 PetscFunctionReturn(0);
1209 }
1210
1211 /*@
1212 KSPSetApplicationContext - Sets the optional user-defined context for the linear solver.
1213
1214 Logically Collective on ksp
1215
1216 Input Parameters:
1217 + ksp - the KSP context
1218 - usrP - optional user context
1219
1220 Fortran Notes:
1221 To use this from Fortran you must write a Fortran interface definition for this
1222 function that tells Fortran the Fortran derived data type that you are passing in as the ctx argument.
1223
1224 Level: intermediate
1225
1226 .seealso: KSPGetApplicationContext()
1227 @*/
KSPSetApplicationContext(KSP ksp,void * usrP)1228 PetscErrorCode KSPSetApplicationContext(KSP ksp,void *usrP)
1229 {
1230 PetscErrorCode ierr;
1231 PC pc;
1232
1233 PetscFunctionBegin;
1234 PetscValidHeaderSpecific(ksp,KSP_CLASSID,1);
1235 ksp->user = usrP;
1236 ierr = KSPGetPC(ksp,&pc);CHKERRQ(ierr);
1237 ierr = PCSetApplicationContext(pc,usrP);CHKERRQ(ierr);
1238 PetscFunctionReturn(0);
1239 }
1240
1241 /*@
1242 KSPGetApplicationContext - Gets the user-defined context for the linear solver.
1243
1244 Not Collective
1245
1246 Input Parameter:
1247 . ksp - KSP context
1248
1249 Output Parameter:
1250 . usrP - user context
1251
1252 Fortran Notes:
1253 To use this from Fortran you must write a Fortran interface definition for this
1254 function that tells Fortran the Fortran derived data type that you are passing in as the ctx argument.
1255
1256 Level: intermediate
1257
1258 .seealso: KSPSetApplicationContext()
1259 @*/
KSPGetApplicationContext(KSP ksp,void * usrP)1260 PetscErrorCode KSPGetApplicationContext(KSP ksp,void *usrP)
1261 {
1262 PetscFunctionBegin;
1263 PetscValidHeaderSpecific(ksp,KSP_CLASSID,1);
1264 *(void**)usrP = ksp->user;
1265 PetscFunctionReturn(0);
1266 }
1267
1268 #include <petsc/private/pcimpl.h>
1269
1270 /*@
1271 KSPCheckSolve - Checks if the PCSetUp() or KSPSolve() failed and set the error flag for the outer PC. A KSP_DIVERGED_ITS is
1272 not considered a failure in this context
1273
1274 Collective on ksp
1275
1276 Input Parameter:
1277 + ksp - the linear solver (KSP) context.
1278 . pc - the preconditioner context
1279 - vec - a vector that will be initialized with Inf to indicate lack of convergence
1280
1281 Notes: this may be called by a subset of the processes in the PC
1282
1283 Level: developer
1284
1285 Developer Note: this is used to manage returning from preconditioners whose inner KSP solvers have failed in some way
1286
1287 .seealso: KSPCreate(), KSPSetType(), KSP, KSPCheckNorm(), KSPCheckDot()
1288 @*/
KSPCheckSolve(KSP ksp,PC pc,Vec vec)1289 PetscErrorCode KSPCheckSolve(KSP ksp,PC pc,Vec vec)
1290 {
1291 PetscErrorCode ierr;
1292 PCFailedReason pcreason;
1293 PC subpc;
1294
1295 PetscFunctionBegin;
1296 ierr = KSPGetPC(ksp,&subpc);CHKERRQ(ierr);
1297 ierr = PCGetFailedReason(subpc,&pcreason);CHKERRQ(ierr);
1298 if (pcreason || (ksp->reason < 0 && ksp->reason != KSP_DIVERGED_ITS)) {
1299 if (pc->erroriffailure) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_NOT_CONVERGED,"Detected not converged in KSP inner solve: KSP reason %s PC reason %s",KSPConvergedReasons[ksp->reason],PCFailedReasons[pcreason]);
1300 else {
1301 ierr = PetscInfo2(ksp,"Detected not converged in KSP inner solve: KSP reason %s PC reason %s\n",KSPConvergedReasons[ksp->reason],PCFailedReasons[pcreason]);CHKERRQ(ierr);
1302 pc->failedreason = PC_SUBPC_ERROR;
1303 if (vec) {
1304 ierr = VecSetInf(vec);CHKERRQ(ierr);
1305 }
1306 }
1307 }
1308 PetscFunctionReturn(0);
1309 }
1310