!&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&& 
!
!                            How to compile
!                           ===============
!
! ifort + MKL dynamic lp64:
! ---------------------------
!   ifort -o dsyerv_ifort_mkl.x dsyerv_check.F90  eispack.F -L/mnt/apps/intel/mkl/lib/intel64 -lmkl_lapack95_lp64 -lmkl_blas95_lp64 -lmkl_intel_lp64 -lmkl_intel_thread -lmkl_core -liomp5
!
! gfortran + MKL dynamic lp64 :
! -----------------------------
!   gfortran -o dsyerv_gfortran_mkl.x dsyerv_check.F90  eispack.F -L/mnt/apps/intel/mkl/lib/intel64 -lmkl_lapack95_lp64 -lmkl_blas95_lp64 -lmkl_intel_lp64 -lmkl_intel_thread -lmkl_core -liomp5
!
! pgf90 + MKL dynamic lp64 
! ------------------------
!   pgf90 -o dsyerv_pgf90_mkl.x  dsyerv_check.F90  eispack.F -L/mnt/apps/intel/mkl/lib/intel64 -lmkl_lapack95_lp64 -lmkl_blas95_lp64 -lmkl_intel_lp64 -lmkl_intel_thread -lmkl_core -liomp5
!
! ifort + MKL dynamic ilp64
! -------------------------
!   ifort -o dsyerv_ifort_i8_mkl.x -i8   dsyerv_check.F90  eispack.F -L/mnt/apps/intel/mkl/lib/intel64  -lmkl_lapack95_ilp64 -lmkl_blas95_ilp64 -lmkl_intel_ilp64 -lmkl_intel_thread -lmkl_core -liomp5
!
! gfortran + MKL dynamic ilp64 : 
! ------------------------------
!  gfortran -o dsyerv_gfortran_i8_mkl.x -fdefault-integer-8  dsyerv_check.F90  eispack.F -L/mnt/apps/intel/mkl/lib/intel64  -lmkl_lapack95_ilp64 -lmkl_blas95_ilp64 -lmkl_intel_ilp64 -lmkl_intel_thread -lmkl_core -liomp5
!
! ifort + dynamic GNU lp64 libs:
! ------------------------------
!   ifort -o dsyerv_ifort_GNUlibs.x   dsyerv_check.F90  eispack.F  -L/usr/lib64 -lblas -llapack
!
! gfortran + dynamic GNU libs lp64:
! ----------------------------------
!  gfortran -o dsyerv_gfortran_GNUlibs.x  dsyerv_check.F90  eispack.F  -L/usr/lib64 -lblas -llapack
!
! pgf90 +  GNU dynamic lp64 
! ------------------------
!   pgf90 -o dsyerv_pgf90_GNUlibs.x  dsyerv_check.F90 eispack.F  -L/usr/lib64 -lblas -llapack
!
!
!&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&& 
module diagmod
!-------------------------------------------------------------------------------------------
!  Module for the standalone testing program of diagonalization routines,
! utils/diag.F
!
! Written by Miro Ilias, August 2014
!-------------------------------------------------------------------------------------------
  integer :: N=0,NZ=0, print_level
  integer :: lu=12, LUPRI=6
  real*8, allocatable :: Ar(:,:), ASr(:,:),Ai(:,:),ASi(:,:),  & 
                         Aj(:,:),ASj(:,:), Ak(:,:), ASk(:,:)
  character*50 :: title_text
  ! hardcoded input matrix name
  character*50 :: matrix_file_name = "Jz_SS_matrix.fermirp2-2"
contains

  subroutine read_matrix
  ! allocates space and reads the matrix for the diagonalization testing from the file
    integer :: i,j, ii,jj
   
    open (lu,file=trim(matrix_file_name), access="sequential", status="old")  
    read(lu,*) N,NZ
    print *,'read N=',N,' NZ=',NZ

    ! allocates space for matrices
    call allocate_space
    ! read matrix elements from its open file
    do i=1,N
    do j=1,N
     if (NZ==1) then
       read(lu,*) ii,jj,Ar(ii,jj)
       ASr(ii,jj)=Ar(ii,jj)
     else if (NZ==2) then
       read(lu,*) ii,jj,Ar(ii,jj),Ai(ii,jj)
       ASr(ii,jj)=Ar(ii,jj); ASi(ii,jj)=Ai(ii,jj)
     else if (NZ==4) then
       read(lu,*) ii,jj,Ar(ii,jj),Ai(ii,jj),Aj(ii,jj),Ak(ii,jj)
       ASr(ii,jj)=Ar(ii,jj); ASi(ii,jj)=Ai(ii,jj)
       ASj(ii,jj)=Aj(ii,jj); ASk(ii,jj)=Ak(ii,jj)
     endif
    enddo
    enddo
    close(lu,status="keep")
    print *,"matrix for diagonalization was read"
 
  end subroutine read_matrix

  subroutine allocate_space
    if (N<=1.or.NZ<=0) then 
       print *,'allocate_space: N=',N,' NZ=',NZ
       stop 5
    endif
    if (NZ==1) then
      allocate( Ar(N,N), stat=ierr )
      allocate( ASr(N,N), stat=ierr )
    else if (NZ==2) then
      allocate( Ar(N,N), Ai(N,N), stat=ierr )
      allocate( ASr(N,N), ASi(N,N), stat=ierr )
    else if (NZ==4) then
      allocate( Ar(N,N), Ai(N,N), Aj(N,N), Ak(N,N), stat=ierr )
      allocate( ASr(N,N), ASi(N,N), ASj(N,N), ASk(N,N), stat=ierr )
    else 
       print *,'wrong value of NZ=',NZ
    endif
  end subroutine allocate_space

