Commit d77d494e authored by Charles Bouillaguet's avatar Charles Bouillaguet
Browse files

toy bl

parent 21cecf2b
# toy implementation of the block-lanczos algorithm
# given a (sparse) symetric matrix A, computes a (block of) vector x s.t.
# x * A * x == 0
# x * M * M.t * x == 0
n = 16
ncols = 7*n + n // 2
N = 8*n
M = matrix.random(GF(2), N, N)
for j in range(ncols, N):
M[:, j] = 0
# gaussian elimination --> independent columns d_i et w_i
def elimination(M_):
d = matrix.zero(GF(2), n)
# first pass: find dependent columns (compute d)
M = matrix(M_)
for j in range(n):
pivot = -1
for i in range(n):
if d[i, i] == 0 and M[i, j] == 1:
pivot = i
break
if pivot < 0:
continue
d[j, j] = 1
M.swap_rows(pivot, j)
for i in range(n):
if i != j and M[i, j] == 1:
M[i] += M[j]
# second pass: compute winv
M = d * M_ * d
winv = matrix(d)
for j in range(n):
pivot = -1
for i in range(j, n):
if M[i, j] == 1:
pivot = i
break
if pivot < 0:
continue
M.swap_rows(pivot, j)
winv.swap_rows(pivot, j)
assert M[j, j] == 1
for i in range(n):
if i != j and M[i, j] == 1:
M[i] += M[j]
winv[i] += winv[j]
return winv, d
# starting point
v = matrix.random(GF(2), N, n)
p = matrix.zero(GF(2), N, n)
W = []
while True:
Mv = M.transpose()*v # block
Av = M*Mv # block
vtAv = v.transpose()*Av # n * n
vtAAv = Av.transpose()*Av # n * n
winv, d = elimination(vtAv)
#assert winv == winv*d == d*winv
#assert d == winv * vtAv*d
#assert winv + winv.transpose() == 0
# checking orthogonality w.r.t. previous iterations
wi = v*d
for wj in W:
assert wj.transpose() * M * M.transpose() * wi == 0
W.append(wi)
if d == 0:
break
# next iterate
dbar = matrix.identity(GF(2), n) - d
c = winv * (vtAAv*d + vtAv*dbar) # n*n
nextv = Av*d + v*dbar - v*c - p*vtAv*d # block
p = v*winv + p*dbar # block
v = nextv
assert v != 0
assert v.transpose() * M == 0
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment