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