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