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 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;

 39:   PetscFunctionBegin;
 46:   PetscCall(PCGetOperators(ksp->pc, &Amat, NULL));
 47:   if (!ksp->guess_zero) {
 48:     PetscCall(KSP_MatMult(ksp, Amat, vsoln, vt1));
 49:     PetscCall(VecWAXPY(vt2, -1.0, vt1, vb));
 50:   } else PetscCall(VecCopy(vb, vt2));
 51:   if (ksp->pc_side == PC_RIGHT) {
 52:     PetscCall(VecCopy(vt2, vres));
 53:   } else if (ksp->pc_side == PC_LEFT) {
 54:     PetscCall(KSP_PCApply(ksp, vt2, vres));
 55:   } else if (ksp->pc_side == PC_SYMMETRIC) {
 56:     PetscCall(PCApplySymmetricLeft(ksp->pc, vt2, vres));
 57:   } else SETERRQ(PetscObjectComm((PetscObject)ksp), PETSC_ERR_SUP, "Invalid preconditioning side %d", (int)ksp->pc_side);
 58:   /* 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 */
 59:   PetscCall(VecFlag(vres, ksp->reason == KSP_DIVERGED_PC_FAILED));
 60:   PetscFunctionReturn(PETSC_SUCCESS);
 61: }

 63: /*@
 64:   KSPUnwindPreconditioner - Unwinds the preconditioning in the solution. That is,
 65:   takes solution to the preconditioned problem and gets the solution to the
 66:   original problem from it.

 68:   Collective

 70:   Input Parameters:
 71: + ksp   - iterative context
 72: . vsoln - solution vector
 73: - vt1   - temporary work vector

 75:   Output Parameter:
 76: . vsoln - contains solution on output

 78:   Level: advanced

 80:   Note:
 81:   If preconditioning either symmetrically or on the right, this routine solves
 82:   for the correction to the unpreconditioned problem.  If preconditioning on
 83:   the left, nothing is done.

 85: .seealso: [](ch_ksp), `KSP`, `KSPSetPCSide()`
 86: @*/
 87: PetscErrorCode KSPUnwindPreconditioner(KSP ksp, Vec vsoln, Vec vt1)
 88: {
 89:   PetscFunctionBegin;
 93:   if (ksp->pc_side == PC_RIGHT) {
 94:     PetscCall(KSP_PCApply(ksp, vsoln, vt1));
 95:     PetscCall(VecCopy(vt1, vsoln));
 96:   } else if (ksp->pc_side == PC_SYMMETRIC) {
 97:     PetscCall(PCApplySymmetricRight(ksp->pc, vsoln, vt1));
 98:     PetscCall(VecCopy(vt1, vsoln));
 99:   }
100:   PetscFunctionReturn(PETSC_SUCCESS);
101: }