LAPACK: Czy operacje na macierzach ze spakowanymi zasobami są szybsze?

Question

Chcę tridiagonalizować prawdziwą macierz symetryczną za pomocą Fortran i LAPACK. LAPACK zasadniczo udostępnia dwie procedury, jedną działającą na pełnej macierzy, a drugą na macierzy w zapakowanym magazynie. Podczas gdy ta ostatnia z pewnością zużywa mniej pamięci, zastanawiałem się, czy cokolwiek można powiedzieć o różnicy prędkości?

Answer 1

Jest to oczywiście pytanie empiryczne: ale ogólnie rzecz biorąc nic nie przychodzi za darmo, a mniej pamięci/więcej czasu pracy jest dość powszechnym kompromisem.

W tym przypadku indeksowanie danych jest bardziej skomplikowane dla zapakowanego przypadku, więc podczas przechodzenia przez macierz koszt pobrania danych jest nieco wyższy. (Skomplikowanie tego obrazu jest takie, że w przypadku macierzy symetrycznych procedury lapack również przyjmują pewien rodzaj upakowania - że macie tylko górny lub dolny komponent macierzy dostępnej).

Wcześniej borykałem się z problemem własnym już dziś, więc wykorzystam to jako punkt odniesienia dla pomiarów; próby z prostym symetrycznym przypadku badanej (Herdon matrycy z http://people.sc.fsu.edu/~jburkardt/m_src/test_mat/test_mat.html) i porównanie ssyevd z 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 

Jest o różnicę 18%, co trzeba przyznać, jest większa niż oczekiwano (również nieco większy błąd w zapakowanej skrzynce?). Jest to w przypadku MKL firmy Intel. Różnica w wydajności zależy oczywiście od twojej macierzy, oczywiście, jak eriktous wskazuje i od problemu, który robisz; im bardziej losowy dostęp do macierzy musisz zrobić, tym gorsze będą koszty ogólne. Kod, którego użyłem, jest następujący:

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