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

NAG AD Library Introduction
Example description
!   E05UC_A1W_F Example Program Text
!   Mark 30.2 Release. NAG Copyright 2024.

    Module e05uc_a1w_fe_mod

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

!     .. Use Statements ..
      Use nagad_library, Only: cos, nagad_a1w_w_rtype, sin, sqrt,              &
                               Operator (**), Assignment (=), Operator (>=),   &
                               Operator (+), Operator (-), Operator (*)
      Use nag_library, Only: nag_wp
!     .. Implicit None Statement ..
      Implicit None
!     .. Accessibility Statements ..
      Private
      Public                           :: mystart, schwefel_confun,            &
                                          schwefel_obj
    Contains
      Subroutine schwefel_obj(ad_handle,mode,n,x,objf,objgrd,nstate,iuser,     &
        ruser)

!       .. 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) :: objf
        Integer, Intent (Inout)        :: mode
        Integer, Intent (In)           :: n, nstate
!       .. Array Arguments ..
        Type (nagad_a1w_w_rtype), Intent (Inout) :: objgrd(n), ruser(*)
        Type (nagad_a1w_w_rtype), Intent (In) :: x(n)
        Integer, Intent (Inout)        :: iuser(*)
!       .. Local Scalars ..
        Integer                        :: i
!       .. Executable Statements ..
        If (nstate==1) Then
!         This is the first call.
!         Take any special action here if desired.
          Continue
        End If
        If (mode==0 .Or. mode==2) Then
!         Evaluate the objective function.
          objf = 0.0_nag_wp
          Do i = 1, n
            If (x(i)>=0.0_nag_wp) Then
              objf = objf + x(i)*sin(ruser(1)*sqrt(x(i)))
            Else
              objf = objf + x(i)*sin(ruser(1)*sqrt(-x(i)))
            End If
          End Do
        End If

        If (mode==1 .Or. mode==2) Then
!         Calculate the gradient of the objective function.
          Do i = 1, n
            If (x(i)>=0.0_nag_wp) Then
              objgrd(i) = sin(ruser(1)*sqrt(x(i))) +                           &
                0.5_nag_wp*ruser(1)*sqrt(x(i))*cos(ruser(1)*sqrt(x(i)))
            Else
              objgrd(i) = sin(ruser(1)*sqrt(-x(i))) +                          &
                0.5_nag_wp*ruser(1)*sqrt(-x(i))*cos(ruser(1)*sqrt(-x(i)))
            End If
          End Do
        End If

        Return
      End Subroutine schwefel_obj
      Subroutine schwefel_confun(ad_handle,mode,ncnln,n,ldcjsl,needc,x,c,cjsl, &
        nstate,iuser,ruser)

!       .. Use Statements ..
        Use iso_c_binding, Only: c_ptr
!       .. Scalar Arguments ..
        Type (c_ptr), Intent (Inout)   :: ad_handle
        Integer, Intent (In)           :: ldcjsl, n, ncnln, nstate
        Integer, Intent (Inout)        :: mode
!       .. Array Arguments ..
        Type (nagad_a1w_w_rtype), Intent (Out) :: c(ncnln)
        Type (nagad_a1w_w_rtype), Intent (Inout) :: cjsl(ldcjsl,n), ruser(*)
        Type (nagad_a1w_w_rtype), Intent (In) :: x(n)
        Integer, Intent (Inout)        :: iuser(*)
        Integer, Intent (In)           :: needc(ncnln)
!       .. Local Scalars ..
        Type (nagad_a1w_w_rtype)       :: t1, t2
        Integer                        :: k
        Logical                        :: evalc, evalcjsl
!       .. Executable Statements ..
        If (nstate==1) Then
!         This is the first call.
!         Take any special action here if desired.
          Continue
        End If

!       mode: what is required - constraints, derivatives, or both?
        evalc = (mode==0 .Or. mode==2)
        evalcjsl = (mode==1 .Or. mode==2)

loop_constraints: Do k = 1, ncnln

          If (needc(k)<=0) Then
            Cycle loop_constraints
          End If

          If (evalc) Then
!           Constraint values are required.
            Select Case (k)
            Case (1)
              c(k) = ruser(2)*x(1)**2 - ruser(3)*x(2)**2 + ruser(4)*x(1)*x(2)
            Case (2)
              c(k) = cos((x(1)*ruser(5))**2+(x(2)*ruser(6)))
            Case Default
!             This constraint is not coded (there are only two).
!             Terminate.
              mode = -1
              Exit loop_constraints
            End Select
          End If

          If (evalcjsl) Then
