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

#define MAX 20


void            dataentry(void);
void            reserve(void);
void            initialize(void);
void            inputdim(void);
void            printresults(void);
void            erroranalysis(void);

void            deter(long *nd, double w[][MAX], double *ddet);

double          aa,
                ah,
                aq,
                temp;
double          a[MAX][MAX],
                b[MAX][MAX],
                x[MAX][MAX];

double          dd[MAX][MAX];
double          ac[MAX][MAX],
                det;

long
                ii,
                ik,
                imethod;
long
                jj,
                kk,
                n,
                nv;


char            filename[5][FILENAME_MAX];
FILE           *pFile[5];

//   Variables
//
//   n is the dimension of the input matrix.
//
//       a is the input matrix.
//
//       x is initialized as the identity matrix.
//    It is transformed to the inverse matrix.
//
//   b is the product of matrix a and the inverse matrix x.
//     It should equal the identity matrix.
//

int             main(void)
{
    long i;
    long j;
    long k;

    printf("\n\n Matrix Inversion via Gaussian Elimination with Partial Pivoting. \n");
    printf("\n Version 1.1,  October 27, 1999 \n");
    printf("\n by Tom Irvine ");
    printf("\n Email:  tomirvine@aol.com \n");

    inputdim();
    initialize();
    dataentry();
    reserve();

    deter(&n, a, &det);

    if (fabs(det) <= 1.0e-200) {
        printf("\n\n Error:  determinant of input matrix  = %12.5e \n", det);
        exit(1);
    } else {
        printf("\n\n Determinant of input matrix  = %12.5e \n", det);
    }

/****** begin forward elimination  *********/

    for (i = 0; i < n; ++i) {
        aq = fabs(a[i][i]);

        for (jj = i + 1; jj < n; ++jj) {
            ah = fabs(a[jj][i]);

            if (aq < ah) {
                aq = fabs(a[jj][i]);

                for (kk = 0; kk < n; ++kk) {
                    temp = a[jj][kk];
                    a[jj][kk] = a[i][kk];
                    a[i][kk] = temp;
                }
                for (nv = 0; nv < n; nv++) {
                    temp = x[jj][nv];
                    x[jj][nv] = x[i][nv];
                    x[i][nv] = temp;
                }
            }
        }
        ik = i + 1;
        for (ii = ik; ii < n; ++ii) {
            if (fabs(a[i][i]) < 1.0e-200) {
                printf("\n Error:  zero in diagonal.\n");
                exit(1);
            }
            aa = a[ii][i] / a[i][i];

            for (nv = 0; nv < n; nv++) {
                x[ii][nv] += -x[i][nv] * a[ii][i] / a[i][i];
            }
            for (j = i; j < n; ++j) {
                a[ii][j] += -aa * a[i][j];
                if (fabs(a[ii][j]) < 1.0e-10) {
                    a[ii][j] = 0.;
                }
            }
        }
    }


/********* back substitution ***********************************/

    for (nv = 0; nv < n; nv++) {
        x[n - 1][nv] /= a[n - 1][n - 1];

        for (i = n - 2; i >= 0; i = i - 1) {
            for (k = i + 1; k < n; ++k) {
                x[i][nv] += (-x[k][nv] * a[i][k]);
            }
            x[i][nv] /= a[i][i];
        }
    }

/************* end of gaussian elimination ***********************/

    printresults();

    deter(&n, x, &det);
    printf("\n\n Determinant of inverse matrix  = %12.5e \n", det);

    erroranalysis();

    deter(&n, b, &det);
    printf("\n Determinant of error matrix  = %12.5e \n\n", det);

    if (imethod == 1) {
        printf("\n Input data written to file:    inv_data.in ");
    }
    printf("\n Output data written to file:   inv_data.out");
    printf("\n Error matrix written to file:  inv_err.out \n");

/********************** end of main program **************************/
    return 0;
}
void            deter(long *nd, double w[][MAX], double *ddet)
{
/********* calculate determinant *********/

    long            i,
                    j,
                    k;
    double          a_a;

    for (i = 0; i < *nd; i++) {
        for (j = 0; j < *nd; j++) {
            dd[i][j] = w[i][j];
        }
    }

    *ddet = 1.;

    for (i = 0; i < *nd; i++) {
        for (k = i + 1; k < *nd; k++) {
            if (fabs(dd[i][i]) < 1.0e-50) {
                *ddet = 0;
                return;
            }
            a_a = dd[k][i] / dd[i][i];

            for (j = i; j < *nd; j++) {
                dd[k][j] -= a_a * dd[i][j];
            }
        }
        *ddet *= dd[i][i];
    }
}

