!=======================================================================
! ch08_jacobi_eigen.f90
!
! Find all eigenvalues and eigenvectors of a real symmetric matrix.
!
! Origin:
!   Based on Jacobi.f90 by Alexander Godunov, December 2009.
!   The cyclic, thresholded Jacobi rotations and the original formulas
!   for c and s are retained. This version limits the number of sweeps
!   and checks the resulting eigenpairs against the original matrix.
!
! Example:
!   A = [ 1  2  3; 2  2 -2; 3 -2  4 ]. Its eigenvalues are approximately
!   -2.54138127, 3.54138127, and 6.00000000.
!
! The columns of X are the eigenvectors associated with the diagonal
! entries of the transformed matrix, in the order produced by Jacobi.
! Eigenvectors may differ by an overall sign. This method is useful
! for a small symmetric matrix; library routines are preferable for
! larger problems.
!=======================================================================
program ch08_jacobi_eigen
implicit none

integer, parameter :: n = 3, max_sweeps = 100
! As in Jacobi.f90, abserr bounds the SUM OF SQUARED off-diagonal
! entries, not the norm itself. The smaller value supports reporting
! eigenvectors and residuals to useful precision.
double precision, parameter :: abserr = 1.0d-18
double precision :: a(n,n), a_work(n,n), x(n,n)
double precision :: offdiag_norm, max_residual, max_orthogonality
double precision :: value, defect
integer :: i, j, k, sweeps, rotations, info

! Each implied-DO fills one row of the symmetric example matrix.
data (a(1,j),j=1,n) / 1.0d0,  2.0d0,  3.0d0 /
data (a(2,j),j=1,n) / 2.0d0,  2.0d0, -2.0d0 /
data (a(3,j),j=1,n) / 3.0d0, -2.0d0,  4.0d0 /

write (*,'(a)') 'Eigenvalues and eigenvectors (Jacobi method)'
write (*,'(/,a)') 'Original symmetric matrix A:'
do i = 1,n
   write (*,'(3f12.6)') (a(i,j),j=1,n)
end do
write (*,'(a,es12.4)') 'Squared off-diagonal tolerance = ',abserr

! Jacobi overwrites its matrix with a nearly diagonal similar matrix.
! Preserve A so that the eigenpairs can be checked independently.
a_work = a
call jacobi(a_work,x,n,abserr,max_sweeps,sweeps,rotations, &
            offdiag_norm,info)
if (info == 1) then
   write (*,'(/,a)') 'Invalid input: A must be symmetric and controls positive.'
   stop 1
end if

write (*,'(/,a)') 'Eigenvalues and corresponding unit eigenvectors:'
do j = 1,n
   write (*,'(a,i0,a,f12.6)') 'lambda(',j,') = ',a_work(j,j)
   do i = 1,n
      write (*,'(a,i0,a,f12.6)') '  v(',i,') = ',x(i,j)
   end do
end do

! Check every component of A*X-X*D using the original matrix.
! The maximum absolute component is easy to compare with the input
! scale. For orthogonality, compare every column dot product with I.
max_residual = 0.0d0
max_orthogonality = 0.0d0
do j = 1,n
   do i = 1,n
      value = 0.0d0
      do k = 1,n
         value = value + a(i,k)*x(k,j)
      end do
      max_residual = max(max_residual,abs(value-a_work(j,j)*x(i,j)))

      value = 0.0d0
      do k = 1,n
         value = value + x(k,i)*x(k,j)
      end do
      defect = 0.0d0
      if (i == j) defect = 1.0d0
      max_orthogonality = max(max_orthogonality,abs(value-defect))
   end do
end do

write (*,'(/,a,i0)') 'Sweeps = ',sweeps
write (*,'(a,i0)') 'Rotations = ',rotations
write (*,'(a,es12.4)') 'Off-diagonal norm of final matrix = ',offdiag_norm
write (*,'(a,es12.4)') 'Maximum |A*X-X*D| component = ',max_residual
write (*,'(a,es12.4)') 'Maximum |X^T*X-I| component = ',max_orthogonality
if (info == 2) then
   write (*,'(a)') 'Sweep limit reached before the tolerance was met.'
   stop 1
end if

contains

