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