# ## MIT License

# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to deal
# in the Software without restriction, including without limitation the rights
# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
# copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
#
# The above copyright notice and this permission notice shall be included in
# all copies or substantial portions of the Software.
#
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
# THE SOFTWARE.
#

# In[1]:
import numpy as np
import numpy.linalg as la
import matplotlib.pyplot as pt
# Let us make a matrix with a defined set of eigenvalues and eigenvectors, given by `eigvals` and `eigvecs`.
# In[2]:
np.random.seed(40)
# Generate matrix with eigenvalues 1...25
n = 25
eigvals = np.linspace(1., n, n)
eigvecs = np.random.randn(n, n)
print(eigvals)
A = la.solve(eigvecs, np.dot(np.diag(eigvals), eigvecs))
print(la.eig(A)[0])
# ## Initialization
# Set up $Q$ and $H$:
# In[3]:
Q = np.zeros((n, n))
H = np.zeros((n, n))
k = 0
# Pick a starting vector, normalize it
# In[4]:
x0 = np.random.randn(n)
x0 = x0/la.norm(x0)
# Poke it into the first column of Q
Q[:, k] = x0
del x0
# Make a list to save arrays of Ritz values:
# In[5]:
ritz_values = []
# ## Algorithm
# Carry out one iteration of Arnoldi iteration.
#
# Run this cell in-place (Ctrl-Enter) until H is filled.
# In[30]:
print(k)
u = A @ Q[:, k]
# Carry out Gram-Schmidt on u against Q
for j in range(k+1):
qj = Q[:, j]
H[j,k] = qj @ u
u = u - H[j,k]*qj
if k+1 < n:
H[k+1, k] = la.norm(u)
Q[:, k+1] = u/H[k+1, k]
k += 1
pt.spy(H)
ritz_values.append(la.eig(H)[0])
# Check that $Q^T A Q =H$:
# In[31]:
la.norm(Q.T @ A @ Q - H)/ la.norm(A)
# Check that Q is orthogonal:
# In[32]:
la.norm(Q.T @ Q - np.eye(n))
# ## Plot convergence of Ritz values
# Enable the Ritz value collection above to make this work.
# In[33]:
for i, rv in enumerate(ritz_values):
pt.plot([i] * len(rv), rv, "x")
# In[ ]: