Fortran的矩阵乘法如何编写更快

10天前浏览652

Fortran是一门古老的适用于数值计算的高级语言，有着悠久的历史。1981年以前，Fortran是世界上最流行的编程语言。各个领域的科学家，工程人员编写了不计其数的Fortran代码用于计算航空，大坝，大气等领域的科学计算。在知名科幻小说三体中，叶文洁第一次进入红岸基地，就对Fortran语言产生了兴趣：

...

叶文洁第一次走进主机房时，看到一排阴极射线管显示屏，她惊奇地发现，屏幕上竟滚动着排排程序代码，可以通过键盘随意-进行编辑和调试。而她在大学里使用计算机时，代码都写在一张张打格的程序纸上，再通过打字机噼噼啪啪地打到纸带上。她听说过从键盘和屏幕输入这回事，现在竟然真-的看到了。但更令她吃惊的是这里的软件技术。她知道了一种叫FORTRAN的东西（注：第一代计算机高级语言），竟能用接近自然语言的代码编写程序。能将数学公-式直接写到代码里，它的编程效率比机器码汇编不知高了多少倍

...

在有限元计算领域，Ansys，Abaqus，Adina和LS-Dyna以及近年开源的Radioss的求解器早期均采用Fortran编写，现在是否依然采用Fortran我们不得而知，但是可以猜测的是，即使在今天，这些商业软件中Fortran依然占据着不小的比重，比如，在Ansys/apdl和LS-DYNA中，有时候计算时的报错提示就是Fortran编译器的报错。实际上，在2024年的今天，高性能计算可选的语言实际上并不多，大部分时候基本上还是Fortran，C/C++。其他可能还有诸如Julia等。

在使用Fortran编写数值计算程序时，矩阵乘法是最常见的内容之一，在Fortran中，矩阵乘法有多种方式可以实现，自带的库函数matmul就能直接实现矩阵乘法。在实际中，对于较大的矩阵，通常也并不直接采用matmul，可能会自行编写或者调用MKL库的DGEMM函数完成。本文不讨论调库操作，仅对自行编写矩阵乘法时的效率进行讨论。

矩阵乘法的具体算法：

根据该算法，我们可以编写以下Fortran代码：

program mainimplicit noneinteger,parameter::n=3000real(8),allocatable::a(:,:),b(:,:),c(:,:)integer::i,j,kreal(8)::t1,t2allocate(a(n,n))allocate(b(n,n))allocate(c(n,n))a=1.0b=2.0c=0.0call CPU_TIME(t1)do i=1,ndo j=1,ndo k=1,nc(i,j)=c(i,j)+a(i,k)*b(k,j)enddoenddoenddocall CPU_TIME(t2)write(*,*)"time",t2-t1write(*,*)c(n,n)deallocate(a)deallocate(b)deallocate(c)end program main

采用gfortran编译-O3优化，运行：

3000阶稠密矩阵所花时间为99.78s，还是比较慢的。

接下来，注意到在Fortran中，数组以列优先存储，即a(1,1)的后一个元素是a(2,1)，因此实际上应尽量使得循环中相邻两次循环取的地址连续，即内存循环变量尽量放在数组的第一个下标。而在上面的程序中，仅b的下标k为最内层循环变量，而a和c的第一下标均为最外层循环变量。因此可以考虑将循环变量交换顺序，得到以下代码：

program mainimplicit noneinteger,parameter::n=3000real(8),allocatable::a(:,:),b(:,:),c(:,:)integer::i,j,kreal(8)::t1,t2allocate(a(n,n))allocate(b(n,n))allocate(c(n,n))a=1.0b=2.0c=0.0call CPU_TIME(t1)do j=1,n            !这里交换了i，j和k的顺序，顺序采用j-k-ido k=1,ndo i=1,nc(i,j)=c(i,j)+a(i,k)*b(k,j)enddoenddoenddocall CPU_TIME(t2)write(*,*)"time",t2-t1write(*,*)c(n,n)deallocate(a)deallocate(b)deallocate(c)end program main