20 #include "../base/base_uses.f90"
24 CHARACTER(len=*),
PARAMETER,
PRIVATE :: moduleN =
'fparser'
36 INTEGER,
PARAMETER,
PRIVATE :: is = selected_int_kind(1)
37 INTEGER(is),
PARAMETER :: cimmed = 1, &
61 CHARACTER(LEN=1),
DIMENSION(cAdd:cPow),
PARAMETER :: ops = (/
'+', &
66 CHARACTER(LEN=5),
DIMENSION(cAbs:cErfc),
PARAMETER :: funcs = (/
'abs ', &
84 INTEGER(is),
DIMENSION(:),
POINTER :: bytecode => null()
85 INTEGER :: bytecodesize = -1
86 REAL(rn),
DIMENSION(:),
POINTER :: immed => null()
87 INTEGER :: immedsize = -1
88 REAL(rn),
DIMENSION(:),
POINTER :: stack => null()
89 INTEGER :: stacksize = -1, &
92 TYPE(tcomp),
DIMENSION(:),
POINTER :: comp => null()
93 INTEGER,
DIMENSION(:),
ALLOCATABLE :: ipos
110 IF (
ASSOCIATED(comp(i)%ByteCode))
THEN
111 DEALLOCATE (comp(i)%ByteCode)
113 IF (
ASSOCIATED(comp(i)%Immed))
THEN
114 DEALLOCATE (comp(i)%Immed)
116 IF (
ASSOCIATED(comp(i)%Stack))
THEN
117 DEALLOCATE (comp(i)%Stack)
133 INTEGER,
INTENT(in) :: n
148 INTEGER,
INTENT(in) :: i
149 CHARACTER(LEN=*),
INTENT(in) :: funcstr
150 CHARACTER(LEN=*),
DIMENSION(:),
INTENT(in) :: var
152 CHARACTER(LEN=LEN(FuncStr)) :: func
157 IF (i < 1 .OR. i >
SIZE(comp))
THEN
158 cpabort(
"Function number is out of range")
160 IF (
SIZE(var) > huge(0_is))
THEN
161 cpabort(
"Too many variables")
163 ALLOCATE (ipos(len_trim(funcstr)))
165 CALL replace(
'**',
'^ ', func)
166 CALL removespaces(func)
167 CALL checksyntax(func, funcstr, var)
169 CALL compile(i, func, var)
183 INTEGER,
INTENT(in) :: i
184 REAL(rn),
DIMENSION(:),
INTENT(in) :: val
187 REAL(rn),
PARAMETER :: zero = 0._rn
189 INTEGER ::
dp, ip, ipow, sp
201 DO ip = 1, comp(i)%ByteCodeSize
202 SELECT CASE (comp(i)%ByteCode(ip))
204 CASE (cimmed); sp = sp + 1; comp(i)%Stack(sp) = comp(i)%Immed(
dp);
dp =
dp + 1
205 CASE (cneg); comp(i)%Stack(sp) = -comp(i)%Stack(sp)
206 CASE (cadd); comp(i)%Stack(sp - 1) = comp(i)%Stack(sp - 1) + comp(i)%Stack(sp); sp = sp - 1
207 CASE (csub); comp(i)%Stack(sp - 1) = comp(i)%Stack(sp - 1) - comp(i)%Stack(sp); sp = sp - 1
208 CASE (cmul); comp(i)%Stack(sp - 1) = comp(i)%Stack(sp - 1)*comp(i)%Stack(sp); sp = sp - 1
209 CASE (cdiv);
IF (comp(i)%Stack(sp) == 0._rn) then;
evalerrtype = 1; res = zero; return;
END IF
210 comp(i)%Stack(sp - 1) = comp(i)%Stack(sp - 1)/comp(i)%Stack(sp); sp = sp - 1
213 IF (comp(i)%Stack(sp - 1) < 0.0_rn)
THEN
214 ipow = floor(comp(i)%Stack(sp))
215 IF (mod(comp(i)%Stack(sp), real(ipow, kind=rn)) == 0.0_rn)
THEN
216 comp(i)%Stack(sp - 1) = comp(i)%Stack(sp - 1)**ipow
218 cpabort(
"Negative floating-point value raised to a real power!")
221 comp(i)%Stack(sp - 1) = comp(i)%Stack(sp - 1)**comp(i)%Stack(sp)
224 CASE (cabs); comp(i)%Stack(sp) = abs(comp(i)%Stack(sp))
225 CASE (cexp); comp(i)%Stack(sp) = exp(comp(i)%Stack(sp))
226 CASE (clog10);
IF (comp(i)%Stack(sp) <= 0._rn) then;
evalerrtype = 3; res = zero; return;
END IF
227 comp(i)%Stack(sp) = log10(comp(i)%Stack(sp))
228 CASE (clog);
IF (comp(i)%Stack(sp) <= 0._rn) then;
evalerrtype = 3; res = zero; return;
END IF
229 comp(i)%Stack(sp) = log(comp(i)%Stack(sp))
230 CASE (csqrt);
IF (comp(i)%Stack(sp) < 0._rn) then;
evalerrtype = 3; res = zero; return;
END IF
231 comp(i)%Stack(sp) = sqrt(comp(i)%Stack(sp))
232 CASE (csinh); comp(i)%Stack(sp) = sinh(comp(i)%Stack(sp))
233 CASE (ccosh); comp(i)%Stack(sp) = cosh(comp(i)%Stack(sp))
234 CASE (ctanh); comp(i)%Stack(sp) = tanh(comp(i)%Stack(sp))
235 CASE (csin); comp(i)%Stack(sp) = sin(comp(i)%Stack(sp))
236 CASE (ccos); comp(i)%Stack(sp) = cos(comp(i)%Stack(sp))
237 CASE (ctan); comp(i)%Stack(sp) = tan(comp(i)%Stack(sp))
238 CASE (casin);
IF ((comp(i)%Stack(sp) < -1._rn) .OR. (comp(i)%Stack(sp) > 1._rn))
THEN
240 comp(i)%Stack(sp) = asin(comp(i)%Stack(sp))
241 CASE (cacos);
IF ((comp(i)%Stack(sp) < -1._rn) .OR. (comp(i)%Stack(sp) > 1._rn))
THEN
243 comp(i)%Stack(sp) = acos(comp(i)%Stack(sp))
244 CASE (catan); comp(i)%Stack(sp) = atan(comp(i)%Stack(sp))
245 CASE (cerf); comp(i)%Stack(sp) = erf(comp(i)%Stack(sp))
246 CASE (cerfc); comp(i)%Stack(sp) = erfc(comp(i)%Stack(sp))
247 CASE DEFAULT; sp = sp + 1; comp(i)%Stack(sp) = val(comp(i)%ByteCode(ip) - varbegin + 1)
251 res = comp(i)%Stack(1)
260 SUBROUTINE checksyntax(Func, FuncStr, Var)
264 CHARACTER(LEN=*),
INTENT(in) :: func, funcstr
265 CHARACTER(LEN=*),
DIMENSION(:),
INTENT(in) :: var
267 INTEGER :: ib, in, j, lfunc, parcnt
268 CHARACTER(LEN=1) :: c
281 lfunc = len_trim(func)
283 IF (j > lfunc)
CALL parseerrmsg(j, funcstr)
288 IF (c ==
'-' .OR. c ==
'+')
THEN
290 IF (j > lfunc)
CALL parseerrmsg(j, funcstr,
'Missing operand')
292 IF (any(c == ops))
CALL parseerrmsg(j, funcstr,
'Multiple operators')
294 n = mathfunctionindex(func(j:))
296 j = j + len_trim(funcs(n))
297 IF (j > lfunc)
CALL parseerrmsg(j, funcstr,
'Missing function argument')
299 IF (c /=
'(')
CALL parseerrmsg(j, funcstr,
'Missing opening parenthesis')
306 IF (scan(c,
'0123456789.') > 0)
THEN
307 r = realnum(func(j:), ib, in, err)
308 IF (err)
CALL parseerrmsg(j, funcstr,
'Invalid number format: '//func(j + ib - 1:j + in - 2))
313 n = variableindex(func(j:), var, ib, in)
314 IF (n == 0)
CALL parseerrmsg(j, funcstr,
'Invalid element: '//func(j + ib - 1:j + in - 2))
321 IF (parcnt < 0)
CALL parseerrmsg(j, funcstr,
'Mismatched parenthesis')
322 IF (func(j - 1:j - 1) ==
'(')
CALL parseerrmsg(j - 1, funcstr,
'Empty parentheses')
331 IF (any(c == ops))
THEN
332 IF (j + 1 > lfunc)
CALL parseerrmsg(j, funcstr)
333 IF (any(func(j + 1:j + 1) == ops))
CALL parseerrmsg(j + 1, funcstr,
'Multiple operators')
335 CALL parseerrmsg(j, funcstr,
'Missing operator')
343 IF (parcnt > 0)
CALL parseerrmsg(j, funcstr,
'Missing )')
344 END SUBROUTINE checksyntax
350 FUNCTION evalerrmsg()
RESULT(msg)
354 CHARACTER(LEN=*),
DIMENSION(4),
PARAMETER :: m = (/
'Division by zero ', &
355 'Argument of SQRT negative ',
'Argument of LOG negative ', &
356 'Argument of ASIN or ACOS illegal'/)
357 CHARACTER(LEN=LEN(m)) :: msg
366 END FUNCTION evalerrmsg
374 SUBROUTINE parseerrmsg(j, FuncStr, Msg)
378 INTEGER,
INTENT(in) :: j
379 CHARACTER(LEN=*),
INTENT(in) :: funcstr
380 CHARACTER(LEN=*),
INTENT(in),
OPTIONAL :: msg
382 CHARACTER(LEN=default_string_length) :: message
387 IF (
PRESENT(msg))
THEN
388 WRITE (unit=message, fmt=
"(A)")
"Syntax error in function string: "//msg
390 WRITE (unit=message, fmt=
"(A)")
"Syntax error in function string"
392 WRITE (*,
'(/,T2,A)') trim(funcstr)
393 IF ((j > lbound(ipos, dim=1)) .AND. (j <= ubound(ipos, dim=1)))
THEN
394 WRITE (*,
'(A)') repeat(
" ", ipos(j))//
"?"
396 WRITE (*,
'(A)') repeat(
" ",
SIZE(ipos) + 1)//
"?"
398 cpabort(trim(message))
400 END SUBROUTINE parseerrmsg
407 FUNCTION operatorindex(c)
RESULT(n)
411 CHARACTER(LEN=1),
INTENT(in) :: c
420 IF (c == ops(j))
THEN
425 END FUNCTION operatorindex
432 FUNCTION mathfunctionindex(str)
RESULT(n)
436 CHARACTER(LEN=*),
INTENT(in) :: str
439 CHARACTER(LEN=LEN(Funcs)) :: fun
447 k = min(len_trim(funcs(j)), len(str))
448 CALL lowcase(str(1:k), fun)
449 IF (fun == funcs(j))
THEN
454 END FUNCTION mathfunctionindex
464 FUNCTION variableindex(str, Var, ibegin, inext)
RESULT(n)
468 CHARACTER(LEN=*),
INTENT(in) :: str
469 CHARACTER(LEN=*),
DIMENSION(:),
INTENT(in) :: var
470 INTEGER,
INTENT(out),
OPTIONAL :: ibegin, inext
473 INTEGER :: ib, in, j, lstr
486 IF (str(ib:ib) /=
' ')
EXIT
489 IF (scan(str(in:in),
'+-*/^) ') > 0)
EXIT
492 IF (str(ib:in - 1) == var(j))
THEN
498 IF (
PRESENT(ibegin)) ibegin = ib
499 IF (
PRESENT(inext)) inext = in
500 END FUNCTION variableindex
506 SUBROUTINE removespaces(str)
510 CHARACTER(LEN=*),
INTENT(inout) :: str
517 ipos(:) = (/(k, k=1, lstr)/)
519 DO WHILE (str(k:lstr) /=
' ')
520 IF (str(k:k) ==
' ')
THEN
521 str(k:lstr) = str(k + 1:lstr)//
' '
522 ipos(k:lstr) = (/ipos(k + 1:lstr), 0/)
527 END SUBROUTINE removespaces
535 SUBROUTINE replace(ca, cb, str)
539 CHARACTER(LEN=*),
INTENT(in) :: ca
540 CHARACTER(LEN=LEN(ca)),
INTENT(in) :: cb
541 CHARACTER(LEN=*),
INTENT(inout) :: str
549 DO j = 1, len_trim(str) - lca + 1
550 IF (str(j:j + lca - 1) == ca) str(j:j + lca - 1) = cb
552 END SUBROUTINE replace
560 SUBROUTINE compile(i, F, Var)
564 INTEGER,
INTENT(in) :: i
565 CHARACTER(LEN=*),
INTENT(in) :: f
566 CHARACTER(LEN=*),
DIMENSION(:),
INTENT(in) :: var
573 IF (
ASSOCIATED(comp(i)%ByteCode))
DEALLOCATE (comp(i)%ByteCode, &
576 comp(i)%ByteCodeSize = 0
577 comp(i)%ImmedSize = 0
578 comp(i)%StackSize = 0
580 CALL compilesubstr(i, f, 1, len_trim(f), var)
581 ALLOCATE (comp(i)%ByteCode(comp(i)%ByteCodeSize), &
582 comp(i)%Immed(comp(i)%ImmedSize), &
583 comp(i)%Stack(comp(i)%StackSize))
584 comp(i)%ByteCodeSize = 0
585 comp(i)%ImmedSize = 0
586 comp(i)%StackSize = 0
588 CALL compilesubstr(i, f, 1, len_trim(f), var)
590 END SUBROUTINE compile
597 SUBROUTINE addcompiledbyte(i, b)
601 INTEGER,
INTENT(in) :: i
602 INTEGER(is),
INTENT(in) :: b
608 comp(i)%ByteCodeSize = comp(i)%ByteCodeSize + 1
609 IF (
ASSOCIATED(comp(i)%ByteCode)) comp(i)%ByteCode(comp(i)%ByteCodeSize) = b
610 END SUBROUTINE addcompiledbyte
619 FUNCTION mathitemindex(i, F, Var)
RESULT(n)
623 INTEGER,
INTENT(in) :: i
624 CHARACTER(LEN=*),
INTENT(in) :: f
625 CHARACTER(LEN=*),
DIMENSION(:),
INTENT(in) :: var
635 IF (scan(f(1:1),
'0123456789.') > 0)
THEN
636 comp(i)%ImmedSize = comp(i)%ImmedSize + 1
637 IF (
ASSOCIATED(comp(i)%Immed)) comp(i)%Immed(comp(i)%ImmedSize) = realnum(f)
640 n = variableindex(f, var)
641 IF (n > 0) n = varbegin + n - 1_is
643 END FUNCTION mathitemindex
652 FUNCTION completelyenclosed(F, b, e)
RESULT(res)
656 CHARACTER(LEN=*),
INTENT(in) :: f
657 INTEGER,
INTENT(in) :: b, e
667 IF (f(b:b) ==
'(' .AND. f(e:e) ==
')')
THEN
670 IF (f(j:j) ==
'(')
THEN
672 ELSEIF (f(j:j) ==
')')
THEN
677 IF (k == 0) res = .true.
679 END FUNCTION completelyenclosed
689 RECURSIVE SUBROUTINE compilesubstr(i, F, b, e, Var)
693 INTEGER,
INTENT(in) :: i
694 CHARACTER(LEN=*),
INTENT(in) :: f
695 INTEGER,
INTENT(in) :: b, e
696 CHARACTER(LEN=*),
DIMENSION(:),
INTENT(in) :: var
698 CHARACTER(LEN=*),
PARAMETER :: &
699 calpha =
'abcdefghijklmnopqrstuvwxyz'//
'ABCDEFGHIJKLMNOPQRSTUVWXYZ'
701 INTEGER :: b2, io, j, k
712 IF (f(b:b) ==
'+')
THEN
714 CALL compilesubstr(i, f, b + 1, e, var)
716 ELSEIF (completelyenclosed(f, b, e))
THEN
718 CALL compilesubstr(i, f, b + 1, e - 1, var)
720 ELSEIF (scan(f(b:b), calpha) > 0)
THEN
721 n = mathfunctionindex(f(b:e))
723 b2 = b + index(f(b:e),
'(') - 1
724 IF (completelyenclosed(f, b2, e))
THEN
726 CALL compilesubstr(i, f, b2 + 1, e - 1, var)
727 CALL addcompiledbyte(i, n)
731 ELSEIF (f(b:b) ==
'-')
THEN
732 IF (completelyenclosed(f, b + 1, e))
THEN
734 CALL compilesubstr(i, f, b + 2, e - 1, var)
735 CALL addcompiledbyte(i, cneg)
737 ELSEIF (scan(f(b + 1:b + 1), calpha) > 0)
THEN
738 n = mathfunctionindex(f(b + 1:e))
740 b2 = b + index(f(b + 1:e),
'(')
741 IF (completelyenclosed(f, b2, e))
THEN
743 CALL compilesubstr(i, f, b2 + 1, e - 1, var)
744 CALL addcompiledbyte(i, n)
745 CALL addcompiledbyte(i, cneg)
757 IF (f(j:j) ==
')')
THEN
759 ELSEIF (f(j:j) ==
'(')
THEN
762 IF (k == 0 .AND. f(j:j) == ops(io) .AND. isbinaryop(j, f))
THEN
763 IF (any(f(j:j) == ops(cmul:cpow)) .AND. f(b:b) ==
'-')
THEN
765 CALL compilesubstr(i, f, b + 1, e, var)
766 CALL addcompiledbyte(i, cneg)
770 CALL compilesubstr(i, f, b, j - 1, var)
771 CALL compilesubstr(i, f, j + 1, e, var)
772 CALL addcompiledbyte(i, operatorindex(ops(io)))
773 comp(i)%StackPtr = comp(i)%StackPtr - 1
783 IF (f(b:b) ==
'-') b2 = b2 + 1
784 n = mathitemindex(i, f(b2:e), var)
786 CALL addcompiledbyte(i, n)
787 comp(i)%StackPtr = comp(i)%StackPtr + 1
788 IF (comp(i)%StackPtr > comp(i)%StackSize) comp(i)%StackSize = comp(i)%StackSize + 1
789 IF (b2 > b)
CALL addcompiledbyte(i, cneg)
790 END SUBROUTINE compilesubstr
798 FUNCTION isbinaryop(j, F)
RESULT(res)
804 INTEGER,
INTENT(in) :: j
805 CHARACTER(LEN=*),
INTENT(in) :: f
809 LOGICAL :: dflag, pflag
817 IF (f(j:j) ==
'+' .OR. f(j:j) ==
'-')
THEN
820 ELSEIF (scan(f(j - 1:j - 1),
'+-*/^(') > 0)
THEN
822 ELSEIF (scan(f(j + 1:j + 1),
'0123456789') > 0 .AND. &
823 scan(f(j - 1:j - 1),
'eEdD') > 0)
THEN
824 dflag = .false.; pflag = .false.
828 IF (scan(f(k:k),
'0123456789') > 0)
THEN
830 ELSEIF (f(k:k) ==
'.')
THEN
840 IF (dflag .AND. (k == 1 .OR. scan(f(k:k),
'+-*/^(') > 0)) res = .false.
843 END FUNCTION isbinaryop
853 FUNCTION realnum(str, ibegin, inext, error)
RESULT(res)
857 CHARACTER(LEN=*),
INTENT(in) :: str
858 INTEGER,
INTENT(out),
OPTIONAL :: ibegin, inext
859 LOGICAL,
INTENT(out),
OPTIONAL :: error
862 INTEGER :: ib, in, istat
863 LOGICAL :: bflag, dinexp, dinman, eflag, err, &
881 bflag = .true.; inman = .false.; pflag = .false.; eflag = .false.; inexp = .false.
882 dinman = .false.; dinexp = .false.
885 DO WHILE (in <= len_trim(str))
886 SELECT CASE (str(in:in))
889 IF (inman .OR. eflag .OR. inexp)
EXIT
892 inman = .true.; bflag = .false.
894 inexp = .true.; eflag = .false.
900 inman = .true.; bflag = .false.
902 inexp = .true.; eflag = .false.
904 IF (inman) dinman = .true.
905 IF (inexp) dinexp = .true.
909 inman = .true.; bflag = .false.
910 ELSEIF (inman .AND. .NOT. pflag)
THEN
915 CASE (
'e',
'E',
'd',
'D')
917 eflag = .true.; inman = .false.
926 err = (ib > in - 1) .OR. (.NOT. dinman) .OR. ((eflag .OR. inexp) .AND. .NOT. dinexp)
930 READ (str(ib:in - 1), *, iostat=istat) res
933 IF (
PRESENT(ibegin)) ibegin = ib
934 IF (
PRESENT(inext)) inext = in
935 IF (
PRESENT(error)) error = err
943 SUBROUTINE lowcase(str1, str2)
947 CHARACTER(LEN=*),
INTENT(in) :: str1
948 CHARACTER(LEN=*),
INTENT(out) :: str2
950 CHARACTER(LEN=*),
PARAMETER :: lc =
'abcdefghijklmnopqrstuvwxyz', &
951 uc =
'ABCDEFGHIJKLMNOPQRSTUVWXYZ'
958 DO j = 1, len_trim(str1)
959 k = index(uc, str1(j:j))
960 IF (k > 0) str2(j:j) = lc(k:k)
962 END SUBROUTINE lowcase
975 FUNCTION evalfd(id_fun, ipar, vals, h, err)
RESULT(derivative)
976 INTEGER,
INTENT(IN) :: id_fun, ipar
977 REAL(kind=rn),
DIMENSION(:),
INTENT(INOUT) :: vals
978 REAL(kind=rn),
INTENT(IN) :: h
979 REAL(kind=rn),
INTENT(OUT) :: err
980 REAL(kind=rn) :: derivative
982 INTEGER,
PARAMETER :: ntab = 10
983 REAL(kind=rn),
PARAMETER :: big_error = 1.0e30_rn, con = 1.4_rn, &
984 con2 = con*con, safe = 2.0_rn
987 REAL(kind=rn) :: a(ntab, ntab), errt,
fac, funcm, funcp, &
990 derivative = huge(0.0_rn)
994 vals(ipar) = xval + hh
995 funcp =
evalf(id_fun, vals)
996 vals(ipar) = xval - hh
997 funcm =
evalf(id_fun, vals)
998 a(1, 1) = (funcp - funcm)/(2.0_rn*hh)
1002 vals(ipar) = xval + hh
1003 funcp =
evalf(id_fun, vals)
1004 vals(ipar) = xval - hh
1005 funcm =
evalf(id_fun, vals)
1006 a(1, i) = (funcp - funcm)/(2.0_rn*hh)
1009 a(j, i) = (a(j - 1, i)*
fac - a(j - 1, i - 1))/(
fac - 1.0_rn)
1011 errt = max(abs(a(j, i) - a(j - 1, i)), abs(a(j, i) - a(j - 1, i - 1)))
1012 IF (errt .LE. err)
THEN
1014 derivative = a(j, i)
1017 IF (abs(a(i, i) - a(i - 1, i - 1)) .GE. safe*err)
RETURN
1020 cpabort(
"DX provided equals zero!")
static GRID_HOST_DEVICE double fac(const int i)
Factorial function, e.g. fac(5) = 5! = 120.
This public domain function parser module is intended for applications where a set of mathematical ex...
real(rn) function, public evalf(i, Val)
...
integer, public evalerrtype
real(kind=rn) function, public evalfd(id_fun, ipar, vals, h, err)
Evaluates derivatives.
subroutine, public finalizef()
...
subroutine, public initf(n)
...
subroutine, public parsef(i, FuncStr, Var)
Parse ith function string FuncStr and compile it into bytecode.
Defines the basic variable types.
integer, parameter, public dp
integer, parameter, public default_string_length