FORTRAN在计算上的优势?

Posted 2023-03-20

人言FORTRAN在工程计算方面有优势，现在工程计算软件还多用FORTRAN开发，请问FORTRAN在计算上到底有什么优势（和C/C++比较），究竟好在哪？最好谈得深入点，谢谢。
1楼的说了半天，根本没解释我关心的“为什么适合工程计算”;或者说，fortran能写的工程计算软件c/c++也能写，和fortran相比有什么差别，为什么会有这样的差别
谢谢gzlu从库的角度分析fortran的优势，还有其他解释吗？

优势如下：

1、FORTRAN历史悠久

自成立以来，距第一个FORTRAN计划诞生已有近50年的时间。长期以来，FORTRAN一直是科学计算语言的唯一选择。相比之下，C / C ++和Matlab迟了将近20年。

2、FORTRAN具有严格的语法要求

FORTRAN更适合严格的科学计算领域。就像C不会检查数组越界一样，如果访问了错误的地址，则计算结果是不可预测的。

3、FORTRAN可以直接对数组和复数进行运算

像C来完成两个矩阵的加法执行两个循环，或调用具有5个形式参数的函数一样，尽管C ++可以完成C = A + B的简单形式，但是基于矩阵类的定义及其成员函数，在重载运算符的基础上，矩阵的初始化成为矩阵类的初始化，增加了工作量。

而FORTRAN可以直接完成矩阵C = A + B。

4、FORTRAN是编译语言

与Matlab相比，FORTRAN是一种编译语言，Matlab与解释型语言类似，并且循环效率更高，因此，如果Matlab要使用很多循环，则必须调用C / C ++或FORTRAN程序。

扩展资料：

Fortran是矩阵运算最方便的语言，通常用于高能物理学，气象学，航空航天等领域的方程计算问题，计算矩阵（array）的计算非常有效。

MKL和IMSL是著名的数学库，Fortran在某些领域仍然受欢迎的原因是，一方面，Fortran具有强大的数学运算能力和高效率。 另一方面，这是因为前辈对于许多问题都有成熟的代码，所以直接使用即可，而不是先考虑算法的答案，然后再对其进行改进。

矩阵越大，方程越复杂，Fortran和matlab这两种语言之间的差异就越大。

参考技术A 我是半年前学的FORTRAN 可以说是刚入门
不过我的老师是个牛人，听他说的这些：
FORTRAN作为为一专门的科学计算语言，其地位到目前为止还是不可动摇的。现今常用与科学计算方面的语言还有，C/C++语言,Matlab语言，FORTRAN与它们比起来有着得天独厚的优越性。
1、FORTRAN历史悠久，底子后，从第一个FORTRAN程序诞生到现在接近50年啦，很长一段时间内是科学计算语言的唯一选择，期间积累的大量的正确、可靠的程序，尽管FORTRAN标准该了多次，但由于其向下兼容，很多程序是招值即来，来之能战。现比之下，C/C++，Matlab晚了近20年
2、FORTRAN语法要求严格，更适合严谨的科学计算领域。像C中并不对数组越界进行检查，如果访问到错误的地址其计算结果是无法预测的。
3、FORTRAN可直接对数组和复数进行运算。像C完成两个矩阵的加要进行两次循环，或者是调用一个带5个形参的函数，C++虽然可以完成C=A+B的简单形式，但是它建立在定义矩阵类及其成员函数，再对运算符进行重载的基础上，前期工作繁琐，而且对矩阵的初始化变成了对矩阵类的初始化
，增加的工作量。而FORTRAN可以直接完成矩阵的C=A+B
4、与Matlab相比，FORTRAN是编译语言，Matlab类似与解释语言，循环效率高的多，所以如果Matlab要大量使用循环就不得不调用C/C++或FORTRAN程序。
5、并行计算中FORTRAN的地位是不可替代的，在并行计算领域FORTRAN是独领风骚本回答被提问者采纳 参考技术B 它只针对科学运算
FORTRAN语言中文叫公式翻译语言，是IBM公司发明的计算机高级语言。是我国六七十年代流行的两大编程语言之一,用于较大的机种。七十年代末起，用于PC。

FORTRAN语言逻辑性强，程序结构清晰，语法语义简捷好懂，特别适合用于科学计算，数据采集处理，调用绘图库（例如GKS,DISPLA等）可以绘图 。大型 MainFrame 计算机，DEC 计算机 等都用 FORTRAN。

从FORTRAN90开始，加入了可视化。现在叫Compag Visual FORTRAN。 可以在 MS Visual Studio 里工作。

凡C,C++ 能干的活，FORTRAN都能做。

