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

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

  IF (command_argument_count() < 2) THEN
     call print_usage()
     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, '(A)') "Couldn't parse kind number argument"
     call print_usage()
     STOP
  END IF

  call get_command_argument(2, impl_arg)

  IF      (trim(impl_arg) .eq. "naiv") THEN
     multype = 1
  ELSE IF (trim(impl_arg) .eq. "bett") THEN
     multype = 2
  ELSE IF (trim(impl_arg) .eq. "dot")  THEN
     multype = 3
  ELSE IF (trim(impl_arg) .eq. "mat")  THEN
     multype = 4
  ELSE IF (trim(impl_arg) .eq. "bett2")  THEN
     multype = 5
  ELSE IF (trim(impl_arg) .eq. "block")  THEN
     multype = 6
  ELSE
     write (error_unit, '(A)') "Unrecognized implementation argument"
     call print_usage()
     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, '(A,I6)')            &
                "wrong kind for real:", real_kind
           STOP

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

     dim = dim * 2
     
  END DO

CONTAINS

  SUBROUTINE print_usage()
    
    write (*, '(A)')                                             &
          "Usage: mull KIND IMPLEMENTATION" // char(10) //       &
          "where KIND is one of: 4, 8, 16" // char(10) //        &
          "      IMPLEMENTATION is one of: " //                  &
          "naiv, bett, dot, mat, bett2, block"
    
  END SUBROUTINE print_usage

  
  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 (4)
          res = matmul(mat1, mat2)
       CASE (5)
          res = bett2mull(mat1, mat2)
       CASE default
          res = blockmull(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 (4)
          res = matmul(mat1, mat2)
       CASE (5)
          res = bett2mull(mat1, mat2)
       CASE default
          res = blockmull(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 (4)
          res = matmul(mat1, mat2)
       CASE (5)
          res = bett2mull(mat1, mat2)
       CASE default
          res = blockmull(mat1, mat2)
       
    END SELECT

    call cpu_time(end)

    time = end - start
    
  END FUNCTION measure_16

END PROGRAM mul