Fortran Coder

标题: 关于FORTRAN的EXP函数 [打印本页]

作者: Kieran 时间: 2021-5-8 09:18
标题: 关于FORTRAN的EXP函数
大家好，

我想咨询下关于精度的问题。

1.如果我想把变量a设置成双精度浮点数，是否可以如下设置？
INTEGER, PARAMETER :: dp=SELECTED_REAL_KIND(15,14)
REAL (KIND=dp) :: a
我是想问，15位有效小数位，取值在10的-14次方到10的14次方的范围。这样设置是否可以，确定变量a为双精度浮点数呢？

2.如果想把整型变量也设置成双精度的话，SELECTED_REAL_KIND中的两个参数要如何设置呢？

3.我按上面的方式设置了双精度的浮点数和整数，并将它们放在了EXP函数里运算，但这个函数却返回Infinity类型的数据。请问要如何解决这个问题呢？这是不是由双精度的浮点数和整数设置不合理导致的呢？

谢谢大家啦。

作者: kyra 时间: 2021-5-8 10:54
1. 可以确定。
2. 整形变量是无误差的，不存在双精度，只存在长度（可以用 SELECTED_INT_KIND）
3. 上代码。

作者: Kieran 时间: 2021-5-8 11:43
谢谢你的回复。

下面是我出问题的代码部分。

[Fortran] 纯文本查看 复制代码

INTEGER, PARAMETER                      :: dp = SELECTED_REAL_KIND(15,14)

INTEGER                                                  :: i, j

REAL (KIND=dp), ALLOCATABLE   :: ev(:,:)

REAL (KIND=dp)                                   :: fermi

REAL (KIND=dp)                                   :: kb

REAL (KIND=dp)                                   :: te

ei = (ev(i,j) - fermi) * (1.0d0 /(EXP((ev(i,j) - fermi)/kb / te) + 1.0d0))&

         + kb * te * DLOG(1.0d0 + EXP(0.0d0 - (ev(i,j) - fermi) /kb /te))

write (unit=*, fmt=*) 'kb', kb

write (unit=*, fmt=*) 'fermi', fermi

write (unit=*, fmt=*) 'te', te

write (unit=*, fmt=*) 'ev', ev(i,j)

write (unit=*, fmt=*) '(ev(i,j) - fermi)', (ev(i,j) - fermi)

write (unit=*, fmt=*) 'f_function', (1.0d0 /(EXP((ev(i,j) - fermi)/kb / te) + 1.0d0))

write (unit=*, fmt=*) 'EXP(0.0d0 - (ev(i,j) - fermi) /kb /te)', EXP(0.0d0 - (ev(i,j) - fermi) /kb /te)

write (unit=*, fmt=*) 'DLOG_function', DLOG(1.0d0 + EXP(0.0d0 - (ev(i,j) - fermi) /kb /te))

write (unit=*, fmt=*) 'ei=', ei

其中，ev(i,j)是从主程序传递进来的一个二维双精度数组。

编译并运行程序后，出现下面的结果。
kb  8.617333262145000E-005
fermi -0.531700000000000
te 30.0000000000000
ev  -20.4013773387809
(ev(i,j) - fermi)  -17.9840108474454
f_function 1.00000000000000
EXP(0.0d0 - (ev(i,j) - fermi) /kb /te)             Infinity
DLOG_function             Infinity
ei=             Infinity

因为EXP函数的结果是Infinity类型的数据，导致后面的DLOG函数结果和ei变量也都变成了Infinity类型的数据。

请问，1. 我对于双精度浮点数据类型的设定是否有误呢？
2. 要如何调整代码，让EXP函数的输出结果合理呢？

谢谢啦。

作者: kyra 时间: 2021-5-8 14:19
EXP((ev(i,j) - fermi)/kb / te)
=exp(-17.9840108474454/8.617333262145000E-005/30)
=exp(-6956.52573015480)
=6.5957E-3022

