NAG Library Manual, Mark 30.3
Interfaces:  FL   CL   CPP   AD 

NAG AD Library Introduction
Example description
!   E05JC_A1W_F Example Program Text
!   Mark 30.3 Release. NAG Copyright 2024.

    Module e05jc_a1w_fe_mod

!     E05JC_A1W_F Example Program Module:
!            Parameters and User-defined Routines

!     .. Use Statements ..
      Use nagad_library, Only: exp, nagad_a1w_w_rtype, Assignment (=),         &
                               Operator (**), Operator (+), Operator (*),      &
                               Operator (-), Operator (/)
      Use nag_library, Only: nag_wp, x04baf
!     .. Implicit None Statement ..
      Implicit None
!     .. Accessibility Statements ..
      Private
      Public                           :: monit, objfun
!     .. Parameters ..
      Integer, Parameter, Public       :: lcomm = 100, nin = 5, ninopt = 7,    &
                                          nout = 6
!     .. Local Scalars ..
      Logical, Public, Save            :: plot
    Contains
      Subroutine outbox(boxl,boxu)

!       Displays edges of box with bounds BOXL and BOXU in format suitable
!       for plotting.

!       .. Array Arguments ..
        Real (Kind=nag_wp), Intent (In) :: boxl(2), boxu(2)
!       .. Local Scalars ..
        Character (80)                 :: rec
!       .. Executable Statements ..
        Write (rec,99999) boxl(1), boxl(2)
        Call x04baf(nout,rec)
        Write (rec,99999) boxl(1), boxu(2)
        Call x04baf(nout,rec)
        Write (rec,'()')
        Call x04baf(nout,rec)
        Write (rec,99999) boxl(1), boxl(2)
        Call x04baf(nout,rec)
        Write (rec,99999) boxu(1), boxl(2)
        Call x04baf(nout,rec)
        Write (rec,'()')
        Call x04baf(nout,rec)
        Write (rec,99999) boxl(1), boxu(2)
        Call x04baf(nout,rec)
        Write (rec,99999) boxu(1), boxu(2)
        Call x04baf(nout,rec)
        Write (rec,'()')
        Call x04baf(nout,rec)
        Write (rec,99999) boxu(1), boxl(2)
        Call x04baf(nout,rec)
        Write (rec,99999) boxu(1), boxu(2)
        Call x04baf(nout,rec)
        Write (rec,'()')
        Call x04baf(nout,rec)

        Return

99999   Format (F20.15,1X,F20.15)
      End Subroutine outbox
      Subroutine objfun(ad_handle,n,x,f,nstate,iuser,ruser,inform)

!       Routine to evaluate E05JB_A1W_F objective function.

!       .. Use Statements ..
        Use iso_c_binding, Only: c_ptr
!       .. Scalar Arguments ..
        Type (c_ptr), Intent (Inout)   :: ad_handle
        Type (nagad_a1w_w_rtype), Intent (Out) :: f
        Integer, Intent (Out)          :: inform
        Integer, Intent (In)           :: n, nstate
!       .. Array Arguments ..
        Type (nagad_a1w_w_rtype), Intent (Inout) :: ruser(*)
        Type (nagad_a1w_w_rtype), Intent (In) :: x(n)
        Integer, Intent (Inout)        :: iuser(*)
!       .. Local Scalars ..
        Type (nagad_a1w_w_rtype)       :: x1, x2
        Character (80)                 :: rec
!       .. Executable Statements ..

!       This is a two-dimensional objective function.
!       As an example of using the inform mechanism,
!       terminate if any other problem size is supplied.

        If (n/=2) Then
          inform = -1
        Else
          inform = 0

          If (inform>=0) Then

!           If INFORM >= 0 then we're prepared to evaluate OBJFUN
!           at the current X

            If (nstate==1) Then

