[ndeg,coeff,root,ierr,work]=cpqr79(ndeg,coeff,root,ierr,work);

function [ndeg,coeff,root,ierr,work]=cpqr79(ndeg,coeff,root,ierr,work);
persistent c k kad khi khr kj km1 kwi kwr scalemlv ; 

if isempty(km1), km1=0; end;
%***BEGIN PROLOGUE  CPQR79
%***PURPOSE  Find the zeros of a polynomial with complex coefficients.
%***LIBRARY   SLATEC
%***CATEGORY  F1A1B
%***TYPE      COMPLEX (RPQR79-S, CPQR79-C)
%***KEYWORDS  COMPLEX POLYNOMIAL, POLYNOMIAL ROOTS, POLYNOMIAL ZEROS
%***AUTHOR  Vandevender, W. H., (SNLA)
%***DESCRIPTION
%
%   Abstract
%       This routine computes all zeros of a polynomial of degree NDEG
%       with complex coefficients by computing the eigenvalues of the
%       companion matrix.
%
%   Description of Parameters
%       The user must dimension all arrays appearing in the call list
%            COEFF(NDEG+1), ROOT(NDEG), WORK(2*NDEG*(NDEG+1))
%
%    --Input--
%      NDEG    degree of polynomial
%
%      COEFF   COMPLEX coefficients in descending order.  i.e.,
%              P(Z)= COEFF(1)*(Z**NDEG) + COEFF(NDEG)*Z + COEFF(NDEG+1)
%
%      WORK    REAL work array of dimension at least 2*NDEG*(NDEG+1)
%
%   --Output--
%      ROOT    COMPLEX vector of roots
%
%      IERR    Output Error Code
%           - Normal Code
%          0  means the roots were computed.
%           - Abnormal Codes
%          1  more than 30 QR iterations on some eigenvalue of the
%             companion matrix
%          2  COEFF(1)=0.0
%          3  NDEG is invalid (less than or equal to 0)
%
%***REFERENCES  (NONE)
%***ROUTINES CALLED  COMQR, XERMSG
%***REVISION HISTORY  (YYMMDD)
%   791201  DATE WRITTEN
%   890531  Changed all specific intrinsics to generic.  (WRB)
%   890531  REVISION DATE from Version 3.2
%   891214  Prologue converted to Version 4.0 format.  (BAB)
%   900315  CALLs to XERROR changed to CALLs to XERMSG.  (THJ)
%   900326  Removed duplicate information from DESCRIPTION section.
%           (WRB)
%   911010  Code reworked and simplified.  (RWC and WRB)
%***end PROLOGUE  CPQR79
coeff_shape=size(coeff);coeff=reshape(coeff,1,[]);
root_shape=size(root);root=reshape(root,1,[]);
if isempty(scalemlv), scalemlv=0; end;
if isempty(c), c=0; end;
work_shape=size(work);work=reshape(work,1,[]);
if isempty(k), k=0; end;
if isempty(khr), khr=0; end;
if isempty(khi), khi=0; end;
if isempty(kwr), kwr=0; end;
if isempty(kwi), kwi=0; end;
if isempty(kad), kad=0; end;
if isempty(kj), kj=0; end;
%***FIRST EXECUTABLE STATEMENT  CPQR79
ierr = 0;
if( abs(coeff(1))==0.0 )
ierr = 2;
xermsg('SLATEC','CPQR79','LEADING COEFFICIENT IS ZERO.',2,1);
coeff_shape=zeros(coeff_shape);coeff_shape(:)=coeff(1:numel(coeff_shape));coeff=coeff_shape;
root_shape=zeros(root_shape);root_shape(:)=root(1:numel(root_shape));root=root_shape;
work_shape=zeros(work_shape);work_shape(:)=work(1:numel(work_shape));work=work_shape;
return;
end;
%
if( ndeg<=0 )
ierr = 3;
xermsg('SLATEC','CPQR79','DEGREE INVALID.',3,1);
coeff_shape=zeros(coeff_shape);coeff_shape(:)=coeff(1:numel(coeff_shape));coeff=coeff_shape;
root_shape=zeros(root_shape);root_shape(:)=root(1:numel(root_shape));root=root_shape;
work_shape=zeros(work_shape);work_shape(:)=work(1:numel(work_shape));work=work_shape;
return;
end;
%
if( ndeg==1 )
root(1) = -coeff(2)./coeff(1);
coeff_shape=zeros(coeff_shape);coeff_shape(:)=coeff(1:numel(coeff_shape));coeff=coeff_shape;
root_shape=zeros(root_shape);root_shape(:)=root(1:numel(root_shape));root=root_shape;
work_shape=zeros(work_shape);work_shape(:)=work(1:numel(work_shape));work=work_shape;
return;
end;
%
scalemlv = 1.0e0./coeff(1);
khr = 1;
khi = fix(khr + ndeg.*ndeg);
kwr = fix(khi + khi - khr);
kwi = fix(kwr + ndeg);
%
for k = 1 : kwr;
work(k) = 0.0e0;
end; k = fix(kwr+1);
%
for k = 1 : ndeg;
kad =fix((k-1).*ndeg + 1);
c = scalemlv.*coeff(k+1);
work(kad) = -real(c);
kj = fix(khi + kad - 1);
work(kj) = -imag(c);
if( k~=ndeg )
work(kad+k) = 1.0e0;
end;
end; k = fix(ndeg+1);
%
ndeg_orig=ndeg; ndeg_orig=ndeg;    [ndeg,dumvar2,dumvar3,dumvar4,dumvar5,dumvar6,dumvar7,dumvar8,ierr]=comqr(ndeg,ndeg,1,ndeg,work(sub2ind(size(work),max(khr,1)):end),work(sub2ind(size(work),max(khi,1)):end),work(sub2ind(size(work),max(kwr,1)):end),work(sub2ind(size(work),max(kwi,1)):end),ierr);    ndeg(dumvar4~=ndeg_orig)=dumvar4(dumvar4~=ndeg_orig); ndeg(dumvar2~=ndeg_orig)=dumvar2(dumvar2~=ndeg_orig);   dumvar5i=find((work(sub2ind(size(work),max(khr,1)):end))~=(dumvar5));dumvar6i=find((work(sub2ind(size(work),max(khi,1)):end))~=(dumvar6));dumvar7i=find((work(sub2ind(size(work),max(kwr,1)):end))~=(dumvar7));dumvar8i=find((work(sub2ind(size(work),max(kwi,1)):end))~=(dumvar8));   work(khr-1+dumvar5i)=dumvar5(dumvar5i); work(khi-1+dumvar6i)=dumvar6(dumvar6i); work(kwr-1+dumvar7i)=dumvar7(dumvar7i); work(kwi-1+dumvar8i)=dumvar8(dumvar8i); 
%
if( ierr~=0 )
ierr = 1;
xermsg('SLATEC','CPQR79','NO CONVERGENCE IN 30 QR ITERATIONS.',1,1);
coeff_shape=zeros(coeff_shape);coeff_shape(:)=coeff(1:numel(coeff_shape));coeff=coeff_shape;
root_shape=zeros(root_shape);root_shape(:)=root(1:numel(root_shape));root=root_shape;
work_shape=zeros(work_shape);work_shape(:)=work(1:numel(work_shape));work=work_shape;
return;
end;
%
for k = 1 : ndeg;
km1 = fix(k - 1);
root(k) = complex(work(kwr+km1),work(kwi+km1));
end; k = fix(ndeg+1);
coeff_shape=zeros(coeff_shape);coeff_shape(:)=coeff(1:numel(coeff_shape));coeff=coeff_shape;
root_shape=zeros(root_shape);root_shape(:)=root(1:numel(root_shape));root=root_shape;
work_shape=zeros(work_shape);work_shape(:)=work(1:numel(work_shape));work=work_shape;
end
%DECK CPROC