Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2024 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 :
8 : ! **************************************************************************************************
9 : !> \brief This public domain function parser module is intended for applications
10 : !> where a set of mathematical expressions is specified at runtime and is
11 : !> then evaluated for a large number of variable values. This is done by
12 : !> compiling the set of function strings into byte code, which is interpreted
13 : !> very efficiently for the various variable values.
14 : !>
15 : !> \author Roland Schmehl <Roland.Schmehl@mach.uni-karlsruhe.de>
16 : ! **************************************************************************************************
17 : MODULE fparser
18 : USE kinds, ONLY: default_string_length,&
19 : rn => dp
20 : #include "../base/base_uses.f90"
21 :
22 : IMPLICIT NONE
23 :
24 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'fparser'
25 :
26 : !------- -------- --------- --------- --------- --------- --------- --------- -------
27 : PUBLIC :: initf, & ! Initialize function parser for n functions
28 : parsef, & ! Parse single function string
29 : evalf, & ! Evaluate single function
30 : finalizef, & ! Finalize the function parser
31 : evalfd
32 : INTEGER, PUBLIC :: EvalErrType ! =0: no error occurred, >0: evaluation error
33 : !------- -------- --------- --------- --------- --------- --------- --------- -------
34 : PRIVATE
35 : SAVE
36 : INTEGER, PARAMETER, PRIVATE :: is = SELECTED_INT_KIND(1) !Data type of bytecode
37 : INTEGER(is), PARAMETER :: cImmed = 1, &
38 : cNeg = 2, &
39 : cAdd = 3, &
40 : cSub = 4, &
41 : cMul = 5, &
42 : cDiv = 6, &
43 : cPow = 7, &
44 : cAbs = 8, &
45 : cExp = 9, &
46 : cLog10 = 10, &
47 : cLog = 11, &
48 : cSqrt = 12, &
49 : cSinh = 13, &
50 : cCosh = 14, &
51 : cTanh = 15, &
52 : cSin = 16, &
53 : cCos = 17, &
54 : cTan = 18, &
55 : cAsin = 19, &
56 : cAcos = 20, &
57 : cAtan = 21, &
58 : cErf = 22, &
59 : cErfc = 23, &
60 : VarBegin = 24
61 : CHARACTER(LEN=1), DIMENSION(cAdd:cPow), PARAMETER :: Ops = (/'+', &
62 : '-', &
63 : '*', &
64 : '/', &
65 : '^'/)
66 : CHARACTER(LEN=5), DIMENSION(cAbs:cErfc), PARAMETER :: Funcs = (/'abs ', &
67 : 'exp ', &
68 : 'log10', &
69 : 'log ', &
70 : 'sqrt ', &
71 : 'sinh ', &
72 : 'cosh ', &
73 : 'tanh ', &
74 : 'sin ', &
75 : 'cos ', &
76 : 'tan ', &
77 : 'asin ', &
78 : 'acos ', &
79 : 'atan ', &
80 : 'erf ', &
81 : 'erfc '/)
82 : ! **************************************************************************************************
83 : TYPE tComp
84 : INTEGER(is), DIMENSION(:), POINTER :: ByteCode => NULL()
85 : INTEGER :: ByteCodeSize = -1
86 : REAL(rn), DIMENSION(:), POINTER :: Immed => NULL()
87 : INTEGER :: ImmedSize = -1
88 : REAL(rn), DIMENSION(:), POINTER :: Stack => NULL()
89 : INTEGER :: StackSize = -1, &
90 : StackPtr = -1
91 : END TYPE tComp
92 : TYPE(tComp), DIMENSION(:), POINTER :: Comp => NULL() ! Bytecode
93 : INTEGER, DIMENSION(:), ALLOCATABLE :: ipos ! Associates function strings
94 : !
95 : CONTAINS
96 : !
97 : ! **************************************************************************************************
98 : !> \brief ...
99 : ! **************************************************************************************************
100 38724 : SUBROUTINE finalizef()
101 : !----- -------- --------- --------- --------- --------- --------- --------- -------
102 : ! Finalize function parser
103 : !----- -------- --------- --------- --------- --------- --------- --------- -------
104 :
105 : INTEGER :: i
106 :
107 : !----- -------- --------- --------- --------- --------- --------- --------- -------
108 :
109 64335 : DO i = 1, SIZE(Comp)
110 25611 : IF (ASSOCIATED(Comp(i)%ByteCode)) THEN
111 25611 : DEALLOCATE (Comp(i)%ByteCode)
112 : END IF
113 25611 : IF (ASSOCIATED(Comp(i)%Immed)) THEN
114 25611 : DEALLOCATE (Comp(i)%Immed)
115 : END IF
116 64335 : IF (ASSOCIATED(Comp(i)%Stack)) THEN
117 25611 : DEALLOCATE (Comp(i)%Stack)
118 : END IF
119 : END DO
120 :
121 38724 : DEALLOCATE (Comp)
122 :
123 38724 : END SUBROUTINE finalizef
124 : !
125 : ! **************************************************************************************************
126 : !> \brief ...
127 : !> \param n ...
128 : ! **************************************************************************************************
129 38724 : SUBROUTINE initf(n)
130 : !----- -------- --------- --------- --------- --------- --------- --------- -------
131 : ! Initialize function parser for n functions
132 : !----- -------- --------- --------- --------- --------- --------- --------- -------
133 : INTEGER, INTENT(in) :: n
134 :
135 : ! Number of functions
136 : !----- -------- --------- --------- --------- --------- --------- --------- -------
137 :
138 108742 : ALLOCATE (Comp(n))
139 38724 : END SUBROUTINE initf
140 : !
141 : ! **************************************************************************************************
142 : !> \brief Parse ith function string FuncStr and compile it into bytecode
143 : !> \param i Function identifier
144 : !> \param FuncStr Function string
145 : !> \param Var Array with variable names
146 : ! **************************************************************************************************
147 51222 : SUBROUTINE parsef(i, FuncStr, Var)
148 : INTEGER, INTENT(in) :: i
149 : CHARACTER(LEN=*), INTENT(in) :: FuncStr
150 : CHARACTER(LEN=*), DIMENSION(:), INTENT(in) :: Var
151 :
152 25611 : CHARACTER(LEN=LEN(FuncStr)) :: Func
153 :
154 : ! Function string, local use
155 : !----- -------- --------- --------- --------- --------- --------- --------- -------
156 :
157 25611 : IF (i < 1 .OR. i > SIZE(Comp)) THEN
158 0 : CPABORT("Function number is out of range")
159 : END IF
160 25611 : IF (SIZE(var) > HUGE(0_is)) THEN
161 0 : CPABORT("Too many variables")
162 : END IF
163 76833 : ALLOCATE (ipos(LEN_TRIM(FuncStr))) ! Char. positions in orig. string
164 25611 : Func = FuncStr ! Local copy of function string
165 25611 : CALL Replace('**', '^ ', Func) ! Exponent into 1-Char. format
166 25611 : CALL RemoveSpaces(Func) ! Condense function string
167 25611 : CALL CheckSyntax(Func, FuncStr, Var)
168 25611 : DEALLOCATE (ipos)
169 25611 : CALL Compile(i, Func, Var) ! Compile into bytecode
170 :
171 25611 : END SUBROUTINE parsef
172 : !
173 : ! **************************************************************************************************
174 : !> \brief ...
175 : !> \param i ...
176 : !> \param Val ...
177 : !> \return ...
178 : ! **************************************************************************************************
179 51473349 : FUNCTION evalf(i, Val) RESULT(res)
180 : !----- -------- --------- --------- --------- --------- --------- --------- -------
181 : ! Evaluate bytecode of ith function for the values passed in array Val(:)
182 : !----- -------- --------- --------- --------- --------- --------- --------- -------
183 : INTEGER, INTENT(in) :: i
184 : REAL(rn), DIMENSION(:), INTENT(in) :: Val
185 : REAL(rn) :: res
186 :
187 : REAL(rn), PARAMETER :: zero = 0._rn
188 :
189 : INTEGER :: DP, IP, ipow, SP
190 :
191 : ! Function identifier
192 : ! Variable values
193 : ! Result
194 : ! Instruction pointer
195 : ! Data pointer
196 : ! Stack pointer
197 : !----- -------- --------- --------- --------- --------- --------- --------- -------
198 :
199 51473349 : DP = 1
200 51473349 : SP = 0
201 1002322792 : DO IP = 1, Comp(i)%ByteCodeSize
202 51473339 : SELECT CASE (Comp(i)%ByteCode(IP))
203 :
204 150630263 : CASE (cImmed); SP = SP + 1; Comp(i)%Stack(SP) = Comp(i)%Immed(DP); DP = DP + 1
205 2762941 : CASE (cNeg); Comp(i)%Stack(SP) = -Comp(i)%Stack(SP)
206 97288529 : CASE (cAdd); Comp(i)%Stack(SP - 1) = Comp(i)%Stack(SP - 1) + Comp(i)%Stack(SP); SP = SP - 1
207 49374414 : CASE (cSub); Comp(i)%Stack(SP - 1) = Comp(i)%Stack(SP - 1) - Comp(i)%Stack(SP); SP = SP - 1
208 100721296 : CASE (cMul); Comp(i)%Stack(SP - 1) = Comp(i)%Stack(SP - 1)*Comp(i)%Stack(SP); SP = SP - 1
209 98846062 : CASE (cDiv); IF (Comp(i)%Stack(SP) == 0._rn) THEN; EvalErrType = 1; res = zero; RETURN; END IF
210 98846052 : Comp(i)%Stack(SP - 1) = Comp(i)%Stack(SP - 1)/Comp(i)%Stack(SP); SP = SP - 1
211 : CASE (cPow)
212 : ! Fixing for possible Negative floating-point value raised to a real power
213 100777224 : IF (Comp(i)%Stack(SP - 1) < 0.0_rn) THEN
214 321254 : ipow = FLOOR(Comp(i)%Stack(SP))
215 321254 : IF (MOD(Comp(i)%Stack(SP), REAL(ipow, KIND=rn)) == 0.0_rn) THEN
216 321254 : Comp(i)%Stack(SP - 1) = Comp(i)%Stack(SP - 1)**ipow
217 : ELSE
218 0 : CPABORT("Negative floating-point value raised to a real power!")
219 : END IF
220 : ELSE
221 100455970 : Comp(i)%Stack(SP - 1) = Comp(i)%Stack(SP - 1)**Comp(i)%Stack(SP)
222 : END IF
223 0 : SP = SP - 1
224 0 : CASE (cAbs); Comp(i)%Stack(SP) = ABS(Comp(i)%Stack(SP))
225 2498087 : CASE (cExp); Comp(i)%Stack(SP) = EXP(Comp(i)%Stack(SP))
226 0 : CASE (cLog10); IF (Comp(i)%Stack(SP) <= 0._rn) THEN; EvalErrType = 3; res = zero; RETURN; END IF
227 0 : Comp(i)%Stack(SP) = LOG10(Comp(i)%Stack(SP))
228 0 : CASE (cLog); IF (Comp(i)%Stack(SP) <= 0._rn) THEN; EvalErrType = 3; res = zero; RETURN; END IF
229 0 : Comp(i)%Stack(SP) = LOG(Comp(i)%Stack(SP))
230 0 : CASE (cSqrt); IF (Comp(i)%Stack(SP) < 0._rn) THEN; EvalErrType = 3; res = zero; RETURN; END IF
231 0 : Comp(i)%Stack(SP) = SQRT(Comp(i)%Stack(SP))
232 0 : CASE (cSinh); Comp(i)%Stack(SP) = SINH(Comp(i)%Stack(SP))
233 0 : CASE (cCosh); Comp(i)%Stack(SP) = COSH(Comp(i)%Stack(SP))
234 0 : CASE (cTanh); Comp(i)%Stack(SP) = TANH(Comp(i)%Stack(SP))
235 84816 : CASE (cSin); Comp(i)%Stack(SP) = SIN(Comp(i)%Stack(SP))
236 13206 : CASE (cCos); Comp(i)%Stack(SP) = COS(Comp(i)%Stack(SP))
237 0 : CASE (cTan); Comp(i)%Stack(SP) = TAN(Comp(i)%Stack(SP))
238 0 : CASE (cAsin); IF ((Comp(i)%Stack(SP) < -1._rn) .OR. (Comp(i)%Stack(SP) > 1._rn)) THEN
239 0 : EvalErrType = 4; res = zero; RETURN; END IF
240 0 : Comp(i)%Stack(SP) = ASIN(Comp(i)%Stack(SP))
241 1548 : CASE (cAcos); IF ((Comp(i)%Stack(SP) < -1._rn) .OR. (Comp(i)%Stack(SP) > 1._rn)) THEN
242 0 : EvalErrType = 4; res = zero; RETURN; END IF
243 1548 : Comp(i)%Stack(SP) = ACOS(Comp(i)%Stack(SP))
244 456 : CASE (cAtan); Comp(i)%Stack(SP) = ATAN(Comp(i)%Stack(SP))
245 0 : CASE (cErf); Comp(i)%Stack(SP) = ERF(Comp(i)%Stack(SP))
246 0 : CASE (cErfc); Comp(i)%Stack(SP) = ERFC(Comp(i)%Stack(SP))
247 950849453 : CASE DEFAULT; SP = SP + 1; Comp(i)%Stack(SP) = Val(Comp(i)%ByteCode(IP) - VarBegin + 1)
248 : END SELECT
249 : END DO
250 51473339 : EvalErrType = 0
251 51473339 : res = Comp(i)%Stack(1)
252 51473339 : END FUNCTION evalf
253 : !
254 : ! **************************************************************************************************
255 : !> \brief ...
256 : !> \param Func ...
257 : !> \param FuncStr ...
258 : !> \param Var ...
259 : ! **************************************************************************************************
260 25611 : SUBROUTINE CheckSyntax(Func, FuncStr, Var)
261 : !----- -------- --------- --------- --------- --------- --------- --------- -------
262 : ! Check syntax of function string, returns 0 if syntax is ok
263 : !----- -------- --------- --------- --------- --------- --------- --------- -------
264 : CHARACTER(LEN=*), INTENT(in) :: Func, FuncStr
265 : CHARACTER(LEN=*), DIMENSION(:), INTENT(in) :: Var
266 :
267 : INTEGER :: ib, in, j, lFunc, ParCnt
268 : CHARACTER(LEN=1) :: c
269 : INTEGER(is) :: n
270 : LOGICAL :: err
271 : REAL(rn) :: r
272 :
273 : ! Function string without spaces
274 : ! Original function string
275 : ! Array with variable names
276 : ! Parenthesis counter
277 : !----- -------- --------- --------- --------- --------- --------- --------- -------
278 :
279 25611 : j = 1
280 25611 : ParCnt = 0
281 25611 : lFunc = LEN_TRIM(Func)
282 : step: DO
283 339781 : IF (j > lFunc) CALL ParseErrMsg(j, FuncStr)
284 339781 : c = Func(j:j)
285 : !-- -------- --------- --------- --------- --------- --------- --------- -------
286 : ! Check for valid operand (must appear)
287 : !-- -------- --------- --------- --------- --------- --------- --------- -------
288 339781 : IF (c == '-' .OR. c == '+') THEN ! Check for leading - or +
289 969 : j = j + 1
290 969 : IF (j > lFunc) CALL ParseErrMsg(j, FuncStr, 'Missing operand')
291 969 : c = Func(j:j)
292 5814 : IF (ANY(c == Ops)) CALL ParseErrMsg(j, FuncStr, 'Multiple operators')
293 : END IF
294 339781 : n = MathFunctionIndex(Func(j:))
295 339781 : IF (n > 0) THEN ! Check for math function
296 803 : j = j + LEN_TRIM(Funcs(n))
297 803 : IF (j > lFunc) CALL ParseErrMsg(j, FuncStr, 'Missing function argument')
298 803 : c = Func(j:j)
299 803 : IF (c /= '(') CALL ParseErrMsg(j, FuncStr, 'Missing opening parenthesis')
300 : END IF
301 339781 : IF (c == '(') THEN ! Check for opening parenthesis
302 110337 : ParCnt = ParCnt + 1
303 110337 : j = j + 1
304 110337 : CYCLE step
305 : END IF
306 229444 : IF (SCAN(c, '0123456789.') > 0) THEN ! Check for number
307 70616 : r = RealNum(Func(j:), ib, in, err)
308 70616 : IF (err) CALL ParseErrMsg(j, FuncStr, 'Invalid number format: '//Func(j + ib - 1:j + in - 2))
309 70616 : j = j + in - 1
310 70616 : IF (j > lFunc) EXIT
311 70099 : c = Func(j:j)
312 : ELSE ! Check for variable
313 158828 : n = VariableIndex(Func(j:), Var, ib, in)
314 158828 : IF (n == 0) CALL ParseErrMsg(j, FuncStr, 'Invalid element: '//Func(j + ib - 1:j + in - 2))
315 158828 : j = j + in - 1
316 158828 : IF (j > lFunc) EXIT
317 155646 : c = Func(j:j)
318 : END IF
319 314170 : DO WHILE (c == ')') ! Check for closing parenthesis
320 110337 : ParCnt = ParCnt - 1
321 110337 : IF (ParCnt < 0) CALL ParseErrMsg(j, FuncStr, 'Mismatched parenthesis')
322 110337 : IF (Func(j - 1:j - 1) == '(') CALL ParseErrMsg(j - 1, FuncStr, 'Empty parentheses')
323 110337 : j = j + 1
324 110337 : IF (j > lFunc) EXIT
325 314170 : c = Func(j:j)
326 : END DO
327 : !-- -------- --------- --------- --------- --------- --------- --------- -------
328 : ! Now, we have a legal operand: A legal operator or end of string must follow
329 : !-- -------- --------- --------- --------- --------- --------- --------- -------
330 225745 : IF (j > lFunc) EXIT
331 630873 : IF (ANY(c == Ops)) THEN ! Check for multiple operators
332 203833 : IF (j + 1 > lFunc) CALL ParseErrMsg(j, FuncStr)
333 1222998 : IF (ANY(Func(j + 1:j + 1) == Ops)) CALL ParseErrMsg(j + 1, FuncStr, 'Multiple operators')
334 : ELSE ! Check for next operand
335 0 : CALL ParseErrMsg(j, FuncStr, 'Missing operator')
336 : END IF
337 : !-- -------- --------- --------- --------- --------- --------- --------- -------
338 : ! Now, we have an operand and an operator: the next loop will check for another
339 : ! operand (must appear)
340 : !-- -------- --------- --------- --------- --------- --------- --------- -------
341 228927 : j = j + 1
342 : END DO step
343 25611 : IF (ParCnt > 0) CALL ParseErrMsg(j, FuncStr, 'Missing )')
344 25611 : END SUBROUTINE CheckSyntax
345 : !
346 : ! **************************************************************************************************
347 : !> \brief ...
348 : !> \return ...
349 : ! **************************************************************************************************
350 0 : FUNCTION EvalErrMsg() RESULT(msg)
351 : !----- -------- --------- --------- --------- --------- --------- --------- -------
352 : ! Return error message
353 : !----- -------- --------- --------- --------- --------- --------- --------- -------
354 : CHARACTER(LEN=*), DIMENSION(4), PARAMETER :: m = (/'Division by zero ', &
355 : 'Argument of SQRT negative ', 'Argument of LOG negative ', &
356 : 'Argument of ASIN or ACOS illegal'/)
357 : CHARACTER(LEN=LEN(m)) :: msg
358 :
359 : !----- -------- --------- --------- --------- --------- --------- --------- -------
360 :
361 0 : IF (EvalErrType < 1 .OR. EvalErrType > SIZE(m)) THEN
362 0 : msg = ''
363 : ELSE
364 0 : msg = m(EvalErrType)
365 : END IF
366 0 : END FUNCTION EvalErrMsg
367 : !
368 : ! **************************************************************************************************
369 : !> \brief ...
370 : !> \param j ...
371 : !> \param FuncStr ...
372 : !> \param Msg ...
373 : ! **************************************************************************************************
374 0 : SUBROUTINE ParseErrMsg(j, FuncStr, Msg)
375 : !----- -------- --------- --------- --------- --------- --------- --------- -------
376 : ! Print error message and terminate program
377 : !----- -------- --------- --------- --------- --------- --------- --------- -------
378 : INTEGER, INTENT(in) :: j
379 : CHARACTER(LEN=*), INTENT(in) :: FuncStr
380 : CHARACTER(LEN=*), INTENT(in), OPTIONAL :: Msg
381 :
382 : CHARACTER(LEN=default_string_length) :: message
383 :
384 : ! Original function string
385 : !----- -------- --------- --------- --------- --------- --------- --------- -------
386 :
387 0 : IF (PRESENT(Msg)) THEN
388 0 : WRITE (UNIT=message, FMT="(A)") "Syntax error in function string: "//Msg
389 : ELSE
390 0 : WRITE (UNIT=message, FMT="(A)") "Syntax error in function string"
391 : END IF
392 0 : WRITE (*, '(/,T2,A)') TRIM(FuncStr)
393 0 : IF ((j > LBOUND(ipos, DIM=1)) .AND. (j <= UBOUND(ipos, DIM=1))) THEN
394 0 : WRITE (*, '(A)') REPEAT(" ", ipos(j))//"?"
395 : ELSE
396 0 : WRITE (*, '(A)') REPEAT(" ", SIZE(ipos) + 1)//"?"
397 : END IF
398 0 : CPABORT(TRIM(message))
399 :
400 0 : END SUBROUTINE ParseErrMsg
401 : !
402 : ! **************************************************************************************************
403 : !> \brief ...
404 : !> \param c ...
405 : !> \return ...
406 : ! **************************************************************************************************
407 407666 : FUNCTION OperatorIndex(c) RESULT(n)
408 : !----- -------- --------- --------- --------- --------- --------- --------- -------
409 : ! Return operator index
410 : !----- -------- --------- --------- --------- --------- --------- --------- -------
411 : CHARACTER(LEN=1), INTENT(in) :: c
412 : INTEGER(is) :: n
413 :
414 : INTEGER(is) :: j
415 :
416 : !----- -------- --------- --------- --------- --------- --------- --------- -------
417 :
418 407666 : n = 0
419 1261746 : DO j = cAdd, cPow
420 1261746 : IF (c == Ops(j)) THEN
421 : n = j
422 : EXIT
423 : END IF
424 : END DO
425 407666 : END FUNCTION OperatorIndex
426 : !
427 : ! **************************************************************************************************
428 : !> \brief ...
429 : !> \param str ...
430 : !> \return ...
431 : ! **************************************************************************************************
432 840419 : FUNCTION MathFunctionIndex(str) RESULT(n)
433 : !----- -------- --------- --------- --------- --------- --------- --------- -------
434 : ! Return index of math function beginning at 1st position of string str
435 : !----- -------- --------- --------- --------- --------- --------- --------- -------
436 : CHARACTER(LEN=*), INTENT(in) :: str
437 : INTEGER(is) :: n
438 :
439 : CHARACTER(LEN=LEN(Funcs)) :: fun
440 : INTEGER :: k
441 : INTEGER(is) :: j
442 :
443 : !----- -------- --------- --------- --------- --------- --------- --------- -------
444 :
445 840419 : n = 0
446 14256154 : DO j = cAbs, cErfc ! Check all math functions
447 13418836 : k = MIN(LEN_TRIM(Funcs(j)), LEN(str))
448 13418836 : CALL LowCase(str(1:k), fun)
449 14256154 : IF (fun == Funcs(j)) THEN ! Compare lower case letters
450 : n = j ! Found a matching function
451 : EXIT
452 : END IF
453 : END DO
454 840419 : END FUNCTION MathFunctionIndex
455 : !
456 : ! **************************************************************************************************
457 : !> \brief ...
458 : !> \param str ...
459 : !> \param Var ...
460 : !> \param ibegin ...
461 : !> \param inext ...
462 : !> \return ...
463 : ! **************************************************************************************************
464 476484 : FUNCTION VariableIndex(str, Var, ibegin, inext) RESULT(n)
465 : !----- -------- --------- --------- --------- --------- --------- --------- -------
466 : ! Return index of variable at begin of string str (returns 0 if no variable found)
467 : !----- -------- --------- --------- --------- --------- --------- --------- -------
468 : CHARACTER(LEN=*), INTENT(in) :: str
469 : CHARACTER(LEN=*), DIMENSION(:), INTENT(in) :: Var
470 : INTEGER, INTENT(out), OPTIONAL :: ibegin, inext
471 : INTEGER(is) :: n
472 :
473 : INTEGER :: ib, in, j, lstr
474 :
475 : ! String
476 : ! Array with variable names
477 : ! Index of variable
478 : ! Start position of variable name
479 : ! Position of character after name
480 : !----- -------- --------- --------- --------- --------- --------- --------- -------
481 :
482 476484 : n = 0
483 476484 : lstr = LEN_TRIM(str)
484 476484 : IF (lstr > 0) THEN
485 476484 : DO ib = 1, lstr ! Search for first character in str
486 476484 : IF (str(ib:ib) /= ' ') EXIT ! When lstr>0 at least 1 char in str
487 : END DO
488 2381112 : DO in = ib, lstr ! Search for name terminators
489 2381112 : IF (SCAN(str(in:in), '+-*/^) ') > 0) EXIT
490 : END DO
491 1091487 : DO j = 1, SIZE(Var)
492 1091487 : IF (str(ib:in - 1) == Var(j)) THEN
493 476484 : n = INT(j, KIND=is) ! Variable name found
494 476484 : EXIT
495 : END IF
496 : END DO
497 : END IF
498 476484 : IF (PRESENT(ibegin)) ibegin = ib
499 476484 : IF (PRESENT(inext)) inext = in
500 476484 : END FUNCTION VariableIndex
501 : !
502 : ! **************************************************************************************************
503 : !> \brief ...
504 : !> \param str ...
505 : ! **************************************************************************************************
506 25611 : SUBROUTINE RemoveSpaces(str)
507 : !----- -------- --------- --------- --------- --------- --------- --------- -------
508 : ! Remove Spaces from string, remember positions of characters in old string
509 : !----- -------- --------- --------- --------- --------- --------- --------- -------
510 : CHARACTER(LEN=*), INTENT(inout) :: str
511 :
512 : INTEGER :: k, lstr
513 :
514 : !----- -------- --------- --------- --------- --------- --------- --------- -------
515 :
516 25611 : lstr = LEN_TRIM(str)
517 2432707 : ipos(:) = (/(k, k=1, lstr)/)
518 25611 : k = 1
519 1229159 : DO WHILE (str(k:lstr) /= ' ')
520 1203548 : IF (str(k:k) == ' ') THEN
521 43492 : str(k:lstr) = str(k + 1:lstr)//' ' ! Move 1 character to left
522 1267300 : ipos(k:lstr) = (/ipos(k + 1:lstr), 0/) ! Move 1 element to left
523 43492 : k = k - 1
524 : END IF
525 1203548 : k = k + 1
526 : END DO
527 25611 : END SUBROUTINE RemoveSpaces
528 : !
529 : ! **************************************************************************************************
530 : !> \brief ...
531 : !> \param ca ...
532 : !> \param cb ...
533 : !> \param str ...
534 : ! **************************************************************************************************
535 25611 : SUBROUTINE Replace(ca, cb, str)
536 : !----- -------- --------- --------- --------- --------- --------- --------- -------
537 : ! Replace ALL appearances of character set ca in string str by character set cb
538 : !----- -------- --------- --------- --------- --------- --------- --------- -------
539 : CHARACTER(LEN=*), INTENT(in) :: ca
540 : CHARACTER(LEN=LEN(ca)), INTENT(in) :: cb
541 : CHARACTER(LEN=*), INTENT(inout) :: str
542 :
543 : INTEGER :: j, lca
544 :
545 : ! LEN(ca) must be LEN(cb)
546 : !----- -------- --------- --------- --------- --------- --------- --------- -------
547 :
548 25611 : lca = LEN(ca)
549 1203548 : DO j = 1, LEN_TRIM(str) - lca + 1
550 1203548 : IF (str(j:j + lca - 1) == ca) str(j:j + lca - 1) = cb
551 : END DO
552 25611 : END SUBROUTINE Replace
553 : !
554 : ! **************************************************************************************************
555 : !> \brief ...
556 : !> \param i ...
557 : !> \param F ...
558 : !> \param Var ...
559 : ! **************************************************************************************************
560 25611 : SUBROUTINE Compile(i, F, Var)
561 : !----- -------- --------- --------- --------- --------- --------- --------- -------
562 : ! Compile i-th function string F into bytecode
563 : !----- -------- --------- --------- --------- --------- --------- --------- -------
564 : INTEGER, INTENT(in) :: i
565 : CHARACTER(LEN=*), INTENT(in) :: F
566 : CHARACTER(LEN=*), DIMENSION(:), INTENT(in) :: Var
567 :
568 : ! Function identifier
569 : ! Function string
570 : ! Array with variable names
571 : !----- -------- --------- --------- --------- --------- --------- --------- -------
572 :
573 25611 : IF (ASSOCIATED(Comp(i)%ByteCode)) DEALLOCATE (Comp(i)%ByteCode, &
574 0 : Comp(i)%Immed, &
575 0 : Comp(i)%Stack)
576 25611 : Comp(i)%ByteCodeSize = 0
577 25611 : Comp(i)%ImmedSize = 0
578 25611 : Comp(i)%StackSize = 0
579 25611 : Comp(i)%StackPtr = 0
580 25611 : CALL CompileSubstr(i, F, 1, LEN_TRIM(F), Var) ! Compile string to determine size
581 : ALLOCATE (Comp(i)%ByteCode(Comp(i)%ByteCodeSize), &
582 : Comp(i)%Immed(Comp(i)%ImmedSize), &
583 177740 : Comp(i)%Stack(Comp(i)%StackSize))
584 25611 : Comp(i)%ByteCodeSize = 0
585 25611 : Comp(i)%ImmedSize = 0
586 25611 : Comp(i)%StackSize = 0
587 25611 : Comp(i)%StackPtr = 0
588 25611 : CALL CompileSubstr(i, F, 1, LEN_TRIM(F), Var) ! Compile string into bytecode
589 : !
590 25611 : END SUBROUTINE Compile
591 : !
592 : ! **************************************************************************************************
593 : !> \brief ...
594 : !> \param i ...
595 : !> \param b ...
596 : ! **************************************************************************************************
597 870098 : SUBROUTINE AddCompiledByte(i, b)
598 : !----- -------- --------- --------- --------- --------- --------- --------- -------
599 : ! Add compiled byte to bytecode
600 : !----- -------- --------- --------- --------- --------- --------- --------- -------
601 : INTEGER, INTENT(in) :: i
602 : INTEGER(is), INTENT(in) :: b
603 :
604 : ! Function identifier
605 : ! Value of byte to be added
606 : !----- -------- --------- --------- --------- --------- --------- --------- -------
607 :
608 870098 : Comp(i)%ByteCodeSize = Comp(i)%ByteCodeSize + 1
609 870098 : IF (ASSOCIATED(Comp(i)%ByteCode)) Comp(i)%ByteCode(Comp(i)%ByteCodeSize) = b
610 870098 : END SUBROUTINE AddCompiledByte
611 : !
612 : ! **************************************************************************************************
613 : !> \brief ...
614 : !> \param i ...
615 : !> \param F ...
616 : !> \param Var ...
617 : !> \return ...
618 : ! **************************************************************************************************
619 458888 : FUNCTION MathItemIndex(i, F, Var) RESULT(n)
620 : !----- -------- --------- --------- --------- --------- --------- --------- -------
621 : ! Return math item index, if item is real number, enter it into Comp-structure
622 : !----- -------- --------- --------- --------- --------- --------- --------- -------
623 : INTEGER, INTENT(in) :: i
624 : CHARACTER(LEN=*), INTENT(in) :: F
625 : CHARACTER(LEN=*), DIMENSION(:), INTENT(in) :: Var
626 : INTEGER(is) :: n
627 :
628 : ! Function identifier
629 : ! Function substring
630 : ! Array with variable names
631 : ! Byte value of math item
632 : !----- -------- --------- --------- --------- --------- --------- --------- -------
633 :
634 458888 : n = 0
635 458888 : IF (SCAN(F(1:1), '0123456789.') > 0) THEN ! Check for begin of a number
636 141232 : Comp(i)%ImmedSize = Comp(i)%ImmedSize + 1
637 141232 : IF (ASSOCIATED(Comp(i)%Immed)) Comp(i)%Immed(Comp(i)%ImmedSize) = RealNum(F)
638 : n = cImmed
639 : ELSE ! Check for a variable
640 317656 : n = VariableIndex(F, Var)
641 317656 : IF (n > 0) n = VarBegin + n - 1_is
642 : END IF
643 458888 : END FUNCTION MathItemIndex
644 : !
645 : ! **************************************************************************************************
646 : !> \brief ...
647 : !> \param F ...
648 : !> \param b ...
649 : !> \param e ...
650 : !> \return ...
651 : ! **************************************************************************************************
652 1092166 : FUNCTION CompletelyEnclosed(F, b, e) RESULT(res)
653 : !----- -------- --------- --------- --------- --------- --------- --------- -------
654 : ! Check if function substring F(b:e) is completely enclosed by a pair of parenthesis
655 : !----- -------- --------- --------- --------- --------- --------- --------- -------
656 : CHARACTER(LEN=*), INTENT(in) :: F
657 : INTEGER, INTENT(in) :: b, e
658 : LOGICAL :: res
659 :
660 : INTEGER :: j, k
661 :
662 : ! Function substring
663 : ! First and last pos. of substring
664 : !----- -------- --------- --------- --------- --------- --------- --------- -------
665 :
666 1092166 : res = .FALSE.
667 1092166 : IF (F(b:b) == '(' .AND. F(e:e) == ')') THEN
668 220936 : k = 0
669 3896792 : DO j = b + 1, e - 1
670 3676118 : IF (F(j:j) == '(') THEN
671 261764 : k = k + 1
672 3414354 : ELSEIF (F(j:j) == ')') THEN
673 262026 : k = k - 1
674 : END IF
675 3896792 : IF (k < 0) EXIT
676 : END DO
677 220936 : IF (k == 0) res = .TRUE. ! All opened parenthesis closed
678 : END IF
679 1092166 : END FUNCTION CompletelyEnclosed
680 : !
681 : ! **************************************************************************************************
682 : !> \brief ...
683 : !> \param i ...
684 : !> \param F ...
685 : !> \param b ...
686 : !> \param e ...
687 : !> \param Var ...
688 : ! **************************************************************************************************
689 1087882 : RECURSIVE SUBROUTINE CompileSubstr(i, F, b, e, Var)
690 : !----- -------- --------- --------- --------- --------- --------- --------- -------
691 : ! Compile i-th function string F into bytecode
692 : !----- -------- --------- --------- --------- --------- --------- --------- -------
693 : INTEGER, INTENT(in) :: i
694 : CHARACTER(LEN=*), INTENT(in) :: F
695 : INTEGER, INTENT(in) :: b, e
696 : CHARACTER(LEN=*), DIMENSION(:), INTENT(in) :: Var
697 :
698 : CHARACTER(LEN=*), PARAMETER :: &
699 : calpha = 'abcdefghijklmnopqrstuvwxyz'//'ABCDEFGHIJKLMNOPQRSTUVWXYZ'
700 :
701 : INTEGER :: b2, io, j, k
702 : INTEGER(is) :: n
703 :
704 : ! Function identifier
705 : ! Function substring
706 : ! Begin and end position substring
707 : ! Array with variable names
708 : !----- -------- --------- --------- --------- --------- --------- --------- -------
709 : ! Check for special cases of substring
710 : !----- -------- --------- --------- --------- --------- --------- --------- -------
711 :
712 1087882 : IF (F(b:b) == '+') THEN ! Case 1: F(b:e) = '+...'
713 : ! WRITE(*,*)'1. F(b:e) = "+..."'
714 0 : CALL CompileSubstr(i, F, b + 1, e, Var)
715 628994 : RETURN
716 1087882 : ELSEIF (CompletelyEnclosed(F, b, e)) THEN ! Case 2: F(b:e) = '(...)'
717 : ! WRITE(*,*)'2. F(b:e) = "(...)"'
718 218948 : CALL CompileSubstr(i, F, b + 1, e - 1, Var)
719 218948 : RETURN
720 868934 : ELSEIF (SCAN(F(b:b), calpha) > 0) THEN
721 499132 : n = MathFunctionIndex(F(b:e))
722 499132 : IF (n > 0) THEN
723 2298 : b2 = b + INDEX(F(b:e), '(') - 1
724 2298 : IF (CompletelyEnclosed(F, b2, e)) THEN ! Case 3: F(b:e) = 'fcn(...)'
725 : ! WRITE(*,*)'3. F(b:e) = "fcn(...)"'
726 1606 : CALL CompileSubstr(i, F, b2 + 1, e - 1, Var)
727 1606 : CALL AddCompiledByte(i, n)
728 1606 : RETURN
729 : END IF
730 : END IF
731 369802 : ELSEIF (F(b:b) == '-') THEN
732 1986 : IF (CompletelyEnclosed(F, b + 1, e)) THEN ! Case 4: F(b:e) = '-(...)'
733 : ! WRITE(*,*)'4. F(b:e) = "-(...)"'
734 120 : CALL CompileSubstr(i, F, b + 2, e - 1, Var)
735 120 : CALL AddCompiledByte(i, cNeg)
736 120 : RETURN
737 1866 : ELSEIF (SCAN(F(b + 1:b + 1), calpha) > 0) THEN
738 1506 : n = MathFunctionIndex(F(b + 1:e))
739 1506 : IF (n > 0) THEN
740 0 : b2 = b + INDEX(F(b + 1:e), '(')
741 0 : IF (CompletelyEnclosed(F, b2, e)) THEN ! Case 5: F(b:e) = '-fcn(...)'
742 : ! WRITE(*,*)'5. F(b:e) = "-fcn(...)"'
743 0 : CALL CompileSubstr(i, F, b2 + 1, e - 1, Var)
744 0 : CALL AddCompiledByte(i, n)
745 0 : CALL AddCompiledByte(i, cNeg)
746 0 : RETURN
747 : END IF
748 : END IF
749 : END IF
750 : END IF
751 : !----- -------- --------- --------- --------- --------- --------- --------- -------
752 : ! Check for operator in substring: check only base level (k=0), exclude expr. in ()
753 : !----- -------- --------- --------- --------- --------- --------- --------- -------
754 4017360 : DO io = cAdd, cPow ! Increasing priority +-*/^
755 3558472 : k = 0
756 33589706 : DO j = e, b, -1
757 29980666 : IF (F(j:j) == ')') THEN
758 1964886 : k = k + 1
759 28015780 : ELSEIF (F(j:j) == '(') THEN
760 1964886 : k = k - 1
761 : END IF
762 33130818 : IF (k == 0 .AND. F(j:j) == Ops(io) .AND. IsBinaryOp(j, F)) THEN
763 1090520 : IF (ANY(F(j:j) == Ops(cMul:cPow)) .AND. F(b:b) == '-') THEN ! Case 6: F(b:e) = '-...Op...' with Op > -
764 : ! WRITE(*,*)'6. F(b:e) = "-...Op..." with Op > -'
765 654 : CALL CompileSubstr(i, F, b + 1, e, Var)
766 654 : CALL AddCompiledByte(i, cNeg)
767 654 : RETURN
768 : ELSE ! Case 7: F(b:e) = '...BinOp...'
769 : ! WRITE(*,*)'7. Binary operator',F(j:j)
770 407666 : CALL CompileSubstr(i, F, b, j - 1, Var)
771 407666 : CALL CompileSubstr(i, F, j + 1, e, Var)
772 407666 : CALL AddCompiledByte(i, OperatorIndex(Ops(io)))
773 407666 : Comp(i)%StackPtr = Comp(i)%StackPtr - 1
774 407666 : RETURN
775 : END IF
776 : END IF
777 : END DO
778 : END DO
779 : !----- -------- --------- --------- --------- --------- --------- --------- -------
780 : ! Check for remaining items, i.e. variables or explicit numbers
781 : !----- -------- --------- --------- --------- --------- --------- --------- -------
782 458888 : b2 = b
783 458888 : IF (F(b:b) == '-') b2 = b2 + 1
784 458888 : n = MathItemIndex(i, F(b2:e), Var)
785 : ! WRITE(*,*)'8. AddCompiledByte ',n
786 458888 : CALL AddCompiledByte(i, n)
787 458888 : Comp(i)%StackPtr = Comp(i)%StackPtr + 1
788 458888 : IF (Comp(i)%StackPtr > Comp(i)%StackSize) Comp(i)%StackSize = Comp(i)%StackSize + 1
789 458888 : IF (b2 > b) CALL AddCompiledByte(i, cNeg)
790 : END SUBROUTINE CompileSubstr
791 : !
792 : ! **************************************************************************************************
793 : !> \brief ...
794 : !> \param j ...
795 : !> \param F ...
796 : !> \return ...
797 : ! **************************************************************************************************
798 410138 : FUNCTION IsBinaryOp(j, F) RESULT(res)
799 : !----- -------- --------- --------- --------- --------- --------- --------- -------
800 : ! Check if operator F(j:j) in string F is binary operator
801 : ! Special cases already covered elsewhere: (that is corrected in v1.1)
802 : ! - operator character F(j:j) is first character of string (j=1)
803 : !----- -------- --------- --------- --------- --------- --------- --------- -------
804 : INTEGER, INTENT(in) :: j
805 : CHARACTER(LEN=*), INTENT(in) :: F
806 : LOGICAL :: res
807 :
808 : INTEGER :: k
809 : LOGICAL :: Dflag, Pflag
810 :
811 : ! Position of Operator
812 : ! String
813 : ! Result
814 : !----- -------- --------- --------- --------- --------- --------- --------- -------
815 :
816 410138 : res = .TRUE.
817 410138 : IF (F(j:j) == '+' .OR. F(j:j) == '-') THEN ! Plus or minus sign:
818 139840 : IF (j == 1) THEN ! - leading unary operator ?
819 : res = .FALSE.
820 138656 : ELSEIF (SCAN(F(j - 1:j - 1), '+-*/^(') > 0) THEN ! - other unary operator ?
821 : res = .FALSE.
822 138022 : ELSEIF (SCAN(F(j + 1:j + 1), '0123456789') > 0 .AND. & ! - in exponent of real number ?
823 : SCAN(F(j - 1:j - 1), 'eEdD') > 0) THEN
824 : Dflag = .FALSE.; Pflag = .FALSE.
825 : k = j - 1
826 0 : DO WHILE (k > 1) ! step to the left in mantissa
827 0 : k = k - 1
828 0 : IF (SCAN(F(k:k), '0123456789') > 0) THEN
829 : Dflag = .TRUE.
830 0 : ELSEIF (F(k:k) == '.') THEN
831 0 : IF (Pflag) THEN
832 : EXIT ! * EXIT: 2nd appearance of '.'
833 : ELSE
834 : Pflag = .TRUE. ! * mark 1st appearance of '.'
835 : END IF
836 : ELSE
837 : EXIT ! * all other characters
838 : END IF
839 : END DO
840 0 : IF (Dflag .AND. (k == 1 .OR. SCAN(F(k:k), '+-*/^(') > 0)) res = .FALSE.
841 : END IF
842 : END IF
843 410138 : END FUNCTION IsBinaryOp
844 : !
845 : ! **************************************************************************************************
846 : !> \brief ...
847 : !> \param str ...
848 : !> \param ibegin ...
849 : !> \param inext ...
850 : !> \param error ...
851 : !> \return ...
852 : ! **************************************************************************************************
853 141232 : FUNCTION RealNum(str, ibegin, inext, error) RESULT(res)
854 : !----- -------- --------- --------- --------- --------- --------- --------- -------
855 : ! Get real number from string - Format: [blanks][+|-][nnn][.nnn][e|E|d|D[+|-]nnn]
856 : !----- -------- --------- --------- --------- --------- --------- --------- -------
857 : CHARACTER(LEN=*), INTENT(in) :: str
858 : INTEGER, INTENT(out), OPTIONAL :: ibegin, inext
859 : LOGICAL, INTENT(out), OPTIONAL :: error
860 : REAL(rn) :: res
861 :
862 : INTEGER :: ib, in, istat
863 : LOGICAL :: Bflag, DInExp, DInMan, Eflag, err, &
864 : InExp, InMan, Pflag
865 :
866 : ! String
867 : ! Real number
868 : ! Start position of real number
869 : ! 1st character after real number
870 : ! Error flag
871 : ! .T. at begin of number in str
872 : ! .T. in mantissa of number
873 : ! .T. after 1st '.' encountered
874 : ! .T. at exponent identifier 'eEdD'
875 : ! .T. in exponent of number
876 : ! .T. if at least 1 digit in mant.
877 : ! .T. if at least 1 digit in exp.
878 : ! Local error flag
879 : !----- -------- --------- --------- --------- --------- --------- --------- -------
880 :
881 141232 : Bflag = .TRUE.; InMan = .FALSE.; Pflag = .FALSE.; Eflag = .FALSE.; InExp = .FALSE.
882 141232 : DInMan = .FALSE.; DInExp = .FALSE.
883 141232 : ib = 1
884 141232 : in = 1
885 335706 : DO WHILE (in <= LEN_TRIM(str))
886 264573 : SELECT CASE (str(in:in))
887 : CASE (' ') ! Only leading blanks permitted
888 0 : ib = ib + 1
889 0 : IF (InMan .OR. Eflag .OR. InExp) EXIT
890 : CASE ('+', '-') ! Permitted only
891 23746 : IF (Bflag) THEN
892 : InMan = .TRUE.; Bflag = .FALSE. ! - at beginning of mantissa
893 23746 : ELSEIF (Eflag) THEN
894 : InExp = .TRUE.; Eflag = .FALSE. ! - at beginning of exponent
895 : ELSE
896 : EXIT ! - otherwise STOP
897 : END IF
898 : CASE ('0':'9') ! Mark
899 190596 : IF (Bflag) THEN
900 : InMan = .TRUE.; Bflag = .FALSE. ! - beginning of mantissa
901 49364 : ELSEIF (Eflag) THEN
902 0 : InExp = .TRUE.; Eflag = .FALSE. ! - beginning of exponent
903 : END IF
904 190596 : IF (InMan) DInMan = .TRUE. ! Mantissa contains digit
905 190596 : IF (InExp) DInExp = .TRUE. ! Exponent contains digit
906 : CASE ('.')
907 3878 : IF (Bflag) THEN
908 : Pflag = .TRUE. ! - mark 1st appearance of '.'
909 : InMan = .TRUE.; Bflag = .FALSE. ! mark beginning of mantissa
910 3878 : ELSEIF (InMan .AND. .NOT. Pflag) THEN
911 : Pflag = .TRUE. ! - mark 1st appearance of '.'
912 : ELSE
913 : EXIT ! - otherwise STOP
914 : END IF
915 : CASE ('e', 'E', 'd', 'D') ! Permitted only
916 0 : IF (InMan) THEN
917 : Eflag = .TRUE.; InMan = .FALSE. ! - following mantissa
918 : ELSE
919 : EXIT ! - otherwise STOP
920 : END IF
921 : CASE DEFAULT
922 264573 : EXIT ! STOP at all other characters
923 : END SELECT
924 218220 : in = in + 1
925 : END DO
926 141232 : err = (ib > in - 1) .OR. (.NOT. DInMan) .OR. ((Eflag .OR. InExp) .AND. .NOT. DInExp)
927 : IF (err) THEN
928 0 : res = 0.0_rn
929 : ELSE
930 141232 : READ (str(ib:in - 1), *, IOSTAT=istat) res
931 141232 : err = istat /= 0
932 : END IF
933 141232 : IF (PRESENT(ibegin)) ibegin = ib
934 141232 : IF (PRESENT(inext)) inext = in
935 141232 : IF (PRESENT(error)) error = err
936 141232 : END FUNCTION RealNum
937 : !
938 : ! **************************************************************************************************
939 : !> \brief ...
940 : !> \param str1 ...
941 : !> \param str2 ...
942 : ! **************************************************************************************************
943 13418836 : SUBROUTINE LowCase(str1, str2)
944 : !----- -------- --------- --------- --------- --------- --------- --------- -------
945 : ! Transform upper case letters in str1 into lower case letters, result is str2
946 : !----- -------- --------- --------- --------- --------- --------- --------- -------
947 : CHARACTER(LEN=*), INTENT(in) :: str1
948 : CHARACTER(LEN=*), INTENT(out) :: str2
949 :
950 : CHARACTER(LEN=*), PARAMETER :: lc = 'abcdefghijklmnopqrstuvwxyz', &
951 : uc = 'ABCDEFGHIJKLMNOPQRSTUVWXYZ'
952 :
953 : INTEGER :: j, k
954 :
955 : !----- -------- --------- --------- --------- --------- --------- --------- -------
956 :
957 13418836 : str2 = str1
958 57057287 : DO j = 1, LEN_TRIM(str1)
959 43638451 : k = INDEX(uc, str1(j:j))
960 57057287 : IF (k > 0) str2(j:j) = lc(k:k)
961 : END DO
962 13418836 : END SUBROUTINE LowCase
963 :
964 : ! **************************************************************************************************
965 : !> \brief Evaluates derivatives
966 : !> \param id_fun ...
967 : !> \param ipar ...
968 : !> \param vals ...
969 : !> \param h ...
970 : !> \param err ...
971 : !> \return ...
972 : !> \author Main algorithm from Numerical Recipes
973 : !> Ridders, C.J.F. 1982 - Advances in Engineering Software, Vol.4, no. 2, pp. 75-76
974 : ! **************************************************************************************************
975 2983 : FUNCTION evalfd(id_fun, ipar, vals, h, err) RESULT(derivative)
976 : INTEGER, INTENT(IN) :: id_fun, ipar
977 : REAL(KIND=rn), DIMENSION(:), INTENT(INOUT) :: vals
978 : REAL(KIND=rn), INTENT(IN) :: h
979 : REAL(KIND=rn), INTENT(OUT) :: err
980 : REAL(KIND=rn) :: derivative
981 :
982 : INTEGER, PARAMETER :: ntab = 10
983 : REAL(KIND=rn), PARAMETER :: big_error = 1.0E30_rn, con = 1.4_rn, &
984 : con2 = con*con, safe = 2.0_rn
985 :
986 : INTEGER :: i, j
987 : REAL(KIND=rn) :: a(ntab, ntab), errt, fac, funcm, funcp, &
988 : hh, xval
989 :
990 2983 : derivative = HUGE(0.0_rn)
991 2983 : IF (h /= 0._rn) THEN
992 2983 : xval = vals(ipar)
993 2983 : hh = h
994 2983 : vals(ipar) = xval + hh
995 2983 : funcp = evalf(id_fun, vals)
996 2983 : vals(ipar) = xval - hh
997 2983 : funcm = evalf(id_fun, vals)
998 2983 : a(1, 1) = (funcp - funcm)/(2.0_rn*hh)
999 2983 : err = big_error
1000 7175 : DO i = 2, ntab
1001 7175 : hh = hh/con
1002 7175 : vals(ipar) = xval + hh
1003 7175 : funcp = evalf(id_fun, vals)
1004 7175 : vals(ipar) = xval - hh
1005 7175 : funcm = evalf(id_fun, vals)
1006 7175 : a(1, i) = (funcp - funcm)/(2.0_rn*hh)
1007 7175 : fac = con2
1008 22816 : DO j = 2, i
1009 15641 : a(j, i) = (a(j - 1, i)*fac - a(j - 1, i - 1))/(fac - 1.0_rn)
1010 15641 : fac = con2*fac
1011 15641 : errt = MAX(ABS(a(j, i) - a(j - 1, i)), ABS(a(j, i) - a(j - 1, i - 1)))
1012 22816 : IF (errt .LE. err) THEN
1013 6997 : err = errt
1014 6997 : derivative = a(j, i)
1015 : END IF
1016 : END DO
1017 7175 : IF (ABS(a(i, i) - a(i - 1, i - 1)) .GE. safe*err) RETURN
1018 : END DO
1019 : ELSE
1020 0 : CPABORT("DX provided equals zero!")
1021 : END IF
1022 0 : vals(ipar) = xval
1023 0 : END FUNCTION evalfd
1024 :
1025 0 : END MODULE fparser
|
