************************************************************************* * ISPACK FORTRAN SUBROUTINE LIBRARY FOR SCIENTIFIC COMPUTING * * Copyright (C) 1998 Keiichi Ishioka * * * * This library is free software; you can redistribute it and/or * * modify it under the terms of the GNU Library General Public * * License as published by the Free Software Foundation; either * * version 2 of the License, or (at your option) any later version. * * * * This library is distributed in the hope that it will be useful, * * but WITHOUT ANY WARRANTY; without even the implied warranty of * * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU * * Library General Public License for more details. * * * * You should have received a copy of the GNU Library General Public * * License along with this library; if not, write to the Free * * Software Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA. * ************************************************************************* * COSINE TRANSFORM (MID-POINT) 2000/09/19 K.Ishioka ************************************************************************* ************************************************************************ SUBROUTINE FTTCMI(N,IT,T) IMPLICIT REAL*8(A-H,O-Z) PARAMETER(PI=3.1415926535897932385D0) DIMENSION T(0:N/2-1,12),IT(5) CALL FTTRUI(N,IT,T) N2=N*2 DO I=0,N/2-1 T(I,5)=SIN(2*PI*(2*I+1)/N2) T(I,6)=COS(2*PI*(2*I+1)/N2) T(I,7)=COS(PI*I/N) T(I,8)=SIN(PI*I/N) T(I,9)=SIN(PI*(4*I+1)/N2) T(I,10)=SIN(PI*(4*I+3)/N2) T(I,11)=1/T(I,9) T(I,12)=1/T(I,10) END DO END ************************************************************************ SUBROUTINE FTTCMF(M,N,X,Y,IT,T) IMPLICIT REAL*8(A-H,O-Z) DIMENSION IT(5),T(0:N/2-1,12) DIMENSION X(M,0:N-1),Y(M,0:N-1) DO J=0,N/2-1 J1=2*J J2=N-1-J1 J3=2*J+1 J4=N-1-J3 DO I=1,M Y(I,J1)=X(I,J1)+X(I,J2)+2*T(J,9 )*(X(I,J1)-X(I,J2)) Y(I,J3)=X(I,J3)+X(I,J4)+2*T(J,10)*(X(I,J3)-X(I,J4)) END DO END DO CALL FTTRUF(M,N,Y,X,IT,T) DO I=1,M X(I,0)=Y(I,0) X(I,N-1)=0.5D0*Y(I,1) END DO DO J=N/2-1,1,-1 DO I=1,M X(I,2*J)=T(J,7)*Y(I,2*J)+T(J,8)*Y(I,2*J+1) X(I,2*J-1)=X(I,2*J+1)+T(J,8)*Y(I,2*J)-T(J,7)*Y(I,2*J+1) END DO END DO END ************************************************************************ SUBROUTINE FTTCMB(M,N,X,Y,IT,T) IMPLICIT REAL*8(A-H,O-Z) DIMENSION IT(5),T(0:N/2-1,12) DIMENSION X(M,0:N-1),Y(M,0:N-1) DO I=1,M Y(I,0)=X(I,0) Y(I,1)=2*X(I,N-1) END DO DO J=1,N/2-1 DO I=1,M Y(I,2*J)=T(J,7)*X(I,2*J)-(X(I,2*J+1)-X(I,2*J-1))*T(J,8) Y(I,2*J+1)=T(J,8)*X(I,2*J)+(X(I,2*J+1)-X(I,2*J-1))*T(J,7) END DO END DO CALL FTTRUB(M,N,Y,X,IT,T) DO J=0,N/2-1 DO I=1,M X(I,2*J )=0.25D0*(Y(I,2*J)+Y(I,2*(N/2-1-J)+1)) & +0.125D0*(Y(I,2*J)-Y(I,2*(N/2-1-J)+1))*T(J,11) X(I,2*J+1)=0.25D0*(Y(I,2*J+1)+Y(I,2*(N/2-1-J))) & +0.125D0*(Y(I,2*J+1)-Y(I,2*(N/2-1-J)))*T(J,12) END DO END DO END ************************************************************************* * ISPACK FORTRAN SUBROUTINE LIBRARY FOR SCIENTIFIC COMPUTING * * Copyright (C) 1998 Keiichi Ishioka * * * * This library is free software; you can redistribute it and/or * * modify it under the terms of the GNU Library General Public * * License as published by the Free Software Foundation; either * * version 2 of the License, or (at your option) any later version. * * * * This library is distributed in the hope that it will be useful, * * but WITHOUT ANY WARRANTY; without even the implied warranty of * * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU * * Library General Public License for more details. * * * * You should have received a copy of the GNU Library General Public * * License along with this library; if not, write to the Free * * Software Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA. * ************************************************************************* * COSINE TRANSFORM (TRAPEZOIDAL) 2000/09/19 K.Ishioka ************************************************************************* SUBROUTINE FTTCTI(N,IT,T) IMPLICIT REAL*8(A-H,O-Z) PARAMETER(PI=3.1415926535897932385D0) DIMENSION T(0:N/2-1,6),IT(5) CALL FTTRUI(N,IT,T) N2=N*2 DO I=0,N/2-1 T(I,5)=SIN(2*PI*(2*I+1)/N2) T(I,6)=COS(2*PI*(2*I+1)/N2) END DO END ************************************************************************ SUBROUTINE FTTCTF(M,N,X,Y,IT,T) IMPLICIT REAL*8(A-H,O-Z) DIMENSION IT(5),T(0:N/2-1,6) DIMENSION X(M,0:N),Y(M,0:N-1) DO J=0,N/2-1 J1=2*J J2=N-J1 J3=2*J+1 J4=N-J3 DO I=1,M Y(I,J1)=X(I,J1)+X(I,J2)+2*T(J,4)*(X(I,J1)-X(I,J2)) Y(I,J3)=X(I,J3)+X(I,J4)+2*T(J,5)*(X(I,J3)-X(I,J4)) END DO END DO F=1D0/N DO I=1,M X(I,N)=(X(I,0)-X(I,N)+2*X(I,1)*T(0,6))*F END DO F=2D0/N DO J=1,N/2-1 DO I=1,M X(I,N)=X(I,N)+(X(I,2*J)*T(J,3)+X(I,2*J+1)*T(J,6))*F END DO END DO CALL FTTRUF(M,N,Y,X,IT,T) DO I=1,M X(I,0)=Y(I,0) X(I,1)=X(I,N) X(I,N)=Y(I,1) END DO DO J=1,N/2-1 DO I=1,M X(I,2*J)=Y(I,2*J) X(I,2*J+1)=X(I,2*J-1)+Y(I,2*J+1) END DO END DO END ************************************************************************ SUBROUTINE FTTCTB(M,N,X,Y,IT,T) IMPLICIT REAL*8(A-H,O-Z) DIMENSION IT(5),T(0:N/2-1,6) DIMENSION X(M,0:N),Y(M,0:N-1) CALL FTTCTF(M,N,X,Y,IT,T) DO J=0,N DO I=1,M X(I,J)=(N/2D0)*X(I,J) END DO END DO END ************************************************************************* * ISPACK FORTRAN SUBROUTINE LIBRARY FOR SCIENTIFIC COMPUTING * * Copyright (C) 1998 Keiichi Ishioka * * * * This library is free software; you can redistribute it and/or * * modify it under the terms of the GNU Library General Public * * License as published by the Free Software Foundation; either * * version 2 of the License, or (at your option) any later version. * * * * This library is distributed in the hope that it will be useful, * * but WITHOUT ANY WARRANTY; without even the implied warranty of * * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU * * Library General Public License for more details. * * * * You should have received a copy of the GNU Library General Public * * License along with this library; if not, write to the Free * * Software Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA. * ************************************************************************* ************************************************************************ * DUMP MESSAGES *----------------------------------------------------------------------- SUBROUTINE FTDMSG(CL,CS,CM) CHARACTER CL*1,CS*(*),CM*(*) CHARACTER CSD*6,CMD*53 DATA MMSG,IMSG/20,0/ SAVE CSD=CS CMD=CM IF(CL.EQ.'E') THEN WRITE(6,'(A)') '***** ERROR ('//CSD//') *** '//CMD STOP END IF IF(IMSG.LT.MMSG) THEN IF(CL.EQ.'W') THEN IMSG=IMSG+1 WRITE(*,*) '*** WARNING ('//CSD//') *** '//CMD ELSE IF(CL.EQ.'M') THEN IMSG=IMSG+1 WRITE(*,*) '*** MESSAGE ('//CSD//') *** '//CMD END IF IF(IMSG.EQ.MMSG) THEN WRITE(*,*) '+++ THE FOLLOWING MESSAGES ARE SUPRRESSED.' END IF END IF END ************************************************************************* * ISPACK FORTRAN SUBROUTINE LIBRARY FOR SCIENTIFIC COMPUTING * * Copyright (C) 1998 Keiichi Ishioka * * * * This library is free software; you can redistribute it and/or * * modify it under the terms of the GNU Library General Public * * License as published by the Free Software Foundation; either * * version 2 of the License, or (at your option) any later version. * * * * This library is distributed in the hope that it will be useful, * * but WITHOUT ANY WARRANTY; without even the implied warranty of * * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU * * Library General Public License for more details. * * * * You should have received a copy of the GNU Library General Public * * License along with this library; if not, write to the Free * * Software Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA. * ************************************************************************* ************************************************************************ SUBROUTINE FTTRUI(N,IT,T) IMPLICIT REAL*8(A-H,O-Z) PARAMETER(PI=3.1415926535897932385D0) DIMENSION T(0:N/2-1,4),IT(5) IF(MOD(N,2).NE.0) THEN CALL FTDMSG('E','FTTRUI','N MUST BE EVEN.') END IF L=N/2 CALL FTTZLI(L,IT,T) DO 10 I=0,L-1 T(I,3)=COS(2*PI*I/N) T(I,4)=SIN(2*PI*I/N) 10 CONTINUE END ************************************************************************ SUBROUTINE FTTRUF(M,N,X,Y,IT,T) IMPLICIT REAL*8(A-H,O-Z) DIMENSION X(M,2,0:N/2-1),Y(M,0:N/2-1,2) DIMENSION T(0:N/2-1,4),IT(5) L=N/2 DO 20 J=0,L-1 DO 10 I=1,M Y(I,J,1)=X(I,1,J) Y(I,J,2)=X(I,2,J) 10 CONTINUE 20 CONTINUE CALL FTTZLM(M,L,Y,X,IT,T) R=1D0/N S=R/2 DO 30 I=1,M X(I,1,0)=R*(Y(I,0,1)+Y(I,0,2)) X(I,2,0)=R*(Y(I,0,1)-Y(I,0,2)) 30 CONTINUE DO 50 J=1,L-1 DO 40 I=1,M X(I,1,J)=S*( (Y(I,L-J,1)+Y(I,J,1)) & +T(J,3)*(Y(I,L-J,2)+Y(I,J,2)) & -T(J,4)*(Y(I,L-J,1)-Y(I,J,1))) X(I,2,J)=S*( (Y(I,L-J,2)-Y(I,J,2)) & -T(J,3)*(Y(I,L-J,1)-Y(I,J,1)) & -T(J,4)*(Y(I,L-J,2)+Y(I,J,2))) 40 CONTINUE 50 CONTINUE END ************************************************************************ SUBROUTINE FTTRUB(M,N,X,Y,IT,T) IMPLICIT REAL*8(A-H,O-Z) DIMENSION X(M,2,0:N/2-1),Y(M,0:N/2-1,2) DIMENSION T(0:N/2-1,4),IT(5) L=N/2 DO 10 I=1,M Y(I,0,1)=X(I,1,0)+X(I,2,0) Y(I,0,2)=X(I,1,0)-X(I,2,0) 10 CONTINUE DO 30 J=1,L-1 DO 20 I=1,M Y(I,J,1)= (X(I,1,L-J)+X(I,1,J)) & -T(J,3)*(X(I,2,L-J)+X(I,2,J)) & +T(J,4)*(X(I,1,L-J)-X(I,1,J)) Y(I,J,2)= -(X(I,2,L-J)-X(I,2,J)) & -T(J,3)*(X(I,1,L-J)-X(I,1,J)) & -T(J,4)*(X(I,2,L-J)+X(I,2,J)) 20 CONTINUE 30 CONTINUE CALL FTTZLM(M,L,Y,X,IT,T) DO 50 J=0,L-1 DO 40 I=1,M X(I,1,J)=Y(I,J,1) X(I,2,J)=Y(I,J,2) 40 CONTINUE 50 CONTINUE END ************************************************************************* * ISPACK FORTRAN SUBROUTINE LIBRARY FOR SCIENTIFIC COMPUTING * * Copyright (C) 1998 Keiichi Ishioka * * * * This library is free software; you can redistribute it and/or * * modify it under the terms of the GNU Library General Public * * License as published by the Free Software Foundation; either * * version 2 of the License, or (at your option) any later version. * * * * This library is distributed in the hope that it will be useful, * * but WITHOUT ANY WARRANTY; without even the implied warranty of * * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU * * Library General Public License for more details. * * * * You should have received a copy of the GNU Library General Public * * License along with this library; if not, write to the Free * * Software Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA. * ************************************************************************* * SINE TRANSFORM (MID-POINT) 2000/09/19 K.Ishioka ************************************************************************* ************************************************************************ SUBROUTINE FTTSMI(N,IT,T) IMPLICIT REAL*8(A-H,O-Z) PARAMETER(PI=3.1415926535897932385D0) DIMENSION T(0:N/2-1,12),IT(5) CALL FTTRUI(N,IT,T) N2=N*2 DO I=0,N/2-1 T(I,5)=SIN(2*PI*(2*I+1)/N2) T(I,6)=COS(2*PI*(2*I+1)/N2) T(I,7)=COS(PI*I/N) T(I,8)=SIN(PI*I/N) T(I,9)=SIN(PI*(4*I+1)/N2) T(I,10)=SIN(PI*(4*I+3)/N2) T(I,11)=1/T(I,9) T(I,12)=1/T(I,10) END DO END ************************************************************************ SUBROUTINE FTTSMF(M,N,X,Y,IT,T) IMPLICIT REAL*8(A-H,O-Z) DIMENSION IT(5),T(0:N/2-1,12) DIMENSION X(M,0:N-1),Y(M,0:N-1) DO J=0,N/2-1 J1=2*J J2=N-1-J1 J3=2*J+1 J4=N-1-J3 DO I=1,M Y(I,2*J)=-(X(I,J1)-X(I,J2))+2*T(J,9 )*(X(I,J1)+X(I,J2)) Y(I,2*J+1)=-(X(I,J3)-X(I,J4))+2*T(J,10)*(X(I,J3)+X(I,J4)) END DO END DO CALL FTTRUF(M,N,Y,X,IT,T) DO I=1,M X(I,0)=0.5D0*Y(I,0) X(I,N-1)=-Y(I,1) END DO DO J=1,N/2-1 DO I=1,M X(I,2*J-1)=T(J,7)*Y(I,2*J+1)-T(J,8)*Y(I,2*J) X(I,2*J)=X(I,2*J-2)+T(J,7)*Y(I,2*J)+T(J,8)*Y(I,2*J+1) END DO END DO END ************************************************************************ SUBROUTINE FTTSMB(M,N,X,Y,IT,T) IMPLICIT REAL*8(A-H,O-Z) DIMENSION IT(5),T(0:N/2-1,12) DIMENSION X(M,0:N-1),Y(M,0:N-1) DO I=1,M Y(I,0)=2*X(I,0) Y(I,1)=-X(I,N-1) END DO DO J=1,N/2-1 DO I=1,M Y(I,2*J)=(X(I,2*J)-X(I,2*J-2))*T(J,7)-T(J,8)*X(I,2*J-1) Y(I,2*J+1)=(X(I,2*J)-X(I,2*J-2))*T(J,8)+T(J,7)*X(I,2*J-1) END DO END DO CALL FTTRUB(M,N,Y,X,IT,T) DO J=0,N/2-1 DO I=1,M X(I,2*J )=-0.25D0*(Y(I,2*J)-Y(I,2*(N/2-1-J)+1)) & +0.125D0*(Y(I,2*J)+Y(I,2*(N/2-1-J)+1))*T(J,11) X(I,2*J+1)=-0.25D0*(Y(I,2*J+1)-Y(I,2*(N/2-1-J))) & +0.125D0*(Y(I,2*J+1)+Y(I,2*(N/2-1-J)))*T(J,12) END DO END DO END ************************************************************************* * ISPACK FORTRAN SUBROUTINE LIBRARY FOR SCIENTIFIC COMPUTING * * Copyright (C) 1998 Keiichi Ishioka * * * * This library is free software; you can redistribute it and/or * * modify it under the terms of the GNU Library General Public * * License as published by the Free Software Foundation; either * * version 2 of the License, or (at your option) any later version. * * * * This library is distributed in the hope that it will be useful, * * but WITHOUT ANY WARRANTY; without even the implied warranty of * * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU * * Library General Public License for more details. * * * * You should have received a copy of the GNU Library General Public * * License along with this library; if not, write to the Free * * Software Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA. * ************************************************************************* * SINE TRANSFORM (TRAPEZOIDAL) 2000/09/19 K.Ishioka ************************************************************************* ************************************************************************ SUBROUTINE FTTSTI(N,IT,T) IMPLICIT REAL*8(A-H,O-Z) PARAMETER(PI=3.1415926535897932385D0) DIMENSION T(0:N/2-1,5),IT(5) CALL FTTRUI(N,IT,T) N2=N*2 DO I=0,N/2-1 T(I,5)=SIN(2*PI*(2*I+1)/N2) END DO END ************************************************************************ SUBROUTINE FTTSTF(M,N,X,Y,IT,T) IMPLICIT REAL*8(A-H,O-Z) DIMENSION IT(5),T(0:N/2-1,5) DIMENSION X(M,N),Y(M,0:N-1) DO J=1,N/2-1 J1=2*J J2=N-J1 J3=2*J+1 J4=N-J3 DO I=1,M Y(I,J1)=X(I,J2)-X(I,J1)+2*T(J,4)*(X(I,J2)+X(I,J1)) Y(I,J3)=X(I,J4)-X(I,J3)+2*T(J,5)*(X(I,J4)+X(I,J3)) END DO END DO J=0 J3=2*J+1 J4=N-J3 DO I=1,M Y(I,0)=0 Y(I,1)=X(I,J4)-X(I,J3)+2*T(J,5)*(X(I,J4)+X(I,J3)) END DO CALL FTTRUF(M,N,Y,X,IT,T) DO I=1,M X(I,N)=0 X(I,1)=0.5D0*Y(I,0) * X(I,N-1)=-0.5D0*Y(I,1) END DO DO J=1,N/2-1 DO I=1,M X(I,2*J)=Y(I,2*J+1) X(I,2*J+1)=X(I,2*J-1)+Y(I,2*J) END DO END DO END ************************************************************************ SUBROUTINE FTTSTB(M,N,X,Y,IT,T) IMPLICIT REAL*8(A-H,O-Z) DIMENSION IT(5),T(0:N/2-1,5) DIMENSION X(M,0:N-1),Y(M,0:N-1) CALL FTTSTF(M,N,X,Y,IT,T) DO J=0,N-1 DO I=1,M X(I,J)=(N/2D0)*X(I,J) END DO END DO END ************************************************************************* * ISPACK FORTRAN SUBROUTINE LIBRARY FOR SCIENTIFIC COMPUTING * * Copyright (C) 1998 Keiichi Ishioka * * * * This library is free software; you can redistribute it and/or * * modify it under the terms of the GNU Library General Public * * License as published by the Free Software Foundation; either * * version 2 of the License, or (at your option) any later version. * * * * This library is distributed in the hope that it will be useful, * * but WITHOUT ANY WARRANTY; without even the implied warranty of * * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU * * Library General Public License for more details. * * * * You should have received a copy of the GNU Library General Public * * License along with this library; if not, write to the Free * * Software Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA. * ************************************************************************* ************************************************************************ SUBROUTINE FTTZUI(N,IT,T) IMPLICIT REAL*8(A-H,O-Z) DIMENSION T(0:N-1,2),IT(5) CALL FTTZLI(N,IT,T) END ************************************************************************ SUBROUTINE FTTZUF(M,N,X,Y,IT,T) IMPLICIT REAL*8(A-H,O-Z) DIMENSION X(M*N,2),Y(M*N,2) DIMENSION T(0:N-1,2),IT(5) F=1D0/N DO I=1,M*N X(I,2)=-X(I,2) END DO CALL FTTZLM(M,N,X,Y,IT,T) DO I=1,M*N X(I,1)= X(I,1)*F X(I,2)=-X(I,2)*F END DO END ************************************************************************ SUBROUTINE FTTZUB(M,N,X,Y,IT,T) IMPLICIT REAL*8(A-H,O-Z) DIMENSION X(M*N,2),Y(M*N,2) DIMENSION T(0:N-1,2),IT(5) CALL FTTZLM(M,N,X,Y,IT,T) END ************************************************************************ SUBROUTINE FTTZLI(N,IT,T) IMPLICIT REAL*8(A-H,O-Z) PARAMETER(PI=3.1415926535897932385D0) DIMENSION T(0:N-1,2),IT(5) IF(N.LE.0) THEN CALL FTDMSG('E','FTTZLI','N MUST BE .GT. 0') END IF J=N DO 20 I=5,2,-1 IT(I)=0 10 CONTINUE K=MOD(J,I) IF(K.EQ.0) THEN IT(I)=IT(I)+1 J=J/I GO TO 10 END IF 20 CONTINUE IF(J.NE.1) THEN CALL FTDMSG('E','FTTZLI','N.NE.(2**P)*(3**Q)*(5**R)') END IF IT(1)=MOD(IT(2)+IT(3)+IT(4)+IT(5),2) IF(IT(1).EQ.1.AND.IT(4).GE.1) THEN IT(2)=IT(2)+2 IT(4)=IT(4)-1 IT(1)=0 END IF DO 30 I=0,N-1 T(I,1)=COS(2*PI*I/N) T(I,2)=SIN(2*PI*I/N) 30 CONTINUE END ************************************************************************ SUBROUTINE FTTZLM(M,N,X,Y,IT,T) IMPLICIT REAL*8(A-H,O-Z) DIMENSION X(M*N,2),Y(M*N,2) DIMENSION T(0:N-1,2),IT(5) K=N L=1 IP=1 DO 10 I=1,IT(1) CALL FTTZL1(M,N,X,Y) IP=-IP 10 CONTINUE DO 20 I=1,IT(2) IF(IP.EQ. 1) CALL FTTZL2(M,K,L,X(1,1),X(1,2),Y(1,1),Y(1,2),T) IF(IP.EQ.-1) CALL FTTZL2(M,K,L,Y(1,1),Y(1,2),X(1,1),X(1,2),T) IP=-IP 20 CONTINUE DO 30 I=1,IT(3) IF(IP.EQ. 1) CALL FTTZL3(M,K,L,X(1,1),X(1,2),Y(1,1),Y(1,2),T) IF(IP.EQ.-1) CALL FTTZL3(M,K,L,Y(1,1),Y(1,2),X(1,1),X(1,2),T) IP=-IP 30 CONTINUE DO 40 I=1,IT(4) IF(IP.EQ. 1) CALL FTTZL4(M,K,L,X(1,1),X(1,2),Y(1,1),Y(1,2),T) IF(IP.EQ.-1) CALL FTTZL4(M,K,L,Y(1,1),Y(1,2),X(1,1),X(1,2),T) IP=-IP 40 CONTINUE DO 50 I=1,IT(5) IF(IP.EQ. 1) CALL FTTZL5(M,K,L,X(1,1),X(1,2),Y(1,1),Y(1,2),T) IF(IP.EQ.-1) CALL FTTZL5(M,K,L,Y(1,1),Y(1,2),X(1,1),X(1,2),T) IP=-IP 50 CONTINUE END ************************************************************************ * FT BY FACTOR 1 *----------------------------------------------------------------------- SUBROUTINE FTTZL1(M,N,X,Y) IMPLICIT REAL*8(A-H,O-Z) DIMENSION X(M*N*2),Y(M*N*2) DO 10 I=1,M*N*2 Y(I)=X(I) 10 CONTINUE END ************************************************************************ * FT BY FACTOR 2 *----------------------------------------------------------------------- SUBROUTINE FTTZL2(M,K,L,A,B,C,D,T) IMPLICIT REAL*8(A-H,O-Z) DIMENSION A(M*L,0:K/2-1,0:1),B(M*L,0:K/2-1,0:1) DIMENSION C(M*L,0:1,0:K/2-1),D(M*L,0:1,0:K/2-1) DIMENSION T(0:L-1,0:K-1,2) DO 20 J=0,K/2-1 DO 10 I=1,M*L C(I,0,J)=A(I,J,0)+A(I,J,1) D(I,0,J)=B(I,J,0)+B(I,J,1) C(I,1,J)= T(0,1*J,1)*(A(I,J,0)-A(I,J,1)) & -T(0,1*J,2)*(B(I,J,0)-B(I,J,1)) D(I,1,J)= T(0,1*J,1)*(B(I,J,0)-B(I,J,1)) & +T(0,1*J,2)*(A(I,J,0)-A(I,J,1)) 10 CONTINUE 20 CONTINUE K=K/2 L=L*2 END ************************************************************************ * FT BY FACTOR 3 *----------------------------------------------------------------------- * C1 = COS(PI/3), S1 = SIN(PI/3) *----------------------------------------------------------------------- SUBROUTINE FTTZL3(M,K,L,A,B,C,D,T) IMPLICIT REAL*8(A-H,O-Z) PARAMETER(C1=0.5D0,S1=0.86602540378443864676D0) DIMENSION A(M*L,0:K/3-1,0:2),B(M*L,0:K/3-1,0:2) DIMENSION C(M*L,0:2,0:K/3-1),D(M*L,0:2,0:K/3-1) DIMENSION T(0:L-1,0:K-1,2) DO 20 J=0,K/3-1 DO 10 I=1,M*L C(I,0,J)=A(I,J,0)+(A(I,J,1)+A(I,J,2)) D(I,0,J)=B(I,J,0)+(B(I,J,1)+B(I,J,2)) C(I,1,J)= & T(0,1*J,1)*( & (A(I,J,0)-C1*(A(I,J,1)+A(I,J,2)))-S1*(B(I,J,1)-B(I,J,2))) & -T(0,1*J,2)*( & (B(I,J,0)-C1*(B(I,J,1)+B(I,J,2)))+S1*(A(I,J,1)-A(I,J,2))) D(I,1,J)= & T(0,1*J,1)*( & (B(I,J,0)-C1*(B(I,J,1)+B(I,J,2)))+S1*(A(I,J,1)-A(I,J,2))) & +T(0,1*J,2)*( & (A(I,J,0)-C1*(A(I,J,1)+A(I,J,2)))-S1*(B(I,J,1)-B(I,J,2))) C(I,2,J)= & T(0,2*J,1)*( & (A(I,J,0)-C1*(A(I,J,1)+A(I,J,2)))+S1*(B(I,J,1)-B(I,J,2))) & -T(0,2*J,2)*( & (B(I,J,0)-C1*(B(I,J,1)+B(I,J,2)))-S1*(A(I,J,1)-A(I,J,2))) D(I,2,J)= & T(0,2*J,1)*( & (B(I,J,0)-C1*(B(I,J,1)+B(I,J,2)))-S1*(A(I,J,1)-A(I,J,2))) & +T(0,2*J,2)*( & (A(I,J,0)-C1*(A(I,J,1)+A(I,J,2)))+S1*(B(I,J,1)-B(I,J,2))) 10 CONTINUE 20 CONTINUE K=K/3 L=L*3 END ************************************************************************ * FT BY FACTOR 4 *----------------------------------------------------------------------- SUBROUTINE FTTZL4(M,K,L,A,B,C,D,T) IMPLICIT REAL*8(A-H,O-Z) DIMENSION A(M*L,0:K/4-1,0:3),B(M*L,0:K/4-1,0:3) DIMENSION C(M*L,0:3,0:K/4-1),D(M*L,0:3,0:K/4-1) DIMENSION T(0:L-1,0:K-1,2) DO 20 J=0,K/4-1 DO 10 I=1,M*L C(I,0,J)=(A(I,J,0)+A(I,J,2))+(A(I,J,1)+A(I,J,3)) D(I,0,J)=(B(I,J,0)+B(I,J,2))+(B(I,J,1)+B(I,J,3)) C(I,2,J)= & T(0,2*J,1)*((A(I,J,0)+A(I,J,2))-(A(I,J,1)+A(I,J,3))) & -T(0,2*J,2)*((B(I,J,0)+B(I,J,2))-(B(I,J,1)+B(I,J,3))) D(I,2,J)= & T(0,2*J,1)*((B(I,J,0)+B(I,J,2))-(B(I,J,1)+B(I,J,3))) & +T(0,2*J,2)*((A(I,J,0)+A(I,J,2))-(A(I,J,1)+A(I,J,3))) C(I,1,J)= & T(0,1*J,1)*((A(I,J,0)-A(I,J,2))-(B(I,J,1)-B(I,J,3))) & -T(0,1*J,2)*((B(I,J,0)-B(I,J,2))+(A(I,J,1)-A(I,J,3))) D(I,1,J)= & T(0,1*J,1)*((B(I,J,0)-B(I,J,2))+(A(I,J,1)-A(I,J,3))) & +T(0,1*J,2)*((A(I,J,0)-A(I,J,2))-(B(I,J,1)-B(I,J,3))) C(I,3,J)= & T(0,3*J,1)*((A(I,J,0)-A(I,J,2))+(B(I,J,1)-B(I,J,3))) & -T(0,3*J,2)*((B(I,J,0)-B(I,J,2))-(A(I,J,1)-A(I,J,3))) D(I,3,J)= & T(0,3*J,1)*((B(I,J,0)-B(I,J,2))-(A(I,J,1)-A(I,J,3))) & +T(0,3*J,2)*((A(I,J,0)-A(I,J,2))+(B(I,J,1)-B(I,J,3))) 10 CONTINUE 20 CONTINUE K=K/4 L=L*4 END ************************************************************************ * FT BY FACTOR 5 *----------------------------------------------------------------------- * C1 = COS(PI/5), S1 = SIN(PI/5), C2 = COS(2*PI/5), S2 = SIN(2*PI/5) *----------------------------------------------------------------------- SUBROUTINE FTTZL5(M,K,L,A,B,C,D,T) IMPLICIT REAL*8(A-H,O-Z) PARAMETER(C1=0.80901699437494742410D0,S1=0.58778525229247312917D0) PARAMETER(C2=0.30901699437494742410D0,S2=0.95105651629515357212D0) DIMENSION A(M*L,0:K/5-1,0:4),B(M*L,0:K/5-1,0:4) DIMENSION C(M*L,0:4,0:K/5-1),D(M*L,0:4,0:K/5-1) DIMENSION T(0:L-1,0:K-1,2) DO 20 J=0,K/5-1 DO 10 I=1,M*L C(I,0,J)=A(I,J,0)+(A(I,J,1)+A(I,J,4))+(A(I,J,2)+A(I,J,3)) D(I,0,J)=B(I,J,0)+(B(I,J,1)+B(I,J,4))+(B(I,J,2)+B(I,J,3)) C(I,1,J)= & T(0,1*J,1)*( & (A(I,J,0)+(C2*(A(I,J,1)+A(I,J,4))-C1*(A(I,J,2)+A(I,J,3))) & -(S2*(B(I,J,1)-B(I,J,4))+S1*(B(I,J,2)-B(I,J,3))))) & -T(0,1*J,2)*( & (B(I,J,0)+(C2*(B(I,J,1)+B(I,J,4))-C1*(B(I,J,2)+B(I,J,3))) & +(S2*(A(I,J,1)-A(I,J,4))+S1*(A(I,J,2)-A(I,J,3))))) D(I,1,J)= & T(0,1*J,1)*( & (B(I,J,0)+(C2*(B(I,J,1)+B(I,J,4))-C1*(B(I,J,2)+B(I,J,3))) & +(S2*(A(I,J,1)-A(I,J,4))+S1*(A(I,J,2)-A(I,J,3))))) & +T(0,1*J,2)*( & (A(I,J,0)+(C2*(A(I,J,1)+A(I,J,4))-C1*(A(I,J,2)+A(I,J,3))) & -(S2*(B(I,J,1)-B(I,J,4))+S1*(B(I,J,2)-B(I,J,3))))) C(I,4,J)= & T(0,4*J,1)*( & (A(I,J,0)+(C2*(A(I,J,1)+A(I,J,4))-C1*(A(I,J,2)+A(I,J,3))) & +(S2*(B(I,J,1)-B(I,J,4))+S1*(B(I,J,2)-B(I,J,3))))) & -T(0,4*J,2)*( & (B(I,J,0)+(C2*(B(I,J,1)+B(I,J,4))-C1*(B(I,J,2)+B(I,J,3))) & -(S2*(A(I,J,1)-A(I,J,4))+S1*(A(I,J,2)-A(I,J,3))))) D(I,4,J)= & T(0,4*J,1)*( & (B(I,J,0)+(C2*(B(I,J,1)+B(I,J,4))-C1*(B(I,J,2)+B(I,J,3))) & -(S2*(A(I,J,1)-A(I,J,4))+S1*(A(I,J,2)-A(I,J,3))))) & +T(0,4*J,2)*( & (A(I,J,0)+(C2*(A(I,J,1)+A(I,J,4))-C1*(A(I,J,2)+A(I,J,3))) & +(S2*(B(I,J,1)-B(I,J,4))+S1*(B(I,J,2)-B(I,J,3))))) C(I,2,J)= & T(0,2*J,1)*( & (A(I,J,0)-(C1*(A(I,J,1)+A(I,J,4))-C2*(A(I,J,2)+A(I,J,3))) & -(S1*(B(I,J,1)-B(I,J,4))-S2*(B(I,J,2)-B(I,J,3))))) & -T(0,2*J,2)*( & (B(I,J,0)-(C1*(B(I,J,1)+B(I,J,4))-C2*(B(I,J,2)+B(I,J,3))) & +(S1*(A(I,J,1)-A(I,J,4))-S2*(A(I,J,2)-A(I,J,3))))) D(I,2,J)= & T(0,2*J,1)*( & (B(I,J,0)-(C1*(B(I,J,1)+B(I,J,4))-C2*(B(I,J,2)+B(I,J,3))) & +(S1*(A(I,J,1)-A(I,J,4))-S2*(A(I,J,2)-A(I,J,3))))) & +T(0,2*J,2)*( & (A(I,J,0)-(C1*(A(I,J,1)+A(I,J,4))-C2*(A(I,J,2)+A(I,J,3))) & -(S1*(B(I,J,1)-B(I,J,4))-S2*(B(I,J,2)-B(I,J,3))))) C(I,3,J)= & T(0,3*J,1)*( & (A(I,J,0)-(C1*(A(I,J,1)+A(I,J,4))-C2*(A(I,J,2)+A(I,J,3))) & +(S1*(B(I,J,1)-B(I,J,4))-S2*(B(I,J,2)-B(I,J,3))))) & -T(0,3*J,2)*( & (B(I,J,0)-(C1*(B(I,J,1)+B(I,J,4))-C2*(B(I,J,2)+B(I,J,3))) & -(S1*(A(I,J,1)-A(I,J,4))-S2*(A(I,J,2)-A(I,J,3))))) D(I,3,J)= & T(0,3*J,1)*( & (B(I,J,0)-(C1*(B(I,J,1)+B(I,J,4))-C2*(B(I,J,2)+B(I,J,3))) & -(S1*(A(I,J,1)-A(I,J,4))-S2*(A(I,J,2)-A(I,J,3))))) & +T(0,3*J,2)*( & (A(I,J,0)-(C1*(A(I,J,1)+A(I,J,4))-C2*(A(I,J,2)+A(I,J,3))) & +(S1*(B(I,J,1)-B(I,J,4))-S2*(B(I,J,2)-B(I,J,3))))) 10 CONTINUE 20 CONTINUE K=K/5 L=L*5 END