Coverage for /usr/share/miniconda3/envs/dolfin/lib/python3.8/site-packages/block/iterative/tfqmr.py: 4%
70 statements
« prev ^ index » next coverage.py v7.2.1, created at 2023-03-20 13:03 +0000
1from __future__ import division
2from .common import *
4def tfqmr(B, A, x, b, tolerance, maxiter, progress, relativeconv=False, callback=None):
5 #####
6 # Adapted from PyKrylov (https://github.com/dpo/pykrylov; LGPL license)
7 #####
9 r0 = b - A*x
11 rho = inner(r0,r0)
12 alphas = []
13 betas = []
14 residuals = [sqrt(rho)]
16 if relativeconv:
17 tolerance *= residuals[0]
19 if residuals[-1] < tolerance:
20 return x, residuals, [], []
22 y = r0.copy() # Initial residual vector
23 w = r0.copy()
24 d = 0*b
25 theta = 0.0
26 eta = 0.0
27 k = 0
29 z = B*y
30 u = A*z
31 v = u.copy()
33 while k < maxiter:
35 k += 1
36 progress += 1
37 sigma = inner(r0,v)
38 alpha = rho/sigma
40 # First pass
41 w -= alpha * u
42 d *= theta * theta * eta / alpha
43 d += z
45 residNorm = residuals[-1]
46 theta = norm(w)/residNorm
47 c = 1.0/sqrt(1 + theta*theta)
48 residNorm *= theta * c
49 eta = c * c * alpha
50 x += eta * d
51 m = 2.0 * k - 1.0
52 if residNorm * sqrt(m+1) < tolerance:
53 break
55 # Second pass
56 m += 1
57 y -= alpha * v
58 z = B*y
60 u = A*z
61 w -= alpha * u
62 d *= theta * theta * eta / alpha
63 d += z
64 theta = norm(w)/residNorm
65 c = 1.0/sqrt(1 + theta*theta)
66 residNorm *= theta * c
67 eta = c * c * alpha
68 x += eta * d
70 residual = residNorm * sqrt(m+1)
72 # Call user provided callback with solution
73 if callable(callback):
74 callback(k=k, x=x, r=residual)
76 residuals.append(residual)
77 if residual < tolerance or k >= maxiter:
78 break
80 # Final updates
81 rho_next = inner(r0,w)
82 beta = rho_next/rho
83 rho = rho_next
85 alphas.append(alpha)
86 betas.append(beta)
88 # Update y
89 y *= beta
90 y += w
92 # Partial update of v with current u
93 v *= beta
94 v += u
95 v *= beta
97 # Update u
98 z = B*y
99 u = A*z
101 # Complete update of v
102 v += u
104 return x, residuals, alphas, betas