/* dgedi.c -- double precision solution of general matrix by
 *	      Gaussian elimination. computes determinant and
 *	      inverse of matrix (following dgeco or dgefa).
 *
 *	From "Linpack Users' Guide", Dongarra, et al.,
 */
#include <stdio.h>
#include <math.h>
#include "func.h"

#define max(a,b)	(((a) > (b)) ? (a) : (b))
#define min(a,b)	(((a) < (b)) ? (a) : (b))
 
void dgedi(double **a,int n,int ipvt[],double det[],double work[],int job)
{
	double t,ten;
	int i,j,k,kb,kp1,l,nm1;
	
/* Compute determinant. */
	if ((job / 10) == 0) goto _70;
	det[0] = 1.0;
	det[1] = 0.0;
	ten = 10.0;
	for (i = 0; i < n; i++) {
		if (ipvt[i] != i) det[0] = -det[0];
		det[0] *= a[i][i];
		if (det[0] == 0.0) goto _60;
		while (fabs(det[0]) < 1.0) {
			det[0] *= ten;
			det[1] -= 1.0;
		}
		while (fabs(det[0]) >= ten) {
			det[0] /= ten;
			det[1] += 1.0;
		}
_60:
	}
_70:

/* Compute inverse. */
	if (job == 10) return;
	for (k = 0; k < n; k++) {
		a[k][k] = 1.0 / a[k][k];
		t = -a[k][k];
		/* dscal(k-1,t,a(1,k),1) */
		for (i = 0; i < k; i++)
			a[i][k] *= t;
		kp1 = k + 1;
		if ( n < kp1+1) goto _90;
		for (j = kp1; j < n; j++) {
			t = a[k][j];
			a[k][j] = 0.0;
			/* daxpy(k,t,a(1,k),1,a(1,j),1) */
			for (i = 0; i < k+1; i++)
				a[i][j] += (t * a[i][k]);
		}
_90:
	}

/* Form inverse(u) * inverse(l) */
	nm1 = n - 1;
	if (nm1 < 0) goto _140;
	for (kb = 0; kb < nm1; kb++) {
		k = n - kb - 2;
		kp1 = k + 1;
		for (i = kp1; i < n; i++) {
			work[i] = a[i][k];
			a[i][k] = 0.0;
		}
		for (j = kp1; j < n; j++) {
			t = work[j];
			/* daxpy(n,t,a(1,j),1,a(1,k),1) */
			for (i = 0; i < n; i++)
				a[i][k] += (t * a[i][j]);
		}
		l = ipvt[k];
		if (l != k) { 
			/* dswap(n,a(1,k),1,a(1,l),1) */	
			for (i = 0; i < n; i++) {
				t = a[i][k];
				a[i][k] = a[i][l];
				a[i][l] = t;
			}
		}
	}
_140:
}