LAPACK：打包存储矩阵上的操作是否更快？

Question

我想使用 Fortran 和 LAPACK 对实对称矩阵进行三对角化。 LAPACK 基本上提供了两个例程，一个在完整矩阵上运行，另一个在打包存储中的矩阵上运行。虽然后者肯定使用更少的内存，但我想知道是否可以对速度差异进行说明？

Answer 1

当然，这是一个经验问题：但一般来说，没有什么是免费的，更少的内存/更多的运行时间是一个非常常见的权衡。

在这种情况下，对于打包的情况，数据的索引更加复杂，因此当您遍历矩阵时，获取数据的成本会稍高一些。 （使这张图变得复杂的是，对于对称矩阵，lapack 例程还假设某种包装 - 您只有矩阵的上部或下部组件可用）。

今天早些时候我正在研究一个特征问题，所以我将用它作为测量基准；尝试使用简单的对称测试用例（赫登矩阵，来自 http:// /people.sc.fsu.edu/~jburkardt/m_src/test_mat/test_mat.html ），并将 ssyevd 与sspevd

$ ./eigen2 500
 Generating a Herdon matrix: 
 Unpacked array:
 Eigenvalues L_infty err =   1.7881393E-06
 Packed array:
 Eigenvalues L_infty err =   3.0994415E-06
 Packed time:   2.800000086426735E-002
 Unpacked time:   2.500000037252903E-002

$ ./eigen2 1000
 Generating a Herdon matrix: 
 Unpacked array:
 Eigenvalues L_infty err =   4.5299530E-06
 Packed array:
 Eigenvalues L_infty err =   5.8412552E-06
 Packed time:   0.193900004029274     
 Unpacked time:   0.165000006556511  

$ ./eigen2 2500
 Generating a Herdon matrix: 
 Unpacked array:
 Eigenvalues L_infty err =   6.1988831E-06
 Packed array:
 Eigenvalues L_infty err =   8.4638596E-06
 Packed time:    3.21040010452271     
 Unpacked time:    2.70149993896484 

大约有 18% 的差异，我必须承认这个差异比我预期的要大（包装箱的误差也稍大一些？）。这是Intel的MKL。当然，正如 eriktous 指出的那样，性能差异通常取决于您的矩阵，以及您正在做的问题；对矩阵的随机访问越多，开销就越严重。我使用的代码如下：

program eigens
      implicit none

      integer :: nargs,n  ! problem size 
      real, dimension(:,:), allocatable :: A, B, Z
      real, dimension(:), allocatable :: PA
      real, dimension(:), allocatable :: work
      integer, dimension(:), allocatable :: iwork
      real, dimension(:), allocatable :: eigenvals, expected
      real :: c, p
      integer :: worksize, iworksize
      character(len=100) :: nstr
      integer :: unpackedclock, packedclock 
      double precision :: unpackedtime, packedtime
      integer :: i,j,info

! get filename
      nargs = command_argument_count()
      if (nargs /= 1) then
          print *,'Usage: eigen2 n'
          print *,'       Where n = size of array'
          stop
      endif
      call get_command_argument(1, nstr)
      read(nstr,'(I)') n
      if (n < 4 .or. n > 25000) then
          print *, 'Invalid n ', nstr
          stop
      endif


! Initialize local arrays    

      allocate(A(n,n),B(n,n))
      allocate(eigenvals(n)) 

! calculate the matrix - unpacked

      print *, 'Generating a Herdon matrix: '

      A = 0.
      c = (1.*n * (1.*n + 1.) * (2.*n - 5.))/6.
      forall (i=1:n-1,j=1:n-1)
        A(i,j) = -1.*i*j/c
      endforall
      forall (i=1:n-1)
        A(i,i) = (c - 1.*i*i)/c
        A(i,n) = 1.*i/c
      endforall
      forall (j=1:n-1)
        A(n,j) = 1.*j/c
      endforall
      A(n,n) = -1./c
      B = A

      ! expected eigenvalues
      allocate(expected(n))
      p = 3. + sqrt((4. * n - 3.) * (n - 1.)*3./(n+1.))
      expected(1) = p/(n*(5.-2.*n))
      expected(2) = 6./(p*(n+1.))
      expected(3:n) = 1.

      print *, 'Unpacked array:'
      allocate(work(1),iwork(1))
      call ssyevd('N','U',n,A,n,eigenvals,work,-1,iwork,-1,info)
      worksize = int(work(1))
      iworksize = int(work(1))
      deallocate(work,iwork)
      allocate(work(worksize),iwork(iworksize))

      call tick(unpackedclock)
      call ssyevd('N','U',n,A,n,eigenvals,work,worksize,iwork,iworksize,info)
      unpackedtime = tock(unpackedclock)
      deallocate(work,iwork)

      if (info /= 0) then
           print *, 'Error -- info = ', info
      endif
      print *,'Eigenvalues L_infty err = ', maxval(eigenvals-expected)


      ! pack array

      print *, 'Packed array:'
      allocate(PA(n*(n+1)/2))
      allocate(Z(n,n))
      do i=1,n 
        do j=i,n
           PA(i+(j-1)*j/2) = B(i,j)
        enddo
      enddo

      allocate(work(1),iwork(1))
      call sspevd('N','U',n,PA,eigenvals,Z,n,work,-1,iwork,-1,info)
      worksize = int(work(1))
      iworksize = iwork(1)
      deallocate(work,iwork)
      allocate(work(worksize),iwork(iworksize))

      call tick(packedclock)
      call sspevd('N','U',n,PA,eigenvals,Z,n,work,worksize,iwork,iworksize,info)
      packedtime = tock(packedclock)
      deallocate(work,iwork)
      deallocate(Z,A,B,PA)

      if (info /= 0) then
           print *, 'Error -- info = ', info
      endif
      print *,'Eigenvalues L_infty err = ', &
      maxval(eigenvals-expected)

      deallocate(eigenvals, expected)


      print *,'Packed time: ', packedtime
      print *,'Unpacked time: ', unpackedtime


contains
    subroutine tick(t)
        integer, intent(OUT) :: t

        call system_clock(t)
    end subroutine tick

    ! returns time in seconds from now to time described by t
    real function tock(t)
        integer, intent(in) :: t
        integer :: now, clock_rate

        call system_clock(now,clock_rate)

        tock = real(now - t)/real(clock_rate)
    end function tock

end program eigens