spffte.f

Go to the documentation of this file.

C> @file
C>
C> Perform multiple fast fourier transforms
C> @author IREDELL @date 96-02-20
 
C> This subprogram performs multiple fast fourier transforms
C> between complex amplitudes in fourier space and real values
C> in cyclic physical space.
C> Subprogram SPFFTE must be invoked first with IDIR=0
C> to initialize trigonemetric data.  use subprogram SPFFT1
C> to perform an fft without previous initialization.
C> This version invokes the ibm essl fft.
C>
C> PROGRAM HISTORY LOG:
C> - 1998-12-18  IREDELL
C> - 2012-11-12  MIRVIS -fixing hard-wired types problem on Intel/Linux 
C>
C> @param IMAX     - INTEGER NUMBER OF VALUES IN THE CYCLIC PHYSICAL SPACE
C>                (SEE LIMITATIONS ON IMAX IN REMARKS BELOW.)
C> @param INCW     - INTEGER FIRST DIMENSION OF THE COMPLEX AMPLITUDE ARRAY
C>                (INCW >= IMAX/2+1)
C> @param INCG     - INTEGER FIRST DIMENSION OF THE REAL VALUE ARRAY
C>                (INCG >= IMAX)
C> @param KMAX     - INTEGER NUMBER OF TRANSFORMS TO PERFORM
C> @param[out] W        - COMPLEX(INCW,KMAX) COMPLEX AMPLITUDES IF IDIR>0
C> @param[out] G        - REAL(INCG,KMAX) REAL VALUES IF IDIR<0
C> @param IDIR     - INTEGER DIRECTION FLAG
C> - IDIR=0 TO INITIALIZE TRIGONOMETRIC DATA
C> - IDIR>0 TO TRANSFORM FROM FOURIER TO PHYSICAL SPACE
C> - IDIR<0 TO TRANSFORM FROM PHYSICAL TO FOURIER SPACE
C> @param[out] AFFT       REAL(8) (50000+4*IMAX) AUXILIARY ARRAY IF IDIR<>0
C>
C> SUBPROGRAMS CALLED:
C>   -scrft()        IBM ESSL COMPLEX TO REAL FOURIER TRANSFORM
C>   -dcrft()        IBM ESSL COMPLEX TO REAL FOURIER TRANSFORM
C>   -srcft()        IBM ESSL REAL TO COMPLEX FOURIER TRANSFORM
C>   -drcft()        IBM ESSL REAL TO COMPLEX FOURIER TRANSFORM
C>
C> @note The restrictions on IMAX are that it must be a multiple
C> of 1 to 25 factors of two, up to 2 factors of three,
C> and up to 1 factor of five, seven and eleven.
C>
C> If IDIR=0, then W and G need not contain any valid data.
C> The other parameters must be supplied and cannot change
C> in succeeding calls until the next time it is called with IDIR=0.
C>
C> This subprogram is thread-safe.
C>      
      SUBROUTINE spffte(IMAX,INCW,INCG,KMAX,W,G,IDIR,AFFT)
        IMPLICIT NONE
        INTEGER,INTENT(IN):: IMAX,INCW,INCG,KMAX,IDIR
        REAL,INTENT(INOUT):: W(2*INCW,KMAX)
        REAL,INTENT(INOUT):: G(INCG,KMAX)
        REAL(8),INTENT(INOUT):: AFFT(50000+4*IMAX)
        INTEGER:: INIT,INC2X,INC2Y,N,M,ISIGN,NAUX1,NAUX2,NAUX3
C ==EM==       ^(4)
        REAL:: SCALE
        REAL(8):: AUX2(20000+2*IMAX),AUX3
        INTEGER:: IACR,IARC
C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
        naux1=25000+2*imax
        naux2=20000+2*imax
        naux3=1
        iacr=1
        iarc=1+naux1
C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
C  INITIALIZATION.
C  FILL AUXILIARY ARRAYS WITH TRIGONOMETRIC DATA
        SELECT CASE(idir)
          CASE(0)
            init=1
            inc2x=incw
            inc2y=incg
            n=imax
            m=kmax
            isign=-1
            scale=1.
            IF(digits(1.).LT.digits(1._8)) THEN
              CALL scrft(init,w,inc2x,g,inc2y,n,m,isign,scale,
     &                   afft(iacr),naux1,aux2,naux2,aux3,naux3)
            ELSE
              CALL dcrft(init,w,inc2x,g,inc2y,n,m,isign,scale,
     &                   afft(iacr),naux1,aux2,naux2)
            ENDIF
            init=1
            inc2x=incg
            inc2y=incw
            n=imax
            m=kmax
            isign=+1
            scale=1./imax
            IF(digits(1.).LT.digits(1._8)) THEN
              CALL srcft(init,g,inc2x,w,inc2y,n,m,isign,scale,
     &                   afft(iarc),naux1,aux2,naux2,aux3,naux3)
            ELSE
              CALL drcft(init,g,inc2x,w,inc2y,n,m,isign,scale,
     &                   afft(iarc),naux1,aux2,naux2)
            ENDIF
C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
C  FOURIER TO PHYSICAL TRANSFORM.
          CASE(1:)
            init=0
            inc2x=incw
            inc2y=incg
            n=imax
            m=kmax
            isign=-1
            scale=1.
            IF(digits(1.).LT.digits(1._8)) THEN
              CALL scrft(init,w,inc2x,g,inc2y,n,m,isign,scale,
     &                   afft(iacr),naux1,aux2,naux2,aux3,naux3)
            ELSE
              CALL dcrft(init,w,inc2x,g,inc2y,n,m,isign,scale,
     &                   afft(iacr),naux1,aux2,naux2)
            ENDIF
C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
C  PHYSICAL TO FOURIER TRANSFORM.
          CASE(:-1)
            init=0
            inc2x=incg
            inc2y=incw
            n=imax
            m=kmax
            isign=+1
            scale=1./imax
            IF(digits(1.).LT.digits(1._8)) THEN
              CALL srcft(init,g,inc2x,w,inc2y,n,m,isign,scale,
     &                   afft(iarc),naux1,aux2,naux2,aux3,naux3)
            ELSE
              CALL drcft(init,g,inc2x,w,inc2y,n,m,isign,scale,
     &                   afft(iarc),naux1,aux2,naux2)
            ENDIF
        END SELECT
      END SUBROUTINE