PROGRAM Main
    IMPLICIT NONE
    INTEGER, PARAMETER :: dp = SELECTED_REAL_KIND(12)
    COMPLEX(KIND=dp), ALLOCATABLE :: cw(:,:), cx(:,:), cy(:,:), w(:)
    INTEGER :: size, n
    CHARACTER(LEN=100) :: filename

! Read the file name from the argument and read the size.
    IF (COMMAND_ARGUMENT_COUNT() /= 1) THEN
        WRITE (*,'("Wrong number of arguments")')
        STOP
    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
    END IF

! Convert the size to the largest power of two that will fit, and
! set up the data.

    n = 1
    DO WHILE (n**2 <= size**2/8)
        n = 2*n
    END DO
    size = n
    if (n**2 <= size**2/4) n = 2*n
    WRITE (*,'("Actual FFT sizes: ",I0," and ",I0,"x",I0)') &
        size**2, size, size
    ALLOCATE(cw(size,size),cx(size,size),cy(size,size),w(size**2))
    READ (42) cw
    CLOSE(42)

! Time and check the FFTs in single-dimensional mode.

    CALL Cffti(w,size**2)
    cx = cw
    CALL Times('')
    CALL Cfftf(1,size**2,cx,cy,w,.True.)
    CALL Cfftb(1,size**2,cx,cy,w,.True.)
    CALL Times('1-D FFT ignoring columns')
    WRITE (*,'("Relative error:",ES9.2)') &
        MAXVAL(ABS(cx/size**2-cw))/MAXVAL(ABS(cw))
    cx = cw
    CALL Times('')
    CALL Cfftf(size**2,1,cx,cy,w,.False.)
    CALL Cfftb(size**2,1,cx,cy,w,.False.)
    CALL Times('1-D FFT ignoring rows')
    WRITE (*,'("Relative error:",ES9.2)') &
        MAXVAL(ABS(cx/size**2-cw))/MAXVAL(ABS(cw))

! Time and check the FFTs using TRANSPOSE.

    CALL Cffti(w,size)
    cx = cw
    CALL Times('')
    CALL Cfftf(size,size,cx,cy,w,.True.)
    CALL My_transpose(cx,size)
    CALL Cfftf(size,size,cx,cy,w,.True.)
    CALL My_transpose(cx,size)
    CALL Cfftb(size,size,cx,cy,w,.True.)
    CALL My_transpose(cx,size)
    CALL Cfftb(size,size,cx,cy,w,.True.)
    CALL My_transpose(cx,size)
    CALL Times('2-D FFT  transposing in place')
    WRITE (*,'("Relative error:",ES9.2)') &
        MAXVAL(ABS(cx/size**2-cw))/MAXVAL(ABS(cw))

! Time and check the FFTs without reordering them.

    CALL Cffti(w,size)
    cx = cw
    CALL Times('')
    CALL Cfftf(size,size,cx,cy,w,.True.)
    CALL Cfftf(size,size,cx,cy,w,.False.)
    CALL Cfftb(size,size,cx,cy,w,.True.)
    CALL Cfftb(size,size,cx,cy,w,.False.)
    CALL Times('2-D FFT without transposing')
    WRITE (*,'("Relative error:",ES9.2)') &
        MAXVAL(ABS(cx/size**2-cw))/MAXVAL(ABS(cw))

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," = ",F0.2," seconds, CPU = ",F0.2," seconds")') &
                which,wall,cpu
        END IF
    END SUBROUTINE Times



! This code is transposition in place, using the obvious code for a
! square matrix, and gives lots of scope for tuning.

    SUBROUTINE My_transpose (array, n)
        INTEGER, INTENT(IN) :: n
        COMPLEX(KIND=dp), INTENT(INOUT) :: array(n,n)
        INTEGER :: i, j
        COMPLEX(KIND=dp) :: x

        DO i = 2,size
            DO j = 1,i-1
                x = array(i,j)
                array(i,j) = array(j,i)
                array(j,i) = x
            END DO
        END DO
    END Subroutine My_transpose



! The following code started off as Schwarztrauber's, but has been cut
! down to powers of two only, and converted to Fortran 90.

    SUBROUTINE Cffti (wa, n)
! Generate the array of constants.
        IMPLICIT NONE
        INTEGER, PARAMETER :: dp = SELECTED_REAL_KIND(12)
        INTEGER, INTENT(IN) :: n
        COMPLEX(KIND=dp), INTENT(OUT) :: wa(:)
        INTEGER :: i, k, m
        REAL(KIND=dp) :: f
     
        f = (4.0_dp*ATAN(1.0_dp))/n
        k = 1
        m = n
        DO WHILE (m > 1)
            IF (MOD(m,2) /= 0) THEN
                WRITE (*,*) 'The vector size is not a power of two:', m
                STOP
            END IF
            m = m/2
            DO i = 0, m-1
                wa(k+i) = CMPLX(COS(f*i),SIN(f*i),KIND=dp)
            END DO
            k = k+m
            f = f+f
        END DO
    END SUBROUTINE Cffti



    SUBROUTINE Cfftf (m, n, c, ch, wa, flag)