void            dataentry(void)
{
    long i;
    long j;
    double          anum;

    if (imethod == 1) {
        strcpy(filename[0], "inv_data.in");
        pFile[0] = fopen(filename[0], "w");

        for (i = 0; i < n; i++) {
            for (j = 0; j < n; j++) {
                printf("\n Enter a[%ld][%ld] \n", i + 1, j + 1);
                scanf("%lf", &anum);

                a[i][j] = anum;

                if (j < (n - 1)) {
                    fprintf(pFile[0], "%12.5e \t", a[i][j]);
                } else {
                    fprintf(pFile[0], "%12.5e \n", a[i][j]);
                }
            }
        }
    }
    if (imethod == 2) {
        printf("\n Input filename.  Format is free, but no header lines allowed. \n");

        scanf("%s", filename[1]);
        pFile[1] = fopen(filename[1], "rb");

        if (pFile[1] == NULL) {
            printf("Input file error \n");
            exit(1);
        }
        for (i = 0; i < n; i++) {
            for (j = 0; j < n; j++) {
                fscanf(pFile[1], "%lf", &anum);
                a[i][j] = anum;
            }
        }

        if (n <= 6) {
            printf("\n\n  Input Matrix \n\n");

            for (i = 0; i < n; i++) {
                for (j = 0; j < n; j++) {
                    if (j < (n - 1)) {
                        printf("%12.5e \t", a[i][j]);
                    } else {
                        printf("%12.5e \n", a[i][j]);
                    }
                }
            }
        }
    }
    if (imethod != 1 && imethod != 2) {
        exit(1);
    }
}
void            reserve(void)
{
/****  copy data into reserve matrices  **/
    long i;
    long j;

    for (i = 0; i < n; ++i) {
        for (j = 0; j < n; ++j) {
            ac[i][j] = a[i][j];
        }
    }

}
void            initialize(void)
{
    long i;
    long j;
    for (i = 0; i < MAX; i++) {
        for (j = 0; j < MAX; j++) {
            a[i][j] = 0.;
            ac[i][j] = 0.;
            x[i][j] = 0.;
            b[i][j] = 0.;
        }
        x[i][i] = 1.;
    }

}
void            inputdim(void)
{
    printf("\n The input matrix has dimension n x n. \n\n Enter n \n");
    scanf("%ld", &n);

    if (n > 20) {
        printf("\n Error:  maximum dimension = 20 \n");
        exit(1);
    }
    if (n < 2) {
        printf("\n Error:  minimum dimension = 2 \n");
        exit(1);
    }
    printf("\n\n Select data entry method:  1=keyboard  2=input file \n");

    scanf("%ld", &imethod);
}
void            printresults(void)
{
    long i;
    long j;

    strcpy(filename[2], "inv_data.out");
    pFile[2] = fopen(filename[2], "w");

    if (n <= 6) {
        printf("\n\n Inverse Matrix \n\n");

        for (i = 0; i < n; i++) {
            for (j = 0; j < n; j++) {
                if (j < (n - 1)) {
                    printf("%12.5e \t", x[i][j]);
                } else {
                    printf("%12.5e \n", x[i][j]);
                }

            }
        }

    }
    for (i = 0; i < n; i++) {
        for (j = 0; j < n; j++) {
            if (j < (n - 1)) {
                fprintf(pFile[2], "%12.5e \t", x[i][j]);
            } else {
                fprintf(pFile[2], "%12.5e \n", x[i][j]);
            }

        }
    }


}
void            erroranalysis(void)
{
    long i;
    long j;
    long k;
    for (i = 0; i < n; i++) {
        for (j = 0; j < n; j++) {
            for (k = 0; k < n; k++) {
                b[i][j] += ac[i][k] * x[k][j];
            }
        }
    }


    strcpy(filename[3], "inv_err.out");
    pFile[3] = fopen(filename[3], "w");

    if (n <= 6) {
        printf("\n Error Analysis: \n Input Matrix x Inverse Matrix (should equal identity matrix). \n\n");

        for (i = 0; i < n; i++) {
            for (j = 0; j < n; j++) {
                if (j < (n - 1)) {
                    printf("%12.5e \t", b[i][j]);
                } else {
                    printf("%12.5e \n", b[i][j]);
                }

            }
        }

    }
    for (i = 0; i < n; i++) {
        for (j = 0; j < n; j++) {
            if (j < (n - 1)) {
                fprintf(pFile[3], "%12.5e \t", b[i][j]);
            } else {
                fprintf(pFile[3], "%12.5e \n", b[i][j]);
            }

        }
    }

}