program main !*****************************************************************************80 ! !! condition_test() tests condition(). ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 31 March 2012 ! ! Author: ! ! John Burkardt ! implicit none call timestamp ( ) write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'condition_test():' write ( *, '(a)' ) ' Fortran90 version' write ( *, '(a)' ) ' Test condition().' write ( *, '(a)' ) ' The R8LIB library must also be available.' call condition_linpack_test ( ) call condition_sample1_test ( ) call condition_hager_test ( ) ! ! Terminate. ! write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'CONDITION_TEST:' write ( *, '(a)' ) ' Normal end of execution.' write ( *, '(a)' ) ' ' call timestamp ( ) stop 0 end subroutine condition_linpack_test ( ) !*****************************************************************************80 ! !! condition_linpack_test() tests condition_linpack(). ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 31 March 2012 ! ! Author: ! ! John Burkardt ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) real ( kind = rk ), allocatable :: a(:,:) real ( kind = rk ), allocatable :: a_inverse(:,:) real ( kind = rk ) a_inverse_norm_l1 real ( kind = rk ) a_norm_l1 real ( kind = rk ) alpha real ( kind = rk ) beta real ( kind = rk ) cond real ( kind = rk ) cond_l1 integer i integer info integer n character ( len = 80 ) name integer, allocatable :: pivot(:) real ( kind = rk ) r8mat_norm_l1 integer seed real ( kind = rk ), allocatable :: z(:) write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'CONDITION_LINPACK_TEST' write ( *, '(a)' ) ' CONDITION_LINPACK estimates the L1 condition number' write ( *, '(a)' ) ' of a matrix in general storage mode.' write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' Matrix Order Condition Linpack' write ( *, '(a)' ) ' ' ! ! Combinatorial matrix. ! name = 'Combinatorial' n = 4 allocate ( a(1:n,1:n) ) allocate ( a_inverse(1:n,1:n) ) allocate ( pivot(1:n) ) allocate ( z(1:n) ) alpha = 2.0D+00 beta = 3.0D+00 call combin ( alpha, beta, n, a ) call combin_inverse ( alpha, beta, n, a_inverse ) a_norm_l1 = r8mat_norm_l1 ( n, n, a ) a_inverse_norm_l1 = r8mat_norm_l1 ( n, n, a_inverse ) cond_l1 = a_norm_l1 * a_inverse_norm_l1 call condition_linpack ( n, a, pivot, cond, z ) write ( *, '(2x,a20,2x,i4,2x,g14.6,2x,g14.6)' ) name, n, cond_l1, cond deallocate ( a ) deallocate ( a_inverse ) deallocate ( pivot ) deallocate ( z ) ! ! CONEX1 ! name = 'CONEX1' n = 4 allocate ( a(1:n,1:n) ) allocate ( a_inverse(1:n,1:n) ) allocate ( pivot(1:n) ) allocate ( z(1:n) ) alpha = 100.0D+00 call conex1 ( alpha, a ) call conex1_inverse ( alpha, a_inverse ) a_norm_l1 = r8mat_norm_l1 ( n, n, a ) a_inverse_norm_l1 = r8mat_norm_l1 ( n, n, a_inverse ) cond_l1 = a_norm_l1 * a_inverse_norm_l1 call condition_linpack ( n, a, pivot, cond, z ) write ( *, '(2x,a20,2x,i4,2x,g14.6,2x,g14.6)' ) name, n, cond_l1, cond deallocate ( a ) deallocate ( a_inverse ) deallocate ( pivot ) deallocate ( z ) ! ! CONEX2 ! name = 'CONEX2' n = 3 allocate ( a(1:n,1:n) ) allocate ( a_inverse(1:n,1:n) ) allocate ( pivot(1:n) ) allocate ( z(1:n) ) alpha = 100.0D+00 call conex2 ( alpha, a ) call conex2_inverse ( alpha, a_inverse ) a_norm_l1 = r8mat_norm_l1 ( n, n, a ) a_inverse_norm_l1 = r8mat_norm_l1 ( n, n, a_inverse ) cond_l1 = a_norm_l1 * a_inverse_norm_l1 call condition_linpack ( n, a, pivot, cond, z ) write ( *, '(2x,a20,2x,i4,2x,g14.6,2x,g14.6)' ) name, n, cond_l1, cond deallocate ( a ) deallocate ( a_inverse ) deallocate ( pivot ) deallocate ( z ) ! ! CONEX3 ! name = 'CONEX3' n = 5 allocate ( a(1:n,1:n) ) allocate ( a_inverse(1:n,1:n) ) allocate ( pivot(1:n) ) allocate ( z(1:n) ) call conex3 ( n, a ) call conex3_inverse ( n, a_inverse ) a_norm_l1 = r8mat_norm_l1 ( n, n, a ) a_inverse_norm_l1 = r8mat_norm_l1 ( n, n, a_inverse ) cond_l1 = a_norm_l1 * a_inverse_norm_l1 call condition_linpack ( n, a, pivot, cond, z ) write ( *, '(2x,a20,2x,i4,2x,g14.6,2x,g14.6)' ) name, n, cond_l1, cond deallocate ( a ) deallocate ( a_inverse ) deallocate ( pivot ) deallocate ( z ) ! ! CONEX4 ! name = 'CONEX4' n = 4 allocate ( a(1:n,1:n) ) allocate ( a_inverse(1:n,1:n) ) allocate ( pivot(1:n) ) allocate ( z(1:n) ) call conex4 ( a ) call conex4_inverse ( a_inverse ) a_norm_l1 = r8mat_norm_l1 ( n, n, a ) a_inverse_norm_l1 = r8mat_norm_l1 ( n, n, a_inverse ) cond_l1 = a_norm_l1 * a_inverse_norm_l1 call condition_linpack ( n, a, pivot, cond, z ) write ( *, '(2x,a20,2x,i4,2x,g14.6,2x,g14.6)' ) name, n, cond_l1, cond deallocate ( a ) deallocate ( a_inverse ) deallocate ( pivot ) deallocate ( z ) ! ! KAHAN ! name = 'KAHAN' n = 4 allocate ( a(1:n,1:n) ) allocate ( a_inverse(1:n,1:n) ) allocate ( pivot(1:n) ) allocate ( z(1:n) ) alpha = 0.25D+00 call kahan ( alpha, n, n, a ) call kahan_inverse ( alpha, n, a_inverse ) a_norm_l1 = r8mat_norm_l1 ( n, n, a ) a_inverse_norm_l1 = r8mat_norm_l1 ( n, n, a_inverse ) cond_l1 = a_norm_l1 * a_inverse_norm_l1 call condition_linpack ( n, a, pivot, cond, z ) write ( *, '(2x,a20,2x,i4,2x,g14.6,2x,g14.6)' ) name, n, cond_l1, cond deallocate ( a ) deallocate ( a_inverse ) deallocate ( pivot ) deallocate ( z ) ! ! Random ! seed = 123456789 do i = 1, 5 name = 'RANDOM' n = 4 allocate ( a(1:n,1:n) ) allocate ( a_inverse(1:n,1:n) ) allocate ( pivot(1:n) ) allocate ( z(1:n) ) call random_number ( harvest = a(1:n,1:n) ) a_inverse(1:n,1:n) = a(1:n,1:n) call r8ge_fa ( n, a_inverse, pivot, info ) call r8ge_inverse ( n, a_inverse, pivot ) a_norm_l1 = r8mat_norm_l1 ( n, n, a ) a_inverse_norm_l1 = r8mat_norm_l1 ( n, n, a_inverse ) cond_l1 = a_norm_l1 * a_inverse_norm_l1 call condition_linpack ( n, a, pivot, cond, z ) write ( *, '(2x,a20,2x,i4,2x,g14.6,2x,g14.6)' ) name, n, cond_l1, cond deallocate ( a ) deallocate ( a_inverse ) deallocate ( pivot ) deallocate ( z ) end do return end subroutine condition_sample1_test ( ) !*****************************************************************************80 ! !! condition_sample1_test() tests condition_sample1(). ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 31 March 2012 ! ! Author: ! ! John Burkardt ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) real ( kind = rk ), allocatable :: a(:,:) real ( kind = rk ), allocatable :: a_inverse(:,:) real ( kind = rk ) a_inverse_norm_l1 real ( kind = rk ) a_norm_l1 real ( kind = rk ) alpha real ( kind = rk ) beta real ( kind = rk ) cond real ( kind = rk ) cond_l1 integer i integer info integer j integer m integer :: m_test(3) = (/ 10, 1000, 100000 /) integer n character ( len = 80 ) name integer, allocatable :: pivot(:) real ( kind = rk ) r8mat_norm_l1 integer seed write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'condition_sample1_test' write ( *, '(a)' ) ' condition_sample1 estimates the L1 condition number using sampling' write ( *, '(a)' ) ' for a matrix in general storage,' write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' Matrix Samples Order Condition Estimate' ! ! Combinatorial matrix. ! name = 'Combinatorial' n = 4 allocate ( a(1:n,1:n) ) allocate ( a_inverse(1:n,1:n) ) alpha = 2.0D+00 beta = 3.0D+00 call combin ( alpha, beta, n, a ) call combin_inverse ( alpha, beta, n, a_inverse ) a_norm_l1 = r8mat_norm_l1 ( n, n, a ) a_inverse_norm_l1 = r8mat_norm_l1 ( n, n, a_inverse ) cond_l1 = a_norm_l1 * a_inverse_norm_l1 write ( *, '(a)' ) ' ' do i = 1, 3 m = m_test(i) call condition_sample1 ( n, a, m, cond ) write ( *, '(2x,a20,2x,i8,2x,i4,2x,g14.6,2x,g14.6)' ) name, m, n, cond_l1, cond end do deallocate ( a ) deallocate ( a_inverse ) ! ! CONEX1 ! name = 'CONEX1' n = 4 allocate ( a(1:n,1:n) ) allocate ( a_inverse(1:n,1:n) ) alpha = 100.0D+00 call conex1 ( alpha, a ) call conex1_inverse ( alpha, a_inverse ) a_norm_l1 = r8mat_norm_l1 ( n, n, a ) a_inverse_norm_l1 = r8mat_norm_l1 ( n, n, a_inverse ) cond_l1 = a_norm_l1 * a_inverse_norm_l1 write ( *, '(a)' ) ' ' do i = 1, 3 m = m_test(i) call condition_sample1 ( n, a, m, cond ) write ( *, '(2x,a20,2x,i8,2x,i4,2x,g14.6,2x,g14.6)' ) name, m, n, cond_l1, cond end do deallocate ( a ) deallocate ( a_inverse ) ! ! CONEX2 ! name = 'CONEX2' n = 3 allocate ( a(1:n,1:n) ) allocate ( a_inverse(1:n,1:n) ) alpha = 100.0D+00 call conex2 ( alpha, a ) call conex2_inverse ( alpha, a_inverse ) a_norm_l1 = r8mat_norm_l1 ( n, n, a ) a_inverse_norm_l1 = r8mat_norm_l1 ( n, n, a_inverse ) cond_l1 = a_norm_l1 * a_inverse_norm_l1 write ( *, '(a)' ) ' ' do i = 1, 3 m = m_test(i) call condition_sample1 ( n, a, m, cond ) write ( *, '(2x,a20,2x,i8,2x,i4,2x,g14.6,2x,g14.6)' ) name, m, n, cond_l1, cond end do deallocate ( a ) deallocate ( a_inverse ) ! ! CONEX3 ! name = 'CONEX3' n = 5 allocate ( a(1:n,1:n) ) allocate ( a_inverse(1:n,1:n) ) call conex3 ( n, a ) call conex3_inverse ( n, a_inverse ) a_norm_l1 = r8mat_norm_l1 ( n, n, a ) a_inverse_norm_l1 = r8mat_norm_l1 ( n, n, a_inverse ) cond_l1 = a_norm_l1 * a_inverse_norm_l1 write ( *, '(a)' ) ' ' do i = 1, 3 m = m_test(i) call condition_sample1 ( n, a, m, cond ) write ( *, '(2x,a20,2x,i8,2x,i4,2x,g14.6,2x,g14.6)' ) name, m, n, cond_l1, cond end do deallocate ( a ) deallocate ( a_inverse ) ! ! CONEX4 ! name = 'CONEX4' n = 4 allocate ( a(1:n,1:n) ) allocate ( a_inverse(1:n,1:n) ) call conex4 ( a ) call conex4_inverse ( a_inverse ) a_norm_l1 = r8mat_norm_l1 ( n, n, a ) a_inverse_norm_l1 = r8mat_norm_l1 ( n, n, a_inverse ) cond_l1 = a_norm_l1 * a_inverse_norm_l1 write ( *, '(a)' ) ' ' do i = 1, 3 m = m_test(i) call condition_sample1 ( n, a, m, cond ) write ( *, '(2x,a20,2x,i8,2x,i4,2x,g14.6,2x,g14.6)' ) name, m, n, cond_l1, cond end do deallocate ( a ) deallocate ( a_inverse ) ! ! KAHAN ! name = 'KAHAN' n = 4 allocate ( a(1:n,1:n) ) allocate ( a_inverse(1:n,1:n) ) alpha = 0.25D+00 call kahan ( alpha, n, n, a ) call kahan_inverse ( alpha, n, a_inverse ) a_norm_l1 = r8mat_norm_l1 ( n, n, a ) a_inverse_norm_l1 = r8mat_norm_l1 ( n, n, a_inverse ) cond_l1 = a_norm_l1 * a_inverse_norm_l1 write ( *, '(a)' ) ' ' do i = 1, 3 m = m_test(i) call condition_sample1 ( n, a, m, cond ) write ( *, '(2x,a20,2x,i8,2x,i4,2x,g14.6,2x,g14.6)' ) name, m, n, cond_l1, cond end do deallocate ( a ) deallocate ( a_inverse ) ! ! Random ! seed = 123456789 do j = 1, 5 name = 'RANDOM' n = 4 allocate ( a(1:n,1:n) ) allocate ( a_inverse(1:n,1:n) ) allocate ( pivot(1:n) ) call random_number ( harvest = a(1:n,1:n) ) a_inverse(1:n,1:n) = a(1:n,1:n) call r8ge_fa ( n, a_inverse, pivot, info ) call r8ge_inverse ( n, a_inverse, pivot ) a_norm_l1 = r8mat_norm_l1 ( n, n, a ) a_inverse_norm_l1 = r8mat_norm_l1 ( n, n, a_inverse ) cond_l1 = a_norm_l1 * a_inverse_norm_l1 write ( *, '(a)' ) ' ' do i = 1, 3 m = m_test(i) call condition_sample1 ( n, a, m, cond ) write ( *, '(2x,a20,2x,i8,2x,i4,2x,g14.6,2x,g14.6)' ) name, m, n, cond_l1, cond end do deallocate ( a ) deallocate ( a_inverse ) deallocate ( pivot ) end do return end subroutine condition_hager_test ( ) !*****************************************************************************80 ! !! condition_hager_test() tests condition_hager(). ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 08 April 2012 ! ! Author: ! ! John Burkardt ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) real ( kind = rk ), allocatable :: a(:,:) real ( kind = rk ), allocatable :: a_inverse(:,:) real ( kind = rk ) a_inverse_norm_l1 real ( kind = rk ) a_norm_l1 real ( kind = rk ) alpha real ( kind = rk ) beta real ( kind = rk ) cond real ( kind = rk ) cond_l1 integer i integer info integer n character ( len = 80 ) name integer, allocatable :: pivot(:) real ( kind = rk ) r8mat_norm_l1 integer seed write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'CONDITION_HAGER_TEST' write ( *, '(a)' ) ' CONDITION_HAGER estimates the L1 condition number' write ( *, '(a)' ) ' for a matrix in general storage,' write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' Matrix Order Condition Hager' write ( *, '(a)' ) ' ' ! ! Combinatorial matrix. ! name = 'Combinatorial' n = 4 allocate ( a(1:n,1:n) ) allocate ( a_inverse(1:n,1:n) ) alpha = 2.0D+00 beta = 3.0D+00 call combin ( alpha, beta, n, a ) call combin_inverse ( alpha, beta, n, a_inverse ) a_norm_l1 = r8mat_norm_l1 ( n, n, a ) a_inverse_norm_l1 = r8mat_norm_l1 ( n, n, a_inverse ) cond_l1 = a_norm_l1 * a_inverse_norm_l1 call condition_hager ( n, a, cond ) write ( *, '(2x,a20,2x,i4,2x,g14.6,2x,g14.6)' ) name, n, cond_l1, cond deallocate ( a ) deallocate ( a_inverse ) ! ! CONEX1 ! name = 'CONEX1' n = 4 allocate ( a(1:n,1:n) ) allocate ( a_inverse(1:n,1:n) ) alpha = 100.0D+00 call conex1 ( alpha, a ) call conex1_inverse ( alpha, a_inverse ) a_norm_l1 = r8mat_norm_l1 ( n, n, a ) a_inverse_norm_l1 = r8mat_norm_l1 ( n, n, a_inverse ) cond_l1 = a_norm_l1 * a_inverse_norm_l1 call condition_hager ( n, a, cond ) write ( *, '(2x,a20,2x,i4,2x,g14.6,2x,g14.6)' ) name, n, cond_l1, cond deallocate ( a ) deallocate ( a_inverse ) ! ! CONEX2 ! name = 'CONEX2' n = 3 allocate ( a(1:n,1:n) ) allocate ( a_inverse(1:n,1:n) ) alpha = 100.0D+00 call conex2 ( alpha, a ) call conex2_inverse ( alpha, a_inverse ) a_norm_l1 = r8mat_norm_l1 ( n, n, a ) a_inverse_norm_l1 = r8mat_norm_l1 ( n, n, a_inverse ) cond_l1 = a_norm_l1 * a_inverse_norm_l1 call condition_hager ( n, a, cond ) write ( *, '(2x,a20,2x,i4,2x,g14.6,2x,g14.6)' ) name, n, cond_l1, cond deallocate ( a ) deallocate ( a_inverse ) ! ! CONEX3 ! name = 'CONEX3' n = 5 allocate ( a(1:n,1:n) ) allocate ( a_inverse(1:n,1:n) ) call conex3 ( n, a ) call conex3_inverse ( n, a_inverse ) a_norm_l1 = r8mat_norm_l1 ( n, n, a ) a_inverse_norm_l1 = r8mat_norm_l1 ( n, n, a_inverse ) cond_l1 = a_norm_l1 * a_inverse_norm_l1 call condition_hager ( n, a, cond ) write ( *, '(2x,a20,2x,i4,2x,g14.6,2x,g14.6)' ) name, n, cond_l1, cond deallocate ( a ) deallocate ( a_inverse ) ! ! CONEX4 ! name = 'CONEX4' n = 4 allocate ( a(1:n,1:n) ) allocate ( a_inverse(1:n,1:n) ) call conex4 ( a ) call conex4_inverse ( a_inverse ) a_norm_l1 = r8mat_norm_l1 ( n, n, a ) a_inverse_norm_l1 = r8mat_norm_l1 ( n, n, a_inverse ) cond_l1 = a_norm_l1 * a_inverse_norm_l1 call condition_hager ( n, a, cond ) write ( *, '(2x,a20,2x,i4,2x,g14.6,2x,g14.6)' ) name, n, cond_l1, cond deallocate ( a ) deallocate ( a_inverse ) ! ! KAHAN ! name = 'KAHAN' n = 4 allocate ( a(1:n,1:n) ) allocate ( a_inverse(1:n,1:n) ) alpha = 0.25D+00 call kahan ( alpha, n, n, a ) call kahan_inverse ( alpha, n, a_inverse ) a_norm_l1 = r8mat_norm_l1 ( n, n, a ) a_inverse_norm_l1 = r8mat_norm_l1 ( n, n, a_inverse ) cond_l1 = a_norm_l1 * a_inverse_norm_l1 call condition_hager ( n, a, cond ) write ( *, '(2x,a20,2x,i4,2x,g14.6,2x,g14.6)' ) name, n, cond_l1, cond deallocate ( a ) deallocate ( a_inverse ) ! ! Random ! seed = 123456789 do i = 1, 5 name = 'RANDOM' n = 4 allocate ( a(1:n,1:n) ) allocate ( a_inverse(1:n,1:n) ) allocate ( pivot(1:n) ) call random_number ( harvest = a(1:n,1:n) ) a_inverse(1:n,1:n) = a(1:n,1:n) call r8ge_fa ( n, a_inverse, pivot, info ) call r8ge_inverse ( n, a_inverse, pivot ) a_norm_l1 = r8mat_norm_l1 ( n, n, a ) a_inverse_norm_l1 = r8mat_norm_l1 ( n, n, a_inverse ) cond_l1 = a_norm_l1 * a_inverse_norm_l1 call condition_hager ( n, a, cond ) write ( *, '(2x,a20,2x,i4,2x,g14.6,2x,g14.6)' ) name, n, cond_l1, cond deallocate ( a ) deallocate ( a_inverse ) deallocate ( pivot ) 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.2,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