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

function [nd,bd,nm1,bm1,nm2,bm2,na,aa,x,y,m,a,b,c,d,u,yy]=cprocp(nd,bd,nm1,bm1,nm2,bm2,na,aa,x,y,m,a,b,c,d,u,yy);
persistent am bh crt den ia id iflg j k m1 m2 mm mm2 rt v y1 y2 yh ym ; 

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;
if isempty(mm2), mm2=0; end;
%***BEGIN PROLOGUE  CPROCP
%***SUBSIDIARY
%***PURPOSE  Subsidiary to CBLKTR
%***LIBRARY   SLATEC
%***TYPE      COMPLEX (CPRODP-S, CPROCP-C)
%***AUTHOR  (UNKNOWN)
%***DESCRIPTION
%
% CPROCP 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,U         are work arrays.
% ISGN        determines whether or not a change in sign is made.
%
%***SEE ALSO  CBLKTR
%***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  CPROCP
%
if isempty(v), v=0; end;
if isempty(den), den=0; end;
if isempty(bh), bh=0; end;
if isempty(ym), ym=0; end;
if isempty(am), am=0; end;
if isempty(y1), y1=0; end;
if isempty(y2), y2=0; end;
if isempty(yh), yh=0; end;
if isempty(crt), crt=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,[]);
u_shape=size(u);u=reshape(u,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  CPROCP
for j = 1 : m;
y(j) = x(j);
end; j = fix(m+1);
mm = fix(m - 1);
mm2 = fix(m - 2);
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);
iflg = 1;
%
% BEGIN SOLUTION TO SYSTEM
%
bh = b(m) - crt;
ym = y(m);
den = b(1) - crt;
d(1) = c(1)./den;
u(1) = a(1)./den;
y(1) = y(1)./den;
v = c(m);
if( mm2>=2 )
for j = 2 : mm2;
den = b(j) - crt - a(j).*d(j-1);
d(j) = c(j)./den;
u(j) = -a(j).*u(j-1)./den;
y(j) =(y(j)-a(j).*y(j-1))./den;
bh = bh - v.*u(j-1);
ym = ym - v.*y(j-1);
v = -v.*d(j-1);
end; j = fix(mm2+1);
end;
den = b(m-1) - crt - a(m-1).*d(m-2);
d(m-1) =(c(m-1)-a(m-1).*u(m-2))./den;
y(m-1) =(y(m-1)-a(m-1).*y(m-2))./den;
am = a(m) - v.*d(m-2);
bh = bh - v.*u(m-2);
ym = ym - v.*y(m-2);
den = bh - am.*d(m-1);
if( abs(den)==0 )
y(m) =complex(1.,0.);
else;
y(m) =(ym-am.*y(m-1))./den;
end;
y(m-1) = y(m-1) - d(m-1).*y(m);
for j = 2 : mm;
k = fix(m - j);
y(k) = y(k) - d(k).*y(k+1) - u(k).*y(m);
end; j = fix(mm+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;
%
% MATRIX MULTIPLICATION
%
yh = y(1);
y1 =(b(1)-rt).*y(1) + c(1).*y(2) + a(1).*y(m);
if( mm>=2 )
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) + c(m).*yh;
y(m-1) = y1;
iflg = 1;
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;
u_shape=zeros(u_shape);u_shape(:)=u(1:numel(u_shape));u=u_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 cprocp
%DECK CPROD