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