PROGRAM Main
    IMPLICIT NONE
    INTEGER, PARAMETER :: dp = SELECTED_REAL_KIND(12)
    REAL(KIND=dp), ALLOCATABLE :: data(:)
    REAL(KIND=dp) :: a, b, c, d, e, f
    INTEGER :: size, n
    CHARACTER(LEN=100) :: filename

! Read the file name from the argument and set up the data.
    IF (COMMAND_ARGUMENT_COUNT() /= 1) THEN
        WRITE (*,'("Wrong number of arguments")')
        STOP 1
    END IF
    CALL GET_COMMAND_ARGUMENT(1,filename)
    OPEN (42,FILE=filename,ACTION='READ',FORM='UNFORMATTED')
    READ (42) size
    IF (size < 1 .OR. size > 1000000) THEN
        WRITE (*,'("Invalid value of size: ",I0)') size
        STOP 1
    END IF
    size = size*size
    ALLOCATE(data(size))
    READ (42) data
    CLOSE(42)

! Time and print the use of Fortran's intrinsics.  The loop is because
! it's unreasonably quick otherwise.
    CALL Times('')
    DO n = 1,250
        a = MINVAL(data)
        b = MAXVAL(data)
        c = SUM(data)/size
        d = SQRT(SUM((data-c)**2)/size)
        e = SUM((data-c)**3)/(size*d**3)
        f = SUM((data-c)**4)/(size*d**4)
    END DO
    CALL Times('Intrinsics')
    WRITE (*,'("Min=",F0.2,", max=",F0.2,", mean=",F0.2,", SE=",F0.3,  &
            & ", skew=",F0.2,", kurt=",F0.2)') &
        a,b,c,d,e,f

! Using the level 1 BLAS is very rarely worth the effort, so we don't
! bother.

! Time and print the use of hand-coded calculations.  The loop is because
! it's unreasonably quick otherwise.
    CALL Times('')
    DO n = 1,250
        CALL Stats(data,a,b,c,d,e,f)
    END DO
    CALL Times('Coded')
    WRITE (*,'("Min=",F0.2,", max=",F0.2,", mean=",F0.2,", SE=",F0.3,  &
            & ", skew=",F0.2,", kurt=",F0.2)') &
        a,b,c,d,e,f

CONTAINS

    SUBROUTINE Times (which)
! If which is not empty, print the times since the previous call.
        CHARACTER(LEN=*), INTENT(IN) :: which
        REAL(KIND=dp), SAVE :: last_wall = 0.0_dp, last_cpu = 0.0_dp
        REAL(KIND=dp) :: wall, cpu
        INTEGER :: m, n

        wall = last_wall
        cpu = last_cpu
        CALL SYSTEM_CLOCK(m,n)
        last_wall = m/REAL(n,KIND=dp)
        CALL CPU_TIME(last_cpu)
        IF (LEN(which) > 0) THEN
            wall = last_wall-wall
            cpu = last_cpu-cpu
            WRITE (*,'(A," time = ",F0.2," seconds, CPU = ",F0.2," seconds")') &
                which,wall,cpu
        END IF
    END SUBROUTINE Times

    SUBROUTINE Stats (data, xmin, xmax, mean, se, skew, kurt)
! Calculate minimum, maximum, mean, standard error, skewness and
! kurtosis.
        REAL(KIND=dp), INTENT(IN) :: data(:)
        REAL(KIND=dp), INTENT(OUT) :: xmin, xmax, mean, se, skew, kurt
        REAL(KIND=dp) :: x, y
        INTEGER :: size, i

        size = UBOUND(data,1)
        xmin = HUGE(xmin)
        xmax = -HUGE(xmax)
        mean = 0.0_dp
        se = 0.0_dp
        skew = 0.0_dp
        kurt = 0.0_dp
        DO i = 1,size
            xmin = MIN(xmin,data(i))
            xmax = MAX(xmax,data(i))
            mean = mean+data(i)
        END DO
        mean = mean/size
        DO i = 1,size
            x = data(i)-mean
            y = x**2
            se = se+y
            skew = skew+x*y
            kurt = kurt+y*y
        END DO
        se = sqrt(se/size)
        skew = skew/(size*se**3)
        kurt = kurt/(size*se**4)
    END SUBROUTINE Stats


END PROGRAM Main
