Actual source code: itres.c


  2: #include <petsc/private/kspimpl.h>

  4: /*@
  5:    KSPInitialResidual - Computes the residual. Either b - A*C*u = b - A*x with right
  6:      preconditioning or C*(b - A*x) with left preconditioning; the latter
  7:      residual is often called the "preconditioned residual".

  9:    Collective

 11:    Input Parameters:
 12: +  vsoln    - solution to use in computing residual
 13: .  vt1, vt2 - temporary work vectors
 14: -  vb       - right-hand-side vector

 16:    Output Parameter:
 17: .  vres     - calculated residual

 19:    Level: developer

 21:    Note:
 22:    This routine assumes that an iterative method, designed for
 23: $     A x = b
 24:    will be used with a preconditioner, C, such that the actual problem is either
 25: $     AC u = b (right preconditioning) or
 26: $     CA x = Cb (left preconditioning).
 27:    This means that the calculated residual will be scaled and/or preconditioned;
 28:    the true residual
 29: $     b-Ax
 30:    is returned in the `vt2` temporary.

 32: .seealso: [](ch_ksp), `KSP`, `KSPSolve()`, `KSPMonitor()`
 33: @*/

 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:    Note:
 91:    If preconditioning either symmetrically or on the right, this routine solves
 92:    for the correction to the unpreconditioned problem.  If preconditioning on
 93:    the left, nothing is done.

 95:    Level: advanced

 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: }