``````import numpy as np

def gmres(A, b, x0, tol=1e-6, max_iter=10):
n = len(b)
x = np.copy(x0)
r = b - A @ x
d0 = np.linalg.norm(r)
v1 = r / d0
H = np.zeros((max_iter + 1, max_iter))  # Hessenberg-Matrix
V = np.zeros((n, max_iter + 1))          # Krylov-space
c = np.zeros(max_iter)
s = np.zeros(max_iter)

for k in range(1, max_iter):
V[:, k] = v1
v_tilde = A @ v1
for i in range(k):
H[i, k] = np.dot(V[:, i], v_tilde)
v_tilde -= H[i, k] * V[:, i]

omega = np.linalg.norm(v_tilde)

# Givens rotation
for i in range(k - 1):
h_ = c[i] * H[i, k] + s[i] * H[i + 1, k]
H[i + 1, k] = -s[i] * H[i, k] + c[i] * H[i + 1, k]
H[i, k] = h_

# Calculate Givens rotation parameters
if omega <= abs(H[k, k]):
tk = omega / abs(H[k, k])
c[k] = H[k, k] / (abs(H[k, k]) * np.sqrt(1 + tk**2))
s[k] = tk / np.sqrt(1 + tk**2)
else:
tk = H[k, k] / omega
c[k] = 1 / np.sqrt(1 + tk**2)
s[k] = tk * c[k]

# Update H matrix
H[k, k] = c[k] * H[k, k] + s[k] * omega

# Update residual norm
dk = -s[k] * d0
d0 = c[k] * d0

# Check for convergence
if abs(dk) < tol:
break

# Calculate solution
y = np.linalg.lstsq(H[:k + 1, :k + 1], np.array([d0]))[0]
x = x0 + np.dot(V[:, :k + 1], y)

return x, abs(dk)

# test
n = 3  # Dimension of Matrix
A = np.array([[4, 1, -1], [2, 5, 1], [1, 2, 4]])
b = np.array([8, 3, 11])
x0 = np.zeros_like(b)

solution, residual_norm = gmres(A, b, x0)
print("solution:", solution)
print("residuumnorm:", residual_norm)
``````

My output is always the zero vector.

So I was trying to implement the GMRES method. But I dont get any good solutions. I think there is a problem with the givens rotation. Because when I simply tried to tackle the problem via least square. I would get a good output/solution. Can someone help with the problem?

New contributor

Oschi is a new contributor to this site. Take care in asking for clarification, commenting, and answering.
Check out our Code of Conduct.