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

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

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(mm2), mm2=0; end;
if isempty(gt), gt=zeros(1,3); end;
%***BEGIN PROLOGUE  PROCP
%***SUBSIDIARY
%***PURPOSE  Subsidiary to CBLKTR
%***LIBRARY   SLATEC
%***TYPE      COMPLEX (PRODP-C, PROCP-C)
%***AUTHOR  (UNKNOWN)
%***DESCRIPTION
%
% PROCP applies a sequence of matrix operations to the vector X and
% stores the result in Y (periodic boundary conditions).
%
% 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,W      are working arrays.
% IS         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  PROCP
%
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,[]);
w_shape=size(w);w=reshape(w,1,[]);
if isempty(den), den=0; end;
if isempty(ym), ym=0; end;
if isempty(v), v=0; end;
if isempty(bh), bh=0; end;
if isempty(am), am=0; end;
%***FIRST EXECUTABLE STATEMENT  PROCP
for j = 1 : m;
y(j) = x(j);
w(j) = y(j);
end; j = fix(m+1);
mm = fix(m - 1);
mm2 = fix(m - 2);
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);
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
%
bh = b(m) - rt;
ym = y(m);
den = b(1) - rt;
d(1) = c(1)./den;
u(1) = a(1)./den;
w(1) = y(1)./den;
v = c(m);
if( mm2>=2 )
for j = 2 : mm2;
den = b(j) - rt - a(j).*d(j-1);
d(j) = c(j)./den;
u(j) = -a(j).*u(j-1)./den;
w(j) =(y(j)-a(j).*w(j-1))./den;
bh = bh - v.*u(j-1);
ym = ym - v.*w(j-1);
v = -v.*d(j-1);
end; j = fix(mm2+1);
end;
den = b(m-1) - rt - a(m-1).*d(m-2);
d(m-1) =(c(m-1)-a(m-1).*u(m-2))./den;
w(m-1) =(y(m-1)-a(m-1).*w(m-2))./den;
am = a(m) - v.*d(m-2);
bh = bh - v.*u(m-2);
ym = ym - v.*w(m-2);
den = bh - am.*d(m-1);
if( abs(den)==0 )
w(m) =complex(1.,0.);
else;
w(m) =(ym-am.*w(m-1))./den;
end;
w(m-1) = w(m-1) - d(m-1).*w(m);
for j = 2 : mm;
k = fix(m - j);
w(k) = w(k) - d(k).*w(k+1) - u(k).*w(m);
end; j = fix(mm+1);
if( na<=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;
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;
w_shape=zeros(w_shape);w_shape(:)=w(1:numel(w_shape));w=w_shape;
end %subroutine procp
%DECK PROD