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