/* minversion v1.0.0 Ver.: OK (c) 02/10 - Christos Iosifides */
/*
   This program inverts a nxn matrix with Gaussian elimination

   Name:			minversion
   Version:		    1.0.0
   Copyright:	    Ch Iossif @ 2010
   Author:		    Christos Iosifidis
   Date:			05/02/10 09:51
   Modified:		05/02/10 11:22
   Description:     This program inverts a nxn matrix with Gaussian elimination
   Compile:
            Linux:     gcc lsqr.c -o lsqr -lm
            Windows:   gcc lsqr.c -o lsqr.exe -lm or similar with any ANSI C compiler

  INPUT:

  N =  3

  A =

  1 3 1
  1 1 2
  2 3 4

  OUTPUT:

  A^(-1) =

   2  9 -5
   0 -2  1
  -1 -3  2

  |A| = -1

 Copyright (C) Dec 2009 Ch Iossif <chiossif@yahoo.com>

 This program is free software: you can redistribute it and/or modify
 it under the terms of the GNU General Public License as published by
 the Free Software Foundation, either version 3 of the License, or
 (at your option) any later version.

 This program is distributed in the hope that it will be useful,
 but WITHOUT ANY WARRANTY; without even the implied warranty of
 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 GNU General Public License for more details.

 You should have received a copy of the GNU General Public License
 along with this program.  If not, see <http://www.gnu.org/licenses/>.
*/

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

#define myABS(x) (((x)<0)?-(x):(x))
#define mySQR(x) ((x)*(x))

//#define FP_TYPE double
//#define EPSILON 1.e-15
//#define FP_FORMAT "%16.10lf",

#define FP_TYPE long double
#define EPSILON 1.e-19
#define FP_FORMAT "%24.15Lf",

//#define FP_TYPE __float128
//#define EPSILON 1.e-34
//#define FP_FORMAT "%40.34Lf",(long double)

FP_TYPE mySQRT(FP_TYPE);
void gpl(void);
int M_inversion(FP_TYPE *,FP_TYPE *,FP_TYPE *,int);
int Matrix_inversionG(FP_TYPE *,FP_TYPE *,FP_TYPE *,int);
int Matrix_inversionC(FP_TYPE *,FP_TYPE *,FP_TYPE *,int);
int Gaussian_elimination(FP_TYPE *,FP_TYPE *,FP_TYPE *,int *,int);
int Cholesky_decomposition(FP_TYPE *, FP_TYPE *, int);

int main(void) {
    FP_TYPE *a, *a_inv, s;
    long double t;
    int n, *l;
    register int i, j, k;

    scanf("%d",&n);
    a=malloc(n*n*sizeof(FP_TYPE));
    l=malloc(n*sizeof(int));
    a_inv=malloc(n*n*sizeof(FP_TYPE));

    for (i=0; i<n; i++)
        for (j=0; j<n; j++) {
            scanf("%Lf",&t);
            a[i*n+j]=t;
        }

    printf("\nminversion inverts a nxn matrix with Gaussian elimination\n");
    printf("\nINPUT:\n");
    printf("\nN = %d\n",n);
    printf("\nA =\n");
    for (i=0; i<n; i++) {
        for (j=0; j<n; j++)
            printf(FP_FORMAT a[i*n+j]);
        printf("\n");
    }

    printf("\nMatrix_inversion with Gaussian elimination\n\nInput:\n\nA :\n");
    for (i=0; i<n; i++) {
        for (j=0; j<n; j++)
            printf(FP_FORMAT a[i*n+j]);
        printf("\n");
    }
    if (Matrix_inversionG(a,a_inv,&s,n)) {
        printf("\nOutput:\n\nA^(-1) :\n");
        for (i=0; i<n; i++) {
            for (j=0; j<n; j++)
                printf(FP_FORMAT a_inv[i*n+j]);
            printf("\n");
        }
    } else {
        printf("\nMatrix A is NOT invertible.\n");
    }

    printf("\nMatrix_inversion with Cholenksy decomposition\n\nInput:\n\nA :\n");
    for (i=0; i<n; i++) {
        for (j=0; j<n; j++)
            printf(FP_FORMAT a[i*n+j]);
        printf("\n");
    }
    if (Matrix_inversionC(a,a_inv,&s,n)) {
        printf("\nOutput:\n\nA^(-1) :\n");
        for (i=0; i<n; i++) {
            for (j=0; j<n; j++)
                printf(FP_FORMAT a_inv[i*n+j]);
            printf("\n");
        }
    } else {
        printf("\nMatrix A is NOT a symmetric, positive-definite matrix.\n");
    }

    printf("\nM_inversion with Gaussian elimination\n\nInput:\n\nA :\n");
    for (i=0; i<n; i++) {
        for (j=0; j<n; j++)
            printf(FP_FORMAT a[i*n+j]);
        printf("\n");
    }
    if (M_inversion(a,a_inv,&s,n)) {
        printf("\nOutput:\n\nA^(-1) :\n");
        for (i=0; i<n; i++) {
            for (j=0; j<n; j++)
                printf(FP_FORMAT a_inv[i*n+j]);
            printf("\n");
        }
    } else {
        printf("\nMatrix A is NOT invertible.\n");
    }

    printf("\nOUTPUT:\n");
    printf("\nA^(-1) =\n");
    for (i=0; i<n; i++) {
        for (j=0; j<n; j++)
            printf(FP_FORMAT a_inv[i*n+j]);
        printf("\n");
    }
    printf("\nDeterminant =");
    printf(FP_FORMAT s);
    printf("\n");

    printf("\nminversion v1.0.0 Ver.: OK (c) 02/10 - Christos Iosifides\n");
    gpl();
    return 0;
}

