C*********************************************************************C
C**   Mixed-Radix Complex Fast Fourier Transform                    **C
C**    for Massively Parallel Processor HITACHI SR2201              **C
C**                                                                 **C
C**   All rights reserved, Copyright (C) 1998, Daisuke TAKAHASHI    **C
C**                                                                 **C
C*********************************************************************C
c/*-
c * Copyright (c) 1998 Information Promotion Agency of Japan and D.Takahashi
c * All rights reserved.
c *
c * Redistribution and non-commercial use in source and binary forms,
c * with or without modification, are permitted provided that the following
c * conditions are met:
c *
c * 1. Redistributions of source code must retain the above copyright
c *    notice, this list of conditions and the following disclaimer.
c *
c * 2. Redistributions in binary form must reproduce the above copyright
c *    notice, this list of conditions and the following disclaimer in the
c *    documentation and/or other materials provided with the distribution.
c *
c * 3. All advertising materials mentioning features or use of this
c *    software must display the following acknowledgement:
c *
c *    This product includes software developed by Information Promotion
c *    Agency of Japan and D.Takahashi
c *
c * 4. The name "Information Promotion Agency" or "D.Takahashi" should not be
c *    used to endorse or promote products derived from this software
c *    without specific prior written permission.
c *
c * 5. Any use of the source code or the binary in a commercial product,
c *    whether may it be the origial representation or in some modified form,
c *    is not permitted without specific prior written permission.
c *
c * THIS SOFTWARE IS PROVIDED BY THE AUTHOR AND CONTRIBUTORS ``AS IS'' AND
c * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
c * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
c * ARE DISCLAIMED.  IN NO EVENT SHALL THE AUTHOR OR CONTRIBUTORS BE LIABLE
c * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
c * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
c * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
c * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
c * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
c * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
c * SUCH DAMAGE.
c */
      SUBROUTINE FFT(A,B,WW,WWW,N,IP,ME,NPU,IOPT,IERR)
      IMPLICIT REAL*8 (A-H,O-Z)
      PARAMETER (NDA2=4096)
      COMPLEX*16 A(*),B(*),WW(*),WWW(*)
      COMPLEX*16 W1(NDA2/2+1),W2(NDA2/2),W3(NDA2/2+1)
      DIMENSION IP(3),IP1(3),IP2(3),IP3(3)
      DATA N0/0/
C
      IF (MOD(N,NPU*NPU) .NE. 0) THEN
        IERR=2000
        RETURN
      END IF
      IF (MOD(N,2) .NE. 0 .AND. MOD(N,3) .NE. 0 .AND.
     1    MOD(N,5) .NE. 0) THEN
        IERR=2000
        RETURN
      END IF
C
      NN=N/NPU
      CALL FACTOR(NPU,IP1)
      IP2(1)=(IP(1)-IP1(1))/2
      IP3(1)=IP(1)-IP1(1)-IP2(1)
      IP2(2)=(IP(2)-IP1(2))/2
      IP3(2)=IP(2)-IP1(2)-IP2(2)
      IP2(3)=IP(3)-(IP(3)-IP1(3))/2
      IP3(3)=IP(3)-IP1(3)-IP2(3)
      DO 5 I=1,3
        IF (IP3(I) .LT. IP1(I)) THEN
          IP3(I)=IP1(I)
          IP2(I)=IP(I)-IP1(I)-IP3(I)
        END IF
    5 CONTINUE
C
      NX=(2**IP1(1))*(3**IP1(2))*(5**IP1(3))
      NY=(2**IP2(1))*(3**IP2(2))*(5**IP2(3))
      NZ=(2**IP3(1))*(3**IP3(2))*(5**IP3(3))
C
      IF (N .NE. NX*NY*NZ) THEN
        IERR=2000
        RETURN
      END IF
C
      IF (N .NE. N0) THEN
        CALL SETTBL(W1,NX)
        CALL SETTBL(W2,NY)
        CALL SETTBL(W3,NZ)
        CALL SETTBL2(WW,NY,NZ)
        CALL SETTBLP(WWW,NZ*NY,ME,NPU)
        N0=N
      END IF
