program diagonalize

  implicit none

  integer, parameter :: dp = selected_real_kind(14)
  integer, parameter :: n=1000 ! Matrix dimension
  real(kind=dp), dimension(n,n) :: A
  real(kind=dp), dimension(n) :: lambda

  ! Helper arrays used for constructing random symmetric A with predefined 
  ! eigenvalues
  real(kind=dp), dimension(n) :: D ! Diagonal matrix
  real(kind=dp), dimension(2*n) :: work 
  integer, dimension(4) :: iseed=(/0, 1, 2, 3/)

  ! Work array for dsyev
  integer, parameter :: lwork = 3*n - 1
  real(kind=dp), dimension(lwork) :: work2

  real(kind=dp) :: t0, t1 ! Timing variables

  integer :: i, info
  
  external dlagsy, dsyev

  ! Initialize matrices
  do i=1,n
        D(i) = i
  end do

  call dlagsy(n, n-1, D, A, n, iseed, work, info)

  ! Lapack routine
  call dsyev("V", "U", n, A, n, lambda, work2, lwork, info)

! Check the results
  do i=1,n
     if (ABS(lambda(i) - D(i)) > 0.001) then
        write(*,*) 'Error in diagonalization'
     end if
  end do


end program diagonalize

