ইনভার্স পাওয়ার মেথড ও শিফটেড মেথড

import math

A = [
    [4.0, 1.0, 1.0],
    [1.0, 3.0, 1.0],
    [1.0, 1.0, 2.0],
]
n = 3

def matvec(M, v):
    return [sum(M[i][j] * v[j] for j in range(len(v))) for i in range(len(M))]

def norm(v):
    return math.sqrt(sum(x * x for x in v))

def gauss_solve(M, b):
    # M3/L10-এর গসিয়ান এলিমিনেশন, পার্শিয়াল পিভোটিংসহ, হাতে-লেখা -- কোনো লাইব্রেরি নেই
    n = len(M)
    aug = [row[:] + [b[i]] for i, row in enumerate(M)]
    for col in range(n):
        piv = max(range(col, n), key=lambda r: abs(aug[r][col]))
        aug[col], aug[piv] = aug[piv], aug[col]
        pivot_val = aug[col][col]
        for row in range(col + 1, n):
            factor = aug[row][col] / pivot_val
            for k in range(col, n + 1):
                aug[row][k] -= factor * aug[col][k]
    x = [0.0] * n
    for row in range(n - 1, -1, -1):
        s = aug[row][n] - sum(aug[row][k] * x[k] for k in range(row + 1, n))
        x[row] = s / aug[row][row]
    return x

def shifted(M, s):
    n = len(M)
    return [[M[i][j] - (s if i == j else 0.0) for j in range(n)] for i in range(n)]

def inverse_power(M, shift, v0, iters):
    Ms = shifted(M, shift)
    v = v0[:]
    prev_eig = None
    print(f"{'ইটার':>4} | {'eigenvalue আনুমান':>20} | {'পরিবর্তন':>14}")
    for i in range(1, iters + 1):
        w = gauss_solve(Ms, v)          # (A - shift*I) w = v  সমাধান -- এটাই A^-1 v
        nrm = norm(w)
        v = [x / nrm for x in w]
        Av = matvec(M, v)
        eig_est = sum(v[k] * Av[k] for k in range(len(v)))   # Rayleigh quotient, মূল A-তে
        if prev_eig is None:
            print(f"{i:4d} | {eig_est:20.10f} | {'--':>14}")
        else:
            print(f"{i:4d} | {eig_est:20.10f} | {abs(eig_est - prev_eig):14.10f}")
        prev_eig = eig_est
    return v, prev_eig

print("--- শিফট = ০ (সবচেয়ে ছোট আইগেনভ্যালু টার্গেট) ---")
v1, eig1 = inverse_power(A, 0.0, [1.0, 1.0, 1.0], 10)
print(f"চূড়ান্ত আইগেনভেক্টর: {[round(x, 6) for x in v1]}")

print("--- শিফট = ২.৫ (মধ্যম আইগেনভ্যালুর কাছাকাছি টার্গেট) ---")
v2, eig2 = inverse_power(A, 2.5, [1.0, 1.0, 1.0], 10)
print(f"চূড়ান্ত আইগেনভেক্টর: {[round(x, 6) for x in v2]}")