Coverage for /usr/share/miniconda3/envs/dolfin/lib/python3.8/site-packages/block/iterative/minres.py: 87%
107 statements
« prev ^ index » next coverage.py v7.2.1, created at 2023-03-20 13:03 +0000
1from __future__ import division
2from __future__ import print_function
3from .common import *
5def minres(B, A, x, b, tolerance, maxiter, progress, relativeconv=False, shift=0, callback=None):
6 #####
7 # Adapted from PyKrylov (https://github.com/dpo/pykrylov; LGPL license)
8 #####
10 msg = [' beta1 = 0. The exact solution is x = 0 ', # 0
11 ' A solution to Ax = b was found, given rtol ', # 1
12 ' A least-squares solution was found, given rtol ', # 2
13 ' Reasonable accuracy achieved, given eps ', # 3
14 ' x has converged to an eigenvector ', # 4
15 ' acond has exceeded 0.1/eps ', # 5
16 ' The iteration limit was reached ', # 6
17 ' Aname does not define a symmetric matrix ', # 7
18 ' Mname does not define a symmetric matrix ', # 8
19 ' Mname does not define a pos-def preconditioner ', # 9
20 ' Userdefined callback function returned true ', #10
21 ' beta2 = 0. If M = I, b and x are eigenvectors '] #-1
23 callback_converged = False
24 istop = 0; itn = 0; Anorm = 0.0; Acond = 0.0;
25 rnorm = 0.0; ynorm = 0.0; done = False;
27 #------------------------------------------------------------------
28 # Set up y and v for the first Lanczos vector v1.
29 # y = beta1 P' v1, where P = C**(-1).
30 # v is really P' v1.
31 #------------------------------------------------------------------
32 r1 = b - A*x
33 y = B*r1
34 beta1 = inner(r1,y)
36 if relativeconv:
37 tolerance *= sqrt(beta1)
38# print ("tolerance ", tolerance, beta1, relativeconv)
42 # Test for an indefinite preconditioner.
43 # If b = 0 exactly, stop with x = 0.
44 if beta1 < 0: 44 ↛ 45line 44 didn't jump to line 45, because the condition on line 44 was never true
45 raise ValueError('B does not define a pos-def preconditioner')
46 if beta1 == 0: 46 ↛ 47line 46 didn't jump to line 47, because the condition on line 46 was never true
47 x *= 0
48 return x, [0], [], []
50 beta1 = sqrt(beta1); # Normalize y to get v1 later.
51 residuals = [beta1]
54 # -------------------------------------------------------------------
55 # Initialize other quantities.
56 # ------------------------------------------------------------------
57 oldb = 0.0; beta = beta1; dbar = 0.0; epsln = 0.0
58 qrnorm = beta1; phibar = beta1; rhs1 = beta1; Arnorm = 0.0
59 rhs2 = 0.0; tnorm2 = 0.0; ynorm2 = 0.0
60 cs = -1.0; sn = 0.0
61 w = 0*b
62 w2 = 0*b
63 r2 = r1.copy()
65 # ---------------------------------------------------------------------
66 # Main iteration loop.
67 # --------------------------------------------------------------------
68 while itn < maxiter: 68 ↛ 209line 68 didn't jump to line 209, because the condition on line 68 was never false
69 itn += 1
70 progress += 1
72 # -------------------------------------------------------------
73 # Obtain quantities for the next Lanczos vector vk+1, k=1,2,...
74 # The general iteration is similar to the case k=1 with v0 = 0:
75 #
76 # p1 = Operator * v1 - beta1 * v0,
77 # alpha1 = v1'p1,
78 # q2 = p2 - alpha1 * v1,
79 # beta2^2 = q2'q2,
80 # v2 = (1/beta2) q2.
81 #
82 # Again, y = betak P vk, where P = C**(-1).
83 # .... more description needed.
84 # -------------------------------------------------------------
85 s = 1.0/beta # Normalize previous vector (in y).
86 v = s*y # v = vk if P = I
88 y = A*v
89 if shift: 89 ↛ 90line 89 didn't jump to line 90, because the condition on line 89 was never true
90 y -= shift*v
92 if itn >= 2:
93 y -= (beta/oldb)*r1
95 alfa = inner(v,y) # alphak
96 y -= alfa/beta*r2
97 r1 = r2
98 r2 = y
99 y = B*r2
101 oldb = beta # oldb = betak
102 beta = inner(r2,y) # beta = betak+1^2
103 if beta < 0: 103 ↛ 104line 103 didn't jump to line 104, because the condition on line 103 was never true
104 raise ValueError('B does not define a pos-def preconditioner')
105 beta = sqrt(beta)
106 tnorm2 += alfa**2 + oldb**2 + beta**2
108 if itn==1: # Initialize a few things.
109 if beta/beta1 <= 10*eps: # beta2 = 0 or ~ 0. 109 ↛ 110line 109 didn't jump to line 110, because the condition on line 109 was never true
110 istop = -1 # Terminate later.
112 # tnorm2 = alfa**2 ??
113 gmax = abs(alfa) # alpha1
114 gmin = gmax # alpha1
116 # Apply previous rotation Qk-1 to get
117 # [deltak epslnk+1] = [cs sn][dbark 0 ]
118 # [gbar k dbar k+1] [sn -cs][alfak betak+1].
120 oldeps = epsln
121 delta = cs * dbar + sn * alfa # delta1 = 0 deltak
122 gbar = sn * dbar - cs * alfa # gbar 1 = alfa1 gbar k
124 # Note: There is severe cancellation in the computation of gbar
125 #print ' sn = %21.15e\n dbar = %21.15e\n cs = %21.15e\n alfa = %21.15e\n sn*dbar-cs*alfa = %21.15e\n gbar =%21.15e' % (sn, dbar, cs, alfa, sn*dbar-cs*alfa, gbar)
127 epsln = sn * beta # epsln2 = 0 epslnk+1
128 dbar = - cs * beta # dbar 2 = beta2 dbar k+1
129 root = sqrt(gbar**2 + dbar**2)
130 Arnorm = phibar * root
132 # Compute the next plane rotation Qk
134 gamma = sqrt(gbar**2 + beta**2)
135 gamma = max(gamma, eps)
136 cs = gbar / gamma # ck
137 sn = beta / gamma # sk
138 phi = cs * phibar # phik
139 phibar = sn * phibar # phibark+1
141 # Update x.
143 denom = 1.0/gamma
144 w1 = w2
145 w2 = w
146 w = (v - oldeps*w1 - delta*w2) * denom
147 x += phi*w
149 # Go round again.
151 gmax = max(gmax, gamma)
152 gmin = min(gmin, gamma)
153 z = rhs1 / gamma
154 ynorm2 = z**2 + ynorm2
155 rhs1 = rhs2 - delta*z
156 rhs2 = - epsln*z
158 # Estimate various norms and test for convergence.
160 Anorm = sqrt(tnorm2)
161 ynorm = sqrt(ynorm2)
162 epsa = Anorm * eps
163 epsx = Anorm * ynorm * eps
164 #epsr = Anorm * ynorm * tolerance
165 diag = gbar
166 if diag==0: diag = epsa
168 qrnorm = phibar
169 rnorm = qrnorm
170 test1 = rnorm / (Anorm*ynorm) # ||r|| / (||A|| ||x||)
171 test2 = root / Anorm # ||Ar|| / (||A|| ||r||)
173 # Estimate cond(A).
174 # In this version we look at the diagonals of R in the
175 # factorization of the lower Hessenberg matrix, Q * H = R,
176 # where H is the tridiagonal matrix from Lanczos with one
177 # extra row, beta(k+1) e_k^T.
179 Acond = gmax/gmin
181 # Call user provided callback with solution
182 if callable(callback): 182 ↛ 183line 182 didn't jump to line 183, because the condition on line 182 was never true
183 callback_converged = callback(k=itn, x=x, r=test1)
185 # See if any of the stopping criteria are satisfied.
186 # In rare cases istop is already -1 from above (Abar = const*I)
188 if istop == 0: 188 ↛ 204line 188 didn't jump to line 204, because the condition on line 188 was never false
189 t1 = 1 + test1 # These tests work if rtol < eps
190 t2 = 1 + test2
191 if t2 <= 1: istop = 2
192 if t1 <= 1: istop = 1
194 if itn >= maxiter: istop = 6
195 if Acond >= 0.1/eps: istop = 4
196 if epsx >= beta1: istop = 3
197 # if rnorm <= epsx: istop = 2
198 # if rnorm <= epsr: istop = 1
199 if test2 <= tolerance: istop = 2
200 if test1 <= tolerance: istop = 1
202 if callback_converged: istop = 10
204 residuals.append(test1)
206 if istop != 0:
207 break
209 if istop != 1: 209 ↛ 210line 209 didn't jump to line 210, because the condition on line 209 was never true
210 print(('MinRes:', msg[istop]))
212 return x, residuals, [], []