C
      KEY=1
      CALL FFT2345(A,B,W3,NY,NZ,IP3,IOPT,KEY)
      IF (KEY .GE. 0) THEN
        CALL CTRANST(A,B,WW,NY,NZ,IOPT)
      ELSE
        CALL CTRANST(B,A,WW,NY,NZ,IOPT)
      END IF
      KEY=-KEY
      CALL FFT2345(A,B,W2,NZ,NY,IP2,IOPT,KEY)
      IF (KEY .GE. 0) THEN
        CALL CTRANST(A,B,WWW,NPU,NZ/NPU*NY,IOPT)
        CALL CTRANSP(B,A,NN,NPU)
      ELSE
        CALL CTRANST(B,A,WWW,NPU,NZ/NPU*NY,IOPT)
        CALL CTRANSP(A,B,NN,NPU)
      END IF
      CALL FFT2345(A,B,W1,NZ/NPU*NY,NX,IP1,IOPT,KEY)
      IF (IOPT .EQ. 1) THEN
        IF (KEY .LT. 0) THEN
          DO 10 I=1,NN
            A(I)=B(I)
   10     CONTINUE
        END IF
      ELSE
        DN=1.0D0/DBLE(N)
        IF (KEY .GE. 0) THEN
          DO 20 I=1,NN
            A(I)=A(I)*DN
   20     CONTINUE
        ELSE
          DO 30 I=1,NN
            A(I)=B(I)*DN
   30     CONTINUE
        END IF
      END IF
      IERR=0
      RETURN
      END
      SUBROUTINE FFT2345(A,B,W,NM,N,IP,IOPT,KEY)
      IMPLICIT REAL*8 (A-H,O-Z)
      COMPLEX*16 A(*),B(*),W(*)
      DIMENSION IP(3)
C
      L=N
      M=1
      DO 10 K=1,IP(1)/2
        L=L/4
        IF (KEY .GE. 0) THEN
          CALL FFT4(A,B,W,NM,M,L,IOPT)
        ELSE
          CALL FFT4(B,A,W,NM,M,L,IOPT)
        END IF
        KEY=-KEY
        M=M*4
   10 CONTINUE
      DO 20 K=1,IP(2)
        L=L/3
        IF (KEY .GE. 0) THEN
          CALL FFT3(A,B,W,NM,M,L,IOPT)
        ELSE
          CALL FFT3(B,A,W,NM,M,L,IOPT)
        END IF
        KEY=-KEY
        M=M*3
   20 CONTINUE
      DO 30 K=1,IP(3)
        L=L/5
        IF (KEY .GE. 0) THEN
          CALL FFT5(A,B,W,NM,M,L,IOPT)
        ELSE
          CALL FFT5(B,A,W,NM,M,L,IOPT)
        END IF
        KEY=-KEY
        M=M*5
   30 CONTINUE
      IF (MOD(IP(1),2) .EQ. 1) THEN
        IF (KEY .GE. 0) THEN
          CALL FFT2(A,B,NM,M)
        ELSE
          CALL FFT2(B,A,NM,M)
        END IF
        KEY=-KEY
      END IF
      RETURN
      END
      SUBROUTINE FFT2(A,B,NM,M)
      IMPLICIT REAL*8 (A-H,O-Z)
      DIMENSION A(2,NM*M,*),B(2,NM*M,*)
C
      DO 10 I=1,NM*M
        B(1,I,1)=A(1,I,1)+A(1,I,2)
        B(2,I,1)=A(2,I,1)+A(2,I,2)
        B(1,I,2)=A(1,I,1)-A(1,I,2)
        B(2,I,2)=A(2,I,1)-A(2,I,2)
   10 CONTINUE
      RETURN
      END
      SUBROUTINE FFT3(A,B,W,NM,M,L,IOPT)
      IMPLICIT REAL*8 (A-H,O-Z)
      DIMENSION A(2,NM*M,L,*),B(2,NM*M,3,*),W(2,M,L,*)
