行主排序对于矩阵向量乘法更有效吗



如果M是一个n x m矩阵,vu是向量,那么就索引而言,矩阵向量乘法看起来像u[i] = sum(M[i,j] v_j, 1 <= j <= m)。由于v是一个矢量,其元素可能存储在面向数值计算的语言的连续存储位置中。如果M按行主顺序存储(如在C、Mathematica和Pascal中),则随着j的递增,总和中的后续M[i,j]也存储在连续的存储位置中,从而使迭代非常有效。如果它是按列主顺序存储的(如Fortran、Matlab、R和Julia),则递增j需要移动与外部矩阵步幅相等的多个存储位置,在这种情况下等于n。对于具有许多行的矩阵来说,这似乎效率较低。(对于矩阵-矩阵乘法,问题不会出现,因为在任何一种排序约定下,增加求和索引都需要在一个矩阵的内存中移动一个主要步长。)

与乘法和加法运算相比,在大多数计算机体系结构中,在内存中移动一个单位和移动多个单位之间的差异是明显的还是可以忽略的?(我猜"可以忽略不计",因为在实践中Fortran通常至少和C一样快,但有人能详细说明为什么吗?)

至少在原则上,大多数计算机体系结构中的差异预计会很大。

矩阵向量乘法是一种有内存限制的计算,因为内存的重用率很低。v的所有(N)个分量都被重用来计算u的每个元素,但矩阵(N^2)的每个元素只使用一次。如果我们考虑典型存储器的延迟(参见例如。,https://gist.github.com/hellerbarde/2843375)与执行浮点运算所需的时间(小于1ns)相比,由于(小于)100ns,我们发现大部分时间都花在了从数组加载和存储值上。

我们仍然可以实现它的缓存友好,即尽可能多地具有数据局部性。由于内存是以行的形式加载到缓存中的,因此我们必须尽可能多地使用已加载的缓存行。这就是为什么访问连续的内存区域可以减少从内存加载数据所花费的时间。

为了支持这一点,让我们尝试一个非常简单的代码:

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有两个索引,其中一个索引在内存中是连续的,具体取决于语言。对于矩阵向量,最好的数据局部性是将"快速"索引映射到矩阵行索引。您可以使用"行专业"one_answers"列专业"语言来实现这一点。你只需要根据这个来存储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)由英特尔编译器自动完成,由BLAS实现自动完成(请参阅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

相关内容

  • 没有找到相关文章

最新更新