(git:4180bf5)
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-2026 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 docf
33 INTEGER, PUBLIC :: evalerrtype ! =0: no error occurred, >0: evaluation error
34 !------- -------- --------- --------- --------- --------- --------- --------- -------
35 PRIVATE
36 SAVE
37 INTEGER, PARAMETER, PRIVATE :: is = selected_int_kind(1) !Data type of bytecode
38 INTEGER(is), PARAMETER :: cimmed = 1, &
39 cneg = 2, &
40 cadd = 3, &
41 csub = 4, &
42 cmul = 5, &
43 cdiv = 6, &
44 cpow = 7, &
45 cabs = 8, &
46 cmin = 9, &
47 cmax = 10, &
48 cexp = 11, &
49 clog10 = 12, &
50 clog = 13, &
51 csqrt = 14, &
52 csinh = 15, &
53 ccosh = 16, &
54 ctanh = 17, &
55 csin = 18, &
56 ccos = 19, &
57 ctan = 20, &
58 casin = 21, &
59 cacos = 22, &
60 catan2 = 23, &
61 catan = 24, &
62 cerf = 25, &
63 cerfc = 26, &
64 varbegin = 27
65 CHARACTER(LEN=1), DIMENSION(cAdd:cPow), PARAMETER :: ops = ['+', &
66 '-', &
67 '*', &
68 '/', &
69 '^']
70 CHARACTER(LEN=5), DIMENSION(cAbs:cErfc), PARAMETER :: funcs = ['abs ', &
71 'min ', &
72 'max ', &
73 'exp ', &
74 'log10', &
75 'log ', &
76 'sqrt ', &
77 'sinh ', &
78 'cosh ', &
79 'tanh ', &
80 'sin ', &
81 'cos ', &
82 'tan ', &
83 'asin ', &
84 'acos ', &
85 'atan2', &
86 'atan ', &
87 'erf ', &
88 'erfc ']
89 INTEGER, DIMENSION(cAbs:cErfc), PARAMETER :: funcsargcnt = [1, & ! abs
90 2, & ! min
91 2, & ! max
92 1, & ! exp
93 1, & ! log10
94 1, & ! log
95 1, & ! sqrt
96 1, & ! sinh
97 1, & ! cosh
98 1, & ! tanh
99 1, & ! sin
100 1, & ! cos
101 1, & ! tan
102 1, & ! asin
103 1, & ! acos
104 2, & ! atan2
105 1, & ! atan
106 1, & ! erf
107 1] ! erfc
108! **************************************************************************************************
109 TYPE tcomp
110 INTEGER(is), DIMENSION(:), POINTER :: bytecode => null()
111 INTEGER :: bytecodesize = -1
112 REAL(rn), DIMENSION(:), POINTER :: immed => null()
113 INTEGER :: immedsize = -1
114 REAL(rn), DIMENSION(:), POINTER :: stack => null()
115 INTEGER :: stacksize = -1, &
116 stackptr = -1
117 END TYPE tcomp
118 TYPE(tcomp), DIMENSION(:), POINTER :: comp => null() ! Bytecode
119 INTEGER, DIMENSION(:), ALLOCATABLE :: ipos ! Associates function strings
120 !
121CONTAINS
122 !
123! **************************************************************************************************
124!> \brief ...
125! **************************************************************************************************
126 SUBROUTINE finalizef()
127 !----- -------- --------- --------- --------- --------- --------- --------- -------
128 ! Finalize function parser
129 !----- -------- --------- --------- --------- --------- --------- --------- -------
130
131 INTEGER :: i
132
133!----- -------- --------- --------- --------- --------- --------- --------- -------
134
135 DO i = 1, SIZE(comp)
136 IF (ASSOCIATED(comp(i)%ByteCode)) THEN
137 DEALLOCATE (comp(i)%ByteCode)
138 END IF
139 IF (ASSOCIATED(comp(i)%Immed)) THEN
140 DEALLOCATE (comp(i)%Immed)
141 END IF
142 IF (ASSOCIATED(comp(i)%Stack)) THEN
143 DEALLOCATE (comp(i)%Stack)
144 END IF
145 END DO
146
147 DEALLOCATE (comp)
148
149 END SUBROUTINE finalizef
150 !
151! **************************************************************************************************
152!> \brief ...
153!> \param n ...
154! **************************************************************************************************
155 SUBROUTINE initf(n)
156 !----- -------- --------- --------- --------- --------- --------- --------- -------
157 ! Initialize function parser for n functions
158 !----- -------- --------- --------- --------- --------- --------- --------- -------
159 INTEGER, INTENT(in) :: n
160
161! Number of functions
162!----- -------- --------- --------- --------- --------- --------- --------- -------
163
164 ALLOCATE (comp(n))
165 END SUBROUTINE initf
166 !
167! **************************************************************************************************
168!> \brief Parse ith function string FuncStr and compile it into bytecode
169!> \param i Function identifier
170!> \param FuncStr Function string
171!> \param Var Array with variable names
172! **************************************************************************************************
173 SUBROUTINE parsef(i, FuncStr, Var)
174 INTEGER, INTENT(in) :: i
175 CHARACTER(LEN=*), INTENT(in) :: funcstr
176 CHARACTER(LEN=*), DIMENSION(:), INTENT(in) :: var
177
178 CHARACTER(LEN=LEN(FuncStr)) :: func
179
180! Function string, local use
181!----- -------- --------- --------- --------- --------- --------- --------- -------
182
183 IF (i < 1 .OR. i > SIZE(comp)) THEN
184 cpabort("Function number is out of range")
185 END IF
186 IF (SIZE(var) > huge(0_is)) THEN
187 cpabort("Too many variables")
188 END IF
189 ALLOCATE (ipos(len_trim(funcstr))) ! Char. positions in orig. string
190 func = funcstr ! Local copy of function string
191 CALL replace('**', '^ ', func) ! Exponent into 1-Char. format
192 CALL removespaces(func) ! Condense function string
193 CALL checksyntax(func, funcstr, var)
194 DEALLOCATE (ipos)
195 CALL compile(i, func, var) ! Compile into bytecode
196
197 END SUBROUTINE parsef
198 !
199! **************************************************************************************************
200!> \brief ...
201!> \param i ...
202!> \param Val ...
203!> \return ...
204! **************************************************************************************************
205 FUNCTION evalf(i, Val) RESULT(res)
206 !----- -------- --------- --------- --------- --------- --------- --------- -------
207 ! Evaluate bytecode of ith function for the values passed in array Val(:)
208 !----- -------- --------- --------- --------- --------- --------- --------- -------
209 INTEGER, INTENT(in) :: i
210 REAL(rn), DIMENSION(:), INTENT(in) :: val
211 REAL(rn) :: res
212
213 REAL(rn), PARAMETER :: zero = 0._rn
214
215 INTEGER :: dp, ip, ipow, sp
216
217! Function identifier
218! Variable values
219! Result
220! Instruction pointer
221! Data pointer
222! Stack pointer
223!----- -------- --------- --------- --------- --------- --------- --------- -------
224
225 dp = 1
226 sp = 0
227 DO ip = 1, comp(i)%ByteCodeSize
228 SELECT CASE (comp(i)%ByteCode(ip))
229
230 CASE (cimmed); sp = sp + 1; comp(i)%Stack(sp) = comp(i)%Immed(dp); dp = dp + 1
231 CASE (cneg); comp(i)%Stack(sp) = -comp(i)%Stack(sp)
232 CASE (cadd); comp(i)%Stack(sp - 1) = comp(i)%Stack(sp - 1) + comp(i)%Stack(sp); sp = sp - 1
233 CASE (csub); comp(i)%Stack(sp - 1) = comp(i)%Stack(sp - 1) - comp(i)%Stack(sp); sp = sp - 1
234 CASE (cmul); comp(i)%Stack(sp - 1) = comp(i)%Stack(sp - 1)*comp(i)%Stack(sp); sp = sp - 1
235 CASE (cdiv); IF (comp(i)%Stack(sp) == 0._rn) then; evalerrtype = 1; res = zero; return; END IF
236 comp(i)%Stack(sp - 1) = comp(i)%Stack(sp - 1)/comp(i)%Stack(sp); sp = sp - 1
237 CASE (cpow)
238 ! Fixing for possible Negative floating-point value raised to a real power
239 IF (comp(i)%Stack(sp - 1) < 0.0_rn) THEN
240 ipow = floor(comp(i)%Stack(sp))
241 IF (mod(comp(i)%Stack(sp), real(ipow, kind=rn)) == 0.0_rn) THEN
242 comp(i)%Stack(sp - 1) = comp(i)%Stack(sp - 1)**ipow
243 ELSE
244 cpabort("Negative floating-point value raised to a real power!")
245 END IF
246 ELSE
247 comp(i)%Stack(sp - 1) = comp(i)%Stack(sp - 1)**comp(i)%Stack(sp)
248 END IF
249 sp = sp - 1
250 CASE (cabs); comp(i)%Stack(sp) = abs(comp(i)%Stack(sp))
251 CASE (cmin); comp(i)%Stack(sp - 1) = min(comp(i)%Stack(sp - 1), comp(i)%Stack(sp)); sp = sp - 1
252 CASE (cmax); comp(i)%Stack(sp - 1) = max(comp(i)%Stack(sp - 1), comp(i)%Stack(sp)); sp = sp - 1
253 CASE (cexp); comp(i)%Stack(sp) = exp(comp(i)%Stack(sp))
254 CASE (clog10); IF (comp(i)%Stack(sp) <= 0._rn) then; evalerrtype = 3; res = zero; return; END IF
255 comp(i)%Stack(sp) = log10(comp(i)%Stack(sp))
256 CASE (clog); IF (comp(i)%Stack(sp) <= 0._rn) then; evalerrtype = 3; res = zero; return; END IF
257 comp(i)%Stack(sp) = log(comp(i)%Stack(sp))
258 CASE (csqrt); IF (comp(i)%Stack(sp) < 0._rn) then; evalerrtype = 3; res = zero; return; END IF
259 comp(i)%Stack(sp) = sqrt(comp(i)%Stack(sp))
260 CASE (csinh); comp(i)%Stack(sp) = sinh(comp(i)%Stack(sp))
261 CASE (ccosh); comp(i)%Stack(sp) = cosh(comp(i)%Stack(sp))
262 CASE (ctanh); comp(i)%Stack(sp) = tanh(comp(i)%Stack(sp))
263 CASE (csin); comp(i)%Stack(sp) = sin(comp(i)%Stack(sp))
264 CASE (ccos); comp(i)%Stack(sp) = cos(comp(i)%Stack(sp))
265 CASE (ctan); comp(i)%Stack(sp) = tan(comp(i)%Stack(sp))
266 CASE (casin); IF ((comp(i)%Stack(sp) < -1._rn) .OR. (comp(i)%Stack(sp) > 1._rn)) THEN
267 evalerrtype = 4; res = zero; return; END IF
268 comp(i)%Stack(sp) = asin(comp(i)%Stack(sp))
269 CASE (cacos); IF ((comp(i)%Stack(sp) < -1._rn) .OR. (comp(i)%Stack(sp) > 1._rn)) THEN
270 evalerrtype = 4; res = zero; return; END IF
271 comp(i)%Stack(sp) = acos(comp(i)%Stack(sp))
272 CASE (catan2); comp(i)%Stack(sp - 1) = atan2(comp(i)%Stack(sp - 1), comp(i)%Stack(sp)); sp = sp - 1
273 CASE (catan); comp(i)%Stack(sp) = atan(comp(i)%Stack(sp))
274 CASE (cerf); comp(i)%Stack(sp) = erf(comp(i)%Stack(sp))
275 CASE (cerfc); comp(i)%Stack(sp) = erfc(comp(i)%Stack(sp))
276 CASE DEFAULT; sp = sp + 1; comp(i)%Stack(sp) = val(comp(i)%ByteCode(ip) - varbegin + 1)
277 END SELECT
278 END DO
279 evalerrtype = 0
280 res = comp(i)%Stack(1)
281 END FUNCTION evalf
282 !
283! **************************************************************************************************
284!> \brief ...
285!> \param Func ...
286!> \param FuncStr ...
287!> \param Var ...
288! **************************************************************************************************
289 SUBROUTINE checksyntax(Func, FuncStr, Var)
290 !----- -------- --------- --------- --------- --------- --------- --------- -------
291 ! Check syntax of function string, returns 0 if syntax is ok
292 !----- -------- --------- --------- --------- --------- --------- --------- -------
293 CHARACTER(LEN=*), INTENT(in) :: func, funcstr
294 CHARACTER(LEN=*), DIMENSION(:), INTENT(in) :: var
295
296 INTEGER :: ib, in, j, lfunc, parcnt
297 CHARACTER(LEN=1) :: c
298 INTEGER(is) :: n
299 LOGICAL :: err
300 REAL(rn) :: r
301
302! Function string without spaces
303! Original function string
304! Array with variable names
305! Parenthesis counter
306!----- -------- --------- --------- --------- --------- --------- --------- -------
307
308 j = 1
309 parcnt = 0
310 lfunc = len_trim(func)
311 step: DO
312 IF (j > lfunc) CALL parseerrmsg(j, funcstr)
313 c = func(j:j)
314 !-- -------- --------- --------- --------- --------- --------- --------- -------
315 ! Check for valid operand (must appear)
316 !-- -------- --------- --------- --------- --------- --------- --------- -------
317 IF (c == '-' .OR. c == '+') THEN ! Check for leading - or +
318 j = j + 1
319 IF (j > lfunc) CALL parseerrmsg(j, funcstr, 'Missing operand')
320 c = func(j:j)
321 IF (any(c == ops)) CALL parseerrmsg(j, funcstr, 'Multiple operators')
322 END IF
323 n = mathfunctionindex(func(j:))
324 IF (n > 0) THEN ! Check for math function
325 j = j + len_trim(funcs(n))
326 IF (j > lfunc) CALL parseerrmsg(j, funcstr, 'Missing function argument')
327 c = func(j:j)
328 IF (c /= '(') CALL parseerrmsg(j, funcstr, 'Missing opening parenthesis')
329 CALL checkfuncargcnt(func, funcstr, j, lfunc, n)
330 END IF
331 IF (c == '(') THEN ! Check for opening parenthesis
332 parcnt = parcnt + 1
333 j = j + 1
334 cycle step
335 END IF
336 IF (scan(c, '0123456789.') > 0) THEN ! Check for number
337 r = realnum(func(j:), ib, in, err)
338 IF (err) CALL parseerrmsg(j, funcstr, 'Invalid number format: '//func(j + ib - 1:j + in - 2))
339 j = j + in - 1
340 IF (j > lfunc) EXIT step
341 c = func(j:j)
342 ELSE ! Check for variable
343 n = variableindex(func(j:), var, ib, in)
344 IF (n == 0) CALL parseerrmsg(j, funcstr, 'Invalid element: '//func(j + ib - 1:j + in - 2))
345 j = j + in - 1
346 IF (j > lfunc) EXIT step
347 c = func(j:j)
348 END IF
349 DO WHILE (c == ')') ! Check for closing parenthesis
350 parcnt = parcnt - 1
351 IF (parcnt < 0) CALL parseerrmsg(j, funcstr, 'Mismatched parenthesis')
352 IF (func(j - 1:j - 1) == '(') CALL parseerrmsg(j - 1, funcstr, 'Empty parentheses')
353 j = j + 1
354 IF (j > lfunc) EXIT
355 c = func(j:j)
356 END DO
357 !-- -------- --------- --------- --------- --------- --------- --------- -------
358 ! Now, we have a legal operand: A legal operator or end of string must follow
359 !-- -------- --------- --------- --------- --------- --------- --------- -------
360 IF (j > lfunc) EXIT step
361 IF (any(c == ops)) THEN ! Check for multiple operators
362 IF (j + 1 > lfunc) CALL parseerrmsg(j, funcstr)
363 IF (any(func(j + 1:j + 1) == ops)) CALL parseerrmsg(j + 1, funcstr, 'Multiple operators')
364 ELSE IF (c == ',') THEN ! Check for commas separating function arguments
365 IF (parcnt == 0) CALL parseerrmsg(j, funcstr, 'Comma outside function')
366 IF (func(j + 1:j + 1) == ',') CALL parseerrmsg(j, funcstr, 'Multiple commas')
367 ELSE ! Check for next operand
368 CALL parseerrmsg(j, funcstr, 'Missing operator')
369 END IF
370 !-- -------- --------- --------- --------- --------- --------- --------- -------
371 ! Now, we have an operand and an operator: the next loop will check for another
372 ! operand (must appear)
373 !-- -------- --------- --------- --------- --------- --------- --------- -------
374 j = j + 1
375 END DO step
376 IF (parcnt > 0) CALL parseerrmsg(j, funcstr, 'Missing )')
377 END SUBROUTINE checksyntax
378 !
379! **************************************************************************************************
380!> \brief ...
381!> \param Func ...
382!> \param FuncStr ...
383!> \param b ...
384!> \param e ...
385!> \param FuncId ...
386! **************************************************************************************************
387 SUBROUTINE checkfuncargcnt(Func, FuncStr, b, e, FuncId)
388 !----- -------- --------- --------- --------- --------- --------- --------- -------
389 ! Check argument count of function substring, returns 0 if count is ok
390 !----- -------- --------- --------- --------- --------- --------- --------- -------
391 CHARACTER(LEN=*), INTENT(in) :: func, funcstr
392 INTEGER, INTENT(in) :: b, e
393 INTEGER(is), INTENT(in) :: funcid
394
395 CHARACTER(len=40) :: msg
396 INTEGER :: argcnt, argpos, j, parcnt
397
398! Function string without spaces
399! Original function string
400! Begin and end position substring
401! ID of function for which arguments are parsed
402!----- -------- --------- --------- --------- --------- --------- --------- -------
403
404 argcnt = 1
405 argpos = 0
406 parcnt = 0
407 DO j = b, e
408 IF (func(j:j) == '(') THEN
409 parcnt = parcnt + 1
410 ELSEIF (func(j:j) == ')') THEN
411 parcnt = parcnt - 1
412 IF (parcnt == 0) EXIT
413 ELSEIF (parcnt == 1 .AND. func(j:j) == ',') THEN
414 argpos = j
415 argcnt = argcnt + 1
416 END IF
417 END DO
418 IF (argcnt /= funcsargcnt(funcid)) THEN
419 IF (argcnt < funcsargcnt(funcid)) argpos = j
420 WRITE (msg, '(I0,A,A,A,I0)') argcnt, ' argument(s) in ', trim(funcs(funcid)), ' instead of ', funcsargcnt(funcid)
421 CALL parseerrmsg(argpos, funcstr, msg)
422 END IF
423 END SUBROUTINE checkfuncargcnt
424 !
425! **************************************************************************************************
426!> \brief ...
427!> \return ...
428! **************************************************************************************************
429 FUNCTION evalerrmsg() RESULT(msg)
430 !----- -------- --------- --------- --------- --------- --------- --------- -------
431 ! Return error message
432 !----- -------- --------- --------- --------- --------- --------- --------- -------
433 CHARACTER(LEN=*), DIMENSION(4), PARAMETER :: m = ['Division by zero ', &
434 'Argument of SQRT negative ', 'Argument of LOG negative ', &
435 'Argument of ASIN or ACOS illegal']
436 CHARACTER(LEN=LEN(m)) :: msg
437
438!----- -------- --------- --------- --------- --------- --------- --------- -------
439
440 IF (evalerrtype < 1 .OR. evalerrtype > SIZE(m)) THEN
441 msg = ''
442 ELSE
443 msg = m(evalerrtype)
444 END IF
445 END FUNCTION evalerrmsg
446 !
447! **************************************************************************************************
448!> \brief ...
449!> \param j ...
450!> \param FuncStr ...
451!> \param Msg ...
452! **************************************************************************************************
453 SUBROUTINE parseerrmsg(j, FuncStr, Msg)
454 !----- -------- --------- --------- --------- --------- --------- --------- -------
455 ! Print error message and terminate program
456 !----- -------- --------- --------- --------- --------- --------- --------- -------
457 INTEGER, INTENT(in) :: j
458 CHARACTER(LEN=*), INTENT(in) :: funcstr
459 CHARACTER(LEN=*), INTENT(in), OPTIONAL :: msg
460
461 CHARACTER(LEN=default_string_length) :: message
462
463! Original function string
464!----- -------- --------- --------- --------- --------- --------- --------- -------
465
466 IF (PRESENT(msg)) THEN
467 WRITE (unit=message, fmt="(A)") "Syntax error in function string: "//msg
468 ELSE
469 WRITE (unit=message, fmt="(A)") "Syntax error in function string"
470 END IF
471 WRITE (*, '(/,T2,A)') trim(funcstr)
472 IF ((j > lbound(ipos, dim=1)) .AND. (j <= ubound(ipos, dim=1))) THEN
473 WRITE (*, '(A)') repeat(" ", ipos(j))//"?"
474 ELSE
475 WRITE (*, '(A)') repeat(" ", SIZE(ipos) + 1)//"?"
476 END IF
477 cpabort(trim(message))
478
479 END SUBROUTINE parseerrmsg
480 !
481! **************************************************************************************************
482!> \brief ...
483!> \param c ...
484!> \return ...
485! **************************************************************************************************
486 FUNCTION operatorindex(c) RESULT(n)
487 !----- -------- --------- --------- --------- --------- --------- --------- -------
488 ! Return operator index
489 !----- -------- --------- --------- --------- --------- --------- --------- -------
490 CHARACTER(LEN=1), INTENT(in) :: c
491 INTEGER(is) :: n
492
493 INTEGER(is) :: j
494
495!----- -------- --------- --------- --------- --------- --------- --------- -------
496
497 n = 0
498 DO j = cadd, cpow
499 IF (c == ops(j)) THEN
500 n = j
501 EXIT
502 END IF
503 END DO
504 END FUNCTION operatorindex
505 !
506! **************************************************************************************************
507!> \brief ...
508!> \param str ...
509!> \return ...
510! **************************************************************************************************
511 FUNCTION mathfunctionindex(str) RESULT(n)
512 !----- -------- --------- --------- --------- --------- --------- --------- -------
513 ! Return index of math function beginning at 1st position of string str
514 !----- -------- --------- --------- --------- --------- --------- --------- -------
515 CHARACTER(LEN=*), INTENT(in) :: str
516 INTEGER(is) :: n
517
518 CHARACTER(LEN=LEN(Funcs)) :: fun
519 INTEGER :: k
520 INTEGER(is) :: j
521
522!----- -------- --------- --------- --------- --------- --------- --------- -------
523
524 n = 0
525 DO j = cabs, cerfc ! Check all math functions
526 k = min(len_trim(funcs(j)), len(str))
527 CALL lowcase(str(1:k), fun)
528 IF (fun == funcs(j)) THEN ! Compare lower case letters
529 n = j ! Found a matching function
530 EXIT
531 END IF
532 END DO
533 END FUNCTION mathfunctionindex
534 !
535! **************************************************************************************************
536!> \brief ...
537!> \param str ...
538!> \param Var ...
539!> \param ibegin ...
540!> \param inext ...
541!> \return ...
542! **************************************************************************************************
543 FUNCTION variableindex(str, Var, ibegin, inext) RESULT(n)
544 !----- -------- --------- --------- --------- --------- --------- --------- -------
545 ! Return index of variable at begin of string str (returns 0 if no variable found)
546 !----- -------- --------- --------- --------- --------- --------- --------- -------
547 CHARACTER(LEN=*), INTENT(in) :: str
548 CHARACTER(LEN=*), DIMENSION(:), INTENT(in) :: var
549 INTEGER, INTENT(out), OPTIONAL :: ibegin, inext
550 INTEGER(is) :: n
551
552 INTEGER :: ib, in, j, lstr
553
554! String
555! Array with variable names
556! Index of variable
557! Start position of variable name
558! Position of character after name
559!----- -------- --------- --------- --------- --------- --------- --------- -------
560
561 n = 0
562 lstr = len_trim(str)
563 IF (lstr > 0) THEN
564 DO ib = 1, lstr ! Search for first character in str
565 IF (str(ib:ib) /= ' ') EXIT ! When lstr>0 at least 1 char in str
566 END DO
567 DO in = ib, lstr ! Search for name terminators
568 IF (scan(str(in:in), '+-*/^) ,') > 0) EXIT
569 END DO
570 DO j = 1, SIZE(var)
571 IF (str(ib:in - 1) == var(j)) THEN
572 n = int(j, kind=is) ! Variable name found
573 EXIT
574 END IF
575 END DO
576 END IF
577 IF (PRESENT(ibegin)) ibegin = ib
578 IF (PRESENT(inext)) inext = in
579 END FUNCTION variableindex
580 !
581! **************************************************************************************************
582!> \brief ...
583!> \param str ...
584! **************************************************************************************************
585 SUBROUTINE removespaces(str)
586 !----- -------- --------- --------- --------- --------- --------- --------- -------
587 ! Remove Spaces from string, remember positions of characters in old string
588 !----- -------- --------- --------- --------- --------- --------- --------- -------
589 CHARACTER(LEN=*), INTENT(inout) :: str
590
591 INTEGER :: k, lstr
592
593!----- -------- --------- --------- --------- --------- --------- --------- -------
594
595 lstr = len_trim(str)
596 ipos(:) = [(k, k=1, lstr)]
597 k = 1
598 DO WHILE (str(k:lstr) /= ' ')
599 IF (str(k:k) == ' ') THEN
600 str(k:lstr) = str(k + 1:lstr)//' ' ! Move 1 character to left
601 ipos(k:lstr) = [ipos(k + 1:lstr), 0] ! Move 1 element to left
602 k = k - 1
603 END IF
604 k = k + 1
605 END DO
606 END SUBROUTINE removespaces
607 !
608! **************************************************************************************************
609!> \brief ...
610!> \param ca ...
611!> \param cb ...
612!> \param str ...
613! **************************************************************************************************
614 SUBROUTINE replace(ca, cb, str)
615 !----- -------- --------- --------- --------- --------- --------- --------- -------
616 ! Replace ALL appearances of character set ca in string str by character set cb
617 !----- -------- --------- --------- --------- --------- --------- --------- -------
618 CHARACTER(LEN=*), INTENT(in) :: ca
619 CHARACTER(LEN=LEN(ca)), INTENT(in) :: cb
620 CHARACTER(LEN=*), INTENT(inout) :: str
621
622 INTEGER :: j, lca
623
624! LEN(ca) must be LEN(cb)
625!----- -------- --------- --------- --------- --------- --------- --------- -------
626
627 lca = len(ca)
628 DO j = 1, len_trim(str) - lca + 1
629 IF (str(j:j + lca - 1) == ca) str(j:j + lca - 1) = cb
630 END DO
631 END SUBROUTINE replace
632 !
633! **************************************************************************************************
634!> \brief ...
635!> \param i ...
636!> \param F ...
637!> \param Var ...
638! **************************************************************************************************
639 SUBROUTINE compile(i, F, Var)
640 !----- -------- --------- --------- --------- --------- --------- --------- -------
641 ! Compile i-th function string F into bytecode
642 !----- -------- --------- --------- --------- --------- --------- --------- -------
643 INTEGER, INTENT(in) :: i
644 CHARACTER(LEN=*), INTENT(in) :: f
645 CHARACTER(LEN=*), DIMENSION(:), INTENT(in) :: var
646
647! Function identifier
648! Function string
649! Array with variable names
650!----- -------- --------- --------- --------- --------- --------- --------- -------
651
652 IF (ASSOCIATED(comp(i)%ByteCode)) DEALLOCATE (comp(i)%ByteCode, &
653 comp(i)%Immed, &
654 comp(i)%Stack)
655 comp(i)%ByteCodeSize = 0
656 comp(i)%ImmedSize = 0
657 comp(i)%StackSize = 0
658 comp(i)%StackPtr = 0
659 CALL compilesubstr(i, f, 1, len_trim(f), var) ! Compile string to determine size
660 ALLOCATE (comp(i)%ByteCode(comp(i)%ByteCodeSize), &
661 comp(i)%Immed(comp(i)%ImmedSize), &
662 comp(i)%Stack(comp(i)%StackSize))
663 comp(i)%ByteCodeSize = 0
664 comp(i)%ImmedSize = 0
665 comp(i)%StackSize = 0
666 comp(i)%StackPtr = 0
667 CALL compilesubstr(i, f, 1, len_trim(f), var) ! Compile string into bytecode
668 !
669 END SUBROUTINE compile
670 !
671! **************************************************************************************************
672!> \brief ...
673!> \param i ...
674!> \param b ...
675! **************************************************************************************************
676 SUBROUTINE addcompiledbyte(i, b)
677 !----- -------- --------- --------- --------- --------- --------- --------- -------
678 ! Add compiled byte to bytecode
679 !----- -------- --------- --------- --------- --------- --------- --------- -------
680 INTEGER, INTENT(in) :: i
681 INTEGER(is), INTENT(in) :: b
682
683! Function identifier
684! Value of byte to be added
685!----- -------- --------- --------- --------- --------- --------- --------- -------
686
687 comp(i)%ByteCodeSize = comp(i)%ByteCodeSize + 1
688 IF (ASSOCIATED(comp(i)%ByteCode)) comp(i)%ByteCode(comp(i)%ByteCodeSize) = b
689 END SUBROUTINE addcompiledbyte
690 !
691! **************************************************************************************************
692!> \brief ...
693!> \param i ...
694!> \param F ...
695!> \param Var ...
696!> \return ...
697! **************************************************************************************************
698 FUNCTION mathitemindex(i, F, Var) RESULT(n)
699 !----- -------- --------- --------- --------- --------- --------- --------- -------
700 ! Return math item index, if item is real number, enter it into Comp-structure
701 !----- -------- --------- --------- --------- --------- --------- --------- -------
702 INTEGER, INTENT(in) :: i
703 CHARACTER(LEN=*), INTENT(in) :: f
704 CHARACTER(LEN=*), DIMENSION(:), INTENT(in) :: var
705 INTEGER(is) :: n
706
707! Function identifier
708! Function substring
709! Array with variable names
710! Byte value of math item
711!----- -------- --------- --------- --------- --------- --------- --------- -------
712
713 n = 0
714 IF (scan(f(1:1), '0123456789.') > 0) THEN ! Check for begin of a number
715 comp(i)%ImmedSize = comp(i)%ImmedSize + 1
716 IF (ASSOCIATED(comp(i)%Immed)) comp(i)%Immed(comp(i)%ImmedSize) = realnum(f)
717 n = cimmed
718 ELSE ! Check for a variable
719 n = variableindex(f, var)
720 IF (n > 0) n = varbegin + n - 1_is
721 END IF
722 END FUNCTION mathitemindex
723 !
724! **************************************************************************************************
725!> \brief ...
726!> \param F ...
727!> \param b ...
728!> \param e ...
729!> \return ...
730! **************************************************************************************************
731 FUNCTION completelyenclosed(F, b, e) RESULT(res)
732 !----- -------- --------- --------- --------- --------- --------- --------- -------
733 ! Check if function substring F(b:e) is completely enclosed by a pair of parenthesis
734 !----- -------- --------- --------- --------- --------- --------- --------- -------
735 CHARACTER(LEN=*), INTENT(in) :: f
736 INTEGER, INTENT(in) :: b, e
737 LOGICAL :: res
738
739 INTEGER :: j, k
740
741! Function substring
742! First and last pos. of substring
743!----- -------- --------- --------- --------- --------- --------- --------- -------
744
745 res = .false.
746 IF (f(b:b) == '(' .AND. f(e:e) == ')') THEN
747 k = 0
748 DO j = b + 1, e - 1
749 IF (f(j:j) == '(') THEN
750 k = k + 1
751 ELSEIF (f(j:j) == ')') THEN
752 k = k - 1
753 END IF
754 IF (k < 0) EXIT
755 END DO
756 IF (k == 0) res = .true. ! All opened parenthesis closed
757 END IF
758 END FUNCTION completelyenclosed
759 !
760! **************************************************************************************************
761!> \brief ...
762!> \param i ...
763!> \param F ...
764!> \param b ...
765!> \param e ...
766!> \param Var ...
767! **************************************************************************************************
768 RECURSIVE SUBROUTINE compilesubstr(i, F, b, e, Var)
769 !----- -------- --------- --------- --------- --------- --------- --------- -------
770 ! Compile i-th function string F into bytecode
771 !----- -------- --------- --------- --------- --------- --------- --------- -------
772 INTEGER, INTENT(in) :: i
773 CHARACTER(LEN=*), INTENT(in) :: f
774 INTEGER, INTENT(in) :: b, e
775 CHARACTER(LEN=*), DIMENSION(:), INTENT(in) :: var
776
777 CHARACTER(LEN=*), PARAMETER :: &
778 calpha = 'abcdefghijklmnopqrstuvwxyz'//'ABCDEFGHIJKLMNOPQRSTUVWXYZ'
779
780 INTEGER :: b2, io, j, k
781 INTEGER(is) :: n
782
783! Function identifier
784! Function substring
785! Begin and end position substring
786! Array with variable names
787!----- -------- --------- --------- --------- --------- --------- --------- -------
788! Check for special cases of substring
789!----- -------- --------- --------- --------- --------- --------- --------- -------
790
791 IF (f(b:b) == '+') THEN ! Case 1: F(b:e) = '+...'
792! WRITE(*,*)'1. F(b:e) = "+..."'
793 CALL compilesubstr(i, f, b + 1, e, var)
794 RETURN
795 ELSEIF (completelyenclosed(f, b, e)) THEN ! Case 2: F(b:e) = '(...)'
796! WRITE(*,*)'2. F(b:e) = "(...)"'
797 CALL compilesubstr(i, f, b + 1, e - 1, var)
798 RETURN
799 ELSEIF (scan(f(b:b), calpha) > 0) THEN
800 n = mathfunctionindex(f(b:e))
801 IF (n > 0) THEN
802 b2 = b + index(f(b:e), '(') - 1
803 IF (completelyenclosed(f, b2, e)) THEN ! Case 3: F(b:e) = 'fcn(...)'
804! WRITE(*,*)'3. F(b:e) = "fcn(...)"'
805 CALL compilefuncargssubstr(i, f, b2 + 1, e - 1, var)
806 CALL addcompiledbyte(i, n)
807 RETURN
808 END IF
809 END IF
810 ELSEIF (f(b:b) == '-') THEN
811 IF (completelyenclosed(f, b + 1, e)) THEN ! Case 4: F(b:e) = '-(...)'
812! WRITE(*,*)'4. F(b:e) = "-(...)"'
813 CALL compilesubstr(i, f, b + 2, e - 1, var)
814 CALL addcompiledbyte(i, cneg)
815 RETURN
816 ELSEIF (scan(f(b + 1:b + 1), calpha) > 0) THEN
817 n = mathfunctionindex(f(b + 1:e))
818 IF (n > 0) THEN
819 b2 = b + index(f(b + 1:e), '(')
820 IF (completelyenclosed(f, b2, e)) THEN ! Case 5: F(b:e) = '-fcn(...)'
821! WRITE(*,*)'5. F(b:e) = "-fcn(...)"'
822 CALL compilefuncargssubstr(i, f, b2 + 1, e - 1, var)
823 CALL addcompiledbyte(i, n)
824 CALL addcompiledbyte(i, cneg)
825 RETURN
826 END IF
827 END IF
828 END IF
829 END IF
830 !----- -------- --------- --------- --------- --------- --------- --------- -------
831 ! Check for operator in substring: check only base level (k=0), exclude expr. in ()
832 !----- -------- --------- --------- --------- --------- --------- --------- -------
833 DO io = cadd, cpow ! Increasing priority +-*/^
834 k = 0
835 DO j = e, b, -1
836 IF (f(j:j) == ')') THEN
837 k = k + 1
838 ELSEIF (f(j:j) == '(') THEN
839 k = k - 1
840 END IF
841 IF (k == 0 .AND. f(j:j) == ops(io) .AND. isbinaryop(j, f)) THEN
842 IF (any(f(j:j) == ops(cmul:cpow)) .AND. f(b:b) == '-') THEN ! Case 6: F(b:e) = '-...Op...' with Op > -
843! WRITE(*,*)'6. F(b:e) = "-...Op..." with Op > -'
844 CALL compilesubstr(i, f, b + 1, e, var)
845 CALL addcompiledbyte(i, cneg)
846 RETURN
847 ELSE ! Case 7: F(b:e) = '...BinOp...'
848! WRITE(*,*)'7. Binary operator',F(j:j)
849 CALL compilesubstr(i, f, b, j - 1, var)
850 CALL compilesubstr(i, f, j + 1, e, var)
851 CALL addcompiledbyte(i, operatorindex(ops(io)))
852 comp(i)%StackPtr = comp(i)%StackPtr - 1
853 RETURN
854 END IF
855 END IF
856 END DO
857 END DO
858 !----- -------- --------- --------- --------- --------- --------- --------- -------
859 ! Check for remaining items, i.e. variables or explicit numbers
860 !----- -------- --------- --------- --------- --------- --------- --------- -------
861 b2 = b
862 IF (f(b:b) == '-') b2 = b2 + 1
863 n = mathitemindex(i, f(b2:e), var)
864! WRITE(*,*)'8. AddCompiledByte ',n
865 CALL addcompiledbyte(i, n)
866 comp(i)%StackPtr = comp(i)%StackPtr + 1
867 IF (comp(i)%StackPtr > comp(i)%StackSize) comp(i)%StackSize = comp(i)%StackSize + 1
868 IF (b2 > b) CALL addcompiledbyte(i, cneg)
869 END SUBROUTINE compilesubstr
870 !
871! **************************************************************************************************
872!> \brief ...
873!> \param i ...
874!> \param F ...
875!> \param b ...
876!> \param e ...
877!> \param Var ...
878! **************************************************************************************************
879 RECURSIVE SUBROUTINE compilefuncargssubstr(i, F, b, e, Var)
880 !----- -------- --------- --------- --------- --------- --------- --------- -------
881 ! Compile i-th function arguments substring F(b,e) into bytecode
882 !----- -------- --------- --------- --------- --------- --------- --------- -------
883 INTEGER, INTENT(in) :: i
884 CHARACTER(LEN=*), INTENT(in) :: f
885 INTEGER, INTENT(in) :: b, e
886 CHARACTER(LEN=*), DIMENSION(:), INTENT(in) :: var
887
888 INTEGER :: b2, j, parcnt
889
890! Function identifier
891! Function substring
892! Begin and end position substring
893! Array with variable names
894!----- -------- --------- --------- --------- --------- --------- --------- -------
895! Check for commas for splitting function arguments into substrings
896!----- -------- --------- --------- --------- --------- --------- --------- -------
897
898 parcnt = 0
899 b2 = b
900 DO j = b, e
901 IF (f(j:j) == '(') THEN
902 parcnt = parcnt + 1
903 ELSEIF (f(j:j) == ')') THEN
904 parcnt = parcnt - 1
905 ELSEIF (parcnt == 0 .AND. f(j:j) == ',') THEN
906 CALL compilesubstr(i, f, b2, j - 1, var)
907 b2 = j + 1
908 END IF
909 END DO
910 CALL compilesubstr(i, f, b2, e, var)
911 END SUBROUTINE compilefuncargssubstr
912 !
913! **************************************************************************************************
914!> \brief ...
915!> \param j ...
916!> \param F ...
917!> \return ...
918! **************************************************************************************************
919 FUNCTION isbinaryop(j, F) RESULT(res)
920 !----- -------- --------- --------- --------- --------- --------- --------- -------
921 ! Check if operator F(j:j) in string F is binary operator
922 ! Special cases already covered elsewhere: (that is corrected in v1.1)
923 ! - operator character F(j:j) is first character of string (j=1)
924 !----- -------- --------- --------- --------- --------- --------- --------- -------
925 INTEGER, INTENT(in) :: j
926 CHARACTER(LEN=*), INTENT(in) :: f
927 LOGICAL :: res
928
929 INTEGER :: k
930 LOGICAL :: dflag, pflag
931
932! Position of Operator
933! String
934! Result
935!----- -------- --------- --------- --------- --------- --------- --------- -------
936
937 res = .true.
938 IF (f(j:j) == '+' .OR. f(j:j) == '-') THEN ! Plus or minus sign:
939 IF (j == 1) THEN ! - leading unary operator ?
940 res = .false.
941 ELSEIF (scan(f(j - 1:j - 1), '+-*/^(,') > 0) THEN ! - other unary operator ?
942 res = .false.
943 ELSEIF (scan(f(j + 1:j + 1), '0123456789') > 0 .AND. & ! - in exponent of real number ?
944 scan(f(j - 1:j - 1), 'eEdD') > 0) THEN
945 dflag = .false.; pflag = .false.
946 k = j - 1
947 DO WHILE (k > 1) ! step to the left in mantissa
948 k = k - 1
949 IF (scan(f(k:k), '0123456789') > 0) THEN
950 dflag = .true.
951 ELSEIF (f(k:k) == '.') THEN
952 IF (pflag) THEN
953 EXIT ! * EXIT: 2nd appearance of '.'
954 ELSE
955 pflag = .true. ! * mark 1st appearance of '.'
956 END IF
957 ELSE
958 EXIT ! * all other characters
959 END IF
960 END DO
961 IF (dflag .AND. (k == 1 .OR. scan(f(k:k), '+-*/^(') > 0)) res = .false.
962 END IF
963 END IF
964 END FUNCTION isbinaryop
965 !
966! **************************************************************************************************
967!> \brief ...
968!> \param str ...
969!> \param ibegin ...
970!> \param inext ...
971!> \param error ...
972!> \return ...
973! **************************************************************************************************
974 FUNCTION realnum(str, ibegin, inext, error) RESULT(res)
975 !----- -------- --------- --------- --------- --------- --------- --------- -------
976 ! Get real number from string - Format: [blanks][+|-][nnn][.nnn][e|E|d|D[+|-]nnn]
977 !----- -------- --------- --------- --------- --------- --------- --------- -------
978 CHARACTER(LEN=*), INTENT(in) :: str
979 INTEGER, INTENT(out), OPTIONAL :: ibegin, inext
980 LOGICAL, INTENT(out), OPTIONAL :: error
981 REAL(rn) :: res
982
983 INTEGER :: ib, in, istat
984 LOGICAL :: bflag, dinexp, dinman, eflag, err, &
985 inexp, inman, pflag
986
987! String
988! Real number
989! Start position of real number
990! 1st character after real number
991! Error flag
992! .T. at begin of number in str
993! .T. in mantissa of number
994! .T. after 1st '.' encountered
995! .T. at exponent identifier 'eEdD'
996! .T. in exponent of number
997! .T. if at least 1 digit in mant.
998! .T. if at least 1 digit in exp.
999! Local error flag
1000!----- -------- --------- --------- --------- --------- --------- --------- -------
1001
1002 bflag = .true.; inman = .false.; pflag = .false.; eflag = .false.; inexp = .false.
1003 dinman = .false.; dinexp = .false.
1004 ib = 1
1005 in = 1
1006 DO WHILE (in <= len_trim(str))
1007 SELECT CASE (str(in:in))
1008 CASE (' ') ! Only leading blanks permitted
1009 ib = ib + 1
1010 IF (inman .OR. eflag .OR. inexp) EXIT
1011 CASE ('+', '-') ! Permitted only
1012 IF (bflag) THEN
1013 inman = .true.; bflag = .false. ! - at beginning of mantissa
1014 ELSEIF (eflag) THEN
1015 inexp = .true.; eflag = .false. ! - at beginning of exponent
1016 ELSE
1017 EXIT ! - otherwise STOP
1018 END IF
1019 CASE ('0':'9') ! Mark
1020 IF (bflag) THEN
1021 inman = .true.; bflag = .false. ! - beginning of mantissa
1022 ELSEIF (eflag) THEN
1023 inexp = .true.; eflag = .false. ! - beginning of exponent
1024 END IF
1025 IF (inman) dinman = .true. ! Mantissa contains digit
1026 IF (inexp) dinexp = .true. ! Exponent contains digit
1027 CASE ('.')
1028 IF (bflag) THEN
1029 pflag = .true. ! - mark 1st appearance of '.'
1030 inman = .true.; bflag = .false. ! mark beginning of mantissa
1031 ELSEIF (inman .AND. .NOT. pflag) THEN
1032 pflag = .true. ! - mark 1st appearance of '.'
1033 ELSE
1034 EXIT ! - otherwise STOP
1035 END IF
1036 CASE ('e', 'E', 'd', 'D') ! Permitted only
1037 IF (inman) THEN
1038 eflag = .true.; inman = .false. ! - following mantissa
1039 ELSE
1040 EXIT ! - otherwise STOP
1041 END IF
1042 CASE DEFAULT
1043 EXIT ! STOP at all other characters
1044 END SELECT
1045 in = in + 1
1046 END DO
1047 err = (ib > in - 1) .OR. (.NOT. dinman) .OR. ((eflag .OR. inexp) .AND. .NOT. dinexp)
1048 IF (err) THEN
1049 res = 0.0_rn
1050 ELSE
1051 READ (str(ib:in - 1), *, iostat=istat) res
1052 err = istat /= 0
1053 END IF
1054 IF (PRESENT(ibegin)) ibegin = ib
1055 IF (PRESENT(inext)) inext = in
1056 IF (PRESENT(error)) error = err
1057 END FUNCTION realnum
1058 !
1059! **************************************************************************************************
1060!> \brief ...
1061!> \param str1 ...
1062!> \param str2 ...
1063! **************************************************************************************************
1064 SUBROUTINE lowcase(str1, str2)
1065 !----- -------- --------- --------- --------- --------- --------- --------- -------
1066 ! Transform upper case letters in str1 into lower case letters, result is str2
1067 !----- -------- --------- --------- --------- --------- --------- --------- -------
1068 CHARACTER(LEN=*), INTENT(in) :: str1
1069 CHARACTER(LEN=*), INTENT(out) :: str2
1070
1071 CHARACTER(LEN=*), PARAMETER :: lc = 'abcdefghijklmnopqrstuvwxyz', &
1072 uc = 'ABCDEFGHIJKLMNOPQRSTUVWXYZ'
1073
1074 INTEGER :: j, k
1075
1076!----- -------- --------- --------- --------- --------- --------- --------- -------
1077
1078 str2 = str1
1079 DO j = 1, len_trim(str1)
1080 k = index(uc, str1(j:j))
1081 IF (k > 0) str2(j:j) = lc(k:k)
1082 END DO
1083 END SUBROUTINE lowcase
1084
1085! **************************************************************************************************
1086!> \brief Evaluates derivatives
1087!> \param id_fun ...
1088!> \param ipar ...
1089!> \param vals ...
1090!> \param h ...
1091!> \param err ...
1092!> \return ...
1093!> \author Main algorithm from Numerical Recipes
1094!> Ridders, C.J.F. 1982 - Advances in Engineering Software, Vol.4, no. 2, pp. 75-76
1095! **************************************************************************************************
1096 FUNCTION evalfd(id_fun, ipar, vals, h, err) RESULT(derivative)
1097 INTEGER, INTENT(IN) :: id_fun, ipar
1098 REAL(kind=rn), DIMENSION(:), INTENT(INOUT) :: vals
1099 REAL(kind=rn), INTENT(IN) :: h
1100 REAL(kind=rn), INTENT(OUT) :: err
1101 REAL(kind=rn) :: derivative
1102
1103 INTEGER, PARAMETER :: ntab = 10
1104 REAL(kind=rn), PARAMETER :: big_error = 1.0e30_rn, con = 1.4_rn, &
1105 con2 = con*con, safe = 2.0_rn
1106
1107 INTEGER :: i, j
1108 REAL(kind=rn) :: a(ntab, ntab), errt, fac, funcm, funcp, &
1109 hh, xval
1110
1111 derivative = huge(0.0_rn)
1112 IF (h /= 0._rn) THEN
1113 xval = vals(ipar)
1114 hh = h
1115 vals(ipar) = xval + hh
1116 funcp = evalf(id_fun, vals)
1117 vals(ipar) = xval - hh
1118 funcm = evalf(id_fun, vals)
1119 a(1, 1) = (funcp - funcm)/(2.0_rn*hh)
1120 err = big_error
1121 DO i = 2, ntab
1122 hh = hh/con
1123 vals(ipar) = xval + hh
1124 funcp = evalf(id_fun, vals)
1125 vals(ipar) = xval - hh
1126 funcm = evalf(id_fun, vals)
1127 a(1, i) = (funcp - funcm)/(2.0_rn*hh)
1128 fac = con2
1129 DO j = 2, i
1130 a(j, i) = (a(j - 1, i)*fac - a(j - 1, i - 1))/(fac - 1.0_rn)
1131 fac = con2*fac
1132 errt = max(abs(a(j, i) - a(j - 1, i)), abs(a(j, i) - a(j - 1, i - 1)))
1133 IF (errt <= err) THEN
1134 err = errt
1135 derivative = a(j, i)
1136 END IF
1137 END DO
1138 IF (abs(a(i, i) - a(i - 1, i - 1)) >= safe*err) RETURN
1139 END DO
1140 ELSE
1141 cpabort("DX provided equals zero!")
1142 END IF
1143 vals(ipar) = xval
1144 END FUNCTION evalfd
1145 !
1146! **************************************************************************************************
1147!> \brief ...
1148!> \return ...
1149! **************************************************************************************************
1150 FUNCTION docf() RESULT(doc)
1151 !----- -------- --------- --------- --------- --------- --------- --------- -------
1152 ! Returns documentation for function parser and its available operators and functions
1153 !----- -------- --------- --------- --------- --------- --------- --------- -------
1154 CHARACTER(LEN=:), ALLOCATABLE :: doc
1155
1156 INTEGER :: i
1157
1158!----- -------- --------- --------- --------- --------- --------- --------- -------
1159
1160 doc = "A functional form is specified. Mathematical Operators recognized are "
1161 DO i = cadd, cpow
1162 IF (i > cadd) doc = doc//", "
1163 doc = doc//ops(i)
1164 END DO
1165 doc = doc//" or alternatively **, whereas symbols for brackets must be (). The function"// &
1166 " parser recognizes the Fortran 90 intrinsic functions "
1167 DO i = cabs, cerfc
1168 IF (i > cabs) doc = doc//", "
1169 doc = doc//trim(funcs(i))
1170 END DO
1171 doc = doc//". Parsing for intrinsic functions is not case sensitive."
1172 END FUNCTION docf
1173
1174END MODULE fparser
fac
static GRID_HOST_DEVICE double fac(const int i)
Factorial function, e.g. fac(5) = 5! = 120.
Definition grid_common.h:51
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:174
fparser::evalf
real(rn) function, public evalf(i, val)
...
Definition fparser.F:206
fparser::evalerrtype
integer, public evalerrtype
Definition fparser.F:33
fparser::evalfd
real(kind=rn) function, public evalfd(id_fun, ipar, vals, h, err)
Evaluates derivatives.
Definition fparser.F:1097
fparser::docf
character(len=:) function, allocatable, public docf()
...
Definition fparser.F:1151
fparser::comp
type(tcomp), dimension(:), pointer comp
Definition fparser.F:118
fparser::finalizef
subroutine, public finalizef()
...
Definition fparser.F:127
fparser::initf
subroutine, public initf(n)
...
Definition fparser.F:156
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