! The normal FFT.
        IMPLICIT NONE
        INTEGER, PARAMETER :: dp = SELECTED_REAL_KIND(12)
        INTEGER, INTENT(IN) :: m, n
        COMPLEX(KIND=dp), INTENT(INOUT) :: ch(m,n), c(m,n)
        COMPLEX(KIND=dp), INTENT(IN) :: wa(:)
        LOGICAL, INTENT(IN) :: flag
        INTEGER :: mn, k, l, iw
        LOGICAL :: na

        na = .True.
        iw = 1
        l = 1
        IF (flag) THEN
            mn = n
        ELSE
            mn = m
        END IF
        k = mn
        DO WHILE (l < mn)
            k = k/2
            IF (na) THEN
                IF (flag) THEN
                    CALL Passf_x(m,k,l,c,ch,wa(iw:iw+k-1))
                ELSE
                    CALL Passf_y(k,l,n,c,ch,wa(iw:iw+k-1))
                END IF
            ELSE
                IF (flag) THEN
                    CALL Passf_x(m,k,l,ch,c,wa(iw:iw+k-1))
                ELSE
                    CALL Passf_y(k,l,n,ch,c,wa(iw:iw+k-1))
                END IF
            END IF
            na = .NOT. na
            l = 2*l
            iw = iw+k
        END DO
        IF (.NOT. na) c = ch
    END SUBROUTINE Cfftf



    SUBROUTINE Cfftb (m, n, c, ch, wa, flag)
! The inverse FFT.
        IMPLICIT NONE
        INTEGER, PARAMETER :: dp = SELECTED_REAL_KIND(12)
        INTEGER, INTENT(IN) :: m, n
        COMPLEX(KIND=dp), INTENT(INOUT) :: ch(m,n), c(m,n)
        COMPLEX(KIND=dp), INTENT(IN) :: wa(:)
        LOGICAL, INTENT(IN) :: flag
        INTEGER :: mn, k, l, iw
        LOGICAL :: na

        na = .True.
        iw = 1
        l = 1
        IF (flag) THEN
            mn = n
        ELSE
            mn = m
        END IF
        k = mn
        DO WHILE (l < mn)
           k = k/2
            IF (na) THEN
                IF (flag) THEN
                    CALL Passb_x(m,k,l,c,ch,wa(iw:iw+k-1))
                ELSE
                    CALL Passb_y(k,l,n,c,ch,wa(iw:iw+k-1))
                END IF
            ELSE
                IF (flag) THEN
                    CALL Passb_x(m,k,l,ch,c,wa(iw:iw+k-1))
                ELSE
                    CALL Passb_y(k,l,n,ch,c,wa(iw:iw+k-1))
                END IF
            END IF
            na = .NOT. na
            l = 2*l
            iw = iw+k
        END DO
        IF (.NOT. na) c = ch
    END SUBROUTINE Cfftb



    SUBROUTINE Passf_x (m, ido, l, cc, ch, wa)
! One pass of the normal FFT.
        IMPLICIT NONE
        INTEGER, PARAMETER :: dp = SELECTED_REAL_KIND(12)
        INTEGER, INTENT(IN) :: m, ido, l
        COMPLEX(KIND=dp), INTENT(IN) :: cc(m,ido,2,l), wa(:)
        COMPLEX(KIND=dp), INTENT(OUT) :: ch(m,ido,l,2)
        INTEGER :: i, k

        DO k = 1, l
            DO i = 1, ido
                ch(:,i,k,1) = cc(:,i,1,k)+cc(:,i,2,k)
                ch(:,i,k,2) = CONJG(wa(i))*(cc(:,i,1,k)-cc(:,i,2,k))
            END DO
        END DO
    END SUBROUTINE Passf_x



    SUBROUTINE Passf_y (ido, l, m, cc, ch, wa)
! One pass of the normal FFT.
        IMPLICIT NONE
        INTEGER, PARAMETER :: dp = SELECTED_REAL_KIND(12)
        INTEGER, INTENT(IN) :: ido, l, m
        COMPLEX(KIND=dp), INTENT(IN) :: cc(ido,2,l,m), wa(:)
        COMPLEX(KIND=dp), INTENT(OUT) :: ch(ido,l,2,m)
        INTEGER :: i, k

        DO k = 1, l
            DO i = 1, ido
                ch(i,k,1,:) = cc(i,1,k,:)+cc(i,2,k,:)
                ch(i,k,2,:) = CONJG(wa(i))*(cc(i,1,k,:)-cc(i,2,k,:))
            END DO
        END DO
    END SUBROUTINE Passf_y



    SUBROUTINE Passb_x (m, ido, l, cc, ch, wa)
! One pass of the inverse FFT.
        IMPLICIT NONE
        INTEGER, PARAMETER :: dp = SELECTED_REAL_KIND(12)
        INTEGER, INTENT(IN) :: m, ido, l
        COMPLEX(KIND=dp), INTENT(IN) :: cc(m,ido,2,l), wa(:)
        COMPLEX(KIND=dp), INTENT(OUT) :: ch(m,ido,l,2)
        INTEGER :: i, k

        DO k = 1, l
            DO i = 1, ido
                ch(:,i,k,1) = cc(:,i,1,k)+cc(:,i,2,k)
                ch(:,i,k,2) = wa(i)*(cc(:,i,1,k)-cc(:,i,2,k))
            END DO
        END DO
    END SUBROUTINE Passb_x



    SUBROUTINE Passb_y (ido, l, m, cc, ch, wa)
! One pass of the inverse FFT.
        IMPLICIT NONE
        INTEGER, PARAMETER :: dp = SELECTED_REAL_KIND(12)
        INTEGER, INTENT(IN) :: ido, l, m
        COMPLEX(KIND=dp), INTENT(IN) :: cc(ido,2,l,m), wa(:)
        COMPLEX(KIND=dp), INTENT(OUT) :: ch(ido,l,2,m)
        INTEGER :: i, k

        DO k = 1, l
            DO i = 1, ido
                ch(i,k,1,:) = cc(i,1,k,:)+cc(i,2,k,:)
                ch(i,k,2,:) = wa(i)*(cc(i,1,k,:)-cc(i,2,k,:))
            END DO
        END DO
    END SUBROUTINE Passb_y



END PROGRAM Main