/* dgesl.c -- double precision solution of general matrix by
 *	      Gaussian elimination.
 *
 *	Solves the system A * X = B or TRANS(A) * X = B.
 *
 *	Operates on the LU decomposition (factorization)
 *	generated by dgefa.c or dgeco.c.
 *
 *	From "Linpack Users' Guide", Dongarra, et al.
 *	Translated from FORTRAN by C. Bond, 1994.
 */
#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 dgesl(double **a,int n,int ipvt[],double b[],int job)
{
	double t;
	int i,k,kb,l,nm1;
	
	nm1 = n - 1;
	if (job != 0) goto _50;

/* job = 0, solve a * x = b.
 * First solve l * y = b.
 */
	if (nm1 < 0) goto _30;
	for (k = 0; k < nm1; k++) {
		l = ipvt[k];
		t = b[l];
		if (l == k) goto _10;
		b[l] = b[k];
		b[k] = t;
_10:
/*		saxpy(n-k,t,a(k+1,k),1,b(k+1),1); */
		for (i=k+1;i < n;i++)
			b[i] += (t * a[i][k]); 
	}
_30: 
/* Now solve u * x = y. */
	for (kb = 0; kb < n; kb++) {
		k = n - kb-1;
		b[k] /= a[k][k];
		t = -b[k];
/*		saxpy(k-1,t,a(1,k),1,b(1),1); */
		for (i = 0; i < k ; i++)
			b[i] += (t * a[i][k]);
	}
	goto _100;
_50:
/* job != 0, solve trans(a) * x = b.
 * First solve trans(u) * x = y.
 */
 	for (k = 0; k < n; k++) {
/* 		t = ddot(k-1,a(1,k),1,b(1),1); */
		t = 0;
		for (i = 0; i < k; i++)
			t += (a[i][k] * b[i]);
		b[k] = (b[k] - t)/a[k][k];
 	}

/* Now solve trans(l) * x = y. */
	if (nm1 < 0) goto _90;
	for (kb = 0; kb < nm1; kb++) {
		k = n - 2 - kb;
/*		b[k] = b[k] + ddot(n-k,a(k+1,k),1,b(k+1),1); */
		t = 0;
		for (i = k+1;i < n; i++)
			t += (a[i][k] * b[i]);
		b[k] += t;
		l = ipvt[k];
		if (l == k) goto _70;
		t = b[l];
		b[l] = b[k];
		b[k] = t;
_70:
	}	
_90:
_100:
}