#include <stdio.h>
#include <math.h>
/*Solves Least Squares based on CLR, p. 771*/

#define MAX_SLOTS 10
typedef char boolean;
short num_slots,i,j,num_points,num_functions;
float basisFunc[MAX_SLOTS+1][MAX_SLOTS+1];
float slot[MAX_SLOTS+1];
typedef float row[MAX_SLOTS+2];
row ab[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 y[MAX_SLOTS+1];

printf("enter n (number of points): ");
scanf("%hd",&num_points);
printf("enter m (number of functions): ");
scanf("%hd",&num_functions);
printf("enter %hd x %hd basis-function matrix\n",num_points,
  num_functions);
for (i=1;i<=num_points;i++)
  for (j=1;j<=num_functions;j++)
    scanf("%f",&basisFunc[i][j]);
printf("enter %hd y values\n",num_points);
for (i=1;i<=num_points;i++)
  scanf("%f",&y[i]);

num_slots=num_functions;
/*Load system to be solved*/
/*a is loaded with the product of basisFunc transposed 
  and basicFunc */
for (i=1;i<=num_functions;i++)
  for (j=1;j<=num_functions;j++)
  {
    ab[i][j]=0.0;
    for (k=1;k<=num_points;k++)
      ab[i][j]+=basisFunc[k][i]*basisFunc[k][j];
  }
/*b is loaded with the product of basisFunc transposed
  and y*/
for (i=1;i<=num_functions;i++)
{
  ab[i][num_functions+1]=0.0;
  for (j=1;j<=num_points;j++)
    ab[i][num_functions+1]+=basisFunc[j][i]*y[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("coefficients\n");
for (i=1;i<=num_slots;i++)
  printf("%f\n",slot[i]);
}