[k,n,l,x,c,b,m]=sossol(k,n,l,x,c,b,m);

function [k,n,l,x,c,b,m]=sossol(k,n,l,x,c,b,m);
persistent j jkm kj km km1 kmm1 kn lk np1 xmax ; 

if isempty(xmax), xmax=0; end;
if isempty(j), j=0; end;
if isempty(jkm), jkm=0; end;
if isempty(kj), kj=0; end;
if isempty(km), km=0; end;
if isempty(km1), km1=0; end;
if isempty(kmm1), kmm1=0; end;
if isempty(kn), kn=0; end;
if isempty(lk), lk=0; end;
if isempty(np1), np1=0; end;
%***BEGIN PROLOGUE  SOSSOL
%***SUBSIDIARY
%***PURPOSE  Subsidiary to SOS
%***LIBRARY   SLATEC
%***TYPE      SINGLE PRECISION (SOSSOL-S, DSOSSL-D)
%***AUTHOR  (UNKNOWN)
%***DESCRIPTION
%
%     SOSSOL solves an upper triangular type of linear system by back
%     substitution.
%
%     The matrix C is upper trapezoidal and stored as a linear array by
%     rows. The equations have been normalized so that the diagonal
%     entries of C are understood to be unity. The off diagonal entries
%     and the elements of the constant right hand side vector B have
%     already been stored as the negatives of the corresponding equation
%     values.
%     with each call to SOSSOL a (K-1) by (K-1) triangular system is
%     resolved. For L greater than K, column L of C is included in the
%     right hand side vector.
%
%***SEE ALSO  SOS
%***ROUTINES CALLED  (NONE)
%***REVISION HISTORY  (YYMMDD)
%   801001  DATE WRITTEN
%   890831  Modified array declarations.  (WRB)
%   891214  Prologue converted to Version 4.0 format.  (BAB)
%   900328  Added TYPE section.  (WRB)
%***end PROLOGUE  SOSSOL
%
%
x_shape=size(x);x=reshape(x,1,[]);
c_shape=size(c);c=reshape(c,1,[]);
b_shape=size(b);b=reshape(b,1,[]);
%
%***FIRST EXECUTABLE STATEMENT  SOSSOL
np1 = fix(n + 1);
km1 = fix(k - 1);
lk = fix(km1);
if( l==k )
lk = fix(k);
end;
kn = fix(m);
%
%
for kj = 1 : km1;
kmm1 = fix(k - kj);
km = fix(kmm1 + 1);
xmax = 0.;
kn = fix(kn - np1 + kmm1);
if( km<=lk )
jkm = fix(kn);
%
for j = km : lk;
jkm = fix(jkm + 1);
xmax = xmax + c(jkm).*x(j);
end; j = fix(lk+1);
end;
%
if( l>k )
jkm = fix(kn + l - kmm1);
xmax = xmax + c(jkm).*x(l);
end;
x(kmm1) = xmax + b(kmm1);
end; kj = fix(km1+1);
%
x_shape=zeros(x_shape);x_shape(:)=x(1:numel(x_shape));x=x_shape;
c_shape=zeros(c_shape);c_shape(:)=c(1:numel(c_shape));c=c_shape;
b_shape=zeros(b_shape);b_shape(:)=b(1:numel(b_shape));b=b_shape;
end
%DECK SPBCO