C
      PI=4.0D0*DATAN(1.0D0)
      S60=DSIN(PI/3.0D0)
      IF (IOPT .EQ. 1) THEN
        DO 10 I=1,NM*M
          X0=A(1,I,1,2)+A(1,I,1,3)
          Y0=A(2,I,1,2)+A(2,I,1,3)
          X1=A(1,I,1,1)-0.5D0*X0
          Y1=A(2,I,1,1)-0.5D0*Y0
          X2=S60*(A(2,I,1,2)-A(2,I,1,3))
          Y2=S60*(A(1,I,1,3)-A(1,I,1,2))
          B(1,I,1,1)=A(1,I,1,1)+X0
          B(2,I,1,1)=A(2,I,1,1)+Y0
          B(1,I,2,1)=X1+X2
          B(2,I,2,1)=Y1+Y2
          B(1,I,3,1)=X1-X2
          B(2,I,3,1)=Y1-Y2
   10   CONTINUE
        DO 30 J=2,L
          WR1=W(1,1,J,1)
          WI1=W(2,1,J,1)
          WR2=WR1*WR1-WI1*WI1
          WI2=WR1*WI1+WI1*WR1
          DO 20 I=1,NM*M
            X0=A(1,I,J,2)+A(1,I,J,3)
            Y0=A(2,I,J,2)+A(2,I,J,3)
            X1=A(1,I,J,1)-0.5D0*X0
            Y1=A(2,I,J,1)-0.5D0*Y0
            X2=S60*(A(2,I,J,2)-A(2,I,J,3))
            Y2=S60*(A(1,I,J,3)-A(1,I,J,2))
            B(1,I,1,J)=A(1,I,J,1)+X0
            B(2,I,1,J)=A(2,I,J,1)+Y0
            B(1,I,2,J)=WR1*(X1+X2)-WI1*(Y1+Y2)
            B(2,I,2,J)=WR1*(Y1+Y2)+WI1*(X1+X2)
            B(1,I,3,J)=WR2*(X1-X2)-WI2*(Y1-Y2)
            B(2,I,3,J)=WR2*(Y1-Y2)+WI2*(X1-X2)
   20     CONTINUE
   30   CONTINUE
      ELSE
        DO 40 I=1,NM*M
          X0=A(1,I,1,2)+A(1,I,1,3)
          Y0=A(2,I,1,2)+A(2,I,1,3)
          X1=A(1,I,1,1)-0.5D0*X0
          Y1=A(2,I,1,1)-0.5D0*Y0
          X2=S60*(A(2,I,1,3)-A(2,I,1,2))
          Y2=S60*(A(1,I,1,2)-A(1,I,1,3))
          B(1,I,1,1)=A(1,I,1,1)+X0
          B(2,I,1,1)=A(2,I,1,1)+Y0
          B(1,I,2,1)=X1+X2
          B(2,I,2,1)=Y1+Y2
          B(1,I,3,1)=X1-X2
          B(2,I,3,1)=Y1-Y2
   40   CONTINUE
        DO 60 J=2,L
          WR1=W(1,1,J,1)
          WI1=W(2,1,J,1)
          WR2=WR1*WR1-WI1*WI1
          WI2=WR1*WI1+WI1*WR1
          DO 50 I=1,NM*M
            X0=A(1,I,J,2)+A(1,I,J,3)
            Y0=A(2,I,J,2)+A(2,I,J,3)
            X1=A(1,I,J,1)-0.5D0*X0
            Y1=A(2,I,J,1)-0.5D0*Y0
            X2=S60*(A(2,I,J,3)-A(2,I,J,2))
            Y2=S60*(A(1,I,J,2)-A(1,I,J,3))
            B(1,I,1,J)=A(1,I,J,1)+X0
            B(2,I,1,J)=A(2,I,J,1)+Y0
            B(1,I,2,J)=WR1*(X1+X2)+WI1*(Y1+Y2)
            B(2,I,2,J)=WR1*(Y1+Y2)-WI1*(X1+X2)
            B(1,I,3,J)=WR2*(X1-X2)+WI2*(Y1-Y2)
            B(2,I,3,J)=WR2*(Y1-Y2)-WI2*(X1-X2)
   50     CONTINUE
   60   CONTINUE
      END IF
      RETURN
      END
      SUBROUTINE FFT4(A,B,W,NM,M,L,IOPT)
      IMPLICIT REAL*8 (A-H,O-Z)
      DIMENSION A(2,NM*M,L,*),B(2,NM*M,4,*),W(2,M,L,*)
