[n,a,x]=bksol(n,a,x);

function [n,a,x]=bksol(n,a,x);
persistent j k m nm1 ; 

if isempty(j), j=0; end;
if isempty(k), k=0; end;
if isempty(m), m=0; end;
if isempty(nm1), nm1=0; end;
%***BEGIN PROLOGUE  BKSOL
%***SUBSIDIARY
%***PURPOSE  Subsidiary to BVSUP
%***LIBRARY   SLATEC
%***TYPE      SINGLE PRECISION (BKSOL-S, DBKSOL-D)
%***AUTHOR  Watts, H. A., (SNLA)
%***DESCRIPTION
%
% **********************************************************************
%     Solution of an upper triangular linear system by
%     back-substitution
%
%     The matrix A is assumed to be stored in a linear
%     array proceeding in a row-wise manner. The
%     vector X contains the given constant vector on input
%     and contains the solution on return.
%     The actual diagonal of A is unity while a diagonal
%     scaling matrix is stored there.
% **********************************************************************
%
%***SEE ALSO  BVSUP
%***ROUTINES CALLED  SDOT
%***REVISION HISTORY  (YYMMDD)
%   750601  DATE WRITTEN
%   891214  Prologue converted to Version 4.0 format.  (BAB)
%   900328  Added TYPE section.  (WRB)
%   910722  Updated AUTHOR section.  (ALS)
%***end PROLOGUE  BKSOL
%
a_shape=size(a);a=reshape(a,1,[]);
x_shape=size(x);x=reshape(x,1,[]);
%
%***FIRST EXECUTABLE STATEMENT  BKSOL
m =fix(fix((n.*(n+1))./2));
x(n) = x(n).*a(m);
if( n~=1 )
nm1 = fix(n - 1);
for k = 1 : nm1;
j = fix(n - k);
m = fix(m - k - 1);
x(j) = x(j).*a(m) - sdot(k,a(sub2ind(size(a),max(m+1,1)):end),1,x(sub2ind(size(x),max(j+1,1)):end),1);
end; k = fix(nm1+1);
end;
%
a_shape=zeros(a_shape);a_shape(:)=a(1:numel(a_shape));a=a_shape;
x_shape=zeros(x_shape);x_shape(:)=x(1:numel(x_shape));x=x_shape;
end
%DECK BLKTR1