!             This is the first call to OBJFUN

              Write (rec,'()')
              Call x04baf(nout,rec)
              Write (rec,99999)
              Call x04baf(nout,rec)
            End If

            x1 = x(1)
            x2 = x(2)
            f = ruser(1)*(ruser(2)-x1)**2*exp(-(x1**2)-(x2+ruser(3))**2) -     &
              ruser(4)*(x1/5.0E0_nag_wp-x1**3-x2**5)*exp(-x1**2-x2**2) -       &
              ruser(5)*exp(-(x1+ruser(6))**2-x2**2)
          End If

        End If

        Return

99999   Format (1X,'(OBJFUN was just called for the first time)')
      End Subroutine objfun
      Subroutine monit(n,ncall,xbest,icount,ninit,list,numpts,initpt,nbaskt,   &
        xbaskt,boxl,boxu,nstate,iuser,ruser,inform)

!       Monitoring routine for E05JB_A1W_F.

!       .. Scalar Arguments ..
        Integer, Intent (Out)          :: inform
        Integer, Intent (In)           :: n, nbaskt, ncall, ninit, nstate
!       .. Array Arguments ..
        Type (nagad_a1w_w_rtype), Intent (In) :: boxl(n), boxu(n),             &
                                          list(n,ninit), xbaskt(n,nbaskt),     &
                                          xbest(n)
        Type (nagad_a1w_w_rtype), Intent (Inout) :: ruser(*)
        Integer, Intent (In)           :: icount(6), initpt(n), numpts(n)
        Integer, Intent (Inout)        :: iuser(*)
!       .. Local Scalars ..
        Integer                        :: i
        Character (80)                 :: rec
!       .. Executable Statements ..
        inform = 0

        If (inform>=0) Then

!         We are going to allow the iterations to continue.

          If (nstate==0 .Or. nstate==1) Then

!           When NSTATE==1, MONIT is called for the first time. When
!           NSTATE==0, MONIT is called for the first AND last time.
!           Display a welcome message

            Write (rec,'()')
            Call x04baf(nout,rec)
            Write (rec,99999)
            Call x04baf(nout,rec)
            Write (rec,'()')
            Call x04baf(nout,rec)

            Write (rec,99998)
            Call x04baf(nout,rec)

            Do i = 1, n
              Write (rec,99997)
              Call x04baf(nout,rec)
              Write (rec,99996) i
              Call x04baf(nout,rec)
              Write (rec,99995) numpts(i)
              Call x04baf(nout,rec)
              Write (rec,99994)
              Call x04baf(nout,rec)
              Write (rec,99993) list(i,1:numpts(i))%value
              Call x04baf(nout,rec)
              Write (rec,99992) initpt(i)
              Call x04baf(nout,rec)
            End Do

            If (plot .And. (n==2)) Then
              Write (rec,99991)
              Call x04baf(nout,rec)
              Write (rec,'()')
              Call x04baf(nout,rec)
            End If

          End If

          If (plot .And. (n==2)) Then

!           Display the coordinates of the edges of the current search
!           box

            Call outbox(boxl%value,boxu%value)

          End If

          If (nstate<=0) Then

!           MONIT is called for the last time

            If (plot .And. (n==2)) Then
              Write (rec,99990)
              Call x04baf(nout,rec)
              Write (rec,'()')
              Call x04baf(nout,rec)
            End If

            Write (rec,99989) icount(1)
            Call x04baf(nout,rec)
            Write (rec,99988) 10*((ncall+5)/10)
            Call x04baf(nout,rec)
            Write (rec,99987)
            Call x04baf(nout,rec)
            Write (rec,99986) 10*((icount(2)+5)/10)
            Call x04baf(nout,rec)
            Write (rec,99985) icount(3)
            Call x04baf(nout,rec)
            Write (rec,99984) icount(4)
            Call x04baf(nout,rec)
            Write (rec,99983) icount(5)
            Call x04baf(nout,rec)
            Write (rec,99982) icount(6)
            Call x04baf(nout,rec)
            Write (rec,99981) nbaskt
            Call x04baf(nout,rec)
            Write (rec,99980)
            Call x04baf(nout,rec)

            Do i = 1, n
              Write (rec,99979) i, xbaskt(i,1:nbaskt)%value
              Call x04baf(nout,rec)
            End Do

            Write (rec,99978)
            Call x04baf(nout,rec)
            Write (rec,99977) xbest(1:n)%value
            Call x04baf(nout,rec)

            Write (rec,'()')
            Call x04baf(nout,rec)
            Write (rec,99976)
            Call x04baf(nout,rec)
            Write (rec,'()')
            Call x04baf(nout,rec)
          End If

        End If

        Return

