Lu decomposition

#include <iostream.h>  /*this one is correct but may not be contigious */#include <stdlib.h>typedef float **matrix;    /* By Anil Pedgaonkar */typedef float *vector;void outvector(vector x,int n){   cout << "\n";for(int  i=1; i <=n; i++)   cout << " x" << i << " = " << x[i];}int lu( matrix a, int n){ int k, r, i,j, pos; vector temp;  float pivot,value ;  const float epsilon =0.000001;  if ( a[0][0] != 1.0)  {  for( k=1; k < n; k++)  // go over all further eqns to right{ pivot =abs(a[k][k]);  pos =k;  for ( r= k+1; r <=n ; r++)  // go over all row entirs below{if ( (value =abs(a[r][k])) > pivot){ pos =r; pivot = value;} } // finished scanning the eqn for pivot so r loop is over  if ( pivot < epsilon)   {cout << "\n error in lu pivot zero ";goto stop;}  else   {temp = a[pos];a[pos] = a[k];a[k] = temp;} for ( r = k+1; r <=n; r++)  a[r][k] /= a[k][k];for ( j =k+1; j <= n; j++)for ( i =k+1; i <= n; i++)a[i][j] -= a[k][j] *a[i][k];}   a[0][0] = 1.0; // signals decomposedmatrix   }   return(0);stop:  return(1);  } void fwd(matrix a, int n) { int r, k;   float sum =0.0;   for ( k =1; k <=n ; k++) // go over eqns that is cols   {    sum =0.0;for ( r =1; r <k; r++)   sum += a[r][k]*a[0][r];a[0][k] -= sum;a[0][k] /= a[k][k];   } } void bwd ( matrix a, vector x, int n) { int r,k, i;   float sum = 0.0;   i = (int) (a[n][0]);   x[i] = a[0][n];   for( k =n-1; k >= 1; k--) // go over eqns that is cols{ for ( r =n; r >k ; r--) // go over variables{ i = (int)( a[r][0]);// extract the name of variable   sum =0.0;  sum += a[r][k] * x[i];} a[0][k] -= sum; i = (int)(a[k][0]); x[i] = a[0][k];}  } int main() { int r, k, m, n, state =0; const char space = ' '; vector x; matrix a; cout << "please give number n of eqns and variables"; cin >> n; m = n; x = new float[n+1]; for ( k =0; k <=n; k++) x[k] =0; a = new vector[m+1]; for ( r = 0; r <= m ; r++) {a[r] = new float[n+1];for ( k = 0 ; k <=n; k++)  a[r][k] = r; } for ( r=1; r <=n ; r++) { cout << "give coefficients of lhs of  eqn no " << r<< space;;for (k=1; k <=n; k++)    cin >> a[k][r]; } cout << "give rhs vector "; for ( k =1; k <=n; k++)    cin >> a[0][k]; if ( !(state =lu(a, n)))  {fwd(a, n);bwd(a, x, n);  } outvector(x, n); return (0);}