Actual source code: cgs.c


  2: /*
  3:     Note that for the complex numbers version, the VecDot() arguments
  4:     within the code MUST remain in the order given for correct computation
  5:     of inner products.
  6: */
  7: #include <petsc/private/kspimpl.h>

  9: static PetscErrorCode KSPSetUp_CGS(KSP ksp)
 10: {
 11:   PetscFunctionBegin;
 12:   PetscCall(KSPSetWorkVecs(ksp, 7));
 13:   PetscFunctionReturn(PETSC_SUCCESS);
 14: }

 16: static PetscErrorCode KSPSolve_CGS(KSP ksp)
 17: {
 18:   PetscInt    i;
 19:   PetscScalar rho, rhoold, a, s, b;
 20:   Vec         X, B, V, P, R, RP, T, Q, U, AUQ;
 21:   PetscReal   dp = 0.0;
 22:   PetscBool   diagonalscale;

 24:   PetscFunctionBegin;
 25:   /* not sure what residual norm it does use, should use for right preconditioning */

 27:   PetscCall(PCGetDiagonalScale(ksp->pc, &diagonalscale));
 28:   PetscCheck(!diagonalscale, PetscObjectComm((PetscObject)ksp), PETSC_ERR_SUP, "Krylov method %s does not support diagonal scaling", ((PetscObject)ksp)->type_name);

 30:   X   = ksp->vec_sol;
 31:   B   = ksp->vec_rhs;
 32:   R   = ksp->work[0];
 33:   RP  = ksp->work[1];
 34:   V   = ksp->work[2];
 35:   T   = ksp->work[3];
 36:   Q   = ksp->work[4];
 37:   P   = ksp->work[5];
 38:   U   = ksp->work[6];
 39:   AUQ = V;

 41:   /* Compute initial preconditioned residual */
 42:   PetscCall(KSPInitialResidual(ksp, X, V, T, R, B));

 44:   /* Test for nothing to do */
 45:   if (ksp->normtype != KSP_NORM_NONE) {
 46:     PetscCall(VecNorm(R, NORM_2, &dp));
 47:     KSPCheckNorm(ksp, dp);
 48:     if (ksp->normtype == KSP_NORM_NATURAL) dp *= dp;
 49:   } else dp = 0.0;

 51:   PetscCall(PetscObjectSAWsTakeAccess((PetscObject)ksp));
 52:   ksp->its   = 0;
 53:   ksp->rnorm = dp;
 54:   PetscCall(PetscObjectSAWsGrantAccess((PetscObject)ksp));
 55:   PetscCall(KSPLogResidualHistory(ksp, dp));
 56:   PetscCall(KSPMonitor(ksp, 0, dp));
 57:   PetscCall((*ksp->converged)(ksp, 0, dp, &ksp->reason, ksp->cnvP));
 58:   if (ksp->reason) PetscFunctionReturn(PETSC_SUCCESS);

 60:   /* Make the initial Rp == R */
 61:   PetscCall(VecCopy(R, RP));
 62:   /*  added for Fidap */
 63:   /* Penalize Startup - Isaac Hasbani Trick for CGS
 64:      Since most initial conditions result in a mostly 0 residual,
 65:      we change all the 0 values in the vector RP to the maximum.
 66:   */
 67:   if (ksp->normtype == KSP_NORM_NATURAL) {
 68:     PetscReal    vr0max;
 69:     PetscScalar *tmp_RP = NULL;
 70:     PetscInt     numnp = 0, *max_pos = NULL;
 71:     PetscCall(VecMax(RP, max_pos, &vr0max));
 72:     PetscCall(VecGetArray(RP, &tmp_RP));
 73:     PetscCall(VecGetLocalSize(RP, &numnp));
 74:     for (i = 0; i < numnp; i++) {
 75:       if (tmp_RP[i] == 0.0) tmp_RP[i] = vr0max;
 76:     }
 77:     PetscCall(VecRestoreArray(RP, &tmp_RP));
 78:   }
 79:   /*  end of addition for Fidap */

 81:   /* Set the initial conditions */
 82:   PetscCall(VecDot(R, RP, &rhoold)); /* rhoold = (r,rp)      */
 83:   PetscCall(VecCopy(R, U));
 84:   PetscCall(VecCopy(R, P));
 85:   PetscCall(KSP_PCApplyBAorAB(ksp, P, V, T));

 87:   i = 0;
 88:   do {
 89:     PetscCall(VecDot(V, RP, &s)); /* s <- (v,rp)          */
 90:     KSPCheckDot(ksp, s);
 91:     a = rhoold / s;                    /* a <- rho / s         */
 92:     PetscCall(VecWAXPY(Q, -a, V, U));  /* q <- u - a v         */
 93:     PetscCall(VecWAXPY(T, 1.0, U, Q)); /* t <- u + q           */
 94:     PetscCall(VecAXPY(X, a, T));       /* x <- x + a (u + q)   */
 95:     PetscCall(KSP_PCApplyBAorAB(ksp, T, AUQ, U));
 96:     PetscCall(VecAXPY(R, -a, AUQ)); /* r <- r - a K (u + q) */
 97:     PetscCall(VecDot(R, RP, &rho)); /* rho <- (r,rp)        */
 98:     KSPCheckDot(ksp, rho);
 99:     if (ksp->normtype == KSP_NORM_NATURAL) {
100:       dp = PetscAbsScalar(rho);
101:     } else if (ksp->normtype != KSP_NORM_NONE) {
102:       PetscCall(VecNorm(R, NORM_2, &dp));
103:       KSPCheckNorm(ksp, dp);
104:     } else dp = 0.0;

106:     PetscCall(PetscObjectSAWsTakeAccess((PetscObject)ksp));
107:     ksp->its++;
108:     ksp->rnorm = dp;
109:     PetscCall(PetscObjectSAWsGrantAccess((PetscObject)ksp));
110:     PetscCall(KSPLogResidualHistory(ksp, dp));
111:     PetscCall(KSPMonitor(ksp, i + 1, dp));
112:     PetscCall((*ksp->converged)(ksp, i + 1, dp, &ksp->reason, ksp->cnvP));
113:     if (ksp->reason) break;

115:     b = rho / rhoold;                /* b <- rho / rhoold    */
116:     PetscCall(VecWAXPY(U, b, Q, R)); /* u <- r + b q         */
117:     PetscCall(VecAXPY(Q, b, P));
118:     PetscCall(VecWAXPY(P, b, Q, U));            /* p <- u + b(q + b p)  */
119:     PetscCall(KSP_PCApplyBAorAB(ksp, P, V, Q)); /* v <- K p    */
120:     rhoold = rho;
121:     i++;
122:   } while (i < ksp->max_it);
123:   if (i >= ksp->max_it) ksp->reason = KSP_DIVERGED_ITS;

125:   PetscCall(KSPUnwindPreconditioner(ksp, X, T));
126:   PetscFunctionReturn(PETSC_SUCCESS);
127: }