C
      IF (IOPT .EQ. 1) THEN
        DO 10 I=1,NM*M
          X0=A(1,I,1,1)+A(1,I,1,3)
          Y0=A(2,I,1,1)+A(2,I,1,3)
          X1=A(1,I,1,1)-A(1,I,1,3)
          Y1=A(2,I,1,1)-A(2,I,1,3)
          X2=A(1,I,1,2)+A(1,I,1,4)
          Y2=A(2,I,1,2)+A(2,I,1,4)
          X3=A(2,I,1,2)-A(2,I,1,4)
          Y3=A(1,I,1,4)-A(1,I,1,2)
          B(1,I,1,1)=X0+X2
          B(2,I,1,1)=Y0+Y2
          B(1,I,2,1)=X1+X3
          B(2,I,2,1)=Y1+Y3
          B(1,I,3,1)=X0-X2
          B(2,I,3,1)=Y0-Y2
          B(1,I,4,1)=X1-X3
          B(2,I,4,1)=Y1-Y3
   10   CONTINUE
        DO 30 J=2,L
          WR1=W(1,1,J,1)
          WI1=W(2,1,J,1)
          WR2=WR1*WR1-WI1*WI1
          WI2=WR1*WI1+WI1*WR1
          WR3=WR1*WR2-WI1*WI2
          WI3=WR1*WI2+WI1*WR2
          DO 20 I=1,NM*M
            X0=A(1,I,J,1)+A(1,I,J,3)
            Y0=A(2,I,J,1)+A(2,I,J,3)
            X1=A(1,I,J,1)-A(1,I,J,3)
            Y1=A(2,I,J,1)-A(2,I,J,3)
            X2=A(1,I,J,2)+A(1,I,J,4)
            Y2=A(2,I,J,2)+A(2,I,J,4)
            X3=A(2,I,J,2)-A(2,I,J,4)
            Y3=A(1,I,J,4)-A(1,I,J,2)
            B(1,I,1,J)=X0+X2
            B(2,I,1,J)=Y0+Y2
            B(1,I,2,J)=WR1*(X1+X3)-WI1*(Y1+Y3)
            B(2,I,2,J)=WR1*(Y1+Y3)+WI1*(X1+X3)
            B(1,I,3,J)=WR2*(X0-X2)-WI2*(Y0-Y2)
            B(2,I,3,J)=WR2*(Y0-Y2)+WI2*(X0-X2)
            B(1,I,4,J)=WR3*(X1-X3)-WI3*(Y1-Y3)
            B(2,I,4,J)=WR3*(Y1-Y3)+WI3*(X1-X3)
   20     CONTINUE
   30   CONTINUE
      ELSE
        DO 40 I=1,NM*M
          X0=A(1,I,1,1)+A(1,I,1,3)
          Y0=A(2,I,1,1)+A(2,I,1,3)
          X1=A(1,I,1,1)-A(1,I,1,3)
          Y1=A(2,I,1,1)-A(2,I,1,3)
          X2=A(1,I,1,2)+A(1,I,1,4)
          Y2=A(2,I,1,2)+A(2,I,1,4)
          X3=A(2,I,1,4)-A(2,I,1,2)
          Y3=A(1,I,1,2)-A(1,I,1,4)
          B(1,I,1,1)=X0+X2
          B(2,I,1,1)=Y0+Y2
          B(1,I,2,1)=X1+X3
          B(2,I,2,1)=Y1+Y3
          B(1,I,3,1)=X0-X2
          B(2,I,3,1)=Y0-Y2
          B(1,I,4,1)=X1-X3
          B(2,I,4,1)=Y1-Y3
   40   CONTINUE
        DO 60 J=2,L
          WR1=W(1,1,J,1)
          WI1=W(2,1,J,1)
          WR2=WR1*WR1-WI1*WI1
          WI2=WR1*WI1+WI1*WR1
          WR3=WR1*WR2-WI1*WI2
          WI3=WR1*WI2+WI1*WR2
          DO 50 I=1,NM*M
            X0=A(1,I,J,1)+A(1,I,J,3)
            Y0=A(2,I,J,1)+A(2,I,J,3)
            X1=A(1,I,J,1)-A(1,I,J,3)
            Y1=A(2,I,J,1)-A(2,I,J,3)
            X2=A(1,I,J,2)+A(1,I,J,4)
            Y2=A(2,I,J,2)+A(2,I,J,4)
            X3=A(2,I,J,4)-A(2,I,J,2)
            Y3=A(1,I,J,2)-A(1,I,J,4)
            B(1,I,1,J)=X0+X2
            B(2,I,1,J)=Y0+Y2
            B(1,I,2,J)=WR1*(X1+X3)+WI1*(Y1+Y3)
            B(2,I,2,J)=WR1*(Y1+Y3)-WI1*(X1+X3)
            B(1,I,3,J)=WR2*(X0-X2)+WI2*(Y0-Y2)
            B(2,I,3,J)=WR2*(Y0-Y2)-WI2*(X0-X2)
            B(1,I,4,J)=WR3*(X1-X3)+WI3*(Y1-Y3)
            B(2,I,4,J)=WR3*(Y1-Y3)-WI3*(X1-X3)
   50     CONTINUE
   60   CONTINUE
      END IF
      RETURN
      END
      SUBROUTINE FFT5(A,B,W,NM,M,L,IOPT)
      IMPLICIT REAL*8 (A-H,O-Z)
      DIMENSION A(2,NM*M,L,*),B(2,NM*M,5,*),W(2,M,L,*)