!           Constraint derivatives are required.
            Select Case (k)
            Case (1)
              cjsl(k,1) = 2.0_nag_wp*ruser(2)*x(1) + ruser(4)*x(2)
              cjsl(k,2) = -2.0_nag_wp*ruser(3)*x(2) + ruser(4)*x(1)
            Case (2)
              t1 = x(1)*ruser(5)
              t2 = x(2)*ruser(6)
              cjsl(k,1) = -sin(t1**2+t2)*2.0_nag_wp*t1*ruser(5)
              cjsl(k,2) = -sin(t1**2+t2)*ruser(6)
            End Select
          End If

        End Do loop_constraints

        Return
      End Subroutine schwefel_confun
      Subroutine mystart(ad_handle,npts,quas,n,repeat,bl,bu,iuser,ruser,mode)

!       Sets the initial points.
!       A typical user-defined start procedure.
!       Only nonzero elements of quas need to be specified here.

!       .. Use Statements ..
        Use iso_c_binding, Only: c_ptr
        Use nag_library, Only: g05kgf, g05saf, nag_wp
!       .. Scalar Arguments ..
        Type (c_ptr), Intent (Inout)   :: ad_handle
        Integer, Intent (Inout)        :: mode
        Integer, Intent (In)           :: n, npts
        Logical, Intent (In)           :: repeat
!       .. Array Arguments ..
        Type (nagad_a1w_w_rtype), Intent (In) :: bl(n), bu(n)
        Type (nagad_a1w_w_rtype), Intent (Inout) :: quas(n,npts), ruser(*)
        Integer, Intent (Inout)        :: iuser(*)
!       .. Local Scalars ..
        Real (Kind=nag_wp)             :: rq
        Integer                        :: genid, i, ifail, lstate, subid
!       .. Local Arrays ..
        Real (Kind=nag_wp), Allocatable :: rquas(:)
        Integer, Allocatable           :: state(:)
!       .. Intrinsic Procedures ..
        Intrinsic                      :: real
!       .. Executable Statements ..
        If (repeat) Then
!         Generate a uniform spread of points between bl and bu.
          Do i = 1, npts
            rq = real(i-1,kind=nag_wp)/real(npts-1,kind=nag_wp)
            quas(1:n,i) = bl(1:n) + rq*(bu(1:n)-bl(1:n))
          End Do
        Else
!         Generate a non-repeatable spread of points between bl and bu.
          genid = 2
          subid = 53
          lstate = -1
          Allocate (state(lstate),rquas(n))
          ifail = 0
          Call g05kgf(genid,subid,state,lstate,ifail)
          Deallocate (state)
          Allocate (state(lstate))
          ifail = 0
          Call g05kgf(genid,subid,state,lstate,ifail)
          Do i = 1, npts
            ifail = 0
            Call g05saf(n,state,rquas,ifail)
            quas(1:n,i) = bl(1:n) + (bu(1:n)-bl(1:n))*rquas(1:n)
          End Do
          Deallocate (state,rquas)
        End If
!       Set mode negative to terminate execution for any reason.
        mode = 0
        Return
      End Subroutine mystart
    End Module e05uc_a1w_fe_mod
    Program e05uc_a1w_fe

!     E05UC_A1W_F Example Main Program

!     .. Use Statements ..
      Use e05uc_a1w_fe_mod, Only: mystart, schwefel_confun, schwefel_obj
      Use iso_c_binding, Only: c_ptr
      Use nagad_library, Only: dgemv_a1w, e05uc_a1w_f, e05zk_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, x04caf
!     .. Implicit None Statement ..
      Implicit None
!     .. Parameters ..
      Integer, Parameter               :: n = 2, nclin = 1, ncnln = 2,         &
                                          nin = 5, nout = 6
!     .. Local Scalars ..
      Type (c_ptr)                     :: ad_handle
      Type (nagad_a1w_w_rtype)         :: alpha, beta
      Integer                          :: i, ifail, j, k, l, lda, ldc, ldcjac, &
                                          ldclda, ldobjd, ldr, ldx, liopts,    &
                                          listat, lopts, nb, npts, sda,        &
                                          sdcjac, sdr
      Logical                          :: repeat
!     .. Local Arrays ..
      Type (nagad_a1w_w_rtype), Allocatable :: a(:,:), bl(:), bu(:), c(:,:),   &
                                          cjac(:,:,:), clamda(:,:), objf(:),   &
                                          objgrd(:,:), opts(:), r(:,:,:),      &
                                          work(:), x(:,:)
      Type (nagad_a1w_w_rtype)         :: ruser(6)
      Real (Kind=nag_wp)               :: x_r(6)
      Integer, Allocatable             :: info(:), iopts(:), istate(:,:),      &
                                          iter(:)
      Integer                          :: iuser(1)
