InverseMatrix

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

#define EPS 1.0e-12

void inverse(double **A,unsigned int n, unsigned int *error)
{
    unsigned int i,j,l,k;
    double   maxA,d,e;
    unsigned int *M;
    M=(unsigned int*)malloc(n*sizeof(unsigned int));
    (*error)=0;
    for(i=1; i<=n; i++)
    {
        // Czesciowy wybor elementu glownego wg wzoru (1.49)
        maxA=0.0;
        for(j=i; j<=n; j++)
        {
            d=A[j - 1][i - 1];
            if(fabs(maxA)<fabs(d))
            {
                maxA=d;
                k=j;
            }
        }
        if(fabs(maxA) < EPS)
        {
            (*error)=1;
            free(M);
            return;
        }
        /* Zapisywanie wskaznikow wierszy wystepowania elementu
          ekstremalnego w i-tej iteracji w postaci wektora M[i] */
        M[i - 1]=k;
        A[k - 1][i - 1]=1.0;
        for(j=1; j<=n; j++)
        {
            // Przestawienie i-tego wiersza z k-tym
            d=(double)(A[k - 1][j - 1]/maxA);
            A[k - 1][j - 1]=A[i - 1][j - 1];
            A[i - 1][j - 1]=d;
        }
        // Generacja ciagu macierzy (1.42) wg wzoru rekurencyjnego (1.41)
        for(j=1; j<=n; j++)
            if(j!=i)
            {
                d=A[j - 1][i - 1];
                A[j - 1][i - 1]=0.0;
                for(l=1; l<=n; l++)
                {
                    e=d*A[i - 1][l - 1];
                    A[j - 1][l - 1]-=e;
                }
            }
    }
    // Przestawianie kolumn macierzy zgodnie z wektorem wskaznikow M[i] (1.46)
    for(i=n; i>=1; i--)
    {
        k=M[i - 1];
        if(k!=i)
            for(j=1; j<=n; j++)
            {
                d=A[j - 1][i - 1];
                A[j - 1][i - 1]=A[j - 1][k - 1];
                A[j - 1][k - 1]=d;
            }
    }
    free(M);
}





int main(int argc, char **argv)
{
    FILE *in;
    FILE *out;
    unsigned int i,j,n;
    unsigned int  error;
    double **A;

    if ((in = fopen(argv[1], "rt"))
            == NULL)
    {
        fprintf(stderr, "Cannot open input file.\n");
        return 1;
    }

    if ((out = fopen(argv[2], "wt"))
            == NULL)
    {
        fprintf(stderr, "Cannot open output file.\n");
        return 1;
    }

    while(!feof(in))
    {
        fscanf(in,"%u",&n);
        A=(double**)malloc(n*sizeof(double*));
        for(i=0; i<n; i++)
            A[i]=(double*)malloc(n*sizeof(double));
        for(i=1; i<=n; i++)
            for(j=1; j<=n; j++)
                fscanf(in,"%lf",&A[i - 1][j - 1]);
        inverse(A,n,&error);
        if(error==0)
        {
            for(i=1; i<=n; i++)
            {
                for(j=1; j<=n; j++)
                    fprintf(out,"%.12lf ",A[i - 1][j - 1]);
                fprintf(out,"\n");
            }
        }
        else
            fprintf(out,"Macierz osobliwa \n");
        fprintf(out,"\n");
        for(i=0; i<n; i++)
            free(A[i]);
        free(A);
    }
    fclose(in);
    fclose(out);
    return 0;
}