Actual source code: minres.c
2: #include <petsc/private/kspimpl.h>
4: typedef struct {
5: PetscReal haptol;
6: } KSP_MINRES;
8: static PetscErrorCode KSPSetUp_MINRES(KSP ksp)
9: {
12: KSPSetWorkVecs(ksp,9);
13: return 0;
14: }
16: static PetscErrorCode KSPSolve_MINRES(KSP ksp)
17: {
18: PetscInt i;
19: PetscScalar alpha,beta,ibeta,betaold,eta,c=1.0,ceta,cold=1.0,coold,s=0.0,sold=0.0,soold;
20: PetscScalar rho0,rho1,irho1,rho2,rho3,dp = 0.0;
21: const PetscScalar none = -1.0;
22: PetscReal np;
23: Vec X,B,R,Z,U,V,W,UOLD,VOLD,WOLD,WOOLD;
24: Mat Amat,Pmat;
25: KSP_MINRES *minres = (KSP_MINRES*)ksp->data;
26: PetscBool diagonalscale;
28: PCGetDiagonalScale(ksp->pc,&diagonalscale);
31: X = ksp->vec_sol;
32: B = ksp->vec_rhs;
33: R = ksp->work[0];
34: Z = ksp->work[1];
35: U = ksp->work[2];
36: V = ksp->work[3];
37: W = ksp->work[4];
38: UOLD = ksp->work[5];
39: VOLD = ksp->work[6];
40: WOLD = ksp->work[7];
41: WOOLD = ksp->work[8];
43: PCGetOperators(ksp->pc,&Amat,&Pmat);
45: ksp->its = 0;
47: VecSet(UOLD,0.0); /* u_old <- 0 */
48: VecSet(VOLD,0.0); /* v_old <- 0 */
49: VecSet(W,0.0); /* w <- 0 */
50: VecSet(WOLD,0.0); /* w_old <- 0 */
52: if (!ksp->guess_zero) {
53: KSP_MatMult(ksp,Amat,X,R); /* r <- b - A*x */
54: VecAYPX(R,-1.0,B);
55: } else {
56: VecCopy(B,R); /* r <- b (x is 0) */
57: }
58: KSP_PCApply(ksp,R,Z); /* z <- B*r */
59: VecNorm(Z,NORM_2,&np); /* np <- ||z|| */
60: KSPCheckNorm(ksp,np);
61: VecDot(R,Z,&dp);
62: KSPCheckDot(ksp,dp);
64: if (PetscRealPart(dp) < minres->haptol && np > minres->haptol) {
66: PetscInfo(ksp,"Detected indefinite operator %g tolerance %g\n",(double)PetscRealPart(dp),(double)minres->haptol);
67: ksp->reason = KSP_DIVERGED_INDEFINITE_MAT;
68: return 0;
69: }
71: ksp->rnorm = 0.0;
72: if (ksp->normtype != KSP_NORM_NONE) ksp->rnorm = np;
73: KSPLogResidualHistory(ksp,ksp->rnorm);
74: KSPMonitor(ksp,0,ksp->rnorm);
75: (*ksp->converged)(ksp,0,ksp->rnorm,&ksp->reason,ksp->cnvP); /* test for convergence */
76: if (ksp->reason) return 0;
78: dp = PetscAbsScalar(dp);
79: dp = PetscSqrtScalar(dp);
80: beta = dp; /* beta <- sqrt(r'*z */
81: eta = beta;
83: VecCopy(R,V);
84: VecCopy(Z,U);
85: ibeta = 1.0 / beta;
86: VecScale(V,ibeta); /* v <- r / beta */
87: VecScale(U,ibeta); /* u <- z / beta */
89: i = 0;
90: do {
91: ksp->its = i+1;
93: /* Lanczos */
95: KSP_MatMult(ksp,Amat,U,R); /* r <- A*u */
96: VecDot(U,R,&alpha); /* alpha <- r'*u */
97: KSP_PCApply(ksp,R,Z); /* z <- B*r */
99: VecAXPY(R,-alpha,V); /* r <- r - alpha v */
100: VecAXPY(Z,-alpha,U); /* z <- z - alpha u */
101: VecAXPY(R,-beta,VOLD); /* r <- r - beta v_old */
102: VecAXPY(Z,-beta,UOLD); /* z <- z - beta u_old */
104: betaold = beta;
106: VecDot(R,Z,&dp);
107: KSPCheckDot(ksp,dp);
108: dp = PetscAbsScalar(dp);
109: beta = PetscSqrtScalar(dp); /* beta <- sqrt(r'*z) */
111: /* QR factorisation */
113: coold = cold; cold = c; soold = sold; sold = s;
115: rho0 = cold * alpha - coold * sold * betaold;
116: rho1 = PetscSqrtScalar(rho0*rho0 + beta*beta);
117: rho2 = sold * alpha + coold * cold * betaold;
118: rho3 = soold * betaold;
120: /* Givens rotation */
122: c = rho0 / rho1;
123: s = beta / rho1;
125: /* Update */
127: VecCopy(WOLD,WOOLD); /* w_oold <- w_old */
128: VecCopy(W,WOLD); /* w_old <- w */
130: VecCopy(U,W); /* w <- u */
131: VecAXPY(W,-rho2,WOLD); /* w <- w - rho2 w_old */
132: VecAXPY(W,-rho3,WOOLD); /* w <- w - rho3 w_oold */
133: irho1 = 1.0 / rho1;
134: VecScale(W,irho1); /* w <- w / rho1 */
136: ceta = c * eta;
137: VecAXPY(X,ceta,W); /* x <- x + c eta w */
139: /*
140: when dp is really small we have either convergence or an indefinite operator so compute true
141: residual norm to check for convergence
142: */
143: if (PetscRealPart(dp) < minres->haptol) {
144: PetscInfo(ksp,"Possible indefinite operator %g tolerance %g\n",(double)PetscRealPart(dp),(double)minres->haptol);
145: KSP_MatMult(ksp,Amat,X,VOLD);
146: VecAXPY(VOLD,none,B);
147: VecNorm(VOLD,NORM_2,&np);
148: KSPCheckNorm(ksp,np);
149: } else {
150: /* otherwise compute new residual norm via recurrence relation */
151: np *= PetscAbsScalar(s);
152: }
154: if (ksp->normtype != KSP_NORM_NONE) ksp->rnorm = np;
155: KSPLogResidualHistory(ksp,ksp->rnorm);
156: KSPMonitor(ksp,i+1,ksp->rnorm);
157: (*ksp->converged)(ksp,i+1,ksp->rnorm,&ksp->reason,ksp->cnvP); /* test for convergence */
158: if (ksp->reason) break;
160: if (PetscRealPart(dp) < minres->haptol) {
162: PetscInfo(ksp,"Detected indefinite operator %g tolerance %g\n",(double)PetscRealPart(dp),(double)minres->haptol);
163: ksp->reason = KSP_DIVERGED_INDEFINITE_MAT;
164: break;
165: }
167: eta = -s * eta;
168: VecCopy(V,VOLD);
169: VecCopy(U,UOLD);
170: VecCopy(R,V);
171: VecCopy(Z,U);
172: ibeta = 1.0 / beta;
173: VecScale(V,ibeta); /* v <- r / beta */
174: VecScale(U,ibeta); /* u <- z / beta */
176: i++;
177: } while (i<ksp->max_it);
178: if (i >= ksp->max_it) ksp->reason = KSP_DIVERGED_ITS;
179: return 0;
180: }
182: /*MC
183: KSPMINRES - This code implements the MINRES (Minimum Residual) method.
185: Options Database Keys:
186: see KSPSolve()
188: Level: beginner
190: Notes:
191: The operator and the preconditioner must be symmetric and the preconditioner must
192: be positive definite for this method.
193: Supports only left preconditioning.
195: Reference: Paige & Saunders, 1975.
197: Contributed by: Robert Scheichl: maprs@maths.bath.ac.uk
199: .seealso: KSPCreate(), KSPSetType(), KSPType (for list of available types), KSP, KSPCG, KSPCR
200: M*/
201: PETSC_EXTERN PetscErrorCode KSPCreate_MINRES(KSP ksp)
202: {
203: KSP_MINRES *minres;
205: KSPSetSupportedNorm(ksp,KSP_NORM_PRECONDITIONED,PC_LEFT,3);
206: KSPSetSupportedNorm(ksp,KSP_NORM_NONE,PC_LEFT,1);
207: PetscNewLog(ksp,&minres);
209: /* this parameter is arbitrary; but e-50 didn't work for __float128 in one example */
210: #if defined(PETSC_USE_REAL___FLOAT128)
211: minres->haptol = 1.e-100;
212: #elif defined(PETSC_USE_REAL_SINGLE)
213: minres->haptol = 1.e-25;
214: #else
215: minres->haptol = 1.e-50;
216: #endif
217: ksp->data = (void*)minres;
219: /*
220: Sets the functions that are associated with this data structure
221: (in C++ this is the same as defining virtual functions)
222: */
223: ksp->ops->setup = KSPSetUp_MINRES;
224: ksp->ops->solve = KSPSolve_MINRES;
225: ksp->ops->destroy = KSPDestroyDefault;
226: ksp->ops->setfromoptions = NULL;
227: ksp->ops->buildsolution = KSPBuildSolutionDefault;
228: ksp->ops->buildresidual = KSPBuildResidualDefault;
229: return 0;
230: }