int M_inversion(FP_TYPE *a,FP_TYPE *a_inv,FP_TYPE *det,int n) {
    FP_TYPE *b, d, s;
    int *l, r_value, n_1, m;
    register int i, j, k;

    l=malloc(n*sizeof(int));
    r_value=1;
    d=1.0;
    n_1=n-1;
    for (k=0; k<n_1; k++) {
        m=k;
        for (i=k+1; i<n; i++)
            if (myABS(a[i*n+k])>myABS(a[m*n+k]))
                m=i;
        l[k]=m;
        d*=a[m*n+k];
        if (myABS(a[m*n+k])<EPSILON)
            return 0;
        if (m!=k) {
            d=-d;
            for (j=k; j<n; j++) {
                s=a[m*n+j];
                a[m*n+j]=a[k*n+j];
                a[k*n+j]=s;
            }
        }
        for (i=k+1; i<n; i++) {
            s=a[i*n+k]/a[k*n+k];
            a[i*n+k]=s;
            for (j=k+1; j<n; j++)
                a[i*n+j]-=s*a[k*n+j];
        }
    }
    d*=a[n_1*n+n_1];
    if (myABS(d)<EPSILON)
        return 0;
    *det=d;

    b=malloc(n*sizeof(FP_TYPE));
    for (i=0; i<n; i++) {
        for (j=0; j<n; j++) {
            b[j]=0.0;
        }
        b[i]=1.0;

        for (k=0; k<n_1; k++) {
            m=l[k];
            s=b[m];
            b[m]=b[k];
            b[k]=s;
            for (j=k+1; j<n; j++)
                b[j]-=s*a[j*n+k];
        }
        b[n_1]/=a[n_1*n+n_1];
        for (k=n-2; k>=0; k--) {
            s=0.;
            for (j=k+1; j<n; j++)
                s+=a[k*n+j]*b[j];
            b[k]=(b[k]-s)/a[k*n+k];
        }
        for (j=0; j<n; j++)
            a_inv[j*n+i]=b[j];
    }

    return r_value;
}

int Matrix_inversionG(FP_TYPE *ata,FP_TYPE *ata_inv,FP_TYPE *det,int n) {
    FP_TYPE *tmpa, *b;
    int *l;
    int r_value;
    register int i, j, k;

    tmpa=malloc(n*n*sizeof(FP_TYPE));
    b=malloc(n*sizeof(FP_TYPE));
    l=malloc(n*sizeof(int));
    r_value=1;
    for (i=0; i<n; i++) {
        for (j=0; j<n; j++) {
            for (k=0; k<n; k++)
                tmpa[j*n+k]=ata[j*n+k];
            b[j]=0.0;
        }
        b[i]=1.0;
        r_value*=Gaussian_elimination(tmpa,b,det,l,n);
//        r_value*=Cholesky_decomposition(tmpa, b, n);
        for (j=0; j<n; j++)
            ata_inv[j*n+i]=b[j];
    }
    return r_value;
}