129: /*MC
130:      KSPCGS - This code implements the CGS (Conjugate Gradient Squared) method.

132:    Level: beginner

134:    Notes:
135:    Does not require a symmetric matrix. Does not apply transpose of the matrix.

137:    Supports left and right preconditioning, but not symmetric.

139:    Developer Note:
140:    Has this weird support for doing the convergence test with the natural norm, I assume this works only with
141:    no preconditioning and symmetric positive definite operator.

143:    References:
144: .  * - Sonneveld, 1989.

146: .seealso: [](ch_ksp), `KSPCreate()`, `KSPSetType()`, `KSPType`, `KSP`, `KSPBCGS`, `KSPSetPCSide()`
147: M*/
148: PETSC_EXTERN PetscErrorCode KSPCreate_CGS(KSP ksp)
149: {
150:   PetscFunctionBegin;
151:   ksp->data = (void *)0;

153:   PetscCall(KSPSetSupportedNorm(ksp, KSP_NORM_PRECONDITIONED, PC_LEFT, 3));
154:   PetscCall(KSPSetSupportedNorm(ksp, KSP_NORM_UNPRECONDITIONED, PC_RIGHT, 2));
155:   PetscCall(KSPSetSupportedNorm(ksp, KSP_NORM_NATURAL, PC_LEFT, 2));
156:   PetscCall(KSPSetSupportedNorm(ksp, KSP_NORM_NATURAL, PC_RIGHT, 2));
157:   PetscCall(KSPSetSupportedNorm(ksp, KSP_NORM_NONE, PC_LEFT, 1));
158:   PetscCall(KSPSetSupportedNorm(ksp, KSP_NORM_NONE, PC_RIGHT, 1));

160:   ksp->ops->setup          = KSPSetUp_CGS;
161:   ksp->ops->solve          = KSPSolve_CGS;
162:   ksp->ops->destroy        = KSPDestroyDefault;
163:   ksp->ops->buildsolution  = KSPBuildSolutionDefault;
164:   ksp->ops->buildresidual  = KSPBuildResidualDefault;
165:   ksp->ops->setfromoptions = NULL;
166:   ksp->ops->view           = NULL;
167:   PetscFunctionReturn(PETSC_SUCCESS);
168: }