subroutine calculate_EVD           &
        (NP,NTOTAL,LDZ,AP_single,W_reverse,Z_reverse)
!------------------------------
!  Eigenvalue decomposition
!
!  Employs the LAPACK programs DSPEV  (double precision), SGEMM (single precision)
!-----------------------------------------------------------------------------
      integer :: NP,NTOTAL,LDZ
      real AP_single(1:NP)                        !! INPUT : Matrix A
      real*8 W_reverse(1:NTOTAL)                  !! OUTPUT: Eigenvalues
      real*8 Z_reverse(1:LDZ,1:NTOTAL)            !! OUTPUT: Eigenvectors
!-----------------------------------------------------------------------------

      real*8,dimension(:),allocatable :: AP
      real*8,dimension(:),allocatable :: W
      real*8,dimension(:,:),allocatable :: Z
      real*8,dimension(:),allocatable :: AUX
      real*8,dimension(:),allocatable :: q1,q2,q3,q4

      real :: sum11,sum12,sum13,sum14
      real :: sum22,sum23,sum24
      real :: sum33,sum34
      real :: sum44
      real :: sumx1,sumx2
      real :: sumy1,sumy2
      real :: sumxx,sumyy,sumxy
      integer :: INFO

      EXTERNAL DSPEV
!====================================================================
!-- define dimensions -----------
!
       NAUX=3*NTOTAL
!----------------------------------------
         allocate(AP(1:NP))
         allocate(W(1:NTOTAL))
         allocate(Z(1:LDZ,1:NTOTAL))
         allocate(AUX(1:NAUX))

         DO K=1,NP
          AP(K)=AP_single(K)
           END DO
!====================================================================
          write(*,*) "------------ "    
          write(*,*) "INPUT MATRIX "    
          write(*,*) "------------ "    
           DO K=1,min(10,NP)
            write(*,*) "AP(",K,")=",AP(K)
             END DO

!---- EIGEN-VALUE DECOMPOSITION ------------------
!---------------------------------------------
!!!---------- LAPACK -----------------------------------
!   'V' - eigenvectors and eigenvalues
!   'N' - eigenvalues only
!
!   'U' - upper triangle of A stored
!   'L' - lower triangle of A stored

       CALL DSPEV ('V','U',NTOTAL,AP,W,Z,LDZ,AUX,INFO)
!---------------------------------------------

          write(*,*) "------------ "    
          write(*,*) "OUTPUT EIGENVALUES "    
          write(*,*) "------------ "    
           DO K=1,NTOTAL
           write(*,*) "ORIGINAL W(",K,")=",W(K)
           END DO

      allocate(q1(1:NTOTAL))
      allocate(q2(1:NTOTAL))
      allocate(q3(1:NTOTAL))
      allocate(q4(1:NTOTAL))

          write(*,*) "------------ "    
          write(*,*) "OUTPUT EIGENVECTORS"    
          write(*,*) "------------ "    
           DO j=1,NTOTAL
          q1(j)=Z(j,1)
          q2(j)=Z(j,2)
          q3(j)=Z(j,3)
          q4(j)=Z(j,4)
           END DO

          sum11=0.
          sum12=0.
          sum13=0.
          sum14=0.
          sum22=0.
          sum23=0.
          sum24=0.
          sum33=0.
          sum34=0.
          sum44=0.
           do i=1,NTOTAL
           sum11=sum11+q1(i)*q1(i)
           sum12=sum12+q1(i)*q2(i)
           sum13=sum13+q1(i)*q3(i)
           sum14=sum14+q1(i)*q4(i)
           sum22=sum22+q2(i)*q2(i)
           sum23=sum23+q2(i)*q3(i)
           sum24=sum24+q2(i)*q4(i)
           sum33=sum33+q3(i)*q3(i)
           sum34=sum34+q3(i)*q4(i)
           sum44=sum44+q4(i)*q4(i)
           end do

!-------------------------------------------
           DO j=1,NTOTAL
          q1(j)=Z(j,1)
          q2(j)=Z(j,2)
          q3(j)=Z(j,NTOTAL-1)
          q4(j)=Z(j,NTOTAL)
           END DO
                                                                                                   
          sumxx=0.
          sumyy=0.
          sumxy=0.
          sumx1=0.
          sumx2=0.
          sumy1=0.
          sumy2=0.
           do i=1,NTOTAL
           sumxx=sumxx+q4(i)*q4(i)
           sumyy=sumyy+q3(i)*q3(i)
           sumxy=sumxy+q4(i)*q3(i)
           sumx1=sumx1+q4(i)*q1(i)
           sumy1=sumy1+q3(i)*q1(i)
           sumx2=sumx2+q4(i)*q2(i)
           sumy2=sumy2+q3(i)*q2(i)
           end do


          write(*,*) "sum11=",sum11
          write(*,*) "sum12=",sum12
          write(*,*) "sum13=",sum13
          write(*,*) "sum14=",sum14
          write(*,*) "sum22=",sum22
          write(*,*) "sum23=",sum23
          write(*,*) "sum24=",sum24
          write(*,*) "sum33=",sum33
          write(*,*) "sum34=",sum34
          write(*,*) "sum44=",sum44
          write(*,*) "sumxx=",sumxx
          write(*,*) "sumyy=",sumyy
          write(*,*) "sumxy=",sumxy
          write(*,*) "sumx1=",sumx1
          write(*,*) "sumx2=",sumx2
          write(*,*) "sumy1=",sumy1
          write(*,*) "sumy2=",sumy2

          write(*,100) (q1(I),I=1,5)
 100      FORMAT('Q1=',5F12.6)
          write(*,200) (q2(I),I=1,5)
 200      FORMAT('Q2=',5F12.6)
          write(*,300) (q3(I),I=1,5)
 300      FORMAT('Q3=',5F12.6)
          write(*,400) (q4(I),I=1,5)
 400      FORMAT('Q4=',5F12.6)

      deallocate(q1)
       deallocate(q2)
        deallocate(q3)
         deallocate(q4)
          deallocate(AP)
           deallocate(AUX)

!--- reverse the order of eigenvalues (from largest to smallest, to accomodate for I+A inverse) ---

           do j=1,NTOTAL
            W_reverse(j)=W(NTOTAL-j+1)
             do i=1,LDZ
            Z_reverse(i,j)=Z(i,NTOTAL-j+1)
           end do
          end do

        deallocate(W)
      deallocate(Z)


           END SUBROUTINE calculate_EVD