/************************************
	 "In Place" LU Decomposition 
	 ---------------------------
Input: a - pointer to [n x n] matrix
	   n - matrix dimenssions
Output: is writen to a
		where l[i][i]=1 [i=0...n-1]
Important note: This routine do not 
perform  pivoting. The final result 
depend on the user's partial pivoting
*************************************/
void ludecomp(float **a,int n)
{
	int i,j,k;
	
	for( j=0 ; j<n ; j++ )
	{
		for( i=0 ; i<=j ; i++ )
		{
			for( k=0 ; k<=i-1 ; k++ )
				a[i][j] -= a[i][k]*a[k][j];
		}
		for( i=j+1 ; i<n ; i++ )
		{
			for( k=0 ; k<=j-1 ; k++ )
				a[i][j] -=a [i][k]*a[k][j];
			a[i][j] /= a[j][j];
		}
	}
}