program main !*****************************************************************************80 ! !! jacobi_eigenvalue_test() tests jacobi_eigenvalue(). ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 14 July 2013 ! ! Author: ! ! John Burkardt ! implicit none call timestamp ( ) write ( *, '(a)' ) '' write ( *, '(a)' ) 'jacobi_eigenvalue_test():' write ( *, '(a)' ) ' Fortran90 version' write ( *, '(a)' ) ' Test jacobi_eigenvalue().' call test01 ( ) call test02 ( ) call test03 ( ) ! ! Terminate. ! write ( *, '(a)' ) '' write ( *, '(a)' ) 'jacobi_eigenvalue_test():' write ( *, '(a)' ) ' Normal end of execution.' write ( *, '(a)' ) '' call timestamp ( ) stop 0 end subroutine test01 ( ) !*****************************************************************************80 ! !! TEST01 uses a 4x4 test matrix. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 14 July 2013 ! ! Author: ! ! John Burkardt ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) integer, parameter :: n = 4 real ( kind = rk ), dimension ( n, n ) :: a = reshape ( (/ & 4.0, -30.0, 60.0, -35.0, & -30.0, 300.0, -675.0, 420.0, & 60.0, -675.0, 1620.0, -1050.0, & -35.0, 420.0, -1050.0, 700.0 /), (/ 4, 4 /) ) real ( kind = rk ) d(n) real ( kind = rk ) error_frobenius integer it_max integer it_num integer rot_num real ( kind = rk ) v(n,n) write ( *, '(a)' ) '' write ( *, '(a)' ) 'TEST01' write ( *, '(a)' ) ' For a symmetric matrix A,' write ( *, '(a)' ) ' JACOBI_EIGENVALUE computes the eigenvalues D' write ( *, '(a)' ) ' and eigenvectors V so that A * V = D * V.' call r8mat_print ( n, n, a, ' Input matrix A:' ) it_max = 100 call jacobi_eigenvalue ( n, a, it_max, v, d, it_num, rot_num ) write ( *, '(a)' ) '' write ( *, '(a,i4)' ) ' Number of iterations = ', it_num write ( *, '(a,i4)' ) ' Number of rotations = ', rot_num call r8vec_print ( n, d, ' Eigenvalues D:' ) call r8mat_print ( n, n, v, ' Eigenvector matrix V:' ) ! ! Compute eigentest. ! call r8mat_is_eigen_right ( n, n, a, v, d, error_frobenius ) write ( *, '(a)' ) '' write ( *, '(a,g14.6)' ) & ' Frobenius norm error in eigensystem A*V-D*V = ', & error_frobenius return end subroutine test02 ( ) !*****************************************************************************80 ! !! TEST02 uses a 4x4 test matrix. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 14 July 2013 ! ! Author: ! ! John Burkardt ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) integer, parameter :: n = 4 real ( kind = rk ), dimension ( n, n ) :: a = reshape ( (/ & 4.0, 0.0, 0.0, 0.0, & 0.0, 1.0, 0.0, 0.0, & 0.0, 0.0, 3.0, 0.0, & 0.0, 0.0, 0.0, 2.0 /), (/ 4, 4 /) ) real ( kind = rk ) d(n) real ( kind = rk ) error_frobenius integer it_max integer it_num integer rot_num real ( kind = rk ) v(n,n) write ( *, '(a)' ) '' write ( *, '(a)' ) 'TEST02' write ( *, '(a)' ) ' For a symmetric matrix A,' write ( *, '(a)' ) ' JACOBI_EIGENVALUE computes the eigenvalues D' write ( *, '(a)' ) ' and eigenvectors V so that A * V = D * V.' write ( *, '(a)' ) '' write ( *, '(a)' ) 'As a sanity check, input a diagonal matrix.' call r8mat_print ( n, n, a, ' Input matrix A:' ) it_max = 100 call jacobi_eigenvalue ( n, a, it_max, v, d, it_num, rot_num ) write ( *, '(a)' ) '' write ( *, '(a,i4)' ) ' Number of iterations = ', it_num write ( *, '(a,i4)' ) ' Number of rotations = ', rot_num call r8vec_print ( n, d, ' Eigenvalues D:' ) call r8mat_print ( n, n, v, ' Eigenvector matrix V:' ) ! ! Compute eigentest. ! call r8mat_is_eigen_right ( n, n, a, v, d, error_frobenius ) write ( *, '(a)' ) '' write ( *, '(a,g14.6)' ) & ' Frobenius norm error in eigensystem A*V-D*V = ', & error_frobenius return end subroutine test03 ( ) !*****************************************************************************80 ! !! TEST03 uses a 5x5 test matrix. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 14 July 2013 ! ! Author: ! ! John Burkardt ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) integer, parameter :: n = 5 real ( kind = rk ) a(n,n) real ( kind = rk ) d(n) real ( kind = rk ) error_frobenius integer i integer it_max integer it_num integer j integer rot_num real ( kind = rk ) v(n,n) write ( *, '(a)' ) '' write ( *, '(a)' ) 'TEST03' write ( *, '(a)' ) ' For a symmetric matrix A,' write ( *, '(a)' ) ' JACOBI_EIGENVALUE computes the eigenvalues D' write ( *, '(a)' ) ' and eigenvectors V so that A * V = D * V.' write ( *, '(a)' ) '' write ( *, '(a)' ) ' Use the discretized second derivative matrix.' do j = 1, n do i = 1, n if ( i == j ) then a(i,j) = -2.0D+00 else if ( i == j + 1 .or. i == j - 1 ) then a(i,j) = 1.0D+00 else a(i,j) = 0.0D+00 end if end do end do call r8mat_print ( n, n, a, ' Input matrix A:' ) it_max = 100 call jacobi_eigenvalue ( n, a, it_max, v, d, it_num, rot_num ) write ( *, '(a)' ) '' write ( *, '(a,i4)' ) ' Number of iterations = ', it_num write ( *, '(a,i4)' ) ' Number of rotations = ', rot_num call r8vec_print ( n, d, ' Eigenvalues D:' ) call r8mat_print ( n, n, v, ' Eigenvector matrix V:' ) ! ! Compute eigentest. ! call r8mat_is_eigen_right ( n, n, a, v, d, error_frobenius ) write ( *, '(a)' ) '' write ( *, '(a,g14.6)' ) & ' Frobenius norm error in eigensystem A*V-D*V = ', & error_frobenius return end subroutine timestamp ( ) !*****************************************************************************80 ! !! TIMESTAMP prints the current YMDHMS date as a time stamp. ! ! Example: ! ! 31 May 2001 9:45:54.872 AM ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 18 May 2013 ! ! Author: ! ! John Burkardt ! implicit none character ( len = 8 ) ampm integer d integer h integer m integer mm character ( len = 9 ), parameter, dimension(12) :: month = (/ & 'January ', 'February ', 'March ', 'April ', & 'May ', 'June ', 'July ', 'August ', & 'September', 'October ', 'November ', 'December ' /) integer n integer s integer values(8) integer y call date_and_time ( values = values ) y = values(1) m = values(2) d = values(3) h = values(5) n = values(6) s = values(7) mm = values(8) if ( h < 12 ) then ampm = 'AM' else if ( h == 12 ) then if ( n == 0 .and. s == 0 ) then ampm = 'Noon' else ampm = 'PM' end if else h = h - 12 if ( h < 12 ) then ampm = 'PM' else if ( h == 12 ) then if ( n == 0 .and. s == 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