symmetric_eigendecomposition Subroutine

subroutine symmetric_eigendecomposition(amat, eval, evec)

   !> Symmetric matrix to decompose
   real(wp), intent(in) :: amat(:, :)

   !> Eigenvalues, sorted by descending magnitude
   real(wp), intent(out) :: eval(:)

   !> Eigenvectors, stored column-wise in the same order as the eigenvalues
   real(wp), intent(out) :: evec(:, :)

   integer :: ndim, isweep, ip, iq, jdim, imax
   real(wp) :: anorm, off, theta, tval, cval, sval, skip
   real(wp) :: ajp, ajq, apj, aqj, vjp, vjq
   real(wp), allocatable :: work(:, :)

   ndim = size(amat, 1)

   evec(:, :) = 0.0_wp
   do ip = 1, ndim
      evec(ip, ip) = 1.0_wp
   end do
   eval(:) = 0.0_wp

   anorm = sqrt(sum(amat**2))
   if (anorm <= 0.0_wp) return

   skip = 0.01_wp * epsilon(1.0_wp) * anorm
   work = amat

   do isweep = 1, max_sweeps
      off = 0.0_wp
      do ip = 1, ndim - 1
         off = off + sum(work(ip, ip+1:)**2)
      end do
      if (sqrt(2.0_wp*off) <= epsilon(1.0_wp) * anorm) exit

      do ip = 1, ndim - 1
         do iq = ip + 1, ndim
            if (abs(work(ip, iq)) <= skip) cycle

            theta = 0.5_wp * (work(iq, iq) - work(ip, ip)) / work(ip, iq)
            tval = sign(1.0_wp, theta) / (abs(theta) + sqrt(1.0_wp + theta*theta))
            cval = 1.0_wp / sqrt(1.0_wp + tval*tval)
            sval = tval * cval

            do jdim = 1, ndim
               ajp = work(jdim, ip)
               ajq = work(jdim, iq)
               work(jdim, ip) = cval*ajp - sval*ajq
               work(jdim, iq) = sval*ajp + cval*ajq
            end do
            do jdim = 1, ndim
               apj = work(ip, jdim)
               aqj = work(iq, jdim)
               work(ip, jdim) = cval*apj - sval*aqj
               work(iq, jdim) = sval*apj + cval*aqj
            end do
            work(ip, iq) = 0.0_wp
            work(iq, ip) = 0.0_wp

            do jdim = 1, ndim
               vjp = evec(jdim, ip)
               vjq = evec(jdim, iq)
               evec(jdim, ip) = cval*vjp - sval*vjq
               evec(jdim, iq) = sval*vjp + cval*vjq
            end do
         end do
      end do
   end do

   do ip = 1, ndim
      eval(ip) = work(ip, ip)
   end do

   do ip = 1, ndim - 1
      imax = ip - 1 + maxloc(abs(eval(ip:)), 1)
      if (imax == ip) cycle
      eval([ip, imax]) = eval([imax, ip])
      evec(:, [ip, imax]) = evec(:, [imax, ip])
   end do

end subroutine symmetric_eigendecomposition