(git:374b731)
Loading...
Searching...
No Matches
fparser.F
Go to the documentation of this file.
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! **************************************************************************************************
17MODULE 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 !
95CONTAINS
96 !
97! **************************************************************************************************
98!> \brief ...
99! **************************************************************************************************
100 SUBROUTINE finalizef()
101 !----- -------- --------- --------- --------- --------- --------- --------- -------
102 ! Finalize function parser
103 !----- -------- --------- --------- --------- --------- --------- --------- -------
104
105 INTEGER :: i
106
107!----- -------- --------- --------- --------- --------- --------- --------- -------
108
109 DO i = 1, SIZE(comp)
110 IF (ASSOCIATED(comp(i)%ByteCode)) THEN
111 DEALLOCATE (comp(i)%ByteCode)
112 END IF
113 IF (ASSOCIATED(comp(i)%Immed)) THEN
114 DEALLOCATE (comp(i)%Immed)
115 END IF
116 IF (ASSOCIATED(comp(i)%Stack)) THEN
117 DEALLOCATE (comp(i)%Stack)
118 END IF
119 END DO
120
121 DEALLOCATE (comp)
122
123 END SUBROUTINE finalizef
124 !
125! **************************************************************************************************
126!> \brief ...
127!> \param n ...
128! **************************************************************************************************
129 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 ALLOCATE (comp(n))
139 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 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 CHARACTER(LEN=LEN(FuncStr)) :: func
153
154! Function string, local use
155!----- -------- --------- --------- --------- --------- --------- --------- -------
156
157 IF (i < 1 .OR. i > SIZE(comp)) THEN
158 cpabort("Function number is out of range")
159 END IF
160 IF (SIZE(var) > huge(0_is)) THEN
161 cpabort("Too many variables")
162 END IF
163 ALLOCATE (ipos(len_trim(funcstr))) ! Char. positions in orig. string
164 func = funcstr ! Local copy of function string
165 CALL replace('**', '^ ', func) ! Exponent into 1-Char. format
166 CALL removespaces(func) ! Condense function string
167 CALL checksyntax(func, funcstr, var)
168 DEALLOCATE (ipos)
169 CALL compile(i, func, var) ! Compile into bytecode
170
171 END SUBROUTINE parsef
172 !
173! **************************************************************************************************
174!> \brief ...
175!> \param i ...
176!> \param Val ...
177!> \return ...
178! **************************************************************************************************
179 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 dp = 1
200 sp = 0
201 DO ip = 1, comp(i)%ByteCodeSize
202 SELECT CASE (comp(i)%ByteCode(ip))
203
204 CASE (cimmed); sp = sp + 1; comp(i)%Stack(sp) = comp(i)%Immed(dp); dp = dp + 1
205 CASE (cneg); comp(i)%Stack(sp) = -comp(i)%Stack(sp)
206 CASE (cadd); comp(i)%Stack(sp - 1) = comp(i)%Stack(sp - 1) + comp(i)%Stack(sp); sp = sp - 1
207 CASE (csub); comp(i)%Stack(sp - 1) = comp(i)%Stack(sp - 1) - comp(i)%Stack(sp); sp = sp - 1
208 CASE (cmul); comp(i)%Stack(sp - 1) = comp(i)%Stack(sp - 1)*comp(i)%Stack(sp); sp = sp - 1
209 CASE (cdiv); IF (comp(i)%Stack(sp) == 0._rn) then; evalerrtype = 1; res = zero; return; END IF
210 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 IF (comp(i)%Stack(sp - 1) < 0.0_rn) THEN
214 ipow = floor(comp(i)%Stack(sp))
215 IF (mod(comp(i)%Stack(sp), real(ipow, kind=rn)) == 0.0_rn) THEN
216 comp(i)%Stack(sp - 1) = comp(i)%Stack(sp - 1)**ipow
217 ELSE
218 cpabort("Negative floating-point value raised to a real power!")
219 END IF
220 ELSE
221 comp(i)%Stack(sp - 1) = comp(i)%Stack(sp - 1)**comp(i)%Stack(sp)
222 END IF
223 sp = sp - 1
224 CASE (cabs); comp(i)%Stack(sp) = abs(comp(i)%Stack(sp))
225 CASE (cexp); comp(i)%Stack(sp) = exp(comp(i)%Stack(sp))
226 CASE (clog10); IF (comp(i)%Stack(sp) <= 0._rn) then; evalerrtype = 3; res = zero; return; END IF
227 comp(i)%Stack(sp) = log10(comp(i)%Stack(sp))
228 CASE (clog); IF (comp(i)%Stack(sp) <= 0._rn) then; evalerrtype = 3; res = zero; return; END IF
229 comp(i)%Stack(sp) = log(comp(i)%Stack(sp))
230 CASE (csqrt); IF (comp(i)%Stack(sp) < 0._rn) then; evalerrtype = 3; res = zero; return; END IF
231 comp(i)%Stack(sp) = sqrt(comp(i)%Stack(sp))
232 CASE (csinh); comp(i)%Stack(sp) = sinh(comp(i)%Stack(sp))
233 CASE (ccosh); comp(i)%Stack(sp) = cosh(comp(i)%Stack(sp))
234 CASE (ctanh); comp(i)%Stack(sp) = tanh(comp(i)%Stack(sp))
235 CASE (csin); comp(i)%Stack(sp) = sin(comp(i)%Stack(sp))
236 CASE (ccos); comp(i)%Stack(sp) = cos(comp(i)%Stack(sp))
237 CASE (ctan); comp(i)%Stack(sp) = tan(comp(i)%Stack(sp))
238 CASE (casin); IF ((comp(i)%Stack(sp) < -1._rn) .OR. (comp(i)%Stack(sp) > 1._rn)) THEN
239 evalerrtype = 4; res = zero; return; END IF
240 comp(i)%Stack(sp) = asin(comp(i)%Stack(sp))
241 CASE (cacos); IF ((comp(i)%Stack(sp) < -1._rn) .OR. (comp(i)%Stack(sp) > 1._rn)) THEN
242 evalerrtype = 4; res = zero; return; END IF
243 comp(i)%Stack(sp) = acos(comp(i)%Stack(sp))
244 CASE (catan); comp(i)%Stack(sp) = atan(comp(i)%Stack(sp))
245 CASE (cerf); comp(i)%Stack(sp) = erf(comp(i)%Stack(sp))
246 CASE (cerfc); comp(i)%Stack(sp) = erfc(comp(i)%Stack(sp))
247 CASE DEFAULT; sp = sp + 1; comp(i)%Stack(sp) = val(comp(i)%ByteCode(ip) - varbegin + 1)
248 END SELECT
249 END DO
250 evalerrtype = 0
251 res = comp(i)%Stack(1)
252 END FUNCTION evalf
253 !
254! **************************************************************************************************
255!> \brief ...
256!> \param Func ...
257!> \param FuncStr ...
258!> \param Var ...
259! **************************************************************************************************
260 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 j = 1
280 parcnt = 0
281 lfunc = len_trim(func)
282 step: DO
283 IF (j > lfunc) CALL parseerrmsg(j, funcstr)
284 c = func(j:j)
285 !-- -------- --------- --------- --------- --------- --------- --------- -------
286 ! Check for valid operand (must appear)
287 !-- -------- --------- --------- --------- --------- --------- --------- -------
288 IF (c == '-' .OR. c == '+') THEN ! Check for leading - or +
289 j = j + 1
290 IF (j > lfunc) CALL parseerrmsg(j, funcstr, 'Missing operand')
291 c = func(j:j)
292 IF (any(c == ops)) CALL parseerrmsg(j, funcstr, 'Multiple operators')
293 END IF
294 n = mathfunctionindex(func(j:))
295 IF (n > 0) THEN ! Check for math function
296 j = j + len_trim(funcs(n))
297 IF (j > lfunc) CALL parseerrmsg(j, funcstr, 'Missing function argument')
298 c = func(j:j)
299 IF (c /= '(') CALL parseerrmsg(j, funcstr, 'Missing opening parenthesis')
300 END IF
301 IF (c == '(') THEN ! Check for opening parenthesis
302 parcnt = parcnt + 1
303 j = j + 1
304 cycle step
305 END IF
306 IF (scan(c, '0123456789.') > 0) THEN ! Check for number
307 r = realnum(func(j:), ib, in, err)
308 IF (err) CALL parseerrmsg(j, funcstr, 'Invalid number format: '//func(j + ib - 1:j + in - 2))
309 j = j + in - 1
310 IF (j > lfunc) EXIT
311 c = func(j:j)
312 ELSE ! Check for variable
313 n = variableindex(func(j:), var, ib, in)
314 IF (n == 0) CALL parseerrmsg(j, funcstr, 'Invalid element: '//func(j + ib - 1:j + in - 2))
315 j = j + in - 1
316 IF (j > lfunc) EXIT
317 c = func(j:j)
318 END IF
319 DO WHILE (c == ')') ! Check for closing parenthesis
320 parcnt = parcnt - 1
321 IF (parcnt < 0) CALL parseerrmsg(j, funcstr, 'Mismatched parenthesis')
322 IF (func(j - 1:j - 1) == '(') CALL parseerrmsg(j - 1, funcstr, 'Empty parentheses')
323 j = j + 1
324 IF (j > lfunc) EXIT
325 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 IF (j > lfunc) EXIT
331 IF (any(c == ops)) THEN ! Check for multiple operators
332 IF (j + 1 > lfunc) CALL parseerrmsg(j, funcstr)
333 IF (any(func(j + 1:j + 1) == ops)) CALL parseerrmsg(j + 1, funcstr, 'Multiple operators')
334 ELSE ! Check for next operand
335 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 j = j + 1
342 END DO step
343 IF (parcnt > 0) CALL parseerrmsg(j, funcstr, 'Missing )')
344 END SUBROUTINE checksyntax
345 !
346! **************************************************************************************************
347!> \brief ...
348!> \return ...
349! **************************************************************************************************
350 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 IF (evalerrtype < 1 .OR. evalerrtype > SIZE(m)) THEN
362 msg = ''
363 ELSE
364 msg = m(evalerrtype)
365 END IF
366 END FUNCTION evalerrmsg
367 !
368! **************************************************************************************************
369!> \brief ...
370!> \param j ...
371!> \param FuncStr ...
372!> \param Msg ...
373! **************************************************************************************************
374 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 IF (PRESENT(msg)) THEN
388 WRITE (unit=message, fmt="(A)") "Syntax error in function string: "//msg
389 ELSE
390 WRITE (unit=message, fmt="(A)") "Syntax error in function string"
391 END IF
392 WRITE (*, '(/,T2,A)') trim(funcstr)
393 IF ((j > lbound(ipos, dim=1)) .AND. (j <= ubound(ipos, dim=1))) THEN
394 WRITE (*, '(A)') repeat(" ", ipos(j))//"?"
395 ELSE
396 WRITE (*, '(A)') repeat(" ", SIZE(ipos) + 1)//"?"
397 END IF
398 cpabort(trim(message))
399
400 END SUBROUTINE parseerrmsg
401 !
402! **************************************************************************************************
403!> \brief ...
404!> \param c ...
405!> \return ...
406! **************************************************************************************************
407 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 n = 0
419 DO j = cadd, cpow
420 IF (c == ops(j)) THEN
421 n = j
422 EXIT
423 END IF
424 END DO
425 END FUNCTION operatorindex
426 !
427! **************************************************************************************************
428!> \brief ...
429!> \param str ...
430!> \return ...
431! **************************************************************************************************
432 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 n = 0
446 DO j = cabs, cerfc ! Check all math functions
447 k = min(len_trim(funcs(j)), len(str))
448 CALL lowcase(str(1:k), fun)
449 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 END FUNCTION mathfunctionindex
455 !
456! **************************************************************************************************
457!> \brief ...
458!> \param str ...
459!> \param Var ...
460!> \param ibegin ...
461!> \param inext ...
462!> \return ...
463! **************************************************************************************************
464 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 n = 0
483 lstr = len_trim(str)
484 IF (lstr > 0) THEN
485 DO ib = 1, lstr ! Search for first character in str
486 IF (str(ib:ib) /= ' ') EXIT ! When lstr>0 at least 1 char in str
487 END DO
488 DO in = ib, lstr ! Search for name terminators
489 IF (scan(str(in:in), '+-*/^) ') > 0) EXIT
490 END DO
491 DO j = 1, SIZE(var)
492 IF (str(ib:in - 1) == var(j)) THEN
493 n = int(j, kind=is) ! Variable name found
494 EXIT
495 END IF
496 END DO
497 END IF
498 IF (PRESENT(ibegin)) ibegin = ib
499 IF (PRESENT(inext)) inext = in
500 END FUNCTION variableindex
501 !
502! **************************************************************************************************
503!> \brief ...
504!> \param str ...
505! **************************************************************************************************
506 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 lstr = len_trim(str)
517 ipos(:) = (/(k, k=1, lstr)/)
518 k = 1
519 DO WHILE (str(k:lstr) /= ' ')
520 IF (str(k:k) == ' ') THEN
521 str(k:lstr) = str(k + 1:lstr)//' ' ! Move 1 character to left
522 ipos(k:lstr) = (/ipos(k + 1:lstr), 0/) ! Move 1 element to left
523 k = k - 1
524 END IF
525 k = k + 1
526 END DO
527 END SUBROUTINE removespaces
528 !
529! **************************************************************************************************
530!> \brief ...
531!> \param ca ...
532!> \param cb ...
533!> \param str ...
534! **************************************************************************************************
535 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 lca = len(ca)
549 DO j = 1, len_trim(str) - lca + 1
550 IF (str(j:j + lca - 1) == ca) str(j:j + lca - 1) = cb
551 END DO
552 END SUBROUTINE replace
553 !
554! **************************************************************************************************
555!> \brief ...
556!> \param i ...
557!> \param F ...
558!> \param Var ...
559! **************************************************************************************************
560 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 IF (ASSOCIATED(comp(i)%ByteCode)) DEALLOCATE (comp(i)%ByteCode, &
574 comp(i)%Immed, &
575 comp(i)%Stack)
576 comp(i)%ByteCodeSize = 0
577 comp(i)%ImmedSize = 0
578 comp(i)%StackSize = 0
579 comp(i)%StackPtr = 0
580 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 comp(i)%Stack(comp(i)%StackSize))
584 comp(i)%ByteCodeSize = 0
585 comp(i)%ImmedSize = 0
586 comp(i)%StackSize = 0
587 comp(i)%StackPtr = 0
588 CALL compilesubstr(i, f, 1, len_trim(f), var) ! Compile string into bytecode
589 !
590 END SUBROUTINE compile
591 !
592! **************************************************************************************************
593!> \brief ...
594!> \param i ...
595!> \param b ...
596! **************************************************************************************************
597 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 comp(i)%ByteCodeSize = comp(i)%ByteCodeSize + 1
609 IF (ASSOCIATED(comp(i)%ByteCode)) comp(i)%ByteCode(comp(i)%ByteCodeSize) = b
610 END SUBROUTINE addcompiledbyte
611 !
612! **************************************************************************************************
613!> \brief ...
614!> \param i ...
615!> \param F ...
616!> \param Var ...
617!> \return ...
618! **************************************************************************************************
619 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 n = 0
635 IF (scan(f(1:1), '0123456789.') > 0) THEN ! Check for begin of a number
636 comp(i)%ImmedSize = comp(i)%ImmedSize + 1
637 IF (ASSOCIATED(comp(i)%Immed)) comp(i)%Immed(comp(i)%ImmedSize) = realnum(f)
638 n = cimmed
639 ELSE ! Check for a variable
640 n = variableindex(f, var)
641 IF (n > 0) n = varbegin + n - 1_is
642 END IF
643 END FUNCTION mathitemindex
644 !
645! **************************************************************************************************
646!> \brief ...
647!> \param F ...
648!> \param b ...
649!> \param e ...
650!> \return ...
651! **************************************************************************************************
652 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 res = .false.
667 IF (f(b:b) == '(' .AND. f(e:e) == ')') THEN
668 k = 0
669 DO j = b + 1, e - 1
670 IF (f(j:j) == '(') THEN
671 k = k + 1
672 ELSEIF (f(j:j) == ')') THEN
673 k = k - 1
674 END IF
675 IF (k < 0) EXIT
676 END DO
677 IF (k == 0) res = .true. ! All opened parenthesis closed
678 END IF
679 END FUNCTION completelyenclosed
680 !
681! **************************************************************************************************
682!> \brief ...
683!> \param i ...
684!> \param F ...
685!> \param b ...
686!> \param e ...
687!> \param Var ...
688! **************************************************************************************************
689 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 IF (f(b:b) == '+') THEN ! Case 1: F(b:e) = '+...'
713! WRITE(*,*)'1. F(b:e) = "+..."'
714 CALL compilesubstr(i, f, b + 1, e, var)
715 RETURN
716 ELSEIF (completelyenclosed(f, b, e)) THEN ! Case 2: F(b:e) = '(...)'
717! WRITE(*,*)'2. F(b:e) = "(...)"'
718 CALL compilesubstr(i, f, b + 1, e - 1, var)
719 RETURN
720 ELSEIF (scan(f(b:b), calpha) > 0) THEN
721 n = mathfunctionindex(f(b:e))
722 IF (n > 0) THEN
723 b2 = b + index(f(b:e), '(') - 1
724 IF (completelyenclosed(f, b2, e)) THEN ! Case 3: F(b:e) = 'fcn(...)'
725! WRITE(*,*)'3. F(b:e) = "fcn(...)"'
726 CALL compilesubstr(i, f, b2 + 1, e - 1, var)
727 CALL addcompiledbyte(i, n)
728 RETURN
729 END IF
730 END IF
731 ELSEIF (f(b:b) == '-') THEN
732 IF (completelyenclosed(f, b + 1, e)) THEN ! Case 4: F(b:e) = '-(...)'
733! WRITE(*,*)'4. F(b:e) = "-(...)"'
734 CALL compilesubstr(i, f, b + 2, e - 1, var)
735 CALL addcompiledbyte(i, cneg)
736 RETURN
737 ELSEIF (scan(f(b + 1:b + 1), calpha) > 0) THEN
738 n = mathfunctionindex(f(b + 1:e))
739 IF (n > 0) THEN
740 b2 = b + index(f(b + 1:e), '(')
741 IF (completelyenclosed(f, b2, e)) THEN ! Case 5: F(b:e) = '-fcn(...)'
742! WRITE(*,*)'5. F(b:e) = "-fcn(...)"'
743 CALL compilesubstr(i, f, b2 + 1, e - 1, var)
744 CALL addcompiledbyte(i, n)
745 CALL addcompiledbyte(i, cneg)
746 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 DO io = cadd, cpow ! Increasing priority +-*/^
755 k = 0
756 DO j = e, b, -1
757 IF (f(j:j) == ')') THEN
758 k = k + 1
759 ELSEIF (f(j:j) == '(') THEN
760 k = k - 1
761 END IF
762 IF (k == 0 .AND. f(j:j) == ops(io) .AND. isbinaryop(j, f)) THEN
763 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 CALL compilesubstr(i, f, b + 1, e, var)
766 CALL addcompiledbyte(i, cneg)
767 RETURN
768 ELSE ! Case 7: F(b:e) = '...BinOp...'
769! WRITE(*,*)'7. Binary operator',F(j:j)
770 CALL compilesubstr(i, f, b, j - 1, var)
771 CALL compilesubstr(i, f, j + 1, e, var)
772 CALL addcompiledbyte(i, operatorindex(ops(io)))
773 comp(i)%StackPtr = comp(i)%StackPtr - 1
774 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 b2 = b
783 IF (f(b:b) == '-') b2 = b2 + 1
784 n = mathitemindex(i, f(b2:e), var)
785! WRITE(*,*)'8. AddCompiledByte ',n
786 CALL addcompiledbyte(i, n)
787 comp(i)%StackPtr = comp(i)%StackPtr + 1
788 IF (comp(i)%StackPtr > comp(i)%StackSize) comp(i)%StackSize = comp(i)%StackSize + 1
789 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 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 res = .true.
817 IF (f(j:j) == '+' .OR. f(j:j) == '-') THEN ! Plus or minus sign:
818 IF (j == 1) THEN ! - leading unary operator ?
819 res = .false.
820 ELSEIF (scan(f(j - 1:j - 1), '+-*/^(') > 0) THEN ! - other unary operator ?
821 res = .false.
822 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 DO WHILE (k > 1) ! step to the left in mantissa
827 k = k - 1
828 IF (scan(f(k:k), '0123456789') > 0) THEN
829 dflag = .true.
830 ELSEIF (f(k:k) == '.') THEN
831 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 IF (dflag .AND. (k == 1 .OR. scan(f(k:k), '+-*/^(') > 0)) res = .false.
841 END IF
842 END IF
843 END FUNCTION isbinaryop
844 !
845! **************************************************************************************************
846!> \brief ...
847!> \param str ...
848!> \param ibegin ...
849!> \param inext ...
850!> \param error ...
851!> \return ...
852! **************************************************************************************************
853 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 bflag = .true.; inman = .false.; pflag = .false.; eflag = .false.; inexp = .false.
882 dinman = .false.; dinexp = .false.
883 ib = 1
884 in = 1
885 DO WHILE (in <= len_trim(str))
886 SELECT CASE (str(in:in))
887 CASE (' ') ! Only leading blanks permitted
888 ib = ib + 1
889 IF (inman .OR. eflag .OR. inexp) EXIT
890 CASE ('+', '-') ! Permitted only
891 IF (bflag) THEN
892 inman = .true.; bflag = .false. ! - at beginning of mantissa
893 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 IF (bflag) THEN
900 inman = .true.; bflag = .false. ! - beginning of mantissa
901 ELSEIF (eflag) THEN
902 inexp = .true.; eflag = .false. ! - beginning of exponent
903 END IF
904 IF (inman) dinman = .true. ! Mantissa contains digit
905 IF (inexp) dinexp = .true. ! Exponent contains digit
906 CASE ('.')
907 IF (bflag) THEN
908 pflag = .true. ! - mark 1st appearance of '.'
909 inman = .true.; bflag = .false. ! mark beginning of mantissa
910 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 IF (inman) THEN
917 eflag = .true.; inman = .false. ! - following mantissa
918 ELSE
919 EXIT ! - otherwise STOP
920 END IF
921 CASE DEFAULT
922 EXIT ! STOP at all other characters
923 END SELECT
924 in = in + 1
925 END DO
926 err = (ib > in - 1) .OR. (.NOT. dinman) .OR. ((eflag .OR. inexp) .AND. .NOT. dinexp)
927 IF (err) THEN
928 res = 0.0_rn
929 ELSE
930 READ (str(ib:in - 1), *, iostat=istat) res
931 err = istat /= 0
932 END IF
933 IF (PRESENT(ibegin)) ibegin = ib
934 IF (PRESENT(inext)) inext = in
935 IF (PRESENT(error)) error = err
936 END FUNCTION realnum
937 !
938! **************************************************************************************************
939!> \brief ...
940!> \param str1 ...
941!> \param str2 ...
942! **************************************************************************************************
943 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 str2 = str1
958 DO j = 1, len_trim(str1)
959 k = index(uc, str1(j:j))
960 IF (k > 0) str2(j:j) = lc(k:k)
961 END DO
962 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 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 derivative = huge(0.0_rn)
991 IF (h /= 0._rn) THEN
992 xval = vals(ipar)
993 hh = h
994 vals(ipar) = xval + hh
995 funcp = evalf(id_fun, vals)
996 vals(ipar) = xval - hh
997 funcm = evalf(id_fun, vals)
998 a(1, 1) = (funcp - funcm)/(2.0_rn*hh)
999 err = big_error
1000 DO i = 2, ntab
1001 hh = hh/con
1002 vals(ipar) = xval + hh
1003 funcp = evalf(id_fun, vals)
1004 vals(ipar) = xval - hh
1005 funcm = evalf(id_fun, vals)
1006 a(1, i) = (funcp - funcm)/(2.0_rn*hh)
1007 fac = con2
1008 DO j = 2, i
1009 a(j, i) = (a(j - 1, i)*fac - a(j - 1, i - 1))/(fac - 1.0_rn)
1010 fac = con2*fac
1011 errt = max(abs(a(j, i) - a(j - 1, i)), abs(a(j, i) - a(j - 1, i - 1)))
1012 IF (errt .LE. err) THEN
1013 err = errt
1014 derivative = a(j, i)
1015 END IF
1016 END DO
1017 IF (abs(a(i, i) - a(i - 1, i - 1)) .GE. safe*err) RETURN
1018 END DO
1019 ELSE
1020 cpabort("DX provided equals zero!")
1021 END IF
1022 vals(ipar) = xval
1023 END FUNCTION evalfd
1024
1025END MODULE fparser
fac
static GRID_HOST_DEVICE double fac(const int i)
Factorial function, e.g. fac(5) = 5! = 120.
Definition grid_common.h:48
fparser
This public domain function parser module is intended for applications where a set of mathematical ex...
Definition fparser.F:17
fparser::parsef
subroutine, public parsef(i, funcstr, var)
Parse ith function string FuncStr and compile it into bytecode.
Definition fparser.F:148
fparser::evalf
real(rn) function, public evalf(i, val)
...
Definition fparser.F:180
fparser::evalerrtype
integer, public evalerrtype
Definition fparser.F:32
fparser::evalfd
real(kind=rn) function, public evalfd(id_fun, ipar, vals, h, err)
Evaluates derivatives.
Definition fparser.F:976
fparser::comp
type(tcomp), dimension(:), pointer comp
Definition fparser.F:92
fparser::finalizef
subroutine, public finalizef()
...
Definition fparser.F:101
fparser::initf
subroutine, public initf(n)
...
Definition fparser.F:130
kinds
Defines the basic variable types.
Definition kinds.F:23
kinds::dp
integer, parameter, public dp
Definition kinds.F:34
kinds::default_string_length
integer, parameter, public default_string_length
Definition kinds.F:57
kinds::sp
integer, parameter, public sp
Definition kinds.F:33