Line data Source code
1 : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2 : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
3 : !!!! !!!!
4 : !!!! ParaMonte: Parallel Monte Carlo and Machine Learning Library. !!!!
5 : !!!! !!!!
6 : !!!! Copyright (C) 2012-present, The Computational Data Science Lab !!!!
7 : !!!! !!!!
8 : !!!! This file is part of the ParaMonte library. !!!!
9 : !!!! !!!!
10 : !!!! LICENSE !!!!
11 : !!!! !!!!
12 : !!!! https://github.com/cdslaborg/paramonte/blob/main/LICENSE.md !!!!
13 : !!!! !!!!
14 : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
15 : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
16 :
17 : !> \brief
18 : !> This file contains the implementations of the tests of module [pm_mathRoot](@ref pm_mathRoot).
19 : !>
20 : !> \fintest
21 : !>
22 : !> \author
23 : !> \AmirShahmoradi, Sunday 4:33 PM, September 19, 2021, Dallas, TX
24 :
25 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
26 :
27 : #if CK_ENABLED
28 : #define TYPE_KIND complex(TKC)
29 : #elif RK_ENABLED
30 : #define TYPE_KIND real(TKC)
31 : #else
32 : #error "Unrecognized interface."
33 : #endif
34 : ! Select root-finding method.
35 : #if Brent_ENABLED || Def_ENABLED
36 : #define METHOD brent_type
37 : #elif False_ENABLED
38 : #define METHOD false_type
39 : #elif Secant_ENABLED
40 : #define METHOD secant_type
41 : #elif Newton_ENABLED
42 : #define METHOD newton_type
43 : #elif Halley_ENABLED
44 : #define METHOD halley_type
45 : #elif Ridders_ENABLED
46 : #define METHOD ridders_type
47 : #elif TOMS748_ENABLED
48 : #define METHOD toms748_type
49 : #elif Schroder_ENABLED
50 : #define METHOD schroder_type
51 : #elif Bisection_ENABLED
52 : #define METHOD bisection_type
53 : #else
54 : #error "Unrecognized interface."
55 : #endif
56 : ! Set the default `getFunc()`.
57 : #if Newton_ENABLED || Halley_ENABLED || Schroder_ENABLED
58 : #define GETFUNC getFuncDiff
59 : #endif
60 : type(METHOD), parameter :: method = METHOD()
61 76 : type(display_type) :: disp
62 : integer(IK) :: itry, neval, niter
63 : real(TKC), parameter :: abstol_ref = epsilon(0._TKC)**.8
64 : TYPE_KIND :: roots(2), lb, ub, xmid
65 76 : assertion = .true._LK
66 :
67 7752 : do itry = 1, 100
68 :
69 22800 : roots = getChoice(getRemoved(getRange(-5, 5), 0), 2_IK, unique = .true._LK)
70 45600 : lb = getChoice([minval(roots, 1), minval(roots, 1) - 1])
71 45600 : ub = getChoice([maxval(roots, 1), maxval(roots, 1) + 1])
72 22800 : if (roots .allinrange. [lb, ub]) then
73 22800 : xmid = .5_TKC * sum(roots(1:2))
74 7600 : if (getUnifRand()) then
75 3831 : lb = xmid
76 : else
77 3769 : ub = xmid
78 : end if
79 : end if
80 7600 : niter = getUnifRand(100, 300)
81 :
82 : #if getRoot_ENABLED
83 4000 : call testWith(method)
84 4000 : call testWith(method, abstol_ref)
85 4000 : call testWith(method, abstol_ref, neval)
86 4000 : call testWith(method, abstol_ref, neval, niter)
87 4000 : call testWith(method, abstol_ref, niter = niter)
88 4000 : call testWith(method, neval = neval, niter = niter)
89 4000 : call testWith(method, neval = neval)
90 4000 : call testWith(method, niter = niter)
91 4000 : call testWith()
92 4000 : call testWith(abstol = abstol_ref)
93 4000 : call testWith(abstol = abstol_ref, neval = neval)
94 4000 : call testWith(abstol = abstol_ref, neval = neval, niter = niter)
95 4000 : call testWith(abstol = abstol_ref, niter = niter)
96 4000 : call testWith(neval = neval, niter = niter)
97 4000 : call testWith(neval = neval)
98 4040 : call testWith(niter = niter)
99 : #elif setRoot_ENABLED
100 36 : block
101 : TYPE_KIND :: root, lf, uf
102 3600 : lf = getFunc(lb)
103 3600 : uf = getFunc(ub)
104 3600 : root = getUnifRand(lb - 1, ub + 1) ! this initialization is required in testing of the Newton method.
105 3600 : call setRoot(method, GETFUNC, root, lb, ub, lf, uf, abstol_ref, neval)
106 5375 : assertion = assertion .and. (any(abs(roots - root) < abstol_ref * 1000) .or. getOption(0_IK, neval) < 0_IK)
107 3600 : call report(__LINE__, root, method, abstol_ref * 1000, neval)
108 :
109 3600 : call setRoot(method, GETFUNC, root, lb, ub, lf, uf, abstol_ref, neval, niter)
110 5375 : assertion = assertion .and. (any(abs(roots - root) < abstol_ref * 1000) .or. getOption(0_IK, neval) < 0_IK)
111 3600 : call report(__LINE__, root, method, abstol_ref * 1000, neval, niter)
112 : end block
113 : #else
114 : #error "Unrecognized interface."
115 : #endif
116 :
117 : end do
118 :
119 : contains
120 :
121 1021708 : pure function getFunc(x) result(func)
122 : real(TKC), intent(in) :: x
123 : real(TKC) :: func
124 6130248 : func = x**2 - sum(roots(1:2)) * x + product(roots) ! (x - roots(1)) * (x - roots(2))
125 1021708 : end function
126 :
127 : #if Newton_ENABLED || Halley_ENABLED || Schroder_ENABLED
128 199113 : pure function getFuncDiff(x, order) result(func)
129 : integer(IK), intent(in) :: order
130 : real(TKC), intent(in) :: x
131 : real(TKC) :: func
132 199113 : if (order == 0) func = getFunc(x)
133 290531 : if (order == 1) func = 2._TKC * x - sum(roots(1:2))
134 199113 : if (order == 2) func = 2._TKC
135 199113 : end function
136 : #endif
137 :
138 : #if getRoot_ENABLED
139 64000 : subroutine testWith(method, abstol, neval, niter)
140 : type(METHOD), intent(in), optional :: method
141 : integer(IK), intent(inout), optional :: neval, niter
142 : real(TKC), intent(in), optional :: abstol
143 : real(TKC) :: abstol_def
144 : integer(IK) :: neval_def
145 : TYPE_KIND :: root
146 64000 : if (present(abstol)) then
147 32000 : abstol_def = abstol * 1000._TKC
148 : else
149 32000 : abstol_def = abstol_ref * (abs(lb) + abs(ub)) * 1000._TKC
150 : end if
151 64000 : if (present(method)) then
152 32000 : root = getRoot(method, GETFUNC, lb, ub, abstol, neval, niter)
153 47968 : assertion = assertion .and. (any(abs(roots - root) < abstol_def) .or. getOption(0_IK, neval) < 0_IK)
154 32000 : call report(__LINE__, root, method, abstol, neval, niter)
155 : #if Newton_ENABLED || Halley_ENABLED || Schroder_ENABLED
156 : !call disp%show("[real(TKC) :: lb, ub, getFunc(lb), getFunc(ub), abstol_def]")
157 : !call disp%show( [real(TKC) :: lb, ub, getFunc(lb), getFunc(ub), abstol_def] )
158 : !call disp%show("[roots, getFunc(roots(1)), getFunc(roots(2))]")
159 : !call disp%show( [roots, getFunc(roots(1)), getFunc(roots(2))] )
160 9600 : root = getRoot(method, GETFUNC, lb, ub, abstol, neval, niter, init = getUnifRand(lb - 1, ub + 1))
161 14408 : assertion = assertion .and. (any(abs(roots - root) < abstol_def) .or. getOption(0_IK, neval) < 0_IK)
162 9600 : call report(__LINE__, root, method, abstol, neval, niter)
163 : #endif
164 : else
165 : !print *, "size(roots)", size(roots)
166 : !print *, "roots", roots
167 : !print *, "froots", getFunc(roots(0)), getFunc(roots(1)), getFunc(roots(2))
168 : !print *, "lb, ub", lb, ub
169 32000 : root = getRoot(getFunc, lb, ub, abstol, neval, niter)
170 47968 : assertion = assertion .and. (any(abs(roots - root) < abstol_def) .or. getOption(0_IK, neval) < 0_IK)
171 32000 : call report(__LINE__, root, method, abstol, neval, niter)
172 : end if
173 64000 : root = getRoot(getFunc, lb, ub, abstol, neval_def, 1_IK)
174 134568 : assertion = assertion .and. (any(abs(roots - root) < abstol_def) .or. neval_def <= 0_IK)
175 64000 : call report(__LINE__, root, method, abstol, neval_def, niter)
176 64000 : end subroutine
177 : #endif
178 :
179 144800 : subroutine report(line, root, method, abstol, neval, niter)
180 : type(METHOD), intent(in), optional :: method
181 : integer(IK), intent(in), optional :: neval, niter
182 : real(TKC), intent(in), optional :: abstol
183 : TYPE_KIND, intent(in) :: root
184 : integer, intent(in) :: line
185 144800 : if (test%traceable .and. .not. assertion) then
186 : ! LCOV_EXCL_START
187 : call disp%skip()
188 : call disp%show("roots")
189 : call disp%show( roots )
190 : call disp%show("[lb, root, ub]")
191 : call disp%show( [lb, root, ub] )
192 : call disp%show("getFunc(root)")
193 : call disp%show( getFunc(root) )
194 : call disp%show("abs(root - roots)")
195 : call disp%show( abs(root - roots) )
196 : call disp%show("present(method)")
197 : call disp%show( present(method) )
198 : call disp%show("present(abstol)")
199 : call disp%show( present(abstol) )
200 : if (present(abstol)) then
201 : call disp%show("abstol")
202 : call disp%show( abstol )
203 : end if
204 : call disp%show("present(neval)")
205 : call disp%show( present(neval) )
206 : if (present(neval)) then
207 : call disp%show("neval")
208 : call disp%show( neval )
209 : end if
210 : call disp%show("present(niter)")
211 : call disp%show( present(niter) )
212 : if (present(niter)) then
213 : call disp%show("niter")
214 : call disp%show( niter )
215 : end if
216 : ! LCOV_EXCL_STOP
217 : end if
218 144800 : call test%assert(assertion, SK_"getRoot() must compute a root of the target function correctly.", int(line, IK))
219 144800 : end subroutine
220 :
221 : #undef GETFUNC
222 : #undef METHOD
|
