!   E05SAF Example Program Text
!   Mark 24 Release. NAG Copyright 2012.

    Module e05safe_mod

!     E05SAF 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        ::                                  &
                                         f_target = -4.189828872724337E2_nag_wp
      Real (Kind=nag_wp), Parameter        :: one = 1.0E0_nag_wp
      Real (Kind=nag_wp), Parameter        :: zero = 0.0_nag_wp
      Integer, Parameter                   :: detail_level = 0,                &
                                              iuser_reference = 1,             &
                                              liopts = 100, lopts = 100,       &
                                              ndim = 20, nout = 6,             &
                                              report_freq = 100,               &
                                              ruser_global = 1
      Real (Kind=nag_wp), Parameter        ::                                  &
                                   x_target(1:ndim) = -420.9687463599820_nag_wp
    Contains
      Subroutine objfun_schwefel_objective(ndim,x,f,w)
!       subroutine to calculate the objective using provided workspace

!       .. Use Statements ..
        Use nag_library, Only: ddot
!       .. 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
!       .. 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 E05SAF.

!       .. 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                              :: threadno, work_1, work_2
        Logical                              :: evalobjf, evalobjg
!       .. External Procedures ..
        Integer, External                    :: omp_get_thread_num
!       .. Intrinsic Procedures ..
        Intrinsic                            :: real
!       .. Executable Statements ..
!       Call OMP_GET_NUM_THREAD to get the current thread number.
!       The use of '!$' also allows this example to be compiled in serial
!       without openMP libraries.
        threadno = 0
!$      threadno = omp_get_thread_num()

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

!         Set up any constant global variables required in RUSER and IUSER
!         This may also be done prior to calling E05SAF
          ruser(ruser_global) = -f_target*real(ndim,nag_wp)

!         Set reference index for thread 0 workspace in IUSER
          iuser(iuser_reference+threadno) = ruser_global + 1

        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 reference index for this
!         threads unique section of RUSER workspace in IUSER.

!         Each thead requires 2*ndim real numbers.
          iuser(iuser_reference+threadno) = iuser(iuser_reference) + &
            2*ndim*threadno

        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 is extremely unlikely, and indicates that an error has
!         occurred on the system
          mode = -1
          If (detail_level>=2) Then
            Write (nout,99999)
          End If
          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.
        Case Default
          mode = -1
          Write (nout,99999)
          Go To 100
        End Select

!       Set work_1 and work_2 to point to the section of RUSER reserved
!       for this thread
        work_1 = iuser(iuser_reference+threadno)
        work_2 = work_1 + ndim

        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(work_1))
!         Use constant safely available to all threads.
          objf = objf + ruser(ruser_global)

        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(work_1), &
            ruser(work_2))
        End If

100     Continue
        Return

99999   Format (1X,'** ERROR DETECTED ')

      End Subroutine objfun_schwefel
      Subroutine monmod(ndim,npar,x,xb,fb,xbest,fbest,itt,iuser,ruser,inform)
!       Example monitor function.

!       .. Scalar Arguments ..
        Real (Kind=nag_wp), Intent (In)      :: fb
        Integer, Intent (Inout)              :: inform
        Integer, Intent (In)                 :: ndim, npar
!       .. Array Arguments ..
        Real (Kind=nag_wp), Intent (In)      :: fbest(npar), xb(ndim),         &
                                                xbest(ndim,npar)
        Real (Kind=nag_wp), Intent (Inout)   :: ruser(*), x(ndim,npar)
        Integer, Intent (In)                 :: itt(6)
        Integer, Intent (Inout)              :: iuser(*)
!       .. Local Scalars ..
        Integer                              :: k
!       .. Intrinsic Procedures ..
        Intrinsic                            :: modulo
!       .. 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) Then
            Write (nout,*)
            Write (nout,*) '* Current global optimum candidate:'
            Do k = 1, ndim
              Write (nout,99999) k, xb(k)
            End Do
            Write (nout,*) '* Current global optimum value:'
            Write (nout,99998) fb
          End If
        End If

