Coverage for /usr/share/miniconda3/envs/dolfin/lib/python3.8/site-packages/block/iterative/symmlq.py: 3%
124 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 symmlq(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={
11 -1:' beta2 = 0. If M = I, b and x are eigenvectors',
12 0:' beta1 = 0. The exact solution is x = 0',
13 1:' Requested accuracy achieved, as determined by tolerance',
14 2:' Reasonable accuracy achieved, given eps',
15 3:' x has converged to an eigenvector',
16 4:' acond has exceeded 0.1/eps',
17 5:' The iteration limit was reached',
18 6:' aprod does not define a symmetric matrix',
19 7:' msolve does not define a symmetric matrix',
20 8:' msolve does not define a pos-def preconditioner'}
22 istop = 0
23 w = 0*b
25 # Set up y for the first Lanczos vector v1.
26 # y is really beta1 * P * v1 where P = C^(-1).
27 # y and beta1 will be zero if b = 0.
29 r1 = 1*b
30 y = B*r1
32 beta1 = inner(r1, y)
34 # Test for an indefinite preconditioner.
35 # If b = 0 exactly, stop with x = 0.
37 if beta1 < 0:
38 raise ValueError('B does not define a pos-def preconditioner')
39 if beta1 == 0:
40 x *= 0
41 return x, [0], [], []
43 beta1 = sqrt(beta1)
44 s = 1.0 / beta1
45 v = s * y
47 y = A*v
49 # Set up y for the second Lanczos vector.
50 # Again, y is beta * P * v2 where P = C^(-1).
51 # y and beta will be zero or very small if Abar = I or constant * I.
53 if shift:
54 y -= shift * v
55 alfa = inner(v, y)
56 y -= (alfa / beta1) * r1
58 # Make sure r2 will be orthogonal to the first v.
60 z = inner(v, y)
61 s = inner(v, v)
62 y -= (z / s) * v
63 r2 = y
64 y = B*y
66 oldb = beta1
67 beta = inner(r2, y)
68 if beta < 0:
69 raise ValueError('B does not define a pos-def preconditioner')
71 # Cause termination (later) if beta is essentially zero.
73 beta = sqrt(beta)
74 if beta <= eps:
75 istop = -1
77 # Initialize other quantities.
78 rhs2 = 0
79 tnorm = alfa**2 + beta**2
80 gbar = alfa
81 dbar = beta
83 bstep = 0
84 ynorm2 = 0
85 snprod = 1
87 gmin = gmax = abs(alfa) + eps
88 rhs1 = beta1
90 # ------------------------------------------------------------------
91 # Main iteration loop.
92 # ------------------------------------------------------------------
93 # Estimate various norms and test for convergence.
95 alphas = []
96 betas = []
97 residuals = [beta1/sqrt(tnorm)]
98 itn = 0
100 while True:
101 itn += 1
102 progress += 1
103 anorm = sqrt(tnorm)
104 ynorm = sqrt(ynorm2)
105 epsa = anorm * eps
106 epsx = anorm * ynorm * eps
107 epsr = anorm * ynorm * tolerance
108 diag = gbar
110 if diag == 0: diag = epsa
112 lqnorm = sqrt(rhs1**2 + rhs2**2)
113 qrnorm = snprod * beta1
114 cgnorm = qrnorm * beta / abs(diag)
116 # Estimate Cond(A).
117 # In this version we look at the diagonals of L in the
118 # factorization of the tridiagonal matrix, T = L*Q.
119 # Sometimes, T(k) can be misleadingly ill-conditioned when
120 # T(k+1) is not, so we must be careful not to overestimate acond
122 if lqnorm < cgnorm:
123 acond = gmax / gmin
124 else:
125 acond = gmax / min(gmin, abs(diag))
127 zbar = rhs1 / diag
128 z = (snprod * zbar + bstep) / beta1
130 # See if any of the stopping criteria are satisfied.
131 # In rare cases, istop is already -1 from above
132 # (Abar = const * I).
134 if istop == 0:
135 if itn > maxiter : istop = 5
136 if acond >= 0.1/eps : istop = 4
137 if epsx >= beta1 : istop = 3
138 if cgnorm <= epsx : istop = 2
139 if cgnorm <= epsr : istop = 1
141 residuals.append(cgnorm / anorm / (ynorm or 1))
143 if istop !=0:
144 break
146 # Obtain the current Lanczos vector v = (1 / beta)*y
147 # and set up y for the next iteration.
149 s = 1/beta
150 v = s * y
151 y = A*v
152 if shift:
153 y -= shift * v
154 y -= (beta / oldb) * r1
155 alfa = inner(v, y)
156 y -= (alfa / beta) * r2
157 r1 = r2.copy()
158 r2 = y
159 y = B*y
160 oldb = beta
161 beta = inner(r2, y)
163 alphas.append(alfa)
164 betas.append(beta)
166 if beta < 0:
167 raise ValueError('A does not define a symmetric matrix')
169 beta = sqrt(beta);
170 tnorm = tnorm + alfa**2 + oldb**2 + beta**2;
172 # Compute the next plane rotation for Q.
174 gamma = sqrt(gbar**2 + oldb**2)
175 cs = gbar / gamma
176 sn = oldb / gamma
177 delta = cs * dbar + sn * alfa
178 gbar = sn * dbar - cs * alfa
179 epsln = sn * beta
180 dbar = - cs * beta
182 # Update X.
184 z = rhs1 / gamma
185 s = z*cs
186 t = z*sn
187 x += s*w + t*v
188 w *= sn
189 w -= cs*v
191 # Accumulate the step along the direction b, and go round again.
193 bstep = snprod * cs * z + bstep
194 snprod = snprod * sn
195 gmax = max(gmax, gamma)
196 gmin = min(gmin, gamma)
197 ynorm2 = z**2 + ynorm2
198 rhs1 = rhs2 - delta * z
199 rhs2 = - epsln * z
201 # ------------------------------------------------------------------
202 # End of main iteration loop.
203 # ------------------------------------------------------------------
205 # Move to the CG point if it seems better.
206 # In this version of SYMMLQ, the convergence tests involve
207 # only cgnorm, so we're unlikely to stop at an LQ point,
208 # EXCEPT if the iteration limit interferes.
210 if cgnorm < lqnorm:
211 zbar = rhs1 / diag
212 bstep = snprod * zbar + bstep
213 ynorm = sqrt(ynorm2 + zbar**2)
214 x += zbar * w
216 # Add the step along b.
218 bstep = bstep / beta1
219 y = B*b
220 x += bstep * y
222 if istop != 1:
223 print(('SymmLQ:',msg[istop]))
225 return x, residuals, [], []