MC01XD

The roots of a third order polynomial

Purpose

  To compute the roots of the polynomial

      P(t) = ALPHA + BETA*t + GAMMA*t^2 + DELTA*t^3 .

      SUBROUTINE MC01XD( ALPHA, BETA, GAMMA, DELTA, EVR, EVI, EVQ,
     $                   DWORK, LDWORK, INFO )
C     .. Scalar Arguments ..
      INTEGER           INFO, LDWORK
      DOUBLE PRECISION  ALPHA, BETA, DELTA, GAMMA
C     .. Array Arguments ..
      DOUBLE PRECISION  DWORK(*), EVI(3), EVQ(3), EVR(3)

Input/Output Parameters

  ALPHA   (input) DOUBLE PRECISION
  BETA    (input) DOUBLE PRECISION
  GAMMA   (input) DOUBLE PRECISION
  DELTA   (input) DOUBLE PRECISION
          The coefficients of the polynomial P.

  EVR     (output) DOUBLE PRECISION array, DIMENSION at least 3
  EVI     (output) DOUBLE PRECISION array, DIMENSION at least 3
  EVQ     (output) DOUBLE PRECISION array, DIMENSION at least 3
          On output, the kth root of P will be equal to
          (EVR(K) + i*EVI(K))/EVQ(K) if EVQ(K) .NE. ZERO. Note that
          the quotient may over- or underflow. If P has a degree d
          less than 3, then 3-d computed roots will be infinite.
          EVQ(K) >= 0.

  DWORK   DOUBLE PRECISION array, DIMENSION (LDWORK)
          On exit, if LDWORK = -1 on input, then DWORK(1) returns
          the optimal value of LDWORK.

  LDWORK  INTEGER
          The dimension of the array DWORK.  LDWORK >= 42.

          If LDWORK = -1, an optimal workspace query is assumed; the
          routine only calculates the optimal size of the DWORK
          array, returns this value as the first entry of the DWORK
          array, and no error message is issued by XERBLA.

  INFO    INTEGER
          = 0:  successful exit;
          < 0:  if INFO = -i, the i-th argument had an illegal
                value;
          > 0:  if INFO = j, 1 <= j <= 6, an error occurred during
                the call to one of the LAPACK Library routines DGGEV
                or DGEEV (1 <= j <= 3). If INFO < 3, the values
                returned in EVR(K), EVI(K), and EVQ(K) should be
                correct for K = INFO+1,...,3.

  A matrix pencil is built, whose eigenvalues are the roots of the
  given polynomial, and they are computed using the QZ algorithm.
  However, when the ratio between the largest and smallest (in
  magnitude) polynomial coefficients is relatively big, and either
  ALPHA or DELTA has the largest magnitude, then a standard
  eigenproblem is solved using the QR algorithm, and EVQ(I) are set
  to 1, for I = 1,2,3.

  The algorithm is numerically stable.

  None

Program Text

*     MC01XD EXAMPLE PROGRAM TEXT.
*     Copyright (c) 2002-2017 NICONET e.V.
*
*     .. Parameters ..
      INTEGER          NIN, NOUT
      PARAMETER        ( NIN = 5, NOUT = 6 )
      INTEGER          NEV
      PARAMETER        ( NEV = 3 )
      INTEGER          LDWORK
      PARAMETER        ( LDWORK = 42 )
      DOUBLE PRECISION ZERO
      PARAMETER        ( ZERO = 0.0D0 )
*     .. Local Scalars ..
      INTEGER          I, INFO, J
      DOUBLE PRECISION ALPHA, BETA, DELTA, GAMMA
*     .. Local Arrays ..
      DOUBLE PRECISION DWORK(LDWORK), EVI(NEV), EVQ(NEV), EVR(NEV),
     $                 RT(2)
*     .. External Subroutines ..
      EXTERNAL         MC01XD
*     .. Executable Statements ..
*
      WRITE ( NOUT, FMT = 99999 )
*     Skip the heading in the data file and read the data.
      READ ( NIN, FMT = '()' )
      READ ( NIN, FMT = * ) ALPHA, BETA, GAMMA, DELTA
*     Compute the roots of the polynomial.
      CALL MC01XD( ALPHA, BETA, GAMMA, DELTA, EVR, EVI, EVQ,
     $             DWORK, LDWORK, INFO )
*
      IF ( INFO.NE.0 ) THEN
         WRITE ( NOUT, FMT = 99998 ) INFO
      ELSE
         WRITE ( NOUT, FMT = 99995 )
         WRITE ( NOUT, FMT = 99996 ) ( EVR(J), J = 1, 3 )
*
         WRITE ( NOUT, FMT = 99994 )
         WRITE ( NOUT, FMT = 99996 ) ( EVI(J), J = 1, 3 )
*
         WRITE ( NOUT, FMT = 99993 )
         WRITE ( NOUT, FMT = 99996 ) ( EVQ(J), J = 1, 3 )
*
         WRITE ( NOUT, FMT = 99992 )
         DO 20 I = 1, 3
            IF ( EVQ(I).NE.ZERO ) THEN
               RT(1) = EVR(I)/EVQ(I)
               RT(2) = EVI(I)/EVQ(I)
               WRITE ( NOUT, FMT = 99996 ) ( RT(J), J = 1, 2 )
            END IF
   20    CONTINUE
*
      END IF
      STOP
*
99999 FORMAT (' MC01XD EXAMPLE PROGRAM RESULTS',/1X)
99998 FORMAT (' INFO on exit from MC01XD = ',I2)
99996 FORMAT (20(1X,F8.4))
99995 FORMAT (' Real part of the numerators of the roots')
99994 FORMAT (' Imaginary part of the numerators of the roots')
99993 FORMAT (' Denominators of the roots')
99992 FORMAT (' Roots of the polynomial',/1X)
      END

 MC01XD EXAMPLE PROGRAM DATA
   -3098110792.0746     4649783048.0746     -2327508384.0000     388574551.8120

 MC01XD EXAMPLE PROGRAM RESULTS

 Real part of the numerators of the roots
   1.0185   1.1048   0.5911
 Imaginary part of the numerators of the roots
   0.0000   0.0455  -0.0244
 Denominators of the roots
   0.5110   0.5528   0.2958
 Roots of the polynomial

   1.9931   0.0000
   1.9984   0.0823
   1.9984  -0.0823

Return to index