您现在的位置是:首页 >其他 >C#,数值计算的进化与发展——FORTRAN 77/80/95源程序 转C# 源程序的软件F2C#网站首页其他
C#,数值计算的进化与发展——FORTRAN 77/80/95源程序 转C# 源程序的软件F2C#
1 F2C#FORTRAN 77/80/95源程序 转C# 源程序的软件
1.1 F2C#起源
全世界科学计算领域在超过40年的时间里累积了巨大数量的FORTRAN源程序(尤其以FORTRAN77居多),实际上目前的许多大型科学软件还是以这些代码为基础的。众所周知的原因,FORTRAN在用户界面、图形等方面的缺陷不是可以很好解决的,即使应用了混合式编程也无法实现许多梦想的效果。C#其实是最适合进行科学计算的程序语言,它可以完美解决建模、计算、可视化的全过程,只不过多数程序员或者科研人员并不能领悟其中的奥秘。还有一个原因是大量原来运行稳定的程序,如果要翻译为C#,需要耗费大量的时间与经费。随着CPU速度的不断提升,内存和外存不断扩大,原来限制科学计算的硬件需求已经不是问题。
联高软件在三维图形系统、大级稀疏矩阵求解、图像处理系统、有限元分析软件、流体力学计算等都全面采用了C#作为编程语言,经过多年运行与计算测试,其性能与FORTRAN及C/C++相差无几。联高公司内部就采用FORTRAN77转C#的工具F2C#将大量的FORTRANS77源程序转换为C#源码。我们建议大家采用C#作为科学计算的编程语言。
为了普及应用C#,帮助您把累积多年的FORTRAN代码转换为C#源代码,这里提供翻译FORTRAN 77转为C#的翻译软件。F2C#可以翻译大多数标准FORTRAN77格式的语句、变量定义、循环等,尤其是可以解决变量定义、数组定义、IF语句和DO循环语句的翻译,节省大量的翻译时间。
1.2 F2C#特色功能
(1)F2C#是逐行翻译的,非常方便与FORTRAN 77源程序进行比对,以进行调试和修改;
(2)F2C#可以自动处理传入的参数变量、局部变量及其初始化(DATA);
(3)F2C#处理了数组的上标、下标问题;多维数组的问题;
(4)F2C#自动处理DO(WHILE,UNTIL)循环、条件语句及其跳转、GOTO (1,2,3,4) 分支语句;
(5)F2C#是可以处理复数Complex变量定义及其相关操作;
(6)我们不能承诺F2C#翻译100%准确;已经知道的未处理语句是WRITE&FORMAT及个别GOTO跳转;
1.3 软件特点
(1)以软件形式发布的F2C#处理速度更快;
(2)F2C#软件可以更好地处理烦人的GOTO跳转,未处理的跳转几乎很少了;
(3)可以一次性处理指定目录及其子目录下的所有FORTRAN 77程序;
(4)可以生成翻译结果的参考文档,便于详细了解翻译后应该如何做少量的修改;
1.4 生动的故事
通过算法的改进可解决语言差异带来的性能问题
不可否认FORTRAN的唯一好处是比C/C++,C#具有更快的计算速度;然而,通过算法的改进完全可以解决语言带来的性能问题。
这里有个生动的故事:有限元分析软件现在基于的算法有两种:一种是所谓的直接解法,一种是所谓的迭代解法。因为有限元软件处理的对象大部分是高元方程组,因此直接解法总能求得解,但速度就不能保证很快,虽然迭代法解题的速度很快,但不能保证所有的算法都是收敛的,因此传统的有限元分析软件大都采用直接解法。1982年,前苏联的三位数学物理博士致力于研究有限元分析的迭代算法的收敛性问题,因为,如果能保证迭代法总是收敛的,就可以大幅度提高解题速度。他们采用穷举法,分析迭代法中所有发散的算法,最后总结出500多种导致迭代法发散的原因并加以有效的对症下药,终于在11年后的1993年发明了举世瞩目的FFE(快速有限元法,Fast Finite Element )算法。FFE方法其实就是针对不同的迭代算法总能保证其收敛的改进的迭代法。
1.5 翻译实例
1.5.1 原始FORTRAN程序
SUBROUTINE EXCHNG( IREL, NSP, DS, EX, VX )
C *****************************************************************
C Finds local exchange energy density and potential
C Adapted by J.M.Soler from routine velect of Froyen's
C pseudopotential generation program. Madrid, Jan'97. Version 0.5.
C **** Input ******************************************************
C INTEGER IREL : relativistic-exchange switch (0=no, 1=yes)
C INTEGER NSP : spin-polarizations (1=>unpolarized, 2=>polarized)
C REAL*8 DS(NSP) : total (nsp=1) or spin (nsp=2) electron density
C **** Output *****************************************************
C REAL*8 EX : exchange energy density
C REAL*8 VX(NSP) : (spin-dependent) exchange potential
C **** Units ******************************************************
C Densities in electrons/Bohr**3
C Energies in Hartrees
C *****************************************************************
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
PARAMETER (ZERO=0.D0,ONE=1.D0,PFIVE=.5D0,OPF=1.5D0,C014=0.014D0)
DIMENSION DS(NSP), VX(NSP)
PI=4*ATAN(ONE)
TRD = ONE/3
FTRD = 4*TRD
TFTM = 2**FTRD-2
A0 = (4/(9*PI))**TRD
C X-alpha parameter:
ALP = 2 * TRD
IF (NSP .EQ. 2) THEN
D1 = MAX(DS(1),0.D0)
D2 = MAX(DS(2),0.D0)
D = D1 + D2
IF (D .LE. ZERO) THEN
EX = ZERO
VX(1) = ZERO
VX(2) = ZERO
RETURN
ENDIF
Z = (D1 - D2) / D
FZ = ((1+Z)**FTRD+(1-Z)**FTRD-2)/TFTM
FZP = FTRD*((1+Z)**TRD-(1-Z)**TRD)/TFTM
ELSE
D = DS(1)
IF (D .LE. ZERO) THEN
EX = ZERO
VX(1) = ZERO
RETURN
ENDIF
Z = ZERO
FZ = ZERO
FZP = ZERO
ENDIF
RS = (3 / (4*PI*D) )**TRD
VXP = -(3*ALP/(2*PI*A0*RS))
EXP_VAR = 3*VXP/4
IF (IREL .EQ. 1) THEN
BETA = C014/RS
SB = SQRT(1+BETA*BETA)
ALB = LOG(BETA+SB)
VXP = VXP * (-PFIVE + OPF * ALB / (BETA*SB))
EXP_VAR = EXP_VAR * (ONE-OPF*((BETA*SB-ALB)/BETA**2)**2)
ENDIF
VXF = 2**TRD*VXP
EXF = 2**TRD*EXP_VAR
IF (NSP .EQ. 2) THEN
VX(1) = VXP + FZ*(VXF-VXP) + (1-Z)*FZP*(EXF-EXP_VAR)
VX(2) = VXP + FZ*(VXF-VXP) - (1+Z)*FZP*(EXF-EXP_VAR)
EX = EXP_VAR + FZ*(EXF-EXP_VAR)
ELSE
VX(1) = VXP
EX = EXP_VAR
ENDIF
END
1.5.2 翻译后的C#程序
#region 版权信息
// 代码作者:
// 版权所有:
#endregion
using System;
using System.Text;
using System.Linq;
using System.Collections;
using System.Collections.Generic;
namespace Namespace1
{
public partial class Class1
{
/// <summary>
/// 子程序EXCHNG
/// </summary>
/// <param name="IREL">int型</param>
/// <param name="NSP">int型</param>
/// <param name="DS">double型数组</param>
/// <param name="EX">double型</param>
/// <param name="VX">double型数组</param>
public static void EXCHNG(int IREL, int NSP, out double[] DS, out double EX, out double[] VX)
{
// *****************************************************************
// Finds local exchange energy density and potential
// Adapted by J.M.Soler from routine velect of Froyen's
// pseudopotential generation program. Madrid, Jan'97. Version 0.5.
// **** Input ******************************************************
// INTEGER IREL : relativistic-exchange switch (0=no, 1=yes)
// INTEGER NSP : spin-polarizations (1=>unpolarized, 2=>polarized)
// REAL*8 DS(NSP) : total (nsp=1) or spin (nsp=2) electron density
// **** Output *****************************************************
// REAL*8 EX : exchange energy density
// REAL*8 VX(NSP) : (spin-dependent) exchange potential
// **** Units ******************************************************
// Densities in electrons/Bohr**3
// Energies in Hartrees
// *****************************************************************
// IMPLICIT DOUBLE(A-H,O-Z)
const double ZERO = 0.0;
const double ONE = 1.0;
const double PFIVE = .50;
const double OPF = 1.50;
const double C014 = 0.0140;
DS = new double[NSP + 1];
VX = new double[NSP + 1];
#region 定义本程序内用到的局部变量
double Z = 0.0;
double RS = 0.0;
double ALB = 0.0;
double A0 = 0.0;
double EXF = 0.0;
double FZP = 0.0;
double D = 0.0;
double SB = 0.0;
double TRD = 0.0;
double PI = 0.0;
double FZ = 0.0;
double EXP_VAR = 0.0;
double VXF = 0.0;
double BETA = 0.0;
double TFTM = 0.0;
double FTRD = 0.0;
double VXP = 0.0;
double D1 = 0.0;
double ALP = 0.0;
double D2 = 0.0;
#endregion 局部变量定义结束
PI = 4 * Math.Atan(ONE);
TRD = ONE / 3;
FTRD = 4 * TRD;
TFTM = Math.Pow(2, FTRD) - 2;
A0 = Math.Pow((4 / (9 * PI)), TRD);
// X-alpha parameter:
ALP = 2 * TRD;
if (NSP == 2)
{
D1 = Math.Max(DS[1], 0.0);
D2 = Math.Max(DS[2], 0.0);
D = D1 + D2;
if (D <= ZERO)
{
EX = ZERO;
VX[1] = ZERO;
VX[2] = ZERO;
return;
}
Z = (D1 - D2) / D;
FZ = (Math.Pow((1 + Z), FTRD) + Math.Pow((1 - Z), FTRD) - 2) / TFTM;
FZP = FTRD * (Math.Pow((1 + Z), TRD) - Math.Pow((1 - Z), TRD)) / TFTM;
}
else
{
D = DS[1];
if (D <= ZERO)
{
EX = ZERO;
VX[1] = ZERO;
return;
}
Z = ZERO;
FZ = ZERO;
FZP = ZERO;
}
RS = Math.Pow((3 / (4 * PI * D)), TRD);
VXP = -(3 * ALP / (2 * PI * A0 * RS));
EXP_VAR = 3 * VXP / 4;
if (IREL == 1)
{
BETA = C014 / RS;
SB = Math.Sqrt(1 + BETA * BETA);
ALB = Math.Log(BETA + SB);
VXP = VXP * (-PFIVE + OPF * ALB / (BETA * SB));
EXP_VAR = EXP_VAR * (ONE - OPF * Math.Pow(((BETA * SB - ALB) / Math.Pow(BETA, 2)), 2));
}
VXF = Math.Pow(2, TRD) * VXP;
EXF = Math.Pow(2, TRD) * EXP_VAR;
if (NSP == 2)
{
VX[1] = VXP + FZ * (VXF - VXP) + (1 - Z) * FZP * (EXF - EXP_VAR);
VX[2] = VXP + FZ * (VXF - VXP) - (1 + Z) * FZP * (EXF - EXP_VAR);
EX = EXP_VAR + FZ * (EXF - EXP_VAR);
}
else
{
VX[1] = VXP;
EX = EXP_VAR;
}
}
}
}
直接被纳入project,没有任何错误!
2 FORTRAN77语法概述
FORTRAN是世界上最早出现的高级编程语言,是工程界最常用的编程语言,它在科学计算中(如航空航天、地质勘探、天气预报和建筑工程等领域)发挥着极其重要的作用。经过40多年的发展,伴随着FORTRAN语言多次版本的更新及相应开发系统的出现,其功能不断完善,最新版本的开发系统几乎具备了VC、VB的所有特点,如图形界面编程、数据库等。目前,工科院校开设的计算机编程语言课首选仍然是FORTRAN :<
说实话,从科技发展的趋势来说这不是好事。您可以设想一下,如果需要用鹅毛笔抄写大量的古籍是什么感受!
强烈建议阅读《发掘C#特性赋予科学计算项目以威力》
2.1 FORTRAN77四则运算符
+ - * / ** (其中**表示乘方)
在表达式中按优先级次序由低到高为: +或-→*或/→**→函数→()
2.2 FORTRAN77变量类型
2.2.1 隐含约定:I-N规则
凡是以字母I,J,K,L,M,N六个字母开头的,即认为是整型变量,其它为实型变量。
2.2.2 用类型说明语句确定变量类型:可以改变I-N规则
INTEGER 整型
REAL 实型
DOUBLE PRECISION 双精度实型
COMPLEX 复型,赋值形式为(实部,虚部),如D=(8.76E+0.5,-67.8E-3),C=(3.0,6.3),如果含表达式则用CMPLX,如C=CMPLX(3.0*A,6.0+B)
LOGICAL 逻辑型,逻辑常量有“T”和“F”,“T”表示“.TRUE.”,“F”表示“.FALSE.”
CHARACTER*N 字符型,N为字符串长度,可以在变量名称后重新指定长度,如CHARACTER*8 STR1,STR2*10 ,赋值形式为STR2='I''M A BOY.'
2.2.3 用IMPLICIT语句将某一字母开头的全部变量指定为所需类型
如 IMPLICIT REAL (I,J)
三种定义的优先级别由低到高顺序为:I-N规则→IMPLICIT语句→类型说明语句,因此,在程序中IMPLICIT语句应放在类型说明语句之前。
2.2.4 数组的说明与使用
使用I-N规则时用DIMENSION说明数组,也可在定义变量类型同时说明数组,说明格式为:数组名(下标下界,下标上界),也可省略下标下界,此时默认为1,例:
DIMENSION IA(0:9),ND(80:99),W(3,2),NUM(-1:0),A(0:2,0:1,0:3)
REAL IA(10),ND(80:99)使用隐含DO循环进行数组输入输出操作:例如
WRITE(*,10) ('I=',I,'A=',A(I),I=1,10,2)
10FORMAT(1X,5(A2,I2,1X,A2,I4))
2.2.5 使用DATA语句给数组赋初值
变量表中可出现变量名,数组名,数组元素名,隐含DO循环,但不许出现任何形式的表达式:例如
DATA A,B,C/-1.0,-1.0,-1.0/
DATA A/-1.0/,B/-1.0/,C/-1.0/
DATA A,B,C/3*-1.0/CHARACTER*6 CHN(10)
DATA CHN/10*' '/INTEGER NUM(1000)
DATA (NUM(I),I=1,500)/500*0/,(NUM(I),I=501,1000)/500*1/
2.3 FORTRAN77程序书写规则
程序中的变量名,不分大小写;
变量名称是以字母开头再加上1到5位字母或数字构成,即变更名字串中只有前6位有效;
一行只能写一个语句;
程序的第一个语句固定为PROGRAM 程序名称字符串
某行的第1个字符至第5个字符位为标号区,只能书写语句标号或空着或注释内容;
某行的第1个字符为C或*号时,则表示该行为注释行,其后面的内容为注释内容;
某行的第6个字符位为非空格和非0字符时,则该行为上一行的续行,一个语句最多可有19个续行;
某行的第7至72字符位为语句区,语句区内可以任加空格以求美观;
某行的第73至80字符位为注释区,80字符位以后不能有内容。
2.4 FORTRAN77关系运算符
.GT. 大于
.GE. 天于或等于
.LT. 小于
.LE. 小于或等于
.EQ. 等于
.NE. 不等于 .AND. 逻辑与
.OR. 逻辑或
.NOT. 逻辑非
.EQV. 逻辑等
.NEQV. 逻辑不等
运算符优先级由高到低顺序为:()→**→*或/→+或-→.GT.或.GE.或.LT.或.LE.或.EQ.或.NE.→.NOT.→.AND.→.OR.→.EQV.或.NEQV
2.5 FORTRAN77语句
2.6 FORTRAN77选择判断语句
2.6.1 逻辑IF语句
IF (逻辑表达式) 程序语句
2.6.2 无ELSE块
IF (逻辑表达式) THEN
程序块
END IF
2.6.3 标准选择
IF (逻辑表达式) THEN
程序块1
ELSE
程序块2
END IF
2.6.4 多重选择块
IF (逻辑表达式1) THEN
程序块1
ELSE IF (逻辑表达式2) THEN
程序块2
ELSE IF (逻辑表达式2) THEN
程序块2
......
ELSE IF (逻辑表达式N) THEN
程序块N
ELSE
程序块N+1
END IF
2.7 FORTRAN77循环语句
2.7.1 GO TO语句
标号程序行
程序块
GO TO 标号
2.7.2 DO语句
DO 标号,记数变量=起始值,终止值,步距”语句,如
DO 标号,N=1,100,1
程序块
标号CONTINUE 7.3 DO WHILE 语句
DO 标号,WHILE(PI.EQ.3.14159)
程序块
标号CONTINUE 7.4 DO UNTIL语句
DO 标号,UNTIL (逻辑表达式)”语句,如
DO 标号,UNTIL(PI.GT.3.14159)
程序块
标号CONTINUE
2.8 FORTRAN77内部函数
2.9 FORTRAN77函数与子程序
2.9.1 FORTRAN77语句函数
当函数十分简单,用一条语句足以定义时(允许使用继续行)才用;
应该放在所有可执行语句之前和有关类型说明语句之后,是非执行语句;
只在其所在程序单位中有意义;
语句函数中的虚参就是变量名,不能是常量、表达式或数组元素等;
语句函数定义语句中的表达式可以包含已经定义过的语句函数、外部函数或内部函数。
语句函数通过表达式得一个函数值,此数值类型必须与函数名的类型一致。语句函数的使用同内部函数相同。
语句函数例子:
YMJ(R)=3.14159265*R*R
ZMJ=YMJ(5)
2.9.2 FORTRAN77自定义函数
定义格式:
类型说明 FUNCTION 函数名(虚拟参数1,虚拟参数2,……,虚拟参数N)
程序块(可以含有RETURN)
函数名=函数值
END
调用格式与内部函数相同。
9.3 FORTRAN77子程序
定义格式:
SUBROUTINE 子程序名(虚拟参数1,虚拟参数2,……,虚拟参数N)
程序块(可以含有RETURN)
END
调用格式:
CALL 子程序名(实在参数1,实在参数2,……,实在参数N)
数据块子程序:只是用来给有名公用区中的变量赋初值,格式如下:
BLOCK DATA 子程序名
DATA语句块
END