!     .. Executable Statements ..
      Write (nout,*) 'E05UC_A1W_F Example Program Results'
      Flush (nout)

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

      Read (nin,*) nb, npts
      Read (nin,*) repeat

      lda = nclin

      If (nclin>0) Then
        sda = n
      Else
        sda = 1
      End If

      ldx = n
      ldobjd = n
      ldc = ncnln
      ldcjac = ncnln

      If (ncnln>0) Then
        sdcjac = n
      Else
        sdcjac = 0
      End If

      ldr = n
      sdr = n
      ldclda = n + nclin + ncnln
      listat = n + nclin + ncnln
      liopts = 740
      lopts = 485
      Allocate (a(lda,sda),bl(n+nclin+ncnln),bu(n+nclin+ncnln),x(ldx,nb),      &
        objf(nb),objgrd(ldobjd,nb),iter(nb),c(ldc,nb),cjac(ldcjac,sdcjac,nb),  &
        r(ldr,sdr,nb),clamda(ldclda,nb),istate(listat,nb),iopts(liopts),       &
        opts(lopts),info(nb),work(nclin))

      bl(1:n) = -500.0_nag_wp
      bl(n+1:n+nclin) = -10000.0_nag_wp
      bl(n+nclin+1:n+nclin+1) = -1.0_nag_wp
      bl(n+nclin+1:n+nclin+ncnln) = -0.9_nag_wp
      bu(1:n) = 500.0_nag_wp
      bu(n+1:n+nclin) = 10.0_nag_wp
      bu(n+nclin+1:n+nclin+1) = 500000.0_nag_wp
      bu(n+nclin+2:n+nclin+ncnln) = 0.9_nag_wp

      a(1,1) = 3.0_nag_wp
      a(1,2) = -2.0_nag_wp

      ruser(1) = 1.0_nag_wp
      ruser(2) = 1.0_nag_wp
      ruser(3) = 1.0_nag_wp
      ruser(4) = 3.0_nag_wp
      ruser(5) = 0.005_nag_wp
      ruser(6) = 0.01_nag_wp

!     Create AD tape
      Call x10za_a1w_f

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

!     Initialize the solver.

      ifail = 0
      Call e05zk_a1w_f(ad_handle,'Initialize = E05UCF',iopts,liopts,opts,      &
        lopts,ifail)

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

!     Solve the problem.

      ifail = -1
      Call e05uc_a1w_f(ad_handle,n,nclin,ncnln,a,lda,bl,bu,schwefel_confun,    &
        schwefel_obj,npts,x,ldx,mystart,repeat,nb,objf,objgrd,ldobjd,iter,c,   &
        ldc,cjac,ldcjac,sdcjac,r,ldr,sdr,clamda,ldclda,istate,listat,iopts,    &
        opts,iuser,ruser,info,ifail)

      Select Case (ifail)
      Case (0)
        l = nb
      Case (8)
        l = info(nb)
        Write (nout,99992) iter(nb)
      Case Default
        Go To 100
      End Select

loop: Do i = 1, l
        Write (nout,99999) i
        Write (nout,99998) info(i)
        Write (nout,99997) 'Varbl'
        Do j = 1, n
          Write (nout,99996) 'V', j, istate(j,i), x(j,i)%value,                &
            clamda(j,i)%value
        End Do
        If (nclin>0) Then
          Write (nout,99997) 'L Con'

!         Below is a call to the adjoint version of DGEMV.
!         This performs the matrix vector multiplication A*X
!         (linear constraint values) and puts the result in
!         the first NCLIN locations of WORK.
          alpha = 1.0_nag_wp
          beta = 0.0_nag_wp
          Call dgemv_a1w('N',nclin,n,alpha,a,lda,x(1,i),1,beta,work,1)

          Do k = n + 1, n + nclin
            j = k - n
            Write (nout,99996) 'L', j, istate(k,i), work(j)%value,             &
              clamda(k,i)%value
          End Do
        End If
        If (ncnln>0) Then
          Write (nout,99997) 'N Con'
          Do k = n + nclin + 1, n + nclin + ncnln
            j = k - n - nclin
            Write (nout,99996) 'N', j, istate(k,i), c(j,i)%value,              &
              clamda(k,i)%value
          End Do
        End If
        Write (nout,99995) objf(i)%value
        Write (nout,99994)
        Write (nout,99993)(clamda(k,i)%value,k=1,n+nclin+ncnln)

        If (l==1) Then
          Exit loop
        End If

        Write (nout,*)

        Write (nout,*)                                                         &
          ' ------------------------------------------------------ '

      End Do loop

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

      Call nagad_a1w_inc_derivative(objf(1),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,' dobjf/druser',ifail)
      Write (nout,*)

100   Continue

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

99999 Format (/,1X,'Solution number',I16)
99998 Format (/,1X,'Local minimization exited with code',I5)
99997 Format (/,1X,A,2X,'Istate',3X,'Value',9X,'Lagr Mult',/)
99996 Format (1X,A,2(1X,I3),4X,F12.4,2X,F12.4)
99995 Format (/,1X,'Final objective value = ',1X,F12.4)
99994 Format (/,1X,'QP multipliers')
99993 Format (1X,F12.4)
99992 Format (1X,I16,' starting points converged')
    End Program e05uc_a1w_fe