/*Translated by FOR_C, v3.4.2 (-), on 07/09/115 at 08:30:55 */
/*FOR_C Options SET: ftn=u io=c no=p op=aimnv s=dbov str=l x=f - prototypes */
#include <math.h>
#include "fcrt.h"
#include "sgbsl.h"
#include <stdlib.h>
void /*FUNCTION*/ sgbsl(
float *abd,
long lda,
long n,
long ml,
long mu,
long ipvt[],
float b[],
long job)
{
#define ABD(I_,J_)	(*(abd+(I_)*(lda)+(J_)))
	long int k, kb, l, la, lb, lm, m, nm1;
	float t;
		/* OFFSET Vectors w/subscript range: 1 to dimension */
	float *const B = &b[0] - 1;
	long *const Ipvt = &ipvt[0] - 1;
		/* end of OFFSET VECTORS */
 
	/* Copyright (c) 1996 California Institute of Technology, Pasadena, CA.
	 * ALL RIGHTS RESERVED.
	 * Based on Government Sponsored Research NAS7-03001.
	 *>> 2001-11-04 SGBSL Krogh  Fixes for F77 and conversion to single
	 *--S replaces "?": ?GBSL, ?GBCO, ?GBFA, ?DOT, ?AXPY
	 *      IMPLICIT NONE
	 ****BEGIN PROLOGUE  SGBSL
	 ****PURPOSE  Solve the real band system A*X=B or TRANS(A)*X=B using
	 *            the factors computed by SGBCO or SGBFA.
	 ****LIBRARY   SLATEC (LINPACK)
	 ****CATEGORY  D2A2
	 ****TYPE      DOUBLE PRECISION (SGBSL-S, SGBSL-D, CGBSL-C)
	 ****KEYWORDS  BANDED, LINEAR ALGEBRA, LINPACK, MATRIX, SOLVE
	 ****AUTHOR  Moler, C. B., (U. of New Mexico)
	 ****DESCRIPTION
	 *
	 *     SGBSL solves the double precision band system
	 *     A * X = B  or  TRANS(A) * X = B
	 *     using the factors computed by SGBCO or SGBFA.
	 *
	 *     On Entry
	 *
	 *        ABD     DOUBLE PRECISION(LDA, N)
	 *                the output from SGBCO or SGBFA.
	 *
	 *        LDA     INTEGER
	 *                the leading dimension of the array  ABD .
	 *
	 *        N       INTEGER
	 *                the order of the original matrix.
	 *
	 *        ML      INTEGER
	 *                number of diagonals below the main diagonal.
	 *
	 *        MU      INTEGER
	 *                number of diagonals above the main diagonal.
	 *
	 *        IPVT    INTEGER(N)
	 *                the pivot vector from SGBCO or SGBFA.
	 *
	 *        B       DOUBLE PRECISION(N)
	 *                the right hand side vector.
	 *
	 *        JOB     INTEGER
	 *                = 0         to solve  A*X = B ,
	 *                = nonzero   to solve  TRANS(A)*X = B , where
	 *                            TRANS(A)  is the transpose.
	 *
	 *     On Return
	 *
	 *        B       the solution vector  X .
	 *
	 *     Error Condition
	 *
	 *        A division by zero will occur if the input factor contains a
	 *        zero on the diagonal.  Technically this indicates singularity
	 *        but it is often caused by improper arguments or improper
	 *        setting of LDA .  It will not occur if the subroutines are
	 *        called correctly and if SGBCO has set RCOND .GT. 0.0
	 *        or SGBFA has set INFO .EQ. 0 .
	 *
	 *     To compute  INVERSE(A) * C  where  C  is a matrix
	 *     with  P  columns
	 *           CALL SGBCO(ABD,LDA,N,ML,MU,IPVT,RCOND,Z)
	 *           IF (RCOND is too small) GO TO ...
	 *           DO 10 J = 1, P
	 *              CALL SGBSL(ABD,LDA,N,ML,MU,IPVT,C(1,J),0)
	 *        10 CONTINUE
	 *
	 ****REFERENCES  J. J. Dongarra, J. R. Bunch, C. B. Moler, and G. W.
	 *                 Stewart, LINPACK Users' Guide, SIAM, 1979.
	 ****ROUTINES CALLED  SAXPY, SDOT
	 ****REVISION HISTORY  (YYMMDD)
	 *   780814  DATE WRITTEN
	 *   890531  Changed all specific intrinsics to generic.  (WRB)
	 *   890831  Modified array declarations.  (WRB)
	 *   890831  REVISION DATE from Version 3.2
	 *   891214  Prologue converted to Version 4.0 format.  (BAB)
	 *   900326  Removed duplicate information from DESCRIPTION section.
	 *           (WRB)
	 *   920501  Reformatted the REFERENCES section.  (WRB)
	 ****END PROLOGUE  SGBSL */
 
	/****FIRST EXECUTABLE STATEMENT  SGBSL */
	m = mu + ml + 1;
	nm1 = n - 1;
	if (job != 0)
		goto L_50;
 
	/*        JOB = 0 , SOLVE  A * X = B
	 *        FIRST SOLVE L*Y = B
	 * */
	if (ml == 0)
		goto L_30;
	if (nm1 < 1)
		goto L_30;
	for (k = 1; k <= nm1; k++)
	{
		lm = min( ml, n - k );
		l = Ipvt[k];
		t = B[l];
		if (l == k)
			goto L_10;
		B[l] = B[k];
		B[k] = t;
L_10:
		;
		saxpy( lm, t, &ABD(k - 1,m), 1, &B[k + 1], 1 );
	}
L_30:
	;
 
	/*        NOW SOLVE  U*X = Y
	 * */
	for (kb = 1; kb <= n; kb++)
	{
		k = n + 1 - kb;
		B[k] /= ABD(k - 1,m - 1);
		lm = min( k, m ) - 1;
		la = m - lm;
		lb = k - lm;
		t = -B[k];
		saxpy( lm, t, &ABD(k - 1,la - 1), 1, &B[lb], 1 );
	}
	goto L_100;
L_50:
	;
 
	/*        JOB = NONZERO, SOLVE  TRANS(A) * X = B
	 *        FIRST SOLVE  TRANS(U)*Y = B
	 * */
	for (k = 1; k <= n; k++)
	{
		lm = min( k, m ) - 1;
		la = m - lm;
		lb = k - lm;
		t = sdot( lm, &ABD(k - 1,la - 1), 1, &B[lb], 1 );
		B[k] = (B[k] - t)/ABD(k - 1,m - 1);
	}
 
	/*        NOW SOLVE TRANS(L)*X = Y
	 * */
	if (ml == 0)
		goto L_90;
	if (nm1 < 1)
		goto L_90;
	for (kb = 1; kb <= nm1; kb++)
	{
		k = n - kb;
		lm = min( ml, n - k );
		B[k] += sdot( lm, &ABD(k - 1,m), 1, &B[k + 1], 1 );
		l = Ipvt[k];
		if (l == k)
			goto L_70;
		t = B[l];
		B[l] = B[k];
		B[k] = t;
L_70:
		;
	}
L_90:
	;
L_100:
	;
	return;
#undef	ABD
} /* end of function */