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