C$TEST MFTE
C TO RUN AS A MAIN PROGRAM REMOVE NEXT LINE
      SUBROUTINE MFTE
C***********************************************************************
C
C  EXAMPLE OF USE OF THE PORT PROGRAM MFTCC
C
C***********************************************************************
      REAL A(25,25,25),B(25,25,25),T(50)
      REAL AA(25,25,25),BB(25,25,25)
      REAL RSTAK(1262)
      REAL SUM,FN1
      INTEGER IFX(25)
      INTEGER I,J,K,L,IFB
      DOUBLE PRECISION DSTAK(631)
      COMMON /CSTAK/DSTAK
      EQUIVALENCE (RSTAK(1),DSTAK(1))
C
C  N IS THE DIMENSION OF THE COMPLEX CUBE
C
      N = 25
C
C  GET ONE PLANE OF WORKSPACE
C
      NWK = 2*N*N + 12
      CALL ISTKIN(NWK,3)
C
C  SET REAL AND IMAGINARY PARTS TO YOUR FAVORITE VALUES HERE
C  AA AND BB ARE COPIES TO COMPARE WITH UNNORMALIZED OUTPUT
C  THIS EXAMPLE USES RANDOM DATA, SEE THE UTILITY CHAPTER OF PORT 3.
C
      DO 1 K = 1,N
      DO 1 J = 1,N
      DO 1 I = 1,N
         A(I,J,K)  = UNI(0)
         AA(I,J,K) = A(I,J,K)
         B(I,J,K)  = UNI(0)
         BB(I,J,K) = B(I,J,K)
 1    CONTINUE
C
C  INITIALIZE TRIGONOMETRIC TABLES AND FACTOR N
C
      CALL MFTCI(N,IFX,T)
C
C  OUTER LOOP COMPUTES A FORWARD (SIGN=1.0), THEN BACKWARD (SIGN=-1.0)
C  TRANSFORM. FIRST X, THEN Y, THEN Z.
C
      SIGN = 1.0E0
      N2   = N*N
      NT   = N
C
C
      DO 2 IFB = 1,2
C  X-DIRECTION TRANSFORMS FOR EACH X-Y PLANE
         DO 3 L = 1,N
            CALL MFTCC(NT,NT,A(1,1,L),B(1,1,L),1,NT,
     *                       A(1,1,L),B(1,1,L),1,NT,IFX,T,SIGN)
 3       CONTINUE
C  Y-DIRECTION TRANSFORMS FOR EACH X-Y PLANE
         DO 4 L = 1,N
            CALL MFTCC(NT,NT,A(1,1,L),B(1,1,L),NT,1,
     *                       A(1,1,L),B(1,1,L),NT,1,IFX,T,SIGN)
 4       CONTINUE
C  Z-DIRECTION TRANSFORMS FOR EACH Y-Z PLANE
         DO 5 L = 1,N
            CALL MFTCC(NT,NT,A(L,1,1),B(L,1,1),NT,N2,
     *                       A(L,1,1),B(L,1,1),NT,N2,IFX,T,SIGN)
 5       CONTINUE
C
         SIGN = -1.0E0
 2    CONTINUE
C
C
C  COMPARE INPUT TO UNNORMALIZED OUTPUT FROM FORWARD/BACKWARD FFT
C
      FN1 = 1.0E0/FLOAT(N*N*N)
      DO 6 K = 1,N
      DO 6 J = 1,N
      DO 6 I = 1,N
         AA(I,J,K) = AA(I,J,K) - FN1*A(I,J,K)
         BB(I,J,K) = BB(I,J,K) - FN1*B(I,J,K)
 6    CONTINUE
C
C  ERR IS THE RMS ERROR, SDOT COMPUTES THE SUM OF SQUARES, SEE THE
C  LINEAR ALGEBRA CHAPTER OF PORT 3
C
      ERR = SDOT(N2,AA,1,AA,1) + SDOT(N2,BB,1,BB,1)
      ERR = SQRT(FN1*ERR)
C
C  PRINT RESULTS
C
      IWRITE = I1MACH(2)
      WRITE(IWRITE,1000)N,N,N,ERR
 1000 FORMAT(1X,18H FOR LATTICE SIZE ,I3,2(2H X,I3),9H ERROR = ,E11.5)
      STOP
      END

.