Actual source code: itres.c
1: #include <petsc/private/kspimpl.h>
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".
8: Collective
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
17: Output Parameter:
18: . vres - calculated residual
20: Level: developer
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.
33: .seealso: [](ch_ksp), `KSP`, `KSPSolve()`, `KSPMonitor()`
34: @*/
35: PetscErrorCode KSPInitialResidual(KSP ksp, Vec vsoln, Vec vt1, Vec vt2, Vec vres, Vec vb)
36: {
37: Mat Amat, Pmat;
39: PetscFunctionBegin;
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: if (ksp->reason == KSP_DIVERGED_PC_FAILED) PetscCall(VecSetInf(vres));
72: PetscFunctionReturn(PETSC_SUCCESS);
73: }
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.
80: Collective
82: Input Parameters:
83: + ksp - iterative context
84: . vsoln - solution vector
85: - vt1 - temporary work vector
87: Output Parameter:
88: . vsoln - contains solution on output
90: Level: advanced
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.
97: .seealso: [](ch_ksp), `KSP`, `KSPSetPCSide()`
98: @*/
99: PetscErrorCode KSPUnwindPreconditioner(KSP ksp, Vec vsoln, Vec vt1)
100: {
101: PetscFunctionBegin;
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: }