99999   Format (1X,'*** Begin monitoring information ***')
99998   Format (1X,'Values controlling initial splitting of a box:')
99997   Format (1X,'**')
99996   Format (1X,'In dimension ',I5)
99995   Format (1X,'Extent of initialization list in this dimension =',I5)
99994   Format (1X,'Initialization points in this dimension:')
99993   Format (1X,'LIST(I,1:NUMPTS(I)) =',(6F9.5))
99992   Format (1X,'Initial point in this dimension: LIST(I,',I5,')')
99991   Format (1X,'<Begin displaying search boxes>')
99990   Format (1X,'<End displaying search boxes>')
99989   Format (1X,'Total sub-boxes =',I5)
99988   Format (1X,'Total function evaluations (rounded to nearest 10) =',I5)
99987   Format (1X,'Total function evaluations used in local search (rounded')
99986   Format (3X,'to nearest 10) =',I5)
99985   Format (1X,'Total points used in local search =',I5)
99984   Format (1X,'Total sweeps through levels =',I5)
99983   Format (1X,'Total splits by init. list =',I5)
99982   Format (1X,'Lowest level with nonsplit boxes =',I5)
99981   Format (1X,'Number of candidate minima in the "shopping basket','" =', &
          I5)
99980   Format (1X,'Shopping basket:')
99979   Format (1X,'XBASKT(',I3,',:) =',(6F9.5))
99978   Format (1X,'Best point:')
99977   Format (1X,'XBEST =',(6F9.5))
99976   Format (1X,'*** End monitoring information ***')
      End Subroutine monit
    End Module e05jc_a1w_fe_mod
    Program e05jc_a1w_fe

!     E05JC_A1W_F Example Main Program

!     This program demonstrates the use of routines to set and get
!     values of optional parameters associated with E05JB_A1W_F

!     .. Use Statements ..
      Use e05jc_a1w_fe_mod, Only: lcomm, monit, nin, ninopt, nout, objfun,     &
                                  plot
      Use iso_c_binding, Only: c_ptr
      Use nagad_library, Only: e05ja_a1w_f, e05jb_a1w_f, e05jc_a1w_f,          &
                               e05jd_a1w_f, e05je_a1w_f, e05jf_a1w_f,          &
                               e05jg_a1w_f, e05jh_a1w_f, e05jj_a1w_f,          &
                               e05jk_a1w_f, e05jl_a1w_f,                       &
                               nagad_a1w_get_derivative,                       &
                               nagad_a1w_inc_derivative,                       &
                               nagad_a1w_ir_interpret_adjoint_sparse,          &
                               nagad_a1w_ir_register_variable,                 &
                               nagad_a1w_ir_remove, nagad_a1w_w_rtype,         &
                               x10aa_a1w_f, x10ab_a1w_f, x10za_a1w_f,          &
                               Assignment (=)
      Use nag_library, Only: nag_wp, x04acf, x04baf, x04caf
!     .. Implicit None Statement ..
      Implicit None
!     .. Parameters ..
      Integer, Parameter               :: n = 2
      Character (*), Parameter         :: fname = 'e05jc_a1w_fe.opt'
