@ti.data_oriented
class GMRES_Solver:
    def __init__(self, rows, cols, n_iter=1, tol=1e-6):
        self.rows = rows
        self.cols = cols
        self.n_iter = n_iter
        self.tol = tol  # Convergence tolerance

        # Taichi fields for 2D vectors and matrices
        self.x0 = ti.field(dtype=ti.f32, shape=(rows, cols))
        self.b = ti.field(dtype=ti.f32, shape=(rows, cols))
        self.V = ti.field(dtype=ti.f32, shape=(n_iter + 1, rows, cols))
        self.H = ti.field(dtype=ti.f32, shape=(n_iter + 1, n_iter))
        self.w = ti.field(dtype=ti.f32, shape=(rows, cols))
        self.beta = ti.field(dtype=ti.f32, shape=())
        self.r0 = ti.field(dtype=ti.f32, shape=(rows, cols))
        self.diagonal_inv = ti.field(dtype=ti.f32, shape=(rows, cols))
    # Kernel for vector subtraction: result = v1 - v2
    @ti.kernel
    def subtract_vectors(self, v1: ti.template(), v2: ti.template(), result: ti.template()):
        for i, j in ti.ndrange(self.rows, self.cols):
            result[i, j] = v1[i, j] - v2[i, j]

    # Kernel for vector normalization: v /= norm(v)
    @ti.kernel
    def normalize(self, v: ti.template(), norm: ti.f32):
        for i, j in ti.ndrange(self.rows, self.cols):
            v[i, j] /= norm

    # Kernel for calculating the 2-norm of a vector
    @ti.kernel
    def norm(self, v: ti.template()) -> ti.f32:
        norm_val = 0.0
        for i, j in ti.ndrange(self.rows, self.cols):
            norm_val += v[i, j] * v[i, j]
        return ti.sqrt(norm_val)

    # Kernel to copy one vector to another (manual element-by-element copying)
    @ti.kernel
    def copy_vector(self, src: ti.template(), dst: ti.template(), idx: ti.i32):
        for i, j in ti.ndrange(self.rows, self.cols):
            dst[idx, i, j] = src[i, j]

    # Kernel for dot product with 2D array row access
    @ti.kernel
    def dot_product_row(self, v1: ti.template(), v2: ti.template(), row: ti.i32) -> ti.f32:
        dot = 0.0
        for i, j in ti.ndrange(self.rows, self.cols):
            dot += v1[i, j] * v2[row, i, j]
        return dot

    # Kernel for in-place subtraction: v1 -= scalar * v2
    @ti.kernel
    def subtract_scaled_vector(self, v1: ti.template(), v2: ti.template(), scalar: ti.f32):
        for i, j in ti.ndrange(self.rows, self.cols):
            v1[i, j] -= scalar * v2[i, j]

    # Kernel for dot product
    @ti.kernel
    def dot_product(self, v1: ti.template(), v2: ti.template()) -> ti.f32:
        dot = 0.0
        for i, j in ti.ndrange(self.rows, self.cols):
            dot += v1[i, j] * v2[i, j]
        return dot

    # Helper kernel to extract the j-th row from V and store it in a temporary vector
    @ti.kernel
    def extract_vector_from_V(self, j: ti.i32, temp: ti.template()):
        for i, k in ti.ndrange(self.rows, self.cols):
            temp[i, k] = self.V[j, i, k]

    # Orthogonalization in Python-scope (modified)
    def orthogonalize_in_python(self, j):
        temp_vec = ti.field(dtype=ti.f32, shape=(self.rows, self.cols))

        for i in range(j + 1):
            self.extract_vector_from_V(i, temp_vec)  # Extract V[i] into temp_vec
            self.H[i, j] = self.dot_product(self.w, temp_vec)  # Perform dot product
            self.subtract_scaled_vector(self.w, temp_vec, self.H[i, j])  # Update w

        self.H[j + 1, j] = self.norm(self.w)  # Norm of w
        self.normalize(self.w, self.H[j + 1, j])  # Normalize w
        self.copy_vector(self.w, self.V, j + 1)  # Copy w to V[j + 1]
    # Function to apply Givens rotations and solve the least squares problem
    # Kernel for in-place Givens rotations and solving least squares problem
    @ti.kernel
    def apply_givens_rotations(self, e1: ti.template(), y: ti.template()):
        # Apply Givens rotations and solve the upper triangular system R_n y_n = e1
        for j in range(self.n_iter):
            rho = self.H[j, j]  # diagonal entry of H (R_n)
            sigma = self.H[j + 1, j]  # sub-diagonal entry of H
            c = rho / ti.sqrt(rho**2 + sigma**2)
            s = sigma / ti.sqrt(rho**2 + sigma**2)

            # Modify the matrix H to reflect the Givens rotation
            for k in range(j, self.n_iter):
                temp = c * self.H[j, k] + s * self.H[j + 1, k]
                self.H[j + 1, k] = -s * self.H[j, k] + c * self.H[j + 1, k]
                self.H[j, k] = temp

            # Apply the same rotation to e1 (right-hand side vector)
            temp_e1 = c * e1[j] + s * e1[j + 1]
            e1[j + 1] = -s * e1[j] + c * e1[j + 1]
            e1[j] = temp_e1

        # Step 2: Solve the upper triangular system R_n y_n = e1 (simulating back substitution)
        for i in range(self.n_iter):
            rev_i = self.n_iter - 1 - i  # Reverse index for back substitution
            sum_val = 0.0
            for j in range(rev_i + 1, self.n_iter):
                sum_val += self.H[rev_i, j] * y[j]
            y[rev_i] = (e1[rev_i] - sum_val) / self.H[rev_i, rev_i]
    # Kernel for preconditioned matrix-vector multiplication: v' = D^{-1} v
    @ti.kernel
    def apply_preconditioner(self, v: ti.template(), result: ti.template()):
        for i, j in ti.ndrange(self.rows, self.cols):
            result[i, j] = self.diagonal_inv[i, j] * v[i, j]  # Applying diagonal inverse

    # GMRES solver function with Jacobi preconditioner
    def gmres_solve(self, b, x0, linear_map, diagonal):
        # Set up the preconditioner (diagonal inverse)
        for i in range(self.rows):
            for j in range(self.cols):
                if diagonal[i,j] != 0:
                    self.diagonal_inv[i, j] = 1.0 / diagonal[i, j]  # Inverse of diagonal elements

        # Step 1: r0 = b - linear_map(x0)
        linear_map(x0, self.r0)
        self.subtract_vectors(b, self.r0, self.r0)

        # Apply Jacobi preconditioner to r0: r0 = D^{-1} r0
        self.apply_preconditioner(self.r0, self.r0)

        # Step 2: beta = norm(r0)
        beta_value = self.norm(self.r0)
        self.beta[None] = beta_value
        # Step 3: V[0] = r0 / beta
        self.normalize(self.r0, beta_value)
        self.copy_vector(self.r0, self.V, 0)

        # Step 4: GMRES iteration loop
        for j in range(self.n_iter):
            temp_vec = ti.field(dtype=ti.f32, shape=(self.rows, self.cols))
            self.extract_vector_from_V(j, temp_vec)

            # Apply preconditioned matrix-vector product: w = D^{-1} * linear_map(V[j])
            linear_map(temp_vec, self.w)  # Compute w = linear_map(V[j])
            self.apply_preconditioner(self.w, self.w)  # Apply Jacobi preconditioner

            self.orthogonalize_in_python(j)  # Orthogonalize
        e1 = ti.field(dtype=ti.f32, shape=(self.n_iter + 1))
        e1[0] = self.beta[None]

        y = ti.field(dtype=ti.f32, shape=(self.n_iter))

        # Apply Givens rotations and solve least squares
        self.apply_givens_rotations(e1, y)

        # Compute the new x
        x_new = ti.field(dtype=ti.f32, shape=(self.rows, self.cols))
        for i in range(self.rows):
            for j in range(self.cols):
                sum_val = 0.0
                for k in range(self.n_iter):
                    sum_val += self.V[k, i, j] * y[k]
                x_new[i, j] = x0[i, j] + sum_val

        return x_new
Make it run faster