! E05SBF Example Program Text
! Mark 27.2 Release. NAG Copyright 2021.
Module e05sbfe_smp_mod
! E05SBF Example Program Module:
! Parameters and User-defined Routines
! .. Use Statements ..
Use nag_library, Only: nag_wp
! .. Implicit None Statement ..
Implicit None
! .. Parameters ..
Real (Kind=nag_wp), Parameter :: zero = 0.0_nag_wp
Integer, Parameter :: detail_level = 0, i_liuser = 1
Integer, Parameter :: i_lruser = i_liuser + 1
Integer, Parameter :: liopts = 100, lopts = 100, &
min_lruser = 0, nin = 5, nout = 6, &
report_freq = 100
Integer, Parameter :: i_liuser_req = i_lruser + 1
Integer, Parameter :: i_lruser_req = i_liuser_req + 1
Integer, Parameter :: i_nlabels = i_lruser_req + 1
Integer, Parameter :: i_nlengths = i_nlabels + 1
Integer, Parameter :: i_max_threads = i_nlengths + 1
Integer, Parameter :: i_ndim = i_max_threads + 1
Integer, Parameter :: i_ncon = i_ndim + 1
Integer, Parameter :: i_iuser_start = i_ncon + 1
Integer, Parameter :: i_ruser_start = i_iuser_start + 1
Integer, Parameter :: i_liuser_callback_shared = &
i_ruser_start + 1
Integer, Parameter :: i_liuser_callback_private = &
i_liuser_callback_shared + 1
Integer, Parameter :: i_lruser_callback_shared = &
i_liuser_callback_private + 1
Integer, Parameter :: i_lruser_callback_private = &
i_lruser_callback_shared + 1
Integer, Parameter :: i_liuser_callback_scratch = &
i_lruser_callback_private + 1
Integer, Parameter :: i_lruser_callback_scratch = &
i_liuser_callback_scratch + 1
Integer, Parameter :: i_label_iuser_callback_shared = &
i_lruser_callback_scratch + 1
Integer, Parameter :: i_label_iuser_callback_private = &
i_label_iuser_callback_shared + 1
Integer, Parameter :: i_label_ruser_callback_shared = &
i_label_iuser_callback_private + 1
Integer, Parameter :: i_label_ruser_callback_private = &
i_label_ruser_callback_shared + 1
Integer, Parameter :: i_label_iuser_scratch = &
i_label_ruser_callback_private + 1
Integer, Parameter :: i_label_ruser_scratch = &
i_label_iuser_scratch + 1
Integer, Parameter :: min_liuser = i_label_ruser_scratch + &
1
Contains
Subroutine initialize_callback_arrays(ndim,ncon,iuser,liuser,ruser, &
lruser,ifail)
! Subroutine to ensure sufficient space in user arrays for labels and to set
! default user array values.
! .. Use Statements ..
Use nag_library, Only: x06acf
! .. Implicit None Statement ..
Implicit None
! .. Scalar Arguments ..
Integer, Intent (Out) :: ifail
Integer, Intent (In) :: liuser, lruser, ncon, ndim
! .. Array Arguments ..
Real (Kind=nag_wp), Intent (Inout) :: ruser(lruser)
Integer, Intent (Inout) :: iuser(liuser)
! .. Local Scalars ..
Integer :: max_threads, nlabels, nlengths
! .. Executable Statements ..
ifail = 0
max_threads = x06acf()
nlabels = (ncon+1)*(max_threads+1)
nlengths = 4*(ncon+1)
iuser(3:liuser) = 0
ruser(1:lruser) = 0.0E0_nag_wp
! Set fixed data
iuser(i_liuser) = liuser
iuser(i_lruser) = lruser
iuser(i_ndim) = ndim
iuser(i_ncon) = ncon
iuser(i_max_threads) = max_threads
iuser(i_liuser_req) = min_liuser + nlabels + nlengths + 1
iuser(i_lruser_req) = min_lruser + 2
! Set references to length references.
iuser(i_liuser_callback_shared) = min_liuser + 1
iuser(i_liuser_callback_private) = iuser(i_liuser_callback_shared) + &
ncon + 1
iuser(i_lruser_callback_shared) = iuser(i_liuser_callback_private) + &
ncon + 1
iuser(i_lruser_callback_private) = iuser(i_lruser_callback_shared) + &
ncon + 1
! Set references to workspace references.
! Require ncon+1 'shared' workpsace references.
iuser(i_label_iuser_callback_shared) = iuser(i_lruser_callback_private &
) + ncon + 1
iuser(i_label_ruser_callback_shared) = iuser( &
i_label_iuser_callback_shared) + ncon + 1
! Require max_threads*(ncon+1) 'private' workspace references.
iuser(i_label_iuser_callback_private) = iuser( &
i_label_ruser_callback_shared) + (ncon+1)
iuser(i_label_ruser_callback_private) = iuser( &
i_label_iuser_callback_private) + (ncon+1)*max_threads
! Require max_threads 'scratch' workspace references
iuser(i_label_iuser_scratch) = iuser(i_label_iuser_callback_private) + &
(ncon+1)*max_threads
iuser(i_label_ruser_scratch) = iuser(i_label_iuser_scratch) + &
max_threads
! Set references to first available workspace to be after last label
iuser(i_iuser_start) = iuser(i_label_ruser_scratch) + max_threads + 1
iuser(i_ruser_start) = min_lruser + 1
Return
End Subroutine initialize_callback_arrays
Subroutine calculate_workspace(ndim,ncon,max_threads,min_liuser, &
min_lruser,fileid,liuser_calc,lruser_calc,ifail)
! Subroutine to calculate required sizes of iuser and ruser from
! information in data file for individual objective and constraint
! functions.
! _shared indicates workspace that is available to all threads.
! _private indicates workspace that is private to one thread,
! and maintained between function calls
! _scratch indicates workspace that is private to each thread
! and not maintained between function calls
! .. Implicit None Statement ..
Implicit None
! .. Scalar Arguments ..
Integer, Intent (In) :: fileid, max_threads, min_liuser, &
min_lruser, ncon, ndim
Integer, Intent (Out) :: ifail, liuser_calc, lruser_calc
! .. Local Scalars ..
Integer :: a, astat, b, c, file_funid, funid, &
liuser_private, liuser_scratch, &
liuser_scratch_max, liuser_shared, &
lruser_private, lruser_scratch, &
lruser_scratch_max, lruser_shared
Character (20) :: buffer
! .. Intrinsic Procedures ..
Intrinsic :: max
! .. Executable Statements ..
ifail = 0
! Initialise to the required space for labels and storage of
! parameters in iuser and ruser
liuser_calc = min_liuser + (ncon+1)*(max_threads+6) + 2*max_threads + &
3
lruser_calc = min_lruser + 2
! Load details of required workspace.
! Rewind data file, then skip lines in data file until the
! correct function or constraint data is found.
Rewind (fileid)
liuser_scratch_max = 0
lruser_scratch_max = 0
Do funid = 0, ncon
loop_data_set: Do
Read (fileid,'(A20)',Iostat=astat) buffer
If (astat<0) Then
! End of file. Assume no data requried.
Exit loop_data_set
Else If (funid==0 .And. buffer(1:6)=='OBJFUN') Then
! Heading for objective function detected.
Exit loop_data_set
Else If (buffer(1:6)=='CONFUN') Then
! Heading for constraint detected.
Read (buffer(7:20),*,Iostat=astat) file_funid
If (funid==file_funid) Then
! Correct constraint number detected.
Exit loop_data_set
End If
End If
End Do loop_data_set
If (astat<0) Then
liuser_shared = 0
lruser_shared = 0
liuser_private = 0
lruser_private = 0
liuser_scratch = 0
lruser_scratch = 0
Else
! Read required memory parameters.
! Here this is stored as a quadratic in ndim.
Read (fileid,*) a, b, c
liuser_shared = (a*ndim+b)*ndim + c
Read (fileid,*) a, b, c
liuser_private = (a*ndim+b)*ndim + c
Read (fileid,*) a, b, c
liuser_scratch = (a*ndim+b)*ndim + c
Read (fileid,*) a, b, c
lruser_shared = (a*ndim+b)*ndim + c
Read (fileid,*) a, b, c
lruser_private = (a*ndim+b)*ndim + c
Read (fileid,*) a, b, c
lruser_scratch = (a*ndim+b)*ndim + c
End If
! Each thread requires only 1 length of scratch workspace.
liuser_scratch_max = max(liuser_scratch_max,liuser_scratch)
lruser_scratch_max = max(lruser_scratch_max,lruser_scratch)
liuser_calc = liuser_calc + max_threads*(liuser_private) + &
liuser_shared
lruser_calc = lruser_calc + max_threads*(lruser_private) + &
lruser_shared
End Do
liuser_calc = liuser_calc + max_threads*liuser_scratch_max
lruser_calc = lruser_calc + max_threads*lruser_scratch_max
Return
End Subroutine calculate_workspace
Subroutine load_callback_data(funid,fileid,iuser,ruser,ifail)
! Subroutine to load user data from the data file
! for individual objective and constraint functions.
! This should only be called when NSTATE = 2
! funid = 0 if called in OBJFUN.
! funid = constraint number if called in confun.
! _shared indicates workspace that is available to all threads.
! _private indicates workspace that is private to one thread,
! and maintained between function calls
! _scratch indicates workspace that is private to each thread
! and not maintained between function calls
! .. Implicit None Statement ..
Implicit None
! .. Scalar Arguments ..
Integer, Intent (In) :: fileid, funid
Integer, Intent (Out) :: ifail
! .. Array Arguments ..
Real (Kind=nag_wp), Intent (Inout) :: ruser(*)
Integer, Intent (Inout) :: iuser(*)
! .. Local Scalars ..
Integer :: a, astat, b, c, file_funid, k, &
labi_private, labi_shared, &
labr_private, labr_shared, &
lab_labi_private, lab_labi_shared, &
lab_labr_private, lab_labr_shared, &
liuser_private, liuser_scratch, &
liuser_shared, lruser_private, &
lruser_scratch, lruser_shared, &
max_threads, ncon, ndim
Character (20) :: buffer
! .. Executable Statements ..
ifail = 0
ndim = iuser(i_ndim)
ncon = iuser(i_ncon)
max_threads = iuser(i_max_threads)
! Load details of required workspace.
! Rewind data file, then skip lines in data file until the
! correct function or constraint data is found.
Rewind (fileid)
loop_data_set: Do
Read (fileid,'(A20)',Iostat=astat) buffer
If (astat<0) Then
! End of file. Assume no data requried.
Exit loop_data_set
Else If (funid==0 .And. buffer(1:6)=='OBJFUN') Then
! Heading for objective function detected.
Exit loop_data_set
Else If (buffer(1:6)=='CONFUN') Then
! Heading for constraint detected.
Read (buffer(7:20),*,Iostat=astat) file_funid
If (funid==file_funid) Then
! Correct constraint number detected.
Exit loop_data_set
End If
End If
End Do loop_data_set
If (astat<0) Then
liuser_shared = 0
lruser_shared = 0
liuser_private = 0
lruser_private = 0
liuser_scratch = 0
lruser_scratch = 0
Else
! Read required memory parameters.
! Here this is stored as a quadratic in ndim.
Read (fileid,*) a, b, c
liuser_shared = (a*ndim+b)*ndim + c
Read (fileid,*) a, b, c
liuser_private = (a*ndim+b)*ndim + c
Read (fileid,*) a, b, c
liuser_scratch = (a*ndim+b)*ndim + c
Read (fileid,*) a, b, c
lruser_shared = (a*ndim+b)*ndim + c
Read (fileid,*) a, b, c
lruser_private = (a*ndim+b)*ndim + c
Read (fileid,*) a, b, c
lruser_scratch = (a*ndim+b)*ndim + c
End If
! Each objective or constraint may have an assigned amount of
! 'shared' and 'private' workspace
iuser(iuser(i_liuser_callback_shared)+funid) = liuser_shared
iuser(iuser(i_lruser_callback_shared)+funid) = lruser_shared
iuser(iuser(i_liuser_callback_private)+funid) = liuser_private
iuser(iuser(i_lruser_callback_private)+funid) = lruser_private
! Each thread requires only 1 length of scratch workspace.
If (iuser(i_liuser_callback_scratch)<liuser_scratch) Then
iuser(i_liuser_callback_scratch) = liuser_scratch
End If
If (iuser(i_lruser_callback_scratch)<lruser_scratch) Then
iuser(i_lruser_callback_scratch) = lruser_scratch
End If
! Set up labels for shared workspace to the next available spaces.
lab_labi_shared = iuser(i_label_iuser_callback_shared) + funid
iuser(lab_labi_shared) = iuser(i_iuser_start)
lab_labr_shared = iuser(i_label_ruser_callback_shared) + funid
iuser(lab_labr_shared) = iuser(i_ruser_start)
! Set up labels for private data for thread 0 (Master)
lab_labi_private = iuser(i_label_iuser_callback_private) + &
max_threads*funid
labi_private = iuser(lab_labi_shared) + liuser_shared
iuser(lab_labi_private) = labi_private
lab_labr_private = iuser(i_label_ruser_callback_private) + &
max_threads*funid
labr_private = iuser(lab_labr_shared) + lruser_shared
iuser(lab_labr_private) = labr_private
! Set i_iuser_start and i_ruser_start to be after shared and private
! space requried for individual function.
iuser(i_iuser_start) = labi_private + max_threads*liuser_private
iuser(i_ruser_start) = labr_private + max_threads*lruser_private
! Adjust scratch workspace labels to point after all
! required private and shared workspace
iuser(iuser(i_label_iuser_scratch)) = iuser(i_iuser_start) + 1
iuser(iuser(i_label_ruser_scratch)) = iuser(i_ruser_start) + 1
! Ensure there is sufficient space allocated in iuser and ruser.
! Note that as this is definitely in a serial region,
! no critical section is required.
iuser(i_liuser_req) = iuser(iuser(i_label_iuser_scratch)) + &
iuser(i_liuser_callback_scratch)*max_threads
iuser(i_lruser_req) = iuser(iuser(i_label_ruser_scratch)) + &
iuser(i_lruser_callback_scratch)*max_threads
If (iuser(i_liuser)<iuser(i_liuser_req) .Or. &
iuser(i_lruser)<iuser(i_lruser_req)) Then
ifail = -2
Write (nout,99999) buffer(1:6)
If (funid==0) Then
Write (nout,99998) ncon, max_threads
Else
Write (nout,99998) funid, max_threads
End If
Write (nout,99997) iuser(i_liuser), iuser(i_liuser_req), &
iuser(i_lruser), iuser(i_lruser_req)
Write (nout,99996)
Go To 100
End If
labi_shared = iuser(lab_labi_shared)
labr_shared = iuser(lab_labr_shared)
! Read in any shared data requried for iuser, if any.
Do k = 1, liuser_shared
Read (fileid,*) iuser(labi_shared+k-1)
End Do
! Read in any shared data required for ruser, if any.
Do k = 1, lruser_shared
Read (fileid,*) ruser(labr_shared+k-1)
End Do
100 Continue
Return
99999 Format (1X,'** Info: IUSER or RUSER is insufficient for ',A6, &
' when **')
99998 Format (1X,'** NCON >= ',I10,' and MAX_THREADS = ',I10,'.**')
99997 Format (1X,'** LIUSER = ',I10,'. Require LIUSER > ',I10,'.**',/,1X, &
'** LRUSER = ',I10,'. Require LRUSER > ',I10,'.**')
99996 Format (1X, &
'** The example program will now try to allocate sufficient space.**'&
)
End Subroutine load_callback_data
Subroutine objfun_schwefel_objective(ndim,x,f,w)
! subroutine to calculate the objective using provided workspace
! .. Use Statements ..
Use nag_library, Only: ddot
! .. Implicit None Statement ..
Implicit None
! .. Scalar Arguments ..
Real (Kind=nag_wp), Intent (Out) :: f
Integer, Intent (In) :: ndim
! .. Array Arguments ..
Real (Kind=nag_wp), Intent (Out) :: w(ndim)
Real (Kind=nag_wp), Intent (In) :: x(ndim)
! .. Intrinsic Procedures ..
Intrinsic :: abs, sin, sqrt
! .. Executable Statements ..
w = abs(x)
w = sqrt(w)
w = sin(w)
f = ddot(ndim,x,1,w,1)
Return
End Subroutine objfun_schwefel_objective
Subroutine objfun_schwefel_gradient(ndim,x,g,w1,w2)
! subroutine to calculate the gradient of the objective function
! with provided workspace
! .. Use Statements ..
Use nag_library, Only: dscal
! .. Implicit None Statement ..
Implicit None
! .. Scalar Arguments ..
Integer, Intent (In) :: ndim
! .. Array Arguments ..
Real (Kind=nag_wp), Intent (Out) :: g(ndim), w1(ndim), w2(ndim)
Real (Kind=nag_wp), Intent (In) :: x(ndim)
! .. Intrinsic Procedures ..
Intrinsic :: abs, cos, sin, sqrt
! .. Executable Statements ..
Continue
w1 = sqrt(abs(x))
w2 = cos(w1)
g = sin(w1)
Call dscal(ndim,0.5E0_nag_wp,w1,1)
g = g + w1*w2
Return
End Subroutine objfun_schwefel_gradient
Subroutine objfun_schwefel(mode,ndim,x,objf,vecout,nstate,iuser,ruser)
! Objfun routine for the Schwefel function for E05SBF_SMP.
! .. Use Statements ..
Use nag_library, Only: x06adf
! .. Implicit None Statement ..
Implicit None
! .. Scalar Arguments ..
Real (Kind=nag_wp), Intent (Inout) :: objf
Integer, Intent (Inout) :: mode
Integer, Intent (In) :: ndim, nstate
! .. Array Arguments ..
Real (Kind=nag_wp), Intent (Inout) :: ruser(*), vecout(ndim)
Real (Kind=nag_wp), Intent (In) :: x(ndim)
Integer, Intent (Inout) :: iuser(*)
! .. Local Scalars ..
Integer :: ifail, labi_private, labi_scratch, &
labi_shared, labr_private, &
labr_scratch, labr_shared, &
lab_labi_private, &
lab_labi_private_master, &
lab_labi_scratch, &
lab_labi_scratch_master, &
lab_labr_private, &
lab_labr_private_master, &
lab_labr_scratch, &
lab_labr_scratch_master, threadno
Logical :: evalobjf, evalobjg
! .. Intrinsic Procedures ..
Intrinsic :: real
! .. Executable Statements ..
! Get the current thread number.
threadno = x06adf()
! Test NSTATE to indicate what stage of computation has been reached.
Select Case (nstate)
Case (2)
! OBJFUN is called for the very first time.
! This occurs in serial on thread 0 only.
Call load_callback_data(0,nin,iuser,ruser,ifail)
If (ifail/=0) Then
mode = ifail
End If
! Perform any one time only calculations required.
labr_shared = iuser(iuser(i_label_ruser_callback_shared))
ruser(labr_shared) = real(ndim,nag_wp)*ruser(labr_shared)
Case (3)
! OBJFUN has been called for the first time on a new thread.
! This occurs at the start of the main parallel region.
! This will not occur on thread 0.
! In this example we take the opportunity to set the labels for this
! threads unique section of RUSER workspace in IUSER.
! This could also be done when NSTATE=2
lab_labi_private_master = iuser(i_label_iuser_callback_private)
lab_labi_private = lab_labi_private_master + threadno
labi_private = iuser(lab_labi_private_master) + &
threadno*(iuser(iuser(i_liuser_callback_private)))
iuser(lab_labi_private) = labi_private
lab_labr_private_master = iuser(i_label_ruser_callback_private)
lab_labr_private = lab_labr_private_master + threadno
labr_private = iuser(lab_labr_private_master) + &
threadno*(iuser(iuser(i_lruser_callback_private)))
iuser(lab_labr_private) = labr_private
! Set private scratch labels.
lab_labi_scratch_master = iuser(i_label_iuser_scratch)
lab_labi_scratch = lab_labi_scratch_master + threadno
labi_scratch = iuser(lab_labi_scratch_master) + &
threadno*iuser(i_liuser_callback_scratch)
iuser(lab_labi_scratch) = labi_scratch
lab_labr_scratch_master = iuser(i_label_ruser_scratch)
lab_labr_scratch = lab_labr_scratch_master + threadno
labr_scratch = iuser(lab_labr_scratch_master) + &
threadno*iuser(i_lruser_callback_scratch)
iuser(lab_labr_scratch) = labr_scratch
Case (1)
! OBJFUN is called on entry to a NAG local minimizer.
Case (0)
! This will be the normal value of NSTATE.
Case Default
! This should never happen, and indicates that an error has
! occurred on the system. It is used here to demonstrate how
! errors may be handled in a thread safe manner, in this case
! by performing indicative IO in a critical section.
mode = -1
!$Omp Critical (e05sbfe_objfun_error)
Write (nout,99999) threadno
!$Omp End Critical (e05sbfe_objfun_error)
Go To 100
End Select
! Test MODE to determine whether to calculate OBJF and/or OBJGRD.
evalobjf = .False.
evalobjg = .False.
Select Case (mode)
Case (0,5)
! Only the value of the objective function is needed.
evalobjf = .True.
Case (1,6)
! Only the values of the NDIM gradients are required.
evalobjg = .True.
Case (2,7)
! Both the objective function and the NDIM gradients are required.
evalobjf = .True.
evalobjg = .True.
End Select
! Set labels for shared and private workspaces.
labi_shared = iuser(iuser(i_label_iuser_callback_shared))
lab_labi_private = iuser(i_label_iuser_callback_private) + threadno
labi_private = iuser(lab_labi_private)
labr_shared = iuser(iuser(i_label_ruser_callback_shared))
lab_labr_private = iuser(i_label_ruser_callback_private) + threadno
labr_private = iuser(lab_labr_private)
! set labels for scratch workspaces.
lab_labi_scratch = iuser(i_label_iuser_scratch) + threadno
labi_scratch = iuser(lab_labi_scratch)
lab_labr_scratch = iuser(i_label_ruser_scratch) + threadno
labr_scratch = iuser(lab_labr_scratch)
If (evalobjf) Then
! Evaluate the objective function.
! objf = const + sum( x*sin(abs(sqrt(x))))
! Pass thread-safe workspace to objfun_schwefel_objective
Call objfun_schwefel_objective(ndim,x,objf,ruser(labr_scratch))
! Use value safely available to all threads.
objf = objf + ruser(labr_shared)
End If
If (evalobjg) Then
! Calculate the gradient of the objective function,
! and return the result in VECOUT.
! gradient = sin(sqrt(abs(x)))
! + sign(cos(sqrt(abs(x)))/2.0*sqrt(abs(x)),x)
Call objfun_schwefel_gradient(ndim,x,vecout,ruser(labr_scratch), &
ruser(labr_scratch+ndim))
End If
100 Continue
Return
99999 Format (1X,'** ERROR DETECTED : Thread number =',1X,I10)
End Subroutine objfun_schwefel
Subroutine confun_non_linear(mode,ncon,ndim,ldcj,needc,x,c,cjac,nstate, &
iuser,ruser)
! Subroutine used to supply constraints
! Unless NCON = 0, this will always be
! called before OBJFUN
! .. Use Statements ..
Use nag_library, Only: daxpy, dcopy, dscal, f06bmf, f06fjf, nag_wp, &
x06adf
! .. Implicit None Statement ..
Implicit None
! .. Scalar Arguments ..
Integer, Intent (In) :: ldcj, ncon, ndim, nstate
Integer, Intent (Inout) :: mode
! .. Array Arguments ..
Real (Kind=nag_wp), Intent (Out) :: c(ncon)
Real (Kind=nag_wp), Intent (Inout) :: cjac(ldcj,ndim), ruser(*)
Real (Kind=nag_wp), Intent (In) :: x(ndim)
Integer, Intent (Inout) :: iuser(*)
Integer, Intent (In) :: needc(ncon)
! .. Local Scalars ..
Real (Kind=nag_wp) :: l2xma, scal, sumsq
Integer :: conid, ifail, labi_private, &
labi_scratch, labi_shared, &
labr_private, labr_scratch, &
labr_shared, lab_labi_private, &
lab_labi_private_master, &
lab_labi_scratch, &
lab_labi_scratch_master, &
lab_labr_private, &
lab_labr_private_master, &
lab_labr_scratch, &
lab_labr_scratch_master, lruser_s, &
threadno
Logical :: evalc, evalcjac
! .. Executable Statements ..
threadno = x06adf()
! Test NSTATE to determine whether a new phase is beginning
Select Case (nstate)
Case (2)
! Very first call to confun.
! This occurs in serial on thread 0 only.
Do conid = 1, ncon
! Load details of required workspace for constraint k.
Call load_callback_data(conid,nin,iuser,ruser,ifail)
If (ifail/=0) Then
mode = ifail
Go To 100
End If
End Do
Case (3)
! First call to confun on a new thread.
! Set labels to alloted workspace in user arrays.
Do conid = 1, ncon
lab_labi_private_master = iuser(i_label_iuser_callback_private) + &
conid*iuser(i_max_threads)
lab_labi_private = lab_labi_private_master + threadno
labi_private = iuser(lab_labi_private_master) + &
threadno*iuser(iuser(i_liuser_callback_private)+conid)
iuser(lab_labi_private) = labi_private
lab_labr_private_master = iuser(i_label_ruser_callback_private) + &
conid*iuser(i_max_threads)
lab_labr_private = lab_labr_private_master + threadno
labr_private = iuser(lab_labr_private_master) + &
threadno*iuser(iuser(i_lruser_callback_private)+conid)
iuser(lab_labr_private) = labr_private
End Do
! Set private scratch labels.
lab_labi_scratch_master = iuser(i_label_iuser_scratch)
lab_labi_scratch = lab_labi_scratch_master + threadno
labi_scratch = iuser(lab_labi_scratch_master) + &
threadno*iuser(i_liuser_callback_scratch)
iuser(lab_labi_scratch) = labi_scratch
lab_labr_scratch_master = iuser(i_label_ruser_scratch)
lab_labr_scratch = lab_labr_scratch_master + threadno
labr_scratch = iuser(lab_labr_scratch_master) + &
threadno*iuser(i_lruser_callback_scratch)
iuser(lab_labr_scratch) = labr_scratch
Case (1)
! First call to confun at the start of a local minimzation.
! Set any constant elements of the Jacobian matrix, if any.
Case Default
! Standard call to confun.
End Select
! MODE: are constraints, derivatives, or both are required?
evalc = mode == 0 .Or. mode == 2
evalcjac = mode == 1 .Or. mode == 2
! Assign scratch workspace labels
lab_labi_scratch = iuser(i_label_iuser_scratch) + threadno
labi_scratch = iuser(lab_labi_scratch)
lab_labr_scratch = iuser(i_label_ruser_scratch) + threadno
labr_scratch = iuser(lab_labr_scratch)
loop_constraints: Do conid = 1, ncon
! Only those for which needc is non-zero need be set.
If (needc(conid)<=0) Then
Cycle loop_constraints
End If
! Assign workspace labels.
labi_shared = iuser(iuser(i_label_iuser_callback_shared)+conid)
labr_shared = iuser(iuser(i_label_ruser_callback_shared)+conid)
lruser_s = iuser(iuser(i_lruser_callback_shared)+conid)
lab_labi_private = iuser(i_label_iuser_callback_private) + &
iuser(i_max_threads)*conid + threadno
labi_private = iuser(lab_labi_private)
lab_labr_private = iuser(i_label_ruser_callback_private) + &
iuser(i_max_threads)*conid + threadno
labr_private = iuser(lab_labr_private)
Select Case (conid)
Case (1)
! ||X - SHARED_DATA ||_2
Call dcopy(ndim,x,1,ruser(labr_scratch),1)
Call daxpy(lruser_s,-1.0E0_nag_wp,ruser(labr_shared),1, &
ruser(labr_scratch),1)
scal = 0.0E0_nag_wp
sumsq = 1.0E0_nag_wp
Call f06fjf(ndim,ruser(labr_scratch),1,scal,sumsq)
l2xma = f06bmf(scal,sumsq)
If (evalc) Then
c(conid) = l2xma
! update private constraint counter
iuser(labi_private) = iuser(labi_private) + 1
! safely update shared constraint counter
!$Omp Critical (e05sbfe_smp_update_c_count)
iuser(labi_shared) = iuser(labi_shared) + 1
!$Omp End Critical (e05sbfe_smp_update_c_count)
End If
If (evalcjac) Then
If (l2xma>0.0E0_nag_wp) Then
l2xma = 1.0E0_nag_wp/l2xma
End If
Call dscal(ndim,l2xma,ruser(labr_scratch),1)
Call dcopy(ndim,ruser(labr_scratch),1,cjac(conid,1),ldcj)
! update private constraint jacobian counter
iuser(labi_private+1) = iuser(labi_private+1) + 1
! Safely update shared constraint jacobian counter
!$Omp Critical (e05sbfe_smp_update_cjac_count)
iuser(labi_shared+1) = iuser(labi_shared+1) + 1
!$Omp End Critical (e05sbfe_smp_update_cjac_count)
End If
End Select
End Do loop_constraints
! If an immediate exit is required, return MODE<0
100 Continue
Return
End Subroutine confun_non_linear
Subroutine monmod(ndim,ncon,npar,x,xb,fb,cb,xbest,fbest,cbest,itt,iuser, &
ruser,inform)
! .. Implicit None Statement ..
Implicit None
! .. Scalar Arguments ..
Real (Kind=nag_wp), Intent (In) :: fb
Integer, Intent (Inout) :: inform
Integer, Intent (In) :: ncon, ndim, npar
! .. Array Arguments ..
Real (Kind=nag_wp), Intent (In) :: cb(ncon), cbest(ncon,npar), &
fbest(npar), xb(ndim), &
xbest(ndim,npar)
Real (Kind=nag_wp), Intent (Inout) :: ruser(*), x(ndim,npar)
Integer, Intent (In) :: itt(7)
Integer, Intent (Inout) :: iuser(*)
! .. Local Scalars ..
Integer :: fid, k, labi_con1_private, &
labi_con1_shared, &
lab_labi_con1_private, &
lab_labi_con1_shared
Character (18) :: fname
Character (3) :: leftno
! .. Intrinsic Procedures ..
Intrinsic :: adjustl, modulo, trim
! .. Executable Statements ..
If (detail_level>=2) Then
! Report on the first iteration, and every report_freq iterations.
If (itt(1)==1 .Or. modulo(itt(1),report_freq)==0 .And. inform<1000) &
Then
! To avoid using critical sections about IO,
! allow each thread to open an output file.
! Note, on entry INFORM = THREADNO.
fid = 10 + inform
lab_labi_con1_shared = iuser(i_label_iuser_callback_shared) + 1
labi_con1_shared = iuser(lab_labi_con1_shared)
lab_labi_con1_private = iuser(i_label_iuser_callback_private) + &
iuser(i_max_threads)*1 + inform
labi_con1_private = iuser(lab_labi_con1_private)
Write (leftno,*) inform
Write (fname,99995) trim(adjustl(leftno))
Open (fid,File=fname,Position='APPEND')
Write (fid,*) '* Iteration *', itt(1)
Write (fid,*) '* Total function evaluations *', itt(5)
Write (fid,*) '* Constraint 1 and Jacobian Evaluations *'
Write (fid,*) '* Total :', iuser(labi_con1_shared), &
iuser(labi_con1_private)
Write (fid,*) '* Current global optimum candidate:'
Do k = 1, ndim
Write (fid,99999) k, xb(k)
End Do
Write (fid,*) '* Current global optimum value:'
Write (fid,99998) fb
Write (fid,99997)
Do k = 1, ncon
Write (fid,99996) k, cb(k)
End Do
Close (fid)
End If
End If
! If required set INFORM<0 to force exit
Return
99999 Format (1X,'* xb(',I3,') = ',F9.2)
99998 Format (1X,'* fb = ',F9.5)
99997 Format ('** Current global optimum constraint violations **')
99996 Format (1X,'* cb(',I3,') = ',F9.2)
99995 Format ('e05sbfe_smp_',A,'.r')
End Subroutine monmod
Subroutine display_result(ndim,ncon,xb,fb,cb,itt,inform)
! Display final results in comparison to known global optimum.
! .. Implicit None Statement ..
Implicit None
! .. Scalar Arguments ..
Real (Kind=nag_wp), Intent (In) :: fb
Integer, Intent (In) :: inform, ncon, ndim
! .. Array Arguments ..
Real (Kind=nag_wp), Intent (In) :: cb(ncon), xb(ndim)
Integer, Intent (In) :: itt(7)
! .. Local Scalars ..
Integer :: k
! .. Executable Statements ..
! Display final counters.
Write (nout,*) ' Algorithm Statistics'
Write (nout,*) ' --------------------'
Write (nout,99992) 'Total complete iterations ', itt(1)
Write (nout,99992) 'Complete iterations since improvement ', itt(2)
Write (nout,99992) 'Total particles converged to xb ', itt(3)
Write (nout,99992) 'Total improvements to global optimum ', itt(4)
Write (nout,99992) 'Total function evaluations ', itt(5)
Write (nout,99992) 'Total particles re-initialized ', itt(6)
Write (nout,99992) 'Total constraints violated ', itt(7)
! Display why finalization occurred.
Write (nout,*)
Select Case (inform)
Case (1)
Write (nout,99999) 'Target value achieved'
Case (2)
Write (nout,99999) 'Minimum swarm standard deviation obtained'
Case (3)
Write (nout,99999) 'Sufficient particles converged'
Case (4)
Write (nout,99999) 'No improvement in preset iteration limit'
Case (5)
Write (nout,99999) 'Maximum complete iterations attained'
Case (6)
Write (nout,99999) 'Maximum function evaluations exceeded'
Case (7)
Write (nout,99999) 'Constrained point located'
Case (:-1)
Write (nout,99998) inform
Case Default
End Select
! Display final objective value and location.
Write (nout,*)
Write (nout,99997) fb
Write (nout,99996)
Do k = 1, ndim
Write (nout,99995) k, xb(k)
End Do
Write (nout,99994)
Do k = 1, ncon
Write (nout,99993) k, cb(k)
End Do
Write (nout,*)
Return
99999 Format (2X,'Solution Status : ',A38)
99998 Format (' User termination case : ',I13)
99997 Format (' Achieved objective value : ',F13.3)
99996 Format (' Final position achieved : ')
99995 Format (' xb(',I4,') = ',F12.2)
99994 Format (' Final constraint violations achieved : ')
99993 Format (' cb(',I4,') = ',F12.2)
99992 Format (2X,A40,' :',I13)
End Subroutine display_result
End Module e05sbfe_smp_mod
Program e05sbfe_smp
! E05SBF SMP Example Main Program
! .. Use Statements ..
Use e05sbfe_smp_mod, Only: calculate_workspace, confun_non_linear, &
display_result, initialize_callback_arrays, &
i_liuser_req, i_lruser_req, liopts, lopts, &
min_liuser, min_lruser, monmod, nag_wp, nin, &
nout, objfun_schwefel, zero
Use nag_library, Only: e05sbf, e05zkf, e05zlf, x06acf
! .. Implicit None Statement ..
Implicit None
! .. Local Scalars ..
Real (Kind=nag_wp) :: fb, rvalue
Integer :: astat, ifail, inform, ivalue, k, &
liuser, lruser, max_threads, ncon, &
ndim, npar, optype
Character (16) :: cvalue
Character (80) :: optstr
! .. Local Arrays ..
Real (Kind=nag_wp), Allocatable :: bl(:), bu(:), cb(:), cbest(:,:), &
fbest(:), ruser(:), xb(:), &
xbest(:,:)
Real (Kind=nag_wp) :: opts(lopts)
Integer :: iopts(liopts), itt(7)
Integer, Allocatable :: iuser(:)
! .. Intrinsic Procedures ..
Intrinsic :: max
! .. Executable Statements ..
! Print advisory information.
Write (nout,*) 'E05SBF Example Program Results'
Write (nout,*)
Write (nout,*) 'Minimization of the Schwefel function.'
Write (nout,*) 'Subject to one non-linear constraint.'
Write (nout,*)
! Skip heading in data file
Read (nin,*)
! Read parameters
Read (nin,*) ndim
Read (nin,*) ncon
max_threads = x06acf()
npar = 10*max(max_threads,ndim)
Allocate (xb(ndim),cb(ncon),xbest(ndim,npar),cbest(ncon,npar), &
fbest(npar),bl(ndim+ncon),bu(ndim+ncon),Stat=astat)
xbest = zero
fbest = zero
cbest = zero
! Set problem specific values.
! Set box bounds
bl(1:ndim) = -500.0E0_nag_wp
bu(1:ndim) = 500.0E0_nag_wp
! Read general constraint bounds.
Do k = 1, ncon
Read (nin,*) bl(k+ndim), bu(k+ndim)
End Do
! Calculate size of LIUSER and LRUSER required.
Call calculate_workspace(ndim,ncon,max_threads,min_liuser,min_lruser, &
nin,liuser,lruser,ifail)
! Initialize the option arrays for E05SBF.
ifail = 0
Call e05zkf('Initialize = E05SBF',iopts,liopts,opts,lopts,ifail)
! ------------------------------------------------------------------
Write (nout,*)
Write (nout,*) '1. Solution without using coupled local minimizer'
Write (nout,*)
! Set various options to non-default values if required.
ifail = 0
Write (optstr,99999) 'Constraint Tolerance', 1.0E-5_nag_wp
Call e05zkf(optstr,iopts,liopts,opts,lopts,ifail)
ifail = 0
Write (optstr,99999) 'Constraint Superiority', 1.0E-6_nag_wp
Call e05zkf(optstr,iopts,liopts,opts,lopts,ifail)
ifail = 0
Call e05zkf('Constraint Norm = L2',iopts,liopts,opts,lopts,ifail)
ifail = 0
Call e05zkf('Repeatability = On',iopts,liopts,opts,lopts,ifail)
ifail = 0
Write (optstr,99999) 'Distance Tolerance', 1.0E-3_nag_wp
Call e05zkf(optstr,iopts,liopts,opts,lopts,ifail)
ifail = 0
Call e05zkf('Maximum Iterations Static = 1000',iopts,liopts,opts,lopts, &
ifail)
ifail = 0
Call e05zkf('Maximum Iterations Static Particles = 0',iopts,liopts,opts, &
lopts,ifail)
ifail = 0
Call e05zkf('Constraint Warning = off',iopts,liopts,opts,lopts,ifail)
ifail = 0
Call e05zkf('Repulsion Initialize = 50',iopts,liopts,opts,lopts,ifail)
ifail = 0
Call e05zkf('Repulsion Finalize = 20',iopts,liopts,opts,lopts,ifail)
ifail = 0
Call e05zkf('Repulsion Particles = 1',iopts,liopts,opts,lopts,ifail)
ifail = 0
Call e05zkf('Boundary = Hyperspherical',iopts,liopts,opts,lopts,ifail)
ifail = 0
! It is essential to indicate callback thread safety.
Call e05zkf('SMP Callback Thread Safe = Yes',iopts,liopts,opts,lopts, &
ifail)
ifail = 0
Call e05zkf('SMP Thread Overrun = 10',iopts,liopts,opts,lopts,ifail)
ifail = 0
! Call E05SBF to search for the global optimum.
! This is incorporated into a loop so that, in the unlikely
! event of insufficient space in iuser and ruser, the code
! will automatically try to allocate sufficient space.
loop_no_min: Do
loop_initialize_user_arrays: Do
Allocate (iuser(liuser),ruser(lruser),Stat=astat)
Call initialize_callback_arrays(ndim,ncon,iuser,liuser,ruser,lruser, &
ifail)
If (ifail==0) Then
Exit loop_initialize_user_arrays
Else
liuser = max(iuser(i_liuser_req),liuser)
lruser = max(iuser(i_lruser_req),lruser)
End If
Deallocate (iuser,ruser)
End Do loop_initialize_user_arrays
! Non-zero IFAIL expected on exit here, so use IFAIL=1 (quiet) on entry.
ifail = 1
Call e05sbf(ndim,ncon,npar,xb,fb,cb,bl,bu,xbest,fbest,cbest, &
objfun_schwefel,confun_non_linear,monmod,iopts,opts,iuser,ruser,itt, &
inform,ifail)
! It is essential to test IFAIL on exit.
Select Case (ifail)
Case (0,1)
! No errors, best found optimum at xb returned in fb.
Call display_result(ndim,ncon,xb,fb,cb,itt,inform)
Case (3)
! Exit flag set in OBJFUN, CONFUN or MONMOD and returned in INFORM.
Select Case (inform)
Case (-2)
! Insufficient memory ineither iuser or ruser
liuser = max(iuser(i_liuser_req),liuser)
lruser = max(iuser(i_lruser_req),lruser)
Deallocate (iuser,ruser)
Cycle loop_no_min
Case Default
Call display_result(ndim,ncon,xb,fb,cb,itt,inform)
End Select
Case Default
! An error was detected. Print message since IFAIL=1 on entry.
Write (nout,99998) '** E05SBF returned with an error, IFAIL = ', &
ifail
Continue
End Select
Exit loop_no_min
End Do loop_no_min
! Deallocate user arrays. We will re-initialise them for the second
! solution below.
Deallocate (iuser,ruser)
! ------------------------------------------------------------------
Write (nout,*) '2. Solution using coupled local minimizer E04UCF'
Write (nout,*)
! Set the local minimizer to be E04UCF and set corresponding options.
ifail = 0
Call e05zkf('Local Minimizer = E04UCF',iopts,liopts,opts,lopts,ifail)
ifail = 0
Call e05zkf('Local Interior Major Iterations = 15',iopts,liopts,opts, &
lopts,ifail)
ifail = 0
Call e05zkf('Local Interior Minor Iterations = 5',iopts,liopts,opts, &
lopts,ifail)
ifail = 0
Call e05zkf('Local Exterior Major Iterations = 50',iopts,liopts,opts, &
lopts,ifail)
ifail = 0
Call e05zkf('Local Exterior Minor Iterations = 15',iopts,liopts,opts, &
lopts,ifail)
! Query the option Distance Tolerance
ifail = 0
optstr = 'Distance Tolerance'
Call e05zlf(optstr,ivalue,rvalue,cvalue,optype,iopts,opts,ifail)
! Adjust Distance Tolerance dependent upon its current value
Write (optstr,99999) 'Distance Tolerance', rvalue*10.0_nag_wp
ifail = 0
Call e05zkf(optstr,iopts,liopts,opts,lopts,ifail)
ifail = 0
Write (optstr,99999) 'Local Interior Tolerance', rvalue
Call e05zkf(optstr,iopts,liopts,opts,lopts,ifail)
ifail = 0
Write (optstr,99999) 'Local Exterior Tolerance', rvalue*1.0E-4_nag_wp
Call e05zkf(optstr,iopts,liopts,opts,lopts,ifail)
! Call E05SBF to search for the global optimum using coupled
! local minimizer E04UCF.
! This is incorporated into a loop so that, in the unlikely
! event of insufficient space in iuser and ruser, the code
! will automatically try to allocate sufficient space.
loop_min_e04ucf: Do
loop_initialize_user_arrays_2: Do
Allocate (iuser(liuser),ruser(lruser),Stat=astat)
Call initialize_callback_arrays(ndim,ncon,iuser,liuser,ruser,lruser, &
ifail)
If (ifail==0) Then
Exit loop_initialize_user_arrays_2
Else
liuser = max(iuser(i_liuser_req),liuser)
lruser = max(iuser(i_lruser_req),lruser)
End If
Deallocate (iuser,ruser)
End Do loop_initialize_user_arrays_2
! Non-zero IFAIL expected on exit here, so use IFAIL=1 (quiet) on entry.
ifail = 1
Call e05sbf(ndim,ncon,npar,xb,fb,cb,bl,bu,xbest,fbest,cbest, &
objfun_schwefel,confun_non_linear,monmod,iopts,opts,iuser,ruser,itt, &
inform,ifail)
! It is essential to test IFAIL on exit.
Select Case (ifail)
Case (0,1)
! No errors, best found optimum at xb returned in fb.
Call display_result(ndim,ncon,xb,fb,cb,itt,inform)
Case (3)
! Exit flag set in OBJFUN, CONFUN or MONMOD and returned in INFORM.
Select Case (inform)
Case (-2)
! Insufficient memory ineither iuser or ruser
liuser = max(iuser(i_liuser_req),liuser)
lruser = max(iuser(i_lruser_req),lruser)
Deallocate (iuser,ruser)
Cycle loop_min_e04ucf
Case Default
Call display_result(ndim,ncon,xb,fb,cb,itt,inform)
End Select
Case Default
! An error was detected. Print message since IFAIL=1 on entry.
Write (nout,99998) '** E05SBF returned with an error, IFAIL = ', &
ifail
Continue
End Select
Exit loop_min_e04ucf
End Do loop_min_e04ucf
99999 Format (A,' = ',E32.16)
99998 Format (1X,A,I6)
End Program e05sbfe_smp