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

NAG FL Interface Introduction
Example description
!   E05SBF Example Program Text
!   Mark 30.1 Release. NAG Copyright 2024.

    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 ..

        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
        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
        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