C
      PI=4.0D0*DATAN(1.0D0)
      S72=DSIN(2.0D0*PI/5.0D0)
      S3672=DSIN(PI/5.0D0)/S72
      SQR5Q=DSQRT(5.0D0)*0.25D0
      IF (IOPT .EQ. 1) THEN
        DO 10 I=1,NM*M
          X0=A(1,I,1,2)+A(1,I,1,5)
          Y0=A(2,I,1,2)+A(2,I,1,5)
          X1=A(1,I,1,3)+A(1,I,1,4)
          Y1=A(2,I,1,3)+A(2,I,1,4)
          X2=S72*(A(1,I,1,2)-A(1,I,1,5))
          Y2=S72*(A(2,I,1,2)-A(2,I,1,5))
          X3=S72*(A(1,I,1,3)-A(1,I,1,4))
          Y3=S72*(A(2,I,1,3)-A(2,I,1,4))
          X4=X0+X1
          Y4=Y0+Y1
          X5=SQR5Q*(X0-X1)
          Y5=SQR5Q*(Y0-Y1)
          X6=A(1,I,1,1)-0.25D0*X4
          Y6=A(2,I,1,1)-0.25D0*Y4
          X7=X6+X5
          Y7=Y6+Y5
          X8=X6-X5
          Y8=Y6-Y5
          X9=Y2+S3672*Y3
          Y9=-X2-S3672*X3
          X10=S3672*Y2-Y3
          Y10=X3-S3672*X2
          B(1,I,1,1)=A(1,I,1,1)+X4
          B(2,I,1,1)=A(2,I,1,1)+Y4
          B(1,I,2,1)=X7+X9
          B(2,I,2,1)=Y7+Y9
          B(1,I,3,1)=X8+X10
          B(2,I,3,1)=Y8+Y10
          B(1,I,4,1)=X8-X10
          B(2,I,4,1)=Y8-Y10
          B(1,I,5,1)=X7-X9
          B(2,I,5,1)=Y7-Y9
   10   CONTINUE
        DO 30 J=2,L
          WR1=W(1,1,J,1)
          WI1=W(2,1,J,1)
          WR2=WR1*WR1-WI1*WI1
          WI2=WR1*WI1+WI1*WR1
          WR3=WR1*WR2-WI1*WI2
          WI3=WR1*WI2+WI1*WR2
          WR4=WR2*WR2-WI2*WI2
          WI4=WR2*WI2+WI2*WR2
          DO 20 I=1,NM*M
            X0=A(1,I,J,2)+A(1,I,J,5)
            Y0=A(2,I,J,2)+A(2,I,J,5)
            X1=A(1,I,J,3)+A(1,I,J,4)
            Y1=A(2,I,J,3)+A(2,I,J,4)
            X2=S72*(A(1,I,J,2)-A(1,I,J,5))
            Y2=S72*(A(2,I,J,2)-A(2,I,J,5))
            X3=S72*(A(1,I,J,3)-A(1,I,J,4))
            Y3=S72*(A(2,I,J,3)-A(2,I,J,4))
            X4=X0+X1
            Y4=Y0+Y1
            X5=SQR5Q*(X0-X1)
            Y5=SQR5Q*(Y0-Y1)
            X6=A(1,I,J,1)-0.25D0*X4
            Y6=A(2,I,J,1)-0.25D0*Y4
            X7=X6+X5
            Y7=Y6+Y5
            X8=X6-X5
            Y8=Y6-Y5
            X9=Y2+S3672*Y3
            Y9=-X2-S3672*X3
            X10=S3672*Y2-Y3
            Y10=X3-S3672*X2
            B(1,I,1,J)=A(1,I,J,1)+X4
            B(2,I,1,J)=A(2,I,J,1)+Y4
            B(1,I,2,J)=WR1*(X7+X9)-WI1*(Y7+Y9)
            B(2,I,2,J)=WR1*(Y7+Y9)+WI1*(X7+X9)
            B(1,I,3,J)=WR2*(X8+X10)-WI2*(Y8+Y10)
            B(2,I,3,J)=WR2*(Y8+Y10)+WI2*(X8+X10)
            B(1,I,4,J)=WR3*(X8-X10)-WI3*(Y8-Y10)
            B(2,I,4,J)=WR3*(Y8-Y10)+WI3*(X8-X10)
            B(1,I,5,J)=WR4*(X7-X9)-WI4*(Y7-Y9)
            B(2,I,5,J)=WR4*(Y7-Y9)+WI4*(X7-X9)
   20     CONTINUE
   30   CONTINUE
      ELSE
        DO 40 I=1,NM*M
          X0=A(1,I,1,2)+A(1,I,1,5)
          Y0=A(2,I,1,2)+A(2,I,1,5)
          X1=A(1,I,1,3)+A(1,I,1,4)
          Y1=A(2,I,1,3)+A(2,I,1,4)
          X2=S72*(A(1,I,1,2)-A(1,I,1,5))
          Y2=S72*(A(2,I,1,2)-A(2,I,1,5))
          X3=S72*(A(1,I,1,3)-A(1,I,1,4))
          Y3=S72*(A(2,I,1,3)-A(2,I,1,4))
          X4=X0+X1
          Y4=Y0+Y1
          X5=SQR5Q*(X0-X1)
          Y5=SQR5Q*(Y0-Y1)
          X6=A(1,I,1,1)-0.25D0*X4
          Y6=A(2,I,1,1)-0.25D0*Y4
          X7=X6+X5
          Y7=Y6+Y5
          X8=X6-X5
          Y8=Y6-Y5
          X9=-Y2-S3672*Y3
          Y9=X2+S3672*X3
          X10=Y3-S3672*Y2
          Y10=S3672*X2-X3
          B(1,I,1,1)=A(1,I,1,1)+X4
          B(2,I,1,1)=A(2,I,1,1)+Y4
          B(1,I,2,1)=X7+X9
          B(2,I,2,1)=Y7+Y9
          B(1,I,3,1)=X8+X10
          B(2,I,3,1)=Y8+Y10
          B(1,I,4,1)=X8-X10
          B(2,I,4,1)=Y8-Y10
          B(1,I,5,1)=X7-X9
          B(2,I,5,1)=Y7-Y9
   40   CONTINUE
        DO 60 J=2,L
          WR1=W(1,1,J,1)
          WI1=W(2,1,J,1)
          WR2=WR1*WR1-WI1*WI1
          WI2=WR1*WI1+WI1*WR1
          WR3=WR1*WR2-WI1*WI2
          WI3=WR1*WI2+WI1*WR2
          WR4=WR2*WR2-WI2*WI2
          WI4=WR2*WI2+WI2*WR2
          DO 50 I=1,NM*M
            X0=A(1,I,J,2)+A(1,I,J,5)
            Y0=A(2,I,J,2)+A(2,I,J,5)
            X1=A(1,I,J,3)+A(1,I,J,4)
            Y1=A(2,I,J,3)+A(2,I,J,4)
            X2=S72*(A(1,I,J,2)-A(1,I,J,5))
            Y2=S72*(A(2,I,J,2)-A(2,I,J,5))
            X3=S72*(A(1,I,J,3)-A(1,I,J,4))
            Y3=S72*(A(2,I,J,3)-A(2,I,J,4))
            X4=X0+X1
            Y4=Y0+Y1
            X5=SQR5Q*(X0-X1)
            Y5=SQR5Q*(Y0-Y1)
            X6=A(1,I,J,1)-0.25D0*X4
            Y6=A(2,I,J,1)-0.25D0*Y4
            X7=X6+X5
            Y7=Y6+Y5
            X8=X6-X5
            Y8=Y6-Y5
            X9=-Y2-S3672*Y3
            Y9=X2+S3672*X3
            X10=Y3-S3672*Y2
            Y10=S3672*X2-X3
            B(1,I,1,J)=A(1,I,J,1)+X4
            B(2,I,1,J)=A(2,I,J,1)+Y4
            B(1,I,2,J)=WR1*(X7+X9)+WI1*(Y7+Y9)
            B(2,I,2,J)=WR1*(Y7+Y9)-WI1*(X7+X9)
            B(1,I,3,J)=WR2*(X8+X10)+WI2*(Y8+Y10)
            B(2,I,3,J)=WR2*(Y8+Y10)-WI2*(X8+X10)
            B(1,I,4,J)=WR3*(X8-X10)+WI3*(Y8-Y10)
            B(2,I,4,J)=WR3*(Y8-Y10)-WI3*(X8-X10)
            B(1,I,5,J)=WR4*(X7-X9)+WI4*(Y7-Y9)
            B(2,I,5,J)=WR4*(Y7-Y9)-WI4*(X7-X9)
   50     CONTINUE
   60   CONTINUE
      END IF
      RETURN
      END
      SUBROUTINE SETTBL(W,N)
      IMPLICIT REAL*8 (A-H,O-Z)
      DIMENSION W(2,*)
