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

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

double complex 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 */
void lu_decomposition(double complex B[N][N],
                      double complex L[N][N],
                      double complex U[N][N])
{
    for(int i=0;i<N;i++){
        for(int j=0;j<N;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++){

        for(int j=i;j<N;j++){
            double complex sum = 0.0 + 0.0*I;

            for(int k=0;k<i;k++)
                sum += L[i][k]*U[k][j];

            U[i][j] = B[i][j] - sum;
        }

        for(int j=i+1;j<N;j++){
            double complex sum = 0.0 + 0.0*I;

            for(int k=0;k<i;k++)
                sum += L[j][k]*U[k][i];

            L[j][i] = (B[j][i]-sum)/U[i][i];
        }
    }
}

/* Forward substitution */
void forward_substitution(double complex L[N][N],
                          double complex x[N],
                          double complex y[N])
{
    for(int i=0;i<N;i++){

        double complex sum = 0.0 + 0.0*I;

        for(int j=0;j<i;j++)
            sum += L[i][j]*y[j];

        y[i] = x[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 + 0.0*I;

        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 complex lambda_hat)
{
    double complex B[N][N];
    double complex L[N][N];
    double complex U[N][N];

    double complex x[N] =
    {
        1.0+0.0*I,
        1.0+0.0*I,
        1.0+0.0*I,
        1.0+0.0*I
    };

    double complex y[N];
    double complex x_next[N];

    for(int i=0;i<N;i++){
        for(int j=0;j<N;j++){

            B[i][j] = A[i][j];

            if(i==j)
                B[i][j] -= lambda_hat;
        }
    }

    lu_decomposition(B,L,U);

    printf("\n=====================================\n");
    printf("shift = %.6f %+ .6fi\n",
           creal(lambda_hat), cimag(lambda_hat));
    printf("=====================================\n");

    int iter=0;
    double diff;
    double complex lambda = 0.0 + 0.0*I;

    do{

        forward_substitution(L,x,y);
        backward_substitution(U,y,x_next);

        /* largest magnitude component */
        double max_mag = 0.0;
        int idx = 0;

        for(int i=0;i<N;i++){
            if(cabs(x_next[i]) > max_mag){
                max_mag = cabs(x_next[i]);
                idx = i;
            }
        }

        double complex scale = x_next[idx];

        for(int i=0;i<N;i++)
            x_next[i] /= scale;

        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/scale;

        printf("Iter %2d : lambda = %.8f %+ .8fi\n",
               iter+1,
               creal(lambda),
               cimag(lambda));

        iter++;

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

    printf("\nFinal eigenvalue:\n");
    printf("%.10f %+ .10fi\n",
           creal(lambda),
           cimag(lambda));

    printf("\nEigenvector:\n");

    for(int i=0;i<N;i++){
        printf("(%.6f %+ .6fi)\n",
               creal(x[i]),
               cimag(x[i]));
    }
}

int main()
{
    inverse_iteration(-5.4 + 0.0*I);

    inverse_iteration(-2.1 + 1.7*I);

    inverse_iteration(-2.1 - 1.7*I);

    return 0;
}