// nacteni standardnich knihoven

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

// nastaveni numericke >nuly< 

#define TINY 1.0e-20

// deklarace pouzitych funkci

int printm(float **,int);
int printv(float *,int);
float *vector(long);
float **matrix(long);
void ludcmp(float **,int,int *,float *);
void lubksb(float **,int,int *,float []);
void mprove(float **,float **,int,int [],float [],float []);

// hlavni program

int main() {
  int n,ii,jj,*indx;
  float **a,**alu,*d,*b,*x;
  
  printf("zadej stupen ctvercove matice >A<\n");
  printf("n=");
  scanf("%d",&n);
  d=(float *)malloc((size_t)sizeof(float));
  indx=(int *)malloc((size_t)n*sizeof(float));
  a=matrix(n);
  alu=matrix(n);
  b=vector(n);
  x=vector(n);
  printf("zadej prvky matice >A<\n");
  for (ii=0;ii<n;ii++) {
    for (jj=0;jj<n;jj++) {
      printf("a(%d,%d)=",ii+1,jj+1);
      scanf("%f",a[ii]+jj);
      alu[ii][jj]=a[ii][jj];
    }
  }
  printf("zadej prvky vektoru >b<\n");  
  for (ii=0;ii<n;ii++) {
    printf("b(%d)=",ii+1);
    scanf("%f",&b[ii]);
    x[ii]=b[ii];
  }  
  printf("matice >A<");
  printm(a,n);
  printf("vektor >b<");
  printv(b,n);
  ludcmp(alu,n,indx,d);
  lubksb(alu,n,indx,x);
  printf("LU rozklad matice >A<");
  printm(alu,n);
  printf("reseni >x< soustavy rovnic >Ax=b<");
  printv(x,n);
  printf("zpresnene reseni >x< soustavy rovnic >Ax=b<");
  mprove(a,alu,n,indx,b,x);
  printv(x,n);
  free(a);
  free(alu);
  free(b);
  free(x);
  free(indx);
  free(d);
  
  return 0;
}

// funkce vypisujici matici

int printm(float **a,int n) {
  int ii,jj;

  for (ii=0;ii<n;ii++) {
    printf("\n[");
    for (jj=0;jj<n;jj++)
      printf("%+1.2e ",a[ii][jj]);
    printf("\b]");
  }
  printf("\n");
  
  return 0;
}

// funkce vypisujci vektor

int printv(float *b,int n) {
  int ii;
 
  printf("\n"); 
  for (ii=0;ii<n;ii++)
    printf("[%+1.2e]\n",b[ii]);
    
  return 0;
}  

// funkce alokujici pamet pro vektor

float *vector(long n) {
  float *v;
  v=(float *)malloc((size_t)(n*sizeof(float)));
  if (!v) printf("allocation failure in vector()\n");
  return v;
}

// funkce alokujici pamet pro matici

float **matrix(long n) {
  int ii;
  float **m;
  m=(float **)malloc((size_t)(n*sizeof(float *)));
  for (ii=0;ii<n;ii++)
    m[ii]=(float *)malloc((size_t)n*sizeof(float));
  if (!m) printf("allocation failure in matrix()\n");
  return m;
}

// LU dekompozice matice 

void ludcmp(float **a,int n,int *indx,float *d) {
  int i,imax,j,k;
  float big,dum,sum,temp;
  float *vv;

  vv=vector(n);
  *d=1.0;
  for (i=0;i<n;i++) {
    big=0.0;
      for (j=0;j<n;j++) {
        if ((temp=fabs(a[i][j]))>big) {
          big=temp;
        }
      }
    if (big==0.0) printf("Singular matrix in routine ludcmp\n");
    vv[i]=1.0/big;
  }
  for (j=0;j<n;j++) {
    for (i=0;i<j;i++) {
      sum=a[i][j];
      for (k=0;k<i;k++)
        sum-=a[i][k]*a[k][j];
      a[i][j]=sum;
    }
    big=0.0;
    for (i=j;i<n;i++) {
      sum=a[i][j];
      for (k=0;k<j;k++) 
        sum-=a[i][k]*a[k][j];
      a[i][j]=sum;
      if ((dum=vv[i]*fabs(sum))>=big) {
        big=dum;
        imax=i;
      }
    } 
    if (j!=imax) {
      for (k=0;k<n;k++) {
        dum=a[imax][k];
        a[imax][k]=a[j][k];
        a[j][k]=dum;
      }
      *d=-(*d);
      vv[imax]=vv[j];
    }
    indx[j]=imax;
    if (a[j][j]==0.0)
      a[j][j]=TINY;
    if (j!=n) {
      dum=1.0/(a[j][j]);
      for (i=j+1;i<n;i++)
        a[i][j]*=dum;
    }
  }
  free(vv);
}

// zpetna procedura Gaussovy eliminace z LU dekompozice matice
 
void lubksb(float **a,int n,int *indx,float b[]) {
  int i,ii=-1,ip,j;
  float sum;
  
  for (i=0;i<n;i++) {
    ip=indx[i];
    sum=b[ip];
    b[ip]=b[i];
    if (ii>-1)
      for (j=ii;j<=i-1;j++)
        sum-=a[i][j]*b[j];
    else if (sum)
      ii=i;
    b[i]=sum;
  }
  for (i=n-1;i>=0;i--) {
    sum=b[i];
    for (j=i+1;j<n;j++)
      sum-=a[i][j]*b[j];
    b[i]=sum/a[i][i];
  }
}

// oprava vyslednoho reseni linearnich rovnic

void mprove(float **a,float **alud,int n,int indx[],float b[],float x[]) {
  int j,i;
  double sdp;
  float *rr;

  rr=vector(n);
  for (i=0;i<n;i++) {
    sdp=-b[i];
    for (j=0;j<n;j++)
      sdp+=a[i][j]*x[j];
    rr[i]=sdp;
  }
  lubksb(alud,n,indx,rr);
  for (i=0;i<n;i++)
    x[i]-=rr[i];
  free(rr);
}