rffti1.f 3.4 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110
  1. *DECK RFFTI1
  2. SUBROUTINE RFFTI1 (N, WA, IFAC)
  3. C***BEGIN PROLOGUE RFFTI1
  4. C***PURPOSE Initialize a real and an integer work array for RFFTF1 and
  5. C RFFTB1.
  6. C***LIBRARY SLATEC (FFTPACK)
  7. C***CATEGORY J1A1
  8. C***TYPE SINGLE PRECISION (RFFTI1-S, CFFTI1-C)
  9. C***KEYWORDS FFTPACK, FOURIER TRANSFORM
  10. C***AUTHOR Swarztrauber, P. N., (NCAR)
  11. C***DESCRIPTION
  12. C
  13. C Subroutine RFFTI1 initializes the work arrays WA and IFAC which are
  14. C used in both RFFTF1 and RFFTB1. The prime factorization of N and a
  15. C tabulation of the trigonometric functions are computed and stored in
  16. C IFAC and WA, respectively.
  17. C
  18. C Input Argument
  19. C
  20. C N the length of the sequence to be transformed.
  21. C
  22. C Output Arguments
  23. C
  24. C WA a real work array which must be dimensioned at least N.
  25. C
  26. C IFAC an integer work array which must be dimensioned at least 15.
  27. C
  28. C The same work arrays can be used for both RFFTF1 and RFFTB1 as long
  29. C as N remains unchanged. Different WA and IFAC arrays are required
  30. C for different values of N. The contents of WA and IFAC must not be
  31. C changed between calls of RFFTF1 or RFFTB1.
  32. C
  33. C***REFERENCES P. N. Swarztrauber, Vectorizing the FFTs, in Parallel
  34. C Computations (G. Rodrigue, ed.), Academic Press,
  35. C 1982, pp. 51-83.
  36. C***ROUTINES CALLED (NONE)
  37. C***REVISION HISTORY (YYMMDD)
  38. C 790601 DATE WRITTEN
  39. C 830401 Modified to use SLATEC library source file format.
  40. C 860115 Modified by Ron Boisvert to adhere to Fortran 77 by
  41. C (a) changing dummy array size declarations (1) to (*),
  42. C (b) changing references to intrinsic function FLOAT
  43. C to REAL, and
  44. C (c) changing definition of variable TPI by using
  45. C FORTRAN intrinsic functions instead of DATA
  46. C statements.
  47. C 881128 Modified by Dick Valent to meet prologue standards.
  48. C 890531 Changed all specific intrinsics to generic. (WRB)
  49. C 891214 Prologue converted to Version 4.0 format. (BAB)
  50. C 900131 Routine changed from subsidiary to user-callable. (WRB)
  51. C 920501 Reformatted the REFERENCES section. (WRB)
  52. C***END PROLOGUE RFFTI1
  53. DIMENSION WA(*), IFAC(*), NTRYH(4)
  54. SAVE NTRYH
  55. DATA NTRYH(1),NTRYH(2),NTRYH(3),NTRYH(4)/4,2,3,5/
  56. C***FIRST EXECUTABLE STATEMENT RFFTI1
  57. NL = N
  58. NF = 0
  59. J = 0
  60. 101 J = J+1
  61. IF (J-4) 102,102,103
  62. 102 NTRY = NTRYH(J)
  63. GO TO 104
  64. 103 NTRY = NTRY+2
  65. 104 NQ = NL/NTRY
  66. NR = NL-NTRY*NQ
  67. IF (NR) 101,105,101
  68. 105 NF = NF+1
  69. IFAC(NF+2) = NTRY
  70. NL = NQ
  71. IF (NTRY .NE. 2) GO TO 107
  72. IF (NF .EQ. 1) GO TO 107
  73. DO 106 I=2,NF
  74. IB = NF-I+2
  75. IFAC(IB+2) = IFAC(IB+1)
  76. 106 CONTINUE
  77. IFAC(3) = 2
  78. 107 IF (NL .NE. 1) GO TO 104
  79. IFAC(1) = N
  80. IFAC(2) = NF
  81. TPI = 8.*ATAN(1.)
  82. ARGH = TPI/N
  83. IS = 0
  84. NFM1 = NF-1
  85. L1 = 1
  86. IF (NFM1 .EQ. 0) RETURN
  87. DO 110 K1=1,NFM1
  88. IP = IFAC(K1+2)
  89. LD = 0
  90. L2 = L1*IP
  91. IDO = N/L2
  92. IPM = IP-1
  93. DO 109 J=1,IPM
  94. LD = LD+L1
  95. I = IS
  96. ARGLD = LD*ARGH
  97. FI = 0.
  98. DO 108 II=3,IDO,2
  99. I = I+2
  100. FI = FI+1.
  101. ARG = FI*ARGLD
  102. WA(I-1) = COS(ARG)
  103. WA(I) = SIN(ARG)
  104. 108 CONTINUE
  105. IS = IS+IDO
  106. 109 CONTINUE
  107. L1 = L2
  108. 110 CONTINUE
  109. RETURN
  110. END