**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