C
      PI2=8.0D0*DATAN(1.0D0)
      PX=-PI2/DBLE(N)
      DO 10 I=1,N/2
        W(1,I)=DCOS(PX*DBLE(I-1))
        W(2,I)=DSIN(PX*DBLE(I-1))
   10 CONTINUE
      RETURN
      END
      SUBROUTINE SETTBL2(W,N1,N2)
      IMPLICIT REAL*8 (A-H,O-Z)
      DIMENSION W(2,N1,*)
C
      N=N1*N2
      PI2=8.0D0*DATAN(1.0D0)
      PX=-PI2/DBLE(N)
      DO 20 K=1,N2
        DO 10 J=1,N1
          W(1,J,K)=DCOS(PX*DBLE(J-1)*DBLE(K-1))
          W(2,J,K)=DSIN(PX*DBLE(J-1)*DBLE(K-1))
   10   CONTINUE
   20 CONTINUE
      RETURN
      END
      SUBROUTINE SETTBLP(W,N1,ME,NPU)
      IMPLICIT REAL*8 (A-H,O-Z)
      DIMENSION W(2,*)
C
      N=N1*NPU
      PI2=8.0D0*DATAN(1.0D0)
      PX=-PI2/DBLE(N)
      DO 10 I=1,N1
        W(1,I)=DCOS(PX*DBLE(I-1)*DBLE(ME))
        W(2,I)=DSIN(PX*DBLE(I-1)*DBLE(ME))
   10 CONTINUE
      RETURN
      END
      SUBROUTINE CTRANST(A,B,W,N1,N2,IOPT)
      IMPLICIT REAL*8 (A-H,O-Z)
      COMPLEX*16 A(N1,*),B(N2,*),W(N1,*)
