!   D02AGF Example Program Text
!   Mark 23 Release. NAG Copyright 2011.

    MODULE d02agfe_mod

!      D02AGF Example Program Module:
!             Global Parameters

!      iprint:    set iprint = 1 for output at each Newton iteration.
!      nin:       the input channel number
!      nout:      the output channel number

!      .. Use Statements ..
       USE nag_library, ONLY : nag_wp
!      .. Implicit None Statement ..
       IMPLICIT NONE
!      .. Parameters ..
       INTEGER, PARAMETER                  :: iprint = 0, nin = 5, nout = 6
    END MODULE d02agfe_mod
    MODULE d02agfe_ex1

!      D02AGF Example Program Module:
!             Parameters and User-defined Routines For Example 1

!      n : number of differential equations,
!      n1: number of parameters.

!      .. Use Statements ..
       USE nag_library, ONLY : nag_wp
       USE d02agfe_mod, ONLY : iprint, nout
!      .. Implicit None Statement ..
       IMPLICIT NONE
!      .. Parameters ..
       INTEGER, PARAMETER                  :: n = 2, n1 = 2
    CONTAINS
!dec$ attributes offload : mic :: aux1
       SUBROUTINE aux1(f,y,x,param)

!         .. Implicit None Statement ..
          IMPLICIT NONE
!         .. Scalar Arguments ..
          REAL (KIND=nag_wp), INTENT (IN)     :: x
!         .. Array Arguments ..
          REAL (KIND=nag_wp), INTENT (OUT)    :: f(*)
          REAL (KIND=nag_wp), INTENT (IN)     :: param(*), y(*)
!         .. Executable Statements ..
          f(1) = y(2)
          f(2) = (y(1)**3-y(2))/(2.0_nag_wp*x)
          RETURN
!dec$ attributes offload : mic :: aux1
       END SUBROUTINE aux1
!dec$ attributes offload : mic :: raaux1
       SUBROUTINE raaux1(x0,x1,r,param)

!         .. Implicit None Statement ..
          IMPLICIT NONE
!         .. Scalar Arguments ..
          REAL (KIND=nag_wp), INTENT (OUT)    :: r, x0, x1
!         .. Array Arguments ..
          REAL (KIND=nag_wp), INTENT (IN)     :: param(*)
!         .. Executable Statements ..
          x0 = 0.1_nag_wp
          x1 = 16.0_nag_wp
          r = 16.0_nag_wp
          RETURN
!dec$ attributes offload : mic :: raaux1
       END SUBROUTINE raaux1
!dec$ attributes offload : mic :: bcaux1
       SUBROUTINE bcaux1(g0,g1,param)

!         .. Implicit None Statement ..
          IMPLICIT NONE
!         .. Array Arguments ..
          REAL (KIND=nag_wp), INTENT (OUT)    :: g0(*), g1(*)
          REAL (KIND=nag_wp), INTENT (IN)     :: param(*)
!         .. Local Scalars ..
          REAL (KIND=nag_wp)                  :: z
!         .. Intrinsic Functions ..
          INTRINSIC                              sqrt
!         .. Executable Statements ..
          z = 0.1_nag_wp
          g0(1) = 0.1_nag_wp + param(1)*sqrt(z)*0.1_nag_wp + 0.01_nag_wp*z
          g0(2) = param(1)*0.05_nag_wp/sqrt(z) + 0.01_nag_wp
          g1(1) = 1.0_nag_wp/6.0_nag_wp
          g1(2) = param(2)
          RETURN
!dec$ attributes offload : mic :: bcaux1
       END SUBROUTINE bcaux1
!dec$ attributes offload : mic :: prsol1
       SUBROUTINE prsol1(param,res,n1,err)

!         .. Implicit None Statement ..
          IMPLICIT NONE
!         .. Scalar Arguments ..
          REAL (KIND=nag_wp), INTENT (IN)     :: res
          INTEGER, INTENT (IN)                :: n1
!         .. Array Arguments ..
          REAL (KIND=nag_wp), INTENT (IN)     :: err(n1), param(n1)
!         .. Local Scalars ..
          INTEGER                             :: i
!         .. Executable Statements ..
          IF (iprint/=0) THEN
             WRITE (nout,99999) 'Current parameters   ', (param(i),i=1,n1)
             WRITE (nout,99998) 'Residuals   ', (err(i),i=1,n1)
             WRITE (nout,99998) 'Sum of residuals squared  ', res
             WRITE (nout,*)
          END IF
          RETURN

99999     FORMAT (1X,A,6(E14.6,2X))
99998     FORMAT (1X,A,6(E12.4,1X))
!dec$ attributes offload : mic :: prsol1
       END SUBROUTINE prsol1
    END MODULE d02agfe_ex1
    MODULE d02agfe_ex2

!      D02AGF Example Program Module:
!             Parameters and User-defined Routines For Example 2

!      n : number of differential equations,
!      n1: number of parameters.

!      .. Use Statements ..
       USE nag_library, ONLY : nag_wp
       USE d02agfe_mod, ONLY : iprint, nout
