! time_crout.f90  time solve AX=Y check with AA YY X_soln   

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

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

  if(n<5) then
    print *, 'init_matrix n=', n
  end if
  do i=1,n ! random A and X_soln
    X_soln(i) = udrnrt();
    do j=1,n
      A(i,j) = udrnrt()
      AA(i,j) = A(i,j)
    end do ! j
  end do ! i
  do i=1,n ! compute Y based on A and X_soln
    sum = 0.0
    do j=1,n
      sum = sum + A(i,j)*X_soln(j)
    end do ! j
    Y(i) = sum
    YY(i) = Y(i)
  end do ! i
  if(n<5) then
    do i=1,n
      do j=1,n
        print *, 'A(', i, ',', j, ')=', A(i,j)
      end do ! j
      print *, 'X_soln(', i, ')=', X_soln(i), ', Y(', i, ')=', Y(i)
    end do ! i
  end if
end subroutine init_matrix


program time_crout
  implicit none 
  integer, parameter :: sz = 2048 
  double precision, dimension(sz,sz) :: A
  double precision, dimension(sz,sz) :: AA
  double precision, dimension(sz) :: X
  double precision, dimension(sz) :: X_soln
  double precision, dimension(sz) :: Y
  double precision, dimension(sz) :: YY
  double precision :: sum, sumerr, maxerr
  double precision :: start, now, next
  integer :: i, j, n

  character*20 :: start_date
  character*20 :: start_time
  integer :: clock_count, clock_rate, clock_max

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

  interface
    subroutine init_matrix(n, sz, A, AA, X_soln, Y, YY)
      integer, intent(in) :: n
      integer, intent(in) :: sz
      double precision, dimension(sz,sz), intent(out) :: A
      double precision, dimension(sz,sz), intent(out) :: AA
      double precision, dimension(sz), intent(out) :: X_soln
      double precision, dimension(sz), intent(out) :: Y
      double precision, dimension(sz), intent(out) :: YY
    end subroutine init_matrix
  end interface

  print *, 'time_crout.f90 running'
  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<=2048)
    call init_matrix(n, sz, A, AA, X_soln, Y, YY)
    call system_clock(count=clock_count)
    now = (1.0D0+clock_count)/clock_rate
    print *, 'matrix initialized in', now-start, ' seconds'
    call crout(n, sz, A, Y, X)
    call system_clock(count=clock_count)
    next = (1.0D0+clock_count)/clock_rate
    sumerr = 0.0
    maxerr = 0.0
    do i=1,n
      sum = 0.0
      do j=1,n
        sum = sum + AA(i,j)*X(j)
      end do ! j
      sumerr = sumerr + dabs(sum-YY(i))
      maxerr = dmax1(maxerr, dabs(X_soln(i)-X(i)))
      if(n<5) then
        print *,'X(', i, ')=', X(i), ', err=', X_soln(i)-X(i)
      end if
    end do ! i
    print *,'n=', n, ' sumerr=', sumerr, ', maxerr=', maxerr
    print *, 'crout took ', next-now, ' seconds'
    n = n+n
  end do ! n
end program time_crout