科学计算，特别是有限元要大量使用数组，说白了有限元法就是求解一个大型甚至是超大型的线性方程组，转换到编程语言里就是用数组来描述了。但上下界都是定值的固定数组在实际编程中的用途真是太狭窄了，真正有用的当数动态数组和可调数组了。在c中，动态数组是用内存分配函数malloc来实现的；fortran则是用虚实结合的参数传递实现的。下面就用一个例题来说明：求解常用对数底数E=2.71828……的任意精度值，算法很简单，把E的表达式展开为泰勒级数，也就是一个多项式，程序模拟手算的方式以求解任意精度，结果用数组保存；这样以来数组的大小就是不确定的了，因此便用到动态数组与可调数组。
c程序：
#include
#include
int *count( int x)
{
int *s,*term,i,j,k,l,m,n;
s=(int*)malloc((x+8)*sizeof(int));
term=(int*)malloc((x+8)*sizeof(int));
s[1]=1;term[1]=1;
for(k=2;k<=x+5;k++)
{ s[k]=0;
term[k]=0;
}
for(i=1;i<=x+1;i++)
{ m=0;
for(k=1;k<=x+5;k++)
{ n=10*m+term[k];
term[k]=n/i;
m=n%i;
s[k]=s[k]+term[k];
}
}
for(k=1;k<=x+5;k++)
{ j=x+6-k;
n=s[j]+l;
s[j]=n%10;
l=n/10;
}
return (s);
free(s);free(term);
}
void main()
{ int m,i,*c;
printf("n Input the accuracy:");
scanf("%d",&m);
c=count(m);
printf("nE=%d.",c[1]);
for(i=2;i<=m+1;i++)printf("%d",c);
}
在count函数中，由malloc函数动态分配了数组s[]和term[]的大小，而它们的大小则是由主函数main中输入的精度参数m决定的。在用完后可以用free函数释放掉，而不占用内存。其中s是存储最终结果的，term是存储临时结果的。
下面提供几个fortran程序，以便对照：
fortran77程序：
program main
integer s(10000),term(10000),x
write(*,*)"Input the accuracy:"
read(*,*)x
call count(s,term,x)
write(*,30)(s(i),i=1,x+1)
30 format(1x,' E=',i1,'.',50i1,9(/5x,50i1))
end
subroutine count(s,term,x)
integer s(x+4),term(x+4),x
s(1)=1
term(1)=1
do 1 k=2,x+4
s(k)=0
term(k)=0
1 continue
do 100 i=1,x
m=0
do 10 k=1,x+4
n=10*m+term(k)
term(k)=n/i
m=mod(n,i)
s(k)=s(k)+term(k)
10 continue
100 continue
l=0
do 20 k=1,x+4
j=x+5-k
n=s(j)+l
s(j)=mod(n,10)
l=n/10
20 continue
end
从上面的程序可以看到，由于fortran77不支持动态数组，只能在主程序中预先定义两个数组(当然是大一点的好)，然后把它们传递到count子程序中，用虚实结合的方法完成计算。这样一来，程序就显得“笨拙”一些，定义了两个大型数组，用完后却不能释放，白白浪费了内存。而且上面的程序还存在着一个bug，当输入的精度参数x大于10000的话，程序可能会出现莫名其妙的错误，甚至崩溃。这是由于语言本身的限制造成的，当然我们可以加上一个条件判断语句，当条件不满足时跳出来。刚开始的时候，我还想到用下面的“方法”来计算：
program main
integer x
read x
call count(x)
.
.
.
end
subroutine count(x)
integer x,s(x+4),term(x+4)
.
.
.
end
乍看起来好像没什么错误，数组s和term是在子程序中定义的，但它们的边界中却含有变量，这样的数组称为“自动数组”，而fortran77是不支持自动数组的，所以我用上面的方法来计算的结果就是——编译根本通不过！这个错误困扰我好长一段时间，毕竟它太隐蔽了，直到我查了fortran90的资料后才明白是怎么回事。fortran90明确的规定支持自动数组和动态数组，上面这个fortran77不能编译的程序在fortran90中可
以畅行无阻，只要把空缺的部分补上就可以了，我就不再重复写了。
下面再看一个具有fortran90特色的版本：
program main
integer x
write(*,*)"Input the accuracy:"
read(*,*)x
call count(x)
end
subroutine count(x)
integer x
integer,allocatable::s(:),term(:)
allocate(s(x+4))
allocate(term(x+4))
s(1)=1
term(1)=1
do k=2,x+4
s(k)=0
term(k)=0
enddo
do i=1,x
m=0
do k=1,x+4
n=10*m+term(k)
term(k)=n/i
m=mod(n,i)
s(k)=s(k)+term(k)
enddo
enddo
l=0
do k=1,x+4
j=x+5-k
n=s(j)+l
s(j)=mod(n,10)
l=n/10
enddo
write(*,30)(s(i),i=1,x+1)
30 format(1x,' E=',i1,'.',50i1,9(/5x,50i1))
deallocate(s)
deallocate(term)
end
这个fortran90版本的程序和fortran77版本有个小区别：把输出语句write(*,30)放到count子程序里面了，这是因为数组s在用完后就释放掉了，在主程序中再出现s的时候就是无定义的了。而c语言版本与fortran90又有不同：虽然c语言中也把s释放了，但由于count()函数把s的地址传回了主函数，所以在main()中依然可以使用它的内容。
二、函数
函数这个范围太大，这里只简单的讨论“以函数作为函数的参数”这个问题。以函数作为函数的参数，这是在科学计算中经常用到的，比如说求某个函数的数值积分，就需要把这个函数作为参数传递到子程序中去。例题：求函数4/(1+x*x)在区间[0，1]上的积分，使用数值分析中常用的变步长梯形求积法。
c语言版本：
#include
#include "math.h"
double ffts(double a,double b,double eps,double (*f)())
{ int n,k;
double fa,fb,h,t1,p,s,x,t;
fa=(*f)(a); fb=(*f)(b);
n=1; h=b-a;
t1=h*(fa+fb)*.5;
p=eps+1.0;
while (p>=eps)
{ s=0.0;
for (k=0;k<=n-1;k++)
{ x=a+(k+0.5)*h;
s=s+(*f)(x);
}
t=(t1+h*s)*.5;
p=fabs(t1-t);
t1=t; n=n+n; h=h*.5;
}
return(t);
}
double fftsf(double x)
{ double y;
y=4.0/(1.0+x*x);
return(y);
}
void main()
{ double a,b,eps,t,fftsf(double);
a=0.0; b=1.0; eps=1e-13;
t=ffts(a,b,eps,fftsf);
printf("n");
printf("t=%.13en",t);
printf("n");
}
c语言中用指向函数的指针解决了这个问题，即用(*f)(a)，(*f)(b)引用原函数的值，使函数具有“通用性”。
fortran程序：
function ffts(a,b,eps,f)
real*8 ffts,a,b,eps,f
integer n,k
real*8 fa,fb,h,t1,s,x,t
fa=f(a)
fb=f(b)
n=1
h=b-a
t1=h*(fa+fb)*0.5
p=eps+1.0
do while (p.ge.eps)
s=0.0
do k=0,n-1
x=a+(k+0.5)*h
s=s+f(x)
end do
t=(t1+h*s)*0.5
p=abs(t1-t)
t1=t
n=n+n
h=h*0.5
end do
ffts=t
return
end
c==========================
function f(x)
real*8 f,x
f=4.0/(1.0+x*x)
return
end
c============================
EXTERNAL F
real*8 a1,b1,eps1,f,ffts
a1=0.0
b1=1.0
eps1=1e-13
write(*,10)ffts(a1,b1,eps1,f)
10 FORMAT(1X,'S=',g25.17)
end
fortran语言中，子程序的功能真的是很强大，不论是什么样的参数——变量、数组、甚至是函数，它都能以相同的方式“接受”，所以fortran的程序比c语言的要简洁一些。fortran90的函数方面真正有价值的改进我觉得不是很多，比较实用的是interface语句，对于上面的算法，90版本的程序还不如77版本的简练，因此我就不再写了。
总结：感觉fortran77真的是“老了”，用它编程时会受到太多的限制，比如不支持动态数组、不支持递归、函数的返回结果只能是变量等等，而这些恰恰又是科学计算中最有用的东西；不过我们应该看到的是：由于fortran77统治数值计算界几十年了，前人给我们积累了大量的经验，这些原代码经历了无数人的修改，那些所谓的“限制”已经被巧妙的“回避”了，这些程序变得非常稳定，有许多至今还在稳定的运行着；细读这些程序，我时常会感慨，前人在软件硬件如此受限的条件下还能写出这样高质量的程序，身处今天的我们还能说什么呢？尽管我们有pentium4、core2，我们还有Compaq
Visual Fortran，Power Station
Fortran，我们能分配几百M，甚至1G的动态数组，但我们仍要吸取前人在这些程序中留下的精髓，它不论在何时都是一笔宝贵的财富。fortran90带给我们的是一种全新的体验，用fortran90写程序那叫一个“畅快”。畅快之后应该看到的是，fortran90这些崭新的内容有多少是我们能用到的。对我来说，要算是指针、动态数组和递归了，当然也许是因为专业的原因有些功能我用不到吧。从上面几个小程序可以看到，有时候fortran90的程序比fortran77会冗长一些，所以我觉得写程序的时候，能用77写的就尽量用77写，或者说是“写77风格的90程序”(当然不是说写那种打孔编码的书写格式)，毕竟这样会简练些。至于C，我不想做过多的评论，c语言是一个博大精深的语言，也并非我等小辈所能评论的。作为和fortran77同时代的语言，c却有着无与伦比的描述能力，fortran77和fortran90能做到的，c也都能做到，虽然有一些烦琐。在数值计算领域，c语言功不可没。因此，只要你的精力还有一点点宽裕，还是把c也学会吧，“技不压人”嘛。
1562
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
