relibc_openlibm/slatec/cpqr79.f

*DECK CPQR79
      SUBROUTINE CPQR79 (NDEG, COEFF, ROOT, IERR, WORK)
C***BEGIN PROLOGUE  CPQR79
C***PURPOSE  Find the zeros of a polynomial with complex coefficients.
C***LIBRARY   SLATEC
C***CATEGORY  F1A1B
C***TYPE      COMPLEX (RPQR79-S, CPQR79-C)
C***KEYWORDS  COMPLEX POLYNOMIAL, POLYNOMIAL ROOTS, POLYNOMIAL ZEROS
C***AUTHOR  Vandevender, W. H., (SNLA)
C***DESCRIPTION
C
C   Abstract
C       This routine computes all zeros of a polynomial of degree NDEG
C       with complex coefficients by computing the eigenvalues of the
C       companion matrix.
C
C   Description of Parameters
C       The user must dimension all arrays appearing in the call list
C            COEFF(NDEG+1), ROOT(NDEG), WORK(2*NDEG*(NDEG+1))
C
C    --Input--
C      NDEG    degree of polynomial
C
C      COEFF   COMPLEX coefficients in descending order.  i.e.,
C              P(Z)= COEFF(1)*(Z**NDEG) + COEFF(NDEG)*Z + COEFF(NDEG+1)
C
C      WORK    REAL work array of dimension at least 2*NDEG*(NDEG+1)
C
C   --Output--
C      ROOT    COMPLEX vector of roots
C
C      IERR    Output Error Code
C           - Normal Code
C          0  means the roots were computed.
C           - Abnormal Codes
C          1  more than 30 QR iterations on some eigenvalue of the
C             companion matrix
C          2  COEFF(1)=0.0
C          3  NDEG is invalid (less than or equal to 0)
C
C***REFERENCES  (NONE)
C***ROUTINES CALLED  COMQR, XERMSG
C***REVISION HISTORY  (YYMMDD)
C   791201  DATE WRITTEN
C   890531  Changed all specific intrinsics to generic.  (WRB)
C   890531  REVISION DATE from Version 3.2
C   891214  Prologue converted to Version 4.0 format.  (BAB)
C   900315  CALLs to XERROR changed to CALLs to XERMSG.  (THJ)
C   900326  Removed duplicate information from DESCRIPTION section.
C           (WRB)
C   911010  Code reworked and simplified.  (RWC and WRB)
C***END PROLOGUE  CPQR79
      COMPLEX COEFF(*), ROOT(*), SCALE, C
      REAL WORK(*)
      INTEGER NDEG, IERR, K, KHR, KHI, KWR, KWI, KAD, KJ
C***FIRST EXECUTABLE STATEMENT  CPQR79
      IERR = 0
      IF (ABS(COEFF(1)) .EQ. 0.0) THEN
         IERR = 2
         CALL XERMSG ('SLATEC', 'CPQR79',
     +      'LEADING COEFFICIENT IS ZERO.', 2, 1)
         RETURN
      ENDIF
C
      IF (NDEG .LE. 0) THEN
         IERR = 3
         CALL XERMSG ('SLATEC', 'CPQR79', 'DEGREE INVALID.', 3, 1)
         RETURN
      ENDIF
C
      IF (NDEG .EQ. 1) THEN
         ROOT(1) = -COEFF(2)/COEFF(1)
         RETURN
      ENDIF
C
      SCALE = 1.0E0/COEFF(1)
      KHR = 1
      KHI = KHR+NDEG*NDEG
      KWR = KHI+KHI-KHR
      KWI = KWR+NDEG
C
      DO 10 K=1,KWR
         WORK(K) = 0.0E0
   10 CONTINUE
C
      DO 20 K=1,NDEG
         KAD = (K-1)*NDEG+1
         C = SCALE*COEFF(K+1)
         WORK(KAD) = -REAL(C)
         KJ = KHI+KAD-1
         WORK(KJ) = -AIMAG(C)
         IF (K .NE. NDEG) WORK(KAD+K) = 1.0E0
   20 CONTINUE
C
      CALL COMQR (NDEG,NDEG,1,NDEG,WORK(KHR),WORK(KHI),WORK(KWR),
     1   WORK(KWI),IERR)
C
      IF (IERR .NE. 0) THEN
         IERR = 1
         CALL XERMSG ('SLATEC', 'CPQR79',
     +      'NO CONVERGENCE IN 30 QR ITERATIONS.', 1, 1)
         RETURN
      ENDIF
C
      DO 30 K=1,NDEG
         KM1 = K-1
         ROOT(K) = CMPLX(WORK(KWR+KM1),WORK(KWI+KM1))
   30 CONTINUE
      RETURN
      END