Line data Source code
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 : ! **************************************************************************************************
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 : 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 : !
121 : CONTAINS
122 : !
123 : ! **************************************************************************************************
124 : !> \brief ...
125 : ! **************************************************************************************************
126 38884 : SUBROUTINE finalizef()
127 : !----- -------- --------- --------- --------- --------- --------- --------- -------
128 : ! Finalize function parser
129 : !----- -------- --------- --------- --------- --------- --------- --------- -------
130 :
131 : INTEGER :: i
132 :
133 : !----- -------- --------- --------- --------- --------- --------- --------- -------
134 :
135 64631 : DO i = 1, SIZE(Comp)
136 25747 : IF (ASSOCIATED(Comp(i)%ByteCode)) THEN
137 25747 : DEALLOCATE (Comp(i)%ByteCode)
138 : END IF
139 25747 : IF (ASSOCIATED(Comp(i)%Immed)) THEN
140 25747 : DEALLOCATE (Comp(i)%Immed)
141 : END IF
142 64631 : IF (ASSOCIATED(Comp(i)%Stack)) THEN
143 25747 : DEALLOCATE (Comp(i)%Stack)
144 : END IF
145 : END DO
146 :
147 38884 : DEALLOCATE (Comp)
148 :
149 38884 : END SUBROUTINE finalizef
150 : !
151 : ! **************************************************************************************************
152 : !> \brief ...
153 : !> \param n ...
154 : ! **************************************************************************************************
155 38884 : 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 109326 : ALLOCATE (Comp(n))
165 38884 : 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 51494 : 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 25747 : CHARACTER(LEN=LEN(FuncStr)) :: Func
179 :
180 : ! Function string, local use
181 : !----- -------- --------- --------- --------- --------- --------- --------- -------
182 :
183 25747 : IF (i < 1 .OR. i > SIZE(Comp)) THEN
184 0 : CPABORT("Function number is out of range")
185 : END IF
186 25747 : IF (SIZE(var) > HUGE(0_is)) THEN
187 0 : CPABORT("Too many variables")
188 : END IF
189 77241 : ALLOCATE (ipos(LEN_TRIM(FuncStr))) ! Char. positions in orig. string
190 25747 : Func = FuncStr ! Local copy of function string
191 25747 : CALL Replace('**', '^ ', Func) ! Exponent into 1-Char. format
192 25747 : CALL RemoveSpaces(Func) ! Condense function string
193 25747 : CALL CheckSyntax(Func, FuncStr, Var)
194 25747 : DEALLOCATE (ipos)
195 25747 : CALL Compile(i, Func, Var) ! Compile into bytecode
196 :
197 25747 : END SUBROUTINE parsef
198 : !
199 : ! **************************************************************************************************
200 : !> \brief ...
201 : !> \param i ...
202 : !> \param Val ...
203 : !> \return ...
204 : ! **************************************************************************************************
205 51475855 : 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 51475855 : DP = 1
226 51475855 : SP = 0
227 1002333048 : DO IP = 1, Comp(i)%ByteCodeSize
228 51475845 : SELECT CASE (Comp(i)%ByteCode(IP))
229 :
230 150630693 : CASE (cImmed); SP = SP + 1; Comp(i)%Stack(SP) = Comp(i)%Immed(DP); DP = DP + 1
231 2762941 : CASE (cNeg); Comp(i)%Stack(SP) = -Comp(i)%Stack(SP)
232 97290043 : CASE (cAdd); Comp(i)%Stack(SP - 1) = Comp(i)%Stack(SP - 1) + Comp(i)%Stack(SP); SP = SP - 1
233 49374406 : CASE (cSub); Comp(i)%Stack(SP - 1) = Comp(i)%Stack(SP - 1) - Comp(i)%Stack(SP); SP = SP - 1
234 100721850 : CASE (cMul); Comp(i)%Stack(SP - 1) = Comp(i)%Stack(SP - 1)*Comp(i)%Stack(SP); SP = SP - 1
235 98846078 : CASE (cDiv); IF (Comp(i)%Stack(SP) == 0._rn) THEN; EvalErrType = 1; res = zero; RETURN; END IF
236 98846068 : 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 100777220 : IF (Comp(i)%Stack(SP - 1) < 0.0_rn) THEN
240 321252 : ipow = FLOOR(Comp(i)%Stack(SP))
241 321252 : IF (MOD(Comp(i)%Stack(SP), REAL(ipow, KIND=rn)) == 0.0_rn) THEN
242 321252 : Comp(i)%Stack(SP - 1) = Comp(i)%Stack(SP - 1)**ipow
243 : ELSE
244 0 : CPABORT("Negative floating-point value raised to a real power!")
245 : END IF
246 : ELSE
247 100455968 : Comp(i)%Stack(SP - 1) = Comp(i)%Stack(SP - 1)**Comp(i)%Stack(SP)
248 : END IF
249 0 : SP = SP - 1
250 0 : CASE (cAbs); Comp(i)%Stack(SP) = ABS(Comp(i)%Stack(SP))
251 0 : CASE (cMin); Comp(i)%Stack(SP - 1) = MIN(Comp(i)%Stack(SP - 1), Comp(i)%Stack(SP)); SP = SP - 1
252 0 : CASE (cMax); Comp(i)%Stack(SP - 1) = MAX(Comp(i)%Stack(SP - 1), Comp(i)%Stack(SP)); SP = SP - 1
253 2498087 : CASE (cExp); Comp(i)%Stack(SP) = EXP(Comp(i)%Stack(SP))
254 0 : CASE (cLog10); IF (Comp(i)%Stack(SP) <= 0._rn) THEN; EvalErrType = 3; res = zero; RETURN; END IF
255 0 : Comp(i)%Stack(SP) = LOG10(Comp(i)%Stack(SP))
256 0 : CASE (cLog); IF (Comp(i)%Stack(SP) <= 0._rn) THEN; EvalErrType = 3; res = zero; RETURN; END IF
257 0 : Comp(i)%Stack(SP) = LOG(Comp(i)%Stack(SP))
258 0 : CASE (cSqrt); IF (Comp(i)%Stack(SP) < 0._rn) THEN; EvalErrType = 3; res = zero; RETURN; END IF
259 0 : Comp(i)%Stack(SP) = SQRT(Comp(i)%Stack(SP))
260 0 : CASE (cSinh); Comp(i)%Stack(SP) = SINH(Comp(i)%Stack(SP))
261 0 : CASE (cCosh); Comp(i)%Stack(SP) = COSH(Comp(i)%Stack(SP))
262 0 : CASE (cTanh); Comp(i)%Stack(SP) = TANH(Comp(i)%Stack(SP))
263 84816 : CASE (cSin); Comp(i)%Stack(SP) = SIN(Comp(i)%Stack(SP))
264 13206 : CASE (cCos); Comp(i)%Stack(SP) = COS(Comp(i)%Stack(SP))
265 0 : CASE (cTan); Comp(i)%Stack(SP) = TAN(Comp(i)%Stack(SP))
266 0 : CASE (cAsin); IF ((Comp(i)%Stack(SP) < -1._rn) .OR. (Comp(i)%Stack(SP) > 1._rn)) THEN
267 0 : EvalErrType = 4; res = zero; RETURN; END IF
268 0 : Comp(i)%Stack(SP) = ASIN(Comp(i)%Stack(SP))
269 1548 : CASE (cAcos); IF ((Comp(i)%Stack(SP) < -1._rn) .OR. (Comp(i)%Stack(SP) > 1._rn)) THEN
270 0 : EvalErrType = 4; res = zero; RETURN; END IF
271 1548 : Comp(i)%Stack(SP) = ACOS(Comp(i)%Stack(SP))
272 542 : CASE (cAtan2); Comp(i)%Stack(SP - 1) = ATAN2(Comp(i)%Stack(SP - 1), Comp(i)%Stack(SP)); SP = SP - 1
273 472 : CASE (cAtan); Comp(i)%Stack(SP) = ATAN(Comp(i)%Stack(SP))
274 0 : CASE (cErf); Comp(i)%Stack(SP) = ERF(Comp(i)%Stack(SP))
275 0 : CASE (cErfc); Comp(i)%Stack(SP) = ERFC(Comp(i)%Stack(SP))
276 950857203 : CASE DEFAULT; SP = SP + 1; Comp(i)%Stack(SP) = Val(Comp(i)%ByteCode(IP) - VarBegin + 1)
277 : END SELECT
278 : END DO
279 51475845 : EvalErrType = 0
280 51475845 : res = Comp(i)%Stack(1)
281 51475845 : END FUNCTION evalf
282 : !
283 : ! **************************************************************************************************
284 : !> \brief ...
285 : !> \param Func ...
286 : !> \param FuncStr ...
287 : !> \param Var ...
288 : ! **************************************************************************************************
289 25747 : 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 25747 : j = 1
309 25747 : ParCnt = 0
310 25747 : lFunc = LEN_TRIM(Func)
311 : step: DO
312 340187 : IF (j > lFunc) CALL ParseErrMsg(j, FuncStr)
313 340187 : c = Func(j:j)
314 : !-- -------- --------- --------- --------- --------- --------- --------- -------
315 : ! Check for valid operand (must appear)
316 : !-- -------- --------- --------- --------- --------- --------- --------- -------
317 340187 : IF (c == '-' .OR. c == '+') THEN ! Check for leading - or +
318 969 : j = j + 1
319 969 : IF (j > lFunc) CALL ParseErrMsg(j, FuncStr, 'Missing operand')
320 969 : c = Func(j:j)
321 5814 : IF (ANY(c == Ops)) CALL ParseErrMsg(j, FuncStr, 'Multiple operators')
322 : END IF
323 340187 : n = MathFunctionIndex(Func(j:))
324 340187 : IF (n > 0) THEN ! Check for math function
325 825 : j = j + LEN_TRIM(Funcs(n))
326 825 : IF (j > lFunc) CALL ParseErrMsg(j, FuncStr, 'Missing function argument')
327 825 : c = Func(j:j)
328 825 : IF (c /= '(') CALL ParseErrMsg(j, FuncStr, 'Missing opening parenthesis')
329 825 : CALL CheckFuncArgCnt(Func, FuncStr, j, lFunc, n)
330 : END IF
331 340187 : IF (c == '(') THEN ! Check for opening parenthesis
332 110461 : ParCnt = ParCnt + 1
333 110461 : j = j + 1
334 110461 : CYCLE step
335 : END IF
336 229726 : IF (SCAN(c, '0123456789.') > 0) THEN ! Check for number
337 70628 : r = RealNum(Func(j:), ib, in, err)
338 70628 : IF (err) CALL ParseErrMsg(j, FuncStr, 'Invalid number format: '//Func(j + ib - 1:j + in - 2))
339 70628 : j = j + in - 1
340 70628 : IF (j > lFunc) EXIT step
341 70099 : c = Func(j:j)
342 : ELSE ! Check for variable
343 159098 : n = VariableIndex(Func(j:), Var, ib, in)
344 159098 : IF (n == 0) CALL ParseErrMsg(j, FuncStr, 'Invalid element: '//Func(j + ib - 1:j + in - 2))
345 159098 : j = j + in - 1
346 159098 : IF (j > lFunc) EXIT step
347 155894 : c = Func(j:j)
348 : END IF
349 314440 : DO WHILE (c == ')') ! Check for closing parenthesis
350 110461 : ParCnt = ParCnt - 1
351 110461 : IF (ParCnt < 0) CALL ParseErrMsg(j, FuncStr, 'Mismatched parenthesis')
352 110461 : IF (Func(j - 1:j - 1) == '(') CALL ParseErrMsg(j - 1, FuncStr, 'Empty parentheses')
353 110461 : j = j + 1
354 110461 : IF (j > lFunc) EXIT
355 314440 : 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 225993 : IF (j > lFunc) EXIT step
361 631173 : IF (ANY(c == Ops)) THEN ! Check for multiple operators
362 203957 : IF (j + 1 > lFunc) CALL ParseErrMsg(j, FuncStr)
363 1223742 : IF (ANY(Func(j + 1:j + 1) == Ops)) CALL ParseErrMsg(j + 1, FuncStr, 'Multiple operators')
364 22 : ELSE IF (c == ',') THEN ! Check for commas separating function arguments
365 22 : IF (ParCnt == 0) CALL ParseErrMsg(j, FuncStr, 'Comma outside function')
366 22 : IF (Func(j + 1:j + 1) == ',') CALL ParseErrMsg(j, FuncStr, 'Multiple commas')
367 : ELSE ! Check for next operand
368 0 : 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 229197 : j = j + 1
375 : END DO step
376 25747 : IF (ParCnt > 0) CALL ParseErrMsg(j, FuncStr, 'Missing )')
377 25747 : END SUBROUTINE CheckSyntax
378 : !
379 : ! **************************************************************************************************
380 : !> \brief ...
381 : !> \param Func ...
382 : !> \param FuncStr ...
383 : !> \param b ...
384 : !> \param e ...
385 : !> \param FuncId ...
386 : ! **************************************************************************************************
387 825 : 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 825 : ArgCnt = 1
405 825 : ArgPos = 0
406 825 : ParCnt = 0
407 9218 : DO j = b, e
408 10043 : IF (Func(j:j) == '(') THEN
409 1329 : ParCnt = ParCnt + 1
410 7889 : ELSEIF (Func(j:j) == ')') THEN
411 1329 : ParCnt = ParCnt - 1
412 1329 : IF (ParCnt == 0) EXIT
413 6560 : ELSEIF (ParCnt == 1 .AND. Func(j:j) == ',') THEN
414 22 : ArgPos = j
415 22 : ArgCnt = ArgCnt + 1
416 : END IF
417 : END DO
418 825 : IF (ArgCnt /= FuncsArgCnt(FuncId)) THEN
419 0 : IF (ArgCnt < FuncsArgCnt(FuncId)) ArgPos = j
420 0 : WRITE (Msg, '(I0,A,A,A,I0)') ArgCnt, ' argument(s) in ', TRIM(Funcs(FuncId)), ' instead of ', FuncsArgCnt(FuncId)
421 0 : CALL ParseErrMsg(ArgPos, FuncStr, Msg)
422 : END IF
423 825 : END SUBROUTINE CheckFuncArgCnt
424 : !
425 : ! **************************************************************************************************
426 : !> \brief ...
427 : !> \return ...
428 : ! **************************************************************************************************
429 0 : 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 0 : IF (EvalErrType < 1 .OR. EvalErrType > SIZE(m)) THEN
441 0 : msg = ''
442 : ELSE
443 0 : msg = m(EvalErrType)
444 : END IF
445 0 : END FUNCTION EvalErrMsg
446 : !
447 : ! **************************************************************************************************
448 : !> \brief ...
449 : !> \param j ...
450 : !> \param FuncStr ...
451 : !> \param Msg ...
452 : ! **************************************************************************************************
453 0 : 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 0 : IF (PRESENT(Msg)) THEN
467 0 : WRITE (UNIT=message, FMT="(A)") "Syntax error in function string: "//Msg
468 : ELSE
469 0 : WRITE (UNIT=message, FMT="(A)") "Syntax error in function string"
470 : END IF
471 0 : WRITE (*, '(/,T2,A)') TRIM(FuncStr)
472 0 : IF ((j > LBOUND(ipos, DIM=1)) .AND. (j <= UBOUND(ipos, DIM=1))) THEN
473 0 : WRITE (*, '(A)') REPEAT(" ", ipos(j))//"?"
474 : ELSE
475 0 : WRITE (*, '(A)') REPEAT(" ", SIZE(ipos) + 1)//"?"
476 : END IF
477 0 : CPABORT(TRIM(message))
478 :
479 0 : END SUBROUTINE ParseErrMsg
480 : !
481 : ! **************************************************************************************************
482 : !> \brief ...
483 : !> \param c ...
484 : !> \return ...
485 : ! **************************************************************************************************
486 407914 : 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 407914 : n = 0
498 1262082 : DO j = cAdd, cPow
499 1262082 : IF (c == Ops(j)) THEN
500 : n = j
501 : EXIT
502 : END IF
503 : END DO
504 407914 : END FUNCTION OperatorIndex
505 : !
506 : ! **************************************************************************************************
507 : !> \brief ...
508 : !> \param str ...
509 : !> \return ...
510 : ! **************************************************************************************************
511 841657 : 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 841657 : n = 0
525 16798740 : DO j = cAbs, cErfc ! Check all math functions
526 15960294 : k = MIN(LEN_TRIM(Funcs(j)), LEN(str))
527 15960294 : CALL LowCase(str(1:k), fun)
528 16798740 : 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 841657 : END FUNCTION MathFunctionIndex
534 : !
535 : ! **************************************************************************************************
536 : !> \brief ...
537 : !> \param str ...
538 : !> \param Var ...
539 : !> \param ibegin ...
540 : !> \param inext ...
541 : !> \return ...
542 : ! **************************************************************************************************
543 477294 : 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 477294 : n = 0
562 477294 : lstr = LEN_TRIM(str)
563 477294 : IF (lstr > 0) THEN
564 477294 : DO ib = 1, lstr ! Search for first character in str
565 477294 : IF (str(ib:ib) /= ' ') EXIT ! When lstr>0 at least 1 char in str
566 : END DO
567 2384220 : DO in = ib, lstr ! Search for name terminators
568 2384220 : IF (SCAN(str(in:in), '+-*/^) ,') > 0) EXIT
569 : END DO
570 1092801 : DO j = 1, SIZE(Var)
571 1092801 : IF (str(ib:in - 1) == Var(j)) THEN
572 477294 : n = INT(j, KIND=is) ! Variable name found
573 477294 : EXIT
574 : END IF
575 : END DO
576 : END IF
577 477294 : IF (PRESENT(ibegin)) ibegin = ib
578 477294 : IF (PRESENT(inext)) inext = in
579 477294 : END FUNCTION VariableIndex
580 : !
581 : ! **************************************************************************************************
582 : !> \brief ...
583 : !> \param str ...
584 : ! **************************************************************************************************
585 25747 : 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 25747 : lstr = LEN_TRIM(str)
596 2435455 : ipos(:) = [(k, k=1, lstr)]
597 25747 : k = 1
598 1230601 : DO WHILE (str(k:lstr) /= ' ')
599 1204854 : IF (str(k:k) == ' ') THEN
600 43492 : str(k:lstr) = str(k + 1:lstr)//' ' ! Move 1 character to left
601 1267300 : ipos(k:lstr) = [ipos(k + 1:lstr), 0] ! Move 1 element to left
602 43492 : k = k - 1
603 : END IF
604 1204854 : k = k + 1
605 : END DO
606 25747 : END SUBROUTINE RemoveSpaces
607 : !
608 : ! **************************************************************************************************
609 : !> \brief ...
610 : !> \param ca ...
611 : !> \param cb ...
612 : !> \param str ...
613 : ! **************************************************************************************************
614 25747 : 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 25747 : lca = LEN(ca)
628 1204854 : DO j = 1, LEN_TRIM(str) - lca + 1
629 1204854 : IF (str(j:j + lca - 1) == ca) str(j:j + lca - 1) = cb
630 : END DO
631 25747 : END SUBROUTINE Replace
632 : !
633 : ! **************************************************************************************************
634 : !> \brief ...
635 : !> \param i ...
636 : !> \param F ...
637 : !> \param Var ...
638 : ! **************************************************************************************************
639 25747 : 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 25747 : IF (ASSOCIATED(Comp(i)%ByteCode)) DEALLOCATE (Comp(i)%ByteCode, &
653 0 : Comp(i)%Immed, &
654 0 : Comp(i)%Stack)
655 25747 : Comp(i)%ByteCodeSize = 0
656 25747 : Comp(i)%ImmedSize = 0
657 25747 : Comp(i)%StackSize = 0
658 25747 : Comp(i)%StackPtr = 0
659 25747 : 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 178568 : Comp(i)%Stack(Comp(i)%StackSize))
663 25747 : Comp(i)%ByteCodeSize = 0
664 25747 : Comp(i)%ImmedSize = 0
665 25747 : Comp(i)%StackSize = 0
666 25747 : Comp(i)%StackPtr = 0
667 25747 : CALL CompileSubstr(i, F, 1, LEN_TRIM(F), Var) ! Compile string into bytecode
668 : !
669 25747 : END SUBROUTINE Compile
670 : !
671 : ! **************************************************************************************************
672 : !> \brief ...
673 : !> \param i ...
674 : !> \param b ...
675 : ! **************************************************************************************************
676 870954 : 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 870954 : Comp(i)%ByteCodeSize = Comp(i)%ByteCodeSize + 1
688 870954 : IF (ASSOCIATED(Comp(i)%ByteCode)) Comp(i)%ByteCode(Comp(i)%ByteCodeSize) = b
689 870954 : END SUBROUTINE AddCompiledByte
690 : !
691 : ! **************************************************************************************************
692 : !> \brief ...
693 : !> \param i ...
694 : !> \param F ...
695 : !> \param Var ...
696 : !> \return ...
697 : ! **************************************************************************************************
698 459452 : 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 459452 : n = 0
714 459452 : IF (SCAN(F(1:1), '0123456789.') > 0) THEN ! Check for begin of a number
715 141256 : Comp(i)%ImmedSize = Comp(i)%ImmedSize + 1
716 141256 : IF (ASSOCIATED(Comp(i)%Immed)) Comp(i)%Immed(Comp(i)%ImmedSize) = RealNum(F)
717 : n = cImmed
718 : ELSE ! Check for a variable
719 318196 : n = VariableIndex(F, Var)
720 318196 : IF (n > 0) n = VarBegin + n - 1_is
721 : END IF
722 459452 : END FUNCTION MathItemIndex
723 : !
724 : ! **************************************************************************************************
725 : !> \brief ...
726 : !> \param F ...
727 : !> \param b ...
728 : !> \param e ...
729 : !> \return ...
730 : ! **************************************************************************************************
731 1093314 : 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 1093314 : res = .FALSE.
746 1093314 : IF (F(b:b) == '(' .AND. F(e:e) == ')') THEN
747 221184 : k = 0
748 3898776 : DO j = b + 1, e - 1
749 3677854 : IF (F(j:j) == '(') THEN
750 261764 : k = k + 1
751 3416090 : ELSEIF (F(j:j) == ')') THEN
752 262026 : k = k - 1
753 : END IF
754 3898776 : IF (k < 0) EXIT
755 : END DO
756 221184 : IF (k == 0) res = .TRUE. ! All opened parenthesis closed
757 : END IF
758 1093314 : END FUNCTION CompletelyEnclosed
759 : !
760 : ! **************************************************************************************************
761 : !> \brief ...
762 : !> \param i ...
763 : !> \param F ...
764 : !> \param b ...
765 : !> \param e ...
766 : !> \param Var ...
767 : ! **************************************************************************************************
768 1088942 : 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 1088942 : IF (F(b:b) == '+') THEN ! Case 1: F(b:e) = '+...'
792 : ! WRITE(*,*)'1. F(b:e) = "+..."'
793 0 : CALL CompileSubstr(i, F, b + 1, e, Var)
794 629490 : RETURN
795 1088942 : ELSEIF (CompletelyEnclosed(F, b, e)) THEN ! Case 2: F(b:e) = '(...)'
796 : ! WRITE(*,*)'2. F(b:e) = "(...)"'
797 219152 : CALL CompileSubstr(i, F, b + 1, e - 1, Var)
798 219152 : RETURN
799 869790 : ELSEIF (SCAN(F(b:b), calpha) > 0) THEN
800 499964 : n = MathFunctionIndex(F(b:e))
801 499964 : IF (n > 0) THEN
802 2386 : b2 = b + INDEX(F(b:e), '(') - 1
803 2386 : IF (CompletelyEnclosed(F, b2, e)) THEN ! Case 3: F(b:e) = 'fcn(...)'
804 : ! WRITE(*,*)'3. F(b:e) = "fcn(...)"'
805 1650 : CALL CompileFuncArgsSubstr(i, F, b2 + 1, e - 1, Var)
806 1650 : CALL AddCompiledByte(i, n)
807 1650 : RETURN
808 : END IF
809 : END IF
810 369826 : ELSEIF (F(b:b) == '-') THEN
811 1986 : IF (CompletelyEnclosed(F, b + 1, e)) THEN ! Case 4: F(b:e) = '-(...)'
812 : ! WRITE(*,*)'4. F(b:e) = "-(...)"'
813 120 : CALL CompileSubstr(i, F, b + 2, e - 1, Var)
814 120 : CALL AddCompiledByte(i, cNeg)
815 120 : RETURN
816 1866 : ELSEIF (SCAN(F(b + 1:b + 1), calpha) > 0) THEN
817 1506 : n = MathFunctionIndex(F(b + 1:e))
818 1506 : IF (n > 0) THEN
819 0 : b2 = b + INDEX(F(b + 1:e), '(')
820 0 : IF (CompletelyEnclosed(F, b2, e)) THEN ! Case 5: F(b:e) = '-fcn(...)'
821 : ! WRITE(*,*)'5. F(b:e) = "-fcn(...)"'
822 0 : CALL CompileFuncArgsSubstr(i, F, b2 + 1, e - 1, Var)
823 0 : CALL AddCompiledByte(i, n)
824 0 : CALL AddCompiledByte(i, cNeg)
825 0 : 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 4021080 : DO io = cAdd, cPow ! Increasing priority +-*/^
834 3561628 : k = 0
835 33603510 : DO j = e, b, -1
836 29990998 : IF (F(j:j) == ')') THEN
837 1964974 : k = k + 1
838 28026024 : ELSEIF (F(j:j) == '(') THEN
839 1964974 : k = k - 1
840 : END IF
841 33144058 : IF (k == 0 .AND. F(j:j) == Ops(io) .AND. IsBinaryOp(j, F)) THEN
842 1091380 : 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 654 : CALL CompileSubstr(i, F, b + 1, e, Var)
845 654 : CALL AddCompiledByte(i, cNeg)
846 654 : RETURN
847 : ELSE ! Case 7: F(b:e) = '...BinOp...'
848 : ! WRITE(*,*)'7. Binary operator',F(j:j)
849 407914 : CALL CompileSubstr(i, F, b, j - 1, Var)
850 407914 : CALL CompileSubstr(i, F, j + 1, e, Var)
851 407914 : CALL AddCompiledByte(i, OperatorIndex(Ops(io)))
852 407914 : Comp(i)%StackPtr = Comp(i)%StackPtr - 1
853 407914 : 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 459452 : b2 = b
862 459452 : IF (F(b:b) == '-') b2 = b2 + 1
863 459452 : n = MathItemIndex(i, F(b2:e), Var)
864 : ! WRITE(*,*)'8. AddCompiledByte ',n
865 459452 : CALL AddCompiledByte(i, n)
866 459452 : Comp(i)%StackPtr = Comp(i)%StackPtr + 1
867 459452 : IF (Comp(i)%StackPtr > Comp(i)%StackSize) Comp(i)%StackSize = Comp(i)%StackSize + 1
868 459452 : 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 1650 : 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 1650 : ParCnt = 0
899 1650 : b2 = b
900 16786 : DO j = b, e
901 16786 : IF (F(j:j) == '(') THEN
902 1008 : ParCnt = ParCnt + 1
903 14128 : ELSEIF (F(j:j) == ')') THEN
904 1008 : ParCnt = ParCnt - 1
905 13120 : ELSEIF (ParCnt == 0 .AND. F(j:j) == ',') THEN
906 44 : CALL CompileSubstr(i, F, b2, j - 1, Var)
907 44 : b2 = j + 1
908 : END IF
909 : END DO
910 1650 : CALL CompileSubstr(i, F, b2, e, Var)
911 1650 : END SUBROUTINE CompileFuncArgsSubstr
912 : !
913 : ! **************************************************************************************************
914 : !> \brief ...
915 : !> \param j ...
916 : !> \param F ...
917 : !> \return ...
918 : ! **************************************************************************************************
919 410386 : 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 410386 : res = .TRUE.
938 410386 : IF (F(j:j) == '+' .OR. F(j:j) == '-') THEN ! Plus or minus sign:
939 140044 : IF (j == 1) THEN ! - leading unary operator ?
940 : res = .FALSE.
941 138860 : ELSEIF (SCAN(F(j - 1:j - 1), '+-*/^(,') > 0) THEN ! - other unary operator ?
942 : res = .FALSE.
943 138226 : 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 0 : DO WHILE (k > 1) ! step to the left in mantissa
948 0 : k = k - 1
949 0 : IF (SCAN(F(k:k), '0123456789') > 0) THEN
950 : Dflag = .TRUE.
951 0 : ELSEIF (F(k:k) == '.') THEN
952 0 : 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 0 : IF (Dflag .AND. (k == 1 .OR. SCAN(F(k:k), '+-*/^(') > 0)) res = .FALSE.
962 : END IF
963 : END IF
964 410386 : END FUNCTION IsBinaryOp
965 : !
966 : ! **************************************************************************************************
967 : !> \brief ...
968 : !> \param str ...
969 : !> \param ibegin ...
970 : !> \param inext ...
971 : !> \param error ...
972 : !> \return ...
973 : ! **************************************************************************************************
974 141256 : 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 141256 : Bflag = .TRUE.; InMan = .FALSE.; Pflag = .FALSE.; Eflag = .FALSE.; InExp = .FALSE.
1003 141256 : DInMan = .FALSE.; DInExp = .FALSE.
1004 141256 : ib = 1
1005 141256 : in = 1
1006 335802 : DO WHILE (in <= LEN_TRIM(str))
1007 264645 : SELECT CASE (str(in:in))
1008 : CASE (' ') ! Only leading blanks permitted
1009 0 : ib = ib + 1
1010 0 : IF (InMan .OR. Eflag .OR. InExp) EXIT
1011 : CASE ('+', '-') ! Permitted only
1012 23746 : IF (Bflag) THEN
1013 : InMan = .TRUE.; Bflag = .FALSE. ! - at beginning of mantissa
1014 23746 : 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 190644 : IF (Bflag) THEN
1021 : InMan = .TRUE.; Bflag = .FALSE. ! - beginning of mantissa
1022 49388 : ELSEIF (Eflag) THEN
1023 0 : InExp = .TRUE.; Eflag = .FALSE. ! - beginning of exponent
1024 : END IF
1025 190644 : IF (InMan) DInMan = .TRUE. ! Mantissa contains digit
1026 190644 : IF (InExp) DInExp = .TRUE. ! Exponent contains digit
1027 : CASE ('.')
1028 3902 : IF (Bflag) THEN
1029 : Pflag = .TRUE. ! - mark 1st appearance of '.'
1030 : InMan = .TRUE.; Bflag = .FALSE. ! mark beginning of mantissa
1031 3902 : 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 0 : IF (InMan) THEN
1038 : Eflag = .TRUE.; InMan = .FALSE. ! - following mantissa
1039 : ELSE
1040 : EXIT ! - otherwise STOP
1041 : END IF
1042 : CASE DEFAULT
1043 264645 : EXIT ! STOP at all other characters
1044 : END SELECT
1045 218292 : in = in + 1
1046 : END DO
1047 141256 : err = (ib > in - 1) .OR. (.NOT. DInMan) .OR. ((Eflag .OR. InExp) .AND. .NOT. DInExp)
1048 : IF (err) THEN
1049 0 : res = 0.0_rn
1050 : ELSE
1051 141256 : READ (str(ib:in - 1), *, IOSTAT=istat) res
1052 141256 : err = istat /= 0
1053 : END IF
1054 141256 : IF (PRESENT(ibegin)) ibegin = ib
1055 141256 : IF (PRESENT(inext)) inext = in
1056 141256 : IF (PRESENT(error)) error = err
1057 141256 : END FUNCTION RealNum
1058 : !
1059 : ! **************************************************************************************************
1060 : !> \brief ...
1061 : !> \param str1 ...
1062 : !> \param str2 ...
1063 : ! **************************************************************************************************
1064 15960294 : 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 15960294 : str2 = str1
1079 67947678 : DO j = 1, LEN_TRIM(str1)
1080 51987384 : k = INDEX(uc, str1(j:j))
1081 67947678 : IF (k > 0) str2(j:j) = lc(k:k)
1082 : END DO
1083 15960294 : 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 3231 : 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 3231 : derivative = HUGE(0.0_rn)
1112 3231 : IF (h /= 0._rn) THEN
1113 3231 : xval = vals(ipar)
1114 3231 : hh = h
1115 3231 : vals(ipar) = xval + hh
1116 3231 : funcp = evalf(id_fun, vals)
1117 3231 : vals(ipar) = xval - hh
1118 3231 : funcm = evalf(id_fun, vals)
1119 3231 : a(1, 1) = (funcp - funcm)/(2.0_rn*hh)
1120 3231 : err = big_error
1121 7901 : DO i = 2, ntab
1122 7901 : hh = hh/con
1123 7901 : vals(ipar) = xval + hh
1124 7901 : funcp = evalf(id_fun, vals)
1125 7901 : vals(ipar) = xval - hh
1126 7901 : funcm = evalf(id_fun, vals)
1127 7901 : a(1, i) = (funcp - funcm)/(2.0_rn*hh)
1128 7901 : fac = con2
1129 25145 : DO j = 2, i
1130 17244 : a(j, i) = (a(j - 1, i)*fac - a(j - 1, i - 1))/(fac - 1.0_rn)
1131 17244 : fac = con2*fac
1132 17244 : errt = MAX(ABS(a(j, i) - a(j - 1, i)), ABS(a(j, i) - a(j - 1, i - 1)))
1133 25145 : IF (errt <= err) THEN
1134 7713 : err = errt
1135 7713 : derivative = a(j, i)
1136 : END IF
1137 : END DO
1138 7901 : IF (ABS(a(i, i) - a(i - 1, i - 1)) >= safe*err) RETURN
1139 : END DO
1140 : ELSE
1141 0 : CPABORT("DX provided equals zero!")
1142 : END IF
1143 0 : vals(ipar) = xval
1144 0 : END FUNCTION evalfd
1145 : !
1146 : ! **************************************************************************************************
1147 : !> \brief ...
1148 : !> \return ...
1149 : ! **************************************************************************************************
1150 57100 : 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 57100 : doc = doc//"A functional form is specified. Mathematical Operators recognized are "
1161 342600 : DO i = cAdd, cPow
1162 513900 : IF (i > cAdd) doc = doc//", "
1163 342600 : doc = doc//Ops(i)
1164 : END DO
1165 : doc = doc//" or alternatively **, whereas symbols for brackets must be (). The function"// &
1166 57100 : " parser recognizes the Fortran 90 intrinsic functions "
1167 1142000 : DO i = cAbs, cErfc
1168 2112700 : IF (i > cAbs) doc = doc//", "
1169 1142000 : doc = doc//TRIM(Funcs(i))
1170 : END DO
1171 57100 : doc = doc//". Parsing for intrinsic functions is not case sensitive."
1172 57100 : END FUNCTION docf
1173 :
1174 2797900 : END MODULE fparser
|