近十年来 FORTRAN 不如 C,C++ 流行。主要因为FORTRAN编译器价格贵。不像C,C++ 有免费的版本。unix 和 工作站流行推动了C的传播，同时MainFrame，DEC的退役，老互联网BITNET消亡，使FORTRAN失去了主要的硬件基础。微软视窗统治了PC世界，与微软视窗紧密结合的VC++比FORTRAN好用。FORTRAN能干的活，C,C++ 都能干。 参考技术C 单纯从计算机语言上对比，没什么差别。但是，早期大量科学与工程计算程序是用fortran写的，留下了大批可用代码。 参考技术D 我觉得当年之所以这么说，是因为FORTRAN有很多数学上用的专业程序库。其实如果不是正好有所需的程序库之外，它没有什么优势。实在要用它的库，还可以用混合语言调用呢。

如何在 Fortran 中有效地计算矩阵内积？

【中文标题】如何在 Fortran 中有效地计算矩阵内积？

【英文标题】：How to efficiently calculate matrix inner product in Fortran?

【发布时间】：2018-04-08 13:16:49

【问题描述】：

我正在尝试计算类似于 Fortran 中加权矩阵内积的东西。我用于计算内积的当前脚本如下

! --> In
real(kind=8), intent(in), dimension(ni, nj, nk, nVar) :: U1, U2
real(kind=8), intent(in), dimension(ni, nj, nk) :: intW

! --> Out
real(kind=8), intent(out) :: innerProd

! --> Local
integer :: ni, nj, nk, nVar, iVar

! --> Computing inner product
do iVar = 1, nVar
    innerProd = innerProd + sum(U1(:,:,:,iVar)*U2(:,:,:,iVar)*intW)
enddo

但是我发现我目前使用的上述脚本效率不是很高。使用 NumPy 在 Python 中可以执行相同的操作，如下所示，

import numpy as np 
import os

# --> Preventing numpy from multi-threading
os.environ['OPENBLAS_NUM_THREADS'] = '1'
os.environ['MKL_NUM_THREADS'] = '1'   

innerProd = 0

# --> Toy matrices
U1 = np.random.random((ni,nj,nk,nVar))
U2 = np.random.random((ni,nj,nk,nVar))
intW = np.random.random((ni,nj,nk))

# --> Reshaping 
U1 = np.reshape(np.ravel(U1), (ni*nj*nk, nVar))
U2 = np.reshape(np.ravel(U1), (ni*nj*nk, nVar))
intW = np.reshape(np.ravel(intW), (ni*nj*nk))

# --> Calculating inner product
for iVar in range(nVar):
    innerProd = innerProd + np.dot(U1[:, iVar], U2[:, iVar]*intW)

使用 Numpy 的第二种方法似乎比使用 Fortran 的方法快得多。对于ni = nj = nk = nVar = 130的具体情况，两种方法所用时间如下

 fortran_time = 25.8641 s
 numpy_time = 6.8924 s

我尝试使用来自 BLAS 的ddot 改进我的 Fortran 代码，如下所示，

do iVar = 1, nVar
    do k = 1, nk
        do j = 1, nj
            innerProd = innerProd + ddot(ni, U1(:,j,k,iVar), 1, U2(:,j,k,iVar)*intW(:,j,k), 1)
        enddo
    enddo
enddo

但时间并没有明显改善。对于ni = nj = nk = nVar = 130 的情况，上述方法所用的时间是~24s。 （我忘了提到我用'-O2'选项编译了Fortran代码以优化性能）。

不幸的是，Fortran 中没有用于逐元素矩阵乘法的 BLAS 函数。而且我不想在 Fortran 中使用 reshape，因为与 Python 不同，Fortran 中的 reshape 会导致将我的数组复制到一个新数组，从而导致更多的 RAM 使用。

有没有什么办法可以加快 Fortran 的性能，从而接近 Numpy 的性能？

【问题讨论】：

请注意 kind=8 是不可移植且丑陋的。

你确定 numpy 使用单核吗？

@AndrasDeak 是的，我很确定 numpy 只使用一个内核。我使用 htop 检查了核心使用情况。由于多线程，我非常有信心 NumPy 更快，但令我惊讶的是我错了。

我想知道将您的权重数组复制到完整大小并摆脱循环是否会有所帮助。你能试试吗？ （我留下了这个评论，然后以为我很困惑并删除了它，但现在我再次认为这实际上与您的内在产品相对应；但我仍然不确定）。

我已经尝试过你的代码，循环在 0.6 秒内完成。您如何安排问题的时间？

【参考方案1】：

您可能没有按照自己的想法进行计时。这是一个完整的fortran示例

