Actual source code: minres.c
1: /*$Id: minres.c,v 1.18 2001/08/10 18:14:38 bsmith Exp $*/
2: /*
3: This code implements the MINRES (Minimum Residual) method.
4: Reference: Paige & Saunders, 1975.
6: Contributed by: Robert Scheichl: maprs@maths.bath.ac.uk
8: */
9: #include src/sles/ksp/kspimpl.h
11: typedef struct {
12: PetscReal haptol;
13: } KSP_MINRES;
15: #undef __FUNCT__
17: int KSPSetUp_MINRES(KSP ksp)
18: {
23: if (ksp->pc_side == PC_RIGHT) {
24: SETERRQ(2,"No right preconditioning for KSPMINRES");
25: } else if (ksp->pc_side == PC_SYMMETRIC) {
26: SETERRQ(2,"No symmetric preconditioning for KSPMINRES");
27: }
29: KSPDefaultGetWork(ksp,9);
31: return(0);
32: }
35: #undef __FUNCT__
37: int KSPSolve_MINRES(KSP ksp,int *its)
38: {
39: int ierr,i,maxit;
40: PetscScalar alpha,malpha,beta,mbeta,ibeta,betaold,eta,c=1.0,ceta,cold=1.0,coold,s=0.0,sold=0.0,soold;
41: PetscScalar rho0,rho1,irho1,rho2,mrho2,rho3,mrho3,mone = -1.0,zero = 0.0,dp = 0.0;
42: PetscReal np;
43: Vec X,B,R,Z,U,V,W,UOLD,VOLD,WOLD,WOOLD;
44: Mat Amat,Pmat;
45: MatStructure pflag;
46: KSP_MINRES *minres = (KSP_MINRES*)ksp->data;
47: PetscTruth diagonalscale;
50: ierr = PCDiagonalScale(ksp->B,&diagonalscale);
51: if (diagonalscale) SETERRQ1(1,"Krylov method %s does not support diagonal scaling",ksp->type_name);
53: maxit = ksp->max_it;
54: X = ksp->vec_sol;
55: B = ksp->vec_rhs;
56: R = ksp->work[0];
57: Z = ksp->work[1];
58: U = ksp->work[2];
59: V = ksp->work[3];
60: W = ksp->work[4];
61: UOLD = ksp->work[5];
62: VOLD = ksp->work[6];
63: WOLD = ksp->work[7];
64: WOOLD = ksp->work[8];
66: PCGetOperators(ksp->B,&Amat,&Pmat,&pflag);
68: ksp->its = 0;
70: VecSet(&zero,UOLD); /* u_old <- 0 */
71: VecCopy(UOLD,VOLD); /* v_old <- 0 */
72: VecCopy(UOLD,W); /* w <- 0 */
73: VecCopy(UOLD,WOLD); /* w_old <- 0 */
75: if (!ksp->guess_zero) {
76: KSP_MatMult(ksp,Amat,X,R); /* r <- b - A*x */
77: VecAYPX(&mone,B,R);
78: } else {
79: VecCopy(B,R); /* r <- b (x is 0) */
80: }
82: KSP_PCApply(ksp,ksp->B,R,Z); /* z <- B*r */
84: VecDot(R,Z,&dp);
85: if (PetscAbsScalar(dp) < minres->haptol) {
86: PetscLogInfo(ksp,"KSPSolve_MINRES:Detected happy breakdown %g tolerance %gn",PetscAbsScalar(dp),minres->haptol);
87: dp = PetscAbsScalar(dp); /* tiny number, can't use 0.0, cause divided by below */
88: if (dp == 0.0) {
89: ksp->reason = KSP_CONVERGED_ATOL;
90: *its = 0;
91: return(0);
92: }
93: }
95: #if !defined(PETSC_USE_COMPLEX)
96: if (dp < 0.0) SETERRQ(PETSC_ERR_KSP_BRKDWN,"Indefinite preconditioner");
97: #endif
98: dp = PetscSqrtScalar(dp);
99: beta = dp; /* beta <- sqrt(r'*z */
100: eta = beta;
102: VecCopy(R,V);
103: VecCopy(Z,U);
104: ibeta = 1.0 / beta;
105: VecScale(&ibeta,V); /* v <- r / beta */
106: VecScale(&ibeta,U); /* u <- z / beta */
108: VecNorm(Z,NORM_2,&np); /* np <- ||z|| */
110: (*ksp->converged)(ksp,0,np,&ksp->reason,ksp->cnvP); /* test for convergence */
111: if (ksp->reason) {*its = 0; return(0);}
112: KSPLogResidualHistory(ksp,np);
113: KSPMonitor(ksp,0,np); /* call any registered monitor routines */
114: ksp->rnorm = np;
116: for (i=0; i<maxit; i++) {
117: ksp->its = i+1;
119: /* Lanczos */
121: KSP_MatMult(ksp,Amat,U,R); /* r <- A*u */
122: VecDot(U,R,&alpha); /* alpha <- r'*u */
123: KSP_PCApply(ksp,ksp->B,R,Z); /* z <- B*r */
125: malpha = - alpha;
126: VecAXPY(&malpha,V,R); /* r <- r - alpha v */
127: VecAXPY(&malpha,U,Z); /* z <- z - alpha u */
128: mbeta = - beta;
129: VecAXPY(&mbeta,VOLD,R); /* r <- r - beta v_old */
130: VecAXPY(&mbeta,UOLD,Z); /* z <- z - beta u_old */
132: betaold = beta;
134: VecDot(R,Z,&dp);
135: if (PetscAbsScalar(dp) < minres->haptol) {
136: PetscLogInfo(ksp,"KSPSolve_MINRES:Detected happy breakdown %g tolerance %gn",PetscAbsScalar(dp),minres->haptol);
137: dp = PetscAbsScalar(dp); /* tiny number, can we use 0.0? */
138: }
140: #if !defined(PETSC_USE_COMPLEX)
141: if (dp < 0.0) SETERRQ1(PETSC_ERR_KSP_BRKDWN,"Indefinite preconditioner R'Z = %g",dp);
142: #endif
143: beta = PetscSqrtScalar(dp); /* beta <- sqrt(r'*z) */
145: /* QR factorisation */
147: coold = cold; cold = c; soold = sold; sold = s;
149: rho0 = cold * alpha - coold * sold * betaold;
150: rho1 = PetscSqrtScalar(rho0*rho0 + beta*beta);
151: rho2 = sold * alpha + coold * cold * betaold;
152: rho3 = soold * betaold;
154: /* Givens rotation */
156: c = rho0 / rho1;
157: s = beta / rho1;
159: /* Update */
161: VecCopy(WOLD,WOOLD); /* w_oold <- w_old */
162: VecCopy(W,WOLD); /* w_old <- w */
163:
164: VecCopy(U,W); /* w <- u */
165: mrho2 = - rho2;
166: VecAXPY(&mrho2,WOLD,W); /* w <- w - rho2 w_old */
167: mrho3 = - rho3;
168: VecAXPY(&mrho3,WOOLD,W); /* w <- w - rho3 w_oold */
169: irho1 = 1.0 / rho1;
170: VecScale(&irho1,W); /* w <- w / rho1 */
172: ceta = c * eta;
173: VecAXPY(&ceta,W,X); /* x <- x + c eta w */
174: eta = - s * eta;
176: VecCopy(V,VOLD);
177: VecCopy(U,UOLD);
178: VecCopy(R,V);
179: VecCopy(Z,U);
180: ibeta = 1.0 / beta;
181: VecScale(&ibeta,V); /* v <- r / beta */
182: VecScale(&ibeta,U); /* u <- z / beta */
183:
184: np = ksp->rnorm * PetscAbsScalar(s);
186: ksp->rnorm = np;
187: KSPLogResidualHistory(ksp,np);
188: KSPMonitor(ksp,i+1,np);
189: (*ksp->converged)(ksp,i+1,np,&ksp->reason,ksp->cnvP); /* test for convergence */
190: if (ksp->reason) break;
191: }
192: if (i == maxit) {
193: ksp->reason = KSP_DIVERGED_ITS;
194: }
195: *its = ksp->its;
196: return(0);
197: }
199: EXTERN_C_BEGIN
200: #undef __FUNCT__
202: int KSPCreate_MINRES(KSP ksp)
203: {
204: KSP_MINRES *minres;
209: ksp->pc_side = PC_LEFT;
210: ierr = PetscNew(KSP_MINRES,&minres);
211: minres->haptol = 1.e-18;
212: ksp->data = (void*)minres;
214: /*
215: Sets the functions that are associated with this data structure
216: (in C++ this is the same as defining virtual functions)
217: */
218: ksp->ops->setup = KSPSetUp_MINRES;
219: ksp->ops->solve = KSPSolve_MINRES;
220: ksp->ops->destroy = KSPDefaultDestroy;
221: ksp->ops->setfromoptions = 0;
222: ksp->ops->buildsolution = KSPDefaultBuildSolution;
223: ksp->ops->buildresidual = KSPDefaultBuildResidual;
225: return(0);
226: }
227: EXTERN_C_END