The ParaMonte Documentation Website
Current view: top level - kernel - Optimization_mod.f90 (source / functions) Hit Total Coverage
Test: ParaMonte 1.5.1 :: Serial Kernel - Code Coverage Report Lines: 174 185 94.1 %
Date: 2021-01-08 13:03:42 Functions: 7 7 100.0 %
Legend: Lines: hit not hit

          Line data    Source code
       1             : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
       2             : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
       3             : !!!!
       4             : !!!!   MIT License
       5             : !!!!
       6             : !!!!   ParaMonte: plain powerful parallel Monte Carlo library.
       7             : !!!!
       8             : !!!!   Copyright (C) 2012-present, The Computational Data Science Lab
       9             : !!!!
      10             : !!!!   This file is part of the ParaMonte library.
      11             : !!!!
      12             : !!!!   Permission is hereby granted, free of charge, to any person obtaining a
      13             : !!!!   copy of this software and associated documentation files (the "Software"),
      14             : !!!!   to deal in the Software without restriction, including without limitation
      15             : !!!!   the rights to use, copy, modify, merge, publish, distribute, sublicense,
      16             : !!!!   and/or sell copies of the Software, and to permit persons to whom the
      17             : !!!!   Software is furnished to do so, subject to the following conditions:
      18             : !!!!
      19             : !!!!   The above copyright notice and this permission notice shall be
      20             : !!!!   included in all copies or substantial portions of the Software.
      21             : !!!!
      22             : !!!!   THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
      23             : !!!!   EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
      24             : !!!!   MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.
      25             : !!!!   IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM,
      26             : !!!!   DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR
      27             : !!!!   OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE
      28             : !!!!   OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
      29             : !!!!
      30             : !!!!   ACKNOWLEDGMENT
      31             : !!!!
      32             : !!!!   ParaMonte is an honor-ware and its currency is acknowledgment and citations.
      33             : !!!!   As per the ParaMonte library license agreement terms, if you use any parts of
      34             : !!!!   this library for any purposes, kindly acknowledge the use of ParaMonte in your
      35             : !!!!   work (education/research/industry/development/...) by citing the ParaMonte
      36             : !!!!   library as described on this page:
      37             : !!!!
      38             : !!!!       https://github.com/cdslaborg/paramonte/blob/main/ACKNOWLEDGMENT.md
      39             : !!!!
      40             : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
      41             : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
      42             : 
      43             : !>  \brief This module contains numerical optimization procedures.
      44             : !>  \author Amir Shahmoradi
      45             : 
      46             : module Optimization_mod
      47             : 
      48             :     use Constants_mod, only: IK, RK
      49             :     use Err_mod, only: Err_type
      50             : 
      51             :     implicit none
      52             : 
      53             :     character(*), parameter :: MODULE_NAME = "@Optimization_mod"
      54             : 
      55             : #if defined OS_IS_WSL
      56             :     procedure(getFuncMD_proc), pointer  :: getFuncMD_WSL                    !< This madness bypasses the Microsoft Subsystem for Linux Internal Function call GFortran Segmentation Fault error.
      57             :     real(RK), allocatable               :: StartVec_WSL(:), DirVec_WSL(:)   !< This madness bypasses the Microsoft Subsystem for Linux Internal Function call GFortran Segmentation Fault error.
      58             :     integer(IK)                         :: ndim_WSL                         !< This madness bypasses the Microsoft Subsystem for Linux Internal Function call GFortran Segmentation Fault error.
      59             : #endif
      60             : 
      61             : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
      62             : 
      63             :     !> The Brent minimizer class.
      64             :     type :: BrentMinimum_type
      65             :         integer(IK)     :: niter                        !< the number of iterations to reach the minimum of the function
      66             :         real(RK)        :: Bracket(3)                   !< the initial 3 Bracketing points that envelop the minimum
      67             :         real(RK)        :: xtol = sqrt(epsilon(1._RK))  !< the stopping rule tolerance
      68             :         real(RK)        :: xmin                         !< the x-value at the minimum of the function
      69             :         real(RK)        :: fmin                         !< the minimum of the function
      70             :         type(Err_type)  :: Err
      71             :     end type BrentMinimum_type
      72             : 
      73             :     ! overload Brent
      74             : 
      75             :     interface BrentMinimum_type
      76             :         module procedure :: minimizeBrent
      77             :     end interface BrentMinimum_type
      78             : 
      79             :     ! interfaces of the objective functions
      80             : 
      81             :     abstract interface
      82             :         function getFunc1D_proc(x) result(funcVal)
      83             :             import :: RK
      84             :             real(RK)    , intent(in)    :: x
      85             :             real(RK)                    :: funcVal
      86             :         end function getFunc1D_proc
      87             :     end interface
      88             : 
      89             : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
      90             : 
      91             :     !> The Powell minimizer class.
      92             :     type :: PowellMinimum_type
      93             :         integer(IK)             :: niter                        !< The number of iterations to reach the minimum of the function.
      94             :         integer(IK)             :: ndim                         !< The number of dimensions of the function.
      95             :         real(RK)                :: ftol = sqrt(epsilon(1._RK))  !< The stopping rule tolerance for the value of function.
      96             :         real(RK), allocatable   :: xmin(:)                      !< The x-value at the minimum of the function.
      97             :         real(RK), allocatable   :: DirMat(:,:)                  !< An initial `(ndim,ndim)` matrix whose columns contain the initial set of directions (usually the `ndim` unit vectors).
      98             :        !real(RK), allocatable   :: StartVec(:)                  !< An initial vector of size `ndim` representing the start of the search.
      99             :         real(RK)                :: fmin                         !< The minimum of the function.
     100             :         type(Err_type)          :: Err
     101             :     end type PowellMinimum_type
     102             : 
     103             :     ! overload Powell
     104             : 
     105             :     interface PowellMinimum_type
     106             :         module procedure :: minimizePowell
     107             :     end interface PowellMinimum_type
     108             : 
     109             :     abstract interface
     110             :         function getFuncMD_proc(ndim,Point) result(funcVal)
     111             :             import :: RK, IK
     112             :             integer(IK) , intent(in)    :: ndim
     113             :             real(RK)    , intent(in)    :: Point(ndim)
     114             :             real(RK)                    :: funcVal
     115             :         end function getFuncMD_proc
     116             :     end interface
     117             : 
     118             : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     119             : 
     120             :     interface minimize
     121             :         module procedure :: minimizeBrent, minimizePowell
     122             :     end interface minimize
     123             : 
     124             : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     125             : 
     126             : contains
     127             : 
     128             : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     129             : 
     130             :     !> The constructor of the class [BrentMinimum_type](@ref brentminimum_type).
     131             :     !> Compute the minimum of the input 1-dimensional function isolated to
     132             :     !> a fractional precision of about xtol using Brent's method.
     133             :     !>
     134             :     !> @param[in]   getFunc :   The 1-dimensional function which will have to be minimized.
     135             :     !> @param[in]   x0      :   The lower point in the set of optional bracketing triplet of abscissas that
     136             :     !>                          bracket the minimum of the function such that, `x0 < x1 < x2` and
     137             :     !>                          `getFunc(x0) > getFunc(x1) < getFunc(x2)` (**optional**).
     138             :     !> @param[in]   x1      :   The middle point in the set of optional bracketing triplet of abscissas that
     139             :     !>                          bracket the minimum of the function such that, `x0 < x1 < x2` and
     140             :     !>                          `getFunc(x0) > getFunc(x1) < getFunc(x2)` (**optional**).
     141             :     !> @param[in]   x2      :   The upper point in the set of optional bracketing triplet of abscissas that
     142             :     !>                          bracket the minimum of the function such that, `x0 < x1 < x2` and
     143             :     !>                          `getFunc(x0) > getFunc(x1) < getFunc(x2)` (**optional**).
     144             :     !> @param[in]   xtol    :   An optional fractional precision within which is the minimum is returned.
     145             :     !>                          The default value is sqrt(epsilon(1._RK)).
     146             :     !>
     147             :     !> \return
     148             :     !> `BrentMinimum` : An object of class [BrentMinimum_type](@ref brentminimum_type) that contains the minimum of the
     149             :     !>                  function (`xmin`) and the function value at the minimum (`fmin`) as well as other relevant information.
     150          56 :     function minimizeBrent(getFunc, x0, x1, x2, xtol) result(BrentMinimum)
     151             : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
     152             :         !DEC$ ATTRIBUTES DLLEXPORT :: minimizeBrent
     153             : #endif
     154             :         use Constants_mod, only: IK, RK
     155             :         procedure(getFunc1D_proc)           :: getFunc
     156             :         real(RK)    , intent(in), optional  :: x0, x1, x2
     157             :         real(RK)    , intent(in), optional  :: xtol
     158             :         type(BrentMinimum_type)             :: BrentMinimum
     159             : 
     160             :         character(*), parameter             :: PROCEDURE_NAME = MODULE_NAME//"@minimizeBrent"
     161             :         integer(IK) , parameter             :: ITMAX = 1000_IK                  ! the maximum number of iterations
     162             :         real(RK)    , parameter             :: CGOLD = 0.3819660_RK             ! Golden Section switch criterion
     163             :         real(RK)    , parameter             :: ZEPS = sqrt(epsilon(1._RK))**3   ! tiny nonzero value
     164             : 
     165             :         integer(IK) :: iter
     166             :         logical     :: isPresentX0, isPresentX1, isPresentX2
     167          56 :         real(RK)    :: a, b, d, e, etemp, fu, fv, fw, fx, p, q, r, tol1, tol2, u, v, w, x, xm
     168          56 :         real(RK)    :: ax, bx, cx, fa, fb, fc
     169             : 
     170          56 :         BrentMinimum%Err%occurred = .false.
     171             : 
     172          56 :         if (present(xtol)) BrentMinimum%xtol = xtol
     173             : 
     174          56 :         isPresentX0 = present(x0)
     175          56 :         isPresentX1 = present(x1)
     176          56 :         isPresentX2 = present(x2)
     177          56 :         if (isPresentX0 .and. isPresentX1 .and. isPresentX2) then
     178         220 :             BrentMinimum%Bracket = [x0, x1, x2]
     179             :         else
     180           1 :             if (isPresentX0) then
     181           0 :                 ax = x0
     182             :             else
     183           1 :                 ax = 0._RK ! assume an initial starting point of zero. This is the worst cae scenario.
     184             :             end if
     185           1 :             if (isPresentX1) then
     186           0 :                 bx = x1
     187             :             else
     188           1 :                 bx = ax + 1._RK
     189             :             end if
     190             :             call getBracket ( ax = ax &
     191             :                             , bx = bx &
     192             :                             , cx = cx &
     193             :                             , fa = fa &
     194             :                             , fb = fb &
     195             :                             , fc = fc &
     196             :                             , getFunc = getFunc &
     197           1 :                             )
     198           4 :             BrentMinimum%Bracket = [ax, bx, cx]
     199             :         end if
     200             : 
     201          56 :         a = min( BrentMinimum%Bracket(1), BrentMinimum%Bracket(3) )
     202          56 :         b = max( BrentMinimum%Bracket(1), BrentMinimum%Bracket(3) )
     203          56 :         v = BrentMinimum%Bracket(2)
     204          56 :         w = v
     205          56 :         x = v
     206          56 :         e = 0._RK
     207          56 :         fx = getFunc(x)
     208          56 :         fv = fx
     209          56 :         fw = fx
     210        1320 :         do iter = 1, ITMAX
     211        1320 :             BrentMinimum%niter = iter
     212        1320 :             xm = 0.5_RK * (a+b)
     213        1320 :             tol1 = xtol*abs(x) + ZEPS
     214        1320 :             tol2 = 2.0_RK*tol1
     215        1320 :             if (abs(x-xm) <= (tol2-0.5_RK*(b-a))) then
     216          56 :                 BrentMinimum%xmin = x
     217          56 :                 BrentMinimum%fmin = fx
     218          56 :                 return
     219             :             end if
     220        1264 :             if (abs(e) > tol1) then
     221        1194 :                 r = (x-w)*(fx-fv)
     222        1194 :                 q = (x-v)*(fx-fw)
     223        1194 :                 p = (x-v)*q-(x-w)*r
     224        1194 :                 q = 2.0_RK*(q-r)
     225        1194 :                 if (q > 0._RK) p = -p
     226        1194 :                 q = abs(q)
     227        1194 :                 etemp = e
     228        1194 :                 e = d
     229        1194 :                 if (abs(p) >= abs(0.5_RK*q*etemp) .or. p <= q*(a-x) .or. p >= q*(b-x)) then
     230         728 :                     e = merge(a-x,b-x, x >= xm )
     231         728 :                     d = cgold*e
     232             :                 else
     233         466 :                     d = p/q
     234         466 :                     u = x+d
     235         466 :                     if (u-a < tol2 .or. b-u < tol2) d = sign(tol1,xm-x)
     236             :                 end if
     237             :             else
     238          70 :                 e = merge(a-x,b-x, x >= xm )
     239          70 :                 d = cgold*e
     240             :             end if
     241        1264 :             u = merge(x+d,x+sign(tol1,d), abs(d) >= tol1 )
     242        1264 :             fu = getFunc(u)
     243        1264 :             if (fu <= fx) then
     244         357 :                 if (u >= x) then
     245         209 :                     a = x
     246             :                 else
     247         148 :                     b = x
     248             :                 end if
     249         357 :                 call shft(v,w,x,u)
     250         357 :                 call shft(fv,fw,fx,fu)
     251             :             else
     252         907 :                 if (u < x) then
     253         494 :                     a = u
     254             :                 else
     255         413 :                     b = u
     256             :                 end if
     257         907 :                 if (fu <= fw .or. w == x) then
     258         217 :                     v = w
     259         217 :                     fv = fw
     260         217 :                     w = u
     261         217 :                     fw = fu
     262         690 :                 else if (fu <= fv .or. v == x .or. v == w) then
     263         207 :                     v = u
     264         207 :                     fv = fu
     265             :                 end if
     266             :             end if
     267             :         end do
     268             : 
     269             :         ! LCOV_EXCL_START
     270             :         BrentMinimum%Err%occurred = .true.
     271             :         BrentMinimum%Err%msg = PROCEDURE_NAME//": maximum number of iterations exceeded."
     272             :         ! LCOV_EXCL_STOP
     273          56 :         return
     274             : 
     275             :     contains
     276             : 
     277         714 :         subroutine shft(a,b,c,d)
     278             :             implicit none
     279             :             real(RK), intent(out) :: a
     280             :             real(RK), intent(inout) :: b,c
     281             :             real(RK), intent(in) :: d
     282         714 :             a = b
     283         714 :             b = c
     284         714 :             c = d
     285          56 :         end subroutine shft
     286             : 
     287             :     end function minimizeBrent
     288             : 
     289             : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     290             : 
     291          55 :     subroutine getBracket(ax,bx,cx,fa,fb,fc,getFunc)
     292             : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
     293             :         !DEC$ ATTRIBUTES DLLEXPORT :: getBracket
     294             : #endif
     295             :         use Constants_mod, only: IK, RK, TINY_RK
     296             :         use Misc_mod, only : swap
     297             :         implicit none
     298             : 
     299             :         real(RK), intent(inout)     :: ax, bx
     300             :         real(RK), intent(out)       :: cx, fa, fb, fc
     301             :         procedure(getFunc1D_proc)   :: getFunc
     302             : 
     303             :         real(RK), parameter         :: GOLD = 1.618034_RK
     304             :         real(RK), parameter         :: GLIMIT = 100.0_RK
     305             :         real(RK), parameter         :: TINY = TINY_RK ! 1.0e-20_RK
     306          55 :         real(RK)                    :: fu, q, r, u, ulim
     307             : 
     308         110 :         fa = getFunc(ax)
     309          55 :         fb = getFunc(bx)
     310          55 :         if (fb > fa) then
     311          39 :             call swap(ax,bx)
     312          39 :             call swap(fa,fb)
     313             :         end if
     314             : 
     315          55 :         cx = bx + GOLD*(bx-ax)
     316          55 :         fc = getFunc(cx)
     317          64 :         do
     318          64 :             if (fb < fc) return
     319          15 :             r = (bx-ax)*(fb-fc)
     320          15 :             q = (bx-cx)*(fb-fa)
     321          15 :             u = bx-((bx-cx)*q-(bx-ax)*r)/(2._RK*sign(max(abs(q-r),TINY),q-r))
     322          15 :             ulim = bx+GLIMIT*(cx-bx)
     323          15 :             if ((bx-u)*(u-cx) > 0._RK) then
     324          10 :                 fu = getFunc(u)
     325          10 :                 if (fu < fc) then
     326           6 :                     ax = bx
     327           6 :                     fa = fb
     328           6 :                     bx = u
     329           6 :                     fb = fu
     330           6 :                     return
     331           4 :                 else if (fu > fb) then
     332           0 :                     cx = u
     333           0 :                     fc = fu
     334           0 :                     return
     335             :                 end if
     336           4 :                 u = cx+GOLD*(cx-bx)
     337           4 :                 fu = getFunc(u)
     338           5 :             else if ((cx-u)*(u-ulim) > 0._RK) then
     339           5 :                 fu = getFunc(u)
     340           5 :                 if (fu < fc) then
     341           5 :                     bx = cx
     342           5 :                     cx = u
     343           5 :                     u = cx+GOLD*(cx-bx)
     344           5 :                     call shft(fb,fc,fu,getFunc(u))
     345             :                 end if
     346           0 :             else if ((u-ulim)*(ulim-cx) >= 0._RK) then
     347           0 :                 u = ulim
     348           0 :                 fu = getFunc(u)
     349             :             else
     350           0 :                 u = cx+GOLD*(cx-bx)
     351           0 :                 fu = getFunc(u)
     352             :             end if
     353           9 :             call shft(ax,bx,cx,u)
     354           9 :             call shft(fa,fb,fc,fu)
     355             :         end do
     356             : 
     357             :     contains
     358             : 
     359          23 :         subroutine shft(a,b,c,d)
     360             :             real(RK), intent(out) :: a
     361             :             real(RK), intent(inout) :: b,c
     362             :             real(RK), intent(in) :: d
     363          23 :             a = b
     364          23 :             b = c
     365          23 :             c = d
     366          55 :         end subroutine shft
     367             : 
     368             :     end subroutine getBracket
     369             : 
     370             : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     371             : 
     372             :     !> The constructor of the class [PowellMinimum_type](@ref powellminimum_type).
     373           4 :     function minimizePowell(ndim, getFuncMD, StartVec, DirMat, ftol) result(PowellMinimum)
     374             : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
     375             :         !DEC$ ATTRIBUTES DLLEXPORT :: minimizePowell
     376             : #endif
     377             :         use Constants_mod, only: IK, RK, TINY_RK ! tiny = 1.0e-25_RK
     378             :         implicit none
     379             : 
     380             :         procedure(getFuncMD_proc)               :: getFuncMD
     381             :         integer(IK) , intent(in)                :: ndim
     382             :         real(RK)    , intent(in)                :: StartVec(ndim)
     383             :         real(RK)    , intent(in)    , optional  :: DirMat(ndim,ndim)
     384             :         real(RK)    , intent(in)    , optional  :: ftol
     385             :         type(PowellMinimum_type)                :: PowellMinimum
     386             : 
     387             :         character(*), parameter                 :: PROCEDURE_NAME = MODULE_NAME//"@minimizeBrent"
     388             :         integer(IK), parameter                  :: ITMAX = 1000
     389             : 
     390             :         integer(IK)                             :: i, ibig
     391           8 :         real(RK)                                :: del,fp,fptt,t
     392          28 :         real(RK)                                :: pt(ndim), ptt(ndim), xit(ndim)
     393             : 
     394           4 :         PowellMinimum%ndim = ndim
     395           4 :         PowellMinimum%Err%occurred = .false.
     396             : 
     397          12 :         allocate(PowellMinimum%xmin, source = StartVec)
     398             : 
     399           4 :         if (present(DirMat)) then
     400           0 :             PowellMinimum%DirMat = DirMat
     401             :         else
     402          28 :             allocate(PowellMinimum%DirMat(ndim,ndim), source = 0._RK)
     403          12 :             do i = 1, ndim
     404          12 :                 PowellMinimum%DirMat(i,i) = 1._RK
     405             :             end do
     406             :         end if
     407           4 :         if (present(ftol)) PowellMinimum%ftol = ftol
     408             : 
     409           4 :         PowellMinimum%fmin = getFuncMD(ndim,PowellMinimum%xmin)
     410          12 :         pt = PowellMinimum%xmin
     411           4 :         PowellMinimum%niter = 0
     412          10 :         do
     413             : 
     414          22 :             PowellMinimum%niter = PowellMinimum%niter + 1
     415          22 :             fp = PowellMinimum%fmin
     416          22 :             ibig = 0
     417          22 :             del = 0._RK
     418          66 :             do i = 1, ndim
     419         132 :                 xit = PowellMinimum%DirMat(1:ndim,i)
     420          44 :                 fptt = PowellMinimum%fmin
     421          44 :                 call linmin(getFuncMD, ndim, PowellMinimum%xmin, xit, PowellMinimum%fmin, PowellMinimum%Err)
     422          44 :                 if (PowellMinimum%Err%occurred) then
     423             :                 ! LCOV_EXCL_START
     424             :                     PowellMinimum%Err%msg = PROCEDURE_NAME//PowellMinimum%Err%msg
     425             :                     return
     426             :                 end if
     427             :                 ! LCOV_EXCL_STOP
     428         110 :                 if (fptt - PowellMinimum%fmin > del) then
     429          28 :                     del = fptt - PowellMinimum%fmin
     430          28 :                     ibig = i
     431             :                 end if
     432             :             end do
     433             : 
     434          22 :             if ( 2._RK*(fp-PowellMinimum%fmin) <= PowellMinimum%ftol*(abs(fp)+abs(PowellMinimum%fmin)) + TINY_RK ) return
     435             : 
     436          18 :             if (PowellMinimum%niter == ITMAX) then
     437             :             ! LCOV_EXCL_START
     438             :                 PowellMinimum%Err%occurred = .true.
     439             :                 PowellMinimum%Err%msg = PROCEDURE_NAME//": maximum number of iterations exceeded."
     440             :                 return
     441             :             end if
     442             :             ! LCOV_EXCL_STOP
     443             : 
     444          54 :             ptt = 2._RK * PowellMinimum%xmin - pt
     445          54 :             xit = PowellMinimum%xmin - pt
     446          54 :             pt = PowellMinimum%xmin
     447          18 :             fptt = getFuncMD(ndim,ptt)
     448          18 :             if (fptt >= fp) cycle
     449          14 :             t = 2._RK * (fp-2._RK*PowellMinimum%fmin+fptt) * (fp-PowellMinimum%fmin-del)**2 - del*(fp-fptt)**2
     450          14 :             if (t >= 0.0) cycle
     451          10 :             call linmin(getFuncMD, ndim, PowellMinimum%xmin, xit, PowellMinimum%fmin, PowellMinimum%Err)
     452          10 :             if (PowellMinimum%Err%occurred) then
     453             :             ! LCOV_EXCL_START
     454             :                 PowellMinimum%Err%msg = PROCEDURE_NAME//PowellMinimum%Err%msg
     455             :                 return
     456             :             end if
     457             :             ! LCOV_EXCL_STOP
     458          30 :             PowellMinimum%DirMat(1:ndim,ibig) = PowellMinimum%DirMat(1:ndim,ndim)
     459          30 :             PowellMinimum%DirMat(1:ndim,ndim) = xit
     460             : 
     461             :         end do
     462             : 
     463           4 :     end function minimizePowell
     464             : 
     465             : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     466             : 
     467          54 :     subroutine linmin(getFuncMD, ndim, StartVec, DirVec, fmin, Err)
     468             : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
     469             :         !DEC$ ATTRIBUTES DLLEXPORT :: linmin
     470             : #endif
     471           4 :         use Constants_mod, only: IK, RK
     472             :         use Err_mod, only: Err_type
     473             :         implicit none
     474             :         procedure(getFuncMD_proc)       :: getFuncMD
     475             :         integer(IK)     , intent(in)    :: ndim
     476             :         real(RK)        , intent(inout) :: StartVec(ndim), DirVec(ndim)
     477             :         real(RK)        , intent(out)   :: fmin
     478             :         type(Err_type)  , intent(out)   :: Err
     479             :         real(RK)        , parameter     :: XTOL = 1.0e-8_RK
     480          54 :         real(RK)                        :: ax, bx, fa, fb, fx, xx
     481          54 :         type(BrentMinimum_type)         :: BrentMinimum
     482             : #if defined OS_IS_WSL
     483             :         getFuncMD_WSL => getFuncMD
     484             :         DirVec_WSL = DirVec
     485             :         StartVec_WSL = StartVec
     486             :         ndim_WSL = ndim
     487             : #endif
     488          54 :         ax = 0.0
     489          54 :         xx = 1.0
     490          54 :         call getBracket(ax,xx,bx,fa,fx,fb,getFunc1D)
     491          54 :         BrentMinimum = minimizeBrent(getFunc1D, ax, xx, bx, XTOL)
     492          54 :         if (BrentMinimum%Err%occurred) then
     493             :         ! LCOV_EXCL_START
     494             :             Err = BrentMinimum%Err
     495             :             return
     496             :         ! LCOV_EXCL_STOP
     497             :         else
     498          54 :             Err%occurred = .false.
     499             :         end if
     500          54 :         fmin = BrentMinimum%fmin
     501         162 :         DirVec = BrentMinimum%xmin * DirVec
     502         162 :         StartVec = StartVec + DirVec
     503             : #if defined OS_IS_WSL
     504             :         nullify(getFuncMD_WSL)
     505             : #else
     506             :     contains
     507        1482 :         function getFunc1D(x) result(funcVal)
     508             :             implicit none
     509             :             real(RK), intent(in)    :: x
     510             :             real(RK)                :: funcVal
     511             :             real(RK), allocatable   :: xt(:)
     512        5927 :             xt = StartVec + x * DirVec
     513        1482 :             funcVal = getFuncMD(ndim,xt)
     514        1536 :         end function getFunc1D
     515             : #endif
     516             :     end subroutine linmin
     517             : 
     518             : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     519             : 
     520             : #if defined OS_IS_WSL
     521             :         !> \brief
     522             :         !> Bypass the Microsoft Subsystem for Linux Internal Function call GFortran Segmentation Fault error.
     523             :         function getFunc1D(x) result(funcVal)
     524             : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
     525             :         !DEC$ ATTRIBUTES DLLEXPORT :: getFunc1D
     526             : #endif
     527             :             implicit none
     528             :             real(RK), intent(in)    :: x
     529             :             real(RK)                :: funcVal
     530             :             real(RK), allocatable   :: xt(:)
     531             :             xt = StartVec_WSL + x * DirVec_WSL
     532             :             funcVal = getFuncMD_WSL(ndim_WSL,xt)
     533             :         end function getFunc1D
     534             : #endif
     535             : 
     536             : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     537             : 
     538             : end module Optimization_mod ! LCOV_EXCL_LINE

ParaMonte: Plain Powerful Parallel Monte Carlo Library 
The Computational Data Science Lab
© Copyright 2012 - 2021