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);
}