! time_simeq.f90

subroutine init_matrix(n, sz, A, Y)
  implicit none
  integer, intent(in) :: n
  integer, intent(in) :: sz
  double precision, dimension(sz,sz), intent(out) :: A
  double precision, dimension(sz), intent(out) :: Y
  integer :: i, j

  interface
    double precision function udrnrt()
    end function udrnrt
  end interface

  do i=1,n
    Y(i) = udrnrt();
    do j=1,n
      A(i,j) = udrnrt()
    end do ! j
  end do ! i
  if(n<5) then
    print *, 'initial A and Y'
    do i=1,n
      print *, (A(i,j), j=1,n)
      print *, 'Y(', i, ')=', Y(i)
    end do ! i
  end if
end subroutine init_matrix

program time_simeq
  implicit none
  integer, parameter :: sz = 1024 
  double precision, dimension(sz,sz) :: A
  double precision, dimension(sz) :: X
  double precision, dimension(sz) :: Y
  integer :: i, j, n
  double precision :: t, total
  double precision :: start, now, next
  character*20 :: start_date
  character*20 :: start_time
  integer :: clock_count, clock_rate, clock_max

  interface
    subroutine simeq(n, sz, A, Y, X)
      integer, intent(in) :: n
      integer, intent(in) :: sz
      double precision, intent(in), dimension(sz,sz) :: A
      double precision, intent(in), dimension(sz) :: Y
      double precision, intent(out), dimension(sz) :: X
    end subroutine simeq
  end interface

  print *, 'time_simeq'
  call date_and_time(date=start_date, time=start_time);
  print *, 'start_date=', start_date(5:6), '/', start_date(7:8), '/', start_date(1:4)
  print *, 'start_time=', start_time(1:2), ':', start_time(3:4), ':', start_time(5:10)
  call system_clock(count=clock_count, count_rate=clock_rate, count_max=clock_max)
  print *, 'clock_count=', clock_count, ', clock_rate=', clock_rate, ', clock_max=', clock_max
  start = (1.0D0+clock_count)/clock_rate
  
  n = 4
  do while(n<=sz)
    call init_matrix(n, sz, A, Y)
    call system_clock(count=clock_count)
    now = (1.0D0+clock_count)/clock_rate
    print *, 'matrix initialized at ', now-start
    call simeq(n, sz, A, Y, X)
    call system_clock(count=clock_count)
    next = (1.0D0+clock_count)/clock_rate
    total = 0.0
    do i=1,n
      t = 0.0
      do j=1,n
        t = t + A(i,j)*X(j)
      end do ! j
      total = total + abs(t-Y(i))
      if(n<5) then
        print *, 'X(', i, ')=', X(i), 't=', t, 'Y(', i, ')=', Y(i)
      end if
    end do ! i
    print *, 'n=', n, ', total error=', total, ', simeq took ', next-now, ' seconds'
    n = n+n
  end do ! n
end program time_simeq