subroutine dif2 ( m, n, a ) !*****************************************************************************80 ! !! dif2() returns the DIF2 matrix. ! ! Example: ! ! N = 5 ! ! 2 -1 . . . ! -1 2 -1 . . ! . -1 2 -1 . ! . . -1 2 -1 ! . . . -1 2 ! ! Properties: ! ! A is banded, with bandwidth 3. ! ! A is tridiagonal. ! ! Because A is tridiagonal, it has property A (bipartite). ! ! A is a special case of the TRIS or tridiagonal scalar matrix. ! ! A is integral, therefore det ( A ) is integral, and ! det ( A ) * inverse ( A ) is integral. ! ! A is Toeplitz: constant along diagonals. ! ! A is symmetric: A' = A. ! ! Because A is symmetric, it is normal. ! ! Because A is normal, it is diagonalizable. ! ! A is persymmetric: A(I,J) = A(N+1-J,N+1-I). ! ! A is positive definite. ! ! A is an M matrix. ! ! A is weakly diagonally dominant, but not strictly diagonally dominant. ! ! A has an LU factorization A = L * U, without pivoting. ! ! The matrix L is lower bidiagonal with subdiagonal elements: ! ! L(I+1,I) = -I/(I+1) ! ! The matrix U is upper bidiagonal, with diagonal elements ! ! U(I,I) = (I+1)/I ! ! and superdiagonal elements which are all -1. ! ! A has a Cholesky factorization A = L * L', with L lower bidiagonal. ! ! L(I,I) = sqrt ( (I+1) / I ) ! L(I,I-1) = -sqrt ( (I-1) / I ) ! ! The eigenvalues are ! ! LAMBDA(I) = 2 + 2 * COS(I*PI/(N+1)) ! = 4 SIN^2(I*PI/(2*N+2)) ! ! The corresponding eigenvector X(I) has entries ! ! X(I)(J) = sqrt(2/(N+1)) * sin ( I*J*PI/(N+1) ). ! ! Simple linear systems: ! ! x = (1,1,1,...,1,1), A*x=(1,0,0,...,0,1) ! ! x = (1,2,3,...,n-1,n), A*x=(0,0,0,...,0,n+1) ! ! det ( A ) = N + 1. ! ! The value of the determinant can be seen by induction, ! and expanding the determinant across the first row: ! ! det ( A(N) ) = 2 * det ( A(N-1) ) - (-1) * (-1) * det ( A(N-2) ) ! = 2 * N - (N-1) ! = N + 1 ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 15 August 2021 ! ! Author: ! ! John Burkardt ! ! Reference: ! ! Robert Gregory, David Karney, ! A Collection of Matrices for Testing Computational Algorithms, ! Wiley, 1969, ! ISBN: 0882756494, ! LC: QA263.68 ! ! Morris Newman, John Todd, ! Example A8, ! The evaluation of matrix inversion programs, ! Journal of the Society for Industrial and Applied Mathematics, ! Volume 6, Number 4, pages 466-476, 1958. ! ! John Todd, ! Basic Numerical Mathematics, ! Volume 2: Numerical Algebra, ! Birkhauser, 1980, ! ISBN: 0817608117, ! LC: QA297.T58. ! ! Joan Westlake, ! A Handbook of Numerical Matrix Inversion and Solution of ! Linear Equations, ! John Wiley, 1968, ! ISBN13: 978-0471936756, ! LC: QA263.W47. ! ! Input: ! ! integer M, N, the order of the matrix. ! ! Output: ! ! real A(M,N), the matrix. ! implicit none integer, parameter :: rk = kind ( 1.0d+00 ) integer m integer n real ( kind = rk ) a(m,n) integer i integer j do j = 1, n do i = 1, m if ( j == i - 1 ) then a(i,j) = -1.0D+00 else if ( j == i ) then a(i,j) = 2.0D+00 else if ( j == i + 1 ) then a(i,j) = -1.0D+00 else a(i,j) = 0.0D+00 end if end do end do return end subroutine jacobi1 ( n, a, b, x, x_new ) !*****************************************************************************80 ! !! jacobi1() carries out one step of the Jacobi iteration. ! ! Discussion: ! ! The linear system A*x=b is to be solved. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 15 August 2021 ! ! Author: ! ! John Burkardt ! ! Input: ! ! integer N, the order of the matrix. ! ! real A(N,N), the matrix. ! ! real B(N), the right hand side. ! ! real X(N), the current solution estimate. ! ! Output: ! ! real X_NEW(N), the solution estimate updated by ! one step of the Jacobi iteration. ! implicit none integer, parameter :: rk = kind ( 1.0d+00 ) integer n real ( kind = rk ) a(n,n) real ( kind = rk ) b(n) integer i integer j real ( kind = rk ) x(n) real ( kind = rk ) x_new(n) do i = 1, n x_new(i) = b(i) do j = 1, n if ( j /= i ) then x_new(i) = x_new(i) - a(i,j) * x(j) end if end do x_new(i) = x_new(i) / a(i,i) end do return end subroutine r8vec_print ( n, a, title ) !*****************************************************************************80 ! !! r8vec_print() prints an R8VEC. ! ! Discussion: ! ! An R8VEC is a vector of R8's. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 15 August 2021 ! ! Author: ! ! John Burkardt ! ! Input: ! ! integer N, the number of components of the vector. ! ! real A(N), the vector to be printed. ! ! character ( len = * ) TITLE, a title. ! implicit none integer, parameter :: rk = kind ( 1.0d+00 ) integer n real ( kind = rk ) a(n) integer i character ( len = * ) title write ( *, '(a)' ) ' ' write ( *, '(a)' ) trim ( title ) write ( *, '(a)' ) ' ' do i = 1, n write ( *, '(2x,i8,a,1x,g16.8)' ) i, ':', a(i) end do 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: ! ! 15 August 2021 ! ! 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