r3fft.f

FFT源代码汇集（各种FFT源代码）

      SUBROUTINE R3FFT(C,ID,NL,NM,NN,WL,WM,WN,IOPT,ISIG,IORD,IWORK,IERR)
****PURPOSE:
*       THIS ROUTINE PERFORMS A 3-DIMENSIONAL REAL FOURIER TRANSFORM,
*       OF ORDER NL*NM*NN  .
****USAGE:
*       THE USER IS EXPECTED TO PROVIDE THE DATA IN A 3-DIMENSIONAL
*       REAL ARRAY C, DIMENSIONED IN THE CALLING PROGRAM C(ID,NM,NN);
*       ID HAS TO BE AN EVEN INTEGER, EQUAL TO NL+2.
*       FOR OUTPUT DATA ARRENGEMENT SEE NOTES TO R2FFT HERE ABOVE.
*       THIS ROUTINE IS
*       INTENDED FOR REPEATED USAGE, THUS SEPARATE SET-UP AND
*       AND OPERATING CALLS ARE AVAILABLE: THE USER SHOULD IN ANY CASE
*       PERFORM A SET-UP CALL (ISIG=0) PASSING THE PARAMETERS BEFORE
*       PERFORMING AN ACTUAL TRANSFORM ( ISIG= +1 OR -1 ); THE USER CAN
*       CHOOSE WHETHER TO OBTAIN THE RESULTS OF THE DIRECT TRANSFORM
*       IN NATURAL ORDER (ISIG=-1,IORD=1) OR LEAVE THEM IN THE
*       BIT-REVERSED  ORDER( ISIG=-1,IORD=0); THIS CHOICE SAVES
*       SOME COMPUTER TIME, AND IT IS RECOMMENDED IN CASES DISCUSSED
*       IN THE LONG WRITE-UP. ANALOGOUSLY, THE INVERSE TRANSFORM ACCEPTS
*       INPUT ( PLEASE NOTE| ) DATA IN NATURAL ORDER ( ISIG=1,IORD=1),
*       OR DATA ALREADY SUBJECTED TO A BIT-REVERSAL PERMUTATION( ISIG=1
*       IORD=0).
*       A SPECIAL TREATMENT IS AVAILABLE TO SPEED UP THE TRANSFORM OF
*       SMALL MATRICES. THIS TREATMENT IS ACTIVATED BY THE FLAG IOPT. IN
*       THIS CASE THE TABLES FOR THE SECOND DIMENSION ( WM ) ARE LARGER,
*       BUT THE INCREASE IN PERFORMANCE IS SUBSTANTIAL WHEN NM<32.
****ARGUMENTS :
*       INPUT :
*       C : ARRAY TO BE TRANSFORMED; DECLARED COMPLEX C(ID,NM,NN) IN THE
*           CALLING PROGRAM;
*       ID : FIRST DIMENSION OF C IN THE CALLING PROGRAM
*            IT HAS TO BE AN EVEN INTEGER .GE. NL+2.
*       ISIG : OPTION FLAG : ISIG=0 : SET-UP RUN, C NOT USED
*                            ISIG=-1: DIRECT TRANSFORM
*                            ISIG=+1: INVERSE TRANSFORM
*       WL,WM,WN : INTEGER ARRAYS,USED TO HOST TABLES FOR THE TRANSFORMS
*               DIMENSIONED IN THE CALLING PROGRAM AT LEAST (6*NL+14)
*               (4*NM+14) AND (4*NN+14) RESPECTIVELY; IF
*               IOPT=1, WM MUST BE DIMENSIONED AT LEAST 4*NM*(ID/2+1)+14
*               IF ISIG.NE.0, THEY ARE ASSUMED TO HAVE BEEN SET BY A
*               PREVIOUS CALL WITH ISIG=0 AND OTHER ARGUMENTS EQUAL, AND
*               NEVER HAVE BEEN MODIFIED ;
*               WHEN THE CORRESPONDING ORDERS ARE EQUAL, THEY DO NOT
*               NEED TO BE DISTINCT
*       NL : ORDER OF THE TRANSFORM ALONG THE COLUMNS OF C
*            IT HAS TO BE AN EVEN INTEGER.
*       NM : ORDER OF THE TRANSFORM ALONG THE ROWS OF C
*       NN : ORDER OF THE TRANSFORM ALONG THE THIRD DIMENSION OF C
*       IOPT : OPTION FLAG : =0 : NORMAL TREATMENT
*                            =1 : SPECIAL TREATMENT FOR IMPROVING
*                                 VECTORIZATION ON MATRICES WITH
*                                 SMALL NL; REQUIRES LONG WM(SEE);IF
*                                 REQUESTED, MUST BE PRESENT IN BOTH
*                                 THE SET-UP AND TRANSFORM CALLS;
*       IORD : OPTION FLAG : =1 : OUTPUT IN NATURAL ORDER (ISIG=-1)
*                                 INPUT IN NATURAL ORDER  (ISIG=+1)
*                            =0 : OUTPUT IN BIT-REVERSED ORDER(ISIG=-1)
*                                 INPUT IN BIT-REVERSED ORDER(ISIG=+1)
*       IWORK : INTEGER ARRAY, USED AS WORK AREA FOR REORDERING IF
*               IORD=1; IT MUST BE AT LEAST MAX(NL,NM,NN) WORDS LONG.
*
*        OUTPUT :
*       C : TRANSFORMED ARRAY
*       WL, WM, WN : ONLY IF ISIG=0, WL,WM AND WN ARE FILLED WITH THE
*                     APPROPRIATE TABLES
*       IWORK : UNDEFINED
*       IERR  : ERROR CODE : =0 : SUCCESSFUL
*                            =1 : WRONG ID PARAMETER
*                            =2 : PRIME FACTORS DIFFERENT FROM 2,3,5
*                                 ARE PRESENT IN DATA DIMENSIONS
*                            =3 : TABLES NOT CORRECTLY INITIALIZED
*                            =4 : FIRST DIMENSION IS AN ODD NUMBER
*
      COMPLEX C(*)
      INTEGER WL(-14:*),WM(-14:*),WN(-14:*)
      INTEGER  IWORK(*)
