lab2

import numpy as np

def parse_input(filename):
    """Parse input file with matrices A1, y1, A2, y2"""
    systems = []
    with open(filename, 'r') as f:
        lines = [line.strip() for line in f if line.strip()]

    i = 0
    while i < len(lines):
        if lines[i].startswith('A'):
            # Read matrix A
            A = []
            i += 1
            while i < len(lines) and not lines[i].startswith('y') and not lines[i].startswith('A'):
                row = list(map(float, lines[i].split()))
                A.append(row)
                i += 1
            A = np.array(A)
        if i < len(lines) and lines[i].startswith('y'):
            # Read vector y
            i += 1
            y = np.array(list(map(float, lines[i].split())))
            i += 1
            systems.append((A, y))

    return systems


def iterative_solve(A, y, eps=1e-6, max_iter=1000):
    """
    Solve Ax = y using Jacobi iterative method.
    Transform to x = y' - B*x where:
      - Divide each equation i by A[i,i]
      - B[i,j] = A[i,j]/A[i,i] for i != j, B[i,i] = 0
      - y'[i] = y[i]/A[i,i]
    """
    n = len(y)

    # Build B and y_tilde: x = y_tilde - B*x
    B = np.zeros((n, n))
    y_tilde = np.zeros(n)

    for i in range(n):
        if abs(A[i, i]) < 1e-15:
            raise ValueError(f"Zero diagonal element at row {i}, cannot apply Jacobi method")
        y_tilde[i] = y[i] / A[i, i]
        for j in range(n):
            if i != j:
                B[i, j] = A[i, j] / A[i, i]

    print("\n--- Матрица B (итерационная матрица) ---")
    print(np.array2string(B, precision=6, suppress_small=True))

    # Operator norm of B (infinity norm = max row sum of abs values)
    norm_B_inf = np.max(np.sum(np.abs(B), axis=1))
    norm_B_1   = np.max(np.sum(np.abs(B), axis=0))
    norm_B_2   = np.linalg.norm(B, 2)

    print(f"\nКоэффициент сжатия (норма B):")
    print(f"  ||B||_inf (строчная) = {norm_B_inf:.6f}")
    print(f"  ||B||_1   (столбцовая) = {norm_B_1:.6f}")
    print(f"  ||B||_2   (спектральная) = {norm_B_2:.6f}")

    # Choose the smallest norm for convergence check
    norm_B = min(norm_B_inf, norm_B_1, norm_B_2)

    if norm_B >= 1:
        print(f"\n⚠️  Предупреждение: ||B|| = {norm_B:.6f} >= 1, метод может не сходиться!")
    else:
        print(f"\n✓ Метод сходится, используем ||B|| = {norm_B:.6f} < 1")

    # Iterative process starting from x0 = 0
    x = np.zeros(n)

    print(f"\n{'Итерация':<10} {'||x_new - x_old||':<22} Значения x")
    print("-" * 80)

    for k in range(max_iter):
        # x_new = y_tilde - B * x  (Jacobi: use old x for all components)
        x_new = y_tilde - B @ x

        diff = np.linalg.norm(x_new - x, np.inf)

        x_str = "  ".join([f"{v:10.6f}" for v in x_new])
        print(f"{k+1:<10} {diff:<22.8e} [{x_str}]")

        if diff < eps:
            print(f"\n✓ Сходимость достигнута за {k+1} итераций (||Δx|| = {diff:.2e} < ε = {eps:.2e})")
            return x_new, k+1, norm_B

        x = x_new

    print(f"\n⚠️  Достигнуто максимальное число итераций ({max_iter})")
    return x, max_iter, norm_B


def verify_solution(A, y, x):
    """Verify solution by computing residual"""
    residual = A @ x - y
    print(f"\nПроверка (A*x - y): [{', '.join([f'{r:.2e}' for r in residual])}]")
    print(f"||A*x - y||_inf = {np.max(np.abs(residual)):.2e}")


def main():
    filename = 'input22.txt'

    # Ask user for epsilon
    try:
        eps = float(input("Введите точность ε (например, 1e-6): "))
    except ValueError:
        eps = 1e-6
        print(f"Неверный ввод, используется ε = {eps}")

    print(f"\nТочность: ε = {eps}\n")

    systems = parse_input(filename)

    for idx, (A, y) in enumerate(systems, 1):
        print("\n" + "=" * 80)
        print(f"СИСТЕМА {idx}: A{idx} * x = y{idx}")
        print("=" * 80)

        n = len(y)
        print(f"\nРазмерность системы: {n}x{n}")

        print("\nМатрица A:")
        print(np.array2string(A, precision=4, suppress_small=True))

        print("\nВектор y:")
        print(y)

        try:
            x_sol, iters, norm_B = iterative_solve(A, y, eps=eps)

            print(f"\n{'='*40}")
            print(f"ОТВЕТ (система {idx}):")
            for i, xi in enumerate(x_sol):
                print(f"  x{i+1} = {xi:.8f}")

            verify_solution(A, y, x_sol)

            # Cross-check with numpy
            x_exact = np.linalg.solve(A, y)
            print(f"\nТочное решение (numpy.linalg.solve):")
            for i, xi in enumerate(x_exact):
                print(f"  x{i+1} = {xi:.8f}")

            print(f"\nПогрешность итерационного решения: {np.max(np.abs(x_sol - x_exact)):.2e}")

        except ValueError as e:
            print(f"\nОшибка: {e}")


if __name__ == "__main__":
    main()