!=======================================================================
! ch08_power_method.f90
!
! Find the dominant eigenvalue and eigenvector by the power method.
!
! Origin:
!   Based on Power.f90 by Alexander Godunov, December 2009.
!   The iteration y <- A*y / ||A*y|| is retained. The eigenvalue
!   estimate now uses the signed Rayleigh quotient from Chapter 8,
!   and the eigenpair residual is used to test convergence.
!
! Example:
!   A = [ 1  2  3; 2  2 -2; 3 -2  4 ], y0 = [ 1, 1, 1 ].
!   The eigenvalues are approximately -2.54138127, 3.54138127,
!   and 6.00000000. The dominant eigenvalue in magnitude is 6.
!
! The method finds one eigenpair. It generally requires a unique
! dominant eigenvalue in magnitude and an initial vector with a
! nonzero component along its eigenvector. A small change in the
! eigenvalue estimate alone does not verify the eigenvector.
!=======================================================================
program ch08_power_method
implicit none

integer, parameter :: n = 3, max_iter = 1000
double precision, parameter :: eps = 1.0d-7
double precision :: a(n,n), y(n), lambda, residual
integer :: i, j, iterations, 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 /
data (y(i),i=1,n)   / 1.0d0,  1.0d0,  1.0d0 /

write (*,'(a)') 'Power method for the dominant eigenpair'
write (*,'(/,a)') 'Matrix A:'
do i = 1,n
   write (*,'(3f12.6)') (a(i,j),j=1,n)
end do
write (*,'(/,a,3f12.6)') 'Initial vector:',y
write (*,'(a,es12.4)') 'Residual tolerance = ',eps

call power_method(a,y,n,eps,max_iter,lambda,iterations,residual,info)

if (info == 1) then
   write (*,'(/,a)') 'Invalid input: check the initial vector and controls.'
   stop 1
end if
if (info == 2) then
   write (*,'(/,a)') 'A*y is zero; the iteration cannot be normalized.'
   stop 1
end if

write (*,'(/,a,f12.6)') 'Dominant eigenvalue = ',lambda
write (*,'(a)') 'Unit eigenvector:'
do i = 1,n
   write (*,'(a,i0,a,f12.6)') 'y(',i,') = ',y(i)
end do
write (*,'(a,i0)') 'Iterations = ',iterations
write (*,'(a,es12.4)') 'Eigenpair residual ||A*y-lambda*y|| = ',residual
if (info == 3) then
   write (*,'(a)') 'Iteration limit reached before the tolerance was met.'
   stop 1
end if

contains

subroutine power_method(a,y,n,eps,max_iter,lambda,iterations,residual,info)
!=======================================================================
! Approximate the eigenpair of largest eigenvalue magnitude.
!
! Input:
!   a(n,n)     - real coefficient matrix, left unchanged
!   y(n)       - nonzero initial vector, overwritten by a unit vector
!   n          - matrix dimension (n >= 1)
!   eps        - positive tolerance for the eigenpair residual
!   max_iter   - maximum number of normalized power iterations
!
! Output:
!   lambda     - signed Rayleigh quotient y^T*A*y / (y^T*y)
!   iterations - completed power iterations
!   residual   - Euclidean norm of A*y - lambda*y
!   info       - 0 if residual <= eps; 1 for invalid input;
!                2 if A*y is zero; 3 if the iteration limit is reached
!
! Method:
!   Normalize the initial vector. Repeatedly multiply by A and
!   normalize, as in Power.f90. Compute the Rayleigh quotient and
!   residual for the resulting vector. The quotient retains the sign
!   of a negative dominant eigenvalue; ||A*y|| alone cannot do so.
!
! A small residual verifies that the returned pair is an eigenpair.
! It cannot by itself establish that the pair is the dominant one
! if the initial vector misses the dominant eigendirection.
!=======================================================================
implicit none
integer, intent(in) :: n, max_iter
double precision, intent(in) :: a(n,n), eps
double precision, intent(inout) :: y(n)
double precision, intent(out) :: lambda, residual
integer, intent(out) :: iterations, info

double precision :: z(n), ay(n), norm_y, norm_z, residual_sq
integer :: i, k

lambda = 0.0d0
residual = huge(1.0d0)
iterations = 0
info = 1
if (n < 1 .or. eps <= 0.0d0 .or. max_iter < 1) return

! An eigenvector has no preferred length, so start with ||y|| = 1.
norm_y = 0.0d0
do i = 1,n
   norm_y = norm_y + y(i)*y(i)
end do
norm_y = sqrt(norm_y)
if (norm_y == 0.0d0) return
do i = 1,n
   y(i) = y(i)/norm_y
end do

do k = 1,max_iter
   ! Power step: multiply by A, then normalize to prevent growth.
   call matrix_vector(a,y,z,n)
   norm_z = 0.0d0
   do i = 1,n
      norm_z = norm_z + z(i)*z(i)
   end do
   norm_z = sqrt(norm_z)
   if (norm_z == 0.0d0) then
      info = 2
      return
   end if
   do i = 1,n
      y(i) = z(i)/norm_z
   end do

   ! Estimate the signed eigenvalue for this normalized vector.
   call matrix_vector(a,y,ay,n)
   lambda = 0.0d0
   do i = 1,n
      lambda = lambda + y(i)*ay(i)
   end do

   ! Check the full eigenpair, not only the change in eigenvalue.
   residual_sq = 0.0d0
   do i = 1,n
      residual_sq = residual_sq + (ay(i)-lambda*y(i))**2
   end do
   residual = sqrt(residual_sq)
   iterations = k
   if (residual <= eps) then
      info = 0
      return
   end if
end do

info = 3
end subroutine power_method

subroutine matrix_vector(a,v,w,n)
!-----------------------------------------------------------------------
! Compute w=A*v explicitly, without calling a matrix library.
!-----------------------------------------------------------------------
implicit none
integer, intent(in) :: n
double precision, intent(in) :: a(n,n), v(n)
double precision, intent(out) :: w(n)
integer :: i, j

do i = 1,n
   w(i) = 0.0d0
   do j = 1,n
      w(i) = w(i) + a(i,j)*v(j)
   end do
end do
end subroutine matrix_vector

end program ch08_power_method