-
Notifications
You must be signed in to change notification settings - Fork 1
/
modCholesky.jl
63 lines (54 loc) · 1.23 KB
/
modCholesky.jl
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
eta = 1.0e-24#betweeen 0 to 1
function modchol(M, eta=1.0e-24)
# print(rank(M))
#Check for square matrix
m,n = size(M)
if m!=n
throw(DimensionMismatch("M should be a square Matrix"))
end
#Making Matrix sparse matrix
#M = sparse(M)
#initialization
M_prev = M
L = sparse(zeros((n,n)))
beta = maximum(diag(M)) #beta = max_i(1,...,m)M_ii
# print(beta)
for i = 1:m
if diag(M_prev)[i] <= beta*eta
# @show(diag(M_prev)[i], beta*eta)
#skip this elimination step
E_curr = get_E_for_modchol(M_prev, i)
M_curr = M_prev - E_curr
#@show(i, E_curr)
else
#perform usual Cholesky elimination step
L[i,i] = sqrt(M_prev[i,i])
M_curr = sparse(zeros(m,n))
for j = i+1:m
L[j,i] = M_prev[i,j] / L[i,i]
end
for j = i+1:m
for k = i+1:m
M_curr[j,k] = M_prev[j,k] - L[j,i]*L[k,i]
end
end
end
#@show(i,M_prev, M_curr, L)
M_prev = M_curr
end
return L
end
function get_E_for_modchol(M, index)
m,n = size(M)
E = sparse(zeros((m,n)))
for i = index:m
E[index,i] = M[index,i]
E[i,index] = M[i,index]
end
return E
end
#M = [25 15 -5;
# 15 18 0;
# -5 0 11]
#modchol(M)
#get_E_for_modchol(M, index)