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 2761 : 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 2761 : real(RK) :: a, b, d, e, etemp, fu, fv, fw, fx, p, q, r, tol1, tol2, u, v, w, x, xm
168 2761 : real(RK) :: ax, bx, cx, fa, fb, fc
169 :
170 2761 : BrentMinimum%Err%occurred = .false.
171 :
172 2761 : if (present(xtol)) BrentMinimum%xtol = xtol
173 :
174 2761 : isPresentX0 = present(x0)
175 2761 : isPresentX1 = present(x1)
176 2761 : isPresentX2 = present(x2)
177 2761 : if (isPresentX0 .and. isPresentX1 .and. isPresentX2) then
178 11032 : BrentMinimum%Bracket = [x0, x1, x2]
179 : else
180 3 : if (isPresentX0) then
181 0 : ax = x0
182 : else
183 3 : ax = 0._RK ! assume an initial starting point of zero. This is the worst cae scenario.
184 : end if
185 3 : if (isPresentX1) then
186 0 : bx = x1
187 : else
188 3 : 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 3 : )
198 12 : BrentMinimum%Bracket = [ax, bx, cx]
199 : end if
200 :
201 2761 : a = min( BrentMinimum%Bracket(1), BrentMinimum%Bracket(3) )
202 2761 : b = max( BrentMinimum%Bracket(1), BrentMinimum%Bracket(3) )
203 2761 : v = BrentMinimum%Bracket(2)
204 2761 : w = v
205 2761 : x = v
206 2761 : e = 0._RK
207 2761 : fx = getFunc(x)
208 2761 : fv = fx
209 2761 : fw = fx
210 62285 : do iter = 1, ITMAX
211 62285 : BrentMinimum%niter = iter
212 62285 : xm = 0.5_RK * (a+b)
213 62285 : tol1 = xtol*abs(x) + ZEPS
214 62285 : tol2 = 2.0_RK*tol1
215 62285 : if (abs(x-xm) <= (tol2-0.5_RK*(b-a))) then
216 2761 : BrentMinimum%xmin = x
217 2761 : BrentMinimum%fmin = fx
218 2761 : return
219 : end if
220 59524 : if (abs(e) > tol1) then
221 56022 : r = (x-w)*(fx-fv)
222 56022 : q = (x-v)*(fx-fw)
223 56022 : p = (x-v)*q-(x-w)*r
224 56022 : q = 2.0_RK*(q-r)
225 56022 : if (q > 0._RK) p = -p
226 56022 : q = abs(q)
227 56022 : etemp = e
228 56022 : e = d
229 56022 : if (abs(p) >= abs(0.5_RK*q*etemp) .or. p <= q*(a-x) .or. p >= q*(b-x)) then
230 35151 : e = merge(a-x,b-x, x >= xm )
231 35151 : d = cgold*e
232 : else
233 20871 : d = p/q
234 20871 : u = x+d
235 20871 : if (u-a < tol2 .or. b-u < tol2) d = sign(tol1,xm-x)
236 : end if
237 : else
238 3502 : e = merge(a-x,b-x, x >= xm )
239 3502 : d = cgold*e
240 : end if
241 59524 : u = merge(x+d,x+sign(tol1,d), abs(d) >= tol1 )
242 59524 : fu = getFunc(u)
243 59524 : if (fu <= fx) then
244 21682 : if (u >= x) then
245 10082 : a = x
246 : else
247 11600 : b = x
248 : end if
249 21682 : call shft(v,w,x,u)
250 21682 : call shft(fv,fw,fx,fu)
251 : else
252 37842 : if (u < x) then
253 17590 : a = u
254 : else
255 20252 : b = u
256 : end if
257 37842 : if (fu <= fw .or. w == x) then
258 11713 : v = w
259 11713 : fv = fw
260 11713 : w = u
261 11713 : fw = fu
262 26129 : else if (fu <= fv .or. v == x .or. v == w) then
263 6780 : v = u
264 6780 : 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 2761 : return
274 :
275 : contains
276 :
277 43364 : 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 43364 : a = b
283 43364 : b = c
284 43364 : c = d
285 2761 : end subroutine shft
286 :
287 : end function minimizeBrent
288 :
289 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
290 :
291 2758 : 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 2758 : real(RK) :: fu, q, r, u, ulim
307 :
308 5516 : fa = getFunc(ax)
309 2758 : fb = getFunc(bx)
310 2758 : if (fb > fa) then
311 2269 : call swap(ax,bx)
312 2269 : call swap(fa,fb)
313 : end if
314 :
315 2758 : cx = bx + GOLD*(bx-ax)
316 2758 : fc = getFunc(cx)
317 3006 : do
318 3006 : if (fb < fc) return
319 387 : r = (bx-ax)*(fb-fc)
320 387 : q = (bx-cx)*(fb-fa)
321 387 : u = bx-((bx-cx)*q-(bx-ax)*r)/(2._RK*sign(max(abs(q-r),TINY),q-r))
322 387 : ulim = bx+GLIMIT*(cx-bx)
323 387 : if ((bx-u)*(u-cx) > 0._RK) then
324 159 : fu = getFunc(u)
325 159 : if (fu < fc) then
326 139 : ax = bx
327 139 : fa = fb
328 139 : bx = u
329 139 : fb = fu
330 139 : return
331 20 : else if (fu > fb) then
332 0 : cx = u
333 0 : fc = fu
334 0 : return
335 : end if
336 20 : u = cx+GOLD*(cx-bx)
337 20 : fu = getFunc(u)
338 228 : else if ((cx-u)*(u-ulim) > 0._RK) then
339 227 : fu = getFunc(u)
340 227 : if (fu < fc) then
341 227 : bx = cx
342 227 : cx = u
343 227 : u = cx+GOLD*(cx-bx)
344 227 : call shft(fb,fc,fu,getFunc(u))
345 : end if
346 1 : else if ((u-ulim)*(ulim-cx) >= 0._RK) then
347 1 : u = ulim
348 1 : fu = getFunc(u)
349 : else
350 0 : u = cx+GOLD*(cx-bx)
351 0 : fu = getFunc(u)
352 : end if
353 248 : call shft(ax,bx,cx,u)
354 248 : call shft(fa,fb,fc,fu)
355 : end do
356 :
357 : contains
358 :
359 723 : 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 723 : a = b
364 723 : b = c
365 723 : c = d
366 2758 : end subroutine shft
367 :
368 : end subroutine getBracket
369 :
370 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
371 :
372 : !> The constructor of the class [PowellMinimum_type](@ref powellminimum_type).
373 221 : 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 442 : real(RK) :: del,fp,fptt,t
392 1547 : real(RK) :: pt(ndim), ptt(ndim), xit(ndim)
393 :
394 221 : PowellMinimum%ndim = ndim
395 221 : PowellMinimum%Err%occurred = .false.
396 :
397 663 : allocate(PowellMinimum%xmin, source = StartVec)
398 :
399 221 : if (present(DirMat)) then
400 0 : PowellMinimum%DirMat = DirMat
401 : else
402 1547 : allocate(PowellMinimum%DirMat(ndim,ndim), source = 0._RK)
403 663 : do i = 1, ndim
404 663 : PowellMinimum%DirMat(i,i) = 1._RK
405 : end do
406 : end if
407 221 : if (present(ftol)) PowellMinimum%ftol = ftol
408 :
409 221 : PowellMinimum%fmin = getFuncMD(ndim,PowellMinimum%xmin)
410 663 : pt = PowellMinimum%xmin
411 221 : PowellMinimum%niter = 0
412 345 : do
413 :
414 1205 : PowellMinimum%niter = PowellMinimum%niter + 1
415 1205 : fp = PowellMinimum%fmin
416 1205 : ibig = 0
417 1205 : del = 0._RK
418 3615 : do i = 1, ndim
419 7230 : xit = PowellMinimum%DirMat(1:ndim,i)
420 2410 : fptt = PowellMinimum%fmin
421 2410 : call linmin(getFuncMD, ndim, PowellMinimum%xmin, xit, PowellMinimum%fmin, PowellMinimum%Err)
422 2410 : 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 6025 : if (fptt - PowellMinimum%fmin > del) then
429 1383 : del = fptt - PowellMinimum%fmin
430 1383 : ibig = i
431 : end if
432 : end do
433 :
434 1205 : if ( 2._RK*(fp-PowellMinimum%fmin) <= PowellMinimum%ftol*(abs(fp)+abs(PowellMinimum%fmin)) + TINY_RK ) return
435 :
436 984 : 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 2952 : ptt = 2._RK * PowellMinimum%xmin - pt
445 2952 : xit = PowellMinimum%xmin - pt
446 2952 : pt = PowellMinimum%xmin
447 984 : fptt = getFuncMD(ndim,ptt)
448 984 : if (fptt >= fp) cycle
449 443 : t = 2._RK * (fp-2._RK*PowellMinimum%fmin+fptt) * (fp-PowellMinimum%fmin-del)**2 - del*(fp-fptt)**2
450 443 : if (t >= 0.0) cycle
451 345 : call linmin(getFuncMD, ndim, PowellMinimum%xmin, xit, PowellMinimum%fmin, PowellMinimum%Err)
452 345 : 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 1035 : PowellMinimum%DirMat(1:ndim,ibig) = PowellMinimum%DirMat(1:ndim,ndim)
459 1035 : PowellMinimum%DirMat(1:ndim,ndim) = xit
460 :
461 : end do
462 :
463 221 : end function minimizePowell
464 :
465 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
466 :
467 2755 : 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 221 : 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 2755 : real(RK) :: ax, bx, fa, fb, fx, xx
481 2755 : 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 2755 : ax = 0.0
489 2755 : xx = 1.0
490 2755 : call getBracket(ax,xx,bx,fa,fx,fb,getFunc1D)
491 2755 : BrentMinimum = minimizeBrent(getFunc1D, ax, xx, bx, XTOL)
492 2755 : if (BrentMinimum%Err%occurred) then
493 : ! LCOV_EXCL_START
494 : Err = BrentMinimum%Err
495 : return
496 : ! LCOV_EXCL_STOP
497 : else
498 2755 : Err%occurred = .false.
499 : end if
500 2755 : fmin = BrentMinimum%fmin
501 8265 : DirVec = BrentMinimum%xmin * DirVec
502 8265 : StartVec = StartVec + DirVec
503 : #if defined OS_IS_WSL
504 : nullify(getFuncMD_WSL)
505 : #else
506 : contains
507 71112 : function getFunc1D(x) result(funcVal)
508 : implicit none
509 : real(RK), intent(in) :: x
510 : real(RK) :: funcVal
511 : real(RK), allocatable :: xt(:)
512 284445 : xt = StartVec + x * DirVec
513 71112 : funcVal = getFuncMD(ndim,xt)
514 73867 : 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
|
