NAG FL Interface
c02agf (poly_real)

Note: this routine is deprecated. Replaced by c02abf.

Old Code

Real (Kind=nag_wp) :: a(0:n)
Real (Kind=nag_wp) :: z(2,n)
Real (Kind=nag_wp) :: w(2*(n+1))
Logical :: scal
...
scal = .True.
Call c02agf(a, n, scal, z, w, ifail)

New Code

Real (Kind=nag_wp) :: a(0:n)
Real (Kind=nag_wp) :: z(2,n)
Complex (Kind=nag_wp) :: zz(n)
Integer :: itmax, polish
Real (Kind=nag_wp) :: berr(n), cond(n)
Integer :: conv(n)
...
itmax = 30
polish=1
Call c02abf(a, n, itmax, polish, zz, berr, cond, conv, ifail)
! Convert array z from Complex to Real
z(1,1:n) = Real(zz(1:n))
z(2,1:n) = Aimag(zz(1:n))

that argument z is a two-dimensional array of type real(kind=nag_wp) for c02agf, but is a one-dimensional array of type complex(kind=nag_wp) for c02abf. Note also that the roots may be returned in a different order in array z.

Keyword Search:

NAG Library Manual, Mark 28.5

Interfaces: FL CL CPP AD PY MB

NAG FL Interface Introduction

C02 (Zeros) Chapter Contents

C02 (Zeros) Chapter Introduction

c02ag: FL CL CPP AD PY MB

▸▿ Contents

1 Purpose

2 Specification

3 Description

4 References

5 Arguments

6 Error Indicators and Warnings

7 Accuracy

8 Parallelism and Performance

9 Further Comments

▸▿ 10 Example

10.1 Program Text

10.2 Program Data

10.3 Program Results

1 Purpose

c02agf finds all the roots of a real polynomial equation, using a variant of Laguerre's method.

2 Specification

Fortran Interface

Subroutine c02agf (

a, n, scal, z, w, ifail)

Integer, Intent (In)	::	n
Integer, Intent (Inout)	::	ifail
Real (Kind=nag_wp), Intent (In)	::	a(n+1)
Real (Kind=nag_wp), Intent (Out)	::	z(2,n), w(2*(n+1))
Logical, Intent (In)	::	scal

C Header Interface

#include <nag.h>

void	c02agf_ (const double a[], const Integer n, const logical scal, double z[], double w[], Integer *ifail)

The routine may be called by the names c02agf or nagf_zeros_poly_real.

3 Description

c02agf attempts to find all the roots of the

n

th degree real polynomial equation

P (z) = a_{0} z^{n} + a_{1} z^{n - 1} + a_{2} z^{n - 2} + \dots + a_{n - 1} z + a_{n} = 0 .

The roots are located using a modified form of Laguerre's method, originally proposed by Smith (1967).

The method of Laguerre (see Wilkinson (1965)) can be described by the iterative scheme

