
!===================!
 module fp_precision
!===================!
    integer, parameter :: dp = kind(1.0d0)
!=======================!
 end module fp_precision
!=======================!



!===================!
 module print_matrix
!===================!
    use fp_precision
    implicit none
!--------!
 contains
!--------!

    subroutine matprt (nl, nc, mat2, ldmat2, cmat2, mat1, cmat1)
    
    integer, intent(in) :: ldmat2, nc, nl

    real(dp), intent(in)           :: mat2(ldmat2,nc)
    real(dp), intent(in), optional :: mat1(nc)

    character(*), intent(in)           :: cmat2
    character(*), intent(in), optional :: cmat1

! Local data

    integer :: i, ilast, inext, j

    logical :: prt1


    if (present(mat1)) then
       prt1 = .true.
    else
       prt1 = .false.
    end if

    write (6,*) 
    inext = 1
    do i = 1, nc/10
       ilast = inext+9
       write (6,3) (j, j=inext,ilast)
       if (prt1) then
          write (6,1) cmat1, mat1(inext:ilast)
          write (6,*)
       end if
       write (6,1) cmat2, mat2(1,inext:ilast)
       do j = 2, nl
          write (6,2) j, mat2(j,inext:ilast)
       end do
       write (6,*)
       inext = ilast + 1
    end do

    if (inext <= nc) then
       write (6,3) (j, j=inext,nc)
       if (prt1) then
          write (6,1) cmat1, mat1(inext:nc)
          write (6,*)
       end if
       write (6,1) cmat2, mat2(1,inext:nc)
       do j = 2, nl
          write (6,2) j, mat2(j,inext:nc)
       end do
       write (6,*)
    end if

    return

  1 format (1x, a8, 2x, 10f15.7)
  2 format (3x, i5, 3x, 10f15.7)
!  1 format (1x, a8, 2x, 10e15.7)
!  2 format (3x, i5, 3x, 10e15.7)
  3 format (9x, 10i15)

    end subroutine matprt
!=======================!
 end module print_matrix
!=======================!

!==================!
 program test_dsbev
!==================!
    use fp_precision
    use print_matrix
    
    implicit none

    integer :: count, i, info, j, kd, n
    integer :: largest_component(1)

    real(dp), allocatable :: eval(:)
    real(dp), allocatable :: work(:)

    real(dp), allocatable :: a(:,:)
    real(dp), allocatable :: aband(:,:)
    real(dp), allocatable :: evec(:,:)
    real(dp), allocatable :: temp(:,:)

! read banded matrix to be diagonalised 

    read (5,*) n, kd
    write (6,'("Matrix dimension = ",i3, ",   number of bands = ", i3)') n, kd

    allocate (aband(kd+1,n)); aband = 0.0_dp

    do i = 1, kd+1
       read (5,*) aband(i,1:n)
    end do

    call matprt (kd+1, n, aband, kd+1, "Matrix")

	allocate (a(n,n));    a = 0.0_dp
    allocate (eval(n));   eval = 0.0_dp
    allocate (evec(n,n)); evec = 0.0_dp
	allocate (temp(n,n));    temp = 0.0_dp
    allocate (work(3*n)); work = 0.0_dp

    do j = 1, n
       do i = j, min(n,j+kd)
          a(i,j) = aband(1+i-j,j)
          a(j,i) = aband(1+i-j,j)
       end do
    end do
    
! diagonalisation using Lapack

    call dsbev ('v', 'l', n, kd, aband, kd+1, eval, evec, n, work, info)

! phase such that the largest component of each vector is positive

    do j = 1, n
       largest_component = maxloc (abs(evec(:,j)))
       if (evec(largest_component(1),j) < 0.0_dp) evec(:,j) = -evec(:,j)
    end do

! print eveigenvalues and associated eigenvectors

    call matprt (n, n, evec, n, "evec", eval, "eval")

! check orthonormality

	temp = matmul (transpose(evec), evec)
	do j = 1,n
	   temp(j,j) = temp(j,j) - 1.0_dp
	end do
	
	write (6,'(" Orthonormality: largest non-zero element is ",e15.7," at ",2i3)') maxval(abs(temp)), maxloc(abs(temp))

! check eigenvectors

    temp = matmul (a, evec)
    a = matmul (transpose(evec), temp)
    
    do j = 1, n
       a(j,j) = a(j,j) - eval(j)
    end do
    
    call matprt (n, n, a, n, "check")
    
	write (6,'(" Eigenvectors: largest non-zero element is ",e15.7," at ",2i3)') maxval(abs(a)), maxloc(abs(a))
	   
! clean up

    deallocate (aband, eval, evec, work)

    stop
!======================!
 end program test_dsbev
!======================!

