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

function [nd,bd,nm1,bm1,nm2,bm2,na,aa,x,yy,m,a,b,c,d,w,y]=cprod(nd,bd,nm1,bm1,nm2,bm2,na,aa,x,yy,m,a,b,c,d,w,y);
persistent crt den ia id iflg j k m1 m2 mm rt y1 y2 ; 

if isempty(rt), rt=0; end;
if isempty(ia), ia=0; end;
if isempty(id), id=0; end;
if isempty(iflg), iflg=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;
%***BEGIN PROLOGUE  CPROD
%***SUBSIDIARY
%***PURPOSE  Subsidiary to BLKTRI
%***LIBRARY   SLATEC
%***TYPE      SINGLE PRECISION (CPROD-S, CPROC-C)
%***AUTHOR  (UNKNOWN)
%***DESCRIPTION
%
% PROD applies a sequence of matrix operations to the vector X and
% stores the result in YY.   (COMPLEX case)
% AA         array containing scalar multipliers of the vector X.
% ND,NM1,NM2 are the lengths of the arrays BD,BM1,BM2 respectively.
% BD,BM1,BM2 are arrays containing roots of certain B polynomials.
% NA         is the length of the array AA.
% X,YY      The matrix operations are applied to X and the result is YY.
% A,B,C      are arrays which contain the tridiagonal matrix.
% M          is the order of the matrix.
% D,W,Y      are working arrays.
% ISGN       determines whether or not a change in sign is made.
%
%***SEE ALSO  BLKTRI
%***ROUTINES CALLED  (NONE)
%***REVISION HISTORY  (YYMMDD)
%   801001  DATE WRITTEN
%   890531  Changed all specific intrinsics to generic.  (WRB)
%   891214  Prologue converted to Version 4.0 format.  (BAB)
%   900402  Added TYPE section.  (WRB)
%***end PROLOGUE  CPROD
%
if isempty(crt), crt=0; end;
if isempty(den), den=0; end;
if isempty(y1), y1=0; end;
if isempty(y2), y2=0; end;
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,[]);
yy_shape=size(yy);yy=reshape(yy,1,[]);
%***FIRST EXECUTABLE STATEMENT  CPROD
for j = 1 : m;
y(j) = complex(x(j),0.);
end; j = fix(m+1);
mm = fix(m - 1);
id = fix(nd);
m1 = fix(nm1);
m2 = fix(nm2);
ia = fix(na);
while( true );
iflg = 0;
if( id>0 )
crt = bd(id);
id = fix(id - 1);
%
% BEGIN SOLUTION TO SYSTEM
%
d(m) = a(m)./(b(m)-crt);
w(m) = y(m)./(b(m)-crt);
for j = 2 : mm;
k = fix(m - j);
den = b(k+1) - crt - 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) - crt - c(1).*d(2);
if( abs(den)==0 )
y(1) =complex(1.,0.);
else;
y(1) =(y(1)-c(1).*w(2))./den;
end;
for j = 2 : m;
y(j) = w(j) - d(j).*y(j-1);
end; j = fix(m+1);
end;
if( m1<=0 )
if( m2<=0 )
if( ia>0 )
rt = aa(ia);
ia = fix(ia - 1);
iflg = 1;
%
% SCALAR MULTIPLICATION
%
for j = 1 : m;
y(j) = rt.*y(j);
end; j = fix(m+1);
end;
if( iflg>0 )
continue;
end;
break;
else;
rt = bm2(m2);
m2 = fix(m2 - 1);
end;
elseif( m2<=0 ) ;
rt = bm1(m1);
m1 = fix(m1 - 1);
elseif( abs(bm1(m1))<=abs(bm2(m2)) ) ;
rt = bm2(m2);
m2 = fix(m2 - 1);
else;
rt = bm1(m1);
m1 = fix(m1 - 1);
end;
y1 =(b(1)-rt).*y(1) + c(1).*y(2);
if( mm>=2 )
%
% MATRIX MULTIPLICATION
%
for j = 2 : mm;
y2 = a(j).*y(j-1) +(b(j)-rt).*y(j) + c(j).*y(j+1);
y(j-1) = y1;
y1 = y2;
end; j = fix(mm+1);
end;
y(m) = a(m).*y(m-1) +(b(m)-rt).*y(m);
y(m-1) = y1;
iflg = 1;
end;
for j = 1 : m;
yy(j) = real(y(j));
end; j = fix(m+1);
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;
yy_shape=zeros(yy_shape);yy_shape(:)=yy(1:numel(yy_shape));yy=yy_shape;
end %subroutine cprod
%DECK CPRODP