C
      IF (IOPT .EQ. 1) THEN
        DO 20 J=1,N1
          DO 10 I=1,N2
            B(I,J)=A(J,I)*W(J,I)
   10     CONTINUE
   20   CONTINUE
      ELSE
        DO 40 J=1,N1
          DO 30 I=1,N2
            B(I,J)=A(J,I)*DCONJG(W(J,I))
   30     CONTINUE
   40   CONTINUE
      END IF
      RETURN
      END
      SUBROUTINE CTRANSP(A,B,NN,NPU)
      IMPLICIT REAL*8 (A-H,O-Z)
      INCLUDE 'mpif.h'
      COMPLEX*16 A(*),B(*)
C
      NN2=NN/NPU
      CALL MPI_ALLTOALL(A,NN2,MPI_DOUBLE_COMPLEX,
     1                  B,NN2,MPI_DOUBLE_COMPLEX,
     2                  MPI_COMM_WORLD,IERR)
      RETURN
      END
      SUBROUTINE FACTOR(N,IP)
      DIMENSION IP(3)
C
      N2=N
      IF (MOD(N2,2) .NE. 0 .AND. MOD(N2,3) .NE. 0 .AND.
     1    MOD(N2,5) .NE. 0) RETURN
      IP(1)=0
      IP(2)=0
      IP(3)=0
   10 IF (N2 .LE. 1) RETURN
      IF (MOD(N2,2) .EQ. 0) THEN
        IP(1)=IP(1)+1
        N2=N2/2
        GO TO 10
      ELSE IF (MOD(N2,3) .EQ. 0) THEN
        IP(2)=IP(2)+1
        N2=N2/3
        GO TO 10
      ELSE IF (MOD(N2,5) .EQ. 0) THEN
        IP(3)=IP(3)+1
        N2=N2/5
        GO TO 10
      END IF
      RETURN
      END