cffti1.f 3.5 KB

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