Skip to content

Instantly share code, notes, and snippets.

@angellicadearaujo
Created November 7, 2016 12:15
Show Gist options
  • Select an option

  • Save angellicadearaujo/9394bde6d99e8ee6d36fd1dc1eab652e to your computer and use it in GitHub Desktop.

Select an option

Save angellicadearaujo/9394bde6d99e8ee6d36fd1dc1eab652e to your computer and use it in GitHub Desktop.
Algoritmo de Arnoldi
def Arnoldi(A,b):
n=len(A)
# Build S set
q=[[] for j in range(n+1)]
S=[[] for j in range(n+1)]
q[0]=b
for i in range(n):
S[i]=multiply(A,q[i])
tmp=[0 for l in range(len(S[i]))]
for j in range(0,i):
tmp=subtract(tmp,Proj(q[j],S[i]))
q[i+1]=addition(S[i], tmp) # ui = vi - proj(uj, vi)
# Build orthogonal with GramSchmidt
n=len(S)
u=[[] for j in range(n)]
u[0]=S[0]
for i in range(1,n):
tmp=[0 for l in range(len(S[i]))]
for j in range(0,i):
tmp=subtract(tmp,Proj(u[j],S[i]))
u[i]=addition(S[i], tmp)
# Build Hessenber Matrix
n=len(S)
H=[[0 for j in range(n)] for i in range(n)]
i=0
while (i+1) < n:
H[i+1][i]=norm(u[i]) # Subdiagonal
j=0
while j<=i:
H[j][i]=internProduct(u[j],u[j])
H[j][i]=internProduct(S[i],u[j])*(1/float(H[j][i])) # hj,n = qj * vi
j += 1
i+=1
return H
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment