Actual source code: tfqmr.c
1: #include <petsc/private/kspimpl.h>
3: static PetscErrorCode KSPSetUp_TFQMR(KSP ksp)
4: {
5: PetscFunctionBegin;
6: PetscCheck(ksp->pc_side != PC_SYMMETRIC, PetscObjectComm((PetscObject)ksp), PETSC_ERR_SUP, "no symmetric preconditioning for KSPTFQMR");
7: PetscCall(KSPSetWorkVecs(ksp, 9));
8: PetscFunctionReturn(PETSC_SUCCESS);
9: }
11: static PetscErrorCode KSPSolve_TFQMR(KSP ksp)
12: {
13: PetscInt i, m;
14: PetscScalar rho, rhoold, a, s, b, eta, etaold, psiold, cf;
15: PetscReal dp, dpold, w, dpest, tau, psi, cm;
16: Vec X, B, V, P, R, RP, T, T1, Q, U, D, AUQ;
18: PetscFunctionBegin;
19: X = ksp->vec_sol;
20: B = ksp->vec_rhs;
21: R = ksp->work[0];
22: RP = ksp->work[1];
23: V = ksp->work[2];
24: T = ksp->work[3];
25: Q = ksp->work[4];
26: P = ksp->work[5];
27: U = ksp->work[6];
28: D = ksp->work[7];
29: T1 = ksp->work[8];
30: AUQ = V;
32: /* Compute initial preconditioned residual */
33: PetscCall(KSPInitialResidual(ksp, X, V, T, R, B));
35: /* Test for nothing to do */
36: PetscCall(VecNorm(R, NORM_2, &dp));
37: KSPCheckNorm(ksp, dp);
38: PetscCall(PetscObjectSAWsTakeAccess((PetscObject)ksp));
39: if (ksp->normtype != KSP_NORM_NONE) ksp->rnorm = dp;
40: else ksp->rnorm = 0.0;
41: ksp->its = 0;
42: PetscCall(PetscObjectSAWsGrantAccess((PetscObject)ksp));
43: PetscCall(KSPLogResidualHistory(ksp, ksp->rnorm));
44: PetscCall(KSPMonitor(ksp, 0, ksp->rnorm));
45: PetscCall((*ksp->converged)(ksp, 0, ksp->rnorm, &ksp->reason, ksp->cnvP));
46: if (ksp->reason) PetscFunctionReturn(PETSC_SUCCESS);
48: /* Make the initial Rp == R */
49: PetscCall(VecCopy(R, RP));
51: /* Set the initial conditions */
52: etaold = 0.0;
53: psiold = 0.0;
54: tau = dp;
55: dpold = dp;
57: PetscCall(VecDot(R, RP, &rhoold)); /* rhoold = (r,rp) */
58: PetscCall(VecCopy(R, U));
59: PetscCall(VecCopy(R, P));
60: PetscCall(KSP_PCApplyBAorAB(ksp, P, V, T));
61: PetscCall(VecSet(D, 0.0));
63: i = 0;
64: do {
65: PetscCall(PetscObjectSAWsTakeAccess((PetscObject)ksp));
66: ksp->its++;
67: PetscCall(PetscObjectSAWsGrantAccess((PetscObject)ksp));
68: PetscCall(VecDot(V, RP, &s)); /* s <- (v,rp) */
69: KSPCheckDot(ksp, s);
70: a = rhoold / s; /* a <- rho / s */
71: PetscCall(VecWAXPY(Q, -a, V, U)); /* q <- u - a v */
72: PetscCall(VecWAXPY(T, 1.0, U, Q)); /* t <- u + q */
73: PetscCall(KSP_PCApplyBAorAB(ksp, T, AUQ, T1));
74: PetscCall(VecAXPY(R, -a, AUQ)); /* r <- r - a K (u + q) */
75: PetscCall(VecNorm(R, NORM_2, &dp));
76: KSPCheckNorm(ksp, dp);
77: for (m = 0; m < 2; m++) {
78: if (!m) w = PetscSqrtReal(dp * dpold);
79: else w = dp;
80: psi = w / tau;
81: cm = 1.0 / PetscSqrtReal(1.0 + psi * psi);
82: tau = tau * psi * cm;
83: eta = cm * cm * a;
84: cf = psiold * psiold * etaold / a;
85: if (!m) {
86: PetscCall(VecAYPX(D, cf, U));
87: } else {
88: PetscCall(VecAYPX(D, cf, Q));
89: }
90: PetscCall(VecAXPY(X, eta, D));
92: dpest = PetscSqrtReal(2 * i + m + 2.0) * tau;
93: PetscCall(PetscObjectSAWsTakeAccess((PetscObject)ksp));
94: if (ksp->normtype != KSP_NORM_NONE) ksp->rnorm = dpest;
95: else ksp->rnorm = 0.0;
96: PetscCall(PetscObjectSAWsGrantAccess((PetscObject)ksp));
97: PetscCall(KSPLogResidualHistory(ksp, ksp->rnorm));
98: PetscCall(KSPMonitor(ksp, i + 1, ksp->rnorm));
99: PetscCall((*ksp->converged)(ksp, i + 1, ksp->rnorm, &ksp->reason, ksp->cnvP));
100: if (ksp->reason) break;
102: etaold = eta;
103: psiold = psi;
104: }
105: if (ksp->reason) break;
107: PetscCall(VecDot(R, RP, &rho)); /* rho <- (r,rp) */
108: b = rho / rhoold; /* b <- rho / rhoold */
109: PetscCall(VecWAXPY(U, b, Q, R)); /* u <- r + b q */
110: PetscCall(VecAXPY(Q, b, P));
111: PetscCall(VecWAXPY(P, b, Q, U)); /* p <- u + b(q + b p) */
112: PetscCall(KSP_PCApplyBAorAB(ksp, P, V, Q)); /* v <- K p */
114: rhoold = rho;
115: dpold = dp;
117: i++;
118: } while (i < ksp->max_it);
119: if (i >= ksp->max_it) ksp->reason = KSP_DIVERGED_ITS;
121: PetscCall(KSPUnwindPreconditioner(ksp, X, T));
122: PetscFunctionReturn(PETSC_SUCCESS);
123: }
125: /*MC
126: KSPTFQMR - An implementation of transpose-free QMR (quasi minimal residual) {cite}`f:93`.
128: Level: intermediate
130: Notes:
131: Transpose-free QMR is an algorithm somewhat similar to QMR, but unlike QMR it does not need the application
132: of the transpose of the matrix or preconditioner.
134: Supports left and right preconditioning, but not symmetric
136: The "residual norm" computed in this algorithm is actually just an upper bound on the actual residual norm.
137: For left preconditioning it is a bound on the preconditioned residual norm and for right preconditioning
138: it is a bound on the true residual norm.
140: The solver has a two-step inner iteration, each of which computes and updates the solution and the residual norm.
141: Hence the values from `KSPGetResidualHistory()` and `KSPGetIterationNumber()` will differ.
143: .seealso: [](ch_ksp), `KSPCreate()`, `KSPSetType()`, `KSPType`, `KSP`, `KSPTCQMR`
144: M*/
145: PETSC_EXTERN PetscErrorCode KSPCreate_TFQMR(KSP ksp)
146: {
147: PetscFunctionBegin;
148: PetscCall(KSPSetSupportedNorm(ksp, KSP_NORM_PRECONDITIONED, PC_LEFT, 3));
149: PetscCall(KSPSetSupportedNorm(ksp, KSP_NORM_UNPRECONDITIONED, PC_RIGHT, 2));
150: PetscCall(KSPSetSupportedNorm(ksp, KSP_NORM_NONE, PC_LEFT, 1));
151: PetscCall(KSPSetSupportedNorm(ksp, KSP_NORM_NONE, PC_RIGHT, 1));
153: ksp->data = NULL;
154: ksp->ops->setup = KSPSetUp_TFQMR;
155: ksp->ops->solve = KSPSolve_TFQMR;
156: ksp->ops->destroy = KSPDestroyDefault;
157: ksp->ops->buildsolution = KSPBuildSolutionDefault;
158: ksp->ops->buildresidual = KSPBuildResidualDefault;
159: ksp->ops->setfromoptions = NULL;
160: ksp->ops->view = NULL;
161: PetscFunctionReturn(PETSC_SUCCESS);
162: }