1 #include <petsc/private/kspimpl.h> /*I "petscksp.h" I*/
2
3 /*@
4 KSPInitialResidual - Computes the residual. Either b - A*C*u = b - A*x with right
5 preconditioning or C*(b - A*x) with left preconditioning; the latter
6 residual is often called the "preconditioned residual".
7
8 Collective
9
10 Input Parameters:
11 + ksp - the `KSP` solver object
12 . vsoln - solution to use in computing residual
13 . vt1 - temporary work vector
14 . vt2 - temporary work vector
15 - vb - right-hand-side vector
16
17 Output Parameter:
18 . vres - calculated residual
19
20 Level: developer
21
22 Note:
23 This routine assumes that an iterative method, designed for $ A x = b $
24 will be used with a preconditioner, C, such that the actual problem is either
25 .vb
26 AC u = b (right preconditioning) or
27 CA x = Cb (left preconditioning).
28 .ve
29 This means that the calculated residual will be scaled and/or preconditioned;
30 the true residual $ b-Ax $
31 is returned in the `vt2` temporary work vector.
32
33 .seealso: [](ch_ksp), `KSP`, `KSPSolve()`, `KSPMonitor()`
34 @*/
KSPInitialResidual(KSP ksp,Vec vsoln,Vec vt1,Vec vt2,Vec vres,Vec vb)35 PetscErrorCode KSPInitialResidual(KSP ksp, Vec vsoln, Vec vt1, Vec vt2, Vec vres, Vec vb)
36 {
37 Mat Amat, Pmat;
38
39 PetscFunctionBegin;
40 PetscValidHeaderSpecific(ksp, KSP_CLASSID, 1);
41 PetscValidHeaderSpecific(vsoln, VEC_CLASSID, 2);
42 PetscValidHeaderSpecific(vres, VEC_CLASSID, 5);
43 PetscValidHeaderSpecific(vb, VEC_CLASSID, 6);
44 if (!ksp->pc) PetscCall(KSPGetPC(ksp, &ksp->pc));
45 PetscCall(PCGetOperators(ksp->pc, &Amat, &Pmat));
46 if (!ksp->guess_zero) {
47 /* skip right scaling since current guess already has it */
48 PetscCall(KSP_MatMult(ksp, Amat, vsoln, vt1));
49 PetscCall(VecCopy(vb, vt2));
50 PetscCall(VecAXPY(vt2, -1.0, vt1));
51 if (ksp->pc_side == PC_RIGHT) {
52 PetscCall(PCDiagonalScaleLeft(ksp->pc, vt2, vres));
53 } else if (ksp->pc_side == PC_LEFT) {
54 PetscCall(KSP_PCApply(ksp, vt2, vres));
55 PetscCall(PCDiagonalScaleLeft(ksp->pc, vres, vres));
56 } else if (ksp->pc_side == PC_SYMMETRIC) {
57 PetscCall(PCApplySymmetricLeft(ksp->pc, vt2, vres));
58 } else SETERRQ(PetscObjectComm((PetscObject)ksp), PETSC_ERR_SUP, "Invalid preconditioning side %d", (int)ksp->pc_side);
59 } else {
60 PetscCall(VecCopy(vb, vt2));
61 if (ksp->pc_side == PC_RIGHT) {
62 PetscCall(PCDiagonalScaleLeft(ksp->pc, vb, vres));
63 } else if (ksp->pc_side == PC_LEFT) {
64 PetscCall(KSP_PCApply(ksp, vb, vres));
65 PetscCall(PCDiagonalScaleLeft(ksp->pc, vres, vres));
66 } else if (ksp->pc_side == PC_SYMMETRIC) {
67 PetscCall(PCApplySymmetricLeft(ksp->pc, vb, vres));
68 } else SETERRQ(PetscObjectComm((PetscObject)ksp), PETSC_ERR_SUP, "Invalid preconditioning side %d", (int)ksp->pc_side);
69 }
70 /* This may be true only on a subset of MPI ranks; setting it here so it will be detected by the first norm computation in the Krylov method */
71 PetscCall(VecFlag(vres, ksp->reason == KSP_DIVERGED_PC_FAILED));
72 PetscFunctionReturn(PETSC_SUCCESS);
73 }
74
75 /*@
76 KSPUnwindPreconditioner - Unwinds the preconditioning in the solution. That is,
77 takes solution to the preconditioned problem and gets the solution to the
78 original problem from it.
79
80 Collective
81
82 Input Parameters:
83 + ksp - iterative context
84 . vsoln - solution vector
85 - vt1 - temporary work vector
86
87 Output Parameter:
88 . vsoln - contains solution on output
89
90 Level: advanced
91
92 Note:
93 If preconditioning either symmetrically or on the right, this routine solves
94 for the correction to the unpreconditioned problem. If preconditioning on
95 the left, nothing is done.
96
97 .seealso: [](ch_ksp), `KSP`, `KSPSetPCSide()`
98 @*/
KSPUnwindPreconditioner(KSP ksp,Vec vsoln,Vec vt1)99 PetscErrorCode KSPUnwindPreconditioner(KSP ksp, Vec vsoln, Vec vt1)
100 {
101 PetscFunctionBegin;
102 PetscValidHeaderSpecific(ksp, KSP_CLASSID, 1);
103 PetscValidHeaderSpecific(vsoln, VEC_CLASSID, 2);
104 if (!ksp->pc) PetscCall(KSPGetPC(ksp, &ksp->pc));
105 if (ksp->pc_side == PC_RIGHT) {
106 PetscCall(KSP_PCApply(ksp, vsoln, vt1));
107 PetscCall(PCDiagonalScaleRight(ksp->pc, vt1, vsoln));
108 } else if (ksp->pc_side == PC_SYMMETRIC) {
109 PetscCall(PCApplySymmetricRight(ksp->pc, vsoln, vt1));
110 PetscCall(VecCopy(vt1, vsoln));
111 } else {
112 PetscCall(PCDiagonalScaleRight(ksp->pc, vsoln, vsoln));
113 }
114 PetscFunctionReturn(PETSC_SUCCESS);
115 }
116