-
Notifications
You must be signed in to change notification settings - Fork 1
/
L_Solve.py
65 lines (51 loc) · 1.78 KB
/
L_Solve.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
#!/usr/bin/python
from numpy import *
from scipy import linalg as la
######################## L_Solve #########################
## solves (A + sigma I)X = B using results from L_Seed) ##
##########################################################
def L_Solve(
seedSystem: dict, # ouput from L_Seed
B: ndarray, # RHS
σ: float, # shift >=0
tol = 5e-3, # absolute tolerance for CG
verbose= False, # verbose output
p_freq = 5 # print frequency (if verbose)
) -> ndarray:
U = seedSystem["U"]
δ = seedSystem["δ"] + σ
β = seedSystem["β"]
ρ = seedSystem["ρ"]
n = U.shape[0]
t = U.shape[1]
maxit = U.shape[2]
X = zeros((n,t,maxit)) ## approximate soln
P = zeros((n,t,maxit)) ## search directions
R = zeros((n,t,maxit)) ## residual
## coefficients
ω = zeros((t,maxit))
γ = ones((t,maxit))
## initial values
R[:,:,0] = B
P[:,:,0] = B
j = 0
cnvg = False
while j < maxit-1:
γ[:,j] = (δ[:,j] - ω[:,j-1]/γ[:,j-1])**-1
ω[:,j] = (β[:,j+1]*γ[:,j])**2
ρ[:,j+1] = -β[:,j+1]*γ[:,j]*ρ[:,j]
## CG vectors update
X[:,:,j+1] = X[:,:,j] + γ[:,j]*P[:,:,j]
R[:,:,j+1] = ρ[:,j+1]*U[:,:,j+1]
P[:,:,j+1] = R[:,:,j+1] + ω[:,j]*P[:,:,j]
j += 1
# res_norm = amax(abs(R[:,:,j]))
res_norm = max(apply_along_axis(la.norm,0, R[:,:,j]))
if verbose and j % p_freq ==0 : print("Error at step ",j,
" is ", res_norm)
if res_norm <= tol:
cnvg = True
break
if verbose and cnvg: print("Converged after ",j," iterations.")
elif verbose: print("Failed to converge after ",j," iterations.")
return X[:,:,j]