双精度可正常表达的精度范围是 -1.79E+308 到 +1.79E+308，最小精度 2.2250738585072014e-308
即，小于 2.2250738585072014e-308 的数字认为是 0

分母为0，则 1.0d0 /(EXP((ev(i,j) - fermi)/kb / te) 为无穷大。

作者: kyra 时间: 2021-5-8 14:31
试试用双双精度。
dp = SELECTED_REAL_KIND(30,3000)

代码中的 1.d0 改成 1.0_dp ， 0.0d0 改成 0.0_dp （其他常量类似）
DLOG 改成 Log

作者: Kieran 时间: 2021-5-11 19:00

kyra 发表于 2021-5-8 14:31
试试用双双精度。
dp = SELECTED_REAL_KIND(30,3000)

非常感谢你的回复。

我还想请教一个很基础的问题。selected_real_kind（[p，r]）函数中的p和r分别指代的是小数精度p和十进制指数范围r。既然FORTRAN中的双精度浮点数范围是从2.2250738585072008D-308到4.94065645841246544D-324。那是不是说，双精度的浮点数据类型的设定应该是selected_real_kind（16，324），而非selected_real_kind（15，14）呢？我是指最大范围地包括所有的双精度浮点数据。

谢谢啦。

作者: vvt 时间: 2021-5-12 08:56
编译器为你提供几种类型的 real，比如单精度，双精度，四精度（也叫双双精度，有的编译器不提供四精度）
分别用不同的 kind 值来代表，对大部分编译器，上面的三种real类型的kind值是4，8，16（但不一定）

而
selected_real_kind(p,r)
函数的目的，就是请编译器为你选择合适的real类型。
这个函数，会返回至少满足 p 位十进制有效位数，且能表达最大 r 位十进制的 real 类型中“最小”的那一个。

比如 p=5,r=15那么单精度能满足，双精度也能满足，四精度仍然满足。那么编译器返回单精度的kind值（通常是4）
比如 p=9,r=15那么单精度不能满足，双精度能满足，四精度也满足。那么编译器返回双精度的kind值（通常是8）
而 p=10,r=100，也是一样，单精度不满足，双精度四精度满足。所以也返回双精度的kind值（通常是8）

回到你的问题，
selected_real_kind（16，324），因为双精度只有15位有效数字，所以不能满足，编译器会返回四精度的kind值（通常是16）
selected_real_kind（15，14），双精度可以满足，因此返回双精度的kind值（通常是8）

作者: Kieran 时间: 2021-5-18 05:14

vvt 发表于 2021-5-12 08:56
编译器为你提供几种类型的 real，比如单精度，双精度，四精度（也叫双双精度，有的编译器不提供四精度）
...

谢谢你的指点。

也就是说，SELECTED_REAL_KIND函数是把，由p和r指定的浮点数，归类成单精度，双精度或者四精度的类型；而其归类的方式，主要是根据p确定的有效数字，来确定的。至于r这个参数的作用不大。因为（p=16,r=324）和(p=15,r=14)的区别主要是由p这个变量来确定的，r的作用似乎可以忽略不计。

我是否可以这样来理解呢？

作者: vvt 时间: 2021-5-18 08:30
不是的，p 和 r 都起作用。
SELECTED_REAL_KIND 的作用是：“选择一个精度类型，使其至少满足 p，也同时满足 r。如果同时满足的有多个，则选择占内存更少的那个”

比如
p=5 , r = 30 。返回单精度的kind值
但
p=5 , r = 80 。则需要返回双精度的kind值。
（因为单精度有效数字是6-7，双精度大约15。都能满足p。但单精度 r 只能到38，要 r=80，必须至少双精度）

作者: Kieran 时间: 2021-5-19 08:51

vvt 发表于 2021-5-18 08:30
不是的，p 和 r 都起作用。
SELECTED_REAL_KIND 的作用是：“选择一个精度类型，使其至少满足 p，也同时满 ...

非常感谢你的解释，我现在懂了。

欢迎光临 Fortran Coder (http://bbs.fcode.cn/)

Powered by Discuz! X3.2