! C02AGF Example Program Text ! Mark 24 Release. NAG Copyright 2012. Module c02agfe_mod ! C02AGF Example Program Module: ! Parameters ! .. Implicit None Statement .. Implicit None ! .. Parameters .. Integer, Parameter :: nin = 5, nout = 6 Logical, Parameter :: 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