program main

c*****************************************************************************80
c
cc eigs_test() tests eigs().
c
c  Licensing:
c
c    This code is distributed under the MIT license.
c
c  Modified:
c
c    06 June 2024
c
c  Author:
c
c    John Burkardt
c
      implicit none

      double precision, allocatable :: A(:,:)
      double complex, allocatable :: lambda(:)
      integer n

      call timestamp ( )

      write ( *, '(a)' ) ''
      write ( *, '(a)' ) 'eigs_test():'
      write ( *, '(a)' ) '  Fortran77 version'
      write ( *, '(a)' ) '  Test eigs():'

      write ( *, '(a)' ) ''
      write ( *, '(a)' ) '  The Clement1 matrix is symmetric, and has'
      write ( *, '(a)' ) '  integer eigenvalues.'
      n = 5
      allocate ( A(1:n,1:n) )
      call clement1_matrix ( n, A )
      call r8mat_print ( n, n, A, '  The matrix A:' )
      allocate ( lambda(1:n) )
      call eigs ( n, A, lambda )
      call c8vec_print ( n, lambda, '  The eigenvalues LAMBDA' )
      deallocate ( A )
      deallocate ( lambda )

      write ( *, '(a)' ) ''
      write ( *, '(a)' ) '  The Helmert matrix is nonsymmetric,'
      write ( *, '(a)' ) '  orthogonal, with complex eigenvalues'
      write ( *, '(a)' ) '  of norm 1.'
      n = 5
      allocate ( A(1:n,1:n) )
      call helmert_matrix ( n, A )
      call r8mat_print ( n, n, A, '  The matrix A:' )
      allocate ( lambda(1:n) )
      call eigs ( n, A, lambda )
      call c8vec_print ( n, lambda, '  The eigenvalues LAMBDA' )
      deallocate ( A )
      deallocate ( lambda )
c
c  Terminate.
c
      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) 'eigs_test():'
      write ( *, '(a)' ) '  Normal end of execution.'

      call timestamp ( )

      stop 0
      end
      subroutine clement1_matrix ( n, a )

c**********************************************************************72
c
cc clement1() returns the CLEMENT1 matrix.
c
c  Formula:
c
c    if ( J = I + 1 )
c      A(I,J) = sqrt(I*(N-I))
c    else if ( I = J + 1 )
c      A(I,J) = sqrt(J*(N-J))
c    else
c      A(I,J) = 0
c
c  Example:
c
c    N = 5
c
c       .    sqrt(4)    .       .       .
c    sqrt(4)    .    sqrt(6)    .       .
c       .    sqrt(6)    .    sqrt(6)    .
c       .       .    sqrt(6)    .    sqrt(4)
c       .       .       .    sqrt(4)    .
c
c  Properties:
c
c    A is tridiagonal.
c
c    A is banded, with bandwidth 3.
c
c    Because A is tridiagonal, it has property A (bipartite).
c
c    A is symmetric: A' = A.
c
c    Because A is symmetric, it is normal.
c
c    Because A is normal, it is diagonalizable.
c
c    A is persymmetric: A(I,J) = A(N+1-J,N+1-I).
c
c    The diagonal of A is zero.
c
c    A is singular if N is odd.
c
c    About 64 percent of the entries of the inverse of A are zero.
c
c    The eigenvalues are plus and minus the numbers
c
c      N-1, N-3, N-5, ..., (1 or 0).
c
c    If N is even,
c
c      det ( A ) = (-1)^(N/2) * (N-1) * (N+1)^(N/2)
c
c    and if N is odd,
c
c      det ( A ) = 0
c
c  Licensing:
c
c    This code is distributed under the MIT license.
c
c  Modified:
c
c    13 June 2011
c
c  Author:
c
c    John Burkardt
c
c  Reference:
c
c    Paul Clement,
c    A class of triple-diagonal matrices for test purposes,
c    SIAM Review,
c    Volume 1, 1959, pages 50-52.
c
c  Input:
c
c    integer N, the order of the matrix.
c
c  Output:
c
c    double precision A(N,N), the matrix.
c
      implicit none

      integer n

      double precision a(n,n)
      integer i
      integer j

      do j = 1, n
        do i = 1, n

          if ( j .eq. i + 1 ) then
            a(i,j) = sqrt ( dble ( i * ( n - i ) ) )
          else if ( i .eq. j + 1 ) then
            a(i,j) = sqrt ( dble ( j * ( n - j ) ) )
          else
            a(i,j) = 0.0D+00
          end if

        end do
      end do

      return
      end
      subroutine helmert_matrix ( n, a )