*
*
      INTEGER IDERR,FACERR,TBERR,ODDERR
      PARAMETER (IDERR=1,FACERR=2,TBERR=3,ODDERR=4)
*
      IF(ID.LT.NL+2)THEN
        IERR=IDERR
        RETURN
      ENDIF
      NL1=NL/2
      IF(NL1*2.NE.NL)THEN
        IERR=ODDERR
        RETURN
      ENDIF
      IERR=0
 
      NMPN=NM*NN
*
*
      IF(ISIG.EQ.0) THEN
        CALL MFFTP(NM,WM,ID/2*IOPT,IERR)
        IF(IERR.NE.0)RETURN
 
        CALL MFFTRP(NL,WL(4*NL))
        IF(NL1.NE.NM) THEN
          CALL MFFTP(NL1,WL,0,IERR)
          IF(IERR.NE.0)RETURN
        ELSE
          CALL MFFTZ0(WM,1,4*NM+14,WL,1)
        ENDIF
*
        IF(NM.EQ.NN) THEN
          CALL MFFTZ0(WM,1,4*NM+14,WN,1)
        ELSE IF(NN.EQ.NL1) THEN
          CALL MFFTZ0(WL,1,4*NL1+14,WN,1)
        ELSE
          CALL MFFTP(NN,WN,0,IERR)
          IF(IERR.NE.0)RETURN
        ENDIF
        RETURN
*
      ELSE   IF(ISIG.GT.0) THEN
*
        IF(IORD.NE.0) THEN
          CALL MFFTOM(C,ID/2,ID/2*NM,1,NM,NN,NL1+1,WM(NM*3),IWORK)
          CALL MFFTOV(C,ID/2*NM,1,NN,ID/2*NM,WN(NN*3),IWORK)
        ENDIF
*
        CALL MFFTIV(C,ID/2*NM,1,NN,ID/2*NM,WN,IERR)
        IF(IERR.NE.0)RETURN
*
        IF(IOPT.EQ.0) THEN
          CALL MFFTIM(C,ID/2,ID/2*NM,1,NM,NN,NL1+1,WM,IERR)
          IF(IERR.NE.0)RETURN
        ELSE
          CALL MFFTIS(C,ID/2,ID/2*NM,1,NM,NN,NL1+1,WM,IERR)
          IF(IERR.NE.0)RETURN
        ENDIF
*
        CALL MFFTRI(C,1,ID/2,NL1,NMPN,WL(4*NL))
        CALL MFFTOV(C,1,ID/2,NL1,NMPN,WL(NL1*3),IWORK)
        CALL MFFTIV(C,1,ID/2,NL1,NMPN,WL,IERR)
        IF(IERR.NE.0)RETURN
*
*
      ELSE
*
*
        CALL MFFTDV(C,1,ID/2,NL1,NMPN,WL,IERR)
        IF(IERR.NE.0)RETURN
        CALL MFFTOV(C,1,ID/2,NL1,NMPN,WL(NL1*2),IWORK)
        CALL MFFTRD(C,1,ID/2,NL1,NMPN,WL(4*NL))
*
*
        IF(IOPT.EQ.0) THEN
          CALL MFFTDM(C,ID/2,ID/2*NM,1,NM,NN,NL1+1,WM,IERR)
          IF(IERR.NE.0)RETURN
        ELSE
          CALL MFFTDS(C,ID/2,ID/2*NM,1,NM,NN,NL1+1,WM,IERR)
          IF(IERR.NE.0)RETURN
        ENDIF
*
        CALL MFFTDV(C,ID/2*NM,1,NN,ID/2*NM,WN,IERR)
        IF(IERR.NE.0)RETURN
*
        IF(IORD.NE.0) THEN
          CALL MFFTOV(C,ID/2*NM,1,NN,ID/2*NM,WN(NN*2),IWORK)
          CALL MFFTOM(C,ID/2,ID/2*NM,1,NM,NN,NL1+1,WM(NM*2),IWORK)
        ENDIF
*
      ENDIF
*
      END