IEEE 丸めモードで区間演算の代用

数値的に不安定な古典的な例として有名なヒルベルト行列の逆行列計算を、丸めモードを変えて計算してみます。
IEEE7754 の丸めモードの切り捨てと切り上げが、区間演算の下限と上限のようなものとして代用できるというアイデアを試してみます。
近接値への丸めによる計算が、必ずしも切り捨てと切り上げの値に挟まれるわけではないようです。

実行結果

単精度での計算では丸めモードの違いで結果に差が出ます。5x5 の Hilbert 行列で約 0.2% の違いが表れています。

倍精度では違いは出ません。

ソース・プログラム

program Hilbert
   use, intrinsic :: ieee_arithmetic
   implicit none
   integer, parameter :: kd = kind(1.0d0)
   real(kd) :: a(5, 5, 4)
   integer :: i, j

   print *, 'round to nearest'
   call ieee_set_rounding_mode(IEEE_NEAREST)
   call sub_hilbert(a(:, :, 1))
   print *, 'round to zero'
   call ieee_set_rounding_mode(IEEE_TO_ZERO)
   call sub_hilbert(a(:, :, 2))
   print *, 'round to down'
   call ieee_set_rounding_mode(IEEE_DOWN)
   call sub_hilbert(a(:, :, 3))
   print *, 'round to up'
   call ieee_set_rounding_mode(IEEE_UP)
   call sub_hilbert(a(:, :, 4))
   
   print *
   print *, 'absolute difference'
   print '(5f15.6)', ((abs(maxval(a(i, j, :)) - minval(a(i, j, :))), j = 1, 5), i = 1, 5)
   
   print *, 'relative difference'
   print '(5(f14.2, "%"))', ((100 * abs(maxval(a(i, j, :)) - minval(a(i, j, :))) / &
                                    maxval(abs(a(i, j, :))), j = 1, 5), i = 1, 5)
   stop   

contains
   subroutine sub_hilbert(a)
     real(kd), intent(in out) :: a(5, 5)
     integer :: i, j
     forall(i = 1:5, j = 1:5) a(i, j) = 1.0_kd / (i + j - 1)
     call inverse(a)
     print '(5f15.6)', (a(i, :), i = 1, 5)
     return
   end subroutine sub_hilbert

   pure subroutine inverse(a)
     real(kd), intent(in out) :: a(5, 5)
     real(kd) :: d
     integer :: k, i
     do k = 1, 5
       d = a(k, k)
       a(k, k) = 1.0_kd
       a(k, :) = a(k, :) / d
       do i = 1, 5
         if (i == k) cycle
         d = a(i, k)
         a(i, k) = 0.0_kd
         a(i, :) = a(i, :) - d * a(k, :)
       end do
     end do
     return
   end subroutine inverse

end program Hilbert