*     C06FFF Example Program Text
*     Mark 14 Revised. NAG Copyright 1989.
*     .. Parameters ..
      INTEGER          NDIM, NMAX, LWORK
      PARAMETER        (NDIM=2,NMAX=96,LWORK=96)
      INTEGER          NIN, NOUT
      PARAMETER        (NIN=5,NOUT=6)
*     .. Local Scalars ..
      INTEGER          IFAIL, L, N
*     .. Local Arrays ..
      DOUBLE PRECISION WORK(LWORK), X(NMAX), Y(NMAX)
      INTEGER          ND(NDIM)
*     .. External Subroutines ..
      EXTERNAL         C06FFF, C06GCF, READXY, WRITXY
*     .. Executable Statements ..
      WRITE (NOUT,*) 'C06FFF Example Program Results'
*     Skip heading in data file
      READ (NIN,*)
   20 READ (NIN,*,END=40) ND(1), ND(2), L
      N = ND(1)*ND(2)
      IF (N.GE.1 .AND. N.LE.NMAX) THEN
         CALL READXY(NIN,X,Y,ND(1),ND(2))
         WRITE (NOUT,*)
         WRITE (NOUT,*) 'Original data'
         CALL WRITXY(NOUT,X,Y,ND(1),ND(2))
         IFAIL = 1
*
*        Compute transform
         CALL C06FFF(NDIM,L,ND,N,X,Y,WORK,LWORK,IFAIL)
*
         IF (IFAIL.EQ.0) THEN
            WRITE (NOUT,*)
            WRITE (NOUT,99999) 'Discrete Fourier transform of variable '
     +        , L
            CALL WRITXY(NOUT,X,Y,ND(1),ND(2))
*
*           Compute inverse transform
            CALL C06GCF(Y,N,IFAIL)
            CALL C06FFF(NDIM,L,ND,N,X,Y,WORK,LWORK,IFAIL)
            CALL C06GCF(Y,N,IFAIL)
*
            WRITE (NOUT,*)
            WRITE (NOUT,*)
     +        'Original sequence as restored by inverse transform'
            CALL WRITXY(NOUT,X,Y,ND(1),ND(2))
            GO TO 20
         ELSE
            WRITE (NOUT,*)
            WRITE (NOUT,99998) ' ** C06FFF returned with IFAIL = ',
     +        IFAIL
         END IF
      ELSE
         WRITE (NOUT,*) 'Invalid value of N'
      END IF
   40 CONTINUE
*
99999 FORMAT (1X,A,I1)
99998 FORMAT (1X,A,I5)
      END
*
      SUBROUTINE READXY(NIN,X,Y,N1,N2)
*     Read 2-dimensional complex data
*     .. Scalar Arguments ..
      INTEGER          N1, N2, NIN
*     .. Array Arguments ..
      DOUBLE PRECISION X(N1,N2), Y(N1,N2)
*     .. Local Scalars ..
      INTEGER          I, J
*     .. Executable Statements ..
      DO 20 I = 1, N1
         READ (NIN,*) (X(I,J),J=1,N2)
         READ (NIN,*) (Y(I,J),J=1,N2)
   20 CONTINUE
      RETURN
      END
*
      SUBROUTINE WRITXY(NOUT,X,Y,N1,N2)
*     Print 2-dimensional complex data
*     .. Scalar Arguments ..
      INTEGER          N1, N2, NOUT
*     .. Array Arguments ..
      DOUBLE PRECISION X(N1,N2), Y(N1,N2)
*     .. Local Scalars ..
      INTEGER          I, J
*     .. Executable Statements ..
      DO 20 I = 1, N1
         WRITE (NOUT,*)
         WRITE (NOUT,99999) 'Real ', (X(I,J),J=1,N2)
         WRITE (NOUT,99999) 'Imag ', (Y(I,J),J=1,N2)
   20 CONTINUE
      RETURN
*
99999 FORMAT (1X,A,7F10.3,/(6X,7F10.3))
      END