PROGRAM mul
  USE naivmat
  USE bettmat
  USE dotmat
  USE iso_fortran_env, only: error_unit
  IMPLICIT none

  integer, parameter :: seed = 123456
  real :: time
  integer :: dim = 10, multype = 1, real_kind, stat
  character(5) :: kind_arg

  IF (command_argument_count() < 1) THEN
     write (error_unit, *) "Usage: mul KIND", char(10), &
                           "where KIND is one of: 4, 8, 16"
     STOP
  END IF

  call get_command_argument(1, kind_arg)
  
  read (kind_arg, *, iostat = stat) real_kind

  IF (stat .ne. 0) THEN
     write (error_unit, *) "Couldn't parse kind number argument"
     STOP
  END IF     
  
  DO WHILE (dim < 2000)

     SELECT CASE(real_kind)
        CASE (4,1)
           time = measure_4(dim)
        CASE (8,2)
           time = measure_8(dim)
        CASE (16,3)
           time = measure_16(dim)
        CASE default
           write (error_unit, *) "wrong kind for real: ", real_kind
           STOP

     END SELECT
        
     print '(I11," ",ES11.5)', dim, time

     dim = dim * 2
     
  END DO

CONTAINS

  SUBROUTINE init_random_seed()
    integer :: i, n
    
    call random_seed(size = n)
    
    call random_seed(put = (/ ((seed + i) * 37, i = 1, n) /))
    
  END SUBROUTINE init_random_seed
  

  real FUNCTION measure_4(dim) result(time)
    integer, intent(in) :: dim
    real(kind=4), dimension(:,:), allocatable :: mat1, mat2, res
    real :: start, end

    call init_random_seed()

    allocate(mat1(dim,dim))
    allocate(mat2(dim,dim))
    allocate(res(dim,dim))

    call random_number(mat1)
    call random_number(mat2)

    call cpu_time(start)

    SELECT CASE(multype)
       CASE (1)
          res = naivmull(mat1, mat2)
       CASE (2)
          res = bettmull(mat1, mat2)
       CASE (3)
          res = dotmull(mat1, mat2)
       CASE default
          res = matmul(mat1, mat2)
       
    END SELECT

    call cpu_time(end)
    
    time = end - start
    
  END FUNCTION measure_4

  
  real FUNCTION measure_8(dim) result(time)
    integer, intent(in) :: dim
    real(kind=8), dimension(:,:), allocatable :: mat1, mat2, res
    real :: start, end

    call init_random_seed()

    allocate(mat1(dim,dim))
    allocate(mat2(dim,dim))
    allocate(res(dim,dim))

    call random_number(mat1)
    call random_number(mat2)

    call cpu_time(start)

    SELECT CASE(multype)
       CASE (1)
          res = naivmull(mat1, mat2)
       CASE (2)
          res = bettmull(mat1, mat2)
       CASE (3)
          res = dotmull(mat1, mat2)
       CASE default
          res = matmul(mat1, mat2)
       
    END SELECT

    call cpu_time(end)
    
    time = end - start
    
  END FUNCTION measure_8


  real FUNCTION measure_16(dim) result(time)
    integer, intent(in) :: dim
    real(kind=16), dimension(:,:), allocatable :: mat1, mat2, res
    real :: start, end

    call init_random_seed()

    allocate(mat1(dim,dim))
    allocate(mat2(dim,dim))
    allocate(res(dim,dim))

    call random_number(mat1)
    call random_number(mat2)

    call cpu_time(start)

    SELECT CASE(multype)
       CASE (1)
          res = naivmull(mat1, mat2)
       CASE (2)
          res = bettmull(mat1, mat2)
       CASE (3)
          res = dotmull(mat1, mat2)
       CASE default
          res = matmul(mat1, mat2)
       
    END SELECT

    call cpu_time(end)
    
    time = end - start
    
  END FUNCTION measure_16

END PROGRAM mul