int Matrix_inversionC(FP_TYPE *ata,FP_TYPE *ata_inv,FP_TYPE *det,int n) {
    FP_TYPE *tmpa, *b;
    int *l;
    int r_value;
    register int i, j, k;

    tmpa=malloc(n*n*sizeof(FP_TYPE));
    b=malloc(n*sizeof(FP_TYPE));
    l=malloc(n*sizeof(int));
    r_value=1;
    for (i=0; i<n; i++) {
        for (j=0; j<n; j++) {
            for (k=0; k<n; k++)
                tmpa[j*n+k]=ata[j*n+k];
            b[j]=0.0;
        }
        b[i]=1.0;
//        r_value*=Gaussian_elimination(tmpa,b,det,l,n);
        r_value*=Cholesky_decomposition(tmpa, b, n);
        for (j=0; j<n; j++)
            ata_inv[j*n+i]=b[j];
    }
    return r_value;
}

int Gaussian_elimination(FP_TYPE *a, FP_TYPE *b, FP_TYPE *det, int *l, int n) {
    FP_TYPE s, d;
    int n_1, m;
    register int i, j, k;

    d=1.0;
    n_1=n-1;
    for (k=0; k<n_1; k++) {
        m=k;
        for (i=k+1; i<n; i++)
            if (myABS(a[i*n+k])>myABS(a[m*n+k]))
                m=i;
        l[k]=m;
        d*=a[m*n+k];
        if (myABS(a[m*n+k])<EPSILON)
            return 0;
        if (m!=k) {
            d=-d;
            for (j=k; j<n; j++) {
                s=a[m*n+j];
                a[m*n+j]=a[k*n+j];
                a[k*n+j]=s;
            }
        }
        for (i=k+1; i<n; i++) {
            s=a[i*n+k]/a[k*n+k];
            a[i*n+k]=s;
            for (j=k+1; j<n; j++)
                a[i*n+j]-=s*a[k*n+j];
        }
    }
    d*=a[n_1*n+n_1];
    if (myABS(d)<EPSILON)
        return 0;
    *det=d;
    for (k=0; k<n_1; k++) {
        m=l[k];
        s=b[m];
        b[m]=b[k];
        b[k]=s;
        for (i=k+1; i<n; i++)
            b[i]-=s*a[i*n+k];
    }
    b[n_1]/=a[n_1*n+n_1];
    for (k=n-2; k>=0; k--) {
        s=0.;
        for (j=k+1; j<n; j++)
            s+=a[k*n+j]*b[j];
        b[k]=(b[k]-s)/a[k*n+k];
    }
    return 1;
}

int Cholesky_decomposition(FP_TYPE *a, FP_TYPE *b, int n) {
    FP_TYPE s;
    register int i, j, k;

    for (i=0; i<n; i++) {
        s=0.0;
        if (i)
            for (j=0; j<i; j++)
                s+=mySQR(a[i*n+j]);
        if ((a[i*n+i]-s)<EPSILON)
            return 0;
        a[i*n+i]=mySQRT(a[i*n+i]-s);
        if (i<n-1)
            for (j=i+1; j<n; j++) {
                s=0.0;
                if (i)
                    for (k=0; k<i; k++)
                        s+=a[i*n+k]*a[j*n+k];
                a[j*n+i]=(a[j*n+i]-s)/a[i*n+i];
            }
    }
    b[0]/=a[0*n+0];
    for (k=1; k<n; k++) {
        s=0.0;
        for (j=0; j<k; j++)
            s+=a[k*n+j]*b[j];
        b[k]=(b[k]-s)/a[k*n+k];
    }
    b[n-1]/=a[(n-1)*n+n-1];
    for (i=0; i<n-1; i++) {
        s=0.0;
        for (j=n-i-1; j<n; j++)
            s+=a[j*n+n-i-2]*b[j];
        b[n-i-2]=(b[n-i-2]-s)/a[(n-i-2)*n+n-i-2];
    }
    return 1;
}

void gpl(void) {
    puts("\nCopyright (C) Dec 2009 Ch Iossif <chiossif@yahoo.com>");

    puts("\nThis program is free software: you can redistribute it and/or modify");
    puts("it under the terms of the GNU General Public License as published by");
    puts("the Free Software Foundation, either version 3 of the License, or");
    puts("(at your option) any later version.");

    puts("\nThis program is distributed in the hope that it will be useful,");
    puts("but WITHOUT ANY WARRANTY; without even the implied warranty of");
    puts("MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the");
    puts("GNU General Public License for more details.");

    puts("\nYou should have received a copy of the GNU General Public License");
    puts("along with this program.  If not, see <http://www.gnu.org/licenses/>.\n");
}

FP_TYPE mySQRT(FP_TYPE x) {
    FP_TYPE px,a;

    a=px=x;
    x/=2.0;
    while (myABS(x-px)>EPSILON) {
        px=x;
        x=(px+a/px)/2;
    }
    return x;
}