!       If required set INFORM<0 to force exit
        Flush (nout)

        Return
99999   Format (1X,'* xb(',I3,') = ',F9.2)
99998   Format (1X,'* fb = ',F9.5)

      End Subroutine monmod
      Subroutine display_result(ndim,xb,fb,itt,inform,ruser)
!       Display final results in comparison to known global optimum.

!       .. Use Statements ..
        Use nag_library, Only: x04cbf
!       .. Parameters ..
        Integer, Parameter                   :: indent = 0, ncols = 80
        Character (10), Parameter            :: clabs(1:2) = (/                &
                                                'x_target  ','xb        '  /)
        Character (1), Parameter             :: diag =   'N',  labcol =        &
                                                'C',  labrow =   'I',          &
                                                matrix =   'G'
        Character (4), Parameter             :: form = 'f9.2'
!       .. Scalar Arguments ..
        Real (Kind=nag_wp), Intent (In)      :: fb
        Integer, Intent (In)                 :: inform, ndim
!       .. Array Arguments ..
        Real (Kind=nag_wp), Intent (In)      :: ruser(*), xb(ndim)
        Integer, Intent (In)                 :: itt(6)
!       .. Local Scalars ..
        Integer                              :: ifail, ldxcom
        Character (ncols)                    :: title
!       .. Local Arrays ..
        Real (Kind=nag_wp), Allocatable      :: xcom(:,:)
!       .. Executable Statements ..
!       Display final counters.
        Write (nout,99996) 'Total complete iterations             : ', itt(1)
        Write (nout,99996) 'Complete iterations since improvement : ', itt(2)
        Write (nout,99996) 'Total particles converged to xb       : ', itt(3)
        Write (nout,99996) 'Total improvements to global optimum  : ', itt(4)
        Write (nout,99996) 'Total function evaluations            : ', itt(5)
        Write (nout,99996) 'Total particles re-initialized        : ', itt(6)

!       Display why finalization occurred.
        Write (nout,*)
        Select Case (inform)
        Case (1)
          Write (nout,99995) 'Target value achieved'
        Case (2)
          Write (nout,99995) 'Minimum swarm standard deviation obtained'
        Case (3)
          Write (nout,99995) 'Sufficient particles converged'
        Case (4)
          Write (nout,99995) 'No improvement in preset iteration limit'
        Case (5)
          Write (nout,99995) 'Maximum complete iterations attained'
        Case (6)
          Write (nout,99995) 'Maximum function evaluations exceeded'
        Case (:-1)
          Write (nout,99994) inform
        Case Default
        End Select

!       Display final objective value and location.
        Write (nout,*)
        Write (nout,99999) 0.0E0_nag_wp
        Write (nout,99998) fb
        Flush (nout)

        ldxcom = ndim
        Allocate (xcom(ldxcom,2))
        xcom(1:ndim,1) = x_target(1:ndim)
        xcom(1:ndim,2) = xb(1:ndim)

        title = ' Comparison between known and achieved optima.'
        ifail = 0
        Call x04cbf(matrix,diag,ndim,2,xcom,ldxcom,form,title,labrow,clabs, &
          labcol,clabs,ncols,indent,ifail)

        Deallocate (xcom)

        Write (nout,99997)

        Return
99999   Format ('  Known objective optimum  : ',F13.5)
99998   Format ('  Achieved objective value : ',F13.5)
99997   Format ('**--------------------**--------------------**'/)
99996   Format (2X,A40,I16)
99995   Format (2X,A43)
99994   Format ('  User termination case                    :',I16)
      End Subroutine display_result
    End Module e05safe_mod

    Program e05safe
!     E05SAF Example Main Program

!     This example program demonstrates how to use E05SAF in standard
!       execution, and with a selection of coupled local minimizers.
!     The non-default option 'REPEATABILITY ON' is used here, giving
!       repeatable results.

!     .. Use Statements ..
      Use nag_library, Only: e05saf, e05zkf
      Use e05safe_mod, Only: display_result, iuser_reference, liopts, lopts,   &
                             monmod, nag_wp, ndim, nout, objfun_schwefel,      &
                             ruser_global
