在大多数计算机架构中,这种差异预计会很大,至少在原则上是这样。
矩阵向量乘法是一种受内存限制的计算,因为内存的重用率很低。 v 的所有 (N) 个分量都被重用来计算 u 的每个元素,但矩阵 (N^2) 的每个元素仅使用一次。如果我们考虑典型内存的延迟(例如,https://gist.github.com/hellerbarde/2843375 https://gist.github.com/hellerbarde/2843375)(小于)100ns 与执行浮点运算所需的时间(小于 1ns)相比,我们看到大部分时间都花在从数组加载和存储值上。
我们仍然可以实现缓存友好的,即尽可能具有数据局部性。由于内存是以行的形式加载到缓存中的,因此我们必须尽可能地使用已加载的缓存行。这就是为什么访问连续的内存区域可以减少从内存加载数据所花费的时间。
为了支持这一点,让我们尝试一个非常简单的代码:
program mv
integer, parameter :: n=10000
real, allocatable :: M(:,:), v(:), u(:)
real :: start, finish
integer :: i, j
allocate(M(n,n),v(n),u(n))
call random_number(M)
call random_number(v)
u(:)=0.
call cpu_time(start)
do i=1,n
do j=1,n
! non-contiguous order
u(i)=u(i)+M(i,j)*v(j)
! contiguous order
! u(i)=u(i)+M(j,i)*v(j)
enddo
enddo
call cpu_time(finish)
print*,'elapsed time: ',finish-start
end program mv
一些结果:
non-contiguous order contiguous order
gfortran -O0 1. 0.5
gfortran -O3 0.3 0.1
ifort -O0 1.5 0.85
ifort -O3 0.037 0.035
正如您所看到的,差异是在没有优化的情况下编译的显着差异。启用优化 gfortran 仍然显示出显着差异,而使用 ifort 则只有很小的差异。查看编译器报告,编译器似乎交换了循环,从而导致内部循环的连续访问。
然而,我们是否可以说具有行优先排序的语言对于矩阵向量计算更有效?不,我不能这么说。不仅因为编译器可以补偿差异。代码本身并不知道 M 的行和列的所有信息:它基本上知道 M 有两个索引,其中一个(取决于语言)在内存中是连续的。对于矩阵向量,数据局部性的最佳方法是将“快速”索引映射到矩阵行索引。您可以使用“行优先”和“列优先”语言来实现此目的。您只需根据此存储 M 的值即可。举个例子,如果你有“代数”矩阵
[ M11 M12 ]
M = [ ]
[ M21 M22 ]
您将其存储为“计算矩阵”
C ==> M[1,1] = M11 ; M[1,2] = M12 ; M[2,1] = M21 ; M[2,2] = M22
Fortran ==> M[1,1] = M11 ; M[2,1] = M12 ; M[1,2] = M21 ; M[2,2] = M22
这样你就总是在“代数矩阵”行中连续。计算机对初始矩阵一无所知,但我们知道计算矩阵是代数矩阵的转置版本。在这两种情况下,我都会让内部循环迭代连续索引,最终结果将是相同的向量。
在复杂的代码中,如果我已经分配了矩阵并用值填充了矩阵,并且无法决定存储转置矩阵,则“行优先”语言可能会提供最佳性能。但是,交换循环(参见https://en.wikipedia.org/wiki/Loop_interchange https://en.wikipedia.org/wiki/Loop_interchange)由英特尔编译器自动完成,并由 BLAS 实现完成(请参阅http://www.netlib.org/lapack/explore-html/db/d58/sgemv_8f_source.html http://www.netlib.org/lapack/explore-html/db/d58/sgemv_8f_source.html),将差异减少到非常小的差异值。因此,使用 Fortran 您可以选择:
do j=1,n
do i=1,n
u(i)=u(i)+M(i,j)*v(j)
enddo
enddo