!     .. Local Scalars ..
      Type (c_ptr)                     :: ad_handle
      Type (nagad_a1w_w_rtype)         :: loctol, obj, sqtol
      Integer                          :: i, ibdchk, ibound, ifail, iinit, j,  &
                                          loclim, mode, n_r, sdlist, stclim
      Character (3)                    :: lcsrch
      Character (80)                   :: rec
!     .. Local Arrays ..
      Type (nagad_a1w_w_rtype)         :: bl(n), bu(n), comm(lcomm), ruser(6), &
                                          x(n)
      Type (nagad_a1w_w_rtype), Allocatable :: list(:,:)
      Real (Kind=nag_wp)               :: br(n), x_r(6)
      Real (Kind=nag_wp), Allocatable  :: rlist(:,:)
      Integer                          :: initpt(n), iuser(1), numpts(n)
!     .. Intrinsic Procedures ..
      Intrinsic                        :: count, sqrt, trim
!     .. Executable Statements ..
      Write (rec,99992) 'E05JC_A1W_F Example Program Results'

!     Create AD tape
      Call x10za_a1w_f

!     Create AD configuration data object
      ifail = 0
      Call x10aa_a1w_f(ad_handle,ifail)

      Call x04baf(nout,rec)

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

      Read (nin,*) sdlist
      Allocate (list(n,sdlist))
      Allocate (rlist(n,sdlist))

      Read (nin,*) ibound

      If (ibound==0) Then

!       Read in the whole of each bound

        Read (nin,*)(br(i),i=1,n)
        bl = br
        Read (nin,*)(br(i),i=1,n)
        bu = br
      Else If (ibound==3) Then

!       Bounds are uniform: read in only the first entry of each

        Read (nin,*) br(1)
        bl(1) = br(1)
        Read (nin,*) br(1)
        bu(1) = br(1)
      End If

      Read (nin,*) iinit

      If (iinit==3) Then

!       User is specifying the initialization list

        Read (nin,*)(numpts(i),i=1,n)
        Read (nin,*)((rlist(i,j),j=1,numpts(i)),i=1,n)
        list = rlist
        Read (nin,*)(initpt(i),i=1,n)
      End If

!     PLOT determines whether MONIT displays information on
!     the current search box:

      Read (nin,*) plot

      ruser(1) = 3.0_nag_wp
      ruser(2) = 1.0_nag_wp
      ruser(3) = 1.0_nag_wp
      ruser(4) = 10.0_nag_wp
      ruser(5) = 1.0_nag_wp/3.0_nag_wp
      ruser(6) = 1.0_nag_wp


!     The first argument to E05JA_A1W_F is a legacy argument and has no
!     significance.

      ifail = 0
      Call e05ja_a1w_f(0,comm,lcomm,ifail)

!     Open the options file for reading

      mode = 0

      ifail = 0
      Call x04acf(ninopt,fname,mode,ifail)

!     Use E05JC_A1W_F to read some options from the options file

      ifail = 0
      Call e05jc_a1w_f(ad_handle,ninopt,comm,lcomm,ifail)

      Write (rec,'()')
      Call x04baf(nout,rec)

!     Use E05JK_A1W_F to find the value of the integer-valued option
!     'Local Searches Limit'

      ifail = 0
      Call e05jk_a1w_f('Local Searches Limit',loclim,comm,lcomm,ifail)

      Write (rec,99999) loclim
      Call x04baf(nout,rec)

!     Compute the number of free variables, then use E05JF_A1W_F to set the value
!     of the integer-valued option 'Static Limit'

      n_r = count(bl(1:n)%value/=bu(1:n)%value)
      stclim = 4*n_r

      ifail = 0
      Call e05jf_a1w_f(ad_handle,'Static Limit',stclim,comm,lcomm,ifail)

