Actual source code: minres.c
petsc-3.6.4 2016-04-12
2: #include <petsc/private/kspimpl.h>
4: typedef struct {
5: PetscReal haptol;
6: } KSP_MINRES;
10: PetscErrorCode KSPSetUp_MINRES(KSP ksp)
11: {
15: if (ksp->pc_side == PC_RIGHT) SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_SUP,"No right preconditioning for KSPMINRES");
16: else if (ksp->pc_side == PC_SYMMETRIC) SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_SUP,"No symmetric preconditioning for KSPMINRES");
17: KSPSetWorkVecs(ksp,9);
18: return(0);
19: }
24: PetscErrorCode KSPSolve_MINRES(KSP ksp)
25: {
27: PetscInt i;
28: PetscScalar alpha,beta,ibeta,betaold,eta,c=1.0,ceta,cold=1.0,coold,s=0.0,sold=0.0,soold;
29: PetscScalar rho0,rho1,irho1,rho2,mrho2,rho3,mrho3,dp = 0.0;
30: PetscReal np;
31: Vec X,B,R,Z,U,V,W,UOLD,VOLD,WOLD,WOOLD;
32: Mat Amat,Pmat;
33: KSP_MINRES *minres = (KSP_MINRES*)ksp->data;
34: PetscBool diagonalscale;
37: PCGetDiagonalScale(ksp->pc,&diagonalscale);
38: if (diagonalscale) SETERRQ1(PetscObjectComm((PetscObject)ksp),PETSC_ERR_SUP,"Krylov method %s does not support diagonal scaling",((PetscObject)ksp)->type_name);
40: X = ksp->vec_sol;
41: B = ksp->vec_rhs;
42: R = ksp->work[0];
43: Z = ksp->work[1];
44: U = ksp->work[2];
45: V = ksp->work[3];
46: W = ksp->work[4];
47: UOLD = ksp->work[5];
48: VOLD = ksp->work[6];
49: WOLD = ksp->work[7];
50: WOOLD = ksp->work[8];
52: PCGetOperators(ksp->pc,&Amat,&Pmat);
54: ksp->its = 0;
56: VecSet(UOLD,0.0); /* u_old <- 0 */
57: VecCopy(UOLD,VOLD); /* v_old <- 0 */
58: VecCopy(UOLD,W); /* w <- 0 */
59: VecCopy(UOLD,WOLD); /* w_old <- 0 */
61: if (!ksp->guess_zero) {
62: KSP_MatMult(ksp,Amat,X,R); /* r <- b - A*x */
63: VecAYPX(R,-1.0,B);
64: } else {
65: VecCopy(B,R); /* r <- b (x is 0) */
66: }
68: KSP_PCApply(ksp,R,Z); /* z <- B*r */
70: VecDot(R,Z,&dp);
71: if (PetscRealPart(dp) < minres->haptol) {
72: PetscInfo2(ksp,"Detected indefinite operator %g tolerance %g\n",(double)PetscRealPart(dp),(double)minres->haptol);
73: ksp->reason = KSP_DIVERGED_INDEFINITE_MAT;
74: return(0);
75: }
77: dp = PetscAbsScalar(dp);
78: dp = PetscSqrtScalar(dp);
79: beta = dp; /* beta <- sqrt(r'*z */
80: eta = beta;
82: VecCopy(R,V);
83: VecCopy(Z,U);
84: ibeta = 1.0 / beta;
85: VecScale(V,ibeta); /* v <- r / beta */
86: VecScale(U,ibeta); /* u <- z / beta */
88: VecNorm(Z,NORM_2,&np); /* np <- ||z|| */
90: KSPLogResidualHistory(ksp,np);
91: KSPMonitor(ksp,0,np);
92: ksp->rnorm = np;
93: (*ksp->converged)(ksp,0,np,&ksp->reason,ksp->cnvP); /* test for convergence */
94: if (ksp->reason) return(0);
96: i = 0;
97: do {
98: ksp->its = i+1;
100: /* Lanczos */
102: KSP_MatMult(ksp,Amat,U,R); /* r <- A*u */
103: VecDot(U,R,&alpha); /* alpha <- r'*u */
104: KSP_PCApply(ksp,R,Z); /* z <- B*r */
106: VecAXPY(R,-alpha,V); /* r <- r - alpha v */
107: VecAXPY(Z,-alpha,U); /* z <- z - alpha u */
108: VecAXPY(R,-beta,VOLD); /* r <- r - beta v_old */
109: VecAXPY(Z,-beta,UOLD); /* z <- z - beta u_old */
111: betaold = beta;
113: VecDot(R,Z,&dp);
114: if ( PetscRealPart(dp) < minres->haptol) {
115: PetscInfo2(ksp,"Detected indefinite operator %g tolerance %g\n",(double)PetscRealPart(dp),(double)minres->haptol);
116: ksp->reason = KSP_DIVERGED_INDEFINITE_MAT;
117: break;
118: }
120: dp = PetscAbsScalar(dp);
121: beta = PetscSqrtScalar(dp); /* beta <- sqrt(r'*z) */
123: /* QR factorisation */
125: coold = cold; cold = c; soold = sold; sold = s;
127: rho0 = cold * alpha - coold * sold * betaold;
128: rho1 = PetscSqrtScalar(rho0*rho0 + beta*beta);
129: rho2 = sold * alpha + coold * cold * betaold;
130: rho3 = soold * betaold;
132: /* Givens rotation */
134: c = rho0 / rho1;
135: s = beta / rho1;
137: /* Update */
139: VecCopy(WOLD,WOOLD); /* w_oold <- w_old */
140: VecCopy(W,WOLD); /* w_old <- w */
142: VecCopy(U,W); /* w <- u */
143: mrho2 = -rho2;
144: VecAXPY(W,mrho2,WOLD); /* w <- w - rho2 w_old */
145: mrho3 = -rho3;
146: VecAXPY(W,mrho3,WOOLD); /* w <- w - rho3 w_oold */
147: irho1 = 1.0 / rho1;
148: VecScale(W,irho1); /* w <- w / rho1 */
150: ceta = c * eta;
151: VecAXPY(X,ceta,W); /* x <- x + c eta w */
152: eta = -s * eta;
154: VecCopy(V,VOLD);
155: VecCopy(U,UOLD);
156: VecCopy(R,V);
157: VecCopy(Z,U);
158: ibeta = 1.0 / beta;
159: VecScale(V,ibeta); /* v <- r / beta */
160: VecScale(U,ibeta); /* u <- z / beta */
162: np = ksp->rnorm * PetscAbsScalar(s);
164: ksp->rnorm = np;
165: KSPLogResidualHistory(ksp,np);
166: KSPMonitor(ksp,i+1,np);
167: (*ksp->converged)(ksp,i+1,np,&ksp->reason,ksp->cnvP); /* test for convergence */
168: if (ksp->reason) break;
169: i++;
170: } while (i<ksp->max_it);
171: if (i >= ksp->max_it) ksp->reason = KSP_DIVERGED_ITS;
172: return(0);
173: }
175: /*MC
176: KSPMINRES - This code implements the MINRES (Minimum Residual) method.
178: Options Database Keys:
179: . see KSPSolve()
181: Level: beginner
183: Notes: The operator and the preconditioner must be symmetric and the preconditioner must
184: be positive definite for this method.
185: Supports only left preconditioning.
187: Reference: Paige & Saunders, 1975.
189: Contributed by: Robert Scheichl: maprs@maths.bath.ac.uk
191: .seealso: KSPCreate(), KSPSetType(), KSPType (for list of available types), KSP, KSPCG, KSPCR
192: M*/
195: PETSC_EXTERN PetscErrorCode KSPCreate_MINRES(KSP ksp)
196: {
197: KSP_MINRES *minres;
201: KSPSetSupportedNorm(ksp,KSP_NORM_PRECONDITIONED,PC_LEFT,3);
202: PetscNewLog(ksp,&minres);
203: minres->haptol = 1.e-50;
204: ksp->data = (void*)minres;
206: /*
207: Sets the functions that are associated with this data structure
208: (in C++ this is the same as defining virtual functions)
209: */
210: ksp->ops->setup = KSPSetUp_MINRES;
211: ksp->ops->solve = KSPSolve_MINRES;
212: ksp->ops->destroy = KSPDestroyDefault;
213: ksp->ops->setfromoptions = 0;
214: ksp->ops->buildsolution = KSPBuildSolutionDefault;
215: ksp->ops->buildresidual = KSPBuildResidualDefault;
216: return(0);
217: }