1 /*
2 Interface KSP routines that the user calls.
3 */
4
5 #include <petsc/private/kspimpl.h> /*I "petscksp.h" I*/
6 #include <petsc/private/matimpl.h> /*I "petscmat.h" I*/
7 #include <petscdm.h>
8
ObjectView(PetscObject obj,PetscViewer viewer,PetscViewerFormat format)9 PETSC_STATIC_INLINE PetscErrorCode ObjectView(PetscObject obj, PetscViewer viewer, PetscViewerFormat format)
10 {
11 PetscErrorCode ierr;
12
13 ierr = PetscViewerPushFormat(viewer, format);CHKERRQ(ierr);
14 ierr = PetscObjectView(obj, viewer);CHKERRQ(ierr);
15 ierr = PetscViewerPopFormat(viewer);CHKERRQ(ierr);
16 return(0);
17 }
18
19 /*@
20 KSPComputeExtremeSingularValues - Computes the extreme singular values
21 for the preconditioned operator. Called after or during KSPSolve().
22
23 Not Collective
24
25 Input Parameter:
26 . ksp - iterative context obtained from KSPCreate()
27
28 Output Parameters:
29 . emin, emax - extreme singular values
30
31 Options Database Keys:
32 . -ksp_view_singularvalues - compute extreme singular values and print when KSPSolve completes.
33
34 Notes:
35 One must call KSPSetComputeSingularValues() before calling KSPSetUp()
36 (or use the option -ksp_view_eigenvalues) in order for this routine to work correctly.
37
38 Many users may just want to use the monitoring routine
39 KSPMonitorSingularValue() (which can be set with option -ksp_monitor_singular_value)
40 to print the extreme singular values at each iteration of the linear solve.
41
42 Estimates of the smallest singular value may be very inaccurate, especially if the Krylov method has not converged.
43 The largest singular value is usually accurate to within a few percent if the method has converged, but is still not
44 intended for eigenanalysis.
45
46 Disable restarts if using KSPGMRES, otherwise this estimate will only be using those iterations after the last
47 restart. See KSPGMRESSetRestart() for more details.
48
49 Level: advanced
50
51 .seealso: KSPSetComputeSingularValues(), KSPMonitorSingularValue(), KSPComputeEigenvalues(), KSP
52 @*/
KSPComputeExtremeSingularValues(KSP ksp,PetscReal * emax,PetscReal * emin)53 PetscErrorCode KSPComputeExtremeSingularValues(KSP ksp,PetscReal *emax,PetscReal *emin)
54 {
55 PetscErrorCode ierr;
56
57 PetscFunctionBegin;
58 PetscValidHeaderSpecific(ksp,KSP_CLASSID,1);
59 PetscValidScalarPointer(emax,2);
60 PetscValidScalarPointer(emin,3);
61 if (!ksp->calc_sings) SETERRQ(PetscObjectComm((PetscObject)ksp),4,"Singular values not requested before KSPSetUp()");
62
63 if (ksp->ops->computeextremesingularvalues) {
64 ierr = (*ksp->ops->computeextremesingularvalues)(ksp,emax,emin);CHKERRQ(ierr);
65 } else {
66 *emin = -1.0;
67 *emax = -1.0;
68 }
69 PetscFunctionReturn(0);
70 }
71
72 /*@
73 KSPComputeEigenvalues - Computes the extreme eigenvalues for the
74 preconditioned operator. Called after or during KSPSolve().
75
76 Not Collective
77
78 Input Parameters:
79 + ksp - iterative context obtained from KSPCreate()
80 - n - size of arrays r and c. The number of eigenvalues computed (neig) will, in
81 general, be less than this.
82
83 Output Parameters:
84 + r - real part of computed eigenvalues, provided by user with a dimension of at least n
85 . c - complex part of computed eigenvalues, provided by user with a dimension of at least n
86 - neig - actual number of eigenvalues computed (will be less than or equal to n)
87
88 Options Database Keys:
89 . -ksp_view_eigenvalues - Prints eigenvalues to stdout
90
91 Notes:
92 The number of eigenvalues estimated depends on the size of the Krylov space
93 generated during the KSPSolve() ; for example, with
94 CG it corresponds to the number of CG iterations, for GMRES it is the number
95 of GMRES iterations SINCE the last restart. Any extra space in r[] and c[]
96 will be ignored.
97
98 KSPComputeEigenvalues() does not usually provide accurate estimates; it is
99 intended only for assistance in understanding the convergence of iterative
100 methods, not for eigenanalysis. For accurate computation of eigenvalues we recommend using
101 the excellent package SLEPc.
102
103 One must call KSPSetComputeEigenvalues() before calling KSPSetUp()
104 in order for this routine to work correctly.
105
106 Many users may just want to use the monitoring routine
107 KSPMonitorSingularValue() (which can be set with option -ksp_monitor_singular_value)
108 to print the singular values at each iteration of the linear solve.
109
110 Level: advanced
111
112 .seealso: KSPSetComputeSingularValues(), KSPMonitorSingularValue(), KSPComputeExtremeSingularValues(), KSP
113 @*/
KSPComputeEigenvalues(KSP ksp,PetscInt n,PetscReal r[],PetscReal c[],PetscInt * neig)114 PetscErrorCode KSPComputeEigenvalues(KSP ksp,PetscInt n,PetscReal r[],PetscReal c[],PetscInt *neig)
115 {
116 PetscErrorCode ierr;
117
118 PetscFunctionBegin;
119 PetscValidHeaderSpecific(ksp,KSP_CLASSID,1);
120 if (n) PetscValidScalarPointer(r,3);
121 if (n) PetscValidScalarPointer(c,4);
122 if (n<0) SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_ARG_OUTOFRANGE,"Requested < 0 Eigenvalues");
123 PetscValidIntPointer(neig,5);
124 if (!ksp->calc_sings) SETERRQ(PetscObjectComm((PetscObject)ksp),4,"Eigenvalues not requested before KSPSetUp()");
125
126 if (n && ksp->ops->computeeigenvalues) {
127 ierr = (*ksp->ops->computeeigenvalues)(ksp,n,r,c,neig);CHKERRQ(ierr);
128 } else {
129 *neig = 0;
130 }
131 PetscFunctionReturn(0);
132 }
133
134 /*@
135 KSPComputeRitz - Computes the Ritz or harmonic Ritz pairs associated to the
136 smallest or largest in modulus, for the preconditioned operator.
137 Called after KSPSolve().
138
139 Not Collective
140
141 Input Parameters:
142 + ksp - iterative context obtained from KSPCreate()
143 . ritz - PETSC_TRUE or PETSC_FALSE for ritz pairs or harmonic Ritz pairs, respectively
144 . small - PETSC_TRUE or PETSC_FALSE for smallest or largest (harmonic) Ritz values, respectively
145 - nrit - number of (harmonic) Ritz pairs to compute
146
147 Output Parameters:
148 + nrit - actual number of computed (harmonic) Ritz pairs
149 . S - multidimensional vector with Ritz vectors
150 . tetar - real part of the Ritz values
151 - tetai - imaginary part of the Ritz values
152
153 Notes:
154 -For GMRES, the (harmonic) Ritz pairs are computed from the Hessenberg matrix obtained during
155 the last complete cycle, or obtained at the end of the solution if the method is stopped before
156 a restart. Then, the number of actual (harmonic) Ritz pairs computed is less or equal to the restart
157 parameter for GMRES if a complete cycle has been performed or less or equal to the number of GMRES
158 iterations.
159 -Moreover, for real matrices, the (harmonic) Ritz pairs are possibly complex-valued. In such a case,
160 the routine selects the complex (harmonic) Ritz value and its conjugate, and two successive columns of S
161 are equal to the real and the imaginary parts of the associated vectors.
162 -the (harmonic) Ritz pairs are given in order of increasing (harmonic) Ritz values in modulus
163 -this is currently not implemented when PETSc is built with complex numbers
164
165 One must call KSPSetComputeRitz() before calling KSPSetUp()
166 in order for this routine to work correctly.
167
168 Level: advanced
169
170 .seealso: KSPSetComputeRitz(), KSP
171 @*/
KSPComputeRitz(KSP ksp,PetscBool ritz,PetscBool small,PetscInt * nrit,Vec S[],PetscReal tetar[],PetscReal tetai[])172 PetscErrorCode KSPComputeRitz(KSP ksp,PetscBool ritz,PetscBool small,PetscInt *nrit,Vec S[],PetscReal tetar[],PetscReal tetai[])
173 {
174 PetscErrorCode ierr;
175
176 PetscFunctionBegin;
177 PetscValidHeaderSpecific(ksp,KSP_CLASSID,1);
178 if (!ksp->calc_ritz) SETERRQ(PetscObjectComm((PetscObject)ksp),4,"Ritz pairs not requested before KSPSetUp()");
179 if (ksp->ops->computeritz) {ierr = (*ksp->ops->computeritz)(ksp,ritz,small,nrit,S,tetar,tetai);CHKERRQ(ierr);}
180 PetscFunctionReturn(0);
181 }
182 /*@
183 KSPSetUpOnBlocks - Sets up the preconditioner for each block in
184 the block Jacobi, block Gauss-Seidel, and overlapping Schwarz
185 methods.
186
187 Collective on ksp
188
189 Input Parameter:
190 . ksp - the KSP context
191
192 Notes:
193 KSPSetUpOnBlocks() is a routine that the user can optinally call for
194 more precise profiling (via -log_view) of the setup phase for these
195 block preconditioners. If the user does not call KSPSetUpOnBlocks(),
196 it will automatically be called from within KSPSolve().
197
198 Calling KSPSetUpOnBlocks() is the same as calling PCSetUpOnBlocks()
199 on the PC context within the KSP context.
200
201 Level: advanced
202
203 .seealso: PCSetUpOnBlocks(), KSPSetUp(), PCSetUp(), KSP
204 @*/
KSPSetUpOnBlocks(KSP ksp)205 PetscErrorCode KSPSetUpOnBlocks(KSP ksp)
206 {
207 PC pc;
208 PetscErrorCode ierr;
209 PCFailedReason pcreason;
210
211 PetscFunctionBegin;
212 PetscValidHeaderSpecific(ksp,KSP_CLASSID,1);
213 ierr = KSPGetPC(ksp,&pc);CHKERRQ(ierr);
214 ierr = PCSetUpOnBlocks(pc);CHKERRQ(ierr);
215 ierr = PCGetFailedReasonRank(pc,&pcreason);CHKERRQ(ierr);
216 /* TODO: this code was wrong and is still wrong, there is no way to propogate the failure to all processes; their is no code to handle a ksp->reason on only some ranks */
217 if (pcreason) {
218 ksp->reason = KSP_DIVERGED_PC_FAILED;
219 }
220 PetscFunctionReturn(0);
221 }
222
223 /*@
224 KSPSetReusePreconditioner - reuse the current preconditioner, do not construct a new one even if the operator changes
225
226 Collective on ksp
227
228 Input Parameters:
229 + ksp - iterative context obtained from KSPCreate()
230 - flag - PETSC_TRUE to reuse the current preconditioner
231
232 Level: intermediate
233
234 .seealso: KSPCreate(), KSPSolve(), KSPDestroy(), PCSetReusePreconditioner(), KSP
235 @*/
KSPSetReusePreconditioner(KSP ksp,PetscBool flag)236 PetscErrorCode KSPSetReusePreconditioner(KSP ksp,PetscBool flag)
237 {
238 PC pc;
239 PetscErrorCode ierr;
240
241 PetscFunctionBegin;
242 PetscValidHeaderSpecific(ksp,KSP_CLASSID,1);
243 ierr = KSPGetPC(ksp,&pc);CHKERRQ(ierr);
244 ierr = PCSetReusePreconditioner(pc,flag);CHKERRQ(ierr);
245 PetscFunctionReturn(0);
246 }
247
248 /*@
249 KSPGetReusePreconditioner - Determines if the KSP reuses the current preconditioner even if the operator in the preconditioner has changed.
250
251 Collective on ksp
252
253 Input Parameters:
254 . ksp - iterative context obtained from KSPCreate()
255
256 Output Parameters:
257 . flag - the boolean flag
258
259 Level: intermediate
260
261 .seealso: KSPCreate(), KSPSolve(), KSPDestroy(), KSPSetReusePreconditioner(), KSP
262 @*/
KSPGetReusePreconditioner(KSP ksp,PetscBool * flag)263 PetscErrorCode KSPGetReusePreconditioner(KSP ksp,PetscBool *flag)
264 {
265 PetscErrorCode ierr;
266
267 PetscFunctionBegin;
268 PetscValidHeaderSpecific(ksp,KSP_CLASSID,1);
269 PetscValidPointer(flag,2);
270 *flag = PETSC_FALSE;
271 if (ksp->pc) {
272 ierr = PCGetReusePreconditioner(ksp->pc,flag);CHKERRQ(ierr);
273 }
274 PetscFunctionReturn(0);
275 }
276
277 /*@
278 KSPSetSkipPCSetFromOptions - prevents KSPSetFromOptions() from call PCSetFromOptions(). This is used if the same PC is shared by more than one KSP so its options are not resetable for each KSP
279
280 Collective on ksp
281
282 Input Parameters:
283 + ksp - iterative context obtained from KSPCreate()
284 - flag - PETSC_TRUE to skip calling the PCSetFromOptions()
285
286 Level: intermediate
287
288 .seealso: KSPCreate(), KSPSolve(), KSPDestroy(), PCSetReusePreconditioner(), KSP
289 @*/
KSPSetSkipPCSetFromOptions(KSP ksp,PetscBool flag)290 PetscErrorCode KSPSetSkipPCSetFromOptions(KSP ksp,PetscBool flag)
291 {
292 PetscFunctionBegin;
293 PetscValidHeaderSpecific(ksp,KSP_CLASSID,1);
294 ksp->skippcsetfromoptions = flag;
295 PetscFunctionReturn(0);
296 }
297
298 /*@
299 KSPSetUp - Sets up the internal data structures for the
300 later use of an iterative solver.
301
302 Collective on ksp
303
304 Input Parameter:
305 . ksp - iterative context obtained from KSPCreate()
306
307 Level: developer
308
309 .seealso: KSPCreate(), KSPSolve(), KSPDestroy(), KSP
310 @*/
KSPSetUp(KSP ksp)311 PetscErrorCode KSPSetUp(KSP ksp)
312 {
313 PetscErrorCode ierr;
314 Mat A,B;
315 Mat mat,pmat;
316 MatNullSpace nullsp;
317 PCFailedReason pcreason;
318
319 PetscFunctionBegin;
320 PetscValidHeaderSpecific(ksp,KSP_CLASSID,1);
321
322 /* reset the convergence flag from the previous solves */
323 ksp->reason = KSP_CONVERGED_ITERATING;
324
325 if (!((PetscObject)ksp)->type_name) {
326 ierr = KSPSetType(ksp,KSPGMRES);CHKERRQ(ierr);
327 }
328 ierr = KSPSetUpNorms_Private(ksp,PETSC_TRUE,&ksp->normtype,&ksp->pc_side);CHKERRQ(ierr);
329
330 if (ksp->dmActive && !ksp->setupstage) {
331 /* first time in so build matrix and vector data structures using DM */
332 if (!ksp->vec_rhs) {ierr = DMCreateGlobalVector(ksp->dm,&ksp->vec_rhs);CHKERRQ(ierr);}
333 if (!ksp->vec_sol) {ierr = DMCreateGlobalVector(ksp->dm,&ksp->vec_sol);CHKERRQ(ierr);}
334 ierr = DMCreateMatrix(ksp->dm,&A);CHKERRQ(ierr);
335 ierr = KSPSetOperators(ksp,A,A);CHKERRQ(ierr);
336 ierr = PetscObjectDereference((PetscObject)A);CHKERRQ(ierr);
337 }
338
339 if (ksp->dmActive) {
340 DMKSP kdm;
341 ierr = DMGetDMKSP(ksp->dm,&kdm);CHKERRQ(ierr);
342
343 if (kdm->ops->computeinitialguess && ksp->setupstage != KSP_SETUP_NEWRHS) {
344 /* only computes initial guess the first time through */
345 ierr = (*kdm->ops->computeinitialguess)(ksp,ksp->vec_sol,kdm->initialguessctx);CHKERRQ(ierr);
346 ierr = KSPSetInitialGuessNonzero(ksp,PETSC_TRUE);CHKERRQ(ierr);
347 }
348 if (kdm->ops->computerhs) {
349 ierr = (*kdm->ops->computerhs)(ksp,ksp->vec_rhs,kdm->rhsctx);CHKERRQ(ierr);
350 }
351
352 if (ksp->setupstage != KSP_SETUP_NEWRHS) {
353 if (kdm->ops->computeoperators) {
354 ierr = KSPGetOperators(ksp,&A,&B);CHKERRQ(ierr);
355 ierr = (*kdm->ops->computeoperators)(ksp,A,B,kdm->operatorsctx);CHKERRQ(ierr);
356 } else SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_ARG_WRONGSTATE,"You called KSPSetDM() but did not use DMKSPSetComputeOperators() or KSPSetDMActive(ksp,PETSC_FALSE);");
357 }
358 }
359
360 if (ksp->setupstage == KSP_SETUP_NEWRHS) PetscFunctionReturn(0);
361 ierr = PetscLogEventBegin(KSP_SetUp,ksp,ksp->vec_rhs,ksp->vec_sol,0);CHKERRQ(ierr);
362
363 switch (ksp->setupstage) {
364 case KSP_SETUP_NEW:
365 ierr = (*ksp->ops->setup)(ksp);CHKERRQ(ierr);
366 break;
367 case KSP_SETUP_NEWMATRIX: { /* This should be replaced with a more general mechanism */
368 if (ksp->setupnewmatrix) {
369 ierr = (*ksp->ops->setup)(ksp);CHKERRQ(ierr);
370 }
371 } break;
372 default: break;
373 }
374
375 if (!ksp->pc) {ierr = KSPGetPC(ksp,&ksp->pc);CHKERRQ(ierr);}
376 ierr = PCGetOperators(ksp->pc,&mat,&pmat);CHKERRQ(ierr);
377 /* scale the matrix if requested */
378 if (ksp->dscale) {
379 PetscScalar *xx;
380 PetscInt i,n;
381 PetscBool zeroflag = PETSC_FALSE;
382 if (!ksp->pc) {ierr = KSPGetPC(ksp,&ksp->pc);CHKERRQ(ierr);}
383 if (!ksp->diagonal) { /* allocate vector to hold diagonal */
384 ierr = MatCreateVecs(pmat,&ksp->diagonal,NULL);CHKERRQ(ierr);
385 }
386 ierr = MatGetDiagonal(pmat,ksp->diagonal);CHKERRQ(ierr);
387 ierr = VecGetLocalSize(ksp->diagonal,&n);CHKERRQ(ierr);
388 ierr = VecGetArray(ksp->diagonal,&xx);CHKERRQ(ierr);
389 for (i=0; i<n; i++) {
390 if (xx[i] != 0.0) xx[i] = 1.0/PetscSqrtReal(PetscAbsScalar(xx[i]));
391 else {
392 xx[i] = 1.0;
393 zeroflag = PETSC_TRUE;
394 }
395 }
396 ierr = VecRestoreArray(ksp->diagonal,&xx);CHKERRQ(ierr);
397 if (zeroflag) {
398 ierr = PetscInfo(ksp,"Zero detected in diagonal of matrix, using 1 at those locations\n");CHKERRQ(ierr);
399 }
400 ierr = MatDiagonalScale(pmat,ksp->diagonal,ksp->diagonal);CHKERRQ(ierr);
401 if (mat != pmat) {ierr = MatDiagonalScale(mat,ksp->diagonal,ksp->diagonal);CHKERRQ(ierr);}
402 ksp->dscalefix2 = PETSC_FALSE;
403 }
404 ierr = PetscLogEventEnd(KSP_SetUp,ksp,ksp->vec_rhs,ksp->vec_sol,0);CHKERRQ(ierr);
405 ierr = PCSetErrorIfFailure(ksp->pc,ksp->errorifnotconverged);CHKERRQ(ierr);
406 ierr = PCSetUp(ksp->pc);CHKERRQ(ierr);
407 ierr = PCGetFailedReasonRank(ksp->pc,&pcreason);CHKERRQ(ierr);
408 /* TODO: this code was wrong and is still wrong, there is no way to propogate the failure to all processes; their is no code to handle a ksp->reason on only some ranks */
409 if (pcreason) {
410 ksp->reason = KSP_DIVERGED_PC_FAILED;
411 }
412
413 ierr = MatGetNullSpace(mat,&nullsp);CHKERRQ(ierr);
414 if (nullsp) {
415 PetscBool test = PETSC_FALSE;
416 ierr = PetscOptionsGetBool(((PetscObject)ksp)->options,((PetscObject)ksp)->prefix,"-ksp_test_null_space",&test,NULL);CHKERRQ(ierr);
417 if (test) {
418 ierr = MatNullSpaceTest(nullsp,mat,NULL);CHKERRQ(ierr);
419 }
420 }
421 ksp->setupstage = KSP_SETUP_NEWRHS;
422 PetscFunctionReturn(0);
423 }
424
425 /*@C
426 KSPConvergedReasonView - Displays the reason a KSP solve converged or diverged to a viewer
427
428 Collective on ksp
429
430 Parameter:
431 + ksp - iterative context obtained from KSPCreate()
432 - viewer - the viewer to display the reason
433
434 Options Database Keys:
435 + -ksp_converged_reason - print reason for converged or diverged, also prints number of iterations
436 - -ksp_converged_reason ::failed - only print reason and number of iterations when diverged
437
438 Notes:
439 To change the format of the output call PetscViewerPushFormat(viewer,format) before this call. Use PETSC_VIEWER_DEFAULT for the default,
440 use PETSC_VIEWER_FAILED to only display a reason if it fails.
441
442 Level: beginner
443
444 .seealso: KSPCreate(), KSPSetUp(), KSPDestroy(), KSPSetTolerances(), KSPConvergedDefault(),
445 KSPSolveTranspose(), KSPGetIterationNumber(), KSP, KSPGetConvergedReason(), PetscViewerPushFormat(), PetscViewerPopFormat()
446 @*/
KSPConvergedReasonView(KSP ksp,PetscViewer viewer)447 PetscErrorCode KSPConvergedReasonView(KSP ksp, PetscViewer viewer)
448 {
449 PetscErrorCode ierr;
450 PetscBool isAscii;
451 PetscViewerFormat format;
452
453 PetscFunctionBegin;
454 if (!viewer) viewer = PETSC_VIEWER_STDOUT_(PetscObjectComm((PetscObject)ksp));
455 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isAscii);CHKERRQ(ierr);
456 if (isAscii) {
457 ierr = PetscViewerGetFormat(viewer, &format);CHKERRQ(ierr);
458 ierr = PetscViewerASCIIAddTab(viewer,((PetscObject)ksp)->tablevel);CHKERRQ(ierr);
459 if (ksp->reason > 0 && format != PETSC_VIEWER_FAILED) {
460 if (((PetscObject) ksp)->prefix) {
461 ierr = PetscViewerASCIIPrintf(viewer,"Linear %s solve converged due to %s iterations %D\n",((PetscObject) ksp)->prefix,KSPConvergedReasons[ksp->reason],ksp->its);CHKERRQ(ierr);
462 } else {
463 ierr = PetscViewerASCIIPrintf(viewer,"Linear solve converged due to %s iterations %D\n",KSPConvergedReasons[ksp->reason],ksp->its);CHKERRQ(ierr);
464 }
465 } else if (ksp->reason <= 0) {
466 if (((PetscObject) ksp)->prefix) {
467 ierr = PetscViewerASCIIPrintf(viewer,"Linear %s solve did not converge due to %s iterations %D\n",((PetscObject) ksp)->prefix,KSPConvergedReasons[ksp->reason],ksp->its);CHKERRQ(ierr);
468 } else {
469 ierr = PetscViewerASCIIPrintf(viewer,"Linear solve did not converge due to %s iterations %D\n",KSPConvergedReasons[ksp->reason],ksp->its);CHKERRQ(ierr);
470 }
471 if (ksp->reason == KSP_DIVERGED_PC_FAILED) {
472 PCFailedReason reason;
473 ierr = PCGetFailedReason(ksp->pc,&reason);CHKERRQ(ierr);
474 ierr = PetscViewerASCIIPrintf(viewer," PC failed due to %s \n",PCFailedReasons[reason]);CHKERRQ(ierr);
475 }
476 }
477 ierr = PetscViewerASCIISubtractTab(viewer,((PetscObject)ksp)->tablevel);CHKERRQ(ierr);
478 }
479 PetscFunctionReturn(0);
480 }
481
482 /*@
483 KSPConvergedReasonViewFromOptions - Processes command line options to determine if/how a KSPReason is to be viewed.
484
485 Collective on ksp
486
487 Input Parameters:
488 . ksp - the KSP object
489
490 Level: intermediate
491
492 @*/
KSPConvergedReasonViewFromOptions(KSP ksp)493 PetscErrorCode KSPConvergedReasonViewFromOptions(KSP ksp)
494 {
495 PetscViewer viewer;
496 PetscBool flg;
497 PetscViewerFormat format;
498 PetscErrorCode ierr;
499
500 PetscFunctionBegin;
501 ierr = PetscOptionsGetViewer(PetscObjectComm((PetscObject)ksp),((PetscObject)ksp)->options,((PetscObject)ksp)->prefix,"-ksp_converged_reason",&viewer,&format,&flg);CHKERRQ(ierr);
502 if (flg) {
503 ierr = PetscViewerPushFormat(viewer,format);CHKERRQ(ierr);
504 ierr = KSPConvergedReasonView(ksp, viewer);CHKERRQ(ierr);
505 ierr = PetscViewerPopFormat(viewer);CHKERRQ(ierr);
506 ierr = PetscViewerDestroy(&viewer);CHKERRQ(ierr);
507 }
508 PetscFunctionReturn(0);
509 }
510
511 #include <petscdraw.h>
512
KSPViewEigenvalues_Internal(KSP ksp,PetscBool isExplicit,PetscViewer viewer,PetscViewerFormat format)513 static PetscErrorCode KSPViewEigenvalues_Internal(KSP ksp, PetscBool isExplicit, PetscViewer viewer, PetscViewerFormat format)
514 {
515 PetscReal *r, *c;
516 PetscInt n, i, neig;
517 PetscBool isascii, isdraw;
518 PetscMPIInt rank;
519 PetscErrorCode ierr;
520
521 PetscFunctionBegin;
522 ierr = MPI_Comm_rank(PetscObjectComm((PetscObject) ksp), &rank);CHKERRQ(ierr);
523 ierr = PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &isascii);CHKERRQ(ierr);
524 ierr = PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERDRAW, &isdraw);CHKERRQ(ierr);
525 if (isExplicit) {
526 ierr = VecGetSize(ksp->vec_sol,&n);CHKERRQ(ierr);
527 ierr = PetscMalloc2(n, &r, n, &c);CHKERRQ(ierr);
528 ierr = KSPComputeEigenvaluesExplicitly(ksp, n, r, c);CHKERRQ(ierr);
529 neig = n;
530 } else {
531 PetscInt nits;
532
533 ierr = KSPGetIterationNumber(ksp, &nits);CHKERRQ(ierr);
534 n = nits+2;
535 if (!nits) {ierr = PetscViewerASCIIPrintf(viewer, "Zero iterations in solver, cannot approximate any eigenvalues\n");CHKERRQ(ierr);PetscFunctionReturn(0);}
536 ierr = PetscMalloc2(n, &r, n, &c);CHKERRQ(ierr);
537 ierr = KSPComputeEigenvalues(ksp, n, r, c, &neig);CHKERRQ(ierr);
538 }
539 if (isascii) {
540 ierr = PetscViewerASCIIPrintf(viewer, "%s computed eigenvalues\n", isExplicit ? "Explicitly" : "Iteratively");CHKERRQ(ierr);
541 for (i = 0; i < neig; ++i) {
542 if (c[i] >= 0.0) {ierr = PetscViewerASCIIPrintf(viewer, "%g + %gi\n", (double) r[i], (double) c[i]);CHKERRQ(ierr);}
543 else {ierr = PetscViewerASCIIPrintf(viewer, "%g - %gi\n", (double) r[i], -(double) c[i]);CHKERRQ(ierr);}
544 }
545 } else if (isdraw && !rank) {
546 PetscDraw draw;
547 PetscDrawSP drawsp;
548
549 if (format == PETSC_VIEWER_DRAW_CONTOUR) {
550 ierr = KSPPlotEigenContours_Private(ksp,neig,r,c);CHKERRQ(ierr);
551 } else {
552 if (!ksp->eigviewer) {ierr = PetscViewerDrawOpen(PETSC_COMM_SELF,NULL,isExplicit ? "Explicitly Computed Eigenvalues" : "Iteratively Computed Eigenvalues",PETSC_DECIDE,PETSC_DECIDE,400,400,&ksp->eigviewer);CHKERRQ(ierr);}
553 ierr = PetscViewerDrawGetDraw(ksp->eigviewer,0,&draw);CHKERRQ(ierr);
554 ierr = PetscDrawSPCreate(draw,1,&drawsp);CHKERRQ(ierr);
555 ierr = PetscDrawSPReset(drawsp);CHKERRQ(ierr);
556 for (i = 0; i < neig; ++i) {ierr = PetscDrawSPAddPoint(drawsp,r+i,c+i);CHKERRQ(ierr);}
557 ierr = PetscDrawSPDraw(drawsp,PETSC_TRUE);CHKERRQ(ierr);
558 ierr = PetscDrawSPSave(drawsp);CHKERRQ(ierr);
559 ierr = PetscDrawSPDestroy(&drawsp);CHKERRQ(ierr);
560 }
561 }
562 ierr = PetscFree2(r, c);CHKERRQ(ierr);
563 PetscFunctionReturn(0);
564 }
565
KSPViewSingularvalues_Internal(KSP ksp,PetscViewer viewer,PetscViewerFormat format)566 static PetscErrorCode KSPViewSingularvalues_Internal(KSP ksp, PetscViewer viewer, PetscViewerFormat format)
567 {
568 PetscReal smax, smin;
569 PetscInt nits;
570 PetscBool isascii;
571 PetscErrorCode ierr;
572
573 PetscFunctionBegin;
574 ierr = PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &isascii);CHKERRQ(ierr);
575 ierr = KSPGetIterationNumber(ksp, &nits);CHKERRQ(ierr);
576 if (!nits) {ierr = PetscViewerASCIIPrintf(viewer, "Zero iterations in solver, cannot approximate any singular values\n");CHKERRQ(ierr);PetscFunctionReturn(0);}
577 ierr = KSPComputeExtremeSingularValues(ksp, &smax, &smin);CHKERRQ(ierr);
578 if (isascii) {ierr = PetscViewerASCIIPrintf(viewer, "Iteratively computed extreme singular values: max %g min %g max/min %g\n",(double)smax,(double)smin,(double)(smax/smin));CHKERRQ(ierr);}
579 PetscFunctionReturn(0);
580 }
581
KSPViewFinalResidual_Internal(KSP ksp,PetscViewer viewer,PetscViewerFormat format)582 static PetscErrorCode KSPViewFinalResidual_Internal(KSP ksp, PetscViewer viewer, PetscViewerFormat format)
583 {
584 PetscBool isascii;
585 PetscErrorCode ierr;
586
587 PetscFunctionBegin;
588 ierr = PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &isascii);CHKERRQ(ierr);
589 if (ksp->dscale && !ksp->dscalefix) SETERRQ(PetscObjectComm((PetscObject) ksp), PETSC_ERR_ARG_WRONGSTATE, "Cannot compute final scale with -ksp_diagonal_scale except also with -ksp_diagonal_scale_fix");
590 if (isascii) {
591 Mat A;
592 Vec t;
593 PetscReal norm;
594
595 ierr = PCGetOperators(ksp->pc, &A, NULL);CHKERRQ(ierr);
596 ierr = VecDuplicate(ksp->vec_rhs, &t);CHKERRQ(ierr);
597 ierr = KSP_MatMult(ksp, A, ksp->vec_sol, t);CHKERRQ(ierr);
598 ierr = VecAYPX(t, -1.0, ksp->vec_rhs);CHKERRQ(ierr);
599 ierr = VecNorm(t, NORM_2, &norm);CHKERRQ(ierr);
600 ierr = VecDestroy(&t);CHKERRQ(ierr);
601 ierr = PetscViewerASCIIPrintf(viewer, "KSP final norm of residual %g\n", (double) norm);CHKERRQ(ierr);
602 }
603 PetscFunctionReturn(0);
604 }
605
KSPSolve_Private(KSP ksp,Vec b,Vec x)606 static PetscErrorCode KSPSolve_Private(KSP ksp,Vec b,Vec x)
607 {
608 PetscErrorCode ierr;
609 PetscBool flg = PETSC_FALSE,inXisinB=PETSC_FALSE,guess_zero;
610 Mat mat,pmat;
611 MPI_Comm comm;
612 MatNullSpace nullsp;
613 Vec btmp,vec_rhs=NULL;
614
615 PetscFunctionBegin;
616 comm = PetscObjectComm((PetscObject)ksp);
617 if (x && x == b) {
618 if (!ksp->guess_zero) SETERRQ(comm,PETSC_ERR_ARG_INCOMP,"Cannot use x == b with nonzero initial guess");
619 ierr = VecDuplicate(b,&x);CHKERRQ(ierr);
620 inXisinB = PETSC_TRUE;
621 }
622 if (b) {
623 ierr = PetscObjectReference((PetscObject)b);CHKERRQ(ierr);
624 ierr = VecDestroy(&ksp->vec_rhs);CHKERRQ(ierr);
625 ksp->vec_rhs = b;
626 }
627 if (x) {
628 ierr = PetscObjectReference((PetscObject)x);CHKERRQ(ierr);
629 ierr = VecDestroy(&ksp->vec_sol);CHKERRQ(ierr);
630 ksp->vec_sol = x;
631 }
632
633 if (ksp->viewPre) {ierr = ObjectView((PetscObject) ksp, ksp->viewerPre, ksp->formatPre);CHKERRQ(ierr);}
634
635 if (ksp->presolve) {ierr = (*ksp->presolve)(ksp,ksp->vec_rhs,ksp->vec_sol,ksp->prectx);CHKERRQ(ierr);}
636
637 /* reset the residual history list if requested */
638 if (ksp->res_hist_reset) ksp->res_hist_len = 0;
639
640 ierr = PetscLogEventBegin(KSP_Solve,ksp,ksp->vec_rhs,ksp->vec_sol,0);CHKERRQ(ierr);
641
642 if (ksp->guess) {
643 PetscObjectState ostate,state;
644
645 ierr = KSPGuessSetUp(ksp->guess);CHKERRQ(ierr);
646 ierr = PetscObjectStateGet((PetscObject)ksp->vec_sol,&ostate);CHKERRQ(ierr);
647 ierr = KSPGuessFormGuess(ksp->guess,ksp->vec_rhs,ksp->vec_sol);CHKERRQ(ierr);
648 ierr = PetscObjectStateGet((PetscObject)ksp->vec_sol,&state);CHKERRQ(ierr);
649 if (state != ostate) {
650 ksp->guess_zero = PETSC_FALSE;
651 } else {
652 ierr = PetscInfo(ksp,"Using zero initial guess since the KSPGuess object did not change the vector\n");CHKERRQ(ierr);
653 ksp->guess_zero = PETSC_TRUE;
654 }
655 }
656
657 /* KSPSetUp() scales the matrix if needed */
658 ierr = KSPSetUp(ksp);CHKERRQ(ierr);
659 ierr = KSPSetUpOnBlocks(ksp);CHKERRQ(ierr);
660
661 ierr = VecSetErrorIfLocked(ksp->vec_sol,3);CHKERRQ(ierr);
662
663 ierr = PCGetOperators(ksp->pc,&mat,&pmat);CHKERRQ(ierr);
664 /* diagonal scale RHS if called for */
665 if (ksp->dscale) {
666 ierr = VecPointwiseMult(ksp->vec_rhs,ksp->vec_rhs,ksp->diagonal);CHKERRQ(ierr);
667 /* second time in, but matrix was scaled back to original */
668 if (ksp->dscalefix && ksp->dscalefix2) {
669 Mat mat,pmat;
670
671 ierr = PCGetOperators(ksp->pc,&mat,&pmat);CHKERRQ(ierr);
672 ierr = MatDiagonalScale(pmat,ksp->diagonal,ksp->diagonal);CHKERRQ(ierr);
673 if (mat != pmat) {ierr = MatDiagonalScale(mat,ksp->diagonal,ksp->diagonal);CHKERRQ(ierr);}
674 }
675
676 /* scale initial guess */
677 if (!ksp->guess_zero) {
678 if (!ksp->truediagonal) {
679 ierr = VecDuplicate(ksp->diagonal,&ksp->truediagonal);CHKERRQ(ierr);
680 ierr = VecCopy(ksp->diagonal,ksp->truediagonal);CHKERRQ(ierr);
681 ierr = VecReciprocal(ksp->truediagonal);CHKERRQ(ierr);
682 }
683 ierr = VecPointwiseMult(ksp->vec_sol,ksp->vec_sol,ksp->truediagonal);CHKERRQ(ierr);
684 }
685 }
686 ierr = PCPreSolve(ksp->pc,ksp);CHKERRQ(ierr);
687
688 if (ksp->guess_zero) { ierr = VecSet(ksp->vec_sol,0.0);CHKERRQ(ierr);}
689 if (ksp->guess_knoll) { /* The Knoll trick is independent on the KSPGuess specified */
690 ierr = PCApply(ksp->pc,ksp->vec_rhs,ksp->vec_sol);CHKERRQ(ierr);
691 ierr = KSP_RemoveNullSpace(ksp,ksp->vec_sol);CHKERRQ(ierr);
692 ksp->guess_zero = PETSC_FALSE;
693 }
694
695 /* can we mark the initial guess as zero for this solve? */
696 guess_zero = ksp->guess_zero;
697 if (!ksp->guess_zero) {
698 PetscReal norm;
699
700 ierr = VecNormAvailable(ksp->vec_sol,NORM_2,&flg,&norm);CHKERRQ(ierr);
701 if (flg && !norm) ksp->guess_zero = PETSC_TRUE;
702 }
703 if (ksp->transpose_solve) {
704 ierr = MatGetNullSpace(pmat,&nullsp);CHKERRQ(ierr);
705 } else {
706 ierr = MatGetTransposeNullSpace(pmat,&nullsp);CHKERRQ(ierr);
707 }
708 if (nullsp) {
709 ierr = VecDuplicate(ksp->vec_rhs,&btmp);CHKERRQ(ierr);
710 ierr = VecCopy(ksp->vec_rhs,btmp);CHKERRQ(ierr);
711 ierr = MatNullSpaceRemove(nullsp,btmp);CHKERRQ(ierr);
712 vec_rhs = ksp->vec_rhs;
713 ksp->vec_rhs = btmp;
714 }
715 ierr = VecLockReadPush(ksp->vec_rhs);CHKERRQ(ierr);
716 if (ksp->reason == KSP_DIVERGED_PC_FAILED) {
717 ierr = VecSetInf(ksp->vec_sol);CHKERRQ(ierr);
718 }
719 ierr = (*ksp->ops->solve)(ksp);CHKERRQ(ierr);
720
721 ierr = VecLockReadPop(ksp->vec_rhs);CHKERRQ(ierr);
722 if (nullsp) {
723 ksp->vec_rhs = vec_rhs;
724 ierr = VecDestroy(&btmp);CHKERRQ(ierr);
725 }
726
727 ksp->guess_zero = guess_zero;
728
729 if (!ksp->reason) SETERRQ(comm,PETSC_ERR_PLIB,"Internal error, solver returned without setting converged reason");
730 ksp->totalits += ksp->its;
731
732 if (ksp->viewReason) {
733 ierr = PetscViewerPushFormat(ksp->viewerReason,ksp->formatReason);CHKERRQ(ierr);
734 ierr = KSPConvergedReasonView(ksp, ksp->viewerReason);CHKERRQ(ierr);
735 ierr = PetscViewerPopFormat(ksp->viewerReason);CHKERRQ(ierr);
736 }
737 ierr = PCPostSolve(ksp->pc,ksp);CHKERRQ(ierr);
738
739 /* diagonal scale solution if called for */
740 if (ksp->dscale) {
741 ierr = VecPointwiseMult(ksp->vec_sol,ksp->vec_sol,ksp->diagonal);CHKERRQ(ierr);
742 /* unscale right hand side and matrix */
743 if (ksp->dscalefix) {
744 Mat mat,pmat;
745
746 ierr = VecReciprocal(ksp->diagonal);CHKERRQ(ierr);
747 ierr = VecPointwiseMult(ksp->vec_rhs,ksp->vec_rhs,ksp->diagonal);CHKERRQ(ierr);
748 ierr = PCGetOperators(ksp->pc,&mat,&pmat);CHKERRQ(ierr);
749 ierr = MatDiagonalScale(pmat,ksp->diagonal,ksp->diagonal);CHKERRQ(ierr);
750 if (mat != pmat) {ierr = MatDiagonalScale(mat,ksp->diagonal,ksp->diagonal);CHKERRQ(ierr);}
751 ierr = VecReciprocal(ksp->diagonal);CHKERRQ(ierr);
752 ksp->dscalefix2 = PETSC_TRUE;
753 }
754 }
755 ierr = PetscLogEventEnd(KSP_Solve,ksp,ksp->vec_rhs,ksp->vec_sol,0);CHKERRQ(ierr);
756 if (ksp->guess) {
757 ierr = KSPGuessUpdate(ksp->guess,ksp->vec_rhs,ksp->vec_sol);CHKERRQ(ierr);
758 }
759 if (ksp->postsolve) {
760 ierr = (*ksp->postsolve)(ksp,ksp->vec_rhs,ksp->vec_sol,ksp->postctx);CHKERRQ(ierr);
761 }
762
763 ierr = PCGetOperators(ksp->pc,&mat,&pmat);CHKERRQ(ierr);
764 if (ksp->viewEV) {ierr = KSPViewEigenvalues_Internal(ksp, PETSC_FALSE, ksp->viewerEV, ksp->formatEV);CHKERRQ(ierr);}
765 if (ksp->viewEVExp) {ierr = KSPViewEigenvalues_Internal(ksp, PETSC_TRUE, ksp->viewerEVExp, ksp->formatEVExp);CHKERRQ(ierr);}
766 if (ksp->viewSV) {ierr = KSPViewSingularvalues_Internal(ksp, ksp->viewerSV, ksp->formatSV);CHKERRQ(ierr);}
767 if (ksp->viewFinalRes) {ierr = KSPViewFinalResidual_Internal(ksp, ksp->viewerFinalRes, ksp->formatFinalRes);CHKERRQ(ierr);}
768 if (ksp->viewMat) {ierr = ObjectView((PetscObject) mat, ksp->viewerMat, ksp->formatMat);CHKERRQ(ierr);}
769 if (ksp->viewPMat) {ierr = ObjectView((PetscObject) pmat, ksp->viewerPMat, ksp->formatPMat);CHKERRQ(ierr);}
770 if (ksp->viewRhs) {ierr = ObjectView((PetscObject) ksp->vec_rhs, ksp->viewerRhs, ksp->formatRhs);CHKERRQ(ierr);}
771 if (ksp->viewSol) {ierr = ObjectView((PetscObject) ksp->vec_sol, ksp->viewerSol, ksp->formatSol);CHKERRQ(ierr);}
772 if (ksp->view) {ierr = ObjectView((PetscObject) ksp, ksp->viewer, ksp->format);CHKERRQ(ierr);}
773 if (ksp->viewDScale) {ierr = ObjectView((PetscObject) ksp->diagonal, ksp->viewerDScale, ksp->formatDScale);CHKERRQ(ierr);}
774 if (ksp->viewMatExp) {
775 Mat A, B;
776
777 ierr = PCGetOperators(ksp->pc, &A, NULL);CHKERRQ(ierr);
778 if (ksp->transpose_solve) {
779 Mat AT;
780
781 ierr = MatCreateTranspose(A, &AT);CHKERRQ(ierr);
782 ierr = MatComputeOperator(AT, MATAIJ, &B);CHKERRQ(ierr);
783 ierr = MatDestroy(&AT);CHKERRQ(ierr);
784 } else {
785 ierr = MatComputeOperator(A, MATAIJ, &B);CHKERRQ(ierr);
786 }
787 ierr = ObjectView((PetscObject) B, ksp->viewerMatExp, ksp->formatMatExp);CHKERRQ(ierr);
788 ierr = MatDestroy(&B);CHKERRQ(ierr);
789 }
790 if (ksp->viewPOpExp) {
791 Mat B;
792
793 ierr = KSPComputeOperator(ksp, MATAIJ, &B);CHKERRQ(ierr);
794 ierr = ObjectView((PetscObject) B, ksp->viewerPOpExp, ksp->formatPOpExp);CHKERRQ(ierr);
795 ierr = MatDestroy(&B);CHKERRQ(ierr);
796 }
797
798 if (inXisinB) {
799 ierr = VecCopy(x,b);CHKERRQ(ierr);
800 ierr = VecDestroy(&x);CHKERRQ(ierr);
801 }
802 ierr = PetscObjectSAWsBlock((PetscObject)ksp);CHKERRQ(ierr);
803 if (ksp->errorifnotconverged && ksp->reason < 0 && ksp->reason != KSP_DIVERGED_ITS) {
804 if (ksp->reason == KSP_DIVERGED_PC_FAILED) {
805 PCFailedReason reason;
806 ierr = PCGetFailedReason(ksp->pc,&reason);CHKERRQ(ierr);
807 SETERRQ2(comm,PETSC_ERR_NOT_CONVERGED,"KSPSolve has not converged, reason %s PC failed due to %s",KSPConvergedReasons[ksp->reason],PCFailedReasons[reason]);
808 } else SETERRQ1(comm,PETSC_ERR_NOT_CONVERGED,"KSPSolve has not converged, reason %s",KSPConvergedReasons[ksp->reason]);
809 }
810 PetscFunctionReturn(0);
811 }
812
813 /*@
814 KSPSolve - Solves linear system.
815
816 Collective on ksp
817
818 Parameters:
819 + ksp - iterative context obtained from KSPCreate()
820 . b - the right hand side vector
821 - x - the solution (this may be the same vector as b, then b will be overwritten with answer)
822
823 Options Database Keys:
824 + -ksp_view_eigenvalues - compute preconditioned operators eigenvalues
825 . -ksp_view_eigenvalues_explicit - compute the eigenvalues by forming the dense operator and using LAPACK
826 . -ksp_view_mat binary - save matrix to the default binary viewer
827 . -ksp_view_pmat binary - save matrix used to build preconditioner to the default binary viewer
828 . -ksp_view_rhs binary - save right hand side vector to the default binary viewer
829 . -ksp_view_solution binary - save computed solution vector to the default binary viewer
830 (can be read later with src/ksp/tutorials/ex10.c for testing solvers)
831 . -ksp_view_mat_explicit - for matrix-free operators, computes the matrix entries and views them
832 . -ksp_view_preconditioned_operator_explicit - computes the product of the preconditioner and matrix as an explicit matrix and views it
833 . -ksp_converged_reason - print reason for converged or diverged, also prints number of iterations
834 . -ksp_view_final_residual - print 2-norm of true linear system residual at the end of the solution process
835 - -ksp_view - print the ksp data structure at the end of the system solution
836
837 Notes:
838
839 If one uses KSPSetDM() then x or b need not be passed. Use KSPGetSolution() to access the solution in this case.
840
841 The operator is specified with KSPSetOperators().
842
843 Call KSPGetConvergedReason() to determine if the solver converged or failed and
844 why. The number of iterations can be obtained from KSPGetIterationNumber().
845
846 If you provide a matrix that has a MatSetNullSpace() and MatSetTransposeNullSpace() this will use that information to solve singular systems
847 in the least squares sense with a norm minimizing solution.
848 $
849 $ A x = b where b = b_p + b_t where b_t is not in the range of A (and hence by the fundamental theorem of linear algebra is in the nullspace(A') see MatSetNullSpace()
850 $
851 $ KSP first removes b_t producing the linear system A x = b_p (which has multiple solutions) and solves this to find the ||x|| minimizing solution (and hence
852 $ it finds the solution x orthogonal to the nullspace(A). The algorithm is simply in each iteration of the Krylov method we remove the nullspace(A) from the search
853 $ direction thus the solution which is a linear combination of the search directions has no component in the nullspace(A).
854 $
855 $ We recommend always using GMRES for such singular systems.
856 $ If nullspace(A) = nullspace(A') (note symmetric matrices always satisfy this property) then both left and right preconditioning will work
857 $ If nullspace(A) != nullspace(A') then left preconditioning will work but right preconditioning may not work (or it may).
858
859 Developer Note: The reason we cannot always solve nullspace(A) != nullspace(A') systems with right preconditioning is because we need to remove at each iteration
860 the nullspace(AB) from the search direction. While we know the nullspace(A) the nullspace(AB) equals B^-1 times the nullspace(A) but except for trivial preconditioners
861 such as diagonal scaling we cannot apply the inverse of the preconditioner to a vector and thus cannot compute the nullspace(AB).
862
863
864 If using a direct method (e.g., via the KSP solver
865 KSPPREONLY and a preconditioner such as PCLU/PCILU),
866 then its=1. See KSPSetTolerances() and KSPConvergedDefault()
867 for more details.
868
869 Understanding Convergence:
870 The routines KSPMonitorSet(), KSPComputeEigenvalues(), and
871 KSPComputeEigenvaluesExplicitly() provide information on additional
872 options to monitor convergence and print eigenvalue information.
873
874 Level: beginner
875
876 .seealso: KSPCreate(), KSPSetUp(), KSPDestroy(), KSPSetTolerances(), KSPConvergedDefault(),
877 KSPSolveTranspose(), KSPGetIterationNumber(), MatNullSpaceCreate(), MatSetNullSpace(), MatSetTransposeNullSpace(), KSP,
878 KSPConvergedReasonView()
879 @*/
KSPSolve(KSP ksp,Vec b,Vec x)880 PetscErrorCode KSPSolve(KSP ksp,Vec b,Vec x)
881 {
882 PetscErrorCode ierr;
883
884 PetscFunctionBegin;
885 PetscValidHeaderSpecific(ksp,KSP_CLASSID,1);
886 if (b) PetscValidHeaderSpecific(b,VEC_CLASSID,2);
887 if (x) PetscValidHeaderSpecific(x,VEC_CLASSID,3);
888 ksp->transpose_solve = PETSC_FALSE;
889 ierr = KSPSolve_Private(ksp,b,x);CHKERRQ(ierr);
890 PetscFunctionReturn(0);
891 }
892
893 /*@
894 KSPSolveTranspose - Solves the transpose of a linear system.
895
896 Collective on ksp
897
898 Input Parameters:
899 + ksp - iterative context obtained from KSPCreate()
900 . b - right hand side vector
901 - x - solution vector
902
903 Notes:
904 For complex numbers this solve the non-Hermitian transpose system.
905
906 Developer Notes:
907 We need to implement a KSPSolveHermitianTranspose()
908
909 Level: developer
910
911 .seealso: KSPCreate(), KSPSetUp(), KSPDestroy(), KSPSetTolerances(), KSPConvergedDefault(),
912 KSPSolve(), KSP
913 @*/
KSPSolveTranspose(KSP ksp,Vec b,Vec x)914 PetscErrorCode KSPSolveTranspose(KSP ksp,Vec b,Vec x)
915 {
916 PetscErrorCode ierr;
917
918 PetscFunctionBegin;
919 PetscValidHeaderSpecific(ksp,KSP_CLASSID,1);
920 if (b) PetscValidHeaderSpecific(b,VEC_CLASSID,2);
921 if (x) PetscValidHeaderSpecific(x,VEC_CLASSID,3);
922 ksp->transpose_solve = PETSC_TRUE;
923 ierr = KSPSolve_Private(ksp,b,x);CHKERRQ(ierr);
924 PetscFunctionReturn(0);
925 }
926
KSPViewFinalMatResidual_Internal(KSP ksp,Mat B,Mat X,PetscViewer viewer,PetscViewerFormat format,PetscInt shift)927 static PetscErrorCode KSPViewFinalMatResidual_Internal(KSP ksp, Mat B, Mat X, PetscViewer viewer, PetscViewerFormat format, PetscInt shift)
928 {
929 Mat A, R;
930 PetscReal *norms;
931 PetscInt i, N;
932 PetscBool flg;
933 PetscErrorCode ierr;
934
935 PetscFunctionBegin;
936 ierr = PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &flg);CHKERRQ(ierr);
937 if (flg) {
938 ierr = PCGetOperators(ksp->pc, &A, NULL);CHKERRQ(ierr);
939 ierr = MatMatMult(A, X, MAT_INITIAL_MATRIX, PETSC_DEFAULT, &R);CHKERRQ(ierr);
940 ierr = MatAYPX(R, -1.0, B, SAME_NONZERO_PATTERN);
941 ierr = MatGetSize(R, NULL, &N);
942 ierr = PetscMalloc1(N, &norms);CHKERRQ(ierr);
943 ierr = MatGetColumnNorms(R, NORM_2, norms);CHKERRQ(ierr);
944 ierr = MatDestroy(&R);CHKERRQ(ierr);
945 for (i = 0; i < N; ++i) {
946 ierr = PetscViewerASCIIPrintf(viewer, "%s #%D %g\n", i == 0 ? "KSP final norm of residual" : " ", shift + i, (double)norms[i]);CHKERRQ(ierr);
947 }
948 ierr = PetscFree(norms);CHKERRQ(ierr);
949 }
950 PetscFunctionReturn(0);
951 }
952
953 /*@
954 KSPMatSolve - Solves a linear system with multiple right-hand sides stored as a MATDENSE. Unlike KSPSolve(), B and X must be different matrices.
955
956 Input Parameters:
957 + ksp - iterative context
958 - B - block of right-hand sides
959
960 Output Parameter:
961 . X - block of solutions
962
963 Notes:
964 This is a stripped-down version of KSPSolve(), which only handles -ksp_view, -ksp_converged_reason, and -ksp_view_final_residual.
965
966 Level: intermediate
967
968 .seealso: KSPSolve(), MatMatSolve(), MATDENSE, KSPHPDDM, PCBJACOBI, PCASM
969 @*/
KSPMatSolve(KSP ksp,Mat B,Mat X)970 PetscErrorCode KSPMatSolve(KSP ksp, Mat B, Mat X)
971 {
972 Mat A, vB, vX;
973 Vec cb, cx;
974 PetscInt m1, M1, m2, M2, n1, N1, n2, N2, Bbn = PETSC_DECIDE;
975 PetscBool match;
976 PetscErrorCode ierr;
977
978 PetscFunctionBegin;
979 PetscValidHeaderSpecific(ksp, KSP_CLASSID, 1);
980 PetscValidHeaderSpecific(B, MAT_CLASSID, 2);
981 PetscValidHeaderSpecific(X, MAT_CLASSID, 3);
982 PetscCheckSameComm(ksp, 1, B, 2);
983 PetscCheckSameComm(ksp, 1, X, 3);
984 if (!B->assembled) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Not for unassembled matrix");
985 MatCheckPreallocated(X, 3);
986 if (!X->assembled) {
987 ierr = MatSetOption(X, MAT_NO_OFF_PROC_ENTRIES, PETSC_TRUE);CHKERRQ(ierr);
988 ierr = MatAssemblyBegin(X, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
989 ierr = MatAssemblyEnd(X, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
990 }
991 if (B == X) SETERRQ(PetscObjectComm((PetscObject)ksp), PETSC_ERR_ARG_IDN, "B and X must be different matrices");
992 ierr = KSPGetOperators(ksp, &A, NULL);CHKERRQ(ierr);
993 ierr = MatGetLocalSize(A, &m1, NULL);CHKERRQ(ierr);
994 ierr = MatGetLocalSize(B, &m2, &n2);CHKERRQ(ierr);
995 ierr = MatGetSize(A, &M1, NULL);CHKERRQ(ierr);
996 ierr = MatGetSize(B, &M2, &N2);CHKERRQ(ierr);
997 if (m1 != m2 || M1 != M2) SETERRQ4(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Cannot use a block of right-hand sides with (m2,M2) = (%D,%D) for a linear system with (m1,M1) = (%D,%D)", m2, M2, m1, M1);
998 ierr = MatGetLocalSize(X, &m1, &n1);CHKERRQ(ierr);
999 ierr = MatGetSize(X, &M1, &N1);CHKERRQ(ierr);
1000 if (m1 != m2 || M1 != M2 || n1 != n2 || N1 != N2) SETERRQ8(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Incompatible block of right-hand sides (m2,M2)x(n2,N2) = (%D,%D)x(%D,%D) and solutions (m1,M1)x(n1,N1) = (%D,%D)x(%D,%D)", m2, M2, n2, N2, m1, M1, n1, N1);
1001 ierr = PetscObjectBaseTypeCompareAny((PetscObject)B, &match, MATSEQDENSE, MATMPIDENSE, "");CHKERRQ(ierr);
1002 if (!match) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Provided block of right-hand sides not stored in a dense Mat");
1003 ierr = PetscObjectBaseTypeCompareAny((PetscObject)X, &match, MATSEQDENSE, MATMPIDENSE, "");CHKERRQ(ierr);
1004 if (!match) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Provided block of solutions not stored in a dense Mat");
1005 ierr = KSPSetUp(ksp);CHKERRQ(ierr);
1006 ierr = KSPSetUpOnBlocks(ksp);CHKERRQ(ierr);
1007 if (ksp->ops->matsolve) {
1008 if (ksp->guess_zero) {
1009 ierr = MatZeroEntries(X);CHKERRQ(ierr);
1010 }
1011 ierr = PetscLogEventBegin(KSP_MatSolve, ksp, B, X, 0);CHKERRQ(ierr);
1012 ierr = KSPGetMatSolveBlockSize(ksp, &Bbn);CHKERRQ(ierr);
1013 /* by default, do a single solve with all columns */
1014 if (Bbn == PETSC_DECIDE) Bbn = N2;
1015 else if (Bbn < 1) SETERRQ1(PetscObjectComm((PetscObject)ksp), PETSC_ERR_ARG_OUTOFRANGE, "KSPMatSolve() block size %D must be positive", Bbn);
1016 ierr = PetscInfo2(ksp, "KSP type %s solving using blocks of width at most %D\n", ((PetscObject)ksp)->type_name, Bbn);CHKERRQ(ierr);
1017 /* if -ksp_matsolve_block_size is greater than the actual number of columns, do a single solve with all columns */
1018 if (Bbn >= N2) {
1019 ierr = (*ksp->ops->matsolve)(ksp, B, X);CHKERRQ(ierr);
1020 if (ksp->viewFinalRes) {
1021 ierr = KSPViewFinalMatResidual_Internal(ksp, B, X, ksp->viewerFinalRes, ksp->formatFinalRes, 0);CHKERRQ(ierr);
1022 }
1023 if (ksp->viewReason) {
1024 ierr = PetscViewerPushFormat(ksp->viewerReason,PETSC_VIEWER_DEFAULT);CHKERRQ(ierr);
1025 ierr = KSPConvergedReasonView(ksp, ksp->viewerReason);CHKERRQ(ierr);
1026 ierr = PetscViewerPopFormat(ksp->viewerReason);CHKERRQ(ierr);
1027 }
1028 } else {
1029 for (n2 = 0; n2 < N2; n2 += Bbn) {
1030 ierr = MatDenseGetSubMatrix(B, n2, PetscMin(n2+Bbn, N2), &vB);CHKERRQ(ierr);
1031 ierr = MatDenseGetSubMatrix(X, n2, PetscMin(n2+Bbn, N2), &vX);CHKERRQ(ierr);
1032 ierr = (*ksp->ops->matsolve)(ksp, vB, vX);CHKERRQ(ierr);
1033 if (ksp->viewFinalRes) {
1034 ierr = KSPViewFinalMatResidual_Internal(ksp, vB, vX, ksp->viewerFinalRes, ksp->formatFinalRes, n2);CHKERRQ(ierr);
1035 }
1036 if (ksp->viewReason) {
1037 ierr = PetscViewerPushFormat(ksp->viewerReason,PETSC_VIEWER_DEFAULT);CHKERRQ(ierr);
1038 ierr = KSPConvergedReasonView(ksp, ksp->viewerReason);CHKERRQ(ierr);
1039 ierr = PetscViewerPopFormat(ksp->viewerReason);CHKERRQ(ierr);
1040 }
1041 ierr = MatDenseRestoreSubMatrix(B, &vB);CHKERRQ(ierr);
1042 ierr = MatDenseRestoreSubMatrix(X, &vX);CHKERRQ(ierr);
1043 }
1044 }
1045 if (ksp->view) {
1046 ierr = KSPView(ksp, ksp->viewer);CHKERRQ(ierr);
1047 }
1048 ierr = PetscLogEventEnd(KSP_MatSolve, ksp, B, X, 0);CHKERRQ(ierr);
1049 } else {
1050 ierr = PetscInfo1(ksp, "KSP type %s solving column by column\n", ((PetscObject)ksp)->type_name);CHKERRQ(ierr);
1051 for (n2 = 0; n2 < N2; ++n2) {
1052 ierr = MatDenseGetColumnVecRead(B, n2, &cb);CHKERRQ(ierr);
1053 ierr = MatDenseGetColumnVecWrite(X, n2, &cx);CHKERRQ(ierr);
1054 ierr = KSPSolve(ksp, cb, cx);CHKERRQ(ierr);
1055 ierr = MatDenseRestoreColumnVecWrite(X, n2, &cx);CHKERRQ(ierr);
1056 ierr = MatDenseRestoreColumnVecRead(B, n2, &cb);CHKERRQ(ierr);
1057 }
1058 }
1059 PetscFunctionReturn(0);
1060 }
1061
1062 /*@
1063 KSPSetMatSolveBlockSize - Sets the maximum number of columns treated simultaneously in KSPMatSolve().
1064
1065 Logically collective
1066
1067 Input Parameters:
1068 + ksp - iterative context
1069 - bs - block size
1070
1071 Level: advanced
1072
1073 .seealso: KSPMatSolve(), KSPGetMatSolveBlockSize(), -mat_mumps_icntl_27, -matmatmult_Bbn
1074 @*/
KSPSetMatSolveBlockSize(KSP ksp,PetscInt bs)1075 PetscErrorCode KSPSetMatSolveBlockSize(KSP ksp, PetscInt bs)
1076 {
1077 PetscErrorCode ierr;
1078
1079 PetscFunctionBegin;
1080 PetscValidHeaderSpecific(ksp, KSP_CLASSID, 1);
1081 PetscValidLogicalCollectiveInt(ksp, bs, 2);
1082 ierr = PetscTryMethod(ksp, "KSPSetMatSolveBlockSize_C", (KSP, PetscInt), (ksp, bs));CHKERRQ(ierr);
1083 PetscFunctionReturn(0);
1084 }
1085
1086 /*@
1087 KSPGetMatSolveBlockSize - Gets the maximum number of columns treated simultaneously in KSPMatSolve().
1088
1089 Input Parameter:
1090 . ksp - iterative context
1091
1092 Output Parameter:
1093 . bs - block size
1094
1095 Level: advanced
1096
1097 .seealso: KSPMatSolve(), KSPSetMatSolveBlockSize(), -mat_mumps_icntl_27, -matmatmult_Bbn
1098 @*/
KSPGetMatSolveBlockSize(KSP ksp,PetscInt * bs)1099 PetscErrorCode KSPGetMatSolveBlockSize(KSP ksp, PetscInt *bs)
1100 {
1101 PetscErrorCode ierr;
1102
1103 PetscFunctionBegin;
1104 PetscValidHeaderSpecific(ksp, KSP_CLASSID, 1);
1105 *bs = PETSC_DECIDE;
1106 ierr = PetscTryMethod(ksp, "KSPGetMatSolveBlockSize_C", (KSP, PetscInt*), (ksp, bs));CHKERRQ(ierr);
1107 PetscFunctionReturn(0);
1108 }
1109
1110 /*@
1111 KSPResetViewers - Resets all the viewers set from the options database during KSPSetFromOptions()
1112
1113 Collective on ksp
1114
1115 Input Parameter:
1116 . ksp - iterative context obtained from KSPCreate()
1117
1118 Level: beginner
1119
1120 .seealso: KSPCreate(), KSPSetUp(), KSPSolve(), KSPSetFromOptions(), KSP
1121 @*/
KSPResetViewers(KSP ksp)1122 PetscErrorCode KSPResetViewers(KSP ksp)
1123 {
1124 PetscErrorCode ierr;
1125
1126 PetscFunctionBegin;
1127 if (ksp) PetscValidHeaderSpecific(ksp,KSP_CLASSID,1);
1128 if (!ksp) PetscFunctionReturn(0);
1129 ierr = PetscViewerDestroy(&ksp->viewer);CHKERRQ(ierr);
1130 ierr = PetscViewerDestroy(&ksp->viewerPre);CHKERRQ(ierr);
1131 ierr = PetscViewerDestroy(&ksp->viewerReason);CHKERRQ(ierr);
1132 ierr = PetscViewerDestroy(&ksp->viewerMat);CHKERRQ(ierr);
1133 ierr = PetscViewerDestroy(&ksp->viewerPMat);CHKERRQ(ierr);
1134 ierr = PetscViewerDestroy(&ksp->viewerRhs);CHKERRQ(ierr);
1135 ierr = PetscViewerDestroy(&ksp->viewerSol);CHKERRQ(ierr);
1136 ierr = PetscViewerDestroy(&ksp->viewerMatExp);CHKERRQ(ierr);
1137 ierr = PetscViewerDestroy(&ksp->viewerEV);CHKERRQ(ierr);
1138 ierr = PetscViewerDestroy(&ksp->viewerSV);CHKERRQ(ierr);
1139 ierr = PetscViewerDestroy(&ksp->viewerEVExp);CHKERRQ(ierr);
1140 ierr = PetscViewerDestroy(&ksp->viewerFinalRes);CHKERRQ(ierr);
1141 ierr = PetscViewerDestroy(&ksp->viewerPOpExp);CHKERRQ(ierr);
1142 ierr = PetscViewerDestroy(&ksp->viewerDScale);CHKERRQ(ierr);
1143 ksp->view = PETSC_FALSE;
1144 ksp->viewPre = PETSC_FALSE;
1145 ksp->viewReason = PETSC_FALSE;
1146 ksp->viewMat = PETSC_FALSE;
1147 ksp->viewPMat = PETSC_FALSE;
1148 ksp->viewRhs = PETSC_FALSE;
1149 ksp->viewSol = PETSC_FALSE;
1150 ksp->viewMatExp = PETSC_FALSE;
1151 ksp->viewEV = PETSC_FALSE;
1152 ksp->viewSV = PETSC_FALSE;
1153 ksp->viewEVExp = PETSC_FALSE;
1154 ksp->viewFinalRes = PETSC_FALSE;
1155 ksp->viewPOpExp = PETSC_FALSE;
1156 ksp->viewDScale = PETSC_FALSE;
1157 PetscFunctionReturn(0);
1158 }
1159
1160 /*@
1161 KSPReset - Resets a KSP context to the kspsetupcalled = 0 state and removes any allocated Vecs and Mats
1162
1163 Collective on ksp
1164
1165 Input Parameter:
1166 . ksp - iterative context obtained from KSPCreate()
1167
1168 Level: beginner
1169
1170 .seealso: KSPCreate(), KSPSetUp(), KSPSolve(), KSP
1171 @*/
KSPReset(KSP ksp)1172 PetscErrorCode KSPReset(KSP ksp)
1173 {
1174 PetscErrorCode ierr;
1175
1176 PetscFunctionBegin;
1177 if (ksp) PetscValidHeaderSpecific(ksp,KSP_CLASSID,1);
1178 if (!ksp) PetscFunctionReturn(0);
1179 if (ksp->ops->reset) {
1180 ierr = (*ksp->ops->reset)(ksp);CHKERRQ(ierr);
1181 }
1182 if (ksp->pc) {ierr = PCReset(ksp->pc);CHKERRQ(ierr);}
1183 if (ksp->guess) {
1184 KSPGuess guess = ksp->guess;
1185 if (guess->ops->reset) { ierr = (*guess->ops->reset)(guess);CHKERRQ(ierr); }
1186 }
1187 ierr = VecDestroyVecs(ksp->nwork,&ksp->work);CHKERRQ(ierr);
1188 ierr = VecDestroy(&ksp->vec_rhs);CHKERRQ(ierr);
1189 ierr = VecDestroy(&ksp->vec_sol);CHKERRQ(ierr);
1190 ierr = VecDestroy(&ksp->diagonal);CHKERRQ(ierr);
1191 ierr = VecDestroy(&ksp->truediagonal);CHKERRQ(ierr);
1192
1193 ierr = KSPResetViewers(ksp);CHKERRQ(ierr);
1194
1195 ksp->setupstage = KSP_SETUP_NEW;
1196 PetscFunctionReturn(0);
1197 }
1198
1199 /*@
1200 KSPDestroy - Destroys KSP context.
1201
1202 Collective on ksp
1203
1204 Input Parameter:
1205 . ksp - iterative context obtained from KSPCreate()
1206
1207 Level: beginner
1208
1209 .seealso: KSPCreate(), KSPSetUp(), KSPSolve(), KSP
1210 @*/
KSPDestroy(KSP * ksp)1211 PetscErrorCode KSPDestroy(KSP *ksp)
1212 {
1213 PetscErrorCode ierr;
1214 PC pc;
1215
1216 PetscFunctionBegin;
1217 if (!*ksp) PetscFunctionReturn(0);
1218 PetscValidHeaderSpecific((*ksp),KSP_CLASSID,1);
1219 if (--((PetscObject)(*ksp))->refct > 0) {*ksp = NULL; PetscFunctionReturn(0);}
1220
1221 ierr = PetscObjectSAWsViewOff((PetscObject)*ksp);CHKERRQ(ierr);
1222
1223 /*
1224 Avoid a cascading call to PCReset(ksp->pc) from the following call:
1225 PCReset() shouldn't be called from KSPDestroy() as it is unprotected by pc's
1226 refcount (and may be shared, e.g., by other ksps).
1227 */
1228 pc = (*ksp)->pc;
1229 (*ksp)->pc = NULL;
1230 ierr = KSPReset((*ksp));CHKERRQ(ierr);
1231 (*ksp)->pc = pc;
1232 if ((*ksp)->ops->destroy) {ierr = (*(*ksp)->ops->destroy)(*ksp);CHKERRQ(ierr);}
1233
1234 ierr = KSPGuessDestroy(&(*ksp)->guess);CHKERRQ(ierr);
1235 ierr = DMDestroy(&(*ksp)->dm);CHKERRQ(ierr);
1236 ierr = PCDestroy(&(*ksp)->pc);CHKERRQ(ierr);
1237 ierr = PetscFree((*ksp)->res_hist_alloc);CHKERRQ(ierr);
1238 if ((*ksp)->convergeddestroy) {
1239 ierr = (*(*ksp)->convergeddestroy)((*ksp)->cnvP);CHKERRQ(ierr);
1240 }
1241 ierr = KSPMonitorCancel((*ksp));CHKERRQ(ierr);
1242 ierr = PetscViewerDestroy(&(*ksp)->eigviewer);CHKERRQ(ierr);
1243 ierr = PetscHeaderDestroy(ksp);CHKERRQ(ierr);
1244 PetscFunctionReturn(0);
1245 }
1246
1247 /*@
1248 KSPSetPCSide - Sets the preconditioning side.
1249
1250 Logically Collective on ksp
1251
1252 Input Parameter:
1253 . ksp - iterative context obtained from KSPCreate()
1254
1255 Output Parameter:
1256 . side - the preconditioning side, where side is one of
1257 .vb
1258 PC_LEFT - left preconditioning (default)
1259 PC_RIGHT - right preconditioning
1260 PC_SYMMETRIC - symmetric preconditioning
1261 .ve
1262
1263 Options Database Keys:
1264 . -ksp_pc_side <right,left,symmetric>
1265
1266 Notes:
1267 Left preconditioning is used by default for most Krylov methods except KSPFGMRES which only supports right preconditioning.
1268
1269 For methods changing the side of the preconditioner changes the norm type that is used, see KSPSetNormType().
1270
1271 Symmetric preconditioning is currently available only for the KSPQCG method. Note, however, that
1272 symmetric preconditioning can be emulated by using either right or left
1273 preconditioning and a pre or post processing step.
1274
1275 Setting the PC side often affects the default norm type. See KSPSetNormType() for details.
1276
1277 Level: intermediate
1278
1279 .seealso: KSPGetPCSide(), KSPSetNormType(), KSPGetNormType(), KSP
1280 @*/
KSPSetPCSide(KSP ksp,PCSide side)1281 PetscErrorCode KSPSetPCSide(KSP ksp,PCSide side)
1282 {
1283 PetscFunctionBegin;
1284 PetscValidHeaderSpecific(ksp,KSP_CLASSID,1);
1285 PetscValidLogicalCollectiveEnum(ksp,side,2);
1286 ksp->pc_side = ksp->pc_side_set = side;
1287 PetscFunctionReturn(0);
1288 }
1289
1290 /*@
1291 KSPGetPCSide - Gets the preconditioning side.
1292
1293 Not Collective
1294
1295 Input Parameter:
1296 . ksp - iterative context obtained from KSPCreate()
1297
1298 Output Parameter:
1299 . side - the preconditioning side, where side is one of
1300 .vb
1301 PC_LEFT - left preconditioning (default)
1302 PC_RIGHT - right preconditioning
1303 PC_SYMMETRIC - symmetric preconditioning
1304 .ve
1305
1306 Level: intermediate
1307
1308 .seealso: KSPSetPCSide(), KSP
1309 @*/
KSPGetPCSide(KSP ksp,PCSide * side)1310 PetscErrorCode KSPGetPCSide(KSP ksp,PCSide *side)
1311 {
1312 PetscErrorCode ierr;
1313
1314 PetscFunctionBegin;
1315 PetscValidHeaderSpecific(ksp,KSP_CLASSID,1);
1316 PetscValidPointer(side,2);
1317 ierr = KSPSetUpNorms_Private(ksp,PETSC_TRUE,&ksp->normtype,&ksp->pc_side);CHKERRQ(ierr);
1318 *side = ksp->pc_side;
1319 PetscFunctionReturn(0);
1320 }
1321
1322 /*@
1323 KSPGetTolerances - Gets the relative, absolute, divergence, and maximum
1324 iteration tolerances used by the default KSP convergence tests.
1325
1326 Not Collective
1327
1328 Input Parameter:
1329 . ksp - the Krylov subspace context
1330
1331 Output Parameters:
1332 + rtol - the relative convergence tolerance
1333 . abstol - the absolute convergence tolerance
1334 . dtol - the divergence tolerance
1335 - maxits - maximum number of iterations
1336
1337 Notes:
1338 The user can specify NULL for any parameter that is not needed.
1339
1340 Level: intermediate
1341
1342 maximum, iterations
1343
1344 .seealso: KSPSetTolerances(), KSP
1345 @*/
KSPGetTolerances(KSP ksp,PetscReal * rtol,PetscReal * abstol,PetscReal * dtol,PetscInt * maxits)1346 PetscErrorCode KSPGetTolerances(KSP ksp,PetscReal *rtol,PetscReal *abstol,PetscReal *dtol,PetscInt *maxits)
1347 {
1348 PetscFunctionBegin;
1349 PetscValidHeaderSpecific(ksp,KSP_CLASSID,1);
1350 if (abstol) *abstol = ksp->abstol;
1351 if (rtol) *rtol = ksp->rtol;
1352 if (dtol) *dtol = ksp->divtol;
1353 if (maxits) *maxits = ksp->max_it;
1354 PetscFunctionReturn(0);
1355 }
1356
1357 /*@
1358 KSPSetTolerances - Sets the relative, absolute, divergence, and maximum
1359 iteration tolerances used by the default KSP convergence testers.
1360
1361 Logically Collective on ksp
1362
1363 Input Parameters:
1364 + ksp - the Krylov subspace context
1365 . rtol - the relative convergence tolerance, relative decrease in the (possibly preconditioned) residual norm
1366 . abstol - the absolute convergence tolerance absolute size of the (possibly preconditioned) residual norm
1367 . dtol - the divergence tolerance, amount (possibly preconditioned) residual norm can increase before KSPConvergedDefault() concludes that the method is diverging
1368 - maxits - maximum number of iterations to use
1369
1370 Options Database Keys:
1371 + -ksp_atol <abstol> - Sets abstol
1372 . -ksp_rtol <rtol> - Sets rtol
1373 . -ksp_divtol <dtol> - Sets dtol
1374 - -ksp_max_it <maxits> - Sets maxits
1375
1376 Notes:
1377 Use PETSC_DEFAULT to retain the default value of any of the tolerances.
1378
1379 See KSPConvergedDefault() for details how these parameters are used in the default convergence test. See also KSPSetConvergenceTest()
1380 for setting user-defined stopping criteria.
1381
1382 Level: intermediate
1383
1384 convergence, maximum, iterations
1385
1386 .seealso: KSPGetTolerances(), KSPConvergedDefault(), KSPSetConvergenceTest(), KSP
1387 @*/
KSPSetTolerances(KSP ksp,PetscReal rtol,PetscReal abstol,PetscReal dtol,PetscInt maxits)1388 PetscErrorCode KSPSetTolerances(KSP ksp,PetscReal rtol,PetscReal abstol,PetscReal dtol,PetscInt maxits)
1389 {
1390 PetscFunctionBegin;
1391 PetscValidHeaderSpecific(ksp,KSP_CLASSID,1);
1392 PetscValidLogicalCollectiveReal(ksp,rtol,2);
1393 PetscValidLogicalCollectiveReal(ksp,abstol,3);
1394 PetscValidLogicalCollectiveReal(ksp,dtol,4);
1395 PetscValidLogicalCollectiveInt(ksp,maxits,5);
1396
1397 if (rtol != PETSC_DEFAULT) {
1398 if (rtol < 0.0 || 1.0 <= rtol) SETERRQ1(PetscObjectComm((PetscObject)ksp),PETSC_ERR_ARG_OUTOFRANGE,"Relative tolerance %g must be non-negative and less than 1.0",(double)rtol);
1399 ksp->rtol = rtol;
1400 }
1401 if (abstol != PETSC_DEFAULT) {
1402 if (abstol < 0.0) SETERRQ1(PetscObjectComm((PetscObject)ksp),PETSC_ERR_ARG_OUTOFRANGE,"Absolute tolerance %g must be non-negative",(double)abstol);
1403 ksp->abstol = abstol;
1404 }
1405 if (dtol != PETSC_DEFAULT) {
1406 if (dtol < 0.0) SETERRQ1(PetscObjectComm((PetscObject)ksp),PETSC_ERR_ARG_OUTOFRANGE,"Divergence tolerance %g must be larger than 1.0",(double)dtol);
1407 ksp->divtol = dtol;
1408 }
1409 if (maxits != PETSC_DEFAULT) {
1410 if (maxits < 0) SETERRQ1(PetscObjectComm((PetscObject)ksp),PETSC_ERR_ARG_OUTOFRANGE,"Maximum number of iterations %D must be non-negative",maxits);
1411 ksp->max_it = maxits;
1412 }
1413 PetscFunctionReturn(0);
1414 }
1415
1416 /*@
1417 KSPSetInitialGuessNonzero - Tells the iterative solver that the
1418 initial guess is nonzero; otherwise KSP assumes the initial guess
1419 is to be zero (and thus zeros it out before solving).
1420
1421 Logically Collective on ksp
1422
1423 Input Parameters:
1424 + ksp - iterative context obtained from KSPCreate()
1425 - flg - PETSC_TRUE indicates the guess is non-zero, PETSC_FALSE indicates the guess is zero
1426
1427 Options database keys:
1428 . -ksp_initial_guess_nonzero : use nonzero initial guess; this takes an optional truth value (0/1/no/yes/true/false)
1429
1430 Level: beginner
1431
1432 Notes:
1433 If this is not called the X vector is zeroed in the call to KSPSolve().
1434
1435 .seealso: KSPGetInitialGuessNonzero(), KSPSetGuessType(), KSPGuessType, KSP
1436 @*/
KSPSetInitialGuessNonzero(KSP ksp,PetscBool flg)1437 PetscErrorCode KSPSetInitialGuessNonzero(KSP ksp,PetscBool flg)
1438 {
1439 PetscFunctionBegin;
1440 PetscValidHeaderSpecific(ksp,KSP_CLASSID,1);
1441 PetscValidLogicalCollectiveBool(ksp,flg,2);
1442 ksp->guess_zero = (PetscBool) !(int)flg;
1443 PetscFunctionReturn(0);
1444 }
1445
1446 /*@
1447 KSPGetInitialGuessNonzero - Determines whether the KSP solver is using
1448 a zero initial guess.
1449
1450 Not Collective
1451
1452 Input Parameter:
1453 . ksp - iterative context obtained from KSPCreate()
1454
1455 Output Parameter:
1456 . flag - PETSC_TRUE if guess is nonzero, else PETSC_FALSE
1457
1458 Level: intermediate
1459
1460 .seealso: KSPSetInitialGuessNonzero(), KSP
1461 @*/
KSPGetInitialGuessNonzero(KSP ksp,PetscBool * flag)1462 PetscErrorCode KSPGetInitialGuessNonzero(KSP ksp,PetscBool *flag)
1463 {
1464 PetscFunctionBegin;
1465 PetscValidHeaderSpecific(ksp,KSP_CLASSID,1);
1466 PetscValidBoolPointer(flag,2);
1467 if (ksp->guess_zero) *flag = PETSC_FALSE;
1468 else *flag = PETSC_TRUE;
1469 PetscFunctionReturn(0);
1470 }
1471
1472 /*@
1473 KSPSetErrorIfNotConverged - Causes KSPSolve() to generate an error if the solver has not converged.
1474
1475 Logically Collective on ksp
1476
1477 Input Parameters:
1478 + ksp - iterative context obtained from KSPCreate()
1479 - flg - PETSC_TRUE indicates you want the error generated
1480
1481 Options database keys:
1482 . -ksp_error_if_not_converged : this takes an optional truth value (0/1/no/yes/true/false)
1483
1484 Level: intermediate
1485
1486 Notes:
1487 Normally PETSc continues if a linear solver fails to converge, you can call KSPGetConvergedReason() after a KSPSolve()
1488 to determine if it has converged.
1489
1490
1491 .seealso: KSPGetErrorIfNotConverged(), KSP
1492 @*/
KSPSetErrorIfNotConverged(KSP ksp,PetscBool flg)1493 PetscErrorCode KSPSetErrorIfNotConverged(KSP ksp,PetscBool flg)
1494 {
1495 PetscFunctionBegin;
1496 PetscValidHeaderSpecific(ksp,KSP_CLASSID,1);
1497 PetscValidLogicalCollectiveBool(ksp,flg,2);
1498 ksp->errorifnotconverged = flg;
1499 PetscFunctionReturn(0);
1500 }
1501
1502 /*@
1503 KSPGetErrorIfNotConverged - Will KSPSolve() generate an error if the solver does not converge?
1504
1505 Not Collective
1506
1507 Input Parameter:
1508 . ksp - iterative context obtained from KSPCreate()
1509
1510 Output Parameter:
1511 . flag - PETSC_TRUE if it will generate an error, else PETSC_FALSE
1512
1513 Level: intermediate
1514
1515 .seealso: KSPSetErrorIfNotConverged(), KSP
1516 @*/
KSPGetErrorIfNotConverged(KSP ksp,PetscBool * flag)1517 PetscErrorCode KSPGetErrorIfNotConverged(KSP ksp,PetscBool *flag)
1518 {
1519 PetscFunctionBegin;
1520 PetscValidHeaderSpecific(ksp,KSP_CLASSID,1);
1521 PetscValidBoolPointer(flag,2);
1522 *flag = ksp->errorifnotconverged;
1523 PetscFunctionReturn(0);
1524 }
1525
1526 /*@
1527 KSPSetInitialGuessKnoll - Tells the iterative solver to use PCApply(pc,b,..) to compute the initial guess (The Knoll trick)
1528
1529 Logically Collective on ksp
1530
1531 Input Parameters:
1532 + ksp - iterative context obtained from KSPCreate()
1533 - flg - PETSC_TRUE or PETSC_FALSE
1534
1535 Level: advanced
1536
1537 Developer Note: the Knoll trick is not currently implemented using the KSPGuess class
1538
1539 .seealso: KSPGetInitialGuessKnoll(), KSPSetInitialGuessNonzero(), KSPGetInitialGuessNonzero(), KSP
1540 @*/
KSPSetInitialGuessKnoll(KSP ksp,PetscBool flg)1541 PetscErrorCode KSPSetInitialGuessKnoll(KSP ksp,PetscBool flg)
1542 {
1543 PetscFunctionBegin;
1544 PetscValidHeaderSpecific(ksp,KSP_CLASSID,1);
1545 PetscValidLogicalCollectiveBool(ksp,flg,2);
1546 ksp->guess_knoll = flg;
1547 PetscFunctionReturn(0);
1548 }
1549
1550 /*@
1551 KSPGetInitialGuessKnoll - Determines whether the KSP solver is using the Knoll trick (using PCApply(pc,b,...) to compute
1552 the initial guess
1553
1554 Not Collective
1555
1556 Input Parameter:
1557 . ksp - iterative context obtained from KSPCreate()
1558
1559 Output Parameter:
1560 . flag - PETSC_TRUE if using Knoll trick, else PETSC_FALSE
1561
1562 Level: advanced
1563
1564 .seealso: KSPSetInitialGuessKnoll(), KSPSetInitialGuessNonzero(), KSPGetInitialGuessNonzero(), KSP
1565 @*/
KSPGetInitialGuessKnoll(KSP ksp,PetscBool * flag)1566 PetscErrorCode KSPGetInitialGuessKnoll(KSP ksp,PetscBool *flag)
1567 {
1568 PetscFunctionBegin;
1569 PetscValidHeaderSpecific(ksp,KSP_CLASSID,1);
1570 PetscValidBoolPointer(flag,2);
1571 *flag = ksp->guess_knoll;
1572 PetscFunctionReturn(0);
1573 }
1574
1575 /*@
1576 KSPGetComputeSingularValues - Gets the flag indicating whether the extreme singular
1577 values will be calculated via a Lanczos or Arnoldi process as the linear
1578 system is solved.
1579
1580 Not Collective
1581
1582 Input Parameter:
1583 . ksp - iterative context obtained from KSPCreate()
1584
1585 Output Parameter:
1586 . flg - PETSC_TRUE or PETSC_FALSE
1587
1588 Options Database Key:
1589 . -ksp_monitor_singular_value - Activates KSPSetComputeSingularValues()
1590
1591 Notes:
1592 Currently this option is not valid for all iterative methods.
1593
1594 Many users may just want to use the monitoring routine
1595 KSPMonitorSingularValue() (which can be set with option -ksp_monitor_singular_value)
1596 to print the singular values at each iteration of the linear solve.
1597
1598 Level: advanced
1599
1600 .seealso: KSPComputeExtremeSingularValues(), KSPMonitorSingularValue(), KSP
1601 @*/
KSPGetComputeSingularValues(KSP ksp,PetscBool * flg)1602 PetscErrorCode KSPGetComputeSingularValues(KSP ksp,PetscBool *flg)
1603 {
1604 PetscFunctionBegin;
1605 PetscValidHeaderSpecific(ksp,KSP_CLASSID,1);
1606 PetscValidBoolPointer(flg,2);
1607 *flg = ksp->calc_sings;
1608 PetscFunctionReturn(0);
1609 }
1610
1611 /*@
1612 KSPSetComputeSingularValues - Sets a flag so that the extreme singular
1613 values will be calculated via a Lanczos or Arnoldi process as the linear
1614 system is solved.
1615
1616 Logically Collective on ksp
1617
1618 Input Parameters:
1619 + ksp - iterative context obtained from KSPCreate()
1620 - flg - PETSC_TRUE or PETSC_FALSE
1621
1622 Options Database Key:
1623 . -ksp_monitor_singular_value - Activates KSPSetComputeSingularValues()
1624
1625 Notes:
1626 Currently this option is not valid for all iterative methods.
1627
1628 Many users may just want to use the monitoring routine
1629 KSPMonitorSingularValue() (which can be set with option -ksp_monitor_singular_value)
1630 to print the singular values at each iteration of the linear solve.
1631
1632 Level: advanced
1633
1634 .seealso: KSPComputeExtremeSingularValues(), KSPMonitorSingularValue(), KSP
1635 @*/
KSPSetComputeSingularValues(KSP ksp,PetscBool flg)1636 PetscErrorCode KSPSetComputeSingularValues(KSP ksp,PetscBool flg)
1637 {
1638 PetscFunctionBegin;
1639 PetscValidHeaderSpecific(ksp,KSP_CLASSID,1);
1640 PetscValidLogicalCollectiveBool(ksp,flg,2);
1641 ksp->calc_sings = flg;
1642 PetscFunctionReturn(0);
1643 }
1644
1645 /*@
1646 KSPGetComputeEigenvalues - Gets the flag indicating that the extreme eigenvalues
1647 values will be calculated via a Lanczos or Arnoldi process as the linear
1648 system is solved.
1649
1650 Not Collective
1651
1652 Input Parameter:
1653 . ksp - iterative context obtained from KSPCreate()
1654
1655 Output Parameter:
1656 . flg - PETSC_TRUE or PETSC_FALSE
1657
1658 Notes:
1659 Currently this option is not valid for all iterative methods.
1660
1661 Level: advanced
1662
1663 .seealso: KSPComputeEigenvalues(), KSPComputeEigenvaluesExplicitly(), KSP
1664 @*/
KSPGetComputeEigenvalues(KSP ksp,PetscBool * flg)1665 PetscErrorCode KSPGetComputeEigenvalues(KSP ksp,PetscBool *flg)
1666 {
1667 PetscFunctionBegin;
1668 PetscValidHeaderSpecific(ksp,KSP_CLASSID,1);
1669 PetscValidBoolPointer(flg,2);
1670 *flg = ksp->calc_sings;
1671 PetscFunctionReturn(0);
1672 }
1673
1674 /*@
1675 KSPSetComputeEigenvalues - Sets a flag so that the extreme eigenvalues
1676 values will be calculated via a Lanczos or Arnoldi process as the linear
1677 system is solved.
1678
1679 Logically Collective on ksp
1680
1681 Input Parameters:
1682 + ksp - iterative context obtained from KSPCreate()
1683 - flg - PETSC_TRUE or PETSC_FALSE
1684
1685 Notes:
1686 Currently this option is not valid for all iterative methods.
1687
1688 Level: advanced
1689
1690 .seealso: KSPComputeEigenvalues(), KSPComputeEigenvaluesExplicitly(), KSP
1691 @*/
KSPSetComputeEigenvalues(KSP ksp,PetscBool flg)1692 PetscErrorCode KSPSetComputeEigenvalues(KSP ksp,PetscBool flg)
1693 {
1694 PetscFunctionBegin;
1695 PetscValidHeaderSpecific(ksp,KSP_CLASSID,1);
1696 PetscValidLogicalCollectiveBool(ksp,flg,2);
1697 ksp->calc_sings = flg;
1698 PetscFunctionReturn(0);
1699 }
1700
1701 /*@
1702 KSPSetComputeRitz - Sets a flag so that the Ritz or harmonic Ritz pairs
1703 will be calculated via a Lanczos or Arnoldi process as the linear
1704 system is solved.
1705
1706 Logically Collective on ksp
1707
1708 Input Parameters:
1709 + ksp - iterative context obtained from KSPCreate()
1710 - flg - PETSC_TRUE or PETSC_FALSE
1711
1712 Notes:
1713 Currently this option is only valid for the GMRES method.
1714
1715 Level: advanced
1716
1717 .seealso: KSPComputeRitz(), KSP
1718 @*/
KSPSetComputeRitz(KSP ksp,PetscBool flg)1719 PetscErrorCode KSPSetComputeRitz(KSP ksp, PetscBool flg)
1720 {
1721 PetscFunctionBegin;
1722 PetscValidHeaderSpecific(ksp,KSP_CLASSID,1);
1723 PetscValidLogicalCollectiveBool(ksp,flg,2);
1724 ksp->calc_ritz = flg;
1725 PetscFunctionReturn(0);
1726 }
1727
1728 /*@
1729 KSPGetRhs - Gets the right-hand-side vector for the linear system to
1730 be solved.
1731
1732 Not Collective
1733
1734 Input Parameter:
1735 . ksp - iterative context obtained from KSPCreate()
1736
1737 Output Parameter:
1738 . r - right-hand-side vector
1739
1740 Level: developer
1741
1742 .seealso: KSPGetSolution(), KSPSolve(), KSP
1743 @*/
KSPGetRhs(KSP ksp,Vec * r)1744 PetscErrorCode KSPGetRhs(KSP ksp,Vec *r)
1745 {
1746 PetscFunctionBegin;
1747 PetscValidHeaderSpecific(ksp,KSP_CLASSID,1);
1748 PetscValidPointer(r,2);
1749 *r = ksp->vec_rhs;
1750 PetscFunctionReturn(0);
1751 }
1752
1753 /*@
1754 KSPGetSolution - Gets the location of the solution for the
1755 linear system to be solved. Note that this may not be where the solution
1756 is stored during the iterative process; see KSPBuildSolution().
1757
1758 Not Collective
1759
1760 Input Parameters:
1761 . ksp - iterative context obtained from KSPCreate()
1762
1763 Output Parameters:
1764 . v - solution vector
1765
1766 Level: developer
1767
1768 .seealso: KSPGetRhs(), KSPBuildSolution(), KSPSolve(), KSP
1769 @*/
KSPGetSolution(KSP ksp,Vec * v)1770 PetscErrorCode KSPGetSolution(KSP ksp,Vec *v)
1771 {
1772 PetscFunctionBegin;
1773 PetscValidHeaderSpecific(ksp,KSP_CLASSID,1);
1774 PetscValidPointer(v,2);
1775 *v = ksp->vec_sol;
1776 PetscFunctionReturn(0);
1777 }
1778
1779 /*@
1780 KSPSetPC - Sets the preconditioner to be used to calculate the
1781 application of the preconditioner on a vector.
1782
1783 Collective on ksp
1784
1785 Input Parameters:
1786 + ksp - iterative context obtained from KSPCreate()
1787 - pc - the preconditioner object (can be NULL)
1788
1789 Notes:
1790 Use KSPGetPC() to retrieve the preconditioner context.
1791
1792 Level: developer
1793
1794 .seealso: KSPGetPC(), KSP
1795 @*/
KSPSetPC(KSP ksp,PC pc)1796 PetscErrorCode KSPSetPC(KSP ksp,PC pc)
1797 {
1798 PetscErrorCode ierr;
1799
1800 PetscFunctionBegin;
1801 PetscValidHeaderSpecific(ksp,KSP_CLASSID,1);
1802 if (pc) {
1803 PetscValidHeaderSpecific(pc,PC_CLASSID,2);
1804 PetscCheckSameComm(ksp,1,pc,2);
1805 }
1806 ierr = PetscObjectReference((PetscObject)pc);CHKERRQ(ierr);
1807 ierr = PCDestroy(&ksp->pc);CHKERRQ(ierr);
1808 ksp->pc = pc;
1809 ierr = PetscLogObjectParent((PetscObject)ksp,(PetscObject)ksp->pc);CHKERRQ(ierr);
1810 PetscFunctionReturn(0);
1811 }
1812
1813 /*@
1814 KSPGetPC - Returns a pointer to the preconditioner context
1815 set with KSPSetPC().
1816
1817 Not Collective
1818
1819 Input Parameters:
1820 . ksp - iterative context obtained from KSPCreate()
1821
1822 Output Parameter:
1823 . pc - preconditioner context
1824
1825 Level: developer
1826
1827 .seealso: KSPSetPC(), KSP
1828 @*/
KSPGetPC(KSP ksp,PC * pc)1829 PetscErrorCode KSPGetPC(KSP ksp,PC *pc)
1830 {
1831 PetscErrorCode ierr;
1832
1833 PetscFunctionBegin;
1834 PetscValidHeaderSpecific(ksp,KSP_CLASSID,1);
1835 PetscValidPointer(pc,2);
1836 if (!ksp->pc) {
1837 ierr = PCCreate(PetscObjectComm((PetscObject)ksp),&ksp->pc);CHKERRQ(ierr);
1838 ierr = PetscObjectIncrementTabLevel((PetscObject)ksp->pc,(PetscObject)ksp,0);CHKERRQ(ierr);
1839 ierr = PetscLogObjectParent((PetscObject)ksp,(PetscObject)ksp->pc);CHKERRQ(ierr);
1840 ierr = PetscObjectSetOptions((PetscObject)ksp->pc,((PetscObject)ksp)->options);CHKERRQ(ierr);
1841 }
1842 *pc = ksp->pc;
1843 PetscFunctionReturn(0);
1844 }
1845
1846 /*@
1847 KSPMonitor - runs the user provided monitor routines, if they exist
1848
1849 Collective on ksp
1850
1851 Input Parameters:
1852 + ksp - iterative context obtained from KSPCreate()
1853 . it - iteration number
1854 - rnorm - relative norm of the residual
1855
1856 Notes:
1857 This routine is called by the KSP implementations.
1858 It does not typically need to be called by the user.
1859
1860 Level: developer
1861
1862 .seealso: KSPMonitorSet()
1863 @*/
KSPMonitor(KSP ksp,PetscInt it,PetscReal rnorm)1864 PetscErrorCode KSPMonitor(KSP ksp,PetscInt it,PetscReal rnorm)
1865 {
1866 PetscInt i, n = ksp->numbermonitors;
1867 PetscErrorCode ierr;
1868
1869 PetscFunctionBegin;
1870 for (i=0; i<n; i++) {
1871 ierr = (*ksp->monitor[i])(ksp,it,rnorm,ksp->monitorcontext[i]);CHKERRQ(ierr);
1872 }
1873 PetscFunctionReturn(0);
1874 }
1875
1876 /*@C
1877 KSPMonitorSet - Sets an ADDITIONAL function to be called at every iteration to monitor
1878 the residual/error etc.
1879
1880 Logically Collective on ksp
1881
1882 Input Parameters:
1883 + ksp - iterative context obtained from KSPCreate()
1884 . monitor - pointer to function (if this is NULL, it turns off monitoring
1885 . mctx - [optional] context for private data for the
1886 monitor routine (use NULL if no context is desired)
1887 - monitordestroy - [optional] routine that frees monitor context
1888 (may be NULL)
1889
1890 Calling Sequence of monitor:
1891 $ monitor (KSP ksp, PetscInt it, PetscReal rnorm, void *mctx)
1892
1893 + ksp - iterative context obtained from KSPCreate()
1894 . it - iteration number
1895 . rnorm - (estimated) 2-norm of (preconditioned) residual
1896 - mctx - optional monitoring context, as set by KSPMonitorSet()
1897
1898 Options Database Keys:
1899 + -ksp_monitor - sets KSPMonitorDefault()
1900 . -ksp_monitor_true_residual - sets KSPMonitorTrueResidualNorm()
1901 . -ksp_monitor_max - sets KSPMonitorTrueResidualMaxNorm()
1902 . -ksp_monitor_lg_residualnorm - sets line graph monitor,
1903 uses KSPMonitorLGResidualNormCreate()
1904 . -ksp_monitor_lg_true_residualnorm - sets line graph monitor,
1905 uses KSPMonitorLGResidualNormCreate()
1906 . -ksp_monitor_singular_value - sets KSPMonitorSingularValue()
1907 - -ksp_monitor_cancel - cancels all monitors that have
1908 been hardwired into a code by
1909 calls to KSPMonitorSet(), but
1910 does not cancel those set via
1911 the options database.
1912
1913 Notes:
1914 The default is to do nothing. To print the residual, or preconditioned
1915 residual if KSPSetNormType(ksp,KSP_NORM_PRECONDITIONED) was called, use
1916 KSPMonitorDefault() as the monitoring routine, with a ASCII viewer as the
1917 context.
1918
1919 Several different monitoring routines may be set by calling
1920 KSPMonitorSet() multiple times; all will be called in the
1921 order in which they were set.
1922
1923 Fortran Notes:
1924 Only a single monitor function can be set for each KSP object
1925
1926 Level: beginner
1927
1928 .seealso: KSPMonitorDefault(), KSPMonitorLGResidualNormCreate(), KSPMonitorCancel(), KSP
1929 @*/
KSPMonitorSet(KSP ksp,PetscErrorCode (* monitor)(KSP,PetscInt,PetscReal,void *),void * mctx,PetscErrorCode (* monitordestroy)(void **))1930 PetscErrorCode KSPMonitorSet(KSP ksp,PetscErrorCode (*monitor)(KSP,PetscInt,PetscReal,void*),void *mctx,PetscErrorCode (*monitordestroy)(void**))
1931 {
1932 PetscInt i;
1933 PetscErrorCode ierr;
1934 PetscBool identical;
1935
1936 PetscFunctionBegin;
1937 PetscValidHeaderSpecific(ksp,KSP_CLASSID,1);
1938 for (i=0; i<ksp->numbermonitors;i++) {
1939 ierr = PetscMonitorCompare((PetscErrorCode (*)(void))monitor,mctx,monitordestroy,(PetscErrorCode (*)(void))ksp->monitor[i],ksp->monitorcontext[i],ksp->monitordestroy[i],&identical);CHKERRQ(ierr);
1940 if (identical) PetscFunctionReturn(0);
1941 }
1942 if (ksp->numbermonitors >= MAXKSPMONITORS) SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_ARG_OUTOFRANGE,"Too many KSP monitors set");
1943 ksp->monitor[ksp->numbermonitors] = monitor;
1944 ksp->monitordestroy[ksp->numbermonitors] = monitordestroy;
1945 ksp->monitorcontext[ksp->numbermonitors++] = (void*)mctx;
1946 PetscFunctionReturn(0);
1947 }
1948
1949 /*@
1950 KSPMonitorCancel - Clears all monitors for a KSP object.
1951
1952 Logically Collective on ksp
1953
1954 Input Parameters:
1955 . ksp - iterative context obtained from KSPCreate()
1956
1957 Options Database Key:
1958 . -ksp_monitor_cancel - Cancels all monitors that have
1959 been hardwired into a code by calls to KSPMonitorSet(),
1960 but does not cancel those set via the options database.
1961
1962 Level: intermediate
1963
1964 .seealso: KSPMonitorDefault(), KSPMonitorLGResidualNormCreate(), KSPMonitorSet(), KSP
1965 @*/
KSPMonitorCancel(KSP ksp)1966 PetscErrorCode KSPMonitorCancel(KSP ksp)
1967 {
1968 PetscErrorCode ierr;
1969 PetscInt i;
1970
1971 PetscFunctionBegin;
1972 PetscValidHeaderSpecific(ksp,KSP_CLASSID,1);
1973 for (i=0; i<ksp->numbermonitors; i++) {
1974 if (ksp->monitordestroy[i]) {
1975 ierr = (*ksp->monitordestroy[i])(&ksp->monitorcontext[i]);CHKERRQ(ierr);
1976 }
1977 }
1978 ksp->numbermonitors = 0;
1979 PetscFunctionReturn(0);
1980 }
1981
1982 /*@C
1983 KSPGetMonitorContext - Gets the monitoring context, as set by
1984 KSPMonitorSet() for the FIRST monitor only.
1985
1986 Not Collective
1987
1988 Input Parameter:
1989 . ksp - iterative context obtained from KSPCreate()
1990
1991 Output Parameter:
1992 . ctx - monitoring context
1993
1994 Level: intermediate
1995
1996 .seealso: KSPMonitorDefault(), KSPMonitorLGResidualNormCreate(), KSP
1997 @*/
KSPGetMonitorContext(KSP ksp,void ** ctx)1998 PetscErrorCode KSPGetMonitorContext(KSP ksp,void **ctx)
1999 {
2000 PetscFunctionBegin;
2001 PetscValidHeaderSpecific(ksp,KSP_CLASSID,1);
2002 *ctx = (ksp->monitorcontext[0]);
2003 PetscFunctionReturn(0);
2004 }
2005
2006 /*@
2007 KSPSetResidualHistory - Sets the array used to hold the residual history.
2008 If set, this array will contain the residual norms computed at each
2009 iteration of the solver.
2010
2011 Not Collective
2012
2013 Input Parameters:
2014 + ksp - iterative context obtained from KSPCreate()
2015 . a - array to hold history
2016 . na - size of a
2017 - reset - PETSC_TRUE indicates the history counter is reset to zero
2018 for each new linear solve
2019
2020 Level: advanced
2021
2022 Notes:
2023 The array is NOT freed by PETSc so the user needs to keep track of
2024 it and destroy once the KSP object is destroyed.
2025
2026 If 'a' is NULL then space is allocated for the history. If 'na' PETSC_DECIDE or PETSC_DEFAULT then a
2027 default array of length 10000 is allocated.
2028
2029 .seealso: KSPGetResidualHistory(), KSP
2030
2031 @*/
KSPSetResidualHistory(KSP ksp,PetscReal a[],PetscInt na,PetscBool reset)2032 PetscErrorCode KSPSetResidualHistory(KSP ksp,PetscReal a[],PetscInt na,PetscBool reset)
2033 {
2034 PetscErrorCode ierr;
2035
2036 PetscFunctionBegin;
2037 PetscValidHeaderSpecific(ksp,KSP_CLASSID,1);
2038
2039 ierr = PetscFree(ksp->res_hist_alloc);CHKERRQ(ierr);
2040 if (na != PETSC_DECIDE && na != PETSC_DEFAULT && a) {
2041 ksp->res_hist = a;
2042 ksp->res_hist_max = na;
2043 } else {
2044 if (na != PETSC_DECIDE && na != PETSC_DEFAULT) ksp->res_hist_max = na;
2045 else ksp->res_hist_max = 10000; /* like default ksp->max_it */
2046 ierr = PetscCalloc1(ksp->res_hist_max,&ksp->res_hist_alloc);CHKERRQ(ierr);
2047
2048 ksp->res_hist = ksp->res_hist_alloc;
2049 }
2050 ksp->res_hist_len = 0;
2051 ksp->res_hist_reset = reset;
2052 PetscFunctionReturn(0);
2053 }
2054
2055 /*@C
2056 KSPGetResidualHistory - Gets the array used to hold the residual history
2057 and the number of residuals it contains.
2058
2059 Not Collective
2060
2061 Input Parameter:
2062 . ksp - iterative context obtained from KSPCreate()
2063
2064 Output Parameters:
2065 + a - pointer to array to hold history (or NULL)
2066 - na - number of used entries in a (or NULL)
2067
2068 Level: advanced
2069
2070 Notes:
2071 Can only be called after a KSPSetResidualHistory() otherwise a and na are set to zero
2072
2073 The Fortran version of this routine has a calling sequence
2074 $ call KSPGetResidualHistory(KSP ksp, integer na, integer ierr)
2075 note that you have passed a Fortran array into KSPSetResidualHistory() and you need
2076 to access the residual values from this Fortran array you provided. Only the na (number of
2077 residual norms currently held) is set.
2078
2079 .seealso: KSPGetResidualHistory(), KSP
2080
2081 @*/
KSPGetResidualHistory(KSP ksp,PetscReal * a[],PetscInt * na)2082 PetscErrorCode KSPGetResidualHistory(KSP ksp,PetscReal *a[],PetscInt *na)
2083 {
2084 PetscFunctionBegin;
2085 PetscValidHeaderSpecific(ksp,KSP_CLASSID,1);
2086 if (a) *a = ksp->res_hist;
2087 if (na) *na = ksp->res_hist_len;
2088 PetscFunctionReturn(0);
2089 }
2090
2091 /*@C
2092 KSPSetConvergenceTest - Sets the function to be used to determine
2093 convergence.
2094
2095 Logically Collective on ksp
2096
2097 Input Parameters:
2098 + ksp - iterative context obtained from KSPCreate()
2099 . converge - pointer to the function
2100 . cctx - context for private data for the convergence routine (may be null)
2101 - destroy - a routine for destroying the context (may be null)
2102
2103 Calling sequence of converge:
2104 $ converge (KSP ksp, PetscInt it, PetscReal rnorm, KSPConvergedReason *reason,void *mctx)
2105
2106 + ksp - iterative context obtained from KSPCreate()
2107 . it - iteration number
2108 . rnorm - (estimated) 2-norm of (preconditioned) residual
2109 . reason - the reason why it has converged or diverged
2110 - cctx - optional convergence context, as set by KSPSetConvergenceTest()
2111
2112
2113 Notes:
2114 Must be called after the KSP type has been set so put this after
2115 a call to KSPSetType(), or KSPSetFromOptions().
2116
2117 The default convergence test, KSPConvergedDefault(), aborts if the
2118 residual grows to more than 10000 times the initial residual.
2119
2120 The default is a combination of relative and absolute tolerances.
2121 The residual value that is tested may be an approximation; routines
2122 that need exact values should compute them.
2123
2124 In the default PETSc convergence test, the precise values of reason
2125 are macros such as KSP_CONVERGED_RTOL, which are defined in petscksp.h.
2126
2127 Level: advanced
2128
2129 .seealso: KSPConvergedDefault(), KSPGetConvergenceContext(), KSPSetTolerances(), KSP, KSPGetConvergenceTest(), KSPGetAndClearConvergenceTest()
2130 @*/
KSPSetConvergenceTest(KSP ksp,PetscErrorCode (* converge)(KSP,PetscInt,PetscReal,KSPConvergedReason *,void *),void * cctx,PetscErrorCode (* destroy)(void *))2131 PetscErrorCode KSPSetConvergenceTest(KSP ksp,PetscErrorCode (*converge)(KSP,PetscInt,PetscReal,KSPConvergedReason*,void*),void *cctx,PetscErrorCode (*destroy)(void*))
2132 {
2133 PetscErrorCode ierr;
2134
2135 PetscFunctionBegin;
2136 PetscValidHeaderSpecific(ksp,KSP_CLASSID,1);
2137 if (ksp->convergeddestroy) {
2138 ierr = (*ksp->convergeddestroy)(ksp->cnvP);CHKERRQ(ierr);
2139 }
2140 ksp->converged = converge;
2141 ksp->convergeddestroy = destroy;
2142 ksp->cnvP = (void*)cctx;
2143 PetscFunctionReturn(0);
2144 }
2145
2146 /*@C
2147 KSPGetConvergenceTest - Gets the function to be used to determine
2148 convergence.
2149
2150 Logically Collective on ksp
2151
2152 Input Parameter:
2153 . ksp - iterative context obtained from KSPCreate()
2154
2155 Output Parameter:
2156 + converge - pointer to convergence test function
2157 . cctx - context for private data for the convergence routine (may be null)
2158 - destroy - a routine for destroying the context (may be null)
2159
2160 Calling sequence of converge:
2161 $ converge (KSP ksp, PetscInt it, PetscReal rnorm, KSPConvergedReason *reason,void *mctx)
2162
2163 + ksp - iterative context obtained from KSPCreate()
2164 . it - iteration number
2165 . rnorm - (estimated) 2-norm of (preconditioned) residual
2166 . reason - the reason why it has converged or diverged
2167 - cctx - optional convergence context, as set by KSPSetConvergenceTest()
2168
2169 Level: advanced
2170
2171 .seealso: KSPConvergedDefault(), KSPGetConvergenceContext(), KSPSetTolerances(), KSP, KSPSetConvergenceTest(), KSPGetAndClearConvergenceTest()
2172 @*/
KSPGetConvergenceTest(KSP ksp,PetscErrorCode (** converge)(KSP,PetscInt,PetscReal,KSPConvergedReason *,void *),void ** cctx,PetscErrorCode (** destroy)(void *))2173 PetscErrorCode KSPGetConvergenceTest(KSP ksp,PetscErrorCode (**converge)(KSP,PetscInt,PetscReal,KSPConvergedReason*,void*),void **cctx,PetscErrorCode (**destroy)(void*))
2174 {
2175 PetscFunctionBegin;
2176 PetscValidHeaderSpecific(ksp,KSP_CLASSID,1);
2177 if (converge) *converge = ksp->converged;
2178 if (destroy) *destroy = ksp->convergeddestroy;
2179 if (cctx) *cctx = ksp->cnvP;
2180 PetscFunctionReturn(0);
2181 }
2182
2183 /*@C
2184 KSPGetAndClearConvergenceTest - Gets the function to be used to determine convergence. Removes the current test without calling destroy on the test context
2185
2186 Logically Collective on ksp
2187
2188 Input Parameter:
2189 . ksp - iterative context obtained from KSPCreate()
2190
2191 Output Parameter:
2192 + converge - pointer to convergence test function
2193 . cctx - context for private data for the convergence routine
2194 - destroy - a routine for destroying the context
2195
2196 Calling sequence of converge:
2197 $ converge (KSP ksp, PetscInt it, PetscReal rnorm, KSPConvergedReason *reason,void *mctx)
2198
2199 + ksp - iterative context obtained from KSPCreate()
2200 . it - iteration number
2201 . rnorm - (estimated) 2-norm of (preconditioned) residual
2202 . reason - the reason why it has converged or diverged
2203 - cctx - optional convergence context, as set by KSPSetConvergenceTest()
2204
2205 Level: advanced
2206
2207 Notes: This is intended to be used to allow transferring the convergence test (and its context) to another testing object (for example another KSP) and then calling
2208 KSPSetConvergenceTest() on this original KSP. If you just called KSPGetConvergenceTest() followed by KSPSetConvergenceTest() the original context information
2209 would be destroyed and hence the transferred context would be invalid and trigger a crash on use
2210
2211 .seealso: KSPConvergedDefault(), KSPGetConvergenceContext(), KSPSetTolerances(), KSP, KSPSetConvergenceTest(), KSPGetConvergenceTest()
2212 @*/
KSPGetAndClearConvergenceTest(KSP ksp,PetscErrorCode (** converge)(KSP,PetscInt,PetscReal,KSPConvergedReason *,void *),void ** cctx,PetscErrorCode (** destroy)(void *))2213 PetscErrorCode KSPGetAndClearConvergenceTest(KSP ksp,PetscErrorCode (**converge)(KSP,PetscInt,PetscReal,KSPConvergedReason*,void*),void **cctx,PetscErrorCode (**destroy)(void*))
2214 {
2215 PetscFunctionBegin;
2216 PetscValidHeaderSpecific(ksp,KSP_CLASSID,1);
2217 *converge = ksp->converged;
2218 *destroy = ksp->convergeddestroy;
2219 *cctx = ksp->cnvP;
2220 ksp->converged = NULL;
2221 ksp->cnvP = NULL;
2222 ksp->convergeddestroy = NULL;
2223 PetscFunctionReturn(0);
2224 }
2225
2226 /*@C
2227 KSPGetConvergenceContext - Gets the convergence context set with
2228 KSPSetConvergenceTest().
2229
2230 Not Collective
2231
2232 Input Parameter:
2233 . ksp - iterative context obtained from KSPCreate()
2234
2235 Output Parameter:
2236 . ctx - monitoring context
2237
2238 Level: advanced
2239
2240 .seealso: KSPConvergedDefault(), KSPSetConvergenceTest(), KSP
2241 @*/
KSPGetConvergenceContext(KSP ksp,void ** ctx)2242 PetscErrorCode KSPGetConvergenceContext(KSP ksp,void **ctx)
2243 {
2244 PetscFunctionBegin;
2245 PetscValidHeaderSpecific(ksp,KSP_CLASSID,1);
2246 *ctx = ksp->cnvP;
2247 PetscFunctionReturn(0);
2248 }
2249
2250 /*@C
2251 KSPBuildSolution - Builds the approximate solution in a vector provided.
2252 This routine is NOT commonly needed (see KSPSolve()).
2253
2254 Collective on ksp
2255
2256 Input Parameter:
2257 . ctx - iterative context obtained from KSPCreate()
2258
2259 Output Parameter:
2260 Provide exactly one of
2261 + v - location to stash solution.
2262 - V - the solution is returned in this location. This vector is created
2263 internally. This vector should NOT be destroyed by the user with
2264 VecDestroy().
2265
2266 Notes:
2267 This routine can be used in one of two ways
2268 .vb
2269 KSPBuildSolution(ksp,NULL,&V);
2270 or
2271 KSPBuildSolution(ksp,v,NULL); or KSPBuildSolution(ksp,v,&v);
2272 .ve
2273 In the first case an internal vector is allocated to store the solution
2274 (the user cannot destroy this vector). In the second case the solution
2275 is generated in the vector that the user provides. Note that for certain
2276 methods, such as KSPCG, the second case requires a copy of the solution,
2277 while in the first case the call is essentially free since it simply
2278 returns the vector where the solution already is stored. For some methods
2279 like GMRES this is a reasonably expensive operation and should only be
2280 used in truly needed.
2281
2282 Level: advanced
2283
2284 .seealso: KSPGetSolution(), KSPBuildResidual(), KSP
2285 @*/
KSPBuildSolution(KSP ksp,Vec v,Vec * V)2286 PetscErrorCode KSPBuildSolution(KSP ksp,Vec v,Vec *V)
2287 {
2288 PetscErrorCode ierr;
2289
2290 PetscFunctionBegin;
2291 PetscValidHeaderSpecific(ksp,KSP_CLASSID,1);
2292 if (!V && !v) SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_ARG_WRONG,"Must provide either v or V");
2293 if (!V) V = &v;
2294 ierr = (*ksp->ops->buildsolution)(ksp,v,V);CHKERRQ(ierr);
2295 PetscFunctionReturn(0);
2296 }
2297
2298 /*@C
2299 KSPBuildResidual - Builds the residual in a vector provided.
2300
2301 Collective on ksp
2302
2303 Input Parameter:
2304 . ksp - iterative context obtained from KSPCreate()
2305
2306 Output Parameters:
2307 + v - optional location to stash residual. If v is not provided,
2308 then a location is generated.
2309 . t - work vector. If not provided then one is generated.
2310 - V - the residual
2311
2312 Notes:
2313 Regardless of whether or not v is provided, the residual is
2314 returned in V.
2315
2316 Level: advanced
2317
2318 .seealso: KSPBuildSolution()
2319 @*/
KSPBuildResidual(KSP ksp,Vec t,Vec v,Vec * V)2320 PetscErrorCode KSPBuildResidual(KSP ksp,Vec t,Vec v,Vec *V)
2321 {
2322 PetscErrorCode ierr;
2323 PetscBool flag = PETSC_FALSE;
2324 Vec w = v,tt = t;
2325
2326 PetscFunctionBegin;
2327 PetscValidHeaderSpecific(ksp,KSP_CLASSID,1);
2328 if (!w) {
2329 ierr = VecDuplicate(ksp->vec_rhs,&w);CHKERRQ(ierr);
2330 ierr = PetscLogObjectParent((PetscObject)ksp,(PetscObject)w);CHKERRQ(ierr);
2331 }
2332 if (!tt) {
2333 ierr = VecDuplicate(ksp->vec_sol,&tt);CHKERRQ(ierr); flag = PETSC_TRUE;
2334 ierr = PetscLogObjectParent((PetscObject)ksp,(PetscObject)tt);CHKERRQ(ierr);
2335 }
2336 ierr = (*ksp->ops->buildresidual)(ksp,tt,w,V);CHKERRQ(ierr);
2337 if (flag) {ierr = VecDestroy(&tt);CHKERRQ(ierr);}
2338 PetscFunctionReturn(0);
2339 }
2340
2341 /*@
2342 KSPSetDiagonalScale - Tells KSP to symmetrically diagonally scale the system
2343 before solving. This actually CHANGES the matrix (and right hand side).
2344
2345 Logically Collective on ksp
2346
2347 Input Parameter:
2348 + ksp - the KSP context
2349 - scale - PETSC_TRUE or PETSC_FALSE
2350
2351 Options Database Key:
2352 + -ksp_diagonal_scale -
2353 - -ksp_diagonal_scale_fix - scale the matrix back AFTER the solve
2354
2355
2356 Notes:
2357 Scales the matrix by D^(-1/2) A D^(-1/2) [D^(1/2) x ] = D^(-1/2) b
2358 where D_{ii} is 1/abs(A_{ii}) unless A_{ii} is zero and then it is 1.
2359
2360 BE CAREFUL with this routine: it actually scales the matrix and right
2361 hand side that define the system. After the system is solved the matrix
2362 and right hand side remain scaled unless you use KSPSetDiagonalScaleFix()
2363
2364 This should NOT be used within the SNES solves if you are using a line
2365 search.
2366
2367 If you use this with the PCType Eisenstat preconditioner than you can
2368 use the PCEisenstatSetNoDiagonalScaling() option, or -pc_eisenstat_no_diagonal_scaling
2369 to save some unneeded, redundant flops.
2370
2371 Level: intermediate
2372
2373 .seealso: KSPGetDiagonalScale(), KSPSetDiagonalScaleFix(), KSP
2374 @*/
KSPSetDiagonalScale(KSP ksp,PetscBool scale)2375 PetscErrorCode KSPSetDiagonalScale(KSP ksp,PetscBool scale)
2376 {
2377 PetscFunctionBegin;
2378 PetscValidHeaderSpecific(ksp,KSP_CLASSID,1);
2379 PetscValidLogicalCollectiveBool(ksp,scale,2);
2380 ksp->dscale = scale;
2381 PetscFunctionReturn(0);
2382 }
2383
2384 /*@
2385 KSPGetDiagonalScale - Checks if KSP solver scales the matrix and
2386 right hand side
2387
2388 Not Collective
2389
2390 Input Parameter:
2391 . ksp - the KSP context
2392
2393 Output Parameter:
2394 . scale - PETSC_TRUE or PETSC_FALSE
2395
2396 Notes:
2397 BE CAREFUL with this routine: it actually scales the matrix and right
2398 hand side that define the system. After the system is solved the matrix
2399 and right hand side remain scaled unless you use KSPSetDiagonalScaleFix()
2400
2401 Level: intermediate
2402
2403 .seealso: KSPSetDiagonalScale(), KSPSetDiagonalScaleFix(), KSP
2404 @*/
KSPGetDiagonalScale(KSP ksp,PetscBool * scale)2405 PetscErrorCode KSPGetDiagonalScale(KSP ksp,PetscBool *scale)
2406 {
2407 PetscFunctionBegin;
2408 PetscValidHeaderSpecific(ksp,KSP_CLASSID,1);
2409 PetscValidPointer(scale,2);
2410 *scale = ksp->dscale;
2411 PetscFunctionReturn(0);
2412 }
2413
2414 /*@
2415 KSPSetDiagonalScaleFix - Tells KSP to diagonally scale the system
2416 back after solving.
2417
2418 Logically Collective on ksp
2419
2420 Input Parameter:
2421 + ksp - the KSP context
2422 - fix - PETSC_TRUE to scale back after the system solve, PETSC_FALSE to not
2423 rescale (default)
2424
2425 Notes:
2426 Must be called after KSPSetDiagonalScale()
2427
2428 Using this will slow things down, because it rescales the matrix before and
2429 after each linear solve. This is intended mainly for testing to allow one
2430 to easily get back the original system to make sure the solution computed is
2431 accurate enough.
2432
2433 Level: intermediate
2434
2435 .seealso: KSPGetDiagonalScale(), KSPSetDiagonalScale(), KSPGetDiagonalScaleFix(), KSP
2436 @*/
KSPSetDiagonalScaleFix(KSP ksp,PetscBool fix)2437 PetscErrorCode KSPSetDiagonalScaleFix(KSP ksp,PetscBool fix)
2438 {
2439 PetscFunctionBegin;
2440 PetscValidHeaderSpecific(ksp,KSP_CLASSID,1);
2441 PetscValidLogicalCollectiveBool(ksp,fix,2);
2442 ksp->dscalefix = fix;
2443 PetscFunctionReturn(0);
2444 }
2445
2446 /*@
2447 KSPGetDiagonalScaleFix - Determines if KSP diagonally scales the system
2448 back after solving.
2449
2450 Not Collective
2451
2452 Input Parameter:
2453 . ksp - the KSP context
2454
2455 Output Parameter:
2456 . fix - PETSC_TRUE to scale back after the system solve, PETSC_FALSE to not
2457 rescale (default)
2458
2459 Notes:
2460 Must be called after KSPSetDiagonalScale()
2461
2462 If PETSC_TRUE will slow things down, because it rescales the matrix before and
2463 after each linear solve. This is intended mainly for testing to allow one
2464 to easily get back the original system to make sure the solution computed is
2465 accurate enough.
2466
2467 Level: intermediate
2468
2469 .seealso: KSPGetDiagonalScale(), KSPSetDiagonalScale(), KSPSetDiagonalScaleFix(), KSP
2470 @*/
KSPGetDiagonalScaleFix(KSP ksp,PetscBool * fix)2471 PetscErrorCode KSPGetDiagonalScaleFix(KSP ksp,PetscBool *fix)
2472 {
2473 PetscFunctionBegin;
2474 PetscValidHeaderSpecific(ksp,KSP_CLASSID,1);
2475 PetscValidPointer(fix,2);
2476 *fix = ksp->dscalefix;
2477 PetscFunctionReturn(0);
2478 }
2479
2480 /*@C
2481 KSPSetComputeOperators - set routine to compute the linear operators
2482
2483 Logically Collective
2484
2485 Input Arguments:
2486 + ksp - the KSP context
2487 . func - function to compute the operators
2488 - ctx - optional context
2489
2490 Calling sequence of func:
2491 $ func(KSP ksp,Mat A,Mat B,void *ctx)
2492
2493 + ksp - the KSP context
2494 . A - the linear operator
2495 . B - preconditioning matrix
2496 - ctx - optional user-provided context
2497
2498 Notes:
2499 The user provided func() will be called automatically at the very next call to KSPSolve(). It will not be called at future KSPSolve() calls
2500 unless either KSPSetComputeOperators() or KSPSetOperators() is called before that KSPSolve() is called.
2501
2502 To reuse the same preconditioner for the next KSPSolve() and not compute a new one based on the most recently computed matrix call KSPSetReusePreconditioner()
2503
2504 Level: beginner
2505
2506 .seealso: KSPSetOperators(), KSPSetComputeRHS(), DMKSPSetComputeOperators(), KSPSetComputeInitialGuess()
2507 @*/
KSPSetComputeOperators(KSP ksp,PetscErrorCode (* func)(KSP,Mat,Mat,void *),void * ctx)2508 PetscErrorCode KSPSetComputeOperators(KSP ksp,PetscErrorCode (*func)(KSP,Mat,Mat,void*),void *ctx)
2509 {
2510 PetscErrorCode ierr;
2511 DM dm;
2512
2513 PetscFunctionBegin;
2514 PetscValidHeaderSpecific(ksp,KSP_CLASSID,1);
2515 ierr = KSPGetDM(ksp,&dm);CHKERRQ(ierr);
2516 ierr = DMKSPSetComputeOperators(dm,func,ctx);CHKERRQ(ierr);
2517 if (ksp->setupstage == KSP_SETUP_NEWRHS) ksp->setupstage = KSP_SETUP_NEWMATRIX;
2518 PetscFunctionReturn(0);
2519 }
2520
2521 /*@C
2522 KSPSetComputeRHS - set routine to compute the right hand side of the linear system
2523
2524 Logically Collective
2525
2526 Input Arguments:
2527 + ksp - the KSP context
2528 . func - function to compute the right hand side
2529 - ctx - optional context
2530
2531 Calling sequence of func:
2532 $ func(KSP ksp,Vec b,void *ctx)
2533
2534 + ksp - the KSP context
2535 . b - right hand side of linear system
2536 - ctx - optional user-provided context
2537
2538 Notes:
2539 The routine you provide will be called EACH you call KSPSolve() to prepare the new right hand side for that solve
2540
2541 Level: beginner
2542
2543 .seealso: KSPSolve(), DMKSPSetComputeRHS(), KSPSetComputeOperators()
2544 @*/
KSPSetComputeRHS(KSP ksp,PetscErrorCode (* func)(KSP,Vec,void *),void * ctx)2545 PetscErrorCode KSPSetComputeRHS(KSP ksp,PetscErrorCode (*func)(KSP,Vec,void*),void *ctx)
2546 {
2547 PetscErrorCode ierr;
2548 DM dm;
2549
2550 PetscFunctionBegin;
2551 PetscValidHeaderSpecific(ksp,KSP_CLASSID,1);
2552 ierr = KSPGetDM(ksp,&dm);CHKERRQ(ierr);
2553 ierr = DMKSPSetComputeRHS(dm,func,ctx);CHKERRQ(ierr);
2554 PetscFunctionReturn(0);
2555 }
2556
2557 /*@C
2558 KSPSetComputeInitialGuess - set routine to compute the initial guess of the linear system
2559
2560 Logically Collective
2561
2562 Input Arguments:
2563 + ksp - the KSP context
2564 . func - function to compute the initial guess
2565 - ctx - optional context
2566
2567 Calling sequence of func:
2568 $ func(KSP ksp,Vec x,void *ctx)
2569
2570 + ksp - the KSP context
2571 . x - solution vector
2572 - ctx - optional user-provided context
2573
2574 Notes: This should only be used in conjunction with KSPSetComputeRHS(), KSPSetComputeOperators(), otherwise
2575 call KSPSetInitialGuessNonzero() and set the initial guess values in the solution vector passed to KSPSolve().
2576
2577 Level: beginner
2578
2579 .seealso: KSPSolve(), KSPSetComputeRHS(), KSPSetComputeOperators(), DMKSPSetComputeInitialGuess()
2580 @*/
KSPSetComputeInitialGuess(KSP ksp,PetscErrorCode (* func)(KSP,Vec,void *),void * ctx)2581 PetscErrorCode KSPSetComputeInitialGuess(KSP ksp,PetscErrorCode (*func)(KSP,Vec,void*),void *ctx)
2582 {
2583 PetscErrorCode ierr;
2584 DM dm;
2585
2586 PetscFunctionBegin;
2587 PetscValidHeaderSpecific(ksp,KSP_CLASSID,1);
2588 ierr = KSPGetDM(ksp,&dm);CHKERRQ(ierr);
2589 ierr = DMKSPSetComputeInitialGuess(dm,func,ctx);CHKERRQ(ierr);
2590 PetscFunctionReturn(0);
2591 }
2592