program main

!*****************************************************************************80
!
!! power_method_test() tests power_method().
!
!  Licensing:
!
!    This code is distributed under the MIT license.
!
!  Modified:
!
!    25 May 2008
!
!  Author:
!
!    John Burkardt
!
  implicit none

  integer, parameter :: rk = kind ( 1.0D+00 )

  call timestamp ( )
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'POWER_METHOD_TEST'
  write ( *, '(a)' ) '  FORTRAN90 version:'
  write ( *, '(a)' ) '  Test the POWER_METHOD library.'

  call test01 ( )
  call test02 ( )
  call test03 ( )
!
!  Terminate.
!
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'POWER_METHOD_TEST:'
  write ( *, '(a)' ) '  Normal end of execution.'
  write ( *, '(a)' ) ' '
  call timestamp ( )

  stop
end
subroutine test01 ( )

!*****************************************************************************80
!
!! TEST01 uses POWER_METHOD on the Fibonacci2 matrix.
!
!  Discussion:
!
!    This matrix, despite having a single dominant eigenvalue, will generally
!    converge only very slowly under the power method.  This has to do with
!    the fact that the matrix has only 3 eigenvectors.
!
!  Licensing:
!
!    This code is distributed under the MIT license.
!
!  Modified:
!
!    20 July 2008
!
!  Author:
!
!    John Burkardt
!
  implicit none

  integer, parameter :: rk = kind ( 1.0D+00 )

  integer, parameter :: n = 100

  real ( kind = rk ) a(n,n)
  real ( kind = rk ) cos_x1x2
  real ( kind = rk ) ctime
  real ( kind = rk ) ctime1
  real ( kind = rk ) ctime2
  integer i
  integer it_max
  integer it_num
  real ( kind = rk ) lambda
  real ( kind = rk ) phi
  real ( kind = rk ) sin_x1x2
  real ( kind = rk ) tol
  real ( kind = rk ) x(n)
  real ( kind = rk ) x2(n)

  call fibonacci2 ( n, a )

  call random_number ( harvest = x(1:n) )

  it_max = 500
  tol = 0.00000001D+0

  phi = ( 1.0 + sqrt ( 5.0 ) ) / 2.0

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST01'
  write ( *, '(a)' ) '  Use POWER_METHOD on the Fibonacci2 matrix.'
  write ( *, '(a)' ) ' '
  write ( *, '(a,i8)'      ) '  Matrix order N         = ', n
  write ( *, '(a,i8)'      ) '  Maximum iterations     = ', it_max
  write ( *, '(a,g14.6)'   ) '  Error tolerance        = ', tol

  call cpu_time ( ctime1 )

  call power_method ( n, a, x, it_max, tol, lambda, it_num )

  call cpu_time ( ctime2 )
  ctime = ctime2 - ctime1

  write ( *, '(a)' ) ' '
  write ( *, '(a,i8)'      ) '  Number of iterations   = ', it_num
  write ( *, '(a,g14.6)'   ) '  CPU time               = ', ctime
  write ( *, '(a,f14.10)'  ) '  Estimated eigenvalue   = ', lambda
  write ( *, '(a,f14.10)'  ) '  Correct value          = ', phi
  write ( *, '(a,g14.6)'   ) '  ||Error||              = ', abs ( lambda - phi )
!
!  X2 is the exact eigenvector.
!  Computing it "backwards" this way avoids overflow when N is large.
!
  x2(n) = 1.0D+00
  do i = n, 2, -1
    x2(i-1) = x2(i) / phi
  end do

  x2(1:n) = x2(1:n) / sqrt ( sum ( x2(1:n)**2 ) )
!
!  The sine of the angle between X and X2 is a measure of error.
!
  cos_x1x2 = dot_product ( x(1:n), x2(1:n) )
  sin_x1x2 = sqrt ( ( 1.0 - cos_x1x2 ) * ( 1.0 + cos_x1x2 ) )

  write ( *, '(a)' ) ' '
  write ( *, '(a,f14.10)' ) &
    '  Sine of angle between true and estimated vectors = ', sin_x1x2

  return
end
subroutine test02 ( )

