Actual source code: tcqmr.c
petsc-3.5.4 2015-05-23
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 <petsc-private/kspimpl.h>
10: #include <../src/ksp/ksp/impls/tcqmr/tcqmrimpl.h>
14: static PetscErrorCode KSPSolve_TCQMR(KSP ksp)
15: {
16: PetscReal rnorm0,rnorm,dp1,Gamma;
17: PetscScalar theta,ep,cl1,sl1,cl,sl,sprod,tau_n1,f;
18: PetscScalar deltmp,rho,beta,eptmp,ta,s,c,tau_n,delta;
19: PetscScalar dp11,dp2,rhom1,alpha,tmp;
23: ksp->its = 0;
25: KSPInitialResidual(ksp,x,u,v,r,b);
26: VecNorm(r,NORM_2,&rnorm0); /* rnorm0 = ||r|| */
28: (*ksp->converged)(ksp,0,rnorm0,&ksp->reason,ksp->cnvP);
29: if (ksp->reason) return(0);
31: VecSet(um1,0.0);
32: VecCopy(r,u);
33: rnorm = rnorm0;
34: tmp = 1.0/rnorm; VecScale(u,tmp);
35: VecSet(vm1,0.0);
36: VecCopy(u,v);
37: VecCopy(u,v0);
38: VecSet(pvec1,0.0);
39: VecSet(pvec2,0.0);
40: VecSet(p,0.0);
41: theta = 0.0;
42: ep = 0.0;
43: cl1 = 0.0;
44: sl1 = 0.0;
45: cl = 0.0;
46: sl = 0.0;
47: sprod = 1.0;
48: tau_n1= rnorm0;
49: f = 1.0;
50: Gamma = 1.0;
51: rhom1 = 1.0;
53: /*
54: CALCULATE SQUARED LANCZOS vectors
55: */
56: (*ksp->converged)(ksp,ksp->its,rnorm,&ksp->reason,ksp->cnvP);
57: while (!ksp->reason) {
58: KSPMonitor(ksp,ksp->its,rnorm);
59: ksp->its++;
61: KSP_PCApplyBAorAB(ksp,u,y,vtmp); /* y = A*u */
62: VecDot(y,v0,&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: f = 1.0 / dp1;
80: VecScale(up1,f);
81: VecAYPX(p,-beta,z); /* p = f*(z-beta*p) */
82: VecScale(p,f);
83: VecCopy(u,um1);
84: VecCopy(up1,u);
85: beta = beta/Gamma;
86: eptmp = beta;
87: KSP_PCApplyBAorAB(ksp,v,vp1,vtmp);
88: VecAXPY(vp1,-alpha,v);
89: VecAXPY(vp1,-beta,vm1);
90: VecNorm(vp1,NORM_2,&Gamma);
91: VecScale(vp1,1.0/Gamma);
92: VecCopy(v,vm1);
93: VecCopy(vp1,v);
95: /*
96: SOLVE Ax = b
97: */
98: /* Apply last two Given's (Gl-1 and Gl) rotations to (beta,alpha,Gamma) */
99: if (ksp->its > 2) {
100: theta = sl1*beta;
101: eptmp = -cl1*beta;
102: }
103: if (ksp->its > 1) {
104: ep = -cl*eptmp + sl*alpha;
105: deltmp = -sl*eptmp - cl*alpha;
106: }
107: if (PetscAbsReal(Gamma) > PetscAbsScalar(deltmp)) {
108: ta = -deltmp / Gamma;
109: s = 1.0 / PetscSqrtScalar(1.0 + ta*ta);
110: c = s*ta;
111: } else {
112: ta = -Gamma/deltmp;
113: c = 1.0 / PetscSqrtScalar(1.0 + ta*ta);
114: s = c*ta;
115: }
117: delta = -c*deltmp + s*Gamma;
118: tau_n = -c*tau_n1; tau_n1 = -s*tau_n1;
119: VecCopy(vm1,pvec);
120: VecAXPY(pvec,-theta,pvec2);
121: VecAXPY(pvec,-ep,pvec1);
122: VecScale(pvec,1.0/delta);
123: VecAXPY(x,tau_n,pvec);
124: cl1 = cl; sl1 = sl; cl = c; sl = s;
126: VecCopy(pvec1,pvec2);
127: VecCopy(pvec,pvec1);
129: /* Compute the upper bound on the residual norm r (See QMR paper p. 13) */
130: sprod = sprod*PetscAbsScalar(s);
131: rnorm = rnorm0 * PetscSqrtReal((PetscReal)ksp->its+2.0) * PetscRealPart(sprod);
132: (*ksp->converged)(ksp,ksp->its,rnorm,&ksp->reason,ksp->cnvP);
133: if (ksp->its >= ksp->max_it) {
134: if (!ksp->reason) ksp->reason = KSP_DIVERGED_ITS;
135: break;
136: }
137: }
138: KSPMonitor(ksp,ksp->its,rnorm);
139: KSPUnwindPreconditioner(ksp,x,vtmp);
140: return(0);
141: }
145: static PetscErrorCode KSPSetUp_TCQMR(KSP ksp)
146: {
150: if (ksp->pc_side == PC_SYMMETRIC) SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_SUP,"no symmetric preconditioning for KSPTCQMR");
151: KSPSetWorkVecs(ksp,TCQMR_VECS);
152: return(0);
153: }
155: /*MC
156: KSPTCQMR - A variant of QMR (quasi minimal residual) developed by Tony Chan
158: Options Database Keys:
159: . see KSPSolve()
161: Level: beginner
163: Notes: 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: Transpose-free formulations of Lanczos-type methods for nonsymmetric linear systems,
171: Tony F. Chan, Lisette de Pillis, and Henk van der Vorst, Numerical Algorithms,
172: Volume 17, Numbers 1-2 / May, 1998 pp. 51-66.
174: .seealso: KSPCreate(), KSPSetType(), KSPType (for list of available types), KSP, KSPTFQMR
176: M*/
180: PETSC_EXTERN PetscErrorCode KSPCreate_TCQMR(KSP ksp)
181: {
185: KSPSetSupportedNorm(ksp,KSP_NORM_PRECONDITIONED,PC_LEFT,3);
186: KSPSetSupportedNorm(ksp,KSP_NORM_UNPRECONDITIONED,PC_RIGHT,2);
188: ksp->data = (void*)0;
189: ksp->ops->buildsolution = KSPBuildSolutionDefault;
190: ksp->ops->buildresidual = KSPBuildResidualDefault;
191: ksp->ops->setup = KSPSetUp_TCQMR;
192: ksp->ops->solve = KSPSolve_TCQMR;
193: ksp->ops->destroy = KSPDestroyDefault;
194: ksp->ops->setfromoptions = 0;
195: ksp->ops->view = 0;
196: return(0);
197: }