/*	dgtsl.c --- LINPACK solver for general tridiagonal matrix. */
#include <stdio.h>
#include <float.h>
#include <math.h>

/* This is a translation of the FORTRAN routine from LINPACK.
 *
 * INPUT parameters:
 *	n	- order of tridiagonal matrix,
 *	c	- singly subscripted array with subdiagonal in
 *		  positions c(1) to c(n-1),
 *	d	- singly subscripted array with diagonal in
 *		  positions d(0) to d(n-1),
 *	e	- singly subscripted array with superdiagonal in
 *		  positions e(0) to e(n-2),
 *	b	- singly subscripted array with right hand side
 *		  of system,
 *
 * RETURN values:
 *	The routine returns an integer which is 0 on successful
 *	completion and contains the index+1 of the first zero
 *	pivot found during the solution process otherwise. 
 *	b	- solution vector,
 */
  
int dgtsl(double c[],double d[],double e[],double b[],int n)
{
	double t;
	int nm1,k,kb,kp1,nm2,retval;
	
/* Initialize return value. */
	retval = 0;

/* Use C array indexing rules. */
	n--;
		
	c[0] = d[0];
	nm1 = n-1;
	if (nm1 < 0) goto _40;
	d[0] = e[0];
	e[0] = 0.0;
	e[n] = 0.0;
	for (k = 0; k <= nm1; k++) {
		kp1 = k + 1;
/* find the largest of the two rows */
		if (fabs(c[kp1]) < fabs(c[k]))
			goto _10;
/* interchange rows */
		t = c[kp1];
		c[kp1] = c[k];
		c[k] = t;
		t = d[kp1];
		d[kp1] = d[k];
		d[k] = t;
		t = e[kp1];
		e[kp1] = e[k];
		e[k] = t;
		t = b[kp1];
		b[kp1] = b[k];
		b[k] = t;
_10:
/* zero elements */		
		if (c[k] == 0.0) return k+1;
/* .............exit */

		t = -c[kp1] / c[k];
		c[kp1] = d[kp1] + t * d[k];
		d[kp1] = e[kp1] + t * e[k];
		e[kp1] = 0.0;
		b[kp1] = b[kp1] + t * b[k];
	}
_40:
	if (c[n] == 0.0) return n+1;

/* back solve */
	nm2 = n - 2;
	b[n] = b[n] / c[n];
	if (n == 0) goto _80;
		b[nm1] = (b[nm1] - d[nm1] * b[n]) / c[nm1];
		if (nm2 < 0) goto _70;
		for (kb = 0; kb <= nm2; kb++) {
			k = nm2 - kb;
			b[k] = (b[k] - d[k] * b[k+1] - e[k] * b[k+2]) / c[k];
		}
		/* endif */
_70:
	/* endif */
_80:
	return retval;
}