#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#define MAX_SLOTS 10
typedef char boolean;
int num_slots,i,j;
double slot[MAX_SLOTS+1];
typedef double row[MAX_SLOTS+2];
row ab[MAX_SLOTS+2],ab1[MAX_SLOTS+2];

void LU()
//SOLVES SYSTEM OF LINEAR EQUATIONS.
//TAKEN FROM M.J. Quinn,P. 147
//CODE HAS BEEN MODIFIED TO PRELOAD AB
//TO SAVE SPACE BY NOT STORING A AND B
//SOLUTION IS PLACED IN SLOT[*]

{
double amax,temp,c,sum;
int n,i,j,k,k1,ell,ell1,itemp;
row y;
int piv[MAX_SLOTS+1];

n=num_slots;
for (i=1;i<=n;i++)
  piv[i]=i;

for (k=1;k<=n-1;k++)
{
  printf("iteration %d\n",k);
  amax=0.0;
  for (ell1=k;ell1<=n;ell1++)
    if (amax<fabs(ab[ell1][k]))
    {
      amax=fabs(ab[ell1][k]);
      ell=ell1;
    }
  if (amax==0.0)
  {
    printf("singular matrix\n");
    abort();
  }
  printf("pivoting row %d\n",ell);
  itemp=piv[ell];
  piv[ell]=piv[k];
  piv[k]=itemp;
  
  for (k1=1;k1<=k;k1++)
  {
    temp=ab[ell][k1];
    ab[ell][k1]=ab[k][k1];
    ab[k][k1]=temp;
  }

  c=1.0/ab[k][k];
  for (i=k+1;i<=n;i++)
    ab[i][k] *= c;
    
  for (j=k+1;j<=n;j++)
  {
    temp=ab[ell][j];
    ab[ell][j]=ab[k][j];
    ab[k][j]=temp;
    
    for (i=k+1;i<=n;i++)
      ab[i][j] -= ab[i][k]*ab[k][j];
  }

 
  printf("LU matrix\n");
  for (i=1;i<=n;i++)
  {
    for (j=1;j<=n+1;j++)
      printf("%f ",ab[i][j]);
    printf("\n");
  }
}
printf("pivots\n");
for (i=1;i<=n;i++)
  printf("%d %d\n",i,piv[i]);
    
printf("forward substitution\n");
y[1]=ab[piv[1]][n+1];
for (i=2;i<=n;i++)
{
  sum=0.0;
  for (j=1;j<i;j++)
    sum+=ab[i][j]*y[j];
  y[i]=ab[piv[i]][n+1]-sum;
}

printf("intermediate vector y:\n");
for (i=1;i<=n;i++)
  printf("%f\n",y[i]);

printf("back substitution\n");
slot[n]=y[n]/ab[n][n];
for (i=n-1;i>=1;i--)
{
  sum=0.0;
  for (j=i+1;j<=n;j++)
    sum+=ab[i][j]*slot[j];
  slot[i]=(y[i]-sum)/ab[i][i];
}
}

int main()
{
int i,j,k;
double sum,L,U;

printf("enter n: ");
scanf("%d",&num_slots);

printf("enter matrix and column\n");
for (i=1;i<=num_slots;i++)
  for (j=1;j<=num_slots+1;j++)
    scanf("%lf",&ab[i][j]);
for (i=1;i<=num_slots;i++)
  for (j=1;j<=num_slots+1;j++)
    ab1[i][j]=ab[i][j];
    
printf("input:\n");
for (i=1;i<=num_slots;i++)
{
  for (j=1;j<=num_slots+1;j++)
    printf("%f ",ab[i][j]);
  printf("\n");
}

LU();

printf("solution\n");
for (i=1;i<=num_slots;i++)
  printf("%f\n",slot[i]);

printf("result of multiplying L & U\n");
for (i=1;i<=num_slots;i++)
{
  for (j=1;j<=num_slots;j++)
  {
    sum=0.0;
    for (k=1;k<=num_slots;k++)
    {
      if (i==k)
        L=1.0;
      if (i<k)
        L=0.0;
      if (i>k)
        L=ab[i][k];
      if (k==j)
        U=ab[k][j];
      if (k<j)
        U=ab[k][j];
      if (k>j)
        U=0.0;
      sum+=L*U;
    }
    printf("%f ",sum);
  }
  printf("\n");
}
printf("sub. sol. into orig. system:\n");
for (i=1;i<=num_slots;i++)
{
  sum=0.0;
  for (j=1;j<=num_slots;j++)
    sum+=ab1[i][j]*slot[j];
  printf("%f\n",sum);
}
}