1 #include <petsc/private/kspimpl.h>
2
KSPSetUp_TFQMR(KSP ksp)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 }
10
KSPSolve_TFQMR(KSP ksp)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;
17
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;
31
32 /* Compute initial preconditioned residual */
33 PetscCall(KSPInitialResidual(ksp, X, V, T, R, B));
34
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);
47
48 /* Make the initial Rp == R */
49 PetscCall(VecCopy(R, RP));
50
51 /* Set the initial conditions */
52 etaold = 0.0;
53 psiold = 0.0;
54 tau = dp;
55 dpold = dp;
56
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));
62
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));
91
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;
101
102 etaold = eta;
103 psiold = psi;
104 }
105 if (ksp->reason) break;
106
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 */
113
114 rhoold = rho;
115 dpold = dp;
116
117 i++;
118 } while (i < ksp->max_it);
119 if (i >= ksp->max_it) ksp->reason = KSP_DIVERGED_ITS;
120
121 PetscCall(KSPUnwindPreconditioner(ksp, X, T));
122 PetscFunctionReturn(PETSC_SUCCESS);
123 }
124
125 /*MC
126 KSPTFQMR - An implementation of transpose-free QMR (quasi minimal residual) {cite}`f:93`.
127
128 Level: intermediate
129
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.
133
134 Supports left and right preconditioning, but not symmetric
135
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.
139
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.
142
143 .seealso: [](ch_ksp), `KSPCreate()`, `KSPSetType()`, `KSPType`, `KSP`, `KSPTCQMR`
144 M*/
KSPCreate_TFQMR(KSP ksp)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));
152
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 }
163