c*********************************************************************72
c
cc helmert_matrix() returns the HELMERT matrix.
c
c  Formula:
c
c    If I = 1 then
c      A(I,J) = 1 / sqrt ( N )
c    else if J < I then
c      A(I,J) = 1 / sqrt ( I * ( I - 1 ) )
c    else if J = I then
c      A(I,J) = - sqrt (I-1) / sqrt ( I )
c    else
c      A(I,J) = 0
c
c  Discussion:
c
c    The matrix given above by Helmert is the classic example of
c    a family of matrices which are now called Helmertian or
c    Helmert matrices.
c
c    A matrix is a (standard) Helmert matrix if it is orthogonal,
c    and the elements which are above the diagonal and below the
c    first row are zero.
c
c    If the elements of the first row are all strictly positive,
c    then the matrix is a strictly Helmertian matrix.
c
c    It is possible to require in addition that all elements below
c    the diagonal be strictly positive, but in the reference, this
c    condition is discarded as being cumbersome and not useful.
c
c    A Helmert matrix can be regarded as a change of basis matrix
c    between a pair of orthonormal coordinate bases.  The first row
c    gives the coordinates of the first new basis vector in the old
c    basis.  Each later row describes combinations of (an increasingly
c    extensive set of) old basis vectors that constitute the new
c    basis vectors.
c
c    Helmert matrices have important applications in statistics.
c
c  Example:
c
c    N = 5
c
c    0.4472    0.4472    0.4472    0.4472    0.4472
c    0.7071   -0.7071         0         0         0
c    0.4082    0.4082   -0.8165         0         0
c    0.2887    0.2887    0.2887   -0.8660         0
c    0.2236    0.2236    0.2236    0.2236   -0.8944
c
c  Properties:
c
c    A is generally not symmetric: A' /= A.
c
c    A is orthogonal: A' * A = A * A' = I.
c
c    Because A is orthogonal, it is normal: A' * A = A * A'.
c
c    A is not symmetric: A' /= A.
c
c    det ( A ) = (-1)^(N+1)
c
c  Licensing:
c
c    This code is distributed under the MIT license.
c
c  Modified:
c
c    18 June 2011
c
c  Author:
c
c    John Burkardt
c
c  Reference:
c
c    HO Lancaster,
c    The Helmert Matrices,
c    American Mathematical Monthly,
c    Volume 72, 1965, pages 4-12.
c
c  Input:
c
c    integer N, the order of the matrix.
c
c  Output:
c
c    double precision A(N,N), the matrix.
c
      implicit none

      integer n

      double precision a(n,n)
      integer i
      integer j
c
c  A begins with the first row, diagonal, and lower triangle set to 1.
c
      do i = 1, n
        do j = 1, n

          if ( i .eq. 1 ) then
            a(i,j) = 1.0D+00 / sqrt ( dble ( n ) )
          else if ( j .lt. i ) then
            a(i,j) = 1.0D+00 / sqrt ( dble ( i * ( i - 1 ) ) )
          else if ( i .eq. j ) then
            a(i,j) = - sqrt ( dble ( i - 1 ) ) 
     &               / sqrt ( dble ( i ) )
          else
            a(i,j) = 0.0D+00
          end if

        end do
      end do

      return
      end
      subroutine c8vec_print ( n, a, title )

c*********************************************************************72
c
cc c8vec_print() prints a C8VEC, with an optional title.
c
c  Discussion:
c
c    A C8VEC is a vector of C8's
c
c  Licensing:
c
c    This code is distributed under the MIT license. 
c
c  Modified:
c
c    07 December 2008
c
c  Author:
c
c    John Burkardt
c
c  Input:
c
c    integer N, the number of components of the vector.
c
c    double complex A(N), the vector to be printed.
c
c    character*(*) TITLE, a title.
c
      implicit none

      integer n

      double complex a(n)
      integer i
      character*(*) title

      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) trim ( title )
      write ( *, '(a)' ) ' '
      do i = 1, n
        write ( *, '(2x,i8,a,1x,2g14.6)' ) i, ':', a(i)
      end do

      return
      end
      subroutine r8mat_print ( m, n, a, title )

