      program example

C To use this distribution, clip off the main program and the 'BUILD.COM'
C file.  The remainder will be the FT code.  You can fragment it further
C if you want to do so. BTW: The BUILD.COM is for a VAX.  I know. Yuck.
C But this compiles just fine on Ultrix, and most everything else, so don't
C get your jeans in an uproar.

      real        a(10200)
      complex     c(1)
      equivalence (a,c)

C Forward 2d transform

      N = 100
      M = 100

      call real2f(a,N,M, +1)

C This has created a 2-d transform, 102 x 100 reals (51 x 100 complex)
C for an N x M array the transform is stored as a (N/2+1) x M complex
C In X orders go from 0 to N/2 (0 to 50 in our example here: 50 is the
C Nyquist order. In Y the orders go from 0 to M/2, and then down from
C -M/2 to -1.  You can reorder the array to get the zero order in the
C middle, but I do not use it like that, and you are probably better
C off leaving well enough alone adn indexing the array with something
C like 'val at order (i,j) = c(i+1+(N+2)*mod(j+M,M))'

C ** NOTE ** N returns as N/2+1 i.e. the X dimension of the complex array
C generated

C Reverse Transform

      call real2f(a,N,M, -1)

      call exit
      end
C --------------------------- Clip Here -----------------------------
$for/ext CMPLFT.FOR
$for/ext DIPRP.FOR
$for/ext HERMFT.FOR
$for/ext MDFTKD.FOR
$for/ext MDP_ERROR.FOR	! <-- This is a dummy routine
$for/ext R2CFTK.FOR
$for/ext R3CFTK.FOR
$for/ext R4CFTK.FOR
$for/ext R5CFTK.FOR
$for/ext R8CFTK.FOR
$for/ext REAL2F.FOR
$for/ext REAL3F.FOR
$for/ext REALFT.FOR
$for/ext RPCFTK.FOR
$for/ext SRFP.FOR
$library/create ten_eyke_fft.olb *.obj
$Delete/noconfirm *.obj;*
$write sys$output "FFT Library complete"
$exit
$! ------------------------- Clip here _________________________________
      SUBROUTINE CMPLFT (X, Y, N, D)
      INTEGER N
      INTEGER D(5)
      REAL X(1), Y(1)
C
      LOGICAL ERROR
      INTEGER P MAX, P SYM, TWO GRP
      INTEGER FACTOR(15), SYM(15), UN SYM(15)
C
      P MAX = 41
      TWO GRP = 8
C
      IF (N .LE. 1) GO TO 100
      CALL S R FP (N, P MAX, TWO GRP, FACTOR, SYM, P SYM, UN SYM, ERROR)
      IF (ERROR) GO TO 200
C
      CALL MD FT KD (N, FACTOR, D, X, Y)
      CALL D IP RP (N, SYM, P SYM, UN SYM, D, X, Y)
C
  100 CONTINUE
      RETURN
C
  200 CONTINUE
      call mdp_error(13,0)
      return
C
c  300 FORMAT ('0INVALID NUMBER OF POINTS FOR CMPL FT.  N =', I10//)
C
      END
      SUBROUTINE DIPRP (PTS,SYM,P SYM,UN SYM,DIM,X,Y)
C     DOUBLE IN PLACE REORDERING PROGRAMME
C
      INTEGER PTS,P SYM
      REAL X (4000), Y (4000)
      INTEGER SYM(15), UN SYM(15), DIM(5)
      INTEGER*4 K
C
      REAL T
      LOGICAL ONE MOD
      INTEGER MODULO (14)
      INTEGER DK,JJ,KK,LK,MODS,MULT,P UN SYM,TEST
      INTEGER NT, SEP, P, P0, P1, P2, P3, P4, P5, SIZE
C
      INTEGER S (14), U (14)
      INTEGER A,B,C,D,E,F,G,H,I,J,L,M,N
      INTEGER    BS,CS,DS,ES,FS,GS,HS,IS,JS,KS,LS,MS,NS
      INTEGER AL,BL,CL,DL,EL,FL,GL,HL,IL,JL,KL,LL,ML,NL
      EQUIVALENCE           (AL,U(1)),(BS,S(2)),(BL,U(2))
      EQUIVALENCE (CS,S(3)),(CL,U(3)),(DS,S(4)),(DL,U(4))
      EQUIVALENCE (ES,S(5)),(EL,U(5)),(FS,S(6)),(FL,U(6))
      EQUIVALENCE (GS,S(7)),(GL,U(7)),(HS,S(8)),(HL,U(8))
      EQUIVALENCE (IS,S(9)),(IL,U(9)),(JS,S(10)),(JL,U(10))
      EQUIVALENCE (KS,S(11)),(KL,U(11)),(LS,S(12)),(LL,U(12))
      EQUIVALENCE (MS,S(13)),(ML,U(13)),(NS,S(14)),(NL,U(14))
      PARAMETER NEST = 14
C
      NT = DIM(1)
      SEP = DIM(2)
      P2 = DIM(3)
      SIZE = DIM(4) - 1
      P4 = DIM(5)
      IF (SYM(1).EQ.0) GO TO 500
      DO 100 J=1,NEST
      U(J)=1
      S(J)=1
  100 CONTINUE
      N=PTS
      DO 200 J=1,NEST
      IF (SYM(J).EQ.0) GO TO 300
      JJ=NEST+1-J
      U(JJ)=N
      S(JJ)=N/SYM(J)
      N=N/SYM(J)
  200 CONTINUE
  300 CONTINUE
C
      JJ=0
      DO 400 A=1,AL
      DO 400 B=A,BL,BS
      DO 400 C=B,CL,CS
      DO 400 D=C,DL,DS
      DO 400 E=D,EL,ES
      DO 400 F=E,FL,FS
      DO 400 G=F,GL,GS
	DO 400 H=G,HL,HS
	DO 400 I=H,IL,IS
	DO 400 J=I,JL,JS
	DO 400 KK=J,KL,KS
	DO 400 L=KK,LL,LS
	DO 400 M=L,ML,MS
	DO 400 N=M,NL,NS
      JJ=JJ+1
      IF (JJ.GE.N) GO TO 400
      DELTA = (N-JJ)*SEP
      P1 = (JJ-1)*SEP + 1
      DO 350 P0 = P1, NT, P2
      P3 = P0 + SIZE
      DO 350 P = P0, P3, P4
      P5 = P + DELTA
      T = X(P)
      X(P) = X(P5)
      X(P5) = T
      T = Y(P)
      Y(P) = Y(P5)
      Y(P5) = T
  350 CONTINUE
  400 CONTINUE
C
  500 CONTINUE
      IF (UN SYM(1).EQ.0) GO TO 1900
      P UN SYM=PTS/P SYM**2
      MULT=P UN SYM/UN SYM(1)
      TEST=(UN SYM(1)*UN SYM(2)-1)*MULT*P SYM
      LK=MULT
      DK=MULT
      DO 600 K=2,NEST
      IF (UN SYM(K).EQ.0) GO TO 700
      LK=LK*UN SYM(K-1)
      DK=DK/UN SYM(K)
      U(K)=(LK-DK)*P SYM
      MODS=K
  600 CONTINUE
  700 CONTINUE
      ONE MOD=MODS.LT.3
      IF (ONE MOD) GO TO 900
      DO 800 J=3,MODS
      JJ=MODS+3-J
      MODULO(JJ)=U(J)
  800 CONTINUE
  900 CONTINUE
      MODULO(2)=U(2)
      JL=(P UN SYM-3)*P SYM
      MS=P UN SYM*P SYM
C
      DO 1800 J=P SYM,JL,P SYM
      K=J
C
1000  K=K*MULT
      IF (ONE MOD) GO TO 1200
      DO 1100 I=3,MODS
1100  K=K-(K/MODULO(I))*MODULO(I)
C
1200  IF (K.GE.TEST) GO TO 1300
      K=K-(K/MODULO(2))*MODULO(2)
      GO TO 1400
1300  K=K-(K/MODULO(2))*MODULO(2)+MODULO(2)
1400  IF (K.LT.J) GO TO 1000
C
      IF (K.EQ.J) GO TO 1700
      DELTA = (K-J)*SEP
      DO 1600 L=1,P SYM
      DO 1500 M=L,PTS,MS
      P1 = (M+J-1)*SEP + 1
      DO 1500 P0 = P1, NT, P2
      P3 = P0 + SIZE
      DO 1500 JJ = P0, P3, P4
      KK = JJ + DELTA
      T=X(JJ)
      X(JJ)=X(KK)
      X(KK)=T
      T=Y(JJ)
      Y(JJ)=Y(KK)
      Y(KK)=T
 1500 CONTINUE
 1600 CONTINUE
 1700 CONTINUE
 1800 CONTINUE
C
 1900 CONTINUE
      RETURN
      END
      SUBROUTINE HERMFT (X, Y, N, DIM)
      INTEGER N
      INTEGER DIM(5)
      REAL X(1), Y(1)
C
      REAL A, B, C, D, E, F, ANGLE, CO, SI, TWO N, TWO PI
      INTEGER NT, D2, D3, D4, D5
      INTEGER I, J, K, N OVER 2, I0, I1, I2
      INTEGER K1
C
      TWO PI = 6.283185
      TWO N = FLOAT(2*N)
C
      NT = DIM(1)
      D2 = DIM(2)
      D3 = DIM(3)
      D4 = DIM(4) - 1
      D5 = DIM(5)
C
      DO 100 I0 = 1, NT, D3
      I1 = I0 + D4
      DO 100 I = I0, I1, D5
      A = X(I)
      B = Y(I)
      X(I) = A + B
      Y(I) = A - B
  100 CONTINUE
C
      N OVER 2 = N/2 + 1
      IF (N OVER 2 .LT. 2) GO TO 500
      DO 400 I0 = 2, N OVER 2
      ANGLE = TWO PI*FLOAT(I0-1)/TWO N
      CO = COS(ANGLE)
      SI = SIN(ANGLE)
      K = (N + 2 - 2*I0)*D2
      K1 = (I0 - 1)*D2 + 1
      DO 300 I1 = K1, NT, D3
      I2 = I1 + D4
      DO 200 I = I1, I2, D5
      J = I + K
      A = X(I) + X(J)
      B = X(I) - X(J)
      C = Y(I) + Y(J)
      D = Y(I) - Y(J)
      E = B*CO + C*SI
      F = B*SI - C*CO
      X(I) = A + F
      X(J) = A - F
      Y(I) = E + D
      Y(J) = E - D
  200 CONTINUE
  300 CONTINUE
  400 CONTINUE
C
      CALL CMPL FT (X, Y, N, DIM)
C
  500 CONTINUE
      RETURN
C
      END
      SUBROUTINE MDFTKD (N, FACTOR, DIM, X, Y)
C     MULTI-DIMENSIONAL COMPLEX FOURIER TRANSFORM KERNEL DRIVER
C
      INTEGER N
      INTEGER FACTOR (10), DIM(5)
      REAL X(10), Y(10)
C
      INTEGER F, M, P, R, S
C
      S = DIM(2)
      F = 0
      M = N
  100 CONTINUE
      F = F + 1
      P = FACTOR(F)
      IF (P.EQ.0) RETURN
      M = M/P
      R = M*S
      IF (P.GT.8) GO TO 700
      GO TO (100, 200, 300, 400, 500, 800, 700, 600), P
      GO TO 800
C
  200 CONTINUE
      CALL R2 CFTK (N, M, X(1), Y(1), X(R+1), Y(R+1), DIM)
      GO TO 100
C
  300 CONTINUE
      CALL R3 CFTK (N, M, X(1), Y(1), X(R+1), Y(R+1), X(2*R+1), Y(2*R+1)
     ., DIM)
      GO TO 100
C
  400 CONTINUE
      CALL R4 CFTK (N, M, X(1), Y(1), X(R+1), Y(R+1), X(2*R+1), Y(2*R+1)
     ., X(3*R+1), Y(3*R+1), DIM)
      GO TO 100
C
  500 CONTINUE
      CALL R5 CFTK (N, M, X(1), Y(1), X(R+1), Y(R+1), X(2*R+1), Y(2*R+1)
     ., X(3*R+1), Y(3*R+1), X(4*R+1), Y(4*R+1), DIM)
      GO TO 100
C
  600 CONTINUE
      CALL R8 CFTK (N, M, X(1), Y(1), X(R+1), Y(R+1), X(2*R+1), Y(2*R+1)
     ., X(3*R+1), Y(3*R+1), X(4*R+1), Y(4*R+1), X(5*R+1), Y(5*R+1),
     .X(6*R+1), Y(6*R+1), X(7*R+1), Y(7*R+1), DIM)
      GO TO 100
C
  700 CONTINUE
      CALL RP CFTK (N, M, P, R, X, Y, DIM)
      GO TO 100
C
  800 CONTINUE
      CALL MDP_ERROR(13,0)
      RETURN
C  900 FORMAT ('0TRANSFER ERROR DETECTED IN MDFTKD'//)
      END
      subroutine mdp_error
      entry mdp_error_t
      return
      end
C++*******************************************************************
C
C $$ R2CFTK:
C
C      CALL R2CFTK(N,M,X0,Y0,X1,Y1,DIM)
C        N
C        M
C        X0
C        Y0
C        X1
C        Y1
C        DIM
C
C--*******************************************************************
C
C
      SUBROUTINE R2 CFTK (N, M, X0, Y0, X1, Y1, DIM)                
C RADIX 2 MULTI-DIMENSIONAL COMPLEX FOURIER TRANSFORM KERNEL      
C                                                 
      INTEGER N, M, DIM(5)                                    
      REAL X0 (10), Y0 (10), X1 (10), Y1 (10)                     
C                                                 
      LOGICAL FOLD,ZERO                                     
      INTEGER J,K,K0,M2,M OVER 2                              
      INTEGER K1, K2, KK, L, L1, MM2, NT, SIZE, SEP                 
      INTEGER NS                                          
      REAL ANGLE,C,IS,IU,RS,RU,S,TWOPI                          
      REAL FJM1, FM2                                        
      DATA TWO PI/6.283185/                                   
C                                                 
      NT = DIM(1)                                           
      SEP = DIM(2)                                          
      L1 = DIM(3)                                           
      SIZE = DIM(4) - 1                                     
      K2 = DIM(5)                                           
      NS = N*SEP                                          
      M2=M*2                                              
      FM2 = FLOAT(M2)                                       
      M OVER 2=M/2+1                                        
      MM2 = SEP*M2                                          
C                                                 
      FJM1 = -1.0                                           
      DO 600 J=1,M OVER 2                                     
      FOLD=J.GT.1 .AND. 2*J.LT.M+2                              
      K0 = (J-1)*SEP + 1                                    
      FJM1 = FJM1 + 1.0                                     
      ANGLE = TWO PI*FJM1/FM2                                 
      ZERO=ANGLE.EQ.0.0                                     
      IF (ZERO) GO TO 200                                     
      C=COS(ANGLE)                                          
      S=SIN(ANGLE)                                          
      GO TO 200                                           
  100 CONTINUE                                            
      FOLD=.FALSE.                                          
      K0 = (M+1-J)*SEP + 1                                    
      C=-C                                                
  200 CONTINUE                                            
C                                                 
      DO 500 KK = K0, NS, MM2                                 
      DO 440 L = KK, NT, L1                                   
      K1 = L + SIZE                                         
      DO 420 K = L, K1, K2                                    
      RS=X0(K)+X1(K)                                        
      IS=Y0(K)+Y1(K)                                        
      RU=X0(K)-X1(K)                                        
      IU=Y0(K)-Y1(K)                                        
      X0(K)=RS                                            
      Y0(K)=IS                                            
      IF (ZERO) GO TO 300                                     
      X1(K)=RU*C+IU*S                                       
      Y1(K)=IU*C-RU*S                                       
      GO TO 400                                           
  300 CONTINUE                                            
      X1(K)=RU                                            
      Y1(K)=IU                                            
  400 CONTINUE                                            
  420 CONTINUE                                            
  440 CONTINUE                                            
  500 CONTINUE                                            
      IF (FOLD) GO TO 100                                     
  600 CONTINUE                                            
C                                                 
      RETURN                                              
      END                                                 
C++*******************************************************************
C
C $$ R3CFTK:
C
C      CALL R3CFTK(N,M,X0,Y0,X1,Y1,X2,Y2,DIM)
C        N
C        M
C        X0
C        Y0
C        X1
C        Y1
C        X2
C        Y2
C        DIM
C
C--*******************************************************************
C
C
      SUBROUTINE R3 CFTK (N, M, X0, Y0, X1, Y1, X2, Y2, DIM)          
C RADIX 3 MULTI-DIMENSIONAL COMPLEX FOURIER TRANSFORM KERNEL      
C                                                 
      INTEGER N, M, DIM(5)                                    
      REAL X0 (10), Y0 (10), X1 (10), Y1 (10), X2 (10), Y2 (10)       
C                                                 
      LOGICAL FOLD,ZERO                                     
      INTEGER J,K,K0,M3,M OVER 2                              
      INTEGER K1, K2, KK, L, L1, MM3, NT, SIZE, SEP                 
      INTEGER NS                                          
      REAL ANGLE,A,B,C1,C2,S1,S2,T,TWOPI                        
      REAL I0,I1,I2,IA,IB,IS,R0,R1,R2,RA,RB,RS                    
      REAL FJM1, FM3                                        
      DATA TWO PI/6.283185/, A/-0.5/, B/0.86602540/                 
C                                                 
      NT = DIM(1)                                           
      SEP = DIM(2)                                          
      L1 = DIM(3)                                           
      SIZE = DIM(4) - 1                                     
      K2 = DIM(5)                                           
      NS = N*SEP                                          
      M3=M*3                                              
      FM3 = FLOAT(M3)                                       
      MM3 = SEP*M3                                          
      M OVER 2=M/2+1                                        
C                                                 
      FJM1 = -1.0                                           
      DO 600 J=1,M OVER 2                                     
      FOLD=J.GT.1 .AND. 2*J.LT.M+2                              
      K0 = (J-1)*SEP + 1                                    
      FJM1 = FJM1 + 1.0                                     
      ANGLE = TWO PI*FJM1/FM3                                 
      ZERO=ANGLE.EQ.0.0                                     
      IF (ZERO) GO TO 200                                     
      C1=COS(ANGLE)                                         
      S1=SIN(ANGLE)                                         
      C2=C1*C1-S1*S1                                        
      S2=S1*C1+C1*S1                                        
      GO TO 200                                           
  100 CONTINUE                                            
      FOLD=.FALSE.                                          
      K0 = (M+1-J)*SEP + 1                                    
      T=C1*A+S1*B                                           
      S1=C1*B-S1*A                                          
      C1=T                                                
      T=C2*A-S2*B                                           
      S2=-C2*B-S2*A                                         
      C2=T                                                
  200 CONTINUE                                            
C                                                 
      DO 500 KK = K0, NS, MM3                                 
      DO 440 L = KK, NT, L1                                   
      K1 = L + SIZE                                         
      DO 420 K = L, K1, K2                                    
      R0=X0(K)                                            
      I0=Y0(K)                                            
      RS=X1(K)+X2(K)                                        
      IS=Y1(K)+Y2(K)                                        
      X0(K)=R0+RS                                           
      Y0(K)=I0+IS                                           
      RA=R0+RS*A                                          
      IA=I0+IS*A                                          
      RB=(X1(K)-X2(K))*B                                    
      IB=(Y1(K)-Y2(K))*B                                    
      IF (ZERO) GO TO 300                                     
      R1=RA+IB                                            
      I1=IA-RB                                            
      R2=RA-IB                                            
      I2=IA+RB                                            
      X1(K)=R1*C1+I1*S1                                     
      Y1(K)=I1*C1-R1*S1                                     
      X2(K)=R2*C2+I2*S2                                     
      Y2(K)=I2*C2-R2*S2                                     
      GO TO 400                                           
  300 CONTINUE                                            
      X1(K)=RA+IB                                           
      Y1(K)=IA-RB                                           
      X2(K)=RA-IB                                           
      Y2(K)=IA+RB                                           
  400 CONTINUE                                            
  420 CONTINUE                                            
  440 CONTINUE                                            
  500 CONTINUE                                            
      IF (FOLD) GO TO 100                                     
  600 CONTINUE                                            
C                                                 
      RETURN                                              
      END                                                 
C++*******************************************************************
C
C $$ R4CFTK:
C
C      CALL R4CFTK(N,M,X0,Y0,X1,Y1,X2,Y2,X3,Y3,DIM)
C        N
C        M
C        X0
C        Y0
C        X1
C        Y1
C        X2
C        Y2
C        X3
C        Y3
C        DIM
C
C--*******************************************************************
C
C
      SUBROUTINE R4 CFTK (N, M, X0, Y0, X1, Y1, X2, Y2, X3, Y3, DIM)    
C RADIX 4 MULTI-DIMENSIONAL COMPLEX FOURIER TRANSFORM KERNEL      
C                                                 
      INTEGER N, M, DIM(5)                                    
      REAL X0 (10), Y0 (10), X1 (10), Y1 (10)                     
      REAL X2 (10), Y2 (10), X3 (10), Y3 (10)                     
C                                                 
      LOGICAL FOLD,ZERO                                     
      INTEGER J,K,K0,M4,M OVER 2                              
      INTEGER K1, K2, KK, L, L1, MM4, NT, SIZE, SEP                 
      INTEGER NS                                          
      REAL ANGLE,C1,C2,C3,S1,S2,S3,T,TWOPI                        
      REAL I1,I2,I3,IS0,IS1,IU0,IU1,R1,R2,R3,RS0,RS1,RU0,RU1          
      REAL FJM1, FM4                                        
      DATA TWO PI/6.283185/                                   
C                                                 
      NT = DIM(1)                                           
      SEP = DIM(2)                                          
      L1 = DIM(3)                                           
      SIZE = DIM(4) - 1                                     
      K2 = DIM(5)                                           
      NS = N*SEP                                          
      M4=M*4                                              
      FM4 = FLOAT(M4)                                       
      MM4 = SEP*M4                                          
      M OVER 2=M/2+1                                        
C                                                 
      FJM1 = -1.0                                           
      DO 600 J=1,M OVER 2                                     
      FOLD=J.GT.1 .AND. 2*J.LT.M+2                              
      K0 = (J-1)*SEP + 1                                    
      FJM1 = FJM1 + 1.0                                     
      ANGLE = TWO PI*FJM1/FM4                                 
      ZERO=ANGLE.EQ.0.0                                     
      IF (ZERO) GO TO 200                                     
      C1=COS(ANGLE)                                         
      S1=SIN(ANGLE)                                         
      C2=C1*C1-S1*S1                                        
      S2=S1*C1+C1*S1                                        
      C3=C2*C1-S2*S1                                        
      S3=S2*C1+C2*S1                                        
      GO TO 200                                           
  100 CONTINUE                                            
      FOLD=.FALSE.                                          
      K0 = (M+1-J)*SEP + 1                                    
      T=C1                                                
      C1=S1                                               
      S1=T                                                
      C2=-C2                                              
      T=C3                                                
      C3=-S3                                              
      S3=-T                                               
  200 CONTINUE                                            
C                                                 
      DO 500 KK = K0, NS, MM4                                 
      DO 440 L = KK, NT, L1                                   
      K1 = L + SIZE                                         
      DO 420 K = L, K1, K2                                    
      RS0=X0(K)+X2(K)                                       
      IS0=Y0(K)+Y2(K)                                       
      RU0=X0(K)-X2(K)                                       
      IU0=Y0(K)-Y2(K)                                       
      RS1=X1(K)+X3(K)                                       
      IS1=Y1(K)+Y3(K)                                       
      RU1=X1(K)-X3(K)                                       
      IU1=Y1(K)-Y3(K)                                       
      X0(K)=RS0+RS1                                         
      Y0(K)=IS0+IS1                                         
      IF (ZERO) GO TO 300                                     
      R1=RU0+IU1                                          
      I1=IU0-RU1                                          
      R2=RS0-RS1                                          
      I2=IS0-IS1                                          
      R3=RU0-IU1                                          
      I3=IU0+RU1                                          
      X2(K)=R1*C1+I1*S1                                     
      Y2(K)=I1*C1-R1*S1                                     
      X1(K)=R2*C2+I2*S2                                     
      Y1(K)=I2*C2-R2*S2                                     
      X3(K)=R3*C3+I3*S3                                     
      Y3(K)=I3*C3-R3*S3                                     
      GO TO 400                                           
  300 CONTINUE                                            
      X2(K)=RU0+IU1                                         
      Y2(K)=IU0-RU1                                         
      X1(K)=RS0-RS1                                         
      Y1(K)=IS0-IS1                                         
      X3(K)=RU0-IU1                                         
      Y3(K)=IU0+RU1                                         
  400 CONTINUE                                            
  420 CONTINUE
  440 CONTINUE
  500 CONTINUE                                            
      IF (FOLD) GO TO 100                                     
  600 CONTINUE                                            
C                                                 
      RETURN                                              
      END                                                 
C++*******************************************************************
C
C $$ R5CFTK:
C
C      CALL R5CFTK(N,M,X0,Y0,X1,Y1,X2,Y2,X3,Y3,X4,Y4,DIM)
C        N
C        M
C        X0
C        Y0
C        X1
C        Y1
C        X2
C        Y2
C        X3
C        Y3
C        X4
C        Y4
C        DIM
C
C--*******************************************************************
C
C
      SUBROUTINE R5 CFTK (N, M, X0, Y0, X1, Y1, X2, Y2, X3, Y3, X4, Y4, 
     . DIM)                                               
C  RADIX 5 MULTI-DIMENSIONAL COMPLEX FOURIER TRANSFORM KERNEL      
C                                                  
      INTEGER N, M, DIM(5)                                    
      REAL X0 (10), Y0 (10), X1 (10), Y1 (10), X2 (10), Y2 (10)       
      REAL X3 (10), Y3 (10), X4 (10), Y4 (10)                     
C                                                  
      LOGICAL FOLD,ZERO                                     
      INTEGER J,K,K0,M5,M OVER 2                              
      INTEGER K1, K2, KK, L, L1, MM5, NT, SIZE, SEP                 
      INTEGER NS                                          
      REAL ANGLE,A1,A2,B1,B2,C1,C2,C3,C4,S1,S2,S3,S4,T,TWOPI          
      REAL R0,R1,R2,R3,R4,RA1,RA2,RB1,RB2,RS1,RS2,RU1,RU2             
      REAL I0,I1,I2,I3,I4,IA1,IA2,IB1,IB2,IS1,IS2,IU1,IU2             
      REAL FJM1, FM5                                        
      DATA TWO PI/6.283185/, A1/0.30901699/, B1/0.95105652/,          
     .      A2/-0.80901699/, B2/0.58778525/                       
C                                                  
      NT = DIM(1)                                           
      SEP = DIM(2)                                          
      L1 = DIM(3)                                           
      SIZE = DIM(4) - 1                                     
      K2 = DIM(5)                                           
      NS = N*SEP                                          
      M5=M*5                                              
      FM5 = FLOAT(M5)                                       
      MM5 = SEP*M5                                          
      M OVER 2=M/2+1                                        
C                                                  
      FJM1 = -1.0                                           
      DO 600 J=1,M OVER 2                                     
      FOLD=J.GT.1 .AND. 2*J.LT.M+2                              
      K0 = (J-1)*SEP + 1                                    
      FJM1 = FJM1 + 1.0                                     
      ANGLE = TWO PI*FJM1/FM5                                 
      ZERO=ANGLE.EQ.0.0                                     
      IF (ZERO) GO TO 200                                     
      C1=COS(ANGLE)                                         
      S1=SIN(ANGLE)                                         
      C2=C1*C1-S1*S1                                        
      S2=S1*C1+C1*S1                                        
      C3=C2*C1-S2*S1                                        
      S3=S2*C1+C2*S1                                        
      C4=C2*C2-S2*S2                                        
      S4=S2*C2+C2*S2                                        
      GO TO 200                                           
  100 CONTINUE                                            
      FOLD=.FALSE.                                          
      K0 = (M+1-J)*SEP + 1                                    
      T=C1*A1+S1*B1                                         
      S1=C1*B1-S1*A1                                        
      C1=T                                                
      T=C2*A2+S2*B2                                         
      S2=C2*B2-S2*A2                                        
      C2=T                                                
      T=C3*A2-S3*B2                                         
      S3=-C3*B2-S3*A2                                       
      C3=T                                                
      T=C4*A1-S4*B1                                         
      S4=-C4*B1-S4*A1                                       
      C4=T                                                
  200 CONTINUE                                            
C                                                  
      DO 500 KK = K0, NS, MM5                                 
      DO 440 L = KK, NT, L1                                   
      K1 = L + SIZE                                         
      DO 420 K = L, K1, K2                                    
      R0=X0(K)                                            
      I0=Y0(K)                                            
      RS1=X1(K)+X4(K)                                       
      IS1=Y1(K)+Y4(K)                                       
      RU1=X1(K)-X4(K)                                       
      IU1=Y1(K)-Y4(K)                                       
      RS2=X2(K)+X3(K)                                       
      IS2=Y2(K)+Y3(K)                                       
      RU2=X2(K)-X3(K)                                       
      IU2=Y2(K)-Y3(K)                                       
      X0(K)=R0+RS1+RS2                                      
      Y0(K)=I0+IS1+IS2                                      
      RA1=R0+RS1*A1+RS2*A2                                    
      IA1=I0+IS1*A1+IS2*A2                                    
      RA2=R0+RS1*A2+RS2*A1                                    
      IA2=I0+IS1*A2+IS2*A1                                    
      RB1=RU1*B1+RU2*B2                                     
      IB1=IU1*B1+IU2*B2                                     
      RB2=RU1*B2-RU2*B1                                     
      IB2=IU1*B2-IU2*B1                                     
      IF (ZERO) GO TO 300                                     
      R1=RA1+IB1                                          
      I1=IA1-RB1                                          
      R2=RA2+IB2                                          
      I2=IA2-RB2                                          
      R3=RA2-IB2                                          
      I3=IA2+RB2                                          
      R4=RA1-IB1                                          
      I4=IA1+RB1                                          
      X1(K)=R1*C1+I1*S1                                     
      Y1(K)=I1*C1-R1*S1                                     
      X2(K)=R2*C2+I2*S2                                     
      Y2(K)=I2*C2-R2*S2                                     
      X3(K)=R3*C3+I3*S3                                     
      Y3(K)=I3*C3-R3*S3                                     
      X4(K)=R4*C4+I4*S4                                     
      Y4(K)=I4*C4-R4*S4                                     
      GO TO 400                                           
  300 CONTINUE                                            
      X1(K)=RA1+IB1                                         
      Y1(K)=IA1-RB1                                         
      X2(K)=RA2+IB2                                         
      Y2(K)=IA2-RB2                                         
      X3(K)=RA2-IB2                                         
      Y3(K)=IA2+RB2                                         
      X4(K)=RA1-IB1                                         
      Y4(K)=IA1+RB1                                         
  400 CONTINUE                                            
  420 CONTINUE                                            
  440 CONTINUE                                            
  500 CONTINUE                                            
      IF (FOLD) GO TO 100                                     
  600 CONTINUE                                            
C                                                 
      RETURN                                              
      END                                                 
C++*******************************************************************
C
C $$ R8CFTK:
C
C      CALL R8CFTK(N,M,X0,Y0,X1,Y1,X2,Y2,X3,Y3,X4,Y4,X5,Y5,
C                     X6,Y6,X7,Y7,DIM)
C
C        N
C        M
C        X0
C        Y0
C        X1
C        Y1
C        X2
C        Y2
C        X3
C        Y3
C        X4
C        Y4
C        X5
C        Y5
C        X6
C        Y6
C        X7
C        Y7
C        DIM
C
C--*******************************************************************
C
C
      SUBROUTINE R8 CFTK (N, M, X0, Y0, X1, Y1, X2, Y2, X3, Y3, X4, Y4, 
     . X5, Y5, X6, Y6, X7, Y7, DIM)                             
C RADIX 8 MULTI-DIMENSIONAL COMPLEX FOURIER TRANSFORM KERNEL      
C                                                 
      INTEGER N, M, DIM(5)                                    
      REAL X0 (10), Y0 (10), X1 (10), Y1 (10), X2 (10), Y2 (10)       
      REAL X3 (10), Y3 (10), X4 (10), Y4 (10), X5 (10), Y5 (10)       
      REAL X6 (10), Y6 (10), X7 (10), Y7 (10)                     
C                                                 
      LOGICAL FOLD,ZERO                                     
      INTEGER J,K,K0,M8,M OVER 2                              
      INTEGER K1, K2, KK, L, L1, MM8, NT, SIZE, SEP                 
      INTEGER NS                                          
      REAL ANGLE,C1,C2,C3,C4,C5,C6,C7,E,S1,S2,S3,S4,S5,S6,S7,T,TWOPI    
      REAL R1,R2,R3,R4,R5,R6,R7,RS0,RS1,RS2,RS3,RU0,RU1,RU2,RU3       
      REAL I1,I2,I3,I4,I5,I6,I7,IS0,IS1,IS2,IS3,IU0,IU1,IU2,IU3       
      REAL RSS0,RSS1,RSU0,RSU1,RUS0,RUS1,RUU0,RUU1                  
      REAL ISS0,ISS1,ISU0,ISU1,IUS0,IUS1,IUU0,IUU1                  
      REAL FJM1, FM8                                        
      DATA TWO PI/6.283185/, E/0.70710678/                        
C                                                 
      NT = DIM(1)                                           
      SEP = DIM(2)                                          
      L1 = DIM(3)                                           
      SIZE = DIM(4) - 1                                     
      K2 = DIM(5)                                           
      NS = N*SEP                                          
      M8=M*8                                              
      FM8 = FLOAT(M8)                                       
      MM8 = SEP*M8                                          
      M OVER 2=M/2+1                                        
C                                                 
      FJM1 = -1.0                                           
      DO 600 J=1,M OVER 2                                     
      FOLD=J.GT.1 .AND. 2*J.LT.M+2                              
      K0 = (J-1)*SEP + 1                                    
      FJM1 = FJM1 + 1.0                                     
      ANGLE = TWO PI*FJM1/FM8                                 
      ZERO=ANGLE.EQ.0.0                                     
      IF (ZERO) GO TO 200                                     
      C1=COS(ANGLE)                                         
      S1=SIN(ANGLE)                                         
      C2=C1*C1-S1*S1                                        
      S2=S1*C1+C1*S1                                        
      C3=C2*C1-S2*S1                                        
      S3=S2*C1+C2*S1                                        
      C4=C2*C2-S2*S2                                        
      S4=S2*C2+C2*S2                                        
      C5=C4*C1-S4*S1                                        
      S5=S4*C1+C4*S1                                        
      C6=C4*C2-S4*S2                                        
      S6=S4*C2+C4*S2                                        
      C7=C4*C3-S4*S3                                        
      S7=S4*C3+C4*S3                                        
      GO TO 200                                           
  100 CONTINUE                                            
      FOLD=.FALSE.                                          
      K0 = (M+1-J)*SEP + 1                                    
      T=(C1+S1)*E                                           
      S1=(C1-S1)*E                                          
      C1=T                                                
      T=S2                                                
      S2=C2                                               
      C2=T                                                
      T=(-C3+S3)*E                                          
      S3=(C3+S3)*E                                          
      C3=T                                                
      C4=-C4                                              
      T=-(C5+S5)*E                                          
      S5=(-C5+S5)*E                                         
      C5=T                                                
      T=-S6                                               
      S6=-C6                                              
      C6=T                                                
      T=(C7-S7)*E                                           
      S7=-(C7+S7)*E                                         
      C7=T                                                
  200 CONTINUE                                            
C                                                 
      DO 500 KK = K0, NS, MM8                                 
      DO 440 L = KK, NT, L1                                   
      K1 = L + SIZE                                         
      DO 420 K = L, K1, K2                                    
      RS0=X0(K)+X4(K)                                       
      IS0=Y0(K)+Y4(K)                                       
      RU0=X0(K)-X4(K)                                       
      IU0=Y0(K)-Y4(K)                                       
      RS1=X1(K)+X5(K)                                       
      IS1=Y1(K)+Y5(K)                                       
      RU1=X1(K)-X5(K)                                       
      IU1=Y1(K)-Y5(K)                                       
      RS2=X2(K)+X6(K)                                       
      IS2=Y2(K)+Y6(K)                                       
      RU2=X2(K)-X6(K)                                       
      IU2=Y2(K)-Y6(K)                                       
      RS3=X3(K)+X7(K)                                       
      IS3=Y3(K)+Y7(K)                                       
      RU3=X3(K)-X7(K)                                       
      IU3=Y3(K)-Y7(K)                                       
      RSS0=RS0+RS2                                          
      ISS0=IS0+IS2                                          
      RSU0=RS0-RS2                                          
      ISU0=IS0-IS2                                          
      RSS1=RS1+RS3                                          
      ISS1=IS1+IS3                                          
      RSU1=RS1-RS3                                          
      ISU1=IS1-IS3                                          
      RUS0=RU0-IU2                                          
      IUS0=IU0+RU2                                          
      RUU0=RU0+IU2                                          
      IUU0=IU0-RU2                                          
      RUS1=RU1-IU3                                          
      IUS1=IU1+RU3                                          
      RUU1=RU1+IU3                                          
      IUU1=IU1-RU3                                          
      T=(RUS1+IUS1)*E                                       
      IUS1=(IUS1-RUS1)*E                                    
      RUS1=T                                              
      T=(RUU1+IUU1)*E                                       
      IUU1=(IUU1-RUU1)*E                                    
      RUU1=T                                              
      X0(K)=RSS0+RSS1                                       
      Y0(K)=ISS0+ISS1                                       
      IF (ZERO) GO TO 300                                     
      R1=RUU0+RUU1                                          
      I1=IUU0+IUU1                                          
      R2=RSU0+ISU1                                          
      I2=ISU0-RSU1                                          
      R3=RUS0+IUS1                                          
      I3=IUS0-RUS1                                          
      R4=RSS0-RSS1                                          
      I4=ISS0-ISS1                                          
      R5=RUU0-RUU1                                          
      I5=IUU0-IUU1                                          
      R6=RSU0-ISU1                                          
      I6=ISU0+RSU1                                          
      R7=RUS0-IUS1                                          
      I7=IUS0+RUS1                                          
      X4(K)=R1*C1+I1*S1                                     
      Y4(K)=I1*C1-R1*S1                                     
      X2(K)=R2*C2+I2*S2                                     
      Y2(K)=I2*C2-R2*S2                                     
      X6(K)=R3*C3+I3*S3                                     
      Y6(K)=I3*C3-R3*S3                                     
      X1(K)=R4*C4+I4*S4                                     
      Y1(K)=I4*C4-R4*S4                                     
      X5(K)=R5*C5+I5*S5                                     
      Y5(K)=I5*C5-R5*S5                                     
      X3(K)=R6*C6+I6*S6                                     
      Y3(K)=I6*C6-R6*S6                                     
      X7(K)=R7*C7+I7*S7                                     
      Y7(K)=I7*C7-R7*S7                                     
      GO TO 400                                           
  300 CONTINUE                                            
      X4(K)=RUU0+RUU1                                       
      Y4(K)=IUU0+IUU1                                       
      X2(K)=RSU0+ISU1                                       
      Y2(K)=ISU0-RSU1                                       
      X6(K)=RUS0+IUS1                                       
      Y6(K)=IUS0-RUS1                                       
      X1(K)=RSS0-RSS1                                       
      Y1(K)=ISS0-ISS1                                       
      X5(K)=RUU0-RUU1                                       
      Y5(K)=IUU0-IUU1                                       
      X3(K)=RSU0-ISU1                                       
      Y3(K)=ISU0+RSU1                                       
      X7(K)=RUS0-IUS1                                       
      Y7(K)=IUS0+RUS1                                       
  400 CONTINUE                                            
  420 CONTINUE                                            
  440 CONTINUE                                            
  500 CONTINUE                                            
      IF (FOLD) GO TO 100                                     
  600 CONTINUE                                            
C                                                 
      RETURN                                              
      END                                                 
      SUBROUTINE REAL2F(A,NX,NY,JD)
C ************************************** Copyright 1983 P.R. Smith.
      INTEGER ID(5)
      REAL A(1)
C
C  TWO DIMENSIONAL TRANSFORM OF REAL OBJECT
C
      IF(.NOT.(MOD(NX,2).EQ.0.OR.JD.LT.0))THEN
        CALL MDP_ERROR(13,0)
        RETURN
      END IF
C
C  SWITCH TO FORWARD OR REVERSE TRANSFORM
C
      IF(JD)1,2,2
C
C  REAL TRANSFORM
C
2     NXX=NX+2
      NS=NXX*NY
C
C  PACK OUT ROWS OF 'A' WITH ZEROS
C
      DO I=NY,1,-1
        L=I
        LF=NXX*(L-1)
        LN=NX*(L-1)
        DO J=NX,1,-1
          A(LF+J)=A(LN+J)
        END DO
      END DO
C
C  NOW PERFORM THE REAL F.T. AS STEP ONE
      NT=NX/2
      ID(1)=NS
      ID(2)=2
      ID(3)=NS
      ID(4)=NS
      ID(5)=NXX
C
      CALL REALFT(A(1),A(2),NT,ID)
C
C  NOW THE COMPLEX F.T.
C
      ID(1)=NS
      ID(2)=NXX
      ID(3)=NS
      ID(4)=NXX
      ID(5)=2
C
      CALL CMPLFT(A(1),A(2),NY,ID)
C
      NX=NT+1
      RETURN
C
C  HERMITIAN SYMMETRIC TRANSFORM
C
1     CONTINUE
      NXX=2*NX
      NT=NXX-2
       NS=NY*NXX
C
C  COMPLEX TRANSFORM
C
      ID(1)=NS
      ID(2)=NXX
      ID(3)=NS
      ID(4)=NXX
      ID(5)=2
C
C  CONJUGATE ARRAY AND TRANSFORM
C
      DO10I=2,NS,2
10    A(I)=-A(I)
C
      CALL CMPLFT(A(1),A(2),NY,ID)
C
C  COPY LAST ELEMENT OVER AND COMPRESS
C
      DO12I=1,NY
      LY=NXX*(I-1)
      LM=NT*(I-1)
      A(2+LY)=A(NXX-1+LY)
      DO13J=1,NT
13    A(LM+J)=A(LY+J)
12    CONTINUE
C
C  NOW THE HERMITIAN F.T.
C
      NS=NY*NT
      ID(1)=NS
      ID(2)=2
      ID(3)=NS
      ID(4)=NS
      ID(5)=NT
C
      CALL HERMFT(A(1),A(2),NT/2,ID)
C
      NX=NT
      RETURN
      END
      SUBROUTINE REAL3F(A,NX,NY,NZ,JD)
C ************************************** Copyright 1990 P.R. Smith.
      INTEGER ID(5)
      REAL A(*)
C
C  THREE DIMENSIONAL TRANSFORM OF REAL OBJECT
C
      IF(.NOT.(MOD(NX,2).EQ.0.OR.JD.LT.0))THEN
        CALL MDP_ERROR(13,0)
        RETURN
      END IF
C
C  SWITCH TO FORWARD OR REVERSE TRANSFORM
C
      IF(JD)1,2,2
C
C  REAL TRANSFORM
C
2     NXX=NX+2
      NS=NXX*NY*NZ
      NXY=NXX*NY
      NYNZ=NY*NZ
C
C  PACK OUT ROWS OF 'A' WITH ZEROS
C
      DO I=NYNZ,1,-1
        L=I
        LF=NXX*(L-1)
        LN=NX*(L-1)
        DO J=NX,1,-1
          A(LF+J)=A(LN+J)
        END DO
      END DO
C
C  NOW PERFORM THE REAL F.T. AS STEP ONE
C
      NT=NX/2
      ID(1)=NS
      ID(2)=2
      ID(3)=NS
      ID(4)=NS
      ID(5)=NXX
C
      CALL REALFT(A(1),A(2),NT,ID)
C
C  NOW THE COMPLEX F.T. IN "Y"
C
      ID(1)=NXY
      ID(2)=NXX
      ID(3)=NXY
      ID(4)=NXX
      ID(5)=2
C
      DO I=1,NZ
        NZ1=NXY*(I-1)+1
        CALL CMPLFT(A(NZ1),A(NZ1+1),NY,ID)
      END DO
C
C  NOW THE COMPLEX F.T. IN "Z"
C
      ID(1)=NS
      ID(2)=NXY
      ID(3)=NS
      ID(4)=NXY
      ID(5)=2
C
      CALL CMPLFT(A(1),A(2),NZ,ID)
C
      NX=NT+1
      RETURN
C
C  HERMITIAN SYMMETRIC TRANSFORM
C
1     CONTINUE
      NXX=2*NX
      NT=NXX-2
      NS=NY*NXX*NZ
      NXY=NXX*NY
C
C  COMPLEX TRANSFORM "Z"
C
      ID(1)=NS
      ID(2)=NXY
      ID(3)=NS
      ID(4)=NXY
      ID(5)=2
C
C  CONJUGATE ARRAY AND TRANSFORM
C
      DO I=2,NS,2
        A(I)=-A(I)
      END DO
C
      CALL CMPLFT(A(1),A(2),NZ,ID)
C
C  COMPLEX TRANSFORM "Y"
C
      ID(1)=NXY
      ID(2)=NXX
      ID(3)=NXY
      ID(4)=NXX
      ID(5)=2
C
      DO I=1,NZ
        NZ1=NXY*(I-1)+1
        CALL CMPLFT(A(NZ1),A(NZ1+1),NY,ID)
      END DO
C
C  COPY LAST ELEMENT OVER AND COMPRESS
C
      DO I=1,NY*NZ
        LY=NXX*(I-1)
        LM=NT*(I-1)
        A(2+LY)=A(NXX-1+LY)
        DO J=1,NT
          A(LM+J)=A(LY+J)
        END DO
      END DO
C
C  NOW THE HERMITIAN F.T.
C
      NS=NY*NT*NZ
      ID(1)=NS
      ID(2)=2
      ID(3)=NS
      ID(4)=NS
      ID(5)=NT
C
      CALL HERMFT(A(1),A(2),NT/2,ID)
C
      NX=NT
      RETURN
      END
      SUBROUTINE REALFT (EVEN, ODD, N, DIM)
      INTEGER N
      INTEGER DIM(5)
      REAL EVEN(1), ODD(1)
C
      REAL A, B, C, D, E, F, ANGLE, CO, SI, TWO PI, TWO N
      INTEGER NT, D2, D3, D4, D5
      INTEGER I, J, K, L, N OVER 2, I0, I1, I2
C
      TWO PI = 6.283185
      TWO N = FLOAT(2*N)
C
      CALL CMPL FT (EVEN, ODD, N, DIM)
C
      NT = DIM(1)
      D2 = DIM(2)
      D3 = DIM(3)
      D4 = DIM(4) - 1
      D5 = DIM(5)
      N OVER 2 = N/2 + 1
C
      IF (N OVER 2 .LT. 2) GO TO 400
      DO 300 I = 2, N OVER 2
      ANGLE = TWO PI*FLOAT(I-1)/TWO N
      CO = COS(ANGLE)
      SI = SIN(ANGLE)
      I0 = (I - 1)*D2 + 1
      J = (N + 2 - 2*I)*D2
      DO 200 I1 = I0, NT, D3
      I2 = I1 + D4
      DO 100 K = I1, I2, D5
      L = K + J
      A = (EVEN(L) + EVEN(K))/2.0
      C = (EVEN(L) - EVEN(K))/2.0
      B = (ODD(L) + ODD(K))/2.0
      D = (ODD(L) - ODD(K))/2.0
      E = C*SI + B*CO
      F = C*CO - B*SI
      EVEN(K) = A + E
      EVEN(L) = A - E
      ODD(K) = F - D
      ODD(L) = F + D
  100 CONTINUE
  200 CONTINUE
  300 CONTINUE
C
  400 CONTINUE
      IF (N .LT. 1) GO TO 600
      J = N*D2
      DO 500 I1 = 1, NT, D3
      I2 = I1 + D4
      DO 500 K = I1, I2, D5
      L = K + J
      EVEN(L) = EVEN(K) - ODD(K)
      ODD(L) = 0.0
      EVEN(K) = EVEN(K) + ODD(K)
      ODD(K) = 0.0
  500 CONTINUE
C
  600 CONTINUE
      RETURN
      END
      SUBROUTINE RPCFTK (N, M, P, R, X, Y, DIM)
C     RADIX PRIME MULTI-DIMENSIONAL COMPLEX FOURIER TRANSFORM KERNEL
C
      INTEGER N, M, P, R, DIM(5)
      REAL X(R,P), Y(R,P)
C
      LOGICAL FOLD,ZERO
      REAL ANGLE,IS,IU,RS,RU,T,TWOPI,XT,YT
      REAL FU, FP, FJM1, FMP
      INTEGER J,JJ,K0,K,M OVER 2,MP,PM,PP,U,V
      INTEGER K1, K2, KK, L, L1, MMP, NT, SIZE, SEP
      INTEGER NS
C
      REAL AA (20,20), BB (20,20)
      REAL A (40), B (40), C (40), S (40)
      REAL IA (20), IB (20), RA (20), RB (20)
      DATA TWO PI/6.283185/
C
      NT = DIM(1)
      SEP = DIM(2)
      L1 = DIM(3)
      SIZE = DIM(4) - 1
      K2 = DIM(5)
      NS = N*SEP
      M OVER 2=M/2+1
      MP=M*P
      FMP = FLOAT(MP)
      MMP = SEP*MP
      PP=P/2
      PM=P-1
      FP = FLOAT(P)
      FU = 0.0
      DO 100 U=1,PP
      FU = FU + 1.0
      ANGLE = TWO PI*FU/FP
      JJ=P-U
      A(U)=COS(ANGLE)
      B(U)=SIN(ANGLE)
      A(JJ)=A(U)
      B(JJ)=-B(U)
  100 CONTINUE
      DO 300 U=1,PP
      DO 200 V=1,PP
      JJ=U*V-U*V/P*P
      AA(V,U)=A(JJ)
      BB(V,U)=B(JJ)
  200 CONTINUE
  300 CONTINUE
C
      FJM1 = -1.0
      DO 1500 J=1,M OVER 2
      FOLD=J.GT.1 .AND. 2*J.LT.M+2
      K0 = (J-1)*SEP + 1
      FJM1 = FJM1 + 1.0
      ANGLE = TWO PI*FJM1/FMP
      ZERO=ANGLE.EQ.0.0
      IF (ZERO) GO TO 700
      C(1)=COS(ANGLE)
      S(1)=SIN(ANGLE)
      DO 400 U=2,PM
      C(U)=C(U-1)*C(1)-S(U-1)*S(1)
      S(U)=S(U-1)*C(1)+C(U-1)*S(1)
  400 CONTINUE
      GO TO 700
  500 CONTINUE
      FOLD=.FALSE.
      K0 = (M+1-J)*SEP + 1
      DO 600 U=1,PM
      T=C(U)*A(U)+S(U)*B(U)
      S(U)=-S(U)*A(U)+C(U)*B(U)
      C(U)=T
  600 CONTINUE
  700 CONTINUE
C
      DO 1400 KK = K0, NS, MMP
      DO 1340 L = KK, NT, L1
      K1 = L + SIZE
      DO 1320 K = L, K1, K2
      XT=X(K,1)
      YT=Y(K,1)
      RS=X(K,2)+X(K,P)
      IS=Y(K,2)+Y(K,P)
      RU=X(K,2)-X(K,P)
      IU=Y(K,2)-Y(K,P)
      DO 800 U=1,PP
      RA(U)=XT+RS*AA(U,1)
      IA(U)=YT+IS*AA(U,1)
      RB(U)=RU*BB(U,1)
      IB(U)=IU*BB(U,1)
  800 CONTINUE
      XT=XT+RS
      YT=YT+IS
      DO 1000 U=2,PP
      JJ=P-U
      RS=X(K,U+1)+X(K,JJ+1)
      IS=Y(K,U+1)+Y(K,JJ+1)
      RU=X(K,U+1)-X(K,JJ+1)
      IU=Y(K,U+1)-Y(K,JJ+1)
      XT=XT+RS
      YT=YT+IS
      DO 900 V=1,PP
      RA(V)=RA(V)+RS*AA(V,U)
      IA(V)=IA(V)+IS*AA(V,U)
      RB(V)=RB(V)+RU*BB(V,U)
      IB(V)=IB(V)+IU*BB(V,U)
  900 CONTINUE
 1000 CONTINUE
      X(K,1)=XT
      Y(K,1)=YT
      DO 1300 U=1,PP
      JJ=P-U
      IF (ZERO) GO TO 1100
      XT=RA(U)+IB(U)
      YT=IA(U)-RB(U)
      X(K,U+1)=XT*C(U)+YT*S(U)
      Y(K,U+1)=YT*C(U)-XT*S(U)
      XT=RA(U)-IB(U)
      YT=IA(U)+RB(U)
      X(K,JJ+1)=XT*C(JJ)+YT*S(JJ)
      Y(K,JJ+1)=YT*C(JJ)-XT*S(JJ)
      GO TO 1200
 1100 CONTINUE
      X(K,U+1)=RA(U)+IB(U)
      Y(K,U+1)=IA(U)-RB(U)
      X(K,JJ+1)=RA(U)-IB(U)
      Y(K,JJ+1)=IA(U)+RB(U)
 1200 CONTINUE
 1300 CONTINUE
 1320 CONTINUE
 1340 CONTINUE
 1400 CONTINUE
      IF (FOLD) GO TO 500
 1500 CONTINUE
C
      RETURN
      END
      SUBROUTINE SRFP (PTS,PMAX,TWO GRP,FACTOR,SYM,P SYM,UN SYM,ERROR)
C     SYMMETRIZED REORDERING FACTORING PROGRAMME
C
      LOGICAL ERROR
      INTEGER PTS,PMAX,TWO GRP,P SYM
      INTEGER FACTOR (10), SYM (10), UN SYM (10)
C
      INTEGER PP(14), QQ (7)
      INTEGER F,J,JJ,N,NEST,P,P TWO,Q,R
C
      NEST=14
C
      N=PTS
      P SYM=1
      F=2
      P=0
      Q=0
  100 CONTINUE
      IF (N.LE.1) GO TO 500
      DO 200 J=F,PMAX
      IF (N.EQ.(N/J)*J) GO TO 300
  200 CONTINUE
      GO TO 1400
  300 CONTINUE
      IF (2*P+Q.GE.NEST) GO TO 1600
      F=J
      N=N/F
      IF (N.EQ.(N/F)*F) GO TO 400
      Q=Q+1
      QQ(Q)=F
      GO TO 100
  400 CONTINUE
      N=N/F
      P=P+1
      PP(P)=F
      P SYM=P SYM*F
      GO TO 100
C
  500 CONTINUE
      R=1
      IF (Q.EQ.0) R=0
      IF (P.LT.1) GO TO 700
      DO 600 J=1,P
      JJ=P+1-J
      SYM(J)=PP(JJ)
      FACTOR(J)=PP(JJ)
      JJ=P+Q+J
      FACTOR(JJ)=PP(J)
      JJ=P+R+J
      SYM(JJ)=PP(J)
  600 CONTINUE
  700 CONTINUE
      IF (Q.LT.1) GO TO 900
      DO 800 J=1,Q
      JJ=P+J
      UN SYM(J)=QQ(J)
      FACTOR(JJ)=QQ(J)
  800 CONTINUE
      SYM(P+1)=PTS/P SYM**2
  900 CONTINUE
      JJ=2*P+Q
      FACTOR(JJ+1)=0
      P TWO=1
      J=0
 1000 CONTINUE
      J=J+1
      IF (FACTOR(J).EQ.0) GO TO 1200
      IF (FACTOR(J).NE.2) GO TO 1000
      P TWO=P TWO*2
      FACTOR(J)=1
      IF (P TWO.GE.TWO GRP) GO TO 1100
      IF (FACTOR(J+1).EQ.2) GO TO 1000
 1100 CONTINUE
      FACTOR(J)=P TWO
      P TWO=1
      GO TO 1000
 1200 CONTINUE
      IF (P.EQ.0) R=0
      JJ=2*P+R
      SYM(JJ+1)=0
      IF (Q.LE.1) Q=0
      UN SYM(Q+1)=0
      ERROR=.FALSE.
C
 1300 CONTINUE
      RETURN
C
 1400 CONTINUE
      call mdp_error_t(56,'Largest factor exceeds maximum')
      ERROR=.TRUE.
      GO TO 1300
C 1500 FORMAT (24H0LARGEST FACTOR EXCEEDS ,I3,7H.  N = ,I6,1H.)
C
 1600 CONTINUE
      call mdp_error_t(56,'Factor count exceeded')
      ERROR=.TRUE.
      GO TO 1300
C 1700 FORMAT (22H0FACTOR COUNT EXCEEDS ,I3,7H.  N = ,I6,1H.)
      END
