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.

1

Khám phá các thẻ bài đăng