/* dgefa.c -- double precision solution of general matrix by
 *	      Gaussian elimination.
 *
 *	This routine performs the LU factorization. (You should
 *	use dgeco.c if you need a condition estimate.)
 *	
 *	Follow this with a call to dgesl.c, which is the linear
 *	equation solver, or dgedi.c which computes the determinant
 *	and inverse.
 *
 *	From "Linpack Users' Guide", Dongarra, et al.
 *	Translated from FORTRAN by C. Bond, 1994.	
 */
#include <stdio.h>
#include <math.h>

#define max(a,b)	(((a) > (b)) ? (a) : (b))
#define min(a,b)	(((a) < (b)) ? (a) : (b))
 
void dgefa(double **a,int n, int ipvt[],int *info)
{
	double dmax,t;
	int i,j,k,kp1,l,nm1;

	*info = 0;
	nm1 = n - 1;
	if (nm1 < 0) goto _70;
	for (k = 0; k < nm1; k++) {
		kp1 = k + 1;
/* Find l = pivot index. */
		dmax = fabs(a[k][k]);
		l = k;
		for (i = k+1; i < n; i++) {
			if (fabs(a[i][k]) <= dmax) continue;
			l = i;
		}
		ipvt[k] = l;
	
/* Zero pivot implies this column already triangularized. */
		if (a[l][k] == 0.0) goto _40;

/* Interchange if necessary. */
		if (l == k) goto _10;
		t = a[l][k];
		a[l][k] = a[k][k];
		a[k][k] = t;
_10:
		
/* Compute multipliers. */
		if (a[k][k] == 0.0) printf("\n!ERROR. Singular matrix.\m");
		t = -1.0/a[k][k];
		for (i = k+1; i < n; i++)
			a[i][k] *= t;
		
/* Row elimination with column indexing. */
		for (j = kp1; j < n; j++) {
			t = a[l][j];
			if (l == k) goto _20;
			a[l][j] = a[k][j];
			a[k][j] = t;
_20:
			for (i = k+1; i < n; i++ )
				a[i][j] += (t * a[i][k]);
		}
		goto _50;
_40:
		*info = k;
_50:
	}
_70:	
				
	ipvt[n-1] = n-1;
	if (a[n-1][n-1] == 0.0) *info = n-1;
}