!&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
 subroutine test_lapack_dsyevr
 !
 ! tests lapack's routine DSYERV for real only symmetric matrices
 !
   implicit none
   integer :: ierr, MATZ, LLWORK, KLWORK, LILWORK
   integer :: i
   integer :: IL, INFO, IU, LDA, LDZ, LIWORK, LWORK, M, LRA, LRV
   character*1 :: JOBZ, RANGE, UPLO
   real*8 :: ABSTOL, VL, VU
   integer, allocatable :: ISUPPZ(:), IWORK(:)
   real*8, allocatable ::  WORK(:), EIG(:),VEC(:,:)
   real*8 :: DDUM, WORKDUM
   integer :: IWORKDUM,IDUM
   real*8, parameter :: D0=0.0D0

   ABSTOL=0.0d0; JOBZ='V'; RANGE = 'A'; UPLO = 'U'
   DDUM = D0; IDUM = 0; WORKDUM = D0; IWORKDUM = 0

   MATZ = 1 
   IF(MATZ.EQ.1) THEN
        JOBZ='V'
   ELSEIF(MATZ.EQ.0) THEN
        JOBZ='N'
   ELSE
     WRITE(*,*) 'Illegal value of MATZ: ',MATZ
     STOP 2
   ENDIF

   allocate(EIG(N),VEC(N,N), stat=ierr)
    
   ! 1.call to determine size 
   IERR =  0 
   ! FYI: LLWORK=-1; LIWORK=-1 on input means calculate temporary storage
   M = N ; LRA = N ;LRV = N
   CALL DSYEVR(JOBZ,'A','U',N,Ar,LRA,DDUM,DDUM,IDUM,IDUM,ABSTOL,  & 
               M,EIG,VEC,LRV,IDUM,EIG,-1,IWORKDUM,-1,IERR)

   LLWORK = NINT(EIG(1))
   LILWORK = IWORKDUM 

   allocate(WORK(LLWORK), stat=ierr)
   allocate(IWORK(LILWORK), stat=ierr)
   allocate(ISUPPZ(2*N),stat=ierr)

   IERR = 0
   CALL DSYEVR(JOBZ,'A','U',N,Ar,LRA,DDUM,DDUM,IDUM,IDUM,ABSTOL,  &
               M,EIG,VEC,LRV,ISUPPZ,WORK,LLWORK,IWORK,LILWORK,IERR)     
   if (IERR /= 0) then
     print *,'The lapack_dsyerv routine ended with error ! ierr=',ierr
   endif

  ! ... print out eigenvalues
  if (print_level >= 0) then
     print *,"LAPACK DSYEVR eigenvalues:"
     do i=1,N
       print *,i,EIG(i)
     enddo
   endif

   call eigv_check(EIG,VEC,'**** LAPACK DSYEVR ****')

   deallocate(WORK,IWORK,EIG,VEC,ISUPPZ,stat=ierr)

  end subroutine test_lapack_dsyevr

!&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
  subroutine test_dirac_rs
! tests DIRAC's EISPACK routine RS real Hermitian matrix
    real*8, allocatable :: eigval(:),eigvect(:,:),temp1(:),temp2(:)
    integer :: ierr,i

! ... allocate arrays for the subroutine
    allocate(eigval(N),eigvect(N,N),temp1(N),temp2(N),stat=ierr)
    print *,"arrays for RS subroutine allocated"

    matz=1 ! have both eigenvalues and eigenvectors
    call RS(N,N,Ar,eigval,matz,eigvect,temp1,temp2,ierr)
    if (ierr /= 0) then
      print *,'The RS routine ended with error !'
    endif