subroutine jacobi(a,x,n,abserr,max_sweeps,sweeps,rotations,offdiag_norm,info)
!=======================================================================
! Diagonalize a real symmetric matrix by thresholded Jacobi rotations.
!
! Input/output:
!   a(n,n)       - symmetric matrix; overwritten by a nearly diagonal
!                  similar matrix whose diagonal holds the eigenvalues
! Output:
!   x(n,n)       - orthonormal eigenvectors in columns: A_original*x_j
!                  is approximately a(j,j)*x_j
!   sweeps       - completed passes over all off-diagonal pairs
!   rotations    - number of rotations actually performed
!   offdiag_norm - sqrt of the sum of all off-diagonal squares
!   info         - 0 converged; 1 invalid input; 2 sweep limit
! Controls:
!   abserr       - positive bound on SUM of squared off-diagonal entries
!   max_sweeps   - positive maximum number of complete passes
!
! Start X as the identity. For each off-diagonal pair above the
! threshold, rotate rows and columns of A by the same plane angle;
! rotate columns of X as well. Orthogonal similarity preserves the
! eigenvalues. Recompute the off-diagonal sum after each sweep rather
! than subtracting an idealized change that can accumulate roundoff.
!
! The c and s formulas follow Jacobi.f90, including the near-90-degree
! rotation it may choose when two diagonal entries differ greatly.
! This code assumes a small, exactly symmetric input matrix.
!=======================================================================
implicit none
integer, intent(in) :: n, max_sweeps
double precision, intent(in) :: abserr
double precision, intent(inout) :: a(n,n)
double precision, intent(out) :: x(n,n), offdiag_norm
integer, intent(out) :: sweeps, rotations, info

double precision :: b2, bar, beta, coeff, c, s, cs, sc
integer :: i, j, k

x = 0.0d0
do i = 1,n
   x(i,i) = 1.0d0
end do
sweeps = 0
rotations = 0
offdiag_norm = huge(1.0d0)
info = 1
if (n < 1 .or. abserr <= 0.0d0 .or. max_sweeps < 1) return
do i = 1,n
   do j = i+1,n
      if (a(i,j) /= a(j,i)) return
   end do
end do

do
   ! Measure the actual off-diagonal entries at the start of each
   ! sweep. In exact arithmetic b2 decreases after every rotation.
   b2 = 0.0d0
   do i = 1,n
      do j = 1,n
         if (i /= j) b2 = b2 + a(i,j)**2
      end do
   end do
   offdiag_norm = sqrt(b2)
   if (b2 <= abserr) then
      info = 0
      return
   end if
   if (sweeps >= max_sweeps) then
      info = 2
      return
   end if

   ! Original code sets bar to half the mean off-diagonal square.
   ! As significant elements disappear, lower the threshold during
   ! this sweep using the expected reduction of b2.
   bar = 0.5d0*b2/dble(n*n)
   do i = 1,n-1
      do j = i+1,n
         if (a(j,i)**2 <= bar) cycle
         b2 = max(0.0d0,b2-2.0d0*a(j,i)**2)
         bar = 0.5d0*b2/dble(n*n)

         ! Choose the same c and s as in the original Jacobi.f90.
         ! The two rotations below eliminate a(j,i), up to roundoff.
         beta = (a(j,j)-a(i,i))/(2.0d0*a(j,i))
         coeff = 0.5d0*beta/sqrt(1.0d0+beta**2)
         s = sqrt(max(0.5d0+coeff,0.0d0))
         c = sqrt(max(0.5d0-coeff,0.0d0))

         ! First rotate rows i and j of A.
         do k = 1,n
            cs =  c*a(i,k)+s*a(j,k)
            sc = -s*a(i,k)+c*a(j,k)
            a(i,k) = cs
            a(j,k) = sc
         end do

         ! Then rotate columns i and j of A. Apply the same column
         ! rotation to X so its columns become the eigenvectors.
         do k = 1,n
            cs =  c*a(k,i)+s*a(k,j)
            sc = -s*a(k,i)+c*a(k,j)
            a(k,i) = cs
            a(k,j) = sc

            cs =  c*x(k,i)+s*x(k,j)
            sc = -s*x(k,i)+c*x(k,j)
            x(k,i) = cs
            x(k,j) = sc
         end do
         rotations = rotations + 1
      end do
   end do
   sweeps = sweeps + 1
end do
end subroutine jacobi

end program ch08_jacobi_eigen