-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathOrtho.py
103 lines (95 loc) · 3.03 KB
/
Ortho.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
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
from scipy.linalg.blas import dgemm
from numpy import matmul
from scipy.linalg import eigh, norm
from numpy import empty, eye, finfo, diag, max, abs
def OrthoB(U, B, V, BU, norm_BV,
tol_ortho=100*finfo(float).eps,
tol_replace=10*finfo(float).eps,
itmax1=6, itmax2=6):
"""
Stathopoulos, A., & Wu, K. (2002)
A block orthogonalization procedure with constant synchronization requirements
SIAM Journal on Scientific Computing, 23(6), 2165-2182
"""
_, p = BU.shape
_, q = V.shape
UtBU = empty((p, p))
VtBU = empty((q, p))
matmul(V.T, BU, out=VtBU)
for _ in range(itmax1):
U[:] = dgemm(alpha=-1.0, a=V, b=VtBU, c=U, overwrite_c=True, beta=1.0)
BU[:] = B @ U
for _ in range(itmax2):
svqbB(U, BU, tol_replace)
matmul(U.T, BU, out=UtBU)
norm_U = norm(U)
err = norm(UtBU - eye(p)) / (norm(BU) * norm_U)
if err < tol_ortho:
break
matmul(V.T, BU, out=VtBU)
err = norm(VtBU) / (norm_BV * norm_U)
if err < tol_ortho:
break
def svqbB(U, BU, tol):
"""
Stathopoulos, A., & Wu, K. (2002)
A block orthogonalization procedure with constant synchronization requirements
SIAM Journal on Scientific Computing, 23(6), 2165-2182
"""
_, p = U.shape
UtBU = empty((p, p))
matmul(U.T, BU, out=UtBU)
D = diag(diag(UtBU) ** -0.5)
matmul(D, matmul(UtBU, D), out=UtBU) # UtBU[:] = D @ UtBU @ D
Theta, Z = eigh(UtBU)
theta_abs_max = max(abs(Theta))
Theta[Theta < tol * theta_abs_max] = tol * theta_abs_max
# Z[:] = D @ Z @ diag(Theta ** -0.5)
matmul(D, matmul(Z, diag(Theta ** -0.5)), out=Z)
U[:] = matmul(U, Z)
BU[:] = matmul(BU, Z)
def Ortho(U, V,
tol_ortho=100*finfo(float).eps,
tol_replace=10*finfo(float).eps,
itmax1=6, itmax2=6):
"""
Stathopoulos, A., & Wu, K. (2002)
A block orthogonalization procedure with constant synchronization requirements
SIAM Journal on Scientific Computing, 23(6), 2165-2182
"""
_, p = U.shape
_, q = V.shape
UtU = empty((p, p))
VtU = empty((q, p))
matmul(V.T, U, out=VtU)
norm_V = norm(V)
for _ in range(itmax1):
U[:] = dgemm(alpha=-1.0, a=V, b=VtU, c=U, overwrite_c=True, beta=1.0)
for _ in range(itmax2):
svqb(U, tol_replace)
norm_U = norm(U)
matmul(U.T, U, out=UtU)
err = norm(UtU - eye(p)) / (norm_U ** 2)
if err < tol_ortho:
break
matmul(V.T, U, out=VtU)
err = norm(VtU) / (norm_V * norm_U)
if err < tol_ortho:
break
def svqb(U, tol):
"""
Stathopoulos, A., & Wu, K. (2002)
A block orthogonalization procedure with constant synchronization requirements
SIAM Journal on Scientific Computing, 23(6), 2165-2182
"""
_, p = U.shape
UtU = empty((p, p))
matmul(U.T, U, out=UtU)
D = diag(diag(UtU) ** -0.5)
matmul(D, matmul(UtU, D), out=UtU) # UtU[:] = D @ UtU @ D
Theta, Z = eigh(UtU)
theta_abs_max = max(abs(Theta))
Theta[Theta < tol * theta_abs_max] = tol * theta_abs_max
# Z[:] = D @ Z @ diag(Theta ** -0.5)
matmul(D, matmul(Z, diag(Theta ** -0.5)), out=Z)
U[:] = matmul(U, Z)