[nd,bd,nm1,bm1,nm2,bm2,na,aa,x,y,m,a,b,c,d,w,u]=prod(nd,bd,nm1,bm1,nm2,bm2,na,aa,x,y,m,a,b,c,d,w,u);

function [nd,bd,nm1,bm1,nm2,bm2,na,aa,x,y,m,a,b,c,d,w,u]=prod(nd,bd,nm1,bm1,nm2,bm2,na,aa,x,y,m,a,b,c,d,w,u);
persistent den gt ia ibr id j k m1 m2 mm rt ; 

if isempty(den), den=0; end;
if isempty(rt), rt=0; end;
if isempty(ia), ia=0; end;
if isempty(ibr), ibr=0; end;
if isempty(id), id=0; end;
if isempty(j), j=0; end;
if isempty(k), k=0; end;
if isempty(m1), m1=0; end;
if isempty(m2), m2=0; end;
if isempty(mm), mm=0; end;
if isempty(gt), gt=zeros(1,3); end;
%***BEGIN PROLOGUE  PROD
%***SUBSIDIARY
%***PURPOSE  Subsidiary to BLKTRI
%***LIBRARY   SLATEC
%***TYPE      SINGLE PRECISION (PROD-S, PROC-C)
%***AUTHOR  (UNKNOWN)
%***DESCRIPTION
%
% PROD applies a sequence of matrix operations to the vector X and
% stores the result in Y.
%
% BD,BM1,BM2 are arrays containing roots of certain B polynomials.
% ND,NM1,NM2 are the lengths of the arrays BD,BM1,BM2 respectively.
% AA         Array containing scalar multipliers of the vector X.
% NA         is the length of the array AA.
% X,Y        The matrix operations are applied to X and the result is Y.
% A,B,C      are arrays which contain the tridiagonal matrix.
% M          is the order of the matrix.
% D,W,U      are working arrays.
% IS         determines whether or not a change in sign is made.
%
%***SEE ALSO  BLKTRI
%***ROUTINES CALLED  (NONE)
%***REVISION HISTORY  (YYMMDD)
%   801001  DATE WRITTEN
%   891214  Prologue converted to Version 4.0 format.  (BAB)
%   900402  Added TYPE section.  (WRB)
%***end PROLOGUE  PROD
%
a_shape=size(a);a=reshape(a,1,[]);
b_shape=size(b);b=reshape(b,1,[]);
c_shape=size(c);c=reshape(c,1,[]);
x_shape=size(x);x=reshape(x,1,[]);
y_shape=size(y);y=reshape(y,1,[]);
d_shape=size(d);d=reshape(d,1,[]);
w_shape=size(w);w=reshape(w,1,[]);
bd_shape=size(bd);bd=reshape(bd,1,[]);
bm1_shape=size(bm1);bm1=reshape(bm1,1,[]);
bm2_shape=size(bm2);bm2=reshape(bm2,1,[]);
aa_shape=size(aa);aa=reshape(aa,1,[]);
u_shape=size(u);u=reshape(u,1,[]);
%***FIRST EXECUTABLE STATEMENT  PROD
for j = 1 : m;
w(j) = x(j);
y(j) = w(j);
end; j = fix(m+1);
mm = fix(m - 1);
id = fix(nd);
ibr = 0;
m1 = fix(nm1);
m2 = fix(nm2);
ia = fix(na);
while( true );
if( ia>0 )
rt = aa(ia);
if( nd==0 )
rt = -rt;
end;
ia = fix(ia - 1);
%
% SCALAR MULTIPLICATION
%
for j = 1 : m;
y(j) = rt.*w(j);
end; j = fix(m+1);
end;
if( id<=0 )
break;
end;
rt = bd(id);
id = fix(id - 1);
if( id==0 )
ibr = 1;
end;
%
% BEGIN SOLUTION TO SYSTEM
%
d(m) = a(m)./(b(m)-rt);
w(m) = y(m)./(b(m)-rt);
for j = 2 : mm;
k = fix(m - j);
den = b(k+1) - rt - c(k+1).*d(k+2);
d(k+1) = a(k+1)./den;
w(k+1) =(y(k+1)-c(k+1).*w(k+2))./den;
end; j = fix(mm+1);
den = b(1) - rt - c(1).*d(2);
w(1) = 1.;
if( den~=0 )
w(1) =(y(1)-c(1).*w(2))./den;
end;
for j = 2 : m;
w(j) = w(j) - d(j).*w(j-1);
end; j = fix(m+1);
if( na<=0 )
gt(:)=0;
while (1);
if(gt(3)==0)
if(gt(2)==0)
if(gt(1)==0)
if( m1>0 )
if( m2>0 )
if( abs(bm1(m1))<=abs(bm2(m2)) )
gt(1)=1;
continue;
end;
end;
if( ibr<=0 )
if( abs(bm1(m1)-bd(id))<abs(bm1(m1)-rt) )
gt(3)=1;
continue;
end;
end;
rt = rt - bm1(m1);
m1 = fix(m1 - 1);
gt(2)=1;
continue;
elseif( m2<=0 ) ;
gt(3)=1;
continue;
end;
end;
gt(1)=0;
if( ibr<=0 )
if( abs(bm2(m2)-bd(id))<abs(bm2(m2)-rt) )
gt(3)=1;
continue;
end;
end;
rt = rt - bm2(m2);
m2 = fix(m2 - 1);
end;
gt(2)=0;
for j = 1 : m;
y(j) = y(j) + rt.*w(j);
end; j = fix(m+1);
break;
end;
gt(3)=0;
for j = 1 : m;
y(j) = w(j);
end; j = fix(m+1);
ibr = 1;
break;
end;
end;
end;
a_shape=zeros(a_shape);a_shape(:)=a(1:numel(a_shape));a=a_shape;
b_shape=zeros(b_shape);b_shape(:)=b(1:numel(b_shape));b=b_shape;
c_shape=zeros(c_shape);c_shape(:)=c(1:numel(c_shape));c=c_shape;
x_shape=zeros(x_shape);x_shape(:)=x(1:numel(x_shape));x=x_shape;
y_shape=zeros(y_shape);y_shape(:)=y(1:numel(y_shape));y=y_shape;
d_shape=zeros(d_shape);d_shape(:)=d(1:numel(d_shape));d=d_shape;
w_shape=zeros(w_shape);w_shape(:)=w(1:numel(w_shape));w=w_shape;
bd_shape=zeros(bd_shape);bd_shape(:)=bd(1:numel(bd_shape));bd=bd_shape;
bm1_shape=zeros(bm1_shape);bm1_shape(:)=bm1(1:numel(bm1_shape));bm1=bm1_shape;
bm2_shape=zeros(bm2_shape);bm2_shape(:)=bm2(1:numel(bm2_shape));bm2=bm2_shape;
aa_shape=zeros(aa_shape);aa_shape(:)=aa(1:numel(aa_shape));aa=aa_shape;
u_shape=zeros(u_shape);u_shape(:)=u(1:numel(u_shape));u=u_shape;
end %subroutine prod
%DECK PRODP