# Refers to Algorithm 2 in the paper, it's needed to implement the algorithm 4. It doesn't work yet. Refer to the file testing.ipynb for more details. This function down here is just a place holder for now
# Refers to Algorithm 2 in the paper
defArnoldi(A,v0,m):
defArnoldi(A,v,m):
v=v0
beta=norm(v)
beta=norm(v)
v=v/beta
v=v/beta
h=sp.sparse.lil_matrix((m,m))
H=sp.sparse.lil_matrix((m+1,m))
V=sp.sparse.lil_matrix((A.shape[0],m+1))
V[:,0]=v# each column of V is a vector v
forjinrange(m):
forjinrange(m):
# print("j = ", j)
w=A@v
w=A@v
foriinrange(j):
foriinrange(j):
tmp=v.T@w
tmp=v.T@w# tmp is a 1x1 matrix, so it's O(1) in memory
h[i,j]=tmp[0,0]
H[i,j]=tmp[0,0]
w=w-h[i,j]*v
w=w-H[i,j]*v
h[j,j-1]=norm(w)# in the paper the index is referred as h[j+1,j] but since python starts from 0 it's h[j,j-1]
H[j+1,j]=norm(w)
ifh[j,j-1]==0:
ifH[j+1,j]==0:
print("Arnoldi breakdown")
print("Arnoldi breakdown")
m=j
m=j
v=0
v=0
break
break
else:
else:
v=w/h[j,j-1]
ifj<m-1:
returnv,h,m,beta,j# this is not the output required in the paper, I should return the matrix V and the matrix H