!      .. Implicit None Statement ..
       IMPLICIT NONE
!      .. Parameters ..
       INTEGER, PARAMETER                  :: n = 3, n1 = 3
    CONTAINS
!dec$ attributes offload : mic :: aux2
       SUBROUTINE aux2(f,y,x,param)

!         .. Implicit None Statement ..
          IMPLICIT NONE
!         .. Parameters ..
          REAL (KIND=nag_wp), PARAMETER       :: eps = 2.0E-5_nag_wp
!         .. Scalar Arguments ..
          REAL (KIND=nag_wp), INTENT (IN)     :: x
!         .. Array Arguments ..
          REAL (KIND=nag_wp), INTENT (OUT)    :: f(*)
          REAL (KIND=nag_wp), INTENT (IN)     :: param(*), y(*)
!         .. Local Scalars ..
          REAL (KIND=nag_wp)                  :: c, s
!         .. Intrinsic Functions ..
          INTRINSIC                              cos, sin
!         .. Executable Statements ..
          c = cos(y(3))
          s = sin(y(3))
          f(1) = s/c
          f(2) = -(param(1)*s+eps*y(2)*y(2))/(y(2)*c)
          f(3) = -param(1)/(y(2)*y(2))
          RETURN
!dec$ attributes offload : mic :: aux2
       END SUBROUTINE aux2
!dec$ attributes offload : mic :: raaux2
       SUBROUTINE raaux2(x0,x1,r,param)

!         .. Implicit None Statement ..
          IMPLICIT NONE
!         .. Scalar Arguments ..
          REAL (KIND=nag_wp), INTENT (OUT)    :: r, x0, x1
!         .. Array Arguments ..
          REAL (KIND=nag_wp), INTENT (IN)     :: param(*)
!         .. Executable Statements ..
          x0 = 0.0_nag_wp
          x1 = param(2)
          r = param(2)
          RETURN
!dec$ attributes offload : mic :: raaux2
       END SUBROUTINE raaux2
!dec$ attributes offload : mic :: bcaux2
       SUBROUTINE bcaux2(g0,g1,param)

!         .. Implicit None Statement ..
          IMPLICIT NONE
!         .. Array Arguments ..
          REAL (KIND=nag_wp), INTENT (OUT)    :: g0(*), g1(*)
          REAL (KIND=nag_wp), INTENT (IN)     :: param(*)
!         .. Executable Statements ..
          g0(1) = 0.0E0_nag_wp
          g0(2) = 500.0E0_nag_wp
          g0(3) = 0.5E0_nag_wp
          g1(1) = 0.0E0_nag_wp
          g1(2) = 450.0E0_nag_wp
          g1(3) = param(3)
          RETURN
!dec$ attributes offload : mic :: bcaux2
       END SUBROUTINE bcaux2
!dec$ attributes offload : mic :: prsol2
       SUBROUTINE prsol2(param,res,n1,err)

!         .. Implicit None Statement ..
          IMPLICIT NONE
!         .. Scalar Arguments ..
          REAL (KIND=nag_wp), INTENT (IN)     :: res
          INTEGER, INTENT (IN)                :: n1
!         .. Array Arguments ..
          REAL (KIND=nag_wp), INTENT (IN)     :: err(n1), param(n1)
!         .. Local Scalars ..
          INTEGER                             :: i
!         .. Executable Statements ..
          IF (iprint/=0) THEN
             WRITE (nout,99999) 'Current parameters   ', (param(i),i=1,n1)
             WRITE (nout,99998) 'Residuals   ', (err(i),i=1,n1)
             WRITE (nout,99998) 'Sum of residuals squared  ', res
             WRITE (nout,*)
          END IF
          RETURN

99999     FORMAT (1X,A,6(E14.6,2X))
99998     FORMAT (1X,A,6(E12.4,1X))
!dec$ attributes offload : mic :: prsol2
       END SUBROUTINE prsol2
    END MODULE d02agfe_ex2
    PROGRAM d02agfe

!      D02AGF Example Main Program

!      .. Use Statements ..
       USE d02agfe_mod, ONLY : nout
!      .. Implicit None Statement ..
       IMPLICIT NONE
!      .. Executable Statements ..
       WRITE (nout,*) 'D02AGF Example Program Results'

       CALL ex1

       CALL ex2

    CONTAINS
!dec$ attributes offload : mic :: ex1
       SUBROUTINE ex1

!         .. Use Statements ..
          USE nag_library, ONLY : d02agf
          USE d02agfe_mod, ONLY : nin
          USE d02agfe_ex1, ONLY : aux1, bcaux1, n, n1, nag_wp, prsol1, raaux1
!         .. Implicit None Statement ..
          IMPLICIT NONE
!         .. Local Scalars ..
          REAL (KIND=nag_wp)                  :: h, r, soler, x, x1, xx
          INTEGER                             :: i, ifail, j, m1
