1*DECK RFFTF1
2      SUBROUTINE RFFTF1 (N, C, CH, WA, IFAC)
3C***BEGIN PROLOGUE  RFFTF1
4C***PURPOSE  Compute the forward transform of a real, periodic sequence.
5C***LIBRARY   SLATEC (FFTPACK)
6C***CATEGORY  J1A1
7C***TYPE      SINGLE PRECISION (RFFTF1-S, CFFTF1-C)
8C***KEYWORDS  FFTPACK, FOURIER TRANSFORM
9C***AUTHOR  Swarztrauber, P. N., (NCAR)
10C***DESCRIPTION
11C
12C   Subroutine RFFTF1 computes the Fourier coefficients of a real
13C   periodic sequence (Fourier analysis).  The transform is defined
14C   below at output parameter C.
15C
16C   The arrays WA and IFAC which are used by subroutine RFFTB1 must be
17C   initialized by calling subroutine RFFTI1.
18C
19C   Input Arguments
20C
21C   N       the length of the array R to be transformed.  The method
22C           is most efficient when N is a product of small primes.
23C           N may change so long as different work arrays are provided.
24C
25C   C       a real array of length N which contains the sequence
26C           to be transformed.
27C
28C   CH      a real work array of length at least N.
29C
30C   WA      a real work array which must be dimensioned at least N.
31C
32C   IFAC    an integer work array which must be dimensioned at least 15.
33C
34C           The WA and IFAC arrays must be initialized by calling
35C           subroutine RFFTI1, and different WA and IFAC arrays must be
36C           used for each different value of N.  This initialization
37C           does not have to be repeated so long as N remains unchanged.
38C           Thus subsequent transforms can be obtained faster than the
39C           first.  The same WA and IFAC arrays can be used by RFFTF1
40C           and RFFTB1.
41C
42C   Output Argument
43C
44C   C       C(1) = the sum from I=1 to I=N of R(I)
45C
46C           If N is even set L = N/2; if N is odd set L = (N+1)/2
47C
48C             then for K = 2,...,L
49C
50C                C(2*K-2) = the sum from I = 1 to I = N of
51C
52C                     C(I)*COS((K-1)*(I-1)*2*PI/N)
53C
54C                C(2*K-1) = the sum from I = 1 to I = N of
55C
56C                    -C(I)*SIN((K-1)*(I-1)*2*PI/N)
57C
58C           If N is even
59C
60C                C(N) = the sum from I = 1 to I = N of
61C
62C                     (-1)**(I-1)*C(I)
63C
64C   Notes:  This transform is unnormalized since a call of RFFTF1
65C           followed by a call of RFFTB1 will multiply the input
66C           sequence by N.
67C
68C           WA and IFAC contain initialization calculations which must
69C           not be destroyed between calls of subroutine RFFTF1 or
70C           RFFTB1.
71C
72C***REFERENCES  P. N. Swarztrauber, Vectorizing the FFTs, in Parallel
73C                 Computations (G. Rodrigue, ed.), Academic Press,
74C                 1982, pp. 51-83.
75C***ROUTINES CALLED  RADF2, RADF3, RADF4, RADF5, RADFG
76C***REVISION HISTORY  (YYMMDD)
77C   790601  DATE WRITTEN
78C   830401  Modified to use SLATEC library source file format.
79C   860115  Modified by Ron Boisvert to adhere to Fortran 77 by
80C           changing dummy array size declarations (1) to (*).
81C   881128  Modified by Dick Valent to meet prologue standards.
82C   891214  Prologue converted to Version 4.0 format.  (BAB)
83C   900131  Routine changed from subsidiary to user-callable.  (WRB)
84C   920501  Reformatted the REFERENCES section.  (WRB)
85C***END PROLOGUE  RFFTF1
86      DIMENSION CH(*), C(*), WA(*), IFAC(*)
87C***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
145