!     Use E05JH_A1W_F to determine whether the real-valued option
!     'Infinite Bound Size' has been set by us (in which case
!     E05JH_A1W_F returns 1) or whether it holds its default value
!     (E05JH_A1W_F returns 0)

      ifail = 0
      Call e05jh_a1w_f('Infinite Bound Size',ibdchk,comm,lcomm,ifail)

      If (ibdchk==1) Then
        Write (rec,99998)
        Call x04baf(nout,rec)
      Else If (ibdchk==0) Then
        Write (rec,99997)
        Call x04baf(nout,rec)
      End If

!     Use E05JL_A1W_F/E05JG_A1W_F to set the real-valued option
!     'Local Searches Tolerance' to the square root of its default

      ifail = 0
      Call e05jl_a1w_f(ad_handle,'Local Searches Tolerance',loctol,comm,lcomm, &
        ifail)

      sqtol = sqrt(loctol%value)
      ifail = 0
      Call e05jg_a1w_f(ad_handle,'Local Searches Tolerance',sqtol,comm,lcomm,  &
        ifail)

!     Use E05JL_A1W_F to get the new value of 'Local Searches Tolerance'

      ifail = 0
      Call e05jl_a1w_f(ad_handle,'Local Searches Tolerance',loctol,comm,lcomm, &
        ifail)

      Write (rec,99996) loctol%value
      Call x04baf(nout,rec)

!     Use E05JD_A1W_F to set the option 'Minimize' (which is the default)

      ifail = 0
      Call e05jd_a1w_f(ad_handle,'Minimize',comm,lcomm,ifail)

!     Use E05JE_A1W_F to set the option 'Local Searches' to 'On' (also
!     the default)

      lcsrch = 'On'

      ifail = 0
      Call e05je_a1w_f(ad_handle,'Local Searches',lcsrch,comm,lcomm,ifail)

!     Get that value of 'Local Searches' using E05JJ_A1W_F

      ifail = 0
      Call e05jj_a1w_f('Local Searches',lcsrch,comm,lcomm,ifail)

      Write (rec,99995) trim(lcsrch)
      Call x04baf(nout,rec)

!     Register variables to differentiate w.r.t.
      Call nagad_a1w_ir_register_variable(ruser)

!     Solve the problem.

      ifail = 0
      Call e05jb_a1w_f(ad_handle,n,objfun,ibound,iinit,bl,bu,sdlist,list,      &
        numpts,initpt,monit,x,obj,comm,lcomm,iuser,ruser,ifail)

      Write (rec,'()')
      Call x04baf(nout,rec)
      Write (rec,99994) obj%value
      Call x04baf(nout,rec)
      Write (rec,99993)(x(i)%value,i=1,n)
      Call x04baf(nout,rec)

      Write (nout,*)
      Write (nout,*) ' Derivatives calculated: First order adjoints'
      Write (nout,*) ' Computational mode    : algorithmic'

      Call nagad_a1w_inc_derivative(obj,1.0_nag_wp)
      ifail = 0
      Call nagad_a1w_ir_interpret_adjoint_sparse(ifail)
      x_r(1:6) = nagad_a1w_get_derivative(ruser(1:6))
      Write (nout,*)
      Call x04caf('General',' ',6,1,x_r,6,' dobj/druser',ifail)
      Write (nout,*)

!     Remove computational data object and tape
      Call x10ab_a1w_f(ad_handle,ifail)
      Call nagad_a1w_ir_remove

99999 Format (1X,'Option "Local Searches Limit" has the value ',I6,'.')
99998 Format (1X,'Option "Infinite Bound Size" has been set by us.')
99997 Format (1X,'Option "Infinite Bound Size" holds its default value.')
99996 Format (1X,'Option "Local Searches Tolerance" has the value ',E13.5,'.')
99995 Format (1X,'Option "Local Searches" has the value "',A,'".')
99994 Format (1X,'Final objective value =',F11.5)
99993 Format (1X,'Global optimum X =',2F9.5)
99992 Format (1X,A)
    End Program e05jc_a1w_fe