LAPACK: पैक किए गए स्टोरेज मैट्रिक्स पर ऑपरेशन तेजी से हैं?

Question

मैं फोरट्रान और लैपैक का उपयोग करके एक वास्तविक सममित मैट्रिक्स को ट्रिडाइगोनलाइज करना चाहता हूं। LAPACK मूल रूप से दो दिनचर्या प्रदान करता है, एक पूर्ण मैट्रिक्स पर एक ऑपरेटिंग, दूसरा पैक किए गए भंडारण में मैट्रिक्स पर। जबकि उत्तरार्द्ध निश्चित रूप से कम स्मृति का उपयोग करता है, मैं सोच रहा था कि गति अंतर के बारे में कुछ भी कहा जा सकता है?

Answer 1

यह एक अनुभवजन्य प्रश्न है, लेकिन सामान्य रूप से, कुछ भी मुफ्त में नहीं आता है, और कम स्मृति/अधिक रनटाइम एक सुंदर आम व्यापार है।

इस मामले में, डेटा के लिए अनुक्रमण पैकिंग मामले के लिए अधिक जटिल है, इसलिए जब आप मैट्रिक्स को पार करते हैं, तो आपका डेटा प्राप्त करने की लागत थोड़ी अधिक होती है। (इस तस्वीर को जटिल बनाना यह है कि सममित मैट्रिस के लिए, लैपैक रूटीन भी एक निश्चित प्रकार की पैकिंग मानते हैं - कि आपके पास केवल मैट्रिक्स के ऊपरी या निचले घटक हैं)।

मैं आज पहले एक ईजिनप्रोबलेम के साथ गड़बड़ कर रहा था, इसलिए मैं इसे माप बेंचमार्क के रूप में उपयोग करूंगा; एक सरल सममित परीक्षण का मामला (Herdon मैट्रिक्स, http://people.sc.fsu.edu/~jburkardt/m_src/test_mat/test_mat.html से) के साथ की कोशिश कर रहा है, और साथ 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% अंतर है जो मैं मानता चाहिए ssyevd की तुलना में बड़ा मेरी अपेक्षा से (यह भी एक थोड़ा बड़ा के साथ है पैक किए गए मामले के लिए त्रुटि?)। यह इंटेल के एमकेएल के साथ है। प्रदर्शन अंतर सामान्य रूप से आपके मैट्रिक्स पर निर्भर करेगा, निश्चित रूप से, जैसा कि आप बता रहे हैं; आपको जो मैट्रिक्स करना है, उससे अधिक यादृच्छिक पहुंच, ओवरहेड खराब होगा। मेरे द्वारा उपयोग किया जाने वाला कोड निम्नानुसार है:

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