LAPACK: Операции с упакованными матрицами памяти быстрее? - PullRequest
Новая
       22

LAPACK: Операции с упакованными матрицами памяти быстрее?

9 голосов
/ 20 января 2012

Я хочу тридиагонализировать реальную симметричную матрицу, используя Fortran и LAPACK.LAPACK в основном предоставляет две подпрограммы, одна из которых работает на полной матрице, а другая - на матрице в упакованном хранилище.Хотя последний, безусловно, использует меньше памяти, мне было интересно, можно ли что-нибудь сказать о разнице в скорости?

1 Ответ

9 голосов
/ 20 января 2012

Это, конечно, эмпирический вопрос: но в общем, ничего не бывает бесплатно, и меньше памяти / больше времени выполнения - довольно распространенный компромисс.

В этом случае индексация данных является более сложнойдля упакованного случая, когда вы пересекаете матрицу, стоимость получения ваших данных немного выше.(Эта картина усложняется тем, что для симметричных матриц подпрограммы Лапака также предполагают определенный вид упаковки - у вас есть только верхний или нижний компонент доступной матрицы).

Ранее сегодня я возился с собственной проблемой, поэтому буду использовать ее в качестве эталона измерения;попытка с простым симметричным тестовым примером (матрица Хердона, от 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%, которую я должен признатьбольше, чем я ожидал (также с немного большей ошибкой для упакованного случая?).Это с MKL от Intel.Разница в производительности, конечно же, будет зависеть от вашей матрицы в целом, как указывает на нее, и от проблемы, которую вы решаете;чем больше случайного доступа к матрице вы должны сделать, тем хуже будут издержки.Код, который я использовал, выглядит следующим образом: 

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
...