## MC01XD

### The roots of a third order polynomial

[Specification] [Arguments] [Method] [References] [Comments] [Example]

Purpose

```  To compute the roots of the polynomial

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

```
Specification
```      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)

```
Arguments

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.

```
Workspace
```  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.

```
Error Indicator
```  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.

```
Method
```  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.

```
Numerical Aspects
```  The algorithm is numerically stable.

```
```  None
```
Example

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

```
Program Data
``` MC01XD EXAMPLE PROGRAM DATA
-3098110792.0746     4649783048.0746     -2327508384.0000     388574551.8120
```
Program Results
``` 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
```