L (z_{k}) = z_{k + 1} - z_{k} = \frac{- n P (z_{k})}{P^{'} (z_{k}) \pm \sqrt{H (z_{k})}},

where

H (z_{k}) = (n - 1) [(n - 1) {(P^{'} (z_{k}))}^{2} - n P (z_{k}) P^{''} (z_{k})]

and

z_{0}

is specified.

The sign in the denominator is chosen so that the modulus of the Laguerre step at

z_{k}

, viz.

| L (z_{k}) |

, is as small as possible. The method can be shown to be cubically convergent for isolated roots (real or complex) and linearly convergent for multiple roots.

The routine generates a sequence of iterates

z_{1}, z_{2}, z_{3}, \dots

, such that

| P (z_{k + 1}) | < | P (z_{k}) |

and ensures that

z_{k + 1} + L (z_{k + 1})

‘roughly’ lies inside a circular region of radius

| F |

about

z_{k}

known to contain a zero of

P (z)

; that is,

| L (z_{k + 1}) | \leq | F |

, where

F

denotes the Fejér bound (see Marden (1966)) at the point

z_{k}

. Following Smith (1967),

F

is taken to be

\min (B, 1.445 n R)

, where

B

is an upper bound for the magnitude of the smallest zero given by

B = 1.0001 \times \min (\sqrt{n} L (z_{k}), | r_{1} |, {| a_{n} / a_{0} |}^{1 / n}),

r_{1}

is the zero

X

of smaller magnitude of the quadratic equation

\frac{P^{''} (z_{k})}{2 n (n - 1)} X^{2} + \frac{P^{'} (z_{k})}{n} X + \frac{1}{2} P (z_{k}) = 0

and the Cauchy lower bound

R

for the smallest zero is computed (using Newton's Method) as the positive root of the polynomial equation

| a_{0} | z^{n} + | a_{1} | z^{n - 1} + | a_{2} | z^{n - 2} + \dots + | a_{n - 1} | z - | a_{n} | = 0 .

Starting from the origin, successive iterates are generated according to the rule

z_{k + 1} = z_{k} + L (z_{k})

, for

k = 1, 2, 3, \dots

, and

L (z_{k})

is ‘adjusted’ so that

| P (z_{k + 1}) | < | P (z_{k}) |

and

| L (z_{k + 1}) | \leq | F |

. The iterative procedure terminates if

P (z_{k + 1})

is smaller in absolute value than the bound on the rounding error in

P (z_{k + 1})

and the current iterate

z_{p} = z_{k + 1}

is taken to be a zero of

P (z)

(as is its conjugate

{\bar{z}}_{p}

z_{p}

is complex). The deflated polynomial

\tilde{P} (z) = P (z) / (z - z_{p})

of degree

n - 1

z_{p}

is real (

\tilde{P} (z) = P (z) / ((z - z_{p}) (z - {\bar{z}}_{p}))

of degree

n - 2

z_{p}

is complex) is then formed, and the above procedure is repeated on the deflated polynomial until

n < 3

, whereupon the remaining roots are obtained via the ‘standard’ closed formulae for a linear (

n = 1

) or quadratic (

n = 2

) equation.

Note that c02ajf, c02akf and c02alf can be used to obtain the roots of a quadratic, cubic (

n = 3

) and quartic (

n = 4

) polynomial, respectively.

4 References

Marden M (1966) Geometry of polynomials Mathematical Surveys 3 American Mathematical Society, Providence, RI

Smith B T (1967) ZERPOL: a zero finding algorithm for polynomials using Laguerre's method Technical Report Department of Computer Science, University of Toronto, Canada

Thompson K W (1991) Error analysis for polynomial solvers Fortran Journal (Volume 3) 3 10–13

Wilkinson J H (1965) The Algebraic Eigenvalue Problem Oxford University Press, Oxford

5 Arguments

1: $a (n + 1)$ – Real (Kind=nag_wp) array Input: On entry: if a is declared with bounds $(0 : n)$ , $a (i)$ must contain $a_{i}$ (i.e., the coefficient of $z^{n - i}$ ), for $i = 0, 1, \dots, n$ .

Constraint: $a (0) \neq 0.0$ .
2: $n$ – Integer Input: On entry: $n$ , the degree of the polynomial.

Constraint: $n \geq 1$ .
3: $scal$ – Logical Input: On entry: indicates whether or not the polynomial is to be scaled. See Section 9 for advice on when it may be preferable to set $scal = .FALSE.$ and for a description of the scaling strategy.

Suggested value: $scal = .TRUE.$ .
4: $z (2, n)$ – Real (Kind=nag_wp) array Output: On exit: the real and imaginary parts of the roots are stored in $z (1, i)$ and $z (2, i)$ respectively, for $i = 1, 2, \dots, n$ . Complex conjugate pairs of roots are stored in consecutive pairs of elements of z; that is, $z (1, i + 1) = z (1, i)$ and $z (2, i + 1) = - z (2, i)$ .
5: $w (2 \times (n + 1))$ – Real (Kind=nag_wp) array Workspace
6: $ifail$ – Integer Input/Output: On entry: ifail must be set to $0$ , $−1$ or $1$ to set behaviour on detection of an error; these values have no effect when no error is detected.
A value of $0$ causes the printing of an error message and program execution will be halted; otherwise program execution continues. A value of $−1$ means that an error message is printed while a value of $1$ means that it is not.

If halting is not appropriate, the value $−1$ or $1$ is recommended. If message printing is undesirable, then the value $1$ is recommended. Otherwise, the value $0$ is recommended. When the value $- 1$ or $1$ is used it is essential to test the value of ifail on exit.

On exit: $ifail = 0$ unless the routine detects an error or a warning has been flagged (see Section 6).

6 Error Indicators and Warnings

If on entry

ifail = 0

−1

, explanatory error messages are output on the current error message unit (as defined by x04aaf).

Errors or warnings detected by the routine:

$ifail = 1$: On entry, $a (0) = 0.0$ .
Constraint: $a (0) \neq 0.0$ .

On entry, $n = ⟨ value ⟩$ .
Constraint: $n \geq 1$ .

$ifail = 2$: The iterative procedure has failed to converge. This error is very unlikely to occur. If it does, please contact NAG immediately, as some basic assumption for the arithmetic has been violated.

$ifail = 3$: c02agf cannot evaluate $P (z)$ near some of its zeros without overflow. If this message occurs please contact NAG.

c02agf cannot evaluate $P (z)$ near some of its zeros without underflow. If this message occurs please contact NAG.

$ifail = - 99$: An unexpected error has been triggered by this routine. Please contact NAG.
See Section 7 in the Introduction to the NAG Library FL Interface for further information.

$ifail = - 399$: Your licence key may have expired or may not have been installed correctly.
See Section 8 in the Introduction to the NAG Library FL Interface for further information.

$ifail = - 999$: Dynamic memory allocation failed.
See Section 9 in the Introduction to the NAG Library FL Interface for further information.

7 Accuracy

All roots are evaluated as accurately as possible, but because of the inherent nature of the problem complete accuracy cannot be guaranteed. See also Section 10.

8 Parallelism and Performance

Background information to multithreading can be found in the Multithreading documentation.

c02agf is not threaded in any implementation.

9 Further Comments

scal = .TRUE.

, then a scaling factor for the coefficients is chosen as a power of the base

b

of the machine so that the largest coefficient in magnitude approaches

thresh = b^{e_{\max} - p}

. You should note that no scaling is performed if the largest coefficient in magnitude exceeds

thresh

, even if

scal = .TRUE.

. (

b

e_{\max}

and

p

are defined in Chapter X02.)

However, with

scal = .TRUE.

, overflow may be encountered when the input coefficients

a_{0}, a_{1}, a_{2}, \dots, a_{n}

vary widely in magnitude, particularly on those machines for which

b^{(4 p)}

overflows. In such cases, scal should be set to .FALSE. and the coefficients scaled so that the largest coefficient in magnitude does not exceed

b^{(e_{\max} - 2 p)}

Even so, the scaling strategy used by c02agf is sometimes insufficient to avoid overflow and/or underflow conditions. In such cases, you are recommended to scale the independent variable

(z)

so that the disparity between the largest and smallest coefficient in magnitude is reduced. That is, use the routine to locate the zeros of the polynomial

d P (c z)

for some suitable values of

c

and

d

. For example, if the original polynomial was

P (z) = 2^{−100} + 2^{100} z^{20}

, then choosing

c = 2^{−10}

and

d = 2^{100}

, for instance, would yield the scaled polynomial

1 + z^{20}

, which is well-behaved relative to overflow and underflow and has zeros which are

2^{10}

times those of

P (z)

If the routine fails with

ifail = 2

3

, then the real and imaginary parts of any roots obtained before the failure occurred are stored in z in the reverse order in which they were found. Let

n_{R}

denote the number of roots found before the failure occurred. Then

z (1, n)

and

z (2, n)

contain the real and imaginary parts of the first root found,

z (1, n - 1)

and

z (2, n - 1)

contain the real and imaginary parts of the second root found,

\dots, z (1, n - n_{R} + 1)

and

z (2, n - n_{R} + 1)

contain the real and imaginary parts of the

n_{R}

th root found. After the failure has occurred, the remaining

2 \times (n - n_{R})

elements of z contain a large negative number (equal to

−1 / (x02amf () \times \sqrt{2})

10 Example

For this routine two examples are presented. There is a single example program for c02agf, with a main program and the code to solve the two example problems given in the subroutines EX1 and EX2.

Example 1 (EX1)

This example finds the roots of the fifth degree polynomial

z^{5} + 2 z^{4} + 3 z^{3} + 4 z^{2} + 5 z + 6 = 0 .

Example 2 (EX2)

This example solves the same problem as subroutine EX1, but in addition attempts to estimate the accuracy of the computed roots using a perturbation analysis. Further details can be found in Thompson (1991).