[mode,g,mdg,nb,ip,ir,x,n,rnorm]=bndsol(mode,g,mdg,nb,ip,ir,x,n,rnorm);

function [mode,g,mdg,nb,ip,ir,x,n,rnorm]=bndsol(mode,g,mdg,nb,ip,ir,x,n,rnorm);
persistent i i1 i2 ie igo ii iopt irm1 ix j jg l nerr np1 rsq s zero ; 

if isempty(rsq), rsq=0; end;
if isempty(s), s=0; end;
if isempty(zero), zero=0; end;
if isempty(i), i=0; end;
if isempty(i1), i1=0; end;
if isempty(i2), i2=0; end;
if isempty(ie), ie=0; end;
if isempty(ii), ii=0; end;
if isempty(iopt), iopt=0; end;
if isempty(irm1), irm1=0; end;
if isempty(ix), ix=0; end;
if isempty(j), j=0; end;
if isempty(jg), jg=0; end;
if isempty(l), l=0; end;
if isempty(nerr), nerr=0; end;
if isempty(np1), np1=0; end;
if isempty(igo), igo=0; end;
%***BEGIN PROLOGUE  BNDSOL
%***PURPOSE  Solve the least squares problem for a banded matrix using
%            sequential accumulation of rows of the data matrix.
%            Exactly one right-hand side vector is permitted.
%***LIBRARY   SLATEC
%***CATEGORY  D9
%***TYPE      SINGLE PRECISION (BNDSOL-S, DBNDSL-D)
%***KEYWORDS  BANDED MATRIX, CURVE FITTING, LEAST SQUARES
%***AUTHOR  Lawson, C. L., (JPL)
%           Hanson, R. J., (SNLA)
%***DESCRIPTION
%
%     These subroutines solve the least squares problem Ax = b for
%     banded matrices A using sequential accumulation of rows of the
%     data matrix.  Exactly one right-hand side vector is permitted.
%
%     These subroutines are intended for the type of least squares
%     systems that arise in applications such as curve or surface
%     fitting of data.  The least squares equations are accumulated and
%     processed using only part of the data.  This requires a certain
%     user interaction during the solution of Ax = b.
%
%     Specifically, suppose the data matrix (A B) is row partitioned
%     into Q submatrices.  Let (E F) be the T-th one of these
%     submatrices where E = (0 C 0).  Here the dimension of E is MT by N
%     and the dimension of C is MT by NB.  The value of NB is the
%     bandwidth of A.  The dimensions of the leading block of zeros in E
%     are MT by JT-1.
%
%     The user of the subroutine BNDACC provides MT,JT,C and F for
%     T=1,...,Q.  Not all of this data must be supplied at once.
%
%     Following the processing of the various blocks (E F), the matrix
%     (A B) has been transformed to the form (R D) where R is upper
%     triangular and banded with bandwidth NB.  The least squares
%     system Rx = d is then easily solved using back substitution by
%     executing the statement CALL BNDSOL(1,...). The sequence of
%     values for JT must be nondecreasing.  This may require some
%     preliminary interchanges of rows and columns of the matrix A.
%
%     The primary reason for these subroutines is that the total
%     processing can take place in a working array of dimension MU by
%     NB+1.  An acceptable value for MU is
%
%                       MU = MAX(MT + N + 1),
%
%     where N is the number of unknowns.
%
%     Here the maximum is taken over all values of MT for T=1,...,Q.
%     Notice that MT can be taken to be a small as one, showing that
%     MU can be as small as N+2.  The subprogram BNDACC processes the
%     rows more efficiently if MU is large enough so that each new
%     block (C F) has a distinct value of JT.
%
%     The four principle parts of these algorithms are obtained by the
%     following call statements
%
%     CALL BNDACC(...)  Introduce new blocks of data.
%
%     CALL BNDSOL(1,...)Compute solution vector and length of
%                       residual vector.
%
%     CALL BNDSOL(2,...)Given any row vector H solve YR = H for the
%                       row vector Y.
%
%     CALL BNDSOL(3,...)Given any column vector W solve RZ = W for
%                       the column vector Z.
%
%     The dots in the above call statements indicate additional
%     arguments that will be specified in the following paragraphs.
%
%     The user must dimension the array appearing in the call list..
%     G(MDG,NB+1)
%
%     Description of calling sequence for BNDACC..
%
%     The entire set of parameters for BNDACC are
%
%     Input..
%
%     G(*,*)            The working array into which the user will
%                       place the MT by NB+1 block (C F) in rows IR
%                       through IR+MT-1, columns 1 through NB+1.
%                       See descriptions of IR and MT below.
%
%     MDG               The number of rows in the working array
%                       G(*,*).  The value of MDG should be .GE. MU.
%                       The value of MU is defined in the abstract
%                       of these subprograms.
%
%     NB                The bandwidth of the data matrix A.
%
%     IP                Set by the user to the value 1 before the
%                       first call to BNDACC.  Its subsequent value
%                       is controlled by BNDACC to set up for the
%                       next call to BNDACC.
%
%     IR                Index of the row of G(*,*) where the user is
%                       the user to the value 1 before the first call
%                       to BNDACC.  Its subsequent value is controlled
%                       by BNDACC. A value of IR .GT. MDG is considered
%                       an error.
%
%     MT,JT             Set by the user to indicate respectively the
%                       number of new rows of data in the block and
%                       the index of the first nonzero column in that
%                       set of rows (E F) = (0 C 0 F) being processed.
%     Output..
%
%     G(*,*)            The working array which will contain the
%                       processed rows of that part of the data
%                       matrix which has been passed to BNDACC.
%
%     IP,IR             The values of these arguments are advanced by
%                       BNDACC to be ready for storing and processing
%                       a new block of data in G(*,*).
%
%     Description of calling sequence for BNDSOL..
%
%     The user must dimension the arrays appearing in the call list..
%
%     G(MDG,NB+1), X(N)
%
%     The entire set of parameters for BNDSOL are
%
%     Input..
%
%     MODE              Set by the user to one of the values 1, 2, or
%                       3.  These values respectively indicate that
%                       the solution of AX = B, YR = H or RZ = W is
%                       required.
%
%     G(*,*),MDG,       These arguments all have the same meaning and
%      NB,IP,IR         contents as following the last call to BNDACC.
%
%     X(*)              With mode=2 or 3 this array contains,
%                       respectively, the right-side vectors H or W of
%                       the systems YR = H or RZ = W.
%
%     N                 The number of variables in the solution
%                       vector.  If any of the N diagonal terms are
%                       zero the subroutine BNDSOL prints an
%                       appropriate message.  This condition is
%                       considered an error.
%
%     Output..
%
%     X(*)              This array contains the solution vectors X,
%                       Y or Z of the systems AX = B, YR = H or
%                       RZ = W depending on the value of MODE=1,
%                       2 or 3.
%
%     RNORM             If MODE=1 RNORM is the Euclidean length of the
%                       residual vector AX-B.  When MODE=2 or 3 RNORM
%                       is set to zero.
%
%     Remarks..
%
%     To obtain the upper triangular matrix and transformed right-hand
%     side vector D so that the super diagonals of R form the columns
%     of G(*,*), execute the following Fortran statements.
%
%     NBP1=NB+1
%
%     DO 10 J=1, NBP1
%
%  10 G(IR,J) = 0.0E0
%
%     MT=1
%
%     JT=N+1
%
%     CALL BNDACC(G,MDG,NB,IP,IR,MT,JT)
%
%***REFERENCES  C. L. Lawson and R. J. Hanson, Solving Least Squares
%                 Problems, Prentice-Hall, Inc., 1974, Chapter 27.
%***ROUTINES CALLED  XERMSG
%***REVISION HISTORY  (YYMMDD)
%   790101  DATE WRITTEN
%   890531  Changed all specific intrinsics to generic.  (WRB)
%   890831  Modified array declarations.  (WRB)
%   891006  Cosmetic changes to prologue.  (WRB)
%   891006  REVISION DATE from Version 3.2
%   891214  Prologue converted to Version 4.0 format.  (BAB)
%   900315  CALLs to XERROR changed to CALLs to XERMSG.  (THJ)
%   900326  Removed duplicate information from DESCRIPTION section.
%           (WRB)
%   920501  Reformatted the REFERENCES section.  (WRB)
%***end PROLOGUE  BNDSOL
g_shape=size(g);g=reshape([g(:).',zeros(1,ceil(numel(g)./prod([mdg])).*prod([mdg])-numel(g))],mdg,[]);
x_shape=size(x);x=reshape(x,1,[]);
%***FIRST EXECUTABLE STATEMENT  BNDSOL
zero = 0.;
%
igo=1;
while(igo==1);
igo=0;
rnorm = zero;
if( mode==2 )
%                                   ********************* MODE = 2
for j = 1 : n;
s = zero;
if( j~=1 )
i1 = fix(max(1,j-nb+1));
i2 = fix(j - 1);
for i = i1 : i2;
l = fix(j - i + 1 + max(0,i-ip));
s = s + x(i).*g(i,l);
end; i = fix(i2+1);
end;
l = fix(max(0,j-ip));
if( g(j,l+1)==0 )
igo=0;
break;
end;
x(j) =(x(j)-s)./g(j,l+1);
end;
if(igo==0)
break;
end;
g_shape=zeros(g_shape);g_shape(:)=g(1:numel(g_shape));g=g_shape;
x_shape=zeros(x_shape);x_shape(:)=x(1:numel(x_shape));x=x_shape;
return;
elseif( mode~=3 ) ;
%                                   ********************* MODE = 1
%                                   ALG. STEP 26
for j = 1 : n;
x(j) = g(j,nb+1);
end; j = fix(n+1);
rsq = zero;
np1 = fix(n + 1);
irm1 = fix(ir - 1);
if( np1<=irm1 )
for j = np1 : irm1;
rsq = rsq + g(j,nb+1).^2;
end; j = fix(irm1+1);
rnorm = sqrt(rsq);
end;
end;
%                                   ********************* MODE = 3
%                                   ALG. STEP 27
for ii = 1 : n;
i = fix(n + 1 - ii);
%                                   ALG. STEP 28
s = zero;
l = fix(max(0,i-ip));
%                                   ALG. STEP 29
if( i~=n )
%                                   ALG. STEP 30
ie = fix(min(n+1-i,nb));
for j = 2 : ie;
jg = fix(j + l);
ix = fix(i - 1 + j);
s = s + g(i,jg).*x(ix);
end; j = fix(ie+1);
end;
%                                   ALG. STEP 31
if( g(i,l+1)==0 )
igo=0;
break;
end;
x(i) =(x(i)-s)./g(i,l+1);
end;
if(igo==0)
break;
end;
%                                   ALG. STEP 32
g_shape=zeros(g_shape);g_shape(:)=g(1:numel(g_shape));g=g_shape;
x_shape=zeros(x_shape);x_shape(:)=x(1:numel(x_shape));x=x_shape;
return;
end;
%
nerr = 1;
iopt = 2;
[dumvar1,dumvar2,dumvar3,nerr,iopt]=xermsg('SLATEC','BNDSOL',['A ZERO DIAGONAL TERM IS IN THE N BY N UPPER TRIANGULAR ','MATRIX.'],nerr,iopt);
g_shape=zeros(g_shape);g_shape(:)=g(1:numel(g_shape));g=g_shape;
x_shape=zeros(x_shape);x_shape(:)=x(1:numel(x_shape));x=x_shape;
end %subroutine bndsol
%DECK BNFAC