1 | !CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC |
---|
2 | ! |
---|
3 | ! FFTPACK 5.0 |
---|
4 | ! |
---|
5 | ! Authors: Paul N. Swarztrauber and Richard A. Valent |
---|
6 | ! |
---|
7 | ! $Id: mradbg.f,v 1.2 2004/06/15 21:29:20 rodney Exp $ |
---|
8 | ! |
---|
9 | !CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC |
---|
10 | |
---|
11 | SUBROUTINE MRADBG (M,IDO,IP,L1,IDL1,CC,C1,C2,IM1,IN1, & |
---|
12 | & CH,CH2,IM2,IN2,WA) |
---|
13 | REAL CH(IN2,IDO,L1,IP) ,CC(IN1,IDO,IP,L1) , & |
---|
14 | & C1(IN1,IDO,L1,IP) ,C2(IN1,IDL1,IP), & |
---|
15 | & CH2(IN2,IDL1,IP) ,WA(IDO) |
---|
16 | ! |
---|
17 | M1D = (M-1)*IM1+1 |
---|
18 | M2S = 1-IM2 |
---|
19 | TPI=2.*4.*ATAN(1.0) |
---|
20 | ARG = TPI/FLOAT(IP) |
---|
21 | DCP = COS(ARG) |
---|
22 | DSP = SIN(ARG) |
---|
23 | IDP2 = IDO+2 |
---|
24 | NBD = (IDO-1)/2 |
---|
25 | IPP2 = IP+2 |
---|
26 | IPPH = (IP+1)/2 |
---|
27 | IF (IDO .LT. L1) GO TO 103 |
---|
28 | DO 102 K=1,L1 |
---|
29 | DO 101 I=1,IDO |
---|
30 | M2 = M2S |
---|
31 | DO 1001 M1=1,M1D,IM1 |
---|
32 | M2 = M2+IM2 |
---|
33 | CH(M2,I,K,1) = CC(M1,I,1,K) |
---|
34 | 1001 CONTINUE |
---|
35 | 101 CONTINUE |
---|
36 | 102 END DO |
---|
37 | GO TO 106 |
---|
38 | 103 DO 105 I=1,IDO |
---|
39 | DO 104 K=1,L1 |
---|
40 | M2 = M2S |
---|
41 | DO 1004 M1=1,M1D,IM1 |
---|
42 | M2 = M2+IM2 |
---|
43 | CH(M2,I,K,1) = CC(M1,I,1,K) |
---|
44 | 1004 CONTINUE |
---|
45 | 104 CONTINUE |
---|
46 | 105 END DO |
---|
47 | 106 DO 108 J=2,IPPH |
---|
48 | JC = IPP2-J |
---|
49 | J2 = J+J |
---|
50 | DO 107 K=1,L1 |
---|
51 | M2 = M2S |
---|
52 | DO 1007 M1=1,M1D,IM1 |
---|
53 | M2 = M2+IM2 |
---|
54 | CH(M2,1,K,J) = CC(M1,IDO,J2-2,K)+CC(M1,IDO,J2-2,K) |
---|
55 | CH(M2,1,K,JC) = CC(M1,1,J2-1,K)+CC(M1,1,J2-1,K) |
---|
56 | 1007 CONTINUE |
---|
57 | 107 CONTINUE |
---|
58 | 108 END DO |
---|
59 | IF (IDO .EQ. 1) GO TO 116 |
---|
60 | IF (NBD .LT. L1) GO TO 112 |
---|
61 | DO 111 J=2,IPPH |
---|
62 | JC = IPP2-J |
---|
63 | DO 110 K=1,L1 |
---|
64 | DO 109 I=3,IDO,2 |
---|
65 | IC = IDP2-I |
---|
66 | M2 = M2S |
---|
67 | DO 1009 M1=1,M1D,IM1 |
---|
68 | M2 = M2+IM2 |
---|
69 | CH(M2,I-1,K,J) = CC(M1,I-1,2*J-1,K)+CC(M1,IC-1,2*J-2,K) |
---|
70 | CH(M2,I-1,K,JC) = CC(M1,I-1,2*J-1,K)-CC(M1,IC-1,2*J-2,K) |
---|
71 | CH(M2,I,K,J) = CC(M1,I,2*J-1,K)-CC(M1,IC,2*J-2,K) |
---|
72 | CH(M2,I,K,JC) = CC(M1,I,2*J-1,K)+CC(M1,IC,2*J-2,K) |
---|
73 | 1009 CONTINUE |
---|
74 | 109 CONTINUE |
---|
75 | 110 CONTINUE |
---|
76 | 111 END DO |
---|
77 | GO TO 116 |
---|
78 | 112 DO 115 J=2,IPPH |
---|
79 | JC = IPP2-J |
---|
80 | DO 114 I=3,IDO,2 |
---|
81 | IC = IDP2-I |
---|
82 | DO 113 K=1,L1 |
---|
83 | M2 = M2S |
---|
84 | DO 1013 M1=1,M1D,IM1 |
---|
85 | M2 = M2+IM2 |
---|
86 | CH(M2,I-1,K,J) = CC(M1,I-1,2*J-1,K)+CC(M1,IC-1,2*J-2,K) |
---|
87 | CH(M2,I-1,K,JC) = CC(M1,I-1,2*J-1,K)-CC(M1,IC-1,2*J-2,K) |
---|
88 | CH(M2,I,K,J) = CC(M1,I,2*J-1,K)-CC(M1,IC,2*J-2,K) |
---|
89 | CH(M2,I,K,JC) = CC(M1,I,2*J-1,K)+CC(M1,IC,2*J-2,K) |
---|
90 | 1013 CONTINUE |
---|
91 | 113 CONTINUE |
---|
92 | 114 CONTINUE |
---|
93 | 115 END DO |
---|
94 | 116 AR1 = 1. |
---|
95 | AI1 = 0. |
---|
96 | DO 120 L=2,IPPH |
---|
97 | LC = IPP2-L |
---|
98 | AR1H = DCP*AR1-DSP*AI1 |
---|
99 | AI1 = DCP*AI1+DSP*AR1 |
---|
100 | AR1 = AR1H |
---|
101 | DO 117 IK=1,IDL1 |
---|
102 | M2 = M2S |
---|
103 | DO 1017 M1=1,M1D,IM1 |
---|
104 | M2 = M2+IM2 |
---|
105 | C2(M1,IK,L) = CH2(M2,IK,1)+AR1*CH2(M2,IK,2) |
---|
106 | C2(M1,IK,LC) = AI1*CH2(M2,IK,IP) |
---|
107 | 1017 CONTINUE |
---|
108 | 117 CONTINUE |
---|
109 | DC2 = AR1 |
---|
110 | DS2 = AI1 |
---|
111 | AR2 = AR1 |
---|
112 | AI2 = AI1 |
---|
113 | DO 119 J=3,IPPH |
---|
114 | JC = IPP2-J |
---|
115 | AR2H = DC2*AR2-DS2*AI2 |
---|
116 | AI2 = DC2*AI2+DS2*AR2 |
---|
117 | AR2 = AR2H |
---|
118 | DO 118 IK=1,IDL1 |
---|
119 | M2 = M2S |
---|
120 | DO 1018 M1=1,M1D,IM1 |
---|
121 | M2 = M2+IM2 |
---|
122 | C2(M1,IK,L) = C2(M1,IK,L)+AR2*CH2(M2,IK,J) |
---|
123 | C2(M1,IK,LC) = C2(M1,IK,LC)+AI2*CH2(M2,IK,JC) |
---|
124 | 1018 CONTINUE |
---|
125 | 118 CONTINUE |
---|
126 | 119 CONTINUE |
---|
127 | 120 END DO |
---|
128 | DO 122 J=2,IPPH |
---|
129 | DO 121 IK=1,IDL1 |
---|
130 | M2 = M2S |
---|
131 | DO 1021 M1=1,M1D,IM1 |
---|
132 | M2 = M2+IM2 |
---|
133 | CH2(M2,IK,1) = CH2(M2,IK,1)+CH2(M2,IK,J) |
---|
134 | 1021 CONTINUE |
---|
135 | 121 CONTINUE |
---|
136 | 122 END DO |
---|
137 | DO 124 J=2,IPPH |
---|
138 | JC = IPP2-J |
---|
139 | DO 123 K=1,L1 |
---|
140 | M2 = M2S |
---|
141 | DO 1023 M1=1,M1D,IM1 |
---|
142 | M2 = M2+IM2 |
---|
143 | CH(M2,1,K,J) = C1(M1,1,K,J)-C1(M1,1,K,JC) |
---|
144 | CH(M2,1,K,JC) = C1(M1,1,K,J)+C1(M1,1,K,JC) |
---|
145 | 1023 CONTINUE |
---|
146 | 123 CONTINUE |
---|
147 | 124 END DO |
---|
148 | IF (IDO .EQ. 1) GO TO 132 |
---|
149 | IF (NBD .LT. L1) GO TO 128 |
---|
150 | DO 127 J=2,IPPH |
---|
151 | JC = IPP2-J |
---|
152 | DO 126 K=1,L1 |
---|
153 | DO 125 I=3,IDO,2 |
---|
154 | M2 = M2S |
---|
155 | DO 1025 M1=1,M1D,IM1 |
---|
156 | M2 = M2+IM2 |
---|
157 | CH(M2,I-1,K,J) = C1(M1,I-1,K,J)-C1(M1,I,K,JC) |
---|
158 | CH(M2,I-1,K,JC) = C1(M1,I-1,K,J)+C1(M1,I,K,JC) |
---|
159 | CH(M2,I,K,J) = C1(M1,I,K,J)+C1(M1,I-1,K,JC) |
---|
160 | CH(M2,I,K,JC) = C1(M1,I,K,J)-C1(M1,I-1,K,JC) |
---|
161 | 1025 CONTINUE |
---|
162 | 125 CONTINUE |
---|
163 | 126 CONTINUE |
---|
164 | 127 END DO |
---|
165 | GO TO 132 |
---|
166 | 128 DO 131 J=2,IPPH |
---|
167 | JC = IPP2-J |
---|
168 | DO 130 I=3,IDO,2 |
---|
169 | DO 129 K=1,L1 |
---|
170 | M2 = M2S |
---|
171 | DO 1029 M1=1,M1D,IM1 |
---|
172 | M2 = M2+IM2 |
---|
173 | CH(M2,I-1,K,J) = C1(M1,I-1,K,J)-C1(M1,I,K,JC) |
---|
174 | CH(M2,I-1,K,JC) = C1(M1,I-1,K,J)+C1(M1,I,K,JC) |
---|
175 | CH(M2,I,K,J) = C1(M1,I,K,J)+C1(M1,I-1,K,JC) |
---|
176 | CH(M2,I,K,JC) = C1(M1,I,K,J)-C1(M1,I-1,K,JC) |
---|
177 | 1029 CONTINUE |
---|
178 | 129 CONTINUE |
---|
179 | 130 CONTINUE |
---|
180 | 131 END DO |
---|
181 | 132 CONTINUE |
---|
182 | IF (IDO .EQ. 1) RETURN |
---|
183 | DO 133 IK=1,IDL1 |
---|
184 | M2 = M2S |
---|
185 | DO 1033 M1=1,M1D,IM1 |
---|
186 | M2 = M2+IM2 |
---|
187 | C2(M1,IK,1) = CH2(M2,IK,1) |
---|
188 | 1033 CONTINUE |
---|
189 | 133 END DO |
---|
190 | DO 135 J=2,IP |
---|
191 | DO 134 K=1,L1 |
---|
192 | M2 = M2S |
---|
193 | DO 1034 M1=1,M1D,IM1 |
---|
194 | M2 = M2+IM2 |
---|
195 | C1(M1,1,K,J) = CH(M2,1,K,J) |
---|
196 | 1034 CONTINUE |
---|
197 | 134 CONTINUE |
---|
198 | 135 END DO |
---|
199 | IF (NBD .GT. L1) GO TO 139 |
---|
200 | IS = -IDO |
---|
201 | DO 138 J=2,IP |
---|
202 | IS = IS+IDO |
---|
203 | IDIJ = IS |
---|
204 | DO 137 I=3,IDO,2 |
---|
205 | IDIJ = IDIJ+2 |
---|
206 | DO 136 K=1,L1 |
---|
207 | M2 = M2S |
---|
208 | DO 1036 M1=1,M1D,IM1 |
---|
209 | M2 = M2+IM2 |
---|
210 | C1(M1,I-1,K,J) = WA(IDIJ-1)*CH(M2,I-1,K,J)-WA(IDIJ)* & |
---|
211 | & CH(M2,I,K,J) |
---|
212 | C1(M1,I,K,J) = WA(IDIJ-1)*CH(M2,I,K,J)+WA(IDIJ)* & |
---|
213 | & CH(M2,I-1,K,J) |
---|
214 | 1036 CONTINUE |
---|
215 | 136 CONTINUE |
---|
216 | 137 CONTINUE |
---|
217 | 138 END DO |
---|
218 | GO TO 143 |
---|
219 | 139 IS = -IDO |
---|
220 | DO 142 J=2,IP |
---|
221 | IS = IS+IDO |
---|
222 | DO 141 K=1,L1 |
---|
223 | IDIJ = IS |
---|
224 | DO 140 I=3,IDO,2 |
---|
225 | IDIJ = IDIJ+2 |
---|
226 | M2 = M2S |
---|
227 | DO 1040 M1=1,M1D,IM1 |
---|
228 | M2 = M2+IM2 |
---|
229 | C1(M1,I-1,K,J) = WA(IDIJ-1)*CH(M2,I-1,K,J)-WA(IDIJ)* & |
---|
230 | & CH(M2,I,K,J) |
---|
231 | C1(M1,I,K,J) = WA(IDIJ-1)*CH(M2,I,K,J)+WA(IDIJ)* & |
---|
232 | & CH(M2,I-1,K,J) |
---|
233 | 1040 CONTINUE |
---|
234 | 140 CONTINUE |
---|
235 | 141 CONTINUE |
---|
236 | 142 END DO |
---|
237 | 143 RETURN |
---|
238 | END |
---|