/*Translated by FOR_C, v3.4.2 (-), on 07/09/115 at 08:31:45 */
/*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 "sbsol.h"
#include <stdlib.h>
		/* PARAMETER translations */
#define	ZERO	0.0e0
		/* end of PARAMETER translations */
 
void /*FUNCTION*/ sbsol(
long mode,
float *g,
long ldg,
long nb,
long ir,
long jtprev,
float x[],
long n,
float *rnorm,
long *ierr3)
{
#define G(I_,J_)	(*(g+(I_)*(ldg)+(J_)))
	long int i, i1, i2, ie, ii, irm1, ix, j, jg, l, np1;
	float rsq, s;
		/* OFFSET Vectors w/subscript range: 1 to dimension */
	float *const X = &x[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.
	 *>> 1994-11-11 SBSOL  Krogh  Declared all vars.
	 *>> 1994-10-20 SBSOL  Krogh  Changes to use M77CON
	 *>> 1987-11-24 SBSOL  Lawson Initial code.
	 *--S replaces "?": ?BSOL
	 *     Solution of banded least squares problem following use of _BACC to
	 *     accumulate the data.
	 *     ------------------------------------------------------------------
	 *                            Subroutine arguments
	 *
	 *     MODE  [in]  Flag selecting operation to be done.
	 *          = 1     SOLVE R*X=Y   WHERE R AND Y ARE IN THE G( ) ARRAY
	 *                  AND X WILL BE STORED IN THE X( ) ARRAY.
	 *            2     SOLVE (R**T)*X=Y   WHERE R IS IN G( ),
	 *                  Y IS INITIALLY IN X( ), AND X REPLACES Y IN X( ),
	 *            3     SOLVE R*X=Y   WHERE R IS IN G( ).
	 *                  Y IS INITIALLY IN X( ), AND X REPLACES Y IN X( ).
	 *
	 *     G(,)  [in]  Array in which the transformed problem data has
	 *           been left by _BACC.
	 *     LDG   [in]  Leading dimensioning parameter for G(,).
	 *     NB    [in]  Bandwidth of data matrix.
	 *     IR    [in]  Must have value set by previous call to _BACC.
	 *     JTPREV    [in]  Must have value set by previous call to _BACC.
	 *     X()   [inout]  Array of length at least N.
	 *           On entry with MODE = 2 or 3 contains the y vector.
	 *           In all (non-error) cases contains the solution vector on
	 *           return.
	 *     N     [in]  Specifies the dimension of the desired solution
	 *           vector.  Can be set smaller than the max value the data
	 *           would support.
	 *     RNORM  [out]  Set to the euclidean norm of the residual vector
	 *           when MODE = 1.  Set to zero otherwise.
	 *     IERR3  [out]  =0 means no errors detected.
	 *                  =1 means a diagonal element is zero.
	 *     ------------------------------------------------------------------
	 *     C.L.LAWSON AND R.J.HANSON, JET PROPULSION LABORATORY, 1973 JUN 12
	 *     Reference:   'SOLVING LEAST SQUARES PROBLEMS', PRENTICE-HALL, 1974
	 *     Comments in this code refer to Algorithm Step numbers on
	 *     pp 213-217 of the book.
	 *     ------------------------------------------------------------------
	 *     1984 July 11, C. L. Lawson, JPL.  Adapted code from book
	 *     to Fortran 77 for JPL MATH 77 library.
	 *     Prefixing subprogram names with S or D for s.p. or d.p. versions.
	 *     Using generic names for any intrinsic functions.
	 *     July 1987.  Using _HTCC in place of _H12.
	 *     ------------------------------------------------------------------
	 *     Subprograms called: IERM1, IERV1
	 *     ------------------------------------------------------------------ */
	/*     ------------------------------------------------------------------
	 * */
	*ierr3 = 0;
	*rnorm = ZERO;
	if (mode == 1)
	{
		/*                                   ********************* MODE = 1
		 *                                   ALG. STEP 26 */
		for (j = 1; j <= n; j++)
		{
			X[j] = G(nb,j - 1);
		}
		rsq = ZERO;
		np1 = n + 1;
		irm1 = ir - 1;
		if (np1 <= irm1)
		{
			for (j = np1; j <= irm1; j++)
			{
				rsq += SQ(G(nb,j - 1));
			}
			*rnorm = sqrtf( rsq );
		}
	}
 
	if (mode == 2)
		goto L_90;
	/*                                   ********************* MODE = 1 or 3
	 *                                   ALG. STEP 27 */
	for (ii = 1; ii <= n; ii++)
	{
		i = n + 1 - ii;
		/*                                   ALG. STEP 28 */
		s = ZERO;
		l = max( 0, i - jtprev );
		/*                                   ALG. STEP 29 */
		if (i != n)
		{
			/*                                   ALG. STEP 30 */
			ie = min( n + 1 - i, nb );
			for (j = 2; j <= ie; j++)
			{
				jg = j + l;
				ix = i - 1 + j;
				s += G(jg - 1,i - 1)*X[ix];
			}
		}
		/*                                   ALG. STEP 31 */
		if (G(l,i - 1) == ZERO)
			goto L_130;
		X[i] = (X[i] - s)/G(l,i - 1);
	}
	/*                                   ALG. STEP 32 */
	return;
	/*     -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
	 *                                   ********************* MODE = 2 */
L_90:
	for (j = 1; j <= n; j++)
	{
		s = ZERO;
		if (j != 1)
		{
			i1 = max( 1, j - nb + 1 );
			i2 = j - 1;
			for (i = i1; i <= i2; i++)
			{
				l = j - i + 1 + max( 0, i - jtprev );
				s += X[i]*G(l - 1,i - 1);
			}
		}
		l = max( 0, j - jtprev );
		if (G(l,j - 1) == ZERO)
		{
			i = j;
			goto L_130;
		}
		X[j] = (X[j] - s)/G(l,j - 1);
	}
	return;
	/*     -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
L_130:
	;
	*ierr3 = i;
	ierm1( "SBSOL", *ierr3, 0, "Singular matrix", "MODE", mode, ',' );
	ierv1( "Index of zero diag elt", i, '.' );
	return;
#undef	G
} /* end of function */