Actual source code: tcqmr.c
2: /*
3: This file contains an implementation of Tony Chan's transpose-free QMR.
5: Note: The vector dot products in the code have not been checked for the
6: complex numbers version, so most probably some are incorrect.
7: */
9: #include <../src/ksp/ksp/impls/tcqmr/tcqmrimpl.h>
11: static PetscErrorCode KSPSolve_TCQMR(KSP ksp)
12: {
13: PetscReal rnorm0,rnorm,dp1,Gamma;
14: PetscScalar theta,ep,cl1,sl1,cl,sl,sprod,tau_n1,f;
15: PetscScalar deltmp,rho,beta,eptmp,ta,s,c,tau_n,delta;
16: PetscScalar dp11,dp2,rhom1,alpha,tmp;
18: ksp->its = 0;
20: KSPInitialResidual(ksp,x,u,v,r,b);
21: VecNorm(r,NORM_2,&rnorm0); /* rnorm0 = ||r|| */
22: KSPCheckNorm(ksp,rnorm0);
23: if (ksp->normtype != KSP_NORM_NONE) ksp->rnorm = rnorm0;
24: else ksp->rnorm = 0;
25: (*ksp->converged)(ksp,0,ksp->rnorm,&ksp->reason,ksp->cnvP);
26: if (ksp->reason) return 0;
28: VecSet(um1,0.0);
29: VecCopy(r,u);
30: rnorm = rnorm0;
31: tmp = 1.0/rnorm; VecScale(u,tmp);
32: VecSet(vm1,0.0);
33: VecCopy(u,v);
34: VecCopy(u,v0);
35: VecSet(pvec1,0.0);
36: VecSet(pvec2,0.0);
37: VecSet(p,0.0);
38: theta = 0.0;
39: ep = 0.0;
40: cl1 = 0.0;
41: sl1 = 0.0;
42: cl = 0.0;
43: sl = 0.0;
44: sprod = 1.0;
45: tau_n1= rnorm0;
46: f = 1.0;
47: Gamma = 1.0;
48: rhom1 = 1.0;
50: /*
51: CALCULATE SQUARED LANCZOS vectors
52: */
53: if (ksp->normtype != KSP_NORM_NONE) ksp->rnorm = rnorm;
54: else ksp->rnorm = 0;
55: (*ksp->converged)(ksp,ksp->its,ksp->rnorm,&ksp->reason,ksp->cnvP);
56: while (!ksp->reason) {
57: KSPMonitor(ksp,ksp->its,ksp->rnorm);
58: ksp->its++;
60: KSP_PCApplyBAorAB(ksp,u,y,vtmp); /* y = A*u */
61: VecDot(y,v0,&dp11);
62: KSPCheckDot(ksp,dp11);
63: VecDot(u,v0,&dp2);
64: alpha = dp11 / dp2; /* alpha = v0'*y/v0'*u */
65: deltmp = alpha;
66: VecCopy(y,z);
67: VecAXPY(z,-alpha,u); /* z = y - alpha u */
68: VecDot(u,v0,&rho);
69: beta = rho / (f*rhom1);
70: rhom1 = rho;
71: VecCopy(z,utmp); /* up1 = (A-alpha*I)*
72: (z-2*beta*p) + f*beta*
73: beta*um1 */
74: VecAXPY(utmp,-2.0*beta,p);
75: KSP_PCApplyBAorAB(ksp,utmp,up1,vtmp);
76: VecAXPY(up1,-alpha,utmp);
77: VecAXPY(up1,f*beta*beta,um1);
78: VecNorm(up1,NORM_2,&dp1);
79: KSPCheckNorm(ksp,dp1);
80: f = 1.0 / dp1;
81: VecScale(up1,f);
82: VecAYPX(p,-beta,z); /* p = f*(z-beta*p) */
83: VecScale(p,f);
84: VecCopy(u,um1);
85: VecCopy(up1,u);
86: beta = beta/Gamma;
87: eptmp = beta;
88: KSP_PCApplyBAorAB(ksp,v,vp1,vtmp);
89: VecAXPY(vp1,-alpha,v);
90: VecAXPY(vp1,-beta,vm1);
91: VecNorm(vp1,NORM_2,&Gamma);
92: KSPCheckNorm(ksp,Gamma);
93: VecScale(vp1,1.0/Gamma);
94: VecCopy(v,vm1);
95: VecCopy(vp1,v);
97: /*
98: SOLVE Ax = b
99: */
100: /* Apply last two Given's (Gl-1 and Gl) rotations to (beta,alpha,Gamma) */
101: if (ksp->its > 2) {
102: theta = sl1*beta;
103: eptmp = -cl1*beta;
104: }
105: if (ksp->its > 1) {
106: ep = -cl*eptmp + sl*alpha;
107: deltmp = -sl*eptmp - cl*alpha;
108: }
109: if (PetscAbsReal(Gamma) > PetscAbsScalar(deltmp)) {
110: ta = -deltmp / Gamma;
111: s = 1.0 / PetscSqrtScalar(1.0 + ta*ta);
112: c = s*ta;
113: } else {
114: ta = -Gamma/deltmp;
115: c = 1.0 / PetscSqrtScalar(1.0 + ta*ta);
116: s = c*ta;
117: }
119: delta = -c*deltmp + s*Gamma;
120: tau_n = -c*tau_n1; tau_n1 = -s*tau_n1;
121: VecCopy(vm1,pvec);
122: VecAXPY(pvec,-theta,pvec2);
123: VecAXPY(pvec,-ep,pvec1);
124: VecScale(pvec,1.0/delta);
125: VecAXPY(x,tau_n,pvec);
126: cl1 = cl; sl1 = sl; cl = c; sl = s;
128: VecCopy(pvec1,pvec2);
129: VecCopy(pvec,pvec1);
131: /* Compute the upper bound on the residual norm r (See QMR paper p. 13) */
132: sprod = sprod*PetscAbsScalar(s);
133: rnorm = rnorm0 * PetscSqrtReal((PetscReal)ksp->its+2.0) * PetscRealPart(sprod);
134: if (ksp->normtype != KSP_NORM_NONE) ksp->rnorm = rnorm;
135: else ksp->rnorm = 0;
136: (*ksp->converged)(ksp,ksp->its,ksp->rnorm,&ksp->reason,ksp->cnvP);
137: if (ksp->its >= ksp->max_it) {
138: if (!ksp->reason) ksp->reason = KSP_DIVERGED_ITS;
139: break;
140: }
141: }
142: KSPMonitor(ksp,ksp->its,ksp->rnorm);
143: KSPUnwindPreconditioner(ksp,x,vtmp);
144: return 0;
145: }
147: static PetscErrorCode KSPSetUp_TCQMR(KSP ksp)
148: {
150: KSPSetWorkVecs(ksp,TCQMR_VECS);
151: return 0;
152: }
154: /*MC
155: KSPTCQMR - A variant of QMR (quasi minimal residual) developed by Tony Chan
157: Options Database Keys:
158: see KSPSolve()
160: Level: beginner
162: Notes:
163: Supports either left or right preconditioning, but not symmetric
165: The "residual norm" computed in this algorithm is actually just an upper bound on the actual residual norm.
166: That is for left preconditioning it is a bound on the preconditioned residual and for right preconditioning
167: it is a bound on the true residual.
169: References:
170: . * - Tony F. Chan, Lisette de Pillis, and Henk van der Vorst, Transpose free formulations of Lanczos type methods for nonsymmetric linear systems,
171: Numerical Algorithms, Volume 17, 1998.
173: .seealso: KSPCreate(), KSPSetType(), KSPType (for list of available types), KSP, KSPTFQMR
175: M*/
177: PETSC_EXTERN PetscErrorCode KSPCreate_TCQMR(KSP ksp)
178: {
179: KSPSetSupportedNorm(ksp,KSP_NORM_PRECONDITIONED,PC_LEFT,3);
180: KSPSetSupportedNorm(ksp,KSP_NORM_UNPRECONDITIONED,PC_RIGHT,2);
181: KSPSetSupportedNorm(ksp,KSP_NORM_NONE,PC_RIGHT,1);
183: ksp->data = (void*)0;
184: ksp->ops->buildsolution = KSPBuildSolutionDefault;
185: ksp->ops->buildresidual = KSPBuildResidualDefault;
186: ksp->ops->setup = KSPSetUp_TCQMR;
187: ksp->ops->solve = KSPSolve_TCQMR;
188: ksp->ops->destroy = KSPDestroyDefault;
189: ksp->ops->setfromoptions = NULL;
190: ksp->ops->view = NULL;
191: return 0;
192: }