1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189
|
*DECK RADBG
SUBROUTINE RADBG (IDO, IP, L1, IDL1, CC, C1, C2, CH, CH2, WA)
C***BEGIN PROLOGUE RADBG
C***SUBSIDIARY
C***PURPOSE Calculate the fast Fourier transform of subvectors of
C arbitrary length.
C***LIBRARY SLATEC (FFTPACK)
C***TYPE SINGLE PRECISION (RADBG-S)
C***AUTHOR Swarztrauber, P. N., (NCAR)
C***ROUTINES CALLED (NONE)
C***REVISION HISTORY (YYMMDD)
C 790601 DATE WRITTEN
C 830401 Modified to use SLATEC library source file format.
C 860115 Modified by Ron Boisvert to adhere to Fortran 77 by
C (a) changing dummy array size declarations (1) to (*),
C (b) changing references to intrinsic function FLOAT
C to REAL, and
C (c) changing definition of variable TPI by using
C FORTRAN intrinsic function ATAN instead of a DATA
C statement.
C 881128 Modified by Dick Valent to meet prologue standards.
C 890531 Changed all specific intrinsics to generic. (WRB)
C 890831 Modified array declarations. (WRB)
C 891214 Prologue converted to Version 4.0 format. (BAB)
C 900402 Added TYPE section. (WRB)
C***END PROLOGUE RADBG
DIMENSION CH(IDO,L1,*), CC(IDO,IP,*), C1(IDO,L1,*),
+ C2(IDL1,*), CH2(IDL1,*), WA(*)
C***FIRST EXECUTABLE STATEMENT RADBG
TPI = 8.*ATAN(1.)
ARG = TPI/IP
DCP = COS(ARG)
DSP = SIN(ARG)
IDP2 = IDO+2
NBD = (IDO-1)/2
IPP2 = IP+2
IPPH = (IP+1)/2
IF (IDO .LT. L1) GO TO 103
DO 102 K=1,L1
DO 101 I=1,IDO
CH(I,K,1) = CC(I,1,K)
101 CONTINUE
102 CONTINUE
GO TO 106
103 DO 105 I=1,IDO
DO 104 K=1,L1
CH(I,K,1) = CC(I,1,K)
104 CONTINUE
105 CONTINUE
106 DO 108 J=2,IPPH
JC = IPP2-J
J2 = J+J
DO 107 K=1,L1
CH(1,K,J) = CC(IDO,J2-2,K)+CC(IDO,J2-2,K)
CH(1,K,JC) = CC(1,J2-1,K)+CC(1,J2-1,K)
107 CONTINUE
108 CONTINUE
IF (IDO .EQ. 1) GO TO 116
IF (NBD .LT. L1) GO TO 112
DO 111 J=2,IPPH
JC = IPP2-J
DO 110 K=1,L1
CDIR$ IVDEP
DO 109 I=3,IDO,2
IC = IDP2-I
CH(I-1,K,J) = CC(I-1,2*J-1,K)+CC(IC-1,2*J-2,K)
CH(I-1,K,JC) = CC(I-1,2*J-1,K)-CC(IC-1,2*J-2,K)
CH(I,K,J) = CC(I,2*J-1,K)-CC(IC,2*J-2,K)
CH(I,K,JC) = CC(I,2*J-1,K)+CC(IC,2*J-2,K)
109 CONTINUE
110 CONTINUE
111 CONTINUE
GO TO 116
112 DO 115 J=2,IPPH
JC = IPP2-J
CDIR$ IVDEP
DO 114 I=3,IDO,2
IC = IDP2-I
DO 113 K=1,L1
CH(I-1,K,J) = CC(I-1,2*J-1,K)+CC(IC-1,2*J-2,K)
CH(I-1,K,JC) = CC(I-1,2*J-1,K)-CC(IC-1,2*J-2,K)
CH(I,K,J) = CC(I,2*J-1,K)-CC(IC,2*J-2,K)
CH(I,K,JC) = CC(I,2*J-1,K)+CC(IC,2*J-2,K)
113 CONTINUE
114 CONTINUE
115 CONTINUE
116 AR1 = 1.
AI1 = 0.
DO 120 L=2,IPPH
LC = IPP2-L
AR1H = DCP*AR1-DSP*AI1
AI1 = DCP*AI1+DSP*AR1
AR1 = AR1H
DO 117 IK=1,IDL1
C2(IK,L) = CH2(IK,1)+AR1*CH2(IK,2)
C2(IK,LC) = AI1*CH2(IK,IP)
117 CONTINUE
DC2 = AR1
DS2 = AI1
AR2 = AR1
AI2 = AI1
DO 119 J=3,IPPH
JC = IPP2-J
AR2H = DC2*AR2-DS2*AI2
AI2 = DC2*AI2+DS2*AR2
AR2 = AR2H
DO 118 IK=1,IDL1
C2(IK,L) = C2(IK,L)+AR2*CH2(IK,J)
C2(IK,LC) = C2(IK,LC)+AI2*CH2(IK,JC)
118 CONTINUE
119 CONTINUE
120 CONTINUE
DO 122 J=2,IPPH
DO 121 IK=1,IDL1
CH2(IK,1) = CH2(IK,1)+CH2(IK,J)
121 CONTINUE
122 CONTINUE
DO 124 J=2,IPPH
JC = IPP2-J
DO 123 K=1,L1
CH(1,K,J) = C1(1,K,J)-C1(1,K,JC)
CH(1,K,JC) = C1(1,K,J)+C1(1,K,JC)
123 CONTINUE
124 CONTINUE
IF (IDO .EQ. 1) GO TO 132
IF (NBD .LT. L1) GO TO 128
DO 127 J=2,IPPH
JC = IPP2-J
DO 126 K=1,L1
CDIR$ IVDEP
DO 125 I=3,IDO,2
CH(I-1,K,J) = C1(I-1,K,J)-C1(I,K,JC)
CH(I-1,K,JC) = C1(I-1,K,J)+C1(I,K,JC)
CH(I,K,J) = C1(I,K,J)+C1(I-1,K,JC)
CH(I,K,JC) = C1(I,K,J)-C1(I-1,K,JC)
125 CONTINUE
126 CONTINUE
127 CONTINUE
GO TO 132
128 DO 131 J=2,IPPH
JC = IPP2-J
DO 130 I=3,IDO,2
DO 129 K=1,L1
CH(I-1,K,J) = C1(I-1,K,J)-C1(I,K,JC)
CH(I-1,K,JC) = C1(I-1,K,J)+C1(I,K,JC)
CH(I,K,J) = C1(I,K,J)+C1(I-1,K,JC)
CH(I,K,JC) = C1(I,K,J)-C1(I-1,K,JC)
129 CONTINUE
130 CONTINUE
131 CONTINUE
132 CONTINUE
IF (IDO .EQ. 1) RETURN
DO 133 IK=1,IDL1
C2(IK,1) = CH2(IK,1)
133 CONTINUE
DO 135 J=2,IP
DO 134 K=1,L1
C1(1,K,J) = CH(1,K,J)
134 CONTINUE
135 CONTINUE
IF (NBD .GT. L1) GO TO 139
IS = -IDO
DO 138 J=2,IP
IS = IS+IDO
IDIJ = IS
DO 137 I=3,IDO,2
IDIJ = IDIJ+2
DO 136 K=1,L1
C1(I-1,K,J) = WA(IDIJ-1)*CH(I-1,K,J)-WA(IDIJ)*CH(I,K,J)
C1(I,K,J) = WA(IDIJ-1)*CH(I,K,J)+WA(IDIJ)*CH(I-1,K,J)
136 CONTINUE
137 CONTINUE
138 CONTINUE
GO TO 143
139 IS = -IDO
DO 142 J=2,IP
IS = IS+IDO
DO 141 K=1,L1
IDIJ = IS
CDIR$ IVDEP
DO 140 I=3,IDO,2
IDIJ = IDIJ+2
C1(I-1,K,J) = WA(IDIJ-1)*CH(I-1,K,J)-WA(IDIJ)*CH(I,K,J)
C1(I,K,J) = WA(IDIJ-1)*CH(I,K,J)+WA(IDIJ)*CH(I-1,K,J)
140 CONTINUE
141 CONTINUE
142 CONTINUE
143 RETURN
END
|