!     .. Implicit None Statement ..
      Implicit None
!     .. Local Scalars ..
      Real (Kind=nag_wp)                   :: fb
      Integer                              :: ifail, inform, liuser, lruser,   &
                                              max_threads, npar
      Character (80)                       :: optstr
!     .. Local Arrays ..
      Real (Kind=nag_wp)                   :: bl(ndim), bu(ndim), opts(lopts), &
                                              xb(ndim)
      Real (Kind=nag_wp), Allocatable      :: ruser(:)
      Integer                              :: iopts(liopts), itt(6)
      Integer, Allocatable                 :: iuser(:)
!     .. External Procedures ..
      Integer, External                    :: omp_get_max_threads
!     .. Intrinsic Procedures ..
      Intrinsic                            :: max
!     .. Executable Statements ..
!     Print advisory information.
      Write (nout,*) 'E05SAF Example Program Results'
      Write (nout,*)
      Write (nout,*) 'Minimization of the Schwefel function.'
      Write (nout,*)

!     Set problem specific values.
!     Set box bounds.       
      bl(1:ndim) = -500.0_nag_wp
      bu(1:ndim) = 500.0_nag_wp

!     Determine the number of threads to be used in the simulation
      max_threads = 1
!$    max_threads = omp_get_max_threads()

!     Ensure NPAR is sufficient for the number of threads
      npar = max(4000,5*max_threads)

!     Assign sufficient sizes for iuser and ruser.
!     In this example, RUSER will be split into workspace unique to each thread,
!     and IUSER will be used to store the indexes to the first element of 
!     workspace for each thread.

      liuser = max(30,iuser_reference+max_threads)
      lruser = ruser_global + 2*max_threads*ndim
      Allocate (iuser(liuser),ruser(lruser))

      Write (nout,*) '* Solution without using coupled local minimizer *'

!     Call E05ZKF to initialize the option arrays
      ifail = 0
      Call e05zkf('Initialize = E05SAF',iopts,liopts,opts,lopts,ifail)

!     Set the option SMP Callback Thread Safe to indicate the callback functions
!     are indeed threadsafe.  This must be done if using multiple threads.

      ifail = 0
      Call e05zkf('SMP Callback Thread Safe = Yes',iopts,liopts,opts,lopts, &
        ifail)

!     Set various options to non-default values if required.
      ifail = 0
      Call e05zkf('Repeatability = On',iopts,liopts,opts,lopts,ifail)
      ifail = 0
      Call e05zkf('Verify Gradients = Off',iopts,liopts,opts,lopts,ifail)
      ifail = 0
      Call e05zkf('Boundary = Hyperspherical',iopts,liopts,opts,lopts,ifail)
      ifail = 0
      Call e05zkf('Maximum iterations static = 150',iopts,liopts,opts,lopts, &
        ifail)
      ifail = 0
      Call e05zkf('Maximum function evaluations = 50000000',iopts,liopts,opts, &
        lopts,ifail)
      ifail = 0
      Call e05zkf('Repulsion Initialize = 50',iopts,liopts,opts,lopts,ifail)
      ifail = 0
      Call e05zkf('Repulsion Finalize = 30',iopts,liopts,opts,lopts,ifail)

!     Set an objective target.
      ifail = 0
      Write (optstr,99999) 'Target Objective Value', 1.0E0_nag_wp
      Call e05zkf(optstr,iopts,liopts,opts,lopts,ifail)
      ifail = 0
      Write (optstr,99999) 'Target Objective Tolerance', 1.0E-5_nag_wp
      Call e05zkf(optstr,iopts,liopts,opts,lopts,ifail)
      ifail = 0
      Write (optstr,99999) 'Target Objective Safeguard', 1.0E-2_nag_wp
      Call e05zkf(optstr,iopts,liopts,opts,lopts,ifail)

!     Set some additional SMP library options if required.
      ifail = 0
      Call e05zkf('SMP thread overrun = 100',iopts,liopts,opts,lopts,ifail)

      ifail = 0
      Write (optstr,99998) 'SMP Subswarm', max_threads
      Call e05zkf(optstr,iopts,liopts,opts,lopts,ifail)

