website/archive/binaries/mac/src/disembl/Tisean_3.0.1/source_f/slatec/rfftf1.f

   1 *DECK RFFTF1
   2       SUBROUTINE RFFTF1 (N, C, CH, WA, IFAC)
   3 C***BEGIN PROLOGUE  RFFTF1
   4 C***PURPOSE  Compute the forward transform of a real, periodic sequence.
   5 C***LIBRARY   SLATEC (FFTPACK)
   6 C***CATEGORY  J1A1
   7 C***TYPE      SINGLE PRECISION (RFFTF1-S, CFFTF1-C)
   8 C***KEYWORDS  FFTPACK, FOURIER TRANSFORM
   9 C***AUTHOR  Swarztrauber, P. N., (NCAR)
  10 C***DESCRIPTION
  11 C
  12 C   Subroutine RFFTF1 computes the Fourier coefficients of a real
  13 C   periodic sequence (Fourier analysis).  The transform is defined
  14 C   below at output parameter C.
  15 C
  16 C   The arrays WA and IFAC which are used by subroutine RFFTB1 must be
  17 C   initialized by calling subroutine RFFTI1.
  18 C
  19 C   Input Arguments
  20 C
  21 C   N       the length of the array R to be transformed.  The method
  22 C           is most efficient when N is a product of small primes.
  23 C           N may change so long as different work arrays are provided.
  24 C
  25 C   C       a real array of length N which contains the sequence
  26 C           to be transformed.
  27 C
  28 C   CH      a real work array of length at least N.
  29 C
  30 C   WA      a real work array which must be dimensioned at least N.
  31 C
  32 C   IFAC    an integer work array which must be dimensioned at least 15.
  33 C
  34 C           The WA and IFAC arrays must be initialized by calling
  35 C           subroutine RFFTI1, and different WA and IFAC arrays must be
  36 C           used for each different value of N.  This initialization
  37 C           does not have to be repeated so long as N remains unchanged.
  38 C           Thus subsequent transforms can be obtained faster than the
  39 C           first.  The same WA and IFAC arrays can be used by RFFTF1
  40 C           and RFFTB1.
  41 C
  42 C   Output Argument
  43 C
  44 C   C       C(1) = the sum from I=1 to I=N of R(I)
  45 C
  46 C           If N is even set L = N/2; if N is odd set L = (N+1)/2
  47 C
  48 C             then for K = 2,...,L
  49 C
  50 C                C(2*K-2) = the sum from I = 1 to I = N of
  51 C
  52 C                     C(I)*COS((K-1)*(I-1)*2*PI/N)
  53 C
  54 C                C(2*K-1) = the sum from I = 1 to I = N of
  55 C
  56 C                    -C(I)*SIN((K-1)*(I-1)*2*PI/N)
  57 C
  58 C           If N is even
  59 C
  60 C                C(N) = the sum from I = 1 to I = N of
  61 C
  62 C                     (-1)**(I-1)*C(I)
  63 C
  64 C   Notes:  This transform is unnormalized since a call of RFFTF1
  65 C           followed by a call of RFFTB1 will multiply the input
  66 C           sequence by N.
  67 C
  68 C           WA and IFAC contain initialization calculations which must
  69 C           not be destroyed between calls of subroutine RFFTF1 or
  70 C           RFFTB1.
  71 C
  72 C***REFERENCES  P. N. Swarztrauber, Vectorizing the FFTs, in Parallel
  73 C                 Computations (G. Rodrigue, ed.), Academic Press,
  74 C                 1982, pp. 51-83.
  75 C***ROUTINES CALLED  RADF2, RADF3, RADF4, RADF5, RADFG
  76 C***REVISION HISTORY  (YYMMDD)
  77 C   790601  DATE WRITTEN
  78 C   830401  Modified to use SLATEC library source file format.
  79 C   860115  Modified by Ron Boisvert to adhere to Fortran 77 by
  80 C           changing dummy array size declarations (1) to (*).
  81 C   881128  Modified by Dick Valent to meet prologue standards.
  82 C   891214  Prologue converted to Version 4.0 format.  (BAB)
  83 C   900131  Routine changed from subsidiary to user-callable.  (WRB)
  84 C   920501  Reformatted the REFERENCES section.  (WRB)
  85 C***END PROLOGUE  RFFTF1
  86       DIMENSION CH(*), C(*), WA(*), IFAC(*)
  87 C***FIRST EXECUTABLE STATEMENT  RFFTF1
  88       NF = IFAC(2)
  89       NA = 1
  90       L2 = N
  91       IW = N
  92       DO 111 K1=1,NF
  93          KH = NF-K1
  94          IP = IFAC(KH+3)
  95          L1 = L2/IP
  96          IDO = N/L2
  97          IDL1 = IDO*L1
  98          IW = IW-(IP-1)*IDO
  99          NA = 1-NA
 100          IF (IP .NE. 4) GO TO 102
 101          IX2 = IW+IDO
 102          IX3 = IX2+IDO
 103          IF (NA .NE. 0) GO TO 101
 104          CALL RADF4 (IDO,L1,C,CH,WA(IW),WA(IX2),WA(IX3))
 105          GO TO 110
 106   101    CALL RADF4 (IDO,L1,CH,C,WA(IW),WA(IX2),WA(IX3))
 107          GO TO 110
 108   102    IF (IP .NE. 2) GO TO 104
 109          IF (NA .NE. 0) GO TO 103
 110          CALL RADF2 (IDO,L1,C,CH,WA(IW))
 111          GO TO 110
 112   103    CALL RADF2 (IDO,L1,CH,C,WA(IW))
 113          GO TO 110
 114   104    IF (IP .NE. 3) GO TO 106
 115          IX2 = IW+IDO
 116          IF (NA .NE. 0) GO TO 105
 117          CALL RADF3 (IDO,L1,C,CH,WA(IW),WA(IX2))
 118          GO TO 110
 119   105    CALL RADF3 (IDO,L1,CH,C,WA(IW),WA(IX2))
 120          GO TO 110
 121   106    IF (IP .NE. 5) GO TO 108
 122          IX2 = IW+IDO
 123          IX3 = IX2+IDO
 124          IX4 = IX3+IDO
 125          IF (NA .NE. 0) GO TO 107
 126          CALL RADF5 (IDO,L1,C,CH,WA(IW),WA(IX2),WA(IX3),WA(IX4))
 127          GO TO 110
 128   107    CALL RADF5 (IDO,L1,CH,C,WA(IW),WA(IX2),WA(IX3),WA(IX4))
 129          GO TO 110
 130   108    IF (IDO .EQ. 1) NA = 1-NA
 131          IF (NA .NE. 0) GO TO 109
 132          CALL RADFG (IDO,IP,L1,IDL1,C,C,C,CH,CH,WA(IW))
 133          NA = 1
 134          GO TO 110
 135   109    CALL RADFG (IDO,IP,L1,IDL1,CH,CH,CH,C,C,WA(IW))
 136          NA = 0
 137   110    L2 = L1
 138   111 CONTINUE
 139       IF (NA .EQ. 1) RETURN
 140       DO 112 I=1,N
 141          C(I) = CH(I)
 142   112 CONTINUE
 143       RETURN
 144       END