/* 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 (n > 0) {
        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) {
                *info = k;
                continue;
            }

/* Interchange if necessary. */
            if (l != k) {
                t = a[l][k];
                a[l][k] = a[k][k];
                a[k][k] = t;
            }
		
/* Compute multipliers. */
            if (a[k][k] == 0.0) printf("\n!ERROR. Singular matrix.\n");
            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) {
                    a[l][j] = a[k][j];
                    a[k][j] = t;
                }
                for (i = k+1; i < n; i++ )
                    a[i][j] += (t * a[i][k]);
            }
        }
    } 
	ipvt[n-1] = n-1;
	if (a[n-1][n-1] == 0.0) *info = n-1;
}