#include <stdio.h>
#include <math.h>
#include <complex.h>
 
#define N 4
#define EPS 1e-6
#define MAX_ITER 100
 
// MATRIX 1 definition (stored as complex numbers)
double complex A[N][N] = { 
    {1.0 + 0.0*I, 2.0 + 0.0*I, 3.0 + 0.0*I, 5.0 + 0.0*I}, 
    {3.0 + 0.0*I, -5.0 + 0.0*I, 1.0 + 0.0*I, 4.0 + 0.0*I}, 
    {5.0 + 0.0*I, 9.0 + 0.0*I, 2.0 + 0.0*I, -6.0 + 0.0*I}, 
    {1.0 + 0.0*I, 7.0 + 0.0*I, 4.0 + 0.0*I, 1.0 + 0.0*I} 
};
 
// Complex LU Decomposition with Partial Pivoting
void lu_decomposition_complex(double complex B[N][N], double complex L[N][N], double complex U[N][N], int P[N]) {
    double complex 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*I : 0.0 + 0.0*I;
            U[i][j] = 0.0 + 0.0*I;
        }
    }
 
    for (int i = 0; i < N; i++) {
        int pivot_row = i;
        double max_val = cabs(A_temp[i][i]); // cabs computes the magnitude of complex number
        for (int r = i + 1; r < N; r++) {
            if (cabs(A_temp[r][i]) > max_val) {
                max_val = cabs(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 complex 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 complex 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 complex 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];
        }
    }
}
 
// Forward Substitution
void forward_substitution(double complex L[N][N], double complex x[N], double complex y[N], int P[N]) {
    double complex px[N];
    for (int i = 0; i < N; i++) px[i] = x[P[i]];
    for (int i = 0; i < N; i++) {
        double complex sum = 0.0;
        for (int j = 0; j < i; j++) sum += L[i][j] * y[j];
        y[i] = px[i] - sum;
    }
}
 
// Backward Substitution
void backward_substitution(double complex U[N][N], double complex y[N], double complex x[N]) {
    for (int i = N - 1; i >= 0; i--) {
        double complex 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];
    }
}
 
// Complex Inverse Iteration Routine
void find_eigen_complex(double complex lambda_hat) {
    double complex B[N][N], L[N][N], U[N][N];
    double complex y[N], x_next[N];
    int P[N];
 
    // Initial guess vector x^(0) = [1, 1, 1, 1]^T
    double complex x[N] = {1.0 + 0.0*I, 1.0 + 0.0*I, 1.0 + 0.0*I, 1.0 + 0.0*I};
 
    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_complex(B, L, U, P);
 
    printf("\n==================================================================\n");
    printf(" Inverse Iteration with Complex Shift: lambda_hat = %.2f + %.2fi\n", creal(lambda_hat), cimag(lambda_hat));
    printf("==================================================================\n");
 
    int iter = 0;
    double diff;
    double complex lambda = 0.0;
 
    do {
        forward_substitution(L, x, y, P);
        backward_substitution(U, y, x_next);
 
        // Scaling factor based on max absolute value (magnitude)
        double complex max_comp = x_next[0];
        double max_mag = cabs(x_next[0]);
        for (int i = 1; i < N; i++) {
            if (cabs(x_next[i]) > max_mag) {
                max_mag = cabs(x_next[i]);
                max_comp = x_next[i];
            }
        }
 
        // Normalize vector entries
        for (int i = 0; i < N; i++) {
            x_next[i] /= max_comp;
        }
 
        // Compute vector diff convergence
        diff = 0.0;
        for (int i = 0; i < N; i++) {
            diff += cabs(x_next[i] - x[i]);
        }
 
        for (int i = 0; i < N; i++) {
            x[i] = x_next[i];
        }
 
        lambda = lambda_hat + (1.0 / max_comp);
 
        printf("Iter %2d: lambda = %7.4f + %7.4fi\n", iter + 1, creal(lambda), cimag(lambda));
 
        iter++;
    } while (diff > EPS && iter < MAX_ITER);
 
    printf("\n CONVERGED SOLUTION:\n");
    printf(" Eigenvalue (lambda) = %9.4f + %9.4fi\n", creal(lambda), cimag(lambda));
    printf(" Eigenvector (x)    = [%.4f+%.4fi, %.4f+%.4fi, %.4f+%.4fi, %.4f+%.4fi]\n", 
           creal(x[0]), cimag(x[0]), creal(x[1]), cimag(x[1]), 
           creal(x[2]), cimag(x[2]), creal(x[3]), cimag(x[3]));
    printf("==================================================================\n");
}
 
