Fortran 高精度定义

http://epsilong.blog.163.com/blog/static/8569511820095143752176/

2009-06-14 15:07:52|  分类： 程序编程 |举报 |字号 订阅

Fortran跟C一样，也分整型(INTEGER)，实型(REAL)，双精度(REAL*8或REAL(8)或DOUBLE
PRECISION)，这些在科学计算中还是比较重要的，以实型数为例：
一般REAL等价于REAL*4或REAL(4),是单精度的；
而双精度在F77中表示为DOUBLE PRECISION，在F90中可以表示为REAL*8或REAL(8)，在高精
度计算中，双精度的变量是很有必要的，对于一般实数可以表示为小数形式或指数形式，
而双精度都表示成指数形式，但指数E要改成D，如：
REAL:100.0或1e2,双精度下就得表示成1D2
由于Fortran中不需要对每个变量都进行声明，所以有时候会在每个程序或子程序开头做个
说明，如下：
IMPLICIT DOUBLE PRECISION(A-H,O-Z)
代表以A-H以及O-Z字母开头的变量默认（在不声明的情况下）是双精度的，否则则是整型
的，如下：
******************************************************************************
PROGRAM D
IMPLICIT DOUBLE PRECISION(A-H,O-Z)
J1=1D-2
J2=-0.5D-1
x=J1+J2
print *,x
end
******************************************************************************
PROGRAM E
implicit double precision (A-I,O-Z)
double precision a,i,e1,e2
data j2 /0.87450547081842D-3/
data j3 /-0.11886910646016D-4/
data j5 /-0.17242068505339D-5/
data j7 /0.10566966079622D-6/
write(*,*) "please input a"
read(*,*) a
write(*,*) "please input i"
read(*,*) i
e1=(j3*sin(i)/(2*a*j2)-5*j5*sin(i)*(1-7*sin(i)**2/2+21*sin(i)**4/8)&
&/(2*a**3*(2-5*sin(i)**2/2))+35*j7*sin(i)*(1-27*sin(i)**2/4+99&
&*sin(i)**4/8-429*sin(i)**6/64)/(3*a**5*(2-5*sin(i)**2/2)))
e2=-(j3*sin(i)/(2*a*j2)-5*j5*sin(i)*(1-7*sin(i)**2/2+21*sin(i)**4/8)&
&/(2*a**3*(2-5*sin(i)**2/2))+35*j7*sin(i)*(1-27*sin(i)**2/4+99&
&*sin(i)**4/8-429*sin(i)**6/64)/(3*a**5*(2-5*sin(i)**2/2)))
write(*,"(E9.2E3)") e1,e2
stop
end
******************************************************************************
第一个程序输出不是-0.4而是0.000000000000000E+000
第二个程序任意输入a、i，并未得到希望得到的结果，而是输出NAN和NAN，关于NAN这个错
误，有时候函数定义域不符合的时候，运行并不报错而是输出NAN，这个时候检查程序这些
地方是检查的重点，当然，会有其他情况，但我碰到的不多，只好就我所知跟大家交流一
下。
这两个程序都因为J开头的变量不属于默认双精度变量，而用双精度表示给它们赋值了，导
致结果跟预期不一致，在程序中把这些以J开头的变量用REAL*8声明一下，或把
implicit double precision (A-I,O-Z)改为：
implicit double precision (A-J,O-Z)，或把这个语句去掉
就可以得到预期的结果了。
对于数组，可以用DIMENSION定义，但需要注意的是，若在程序头未做声明（implicit
none）时，用DIMENSION定义数组时，当数组名首字母不属于(A-J,O-Z)里时，其值输出时
为整型，当然做了如下声明情况也会如此：（implicit double precision (A-I,O-Z)）
如下：
******************************************************************************
PROGRAM F
dimension m(2)
m(1)=1.5
m(2)=2.5
print *,m(1),m(2)
end
******************************************************************************
输出的结果是“1，2”而不是“1.500000,2.500000”
当把程序中m改为a时，输出“1.500000,2.500000”
所以，比较好的方法是尝试用REAL来定义数组（当然也可以用REAL*8）：
******************************************************************************
PROGRAM G
real m(2)
m(1)=1.5
m(2)=2.5
print *,m(1),m(2)
end
******************************************************************************
另外，要说的是，变量可以不定义而直接赋值，但会出现如上面PROGRAM D-E的问题，所以
建议大家在编程的时候对非整型变量声明一下，尽管麻烦，但不容易出错，有时候正是这
类错误会让初学者困扰好久。
定义变量时，经常会看到两种定义的写法：以REAL为例：
可以有
real m
和 real:: m
第一种方式不可以直接赋值，必须写成这样：
******************************************************************************
PROGRAM H
real m
m=1.0
print *,m
end
******************************************************************************
第二种则可以：
******************************************************************************
PROGRAM I
real:: m=1.0
print *,m
end