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:   PetscCall(VecFlag(vres, ksp->reason == KSP_DIVERGED_PC_FAILED));
     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: }