!*****************************************************************************80
!
!! TEST02 uses POWER_METHOD2 on the Fibonacci2 matrix.
!
!  Discussion:
!
!    This matrix, despite having a single dominant eigenvalue, will generally
!    converge only very slowly under the power method.  This has to do with
!    the fact that the matrix has only 3 eigenvectors.
!
!  Licensing:
!
!    This code is distributed under the MIT license.
!
!  Modified:
!
!    20 July 2008
!
!  Author:
!
!    John Burkardt
!
  implicit none

  integer, parameter :: rk = kind ( 1.0D+00 )

  integer, parameter :: ck = kind ( ( 1.0D+00, 1.0D+00 ) )

  integer, parameter :: n = 10

  real ( kind = rk ) a(n,n)
  real ( kind = rk ) cos_x1x2
  real ( kind = rk ) ctime
  real ( kind = rk ) ctime1
  real ( kind = rk ) ctime2
  integer i
  integer it_max
  integer it_num
  complex ( kind = ck ) lambda
  real ( kind = rk ) phi
  real ( kind = rk ) sin_x1x2
  real ( kind = rk ) tol
  complex ( kind = ck ) v(n)
  real ( kind = rk ) x(n)
  real ( kind = rk ) x2(n)

  call fibonacci2 ( n, a )

  call random_number ( harvest = x(1:n) )

  it_max = 500
  tol = 0.00000001D+0

  phi = ( 1.0 + sqrt ( 5.0 ) ) / 2.0

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST02'
  write ( *, '(a)' ) '  Use POWER_METHOD2 on the Fibonacci2 matrix.'
  write ( *, '(a)' ) ' '
  write ( *, '(a,i8)'      ) '  Matrix order N         = ', n
  write ( *, '(a,i8)'      ) '  Maximum iterations     = ', it_max
  write ( *, '(a,g14.6)'   ) '  Error tolerance        = ', tol

  call cpu_time ( ctime1 )

  call power_method2 ( n, a, x, it_max, tol, lambda, v, it_num )
  x = real ( v, kind = rk )
  x = x / sqrt ( sum ( x(1:n)**2 ) )

  call cpu_time ( ctime2 )
  ctime = ctime2 - ctime1

  write ( *, '(a)' ) ' '
  write ( *, '(a,i8)'      ) '  Number of iterations   = ', it_num
  write ( *, '(a,g14.6)'   ) '  CPU time               = ', ctime
  write ( *, '(a,2f14.10)' ) '  Estimated eigenvalue   = ', lambda
  write ( *, '(a,f14.10)'  ) '  Correct value          = ', phi
  write ( *, '(a,g14.6)'   ) '  ||Error||              = ', abs ( lambda - phi )
!
!  X2 is the exact eigenvector.
!  Computing it "backwards" this way avoids overflow when N is large.
!
  x2(n) = 1.0D+00
  do i = n, 2, -1
    x2(i-1) = x2(i) / phi
  end do

  x2(1:n) = x2(1:n) / sqrt ( sum ( x2(1:n)**2 ) )
!
!  The sine of the angle between X and X2 is a measure of error.
!
  cos_x1x2 = dot_product ( x(1:n), x2(1:n) )
  sin_x1x2 = sqrt ( ( 1.0 - cos_x1x2 ) * ( 1.0 + cos_x1x2 ) )

  write ( *, '(a)' ) ' '
  write ( *, '(a,f14.10)' ) &
    '  Sine of angle between true and estimated vectors = ', sin_x1x2

  return
end
subroutine test03 ( )

!*****************************************************************************80
!
!! TEST03 uses POWER_METHOD2 on the TRIS matrix.
!
!  Licensing:
!
!    This code is distributed under the MIT license.
!
!  Modified:
!
!    20 July 2008
!
!  Author:
!
!    John Burkardt
!
  implicit none

  integer, parameter :: rk = kind ( 1.0D+00 )

  integer, parameter :: ck = kind ( ( 1.0D+00, 1.0D+00 ) )

  integer, parameter :: n = 100

  real ( kind = rk ) a(n,n)
  real ( kind = rk ) alpha
  real ( kind = rk ) beta
  real ( kind = rk ) ctime
  real ( kind = rk ) ctime1
  real ( kind = rk ) ctime2
  real ( kind = rk ) gamma
  integer i
  integer it_max
  integer it_num
  complex ( kind = ck ) lambda
  complex ( kind = ck ) lambda_max
  complex ( kind = ck ) lambda_vec(n)
  real ( kind = rk ) tol
  complex ( kind = ck ) v(n)
  real ( kind = rk ) x(n)
!
!  If ALPHA * GAMMA is negative, this matrix will have complex eigenvalues.
!
  alpha = -1.0D+00
  beta = 10.0D+00
  gamma = +8.0D+0

  call tris ( n, n, alpha, beta, gamma, a )

  call random_number ( harvest = x(1:n) )

  it_max = 4000
  tol = 0.00000001D+0

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST03'
  write ( *, '(a)' ) '  Use POWER_METHOD2 on the TRIS (tridiagonal scalar) matrix.'
  write ( *, '(a)' ) ' '
  write ( *, '(a,i8)'      ) '  Matrix order N         = ', n
  write ( *, '(a,i8)'      ) '  Maximum iterations     = ', it_max
  write ( *, '(a,g14.6)'   ) '  Error tolerance        = ', tol
!
!  Estimate the solution.
!
  call cpu_time ( ctime1 )

  call power_method2 ( n, a, x, it_max, tol, lambda, v, it_num )

  x = real ( v, kind = rk )
  x = x / sqrt ( sum ( x(1:n)**2 ) )

  call cpu_time ( ctime2 )
  ctime = ctime2 - ctime1

  write ( *, '(a)' ) ' '
  write ( *, '(a,i8)'      ) '  Number of iterations   = ', it_num
  write ( *, '(a,g14.6)'   ) '  CPU time               = ', ctime
  write ( *, '(a,2f14.10)' ) '  Estimated eigenvalue   = ', lambda
!
!  Get the exact solution.
!
  call tris_eigenvalues ( n, alpha, beta, gamma, lambda_vec )

  lambda_max = lambda_vec(1)
  do i = 2, n
    if ( abs ( lambda_max ) < abs ( lambda_vec(i) ) ) then
      lambda_max = lambda_vec(i)
    end if
  end do

  write ( *, '(a,2f14.10)' ) '  Correct max eigenvalue = ', lambda_max
  write ( *, '(a,g14.6)'   ) '  ||Error||              = ', abs ( lambda - lambda_max )

  return
end