!         .. Local Arrays ..
          REAL (KIND=nag_wp), ALLOCATABLE     :: c(:,:), e(:), mat(:,:),       &
                                                 param(:), parerr(:),          &
                                                 wspac1(:), wspac2(:),         &
                                                 wspace(:,:)
          REAL (KIND=nag_wp)                  :: dummy(1,1)
!         .. Intrinsic Functions ..
          INTRINSIC                              real
!         .. Executable Statements ..
!         Skip heading in data file
          READ (nin,*)
!         m1: final solution calculated at m1 points in range including end-points,
          READ (nin,*) m1
          ALLOCATE (c(m1,n),e(n),mat(n,n),param(n),parerr(n),wspac1(n), &
             wspac2(n),wspace(n,9))
!         h:      Step size estimate,          param: initial parameter estimates
!         parerr: Newton iteration tolerances, soler: bound on the local error.
          READ (nin,*) h
          READ (nin,*) param(1:n1)
          READ (nin,*) parerr(1:n1)
          READ (nin,*) soler
          e(1:n) = soler
          WRITE (nout,*)
          WRITE (nout,*)
          WRITE (nout,*) 'Case 1'
          WRITE (nout,*)

!         ifail: behaviour on error exit   
!                =0 for hard exit, =1 for quiet-soft, =-1 for noisy-soft
          ifail = 0
          CALL d02agf(h,e,parerr,param,c,n,n1,m1,aux1,bcaux1,raaux1,prsol1, &
             mat,dummy,wspace,wspac1,wspac2,ifail)

          WRITE (nout,*) 'Final parameters'
          WRITE (nout,99999) (param(i),i=1,n1)
          WRITE (nout,*)
          WRITE (nout,*) 'Final solution'
          WRITE (nout,*) 'X-value     Components of solution'
          CALL raaux1(x,x1,r,param)
          h = (x1-x)/real(m1-1,kind=nag_wp)
          xx = x
          DO i = 1, m1
             WRITE (nout,99998) xx, (c(i,j),j=1,n)
             xx = xx + h
          END DO

          RETURN

99999     FORMAT (1X,3E16.6)
99998     FORMAT (1X,F7.2,3E13.4)
!dec$ attributes offload : mic :: ex1
       END SUBROUTINE ex1
!dec$ attributes offload : mic :: ex2
       SUBROUTINE ex2

!         .. Use Statements ..
          USE nag_library, ONLY : d02agf
          USE d02agfe_mod, ONLY : nin
          USE d02agfe_ex2, ONLY : aux2, bcaux2, n, n1, nag_wp, prsol2, raaux2
!         .. Implicit None Statement ..
          IMPLICIT NONE
!         .. Local Scalars ..
          REAL (KIND=nag_wp)                  :: h, r, soler, x, x1, xx
          INTEGER                             :: i, ifail, j, m1
!         .. Local Arrays ..
          REAL (KIND=nag_wp), ALLOCATABLE     :: c(:,:), e(:), mat(:,:),       &
                                                 param(:), parerr(:),          &
                                                 wspac1(:), wspac2(:),         &
                                                 wspace(:,:)
          REAL (KIND=nag_wp)                  :: dummy(1,1)
!         .. Executable Statements ..
          READ (nin,*)
!         Read in problem parameters
!         m1: final solution calculated at m1 points in range including end-points,
          READ (nin,*) m1
          ALLOCATE (c(m1,n),e(n),mat(n,n),param(n),parerr(n),wspac1(n), &
             wspac2(n),wspace(n,9))
          WRITE (nout,*)
          WRITE (nout,*)
          WRITE (nout,*) 'Case 2'
          WRITE (nout,*)
!         h:      Step size estimate,          param: initial parameter estimates
!         parerr: Newton iteration tolerances, soler: bound on the local error.
          READ (nin,*) h
          READ (nin,*) param(1:n1)
          READ (nin,*) parerr(1:n1)
          READ (nin,*) soler
          e(1:n) = soler

!         ifail: behaviour on error exit   
!                =0 for hard exit, =1 for quiet-soft, =-1 for noisy-soft
          ifail = 0
          CALL d02agf(h,e,parerr,param,c,n,n1,m1,aux2,bcaux2,raaux2,prsol2, &
             mat,dummy,wspace,wspac1,wspac2,ifail)

          WRITE (nout,*) 'Final parameters'
          WRITE (nout,99999) (param(i),i=1,n)
          WRITE (nout,*)
          WRITE (nout,*) 'Final solution'
          WRITE (nout,*) 'X-value     Components of solution'
          CALL raaux2(x,x1,r,param)
          h = (x1-x)/5.0E0_nag_wp
          xx = x
          DO i = 1, 6
             WRITE (nout,99998) xx, (c(i,j),j=1,n)
             xx = xx + h
          END DO

          RETURN

99999     FORMAT (1X,3E16.6)
99998     FORMAT (1X,F7.0,3E13.4)
!dec$ attributes offload : mic :: ex2
       END SUBROUTINE ex2
    END PROGRAM d02agfe
