program main
  implicit none
  integer :: n
  real(8), allocatable :: A(:,:)
  real(8), allocatable :: U(:,:)
  real(8), allocatable :: w(:)
  real(8), allocatable :: L(:,:)
  real(8), allocatable :: Uinv(:,:)
  real(8), allocatable :: B(:,:)
  integer :: m
  integer :: info
  integer :: i
  character(64) :: wfmt

  n = 4
  allocate(A(n,n))
  allocate(w(n))
  allocate(U(n,n))
  A(1,:) = (/1,2,3,4/)
  A(2,:) = (/2,2,3,4/)
  A(3,:) = (/3,3,3,4/)
  A(4,:) = (/4,4,4,4/)

  call eig_sym(A, m, w, U, info)

  print*, 'info', info

  write(wfmt,"('(',i0,'(1x,f10.7))')") m
  print*, 'Eigenvalues'
  print wfmt, w(1:m)

  print*, 'Eigenvectors'
  do i = 1, m
    print wfmt, U(i,:)
  enddo

  allocate(B(n,n))
  allocate(L(n,n))
  allocate(Uinv(n,n))
  L = diag(w)
  call inv(U, Uinv, info)
  B = U .mul. L .mul. Uinv

  print*, 'ULU-1'
  do i = 1, n
    print wfmt, B(i,:)
  enddo
end program main