! ... print out eigenvalues
    if (print_level >= 0) then
      print *,"Eispack's RS eigenvalues:"
      do i=1,N
        print *,i,eigval(i)
      enddo
    endif

    call eigv_check(eigval,eigvect,'**** EISPACK RS ****')

    deallocate(eigval,eigvect,temp1,temp2,stat=ierr)
  end subroutine test_dirac_rs

  subroutine eigv_check(eigval,eigvect,routine_name)
  !-----------------------------------------------------------------------------------
  !
  ! Do the diagonalization checks using original Hermitian matrix, eigenvalues and
  ! eigenvectors
  ! 
  !           eigvect^+ * AS * eigvect - eigval  = 0
  !
  !           eigvect^+ * eigvect = I
  !
  !           eigvect * eigvect^+ = I
  !
  !-----------------------------------------------------------------------------------
   implicit none
    real*8, intent(in) :: eigval(N),eigvect(N,N)
    character(*), intent(in) :: routine_name
    real*8, allocatable :: AXr(:,:),AYr(:,:),C(:,:)
    integer :: i,j, ierr, iprint
    real*8 :: C_norm_diag, C_norm_offdiag
    real*8 :: AXr_norm_diag, AXr_norm_offdiag
    real*8 :: AYr_norm_diag, AYr_norm_offdiag
    integer(kind=8) :: maxsize = -1
    
    allocate( AXr(N,N),AYr(N,N),C(N,N),stat=ierr )
 
   ! Do  eigvect^{+)*ASr = AXr
    call dgemm('T','N',N,N,N,1.0d0,eigvect,N,ASr,N,0.0d0,AXr,N)

   ! Do  AXr * eigvect = C
    call dgemm('N','N',N,N,N,1.0d0,AXr,N,eigvect,N,0.0d0,C,N)

   ! Another check:  eigvect^{+)*eigvect = I = AXr
    call dgemm('T','N',N,N,N,1.0d0,eigvect,N,eigvect,N,0.0d0,AXr,N)
    if (print_level>=1) then
      write(LUPRI,*) "Printout of eigvect^{+)*eigvect =? I "
    endif

   ! Another check:  eigvect*eigvect^{+} = I = AYr
    call dgemm('N','T',N,N,N,1.0d0,eigvect,N,eigvect,N,0.0d0,AYr,N)
    if (print_level>=1) then
      write(LUPRI,*) "Printout of eigvect*eigvect^{+} =? I "
    endif

    ! substract eigenvalues ...
    C_norm_diag = 0.0d0; AXr_norm_diag = 0.0d0; AYr_norm_diag = 0.0d0
    do i=1,N
      C(i,i) = C(i,i) -eigval(i); 
      AXr(i,i)=AXr(i,i)-1.0d0; AYr(i,i)=AYr(i,i)-1.0d0
      C_norm_diag = C_norm_diag + DABS(C(i,i))
      AXr_norm_diag = AXr_norm_diag  + DABS(AXr(i,i))
      AYr_norm_diag = AYr_norm_diag  + DABS(AYr(i,i))
    enddo
    C_norm_diag   = C_norm_diag/DFLOAT(N)
    AXr_norm_diag = AXr_norm_diag/DFLOAT(N)
    AYr_norm_diag = AYr_norm_diag/DFLOAT(N)

    C_norm_offdiag = 0.0d0; AXr_norm_offdiag = 0.0d0;  AYr_norm_offdiag = 0.0d0
    do i = 1, N
    do j = 1, N
     if (i /= j) then
       C_norm_offdiag = C_norm_offdiag + DABS(C(i,j))
       AXr_norm_offdiag = AXr_norm_offdiag + DABS(AXr(i,j))
       AYr_norm_offdiag = AYr_norm_offdiag + DABS(AYr(i,j))
     endif
    enddo
    enddo
    C_norm_offdiag = C_norm_offdiag/((DFLOAT(N)*DFLOAT(N))-DFLOAT(N))
    AXr_norm_offdiag = AXr_norm_offdiag/((DFLOAT(N)*DFLOAT(N))-DFLOAT(N))
    AYr_norm_offdiag = AYr_norm_offdiag/((DFLOAT(N)*DFLOAT(N))-DFLOAT(N))

    write(LUPRI,"(2X,A)") routine_name
    write(LUPRI,"(1X,A,D10.4,2X,A,D10.4)") & 
 "U^{+}*A*U - eps ?= 0> norm/diag:",C_norm_diag,  "norm/offdiag:",C_norm_offdiag

    write(LUPRI,"(1X,A,D10.4,2X,A,D10.4)") & 
 "    U^{+}*U - I ?= 0> norm/diag:",AXr_norm_diag,"norm/offdiag:",AXr_norm_offdiag
    write(LUPRI,"(1X,A,D10.4,2X,A,D10.4)") & 
 "    U*U^{+} - I ?= 0> norm/diag:",AYr_norm_diag,"norm/offdiag:",  &
  AYr_norm_offdiag


    deallocate( AXr, AYr, C, stat=ierr )

  end subroutine eigv_check

end module diagmod


Program Test_DIRAC_Diagonalization_Routines
  use diagmod
   write(6,'(2x,a)') "Welcome to the testing program for DIRAC's diagonalization routines !"

  call read_matrix
  call test_dirac_rs
  call test_lapack_dsyevr

End Program Test_DIRAC_Diagonalization_Routines