Try out NAG Library functions

Explore NAG maths and stats routines with interactive demos
Function ID
C02AGF
Name
nagf_zeros_poly_real
Description
Zeros of a polynomial with real coefficients
Keywords
Laguerre's method | root-finding
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.
The example data reflects that shown in the "Example" section of the routine documentation. You can change this here to try alternative inputs. The formatting will need to be kept as it is here, otherwise the program is likely to fail to run correctly.

Please note that incompatible data will however cause the example output to display an error message. These error messages are fully explained in the Routine document
!   C02AGF Example Program Text
!   Mark 26 Release. NAG Copyright 2016.

    Module c02agfe_mod

!     C02AGF Example Program Module:
!            Parameters

!     .. Implicit None Statement ..
      Implicit None
!     .. Accessibility Statements ..
      Private
!     .. Parameters ..
      Integer, Parameter, Public       :: nin = 5, nout = 6
      Logical, Parameter, Public       :: scal = .True.
    End Module c02agfe_mod
    Program c02agfe

!     C02AGF Example Main Program

!     .. Use Statements ..
      Use c02agfe_mod, Only: nout
!     .. Implicit None Statement ..
      Implicit None
!     .. Executable Statements ..
      Write (nout,*) 'C02AGF Example Program Results'

      Call ex1

      Call ex2

    Contains
      Subroutine ex1

!       .. Use Statements ..
        Use nag_library, Only: c02agf, nag_wp
        Use c02agfe_mod, Only: nin, scal
!       .. Local Scalars ..
        Real (Kind=nag_wp)             :: zi, zr
        Integer                        :: i, ifail, n, nroot
!       .. Local Arrays ..
        Real (Kind=nag_wp), Allocatable :: a(:), w(:), z(:,:)
!       .. Intrinsic Procedures ..
        Intrinsic                      :: abs
!       .. Executable Statements ..
        Write (nout,*)
        Write (nout,*)
        Write (nout,*) 'Example 1'

!       Skip heading in data file
        Read (nin,*)
        Read (nin,*)
        Read (nin,*)

        Read (nin,*) n
        Allocate (a(0:n),w(2*(n+1)),z(2,n))

        Read (nin,*)(a(i),i=0,n)

        Write (nout,*)
        Write (nout,99999) 'Degree of polynomial = ', n

        ifail = 0
        Call c02agf(a,n,scal,z,w,ifail)

        Write (nout,99998) 'Computed roots of polynomial'

        nroot = 1

        Do While (nroot<=n)

          zr = z(1,nroot)
          zi = z(2,nroot)
          If (zi==0.0E0_nag_wp) Then
            Write (nout,99997) 'z = ', zr
            nroot = nroot + 1
          Else
            Write (nout,99997) 'z = ', zr, ' +/- ', abs(zi), '*i'
            nroot = nroot + 2
          End If

        End Do

99999   Format (/,1X,A,I4)
99998   Format (/,1X,A,/)
99997   Format (1X,A,1P,E12.4,A,E12.4,A)
      End Subroutine ex1
      Subroutine ex2

!       .. Use Statements ..
        Use nag_library, Only: a02abf, c02agf, nag_wp, x02ajf, x02alf
        Use c02agfe_mod, Only: nin, scal
!       .. Local Scalars ..
        Real (Kind=nag_wp)             :: deltac, deltai, di, eps, epsbar, f,  &
                                          r1, r2, r3, rmax
        Integer                        :: i, ifail, j, jmin, n
!       .. Local Arrays ..
        Real (Kind=nag_wp), Allocatable :: a(:), abar(:), r(:), w(:), z(:,:),  &
                                          zbar(:,:)
        Integer, Allocatable           :: m(:)
!       .. Intrinsic Procedures ..
        Intrinsic                      :: abs, max, min
!       .. Executable Statements ..
        Write (nout,*)
        Write (nout,*)
        Write (nout,*) 'Example 2'

!       Skip heading in data file
        Read (nin,*)
        Read (nin,*)

        Read (nin,*) n
        Allocate (a(0:n),abar(0:n),r(n),w(2*(n+1)),z(2,n),zbar(2,n),m(n))

!       Read in the coefficients of the original polynomial.

        Read (nin,*)(a(i),i=0,n)

!       Compute the roots of the original polynomial.

        ifail = 0
        Call c02agf(a,n,scal,z,w,ifail)

!       Form the coefficients of the perturbed polynomial.

        eps = x02ajf()
        epsbar = 3.0_nag_wp*eps

        Do i = 0, n

          If (a(i)/=0.0_nag_wp) Then
            f = 1.0_nag_wp + epsbar
            epsbar = -epsbar
            abar(i) = f*a(i)
          Else
            abar(i) = 0.0E0_nag_wp
          End If

        End Do

!       Compute the roots of the perturbed polynomial.

        ifail = 0
        Call c02agf(abar,n,scal,zbar,w,ifail)

!       Perform error analysis.

!       Initialize markers to 0 (unmarked).

        m(1:n) = 0

        rmax = x02alf()

!       Loop over all unperturbed roots (stored in Z).

        Do i = 1, n
          deltai = rmax
          r1 = a02abf(z(1,i),z(2,i))

!         Loop over all perturbed roots (stored in ZBAR).

          Do j = 1, n

!           Compare the current unperturbed root to all unmarked
!           perturbed roots.

            If (m(j)==0) Then
              r2 = a02abf(zbar(1,j),zbar(2,j))
              deltac = abs(r1-r2)

              If (deltac<deltai) Then
                deltai = deltac
                jmin = j
              End If

            End If

          End Do

!         Mark the selected perturbed root.

          m(jmin) = 1

!         Compute the relative error.

          If (r1/=0.0E0_nag_wp) Then
            r3 = a02abf(zbar(1,jmin),zbar(2,jmin))
            di = min(r1,r3)
            r(i) = max(deltai/max(di,deltai/rmax),eps)
          Else
            r(i) = 0.0_nag_wp
          End If

        End Do

        Write (nout,*)
        Write (nout,99999) 'Degree of polynomial = ', n
        Write (nout,*)
        Write (nout,*) 'Computed roots of polynomial   ', '  Error estimates'
        Write (nout,*) '                            ',                         &
          '   (machine-dependent)'
        Write (nout,*)

        Do i = 1, n
          Write (nout,99998) 'z = ', z(1,i), z(2,i), '*i', r(i)
        End Do

99999   Format (1X,A,I4)
99998   Format (1X,A,1P,E12.4,Sp,E12.4,A,5X,Ss,E9.1)
      End Subroutine ex2
    End Program c02agfe
The NAG Library
The world’s largest collection of robust, documented, tested and maintained numerical algorithms.
Learn more here or contact us for purchasing information