LAPACK: Các hoạt động trên ma trận lưu trữ đóng gói có nhanh hơn không?

Question

Tôi muốn tridiagonalize một ma trận đối xứng thực sự bằng cách sử dụng Fortran và LAPACK. LAPACK về cơ bản cung cấp hai thói quen, một hoạt động trên ma trận đầy đủ, một hoạt động khác trên ma trận trong bộ nhớ đóng gói. Trong khi sau này chắc chắn sử dụng bộ nhớ ít hơn, tôi đã tự hỏi nếu bất cứ điều gì có thể nói về sự khác biệt tốc độ?

Answer 1

Đó là một câu hỏi thực nghiệm, tất nhiên: nhưng nói chung, không có gì đến miễn phí, và ít bộ nhớ/thời gian chạy hơn là một sự cân bằng khá phổ biến.

Trong trường hợp này, việc lập chỉ mục cho dữ liệu phức tạp hơn đối với trường hợp được đóng gói, vì vậy khi bạn duyệt qua ma trận, chi phí nhận dữ liệu của bạn cao hơn một chút. (Phức tạp hình ảnh này là đối với ma trận đối xứng, các thói quen lapack cũng giả định một loại bao bì nhất định - rằng bạn chỉ có thành phần trên hoặc dưới của ma trận có sẵn).

Tôi đã rối tung xung quanh với một giải pháp bản địa trước đó, vì vậy tôi sẽ sử dụng nó làm điểm chuẩn đo lường; cố gắng với một trường hợp thử nghiệm đối xứng đơn giản (Ma trận Herdon, từ http://people.sc.fsu.edu/~jburkardt/m_src/test_mat/test_mat.html), và so sánh ssyevd với 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 

Có khoảng một sự khác biệt 18%, mà tôi phải thừa nhận là lớn hơn tôi mong đợi (cũng với một hơi lớn lỗi cho trường hợp được đóng gói?). Đây là với MKL của intel. Sự khác biệt về hiệu suất sẽ phụ thuộc vào ma trận của bạn nói chung, tất nhiên, như là điểm nổi bật, và về vấn đề bạn đang làm; sự truy cập ngẫu nhiên hơn vào ma trận bạn phải làm, thì mức phí trên sẽ càng thấp. Mã tôi đã sử dụng như sau:

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