!     Call E05SAF to search for the global optimum.
!     As guaranteed success cannot by default be achieved, IFAIL=1,-1 is
!     recommended on entry.
      ifail = -1
      Call e05saf(ndim,npar,xb,fb,bl,bu,objfun_schwefel,monmod,iopts,opts, &
        iuser,ruser,itt,inform,ifail)
!     It is essential to test IFAIL on exit.
      Select Case (ifail)
      Case (0,1)
!       E05SAF encountered no errors during operation,
!       and will have returned the best found optimum.
        Call display_result(ndim,xb,fb,itt,inform,ruser)
      Case (-999,-399,2:)
!       An error was detected.
        Continue
      Case Default
!       An instruction to exit was received by E05SAF from OBJFUN or MONMOD.
!       The exit flag will have been returned in INFORM.
        Call display_result(ndim,xb,fb,itt,inform,ruser)
      End Select

!     ------------------------------------------------------------------
      Write (nout,*) '* Solution using coupled local minimizer E04JYF *'



!     Set the local minimizer to be E04JYF and set corresponding options.
      ifail = 0
      Call e05zkf('Local Minimizer = E04JYF',iopts,liopts,opts,lopts,ifail)
      ifail = 0
      Call e05zkf('Local Interior Iterations = 30',iopts,liopts,opts,lopts, &
        ifail)
      ifail = 0
      Call e05zkf('Local Exterior Iterations = 200',iopts,liopts,opts,lopts, &
        ifail)
      ifail = 0
      Write (optstr,99999) 'Local Interior Tolerance', 1.0E-4_nag_wp
      Call e05zkf(optstr,iopts,liopts,opts,lopts,ifail)
      ifail = 0
      Write (optstr,99999) 'Local Exterior Tolerance', 1.0E-8_nag_wp
      Call e05zkf(optstr,iopts,liopts,opts,lopts,ifail)

!     Call E05SAF to search for the global optimum.
      ifail = -1
      Call e05saf(ndim,npar,xb,fb,bl,bu,objfun_schwefel,monmod,iopts,opts, &
        iuser,ruser,itt,inform,ifail)

!     It is essential to test IFAIL on exit.
      Select Case (ifail)
      Case (0,1)
!       E05SAF encountered no errors during operation,
!       and will have returned the best found optimum.
        Call display_result(ndim,xb,fb,itt,inform,ruser)
      Case (-999,-399,2:)
!       An error was detected.
        Continue
      Case Default
!       An instruction to exit was received by E05SAF from OBJFUN or MONMOD.
!       The exit flag will have been returned in INFORM.
        Call display_result(ndim,xb,fb,itt,inform,ruser)
      End Select

!     -----------------------------------------------------------------
      Write (nout,*) '* Solution using coupled local minimizer E04UCF *'

!     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 Iterations = 10',iopts,liopts,opts,lopts, &
        ifail)
      ifail = 0
      Call e05zkf('Local Exterior Iterations = 100',iopts,liopts,opts,lopts, &
        ifail)

!     Call E05SAF to search for the global optimum.
      ifail = -1
      Call e05saf(ndim,npar,xb,fb,bl,bu,objfun_schwefel,monmod,iopts,opts, &
        iuser,ruser,itt,inform,ifail)

!     It is essential to test IFAIL on exit.
      Select Case (ifail)
      Case (0,1)
!       E05SAF encountered no errors during operation,
!       and will have returned the best found optimum.
        Call display_result(ndim,xb,fb,itt,inform,ruser)
      Case (-999,-399,2:)
!       An error was detected.
        Continue
      Case Default
!       An instruction to exit was received by E05SAF from OBJFUN or MONMOD.
!       The exit flag will have been returned in INFORM.
        Call display_result(ndim,xb,fb,itt,inform,ruser)
      End Select

99999 Format (A,' = ',E32.16)
99998 Format (A,' = ',I20)
    End Program e05safe
