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

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

{
float amax,temp,c,sum;
short n,i,j,k,k1,l,l1,itemp;
row y;
short 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 %hd\n",k);
  amax=0.0;
  for (l1=k;l1<=n;l1++)
    if (amax<fabs(ab[l1][k]))
    {
      amax=fabs(ab[l1][k]);
      l=l1;
    }
  if (amax==0.0)
  {
    printf("singular matrix\n");
    abort();
  }
  printf("pivoting row %hd\n",l);
  itemp=piv[l];
  piv[l]=piv[k];
  piv[k]=itemp;
  
/* This corrects an error in Lord, J.ACM 30 (1) */
  for (k1=1;k1<=k;k1++)
  {
    temp=ab[l][k1];
    ab[l][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[l][j];
    ab[l][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("%hd %hd\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];
}
}

main()
{
short i,j,k;
float sum,L,U;

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

printf("enter matrix and column\n");
for (i=1;i<=num_slots;i++)
  for (j=1;j<=num_slots+1;j++)
    scanf("%f",&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 and 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("plugging solution into original 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);
}
}