! E05UC_A1W_F Example Program Text
! Mark 28.3 Release. NAG Copyright 2022.
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