int main() {
    printf("=== Complex Inverse Iteration Analysis for Matrix 1 ===\n");
 
    // 1. Finding the known real negative root
    find_eigen_complex(-5.0 + 0.0*I); 
 
    // 2. Finding the complex conjugate pair by supplying complex guesses near the spectrum floor
    find_eigen_complex(-2.0 + 1.5*I); // Targets upper complex plane root
    find_eigen_complex(-2.0 - 1.5*I); // Targets lower complex plane root
 
    return 0;
}

I2luY2x1ZGUgPHN0ZGlvLmg+CiNpbmNsdWRlIDxtYXRoLmg+CiNpbmNsdWRlIDxjb21wbGV4Lmg+CgojZGVmaW5lIE4gNAojZGVmaW5lIEVQUyAxZS02CiNkZWZpbmUgTUFYX0lURVIgMTAwCgovLyBNQVRSSVggMSBkZWZpbml0aW9uIChzdG9yZWQgYXMgY29tcGxleCBudW1iZXJzKQpkb3VibGUgY29tcGxleCBBW05dW05dID0geyAKICAgIHsxLjAgKyAwLjAqSSwgMi4wICsgMC4wKkksIDMuMCArIDAuMCpJLCA1LjAgKyAwLjAqSX0sIAogICAgezMuMCArIDAuMCpJLCAtNS4wICsgMC4wKkksIDEuMCArIDAuMCpJLCA0LjAgKyAwLjAqSX0sIAogICAgezUuMCArIDAuMCpJLCA5LjAgKyAwLjAqSSwgMi4wICsgMC4wKkksIC02LjAgKyAwLjAqSX0sIAogICAgezEuMCArIDAuMCpJLCA3LjAgKyAwLjAqSSwgNC4wICsgMC4wKkksIDEuMCArIDAuMCpJfSAKfTsKCi8vIENvbXBsZXggTFUgRGVjb21wb3NpdGlvbiB3aXRoIFBhcnRpYWwgUGl2b3RpbmcKdm9pZCBsdV9kZWNvbXBvc2l0aW9uX2NvbXBsZXgoZG91YmxlIGNvbXBsZXggQltOXVtOXSwgZG91YmxlIGNvbXBsZXggTFtOXVtOXSwgZG91YmxlIGNvbXBsZXggVVtOXVtOXSwgaW50IFBbTl0pIHsKICAgIGRvdWJsZSBjb21wbGV4IEFfdGVtcFtOXVtOXTsKICAgIGZvciAoaW50IGkgPSAwOyBpIDwgTjsgaSsrKSB7CiAgICAgICAgUFtpXSA9IGk7CiAgICAgICAgZm9yIChpbnQgaiA9IDA7IGogPCBOOyBqKyspIHsKICAgICAgICAgICAgQV90ZW1wW2ldW2pdID0gQltpXVtqXTsKICAgICAgICAgICAgTFtpXVtqXSA9IChpID09IGopID8gMS4wICsgMC4wKkkgOiAwLjAgKyAwLjAqSTsKICAgICAgICAgICAgVVtpXVtqXSA9IDAuMCArIDAuMCpJOwogICAgICAgIH0KICAgIH0KCiAgICBmb3IgKGludCBpID0gMDsgaSA8IE47IGkrKykgewogICAgICAgIGludCBwaXZvdF9yb3cgPSBpOwogICAgICAgIGRvdWJsZSBtYXhfdmFsID0gY2FicyhBX3RlbXBbaV1baV0pOyAvLyBjYWJzIGNvbXB1dGVzIHRoZSBtYWduaXR1ZGUgb2YgY29tcGxleCBudW1iZXIKICAgICAgICBmb3IgKGludCByID0gaSArIDE7IHIgPCBOOyByKyspIHsKICAgICAgICAgICAgaWYgKGNhYnMoQV90ZW1wW3JdW2ldKSA+IG1heF92YWwpIHsKICAgICAgICAgICAgICAgIG1heF92YWwgPSBjYWJzKEFfdGVtcFtyXVtpXSk7CiAgICAgICAgICAgICAgICBwaXZvdF9yb3cgPSByOwogICAgICAgICAgICB9CiAgICAgICAgfQoKICAgICAgICBpZiAocGl2b3Rfcm93ICE9IGkpIHsKICAgICAgICAgICAgaW50IHQgPSBQW2ldOyBQW2ldID0gUFtwaXZvdF9yb3ddOyBQW3Bpdm90X3Jvd10gPSB0OwogICAgICAgICAgICBmb3IgKGludCBrID0gMDsgayA8IE47IGsrKykgewogICAgICAgICAgICAgICAgZG91YmxlIGNvbXBsZXggdG1wID0gQV90ZW1wW2ldW2tdOwogICAgICAgICAgICAgICAgQV90ZW1wW2ldW2tdID0gQV90ZW1wW3Bpdm90X3Jvd11ba107CiAgICAgICAgICAgICAgICBBX3RlbXBbcGl2b3Rfcm93XVtrXSA9IHRtcDsKICAgICAgICAgICAgfQogICAgICAgIH0KCiAgICAgICAgZm9yIChpbnQgaiA9IGk7IGogPCBOOyBqKyspIHsKICAgICAgICAgICAgZG91YmxlIGNvbXBsZXggc3VtID0gMC4wOwogICAgICAgICAgICBmb3IgKGludCBrID0gMDsgayA8IGk7IGsrKykgc3VtICs9IExbaV1ba10gKiBVW2tdW2pdOwogICAgICAgICAgICBVW2ldW2pdID0gQV90ZW1wW2ldW2pdIC0gc3VtOwogICAgICAgIH0KICAgICAgICBmb3IgKGludCBqID0gaSArIDE7IGogPCBOOyBqKyspIHsKICAgICAgICAgICAgZG91YmxlIGNvbXBsZXggc3VtID0gMC4wOwogICAgICAgICAgICBmb3IgKGludCBrID0gMDsgayA8IGk7IGsrKykgc3VtICs9IExbal1ba10gKiBVW2tdW2ldOwogICAgICAgICAgICBMW2pdW2ldID0gKEFfdGVtcFtqXVtpXSAtIHN1bSkgLyBVW2ldW2ldOwogICAgICAgIH0KICAgIH0KfQoKLy8gRm9yd2FyZCBTdWJzdGl0dXRpb24Kdm9pZCBmb3J3YXJkX3N1YnN0aXR1dGlvbihkb3VibGUgY29tcGxleCBMW05dW05dLCBkb3VibGUgY29tcGxleCB4W05dLCBkb3VibGUgY29tcGxleCB5W05dLCBpbnQgUFtOXSkgewogICAgZG91YmxlIGNvbXBsZXggcHhbTl07CiAgICBmb3IgKGludCBpID0gMDsgaSA8IE47IGkrKykgcHhbaV0gPSB4W1BbaV1dOwogICAgZm9yIChpbnQgaSA9IDA7IGkgPCBOOyBpKyspIHsKICAgICAgICBkb3VibGUgY29tcGxleCBzdW0gPSAwLjA7CiAgICAgICAgZm9yIChpbnQgaiA9IDA7IGogPCBpOyBqKyspIHN1bSArPSBMW2ldW2pdICogeVtqXTsKICAgICAgICB5W2ldID0gcHhbaV0gLSBzdW07CiAgICB9Cn0KCi8vIEJhY2t3YXJkIFN1YnN0aXR1dGlvbgp2b2lkIGJhY2t3YXJkX3N1YnN0aXR1dGlvbihkb3VibGUgY29tcGxleCBVW05dW05dLCBkb3VibGUgY29tcGxleCB5W05dLCBkb3VibGUgY29tcGxleCB4W05dKSB7CiAgICBmb3IgKGludCBpID0gTiAtIDE7IGkgPj0gMDsgaS0tKSB7CiAgICAgICAgZG91YmxlIGNvbXBsZXggc3VtID0gMC4wOwogICAgICAgIGZvciAoaW50IGogPSBpICsgMTsgaiA8IE47IGorKykgc3VtICs9IFVbaV1bal0gKiB4W2pdOwogICAgICAgIHhbaV0gPSAoeVtpXSAtIHN1bSkgLyBVW2ldW2ldOwogICAgfQp9CgovLyBDb21wbGV4IEludmVyc2UgSXRlcmF0aW9uIFJvdXRpbmUKdm9pZCBmaW5kX2VpZ2VuX2NvbXBsZXgoZG91YmxlIGNvbXBsZXggbGFtYmRhX2hhdCkgewogICAgZG91YmxlIGNvbXBsZXggQltOXVtOXSwgTFtOXVtOXSwgVVtOXVtOXTsKICAgIGRvdWJsZSBjb21wbGV4IHlbTl0sIHhfbmV4dFtOXTsKICAgIGludCBQW05dOwogICAgCiAgICAvLyBJbml0aWFsIGd1ZXNzIHZlY3RvciB4XigwKSA9IFsxLCAxLCAxLCAxXV5UCiAgICBkb3VibGUgY29tcGxleCB4W05dID0gezEuMCArIDAuMCpJLCAxLjAgKyAwLjAqSSwgMS4wICsgMC4wKkksIDEuMCArIDAuMCpJfTsKCiAgICBmb3IgKGludCBpID0gMDsgaSA8IE47IGkrKykgewogICAgICAgIGZvciAoaW50IGogPSAwOyBqIDwgTjsgaisrKSB7CiAgICAgICAgICAgIEJbaV1bal0gPSBBW2ldW2pdIC0gKGkgPT0gaiA/IGxhbWJkYV9oYXQgOiAwLjApOwogICAgICAgIH0KICAgIH0KCiAgICBsdV9kZWNvbXBvc2l0aW9uX2NvbXBsZXgoQiwgTCwgVSwgUCk7CgogICAgcHJpbnRmKCJcbj09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PVxuIik7CiAgICBwcmludGYoIiBJbnZlcnNlIEl0ZXJhdGlvbiB3aXRoIENvbXBsZXggU2hpZnQ6IGxhbWJkYV9oYXQgPSAlLjJmICsgJS4yZmlcbiIsIGNyZWFsKGxhbWJkYV9oYXQpLCBjaW1hZyhsYW1iZGFfaGF0KSk7CiAgICBwcmludGYoIj09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PVxuIik7CiAgICAKICAgIGludCBpdGVyID0gMDsKICAgIGRvdWJsZSBkaWZmOwogICAgZG91YmxlIGNvbXBsZXggbGFtYmRhID0gMC4wOwoKICAgIGRvIHsKICAgICAgICBmb3J3YXJkX3N1YnN0aXR1dGlvbihMLCB4LCB5LCBQKTsKICAgICAgICBiYWNrd2FyZF9zdWJzdGl0dXRpb24oVSwgeSwgeF9uZXh0KTsKCiAgICAgICAgLy8gU2NhbGluZyBmYWN0b3IgYmFzZWQgb24gbWF4IGFic29sdXRlIHZhbHVlIChtYWduaXR1ZGUpCiAgICAgICAgZG91YmxlIGNvbXBsZXggbWF4X2NvbXAgPSB4X25leHRbMF07CiAgICAgICAgZG91YmxlIG1heF9tYWcgPSBjYWJzKHhfbmV4dFswXSk7CiAgICAgICAgZm9yIChpbnQgaSA9IDE7IGkgPCBOOyBpKyspIHsKICAgICAgICAgICAgaWYgKGNhYnMoeF9uZXh0W2ldKSA+IG1heF9tYWcpIHsKICAgICAgICAgICAgICAgIG1heF9tYWcgPSBjYWJzKHhfbmV4dFtpXSk7CiAgICAgICAgICAgICAgICBtYXhfY29tcCA9IHhfbmV4dFtpXTsKICAgICAgICAgICAgfQogICAgICAgIH0KCiAgICAgICAgLy8gTm9ybWFsaXplIHZlY3RvciBlbnRyaWVzCiAgICAgICAgZm9yIChpbnQgaSA9IDA7IGkgPCBOOyBpKyspIHsKICAgICAgICAgICAgeF9uZXh0W2ldIC89IG1heF9jb21wOwogICAgICAgIH0KCiAgICAgICAgLy8gQ29tcHV0ZSB2ZWN0b3IgZGlmZiBjb252ZXJnZW5jZQogICAgICAgIGRpZmYgPSAwLjA7CiAgICAgICAgZm9yIChpbnQgaSA9IDA7IGkgPCBOOyBpKyspIHsKICAgICAgICAgICAgZGlmZiArPSBjYWJzKHhfbmV4dFtpXSAtIHhbaV0pOwogICAgICAgIH0KCiAgICAgICAgZm9yIChpbnQgaSA9IDA7IGkgPCBOOyBpKyspIHsKICAgICAgICAgICAgeFtpXSA9IHhfbmV4dFtpXTsKICAgICAgICB9CgogICAgICAgIGxhbWJkYSA9IGxhbWJkYV9oYXQgKyAoMS4wIC8gbWF4X2NvbXApOwogICAgICAgIAogICAgICAgIHByaW50ZigiSXRlciAlMmQ6IGxhbWJkYSA9ICU3LjRmICsgJTcuNGZpXG4iLCBpdGVyICsgMSwgY3JlYWwobGFtYmRhKSwgY2ltYWcobGFtYmRhKSk7CgogICAgICAgIGl0ZXIrKzsKICAgIH0gd2hpbGUgKGRpZmYgPiBFUFMgJiYgaXRlciA8IE1BWF9JVEVSKTsKCiAgICBwcmludGYoIlxuIENPTlZFUkdFRCBTT0xVVElPTjpcbiIpOwogICAgcHJpbnRmKCIgRWlnZW52YWx1ZSAobGFtYmRhKSA9ICU5LjRmICsgJTkuNGZpXG4iLCBjcmVhbChsYW1iZGEpLCBjaW1hZyhsYW1iZGEpKTsKICAgIHByaW50ZigiIEVpZ2VudmVjdG9yICh4KSAgICA9IFslLjRmKyUuNGZpLCAlLjRmKyUuNGZpLCAlLjRmKyUuNGZpLCAlLjRmKyUuNGZpXVxuIiwgCiAgICAgICAgICAgY3JlYWwoeFswXSksIGNpbWFnKHhbMF0pLCBjcmVhbCh4WzFdKSwgY2ltYWcoeFsxXSksIAogICAgICAgICAgIGNyZWFsKHhbMl0pLCBjaW1hZyh4WzJdKSwgY3JlYWwoeFszXSksIGNpbWFnKHhbM10pKTsKICAgIHByaW50ZigiPT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09XG4iKTsKfQoKaW50IG1haW4oKSB7CiAgICBwcmludGYoIj09PSBDb21wbGV4IEludmVyc2UgSXRlcmF0aW9uIEFuYWx5c2lzIGZvciBNYXRyaXggMSA9PT1cbiIpOwogICAgCiAgICAvLyAxLiBGaW5kaW5nIHRoZSBrbm93biByZWFsIG5lZ2F0aXZlIHJvb3QKICAgIGZpbmRfZWlnZW5fY29tcGxleCgtNS4wICsgMC4wKkkpOyAKICAgIAogICAgLy8gMi4gRmluZGluZyB0aGUgY29tcGxleCBjb25qdWdhdGUgcGFpciBieSBzdXBwbHlpbmcgY29tcGxleCBndWVzc2VzIG5lYXIgdGhlIHNwZWN0cnVtIGZsb29yCiAgICBmaW5kX2VpZ2VuX2NvbXBsZXgoLTIuMCArIDEuNSpJKTsgLy8gVGFyZ2V0cyB1cHBlciBjb21wbGV4IHBsYW5lIHJvb3QKICAgIGZpbmRfZWlnZW5fY29tcGxleCgtMi4wIC0gMS41KkkpOyAvLyBUYXJnZXRzIGxvd2VyIGNvbXBsZXggcGxhbmUgcm9vdAogICAgCiAgICByZXR1cm4gMDsKfQ==