program test                                                        
    use iso_fortran_env, r8 => real64                               
    implicit none                                                   

    integer, parameter :: ni = 130, nj = 130, nk = 130, nvar = 130  
    real(r8), allocatable :: u1(:,:,:,:), u2(:,:,:,:), w(:,:,:)     
    real(r8) :: sum, t0, t1                                         
    integer :: i,j,k,n                                              

    call cpu_time(t0)                                               
    allocate(u1(ni,nj,nk,nvar))                                     
    allocate(u2(ni,nj,nk,nvar))                                     
    allocate(w(ni,nj,nk))                                           
    call cpu_time(t1)                                               
    write(*,'("allocation time(s):",es15.5)') t1-t0                 

    call cpu_time(t0)                                               
    call random_seed()                                              
    call random_number(u1)                                          
    call random_number(u2)                                          
    call random_number(w)                                           
    call cpu_time(t1)                                               
    write(*,'("random init time (s):",es15.5)') t1-t0               

    sum = 0.0_r8                                                    
    call cpu_time(t0)                                               
    do n = 1, nvar                                                  
        do k = 1, nk                                                
            do j = 1, nj                                            
                do i = 1, ni                                        
                    sum = sum + u1(i,j,k,n)*u2(i,j,k,n)*w(i,j,k)    
                end do                                              
            end do                                                  
        end do                                                      
    end do                                                          
    call cpu_time(t1)                                               
    write(*,'("Sum:",es15.5," time(s):",es15.5)') sum, t1-t0        

end program

还有输出：

$ gfortran -O2 -o inner_product inner_product.f90            
$ time ./inner_product 
allocation time(s):    3.00000E-05
random init time (s):    5.73293E+00
Sum:    3.57050E+07 time(s):    5.69066E-01

real    0m6.465s
user    0m4.634s
sys 0m1.798s

在这个 fortran 代码中计算内积的时间不到 10%。你如何/什么时间是非常重要的。你确定你在 fortran 和 python 版本中的时间是一样的吗？你确定你只计时 inner_product 计算？

【参考方案2】：

这样可以避免制作任何副本。 （注意 blas ddot 方法仍然需要为元素产品制作副本）

   subroutine dot3(n,a,b,c,result)
   implicit none
   real(kind=..) a(*),b(*),c(*),result
   integer i,n
   result=0
   do i=1,n
    result=result+a(i)*b(i)*c(i)
   enddo
   end

dot3 是外部的，意味着 not 在模块/包含结构中。 kind 显然应该匹配主声明。

在主代码中：

  innerprod=0
  do iVar = 1, nVar 
  call dot3(ni*nj*nk, U1(1,1,1,iVar),U2(1,1,1,iVar),intW,result)
  innerProd=innerProd+result
  enddo

【讨论】：

是的，我刚才尝试了这种方法。但是运行时间仍然没有改善。它仍然是~24-25 s。我很惊讶我在 Fortran 中尝试过的所有方法都具有几乎相同的运行时间。

@AdhityaRavi 我不会对此感到惊讶。在金属附近，您只需要遍历内存中的每个项目，乘以三个双精度数，然后总结结果；没有太大的改进空间。我们可以做的是尽量减少开销，例如创建临时数组。只有在高级语言中，根据您的方法可能会出现巨大差异，例如本机 python 循环与 numpy 魔术。但是，如果 4 的差异因素不想消失，无论你做什么，我仍然会怀疑多核恶作剧......

@AndrasDeak 我猜你是对的。我也感觉到 numpy 以某种方式是多线程的，因为我找不到其他解释来说明它快 4 到 5 倍。我尝试使用os.environ['OPENBLAS_NUM_THREADS'] = '1' 和os.environ['MKL_NUM_THREADS'] = '1' 来抑制这种情况。但是，我不确定如何确保 numpy 使用单核。

@AdhityaRavi 众所周知，影响这一点的 envvar 取决于 numpy 使用的 blas 库，所以它可能是其他东西，例如 OMP_NUM_THREADS 我认为。不过，您之前提到的 CPU 负载应该是确定的。无论如何，您始终可以通过从 numpy 版本中删除线程限制开关来检查它是否运行得快 n 倍（以及 CPU 负载是否相应增加）。

请注意，过程是外部的、内部的、模块的还是其他任何对序列关联处理参数的方式没有影响。

【参考方案3】：

比较 Numpy 和 Fortran 代码时，我有同样的观察结果。

差异原来是 BLAS 的版本，我发现使用 netlib 中的 DGEMM 类似于循环，并且比 OpenBLAS 慢大约三倍（请参阅this 答案中的配置文件）。

对我来说最令人惊讶的是，OpenBLAS 提供的代码比仅编译 Fortran 三重嵌套循环要快得多。看来这就是 GotoBLAS 的重点，它是用汇编代码为处理器架构手写的。

即使定时正确、正确排序循环、避免复制并使用每个优化标志（在 gfortran 中），性能仍然比 OpenBLAS 慢大约三倍。我没有尝试过 ifort 或 pgi，但我想知道这是否解释了@kvantour 的赞成评论“我的循环在 0.6 秒内完成”（注意，在某些实现中，内在的 matmul 被 BLAS 取代）。

【讨论】：

我不再从事我发布此问题的项目。但是，这听起来像是一个非常合理的解释。感谢您的回答，我会将此线程转发给仍在使用慢速 fortran 代码的同事，以检查我的情况是否也发生了这种情况，以便我可以关闭;)。