#include <stdio.h>
#include <math.h>

#define N 4
#define EPS 1e-6
#define MAX_ITER 100

// MATRIX 1 definition
double A[N][N] = { 
    {1.0, 2.0, 3.0, 5.0}, 
    {3.0, -5.0, 1.0, 4.0}, 
    {5.0, 9.0, 2.0, -6.0}, 
    {1.0, 7.0, 4.0, 1.0} 
};

// LU Decomposition with Partial Pivoting to prevent numerical breakdown
void lu_decomposition(double B[N][N], double L[N][N], double U[N][N], int P[N]) {
    double A_temp[N][N];
    for (int i = 0; i < N; i++) {
        P[i] = i;
        for (int j = 0; j < N; j++) {
            A_temp[i][j] = B[i][j];
            L[i][j] = (i == j) ? 1.0 : 0.0;
            U[i][j] = 0.0;
        }
    }

    for (int i = 0; i < N; i++) {
        int pivot_row = i;
        double max_val = fabs(A_temp[i][i]);
        for (int r = i + 1; r < N; r++) {
            if (fabs(A_temp[r][i]) > max_val) {
                max_val = fabs(A_temp[r][i]);
                pivot_row = r;
            }
        }

        if (pivot_row != i) {
            int t = P[i]; P[i] = P[pivot_row]; P[pivot_row] = t;
            for (int k = 0; k < N; k++) {
                double tmp = A_temp[i][k];
                A_temp[i][k] = A_temp[pivot_row][k];
                A_temp[pivot_row][k] = tmp;
            }
        }

        for (int j = i; j < N; j++) {
            double sum = 0.0;
            for (int k = 0; k < i; k++) sum += L[i][k] * U[k][j];
            U[i][j] = A_temp[i][j] - sum;
        }
        for (int j = i + 1; j < N; j++) {
            double sum = 0.0;
            for (int k = 0; k < i; k++) sum += L[j][k] * U[k][i];
            L[j][i] = (A_temp[j][i] - sum) / U[i][i];
        }
    }
}

void forward_substitution(double L[N][N], double x[N], double y[N], int P[N]) {
    double px[N];
    for (int i = 0; i < N; i++) px[i] = x[P[i]];
    for (int i = 0; i < N; i++) {
        double sum = 0.0;
        for (int j = 0; j < i; j++) sum += L[i][j] * y[j];
        y[i] = px[i] - sum;
    }
}

void backward_substitution(double U[N][N], double y[N], double x[N]) {
    for (int i = N - 1; i >= 0; i--) {
        double sum = 0.0;
        for (int j = i + 1; j < N; j++) sum += U[i][j] * x[j];
        x[i] = (y[i] - sum) / U[i][i];
    }
}

void inverse_iteration(double lambda_hat) {
    double B[N][N], L[N][N], U[N][N];
    double y[N], x_next[N];
    int P[N];
    double x[N] = {1.0, 1.0, 1.0, 1.0}; // Initial vector

    for (int i = 0; i < N; i++) {
        for (int j = 0; j < N; j++) {
            B[i][j] = A[i][j] - (i == j ? lambda_hat : 0.0);
        }
    }

    lu_decomposition(B, L, U, P);

    printf("\n--- Shifting value (lambda_hat) = %.2f ---\n", lambda_hat);
    
    int iter = 0;
    double diff;
    double lambda = 0.0;

    do {
        forward_substitution(L, x, y, P);
        backward_substitution(U, y, x_next);

        // Find the element with maximum absolute value (keeping its sign)
        double max_val = x_next[0];
        for (int i = 1; i < N; i++) {
            if (fabs(x_next[i]) > fabs(max_val)) {
                max_val = x_next[i];
            }
        }

        // Normalize the vector (Slide 9/10 style)
        for (int i = 0; i < N; i++) {
            x_next[i] /= max_val;
        }

        // Check convergence
        diff = 0.0;
        for (int i = 0; i < N; i++) {
            diff += fabs(x_next[i] - x[i]);
        }

        for (int i = 0; i < N; i++) {
            x[i] = x_next[i];
        }

        // Slide 10 Formula: lambda = lambda_hat + 1 / max_val
        lambda = lambda_hat + (1.0 / max_val);
        iter++;
        
        printf("Iter %2d: Lambda = %10.6f\n", iter, lambda);

    } while (diff > EPS && iter < MAX_ITER);

    printf("Result -> Eigenvalue: %10.6f\n", lambda);
    printf("Result -> Eigenvector: [%.4f, %.4f, %.4f, %.4f]\n", x[0], x[1], x[2], x[3]);
}

int main() {
    printf("=== Matrix 1 Inverse Iteration Solutions ===\n");
    
    // Try running it with a couple of different real shifting guesses
    inverse_iteration(8.5);  // Targets the dominant real eigenvalue
    inverse_iteration(-6.0); // Targets the negative real eigenvalue
    
    return 0;
}