c*********************************************************************72
c
cc r8mat_print() prints an R8MAT.
c
c  Discussion:
c
c    An R8MAT is an array of R8's.
c
c  Licensing:
c
c    This code is distributed under the MIT license.
c
c  Modified:
c
c    20 May 2004
c
c  Author:
c
c    John Burkardt
c
c  Input:
c
c    integer M, the number of rows in A.
c
c    integer N, the number of columns in A.
c
c    double precision A(M,N), the matrix.
c
c    character * ( * ) TITLE, a title.
c
      implicit none

      integer m
      integer n

      double precision a(m,n)
      character * ( * ) title

      call r8mat_print_some ( m, n, a, 1, 1, m, n, title )

      return
      end
      subroutine r8mat_print_some ( m, n, a, ilo, jlo, ihi, jhi,
     &  title )

c*********************************************************************72
c
cc r8mat_print_some() prints some of an R8MAT.
c
c  Discussion:
c
c    An R8MAT is an array of R8's.
c
c  Licensing:
c
c    This code is distributed under the MIT license.
c
c  Modified:
c
c    25 January 2007
c
c  Author:
c
c    John Burkardt
c
c  Input:
c
c    integer M, N, the number of rows and columns.
c
c    double precision A(M,N), an M by N matrix to be printed.
c
c    integer ILO, JLO, the first row and column to print.
c
c    integer IHI, JHI, the last row and column to print.
c
c    character * ( * ) TITLE, a title.
c
      implicit none

      integer incx
      parameter ( incx = 5 )
      integer m
      integer n

      double precision a(m,n)
      character * ( 14 ) ctemp(incx)
      integer i
      integer i2hi
      integer i2lo
      integer ihi
      integer ilo
      integer inc
      integer j
      integer j2
      integer j2hi
      integer j2lo
      integer jhi
      integer jlo
      character * ( * ) title

      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) trim ( title )

      if ( m .le. 0 .or. n .le. 0 ) then
        write ( *, '(a)' ) ' '
        write ( *, '(a)' ) '  (None)'
        return
      end if

      do j2lo = max ( jlo, 1 ), min ( jhi, n ), incx

        j2hi = j2lo + incx - 1
        j2hi = min ( j2hi, n )
        j2hi = min ( j2hi, jhi )

        inc = j2hi + 1 - j2lo

        write ( *, '(a)' ) ' '

        do j = j2lo, j2hi
          j2 = j + 1 - j2lo
          write ( ctemp(j2), '(i7,7x)') j
        end do

        write ( *, '(''  Col   '',5a14)' ) ( ctemp(j), j = 1, inc )
        write ( *, '(a)' ) '  Row'
        write ( *, '(a)' ) ' '

        i2lo = max ( ilo, 1 )
        i2hi = min ( ihi, m )

        do i = i2lo, i2hi

          do j2 = 1, inc

            j = j2lo - 1 + j2

            write ( ctemp(j2), '(g14.6)' ) a(i,j)

          end do

          write ( *, '(i5,a,5a14)' ) i, ':', ( ctemp(j), j = 1, inc )

        end do

      end do

      return
      end
      subroutine timestamp ( )

c*********************************************************************72
c
cc timestamp() prints the YMDHMS date as a timestamp.
c
c  Licensing:
c
c    This code is distributed under the MIT license.
c
c  Modified:
c
c    12 June 2014
c
c  Author:
c
c    John Burkardt
c
      implicit none

      character * ( 8 ) ampm
      integer d
      character * ( 8 ) date
      integer h
      integer m
      integer mm
      character * ( 9 ) month(12)
      integer n
      integer s
      character * ( 10 ) time
      integer y

      save month

      data month /
     &  'January  ', 'February ', 'March    ', 'April    ', 
     &  'May      ', 'June     ', 'July     ', 'August   ', 
     &  'September', 'October  ', 'November ', 'December ' /

      call date_and_time ( date, time )

      read ( date, '(i4,i2,i2)' ) y, m, d
      read ( time, '(i2,i2,i2,1x,i3)' ) h, n, s, mm

      if ( h .lt. 12 ) then
        ampm = 'AM'
      else if ( h .eq. 12 ) then
        if ( n .eq. 0 .and. s .eq. 0 ) then
          ampm = 'Noon'
        else
          ampm = 'PM'
        end if
      else
        h = h - 12
        if ( h .lt. 12 ) then
          ampm = 'PM'
        else if ( h .eq. 12 ) then
          if ( n .eq. 0 .and. s .eq. 0 ) then
            ampm = 'Midnight'
          else
            ampm = 'AM'
          end if
        end if
      end if

      write ( *, 
     &  '(i2,1x,a,1x,i4,2x,i2,a1,i2.2,a1,i2.2,a1,i3.3,1x,a)' ) 
     &  d, trim ( month(m) ), y, h, ':', n, ':', s, '.', mm, 
     &  trim ( ampm )

      return
      end