[4f24798] | 1 | !c Copyright 1981-2016 ECMWF. |
---|
| 2 | !c |
---|
| 3 | !c This software is licensed under the terms of the Apache Licence |
---|
| 4 | !c Version 2.0 which can be obtained at http://www.apache.org/licenses/LICENSE-2.0. |
---|
| 5 | !c |
---|
| 6 | !c In applying this licence, ECMWF does not waive the privileges and immunities |
---|
| 7 | !c granted to it by virtue of its status as an intergovernmental organisation |
---|
| 8 | !c nor does it submit to any jurisdiction. |
---|
| 9 | !c |
---|
| 10 | |
---|
| 11 | SUBROUTINE QPASSM(A,B,C,D,TRIGS,INC1,INC2,INC3,INC4,ILOT,N,IFAC,ILA,IERR) |
---|
| 12 | REAL A(*),B(*),C(*),D(*),TRIGS(*) |
---|
| 13 | !C |
---|
| 14 | !C SUBROUTINE 'QPASSM' - PERFORMS ONE PASS THROUGH DATA AS PART |
---|
| 15 | !C OF MULTIPLE REAL FFT (FOURIER ANALYSIS) ROUTINE |
---|
| 16 | !C |
---|
| 17 | !C A IS FIRST REAL INPUT VECTOR |
---|
| 18 | !C EQUIVALENCE B(1) WITH A(IFAC*ILA*INC1+1) |
---|
| 19 | !C C IS FIRST REAL OUTPUT VECTOR |
---|
| 20 | !C EQUIVALENCE D(1) WITH C(ILA*INC2+1) |
---|
| 21 | !C TRIGS IS A PRECALCULATED LIST OF SINES & COSINES |
---|
| 22 | !C INC1 IS THE ADDRESSING INCREMENT FOR A |
---|
| 23 | !C INC2 IS THE ADDRESSING INCREMENT FOR C |
---|
| 24 | !C INC3 IS THE INCREMENT BETWEEN INPUT VECTORS A |
---|
| 25 | !C INC4 IS THE INCREMENT BETWEEN OUTPUT VECTORS C |
---|
| 26 | !C ILOT IS THE NUMBER OF VECTORS |
---|
| 27 | !C N IS THE LENGTH OF THE VECTORS |
---|
| 28 | !C IFAC IS THE CURRENT FACTOR OF N |
---|
| 29 | !C ILA = N/(PRODUCT OF FACTORS USED SO FAR) |
---|
| 30 | !C IERR IS AN ERROR INDICATOR: |
---|
| 31 | !C 0 - PASS COMPLETED WITHOUT ERROR |
---|
| 32 | !C 1 - ILOT GREATER THAN 64 |
---|
| 33 | !C 2 - IFAC NOT CATERED FOR |
---|
| 34 | !C 3 - IFAC ONLY CATERED FOR IF ILA=N/IFAC |
---|
| 35 | !C |
---|
| 36 | !C----------------------------------------------------------------------- |
---|
| 37 | !C |
---|
| 38 | SAVE SIN36, SIN72, QRT5, SIN60 |
---|
| 39 | |
---|
| 40 | DATA SIN36/0.587785252292473/,SIN72/0.951056516295154/,& |
---|
| 41 | QRT5/0.559016994374947/,SIN60/0.866025403784437/ |
---|
| 42 | |
---|
| 43 | M=N/IFAC |
---|
| 44 | IINK=ILA*INC1 |
---|
| 45 | JINK=ILA*INC2 |
---|
| 46 | IJUMP=(IFAC-1)*IINK |
---|
| 47 | KSTOP=(N-IFAC)/(2*IFAC) |
---|
| 48 | |
---|
| 49 | IBAD=1 |
---|
| 50 | IF (ILOT.GT.512) GO TO 910 |
---|
| 51 | IBASE=0 |
---|
| 52 | JBASE=0 |
---|
| 53 | IGO=IFAC-1 |
---|
| 54 | IF (IGO.EQ.7) IGO=6 |
---|
| 55 | IBAD=2 |
---|
| 56 | IF (IGO.LT.1.OR.IGO.GT.6) GO TO 910 |
---|
| 57 | GO TO (200,300,400,500,600,800),IGO |
---|
| 58 | |
---|
| 59 | ! CODING FOR FACTOR 2 |
---|
| 60 | ! ------------------- |
---|
| 61 | 200 CONTINUE |
---|
| 62 | IA=1 |
---|
| 63 | IB=IA+IINK |
---|
| 64 | JA=1 |
---|
| 65 | JB=JA+(2*M-ILA)*INC2 |
---|
| 66 | |
---|
| 67 | IF (ILA.EQ.M) GO TO 290 |
---|
| 68 | |
---|
| 69 | DO 220 JL=1,ILA |
---|
| 70 | I=IBASE |
---|
| 71 | J=JBASE |
---|
| 72 | !DIR$ IVDEP |
---|
| 73 | !VOCL LOOP,NOVREC |
---|
| 74 | DO 210 IJK=1,ILOT |
---|
| 75 | C(JA+J)=A(IA+I)+A(IB+I) |
---|
| 76 | C(JB+J)=A(IA+I)-A(IB+I) |
---|
| 77 | I=I+INC3 |
---|
| 78 | J=J+INC4 |
---|
| 79 | 210 CONTINUE |
---|
| 80 | IBASE=IBASE+INC1 |
---|
| 81 | JBASE=JBASE+INC2 |
---|
| 82 | 220 CONTINUE |
---|
| 83 | JA=JA+JINK |
---|
| 84 | JINK=2*JINK |
---|
| 85 | JB=JB-JINK |
---|
| 86 | IBASE=IBASE+IJUMP |
---|
| 87 | IJUMP=2*IJUMP+IINK |
---|
| 88 | IF (JA.EQ.JB) GO TO 260 |
---|
| 89 | DO 250 K=ILA,KSTOP,ILA |
---|
| 90 | KB=K+K |
---|
| 91 | C1=TRIGS(KB+1) |
---|
| 92 | S1=TRIGS(KB+2) |
---|
| 93 | JBASE=0 |
---|
| 94 | DO 240 JL=1,ILA |
---|
| 95 | I=IBASE |
---|
| 96 | J=JBASE |
---|
| 97 | !DIR$ IVDEP |
---|
| 98 | !VOCL LOOP,NOVREC |
---|
| 99 | DO 230 IJK=1,ILOT |
---|
| 100 | C(JA+J)=A(IA+I)+(C1*A(IB+I)+S1*B(IB+I)) |
---|
| 101 | C(JB+J)=A(IA+I)-(C1*A(IB+I)+S1*B(IB+I)) |
---|
| 102 | D(JA+J)=(C1*B(IB+I)-S1*A(IB+I))+B(IA+I) |
---|
| 103 | D(JB+J)=(C1*B(IB+I)-S1*A(IB+I))-B(IA+I) |
---|
| 104 | I=I+INC3 |
---|
| 105 | J=J+INC4 |
---|
| 106 | 230 CONTINUE |
---|
| 107 | IBASE=IBASE+INC1 |
---|
| 108 | JBASE=JBASE+INC2 |
---|
| 109 | 240 CONTINUE |
---|
| 110 | IBASE=IBASE+IJUMP |
---|
| 111 | JA=JA+JINK |
---|
| 112 | JB=JB-JINK |
---|
| 113 | 250 CONTINUE |
---|
| 114 | IF (JA.GT.JB) GO TO 900 |
---|
| 115 | 260 CONTINUE |
---|
| 116 | JBASE=0 |
---|
| 117 | DO 280 JL=1,ILA |
---|
| 118 | I=IBASE |
---|
| 119 | J=JBASE |
---|
| 120 | !DIR$ IVDEP |
---|
| 121 | !VOCL LOOP,NOVREC |
---|
| 122 | DO 270 IJK=1,ILOT |
---|
| 123 | C(JA+J)=A(IA+I) |
---|
| 124 | D(JA+J)=-A(IB+I) |
---|
| 125 | I=I+INC3 |
---|
| 126 | J=J+INC4 |
---|
| 127 | 270 CONTINUE |
---|
| 128 | IBASE=IBASE+INC1 |
---|
| 129 | JBASE=JBASE+INC2 |
---|
| 130 | 280 CONTINUE |
---|
| 131 | GO TO 900 |
---|
| 132 | |
---|
| 133 | 290 CONTINUE |
---|
| 134 | Z=1.0/FLOAT(N) |
---|
| 135 | DO 294 JL=1,ILA |
---|
| 136 | I=IBASE |
---|
| 137 | J=JBASE |
---|
| 138 | !DIR$ IVDEP |
---|
| 139 | !VOCL LOOP,NOVREC |
---|
| 140 | DO 292 IJK=1,ILOT |
---|
| 141 | C(JA+J)=Z*(A(IA+I)+A(IB+I)) |
---|
| 142 | C(JB+J)=Z*(A(IA+I)-A(IB+I)) |
---|
| 143 | I=I+INC3 |
---|
| 144 | J=J+INC4 |
---|
| 145 | 292 CONTINUE |
---|
| 146 | IBASE=IBASE+INC1 |
---|
| 147 | JBASE=JBASE+INC2 |
---|
| 148 | 294 CONTINUE |
---|
| 149 | GO TO 900 |
---|
| 150 | |
---|
| 151 | ! CODING FOR FACTOR 3 |
---|
| 152 | ! ------------------- |
---|
| 153 | 300 CONTINUE |
---|
| 154 | IA=1 |
---|
| 155 | IB=IA+IINK |
---|
| 156 | IC=IB+IINK |
---|
| 157 | JA=1 |
---|
| 158 | JB=JA+(2*M-ILA)*INC2 |
---|
| 159 | JC=JB |
---|
| 160 | |
---|
| 161 | IF (ILA.EQ.M) GO TO 390 |
---|
| 162 | |
---|
| 163 | DO 320 JL=1,ILA |
---|
| 164 | I=IBASE |
---|
| 165 | J=JBASE |
---|
| 166 | !DIR$ IVDEP |
---|
| 167 | !VOCL LOOP,NOVREC |
---|
| 168 | DO 310 IJK=1,ILOT |
---|
| 169 | C(JA+J)=A(IA+I)+(A(IB+I)+A(IC+I)) |
---|
| 170 | C(JB+J)=A(IA+I)-0.5*(A(IB+I)+A(IC+I)) |
---|
| 171 | D(JB+J)=SIN60*(A(IC+I)-A(IB+I)) |
---|
| 172 | I=I+INC3 |
---|
| 173 | J=J+INC4 |
---|
| 174 | 310 CONTINUE |
---|
| 175 | IBASE=IBASE+INC1 |
---|
| 176 | JBASE=JBASE+INC2 |
---|
| 177 | 320 CONTINUE |
---|
| 178 | JA=JA+JINK |
---|
| 179 | JINK=2*JINK |
---|
| 180 | JB=JB+JINK |
---|
| 181 | JC=JC-JINK |
---|
| 182 | IBASE=IBASE+IJUMP |
---|
| 183 | IJUMP=2*IJUMP+IINK |
---|
| 184 | IF (JA.EQ.JC) GO TO 360 |
---|
| 185 | DO 350 K=ILA,KSTOP,ILA |
---|
| 186 | KB=K+K |
---|
| 187 | KC=KB+KB |
---|
| 188 | C1=TRIGS(KB+1) |
---|
| 189 | S1=TRIGS(KB+2) |
---|
| 190 | C2=TRIGS(KC+1) |
---|
| 191 | S2=TRIGS(KC+2) |
---|
| 192 | JBASE=0 |
---|
| 193 | DO 340 JL=1,ILA |
---|
| 194 | I=IBASE |
---|
| 195 | J=JBASE |
---|
| 196 | !DIR$ IVDEP |
---|
| 197 | !VOCL LOOP,NOVREC |
---|
| 198 | DO 330 IJK=1,ILOT |
---|
| 199 | A1=(C1*A(IB+I)+S1*B(IB+I))+(C2*A(IC+I)+S2*B(IC+I)) |
---|
| 200 | B1=(C1*B(IB+I)-S1*A(IB+I))+(C2*B(IC+I)-S2*A(IC+I)) |
---|
| 201 | A2=A(IA+I)-0.5*A1 |
---|
| 202 | B2=B(IA+I)-0.5*B1 |
---|
| 203 | A3=SIN60*((C1*A(IB+I)+S1*B(IB+I))-(C2*A(IC+I)+S2*B(IC+I))) |
---|
| 204 | B3=SIN60*((C1*B(IB+I)-S1*A(IB+I))-(C2*B(IC+I)-S2*A(IC+I))) |
---|
| 205 | C(JA+J)=A(IA+I)+A1 |
---|
| 206 | D(JA+J)=B(IA+I)+B1 |
---|
| 207 | C(JB+J)=A2+B3 |
---|
| 208 | D(JB+J)=B2-A3 |
---|
| 209 | C(JC+J)=A2-B3 |
---|
| 210 | D(JC+J)=-(B2+A3) |
---|
| 211 | I=I+INC3 |
---|
| 212 | J=J+INC4 |
---|
| 213 | 330 CONTINUE |
---|
| 214 | IBASE=IBASE+INC1 |
---|
| 215 | JBASE=JBASE+INC2 |
---|
| 216 | 340 CONTINUE |
---|
| 217 | IBASE=IBASE+IJUMP |
---|
| 218 | JA=JA+JINK |
---|
| 219 | JB=JB+JINK |
---|
| 220 | JC=JC-JINK |
---|
| 221 | 350 CONTINUE |
---|
| 222 | IF (JA.GT.JC) GO TO 900 |
---|
| 223 | 360 CONTINUE |
---|
| 224 | JBASE=0 |
---|
| 225 | DO 380 JL=1,ILA |
---|
| 226 | I=IBASE |
---|
| 227 | J=JBASE |
---|
| 228 | !DIR$ IVDEP |
---|
| 229 | !VOCL LOOP,NOVREC |
---|
| 230 | DO 370 IJK=1,ILOT |
---|
| 231 | C(JA+J)=A(IA+I)+0.5*(A(IB+I)-A(IC+I)) |
---|
| 232 | D(JA+J)=-SIN60*(A(IB+I)+A(IC+I)) |
---|
| 233 | C(JB+J)=A(IA+I)-(A(IB+I)-A(IC+I)) |
---|
| 234 | I=I+INC3 |
---|
| 235 | J=J+INC4 |
---|
| 236 | 370 CONTINUE |
---|
| 237 | IBASE=IBASE+INC1 |
---|
| 238 | JBASE=JBASE+INC2 |
---|
| 239 | 380 CONTINUE |
---|
| 240 | GO TO 900 |
---|
| 241 | |
---|
| 242 | 390 CONTINUE |
---|
| 243 | Z=1.0/FLOAT(N) |
---|
| 244 | ZSIN60=Z*SIN60 |
---|
| 245 | DO 394 JL=1,ILA |
---|
| 246 | I=IBASE |
---|
| 247 | J=JBASE |
---|
| 248 | !DIR$ IVDEP |
---|
| 249 | !VOCL LOOP,NOVREC |
---|
| 250 | DO 392 IJK=1,ILOT |
---|
| 251 | C(JA+J)=Z*(A(IA+I)+(A(IB+I)+A(IC+I))) |
---|
| 252 | C(JB+J)=Z*(A(IA+I)-0.5*(A(IB+I)+A(IC+I))) |
---|
| 253 | D(JB+J)=ZSIN60*(A(IC+I)-A(IB+I)) |
---|
| 254 | I=I+INC3 |
---|
| 255 | J=J+INC4 |
---|
| 256 | 392 CONTINUE |
---|
| 257 | IBASE=IBASE+INC1 |
---|
| 258 | JBASE=JBASE+INC2 |
---|
| 259 | 394 CONTINUE |
---|
| 260 | GO TO 900 |
---|
| 261 | |
---|
| 262 | ! CODING FOR FACTOR 4 |
---|
| 263 | ! ------------------- |
---|
| 264 | 400 CONTINUE |
---|
| 265 | IA=1 |
---|
| 266 | IB=IA+IINK |
---|
| 267 | IC=IB+IINK |
---|
| 268 | ID=IC+IINK |
---|
| 269 | JA=1 |
---|
| 270 | JB=JA+(2*M-ILA)*INC2 |
---|
| 271 | JC=JB+2*M*INC2 |
---|
| 272 | JD=JB |
---|
| 273 | |
---|
| 274 | IF (ILA.EQ.M) GO TO 490 |
---|
| 275 | |
---|
| 276 | DO 420 JL=1,ILA |
---|
| 277 | I=IBASE |
---|
| 278 | J=JBASE |
---|
| 279 | !DIR$ IVDEP |
---|
| 280 | !VOCL LOOP,NOVREC |
---|
| 281 | DO 410 IJK=1,ILOT |
---|
| 282 | C(JA+J)=(A(IA+I)+A(IC+I))+(A(IB+I)+A(ID+I)) |
---|
| 283 | C(JC+J)=(A(IA+I)+A(IC+I))-(A(IB+I)+A(ID+I)) |
---|
| 284 | C(JB+J)=A(IA+I)-A(IC+I) |
---|
| 285 | D(JB+J)=A(ID+I)-A(IB+I) |
---|
| 286 | I=I+INC3 |
---|
| 287 | J=J+INC4 |
---|
| 288 | 410 CONTINUE |
---|
| 289 | IBASE=IBASE+INC1 |
---|
| 290 | JBASE=JBASE+INC2 |
---|
| 291 | 420 CONTINUE |
---|
| 292 | JA=JA+JINK |
---|
| 293 | JINK=2*JINK |
---|
| 294 | JB=JB+JINK |
---|
| 295 | JC=JC-JINK |
---|
| 296 | JD=JD-JINK |
---|
| 297 | IBASE=IBASE+IJUMP |
---|
| 298 | IJUMP=2*IJUMP+IINK |
---|
| 299 | IF (JB.EQ.JC) GO TO 460 |
---|
| 300 | DO 450 K=ILA,KSTOP,ILA |
---|
| 301 | KB=K+K |
---|
| 302 | KC=KB+KB |
---|
| 303 | KD=KC+KB |
---|
| 304 | C1=TRIGS(KB+1) |
---|
| 305 | S1=TRIGS(KB+2) |
---|
| 306 | C2=TRIGS(KC+1) |
---|
| 307 | S2=TRIGS(KC+2) |
---|
| 308 | C3=TRIGS(KD+1) |
---|
| 309 | S3=TRIGS(KD+2) |
---|
| 310 | JBASE=0 |
---|
| 311 | DO 440 JL=1,ILA |
---|
| 312 | I=IBASE |
---|
| 313 | J=JBASE |
---|
| 314 | !DIR$ IVDEP |
---|
| 315 | !VOCL LOOP,NOVREC |
---|
| 316 | DO 430 IJK=1,ILOT |
---|
| 317 | A0=A(IA+I)+(C2*A(IC+I)+S2*B(IC+I)) |
---|
| 318 | A2=A(IA+I)-(C2*A(IC+I)+S2*B(IC+I)) |
---|
| 319 | A1=(C1*A(IB+I)+S1*B(IB+I))+(C3*A(ID+I)+S3*B(ID+I)) |
---|
| 320 | A3=(C1*A(IB+I)+S1*B(IB+I))-(C3*A(ID+I)+S3*B(ID+I)) |
---|
| 321 | B0=B(IA+I)+(C2*B(IC+I)-S2*A(IC+I)) |
---|
| 322 | B2=B(IA+I)-(C2*B(IC+I)-S2*A(IC+I)) |
---|
| 323 | B1=(C1*B(IB+I)-S1*A(IB+I))+(C3*B(ID+I)-S3*A(ID+I)) |
---|
| 324 | B3=(C1*B(IB+I)-S1*A(IB+I))-(C3*B(ID+I)-S3*A(ID+I)) |
---|
| 325 | C(JA+J)=A0+A1 |
---|
| 326 | C(JC+J)=A0-A1 |
---|
| 327 | D(JA+J)=B0+B1 |
---|
| 328 | D(JC+J)=B1-B0 |
---|
| 329 | C(JB+J)=A2+B3 |
---|
| 330 | C(JD+J)=A2-B3 |
---|
| 331 | D(JB+J)=B2-A3 |
---|
| 332 | D(JD+J)=-(B2+A3) |
---|
| 333 | I=I+INC3 |
---|
| 334 | J=J+INC4 |
---|
| 335 | 430 CONTINUE |
---|
| 336 | IBASE=IBASE+INC1 |
---|
| 337 | JBASE=JBASE+INC2 |
---|
| 338 | 440 CONTINUE |
---|
| 339 | IBASE=IBASE+IJUMP |
---|
| 340 | JA=JA+JINK |
---|
| 341 | JB=JB+JINK |
---|
| 342 | JC=JC-JINK |
---|
| 343 | JD=JD-JINK |
---|
| 344 | 450 CONTINUE |
---|
| 345 | IF (JB.GT.JC) GO TO 900 |
---|
| 346 | 460 CONTINUE |
---|
| 347 | SIN45=SQRT(0.5) |
---|
| 348 | JBASE=0 |
---|
| 349 | DO 480 JL=1,ILA |
---|
| 350 | I=IBASE |
---|
| 351 | J=JBASE |
---|
| 352 | !DIR$ IVDEP |
---|
| 353 | !VOCL LOOP,NOVREC |
---|
| 354 | DO 470 IJK=1,ILOT |
---|
| 355 | C(JA+J)=A(IA+I)+SIN45*(A(IB+I)-A(ID+I)) |
---|
| 356 | C(JB+J)=A(IA+I)-SIN45*(A(IB+I)-A(ID+I)) |
---|
| 357 | D(JA+J)=-A(IC+I)-SIN45*(A(IB+I)+A(ID+I)) |
---|
| 358 | D(JB+J)=A(IC+I)-SIN45*(A(IB+I)+A(ID+I)) |
---|
| 359 | I=I+INC3 |
---|
| 360 | J=J+INC4 |
---|
| 361 | 470 CONTINUE |
---|
| 362 | IBASE=IBASE+INC1 |
---|
| 363 | JBASE=JBASE+INC2 |
---|
| 364 | 480 CONTINUE |
---|
| 365 | GO TO 900 |
---|
| 366 | |
---|
| 367 | 490 CONTINUE |
---|
| 368 | Z=1.0/FLOAT(N) |
---|
| 369 | DO 494 JL=1,ILA |
---|
| 370 | I=IBASE |
---|
| 371 | J=JBASE |
---|
| 372 | !DIR$ IVDEP |
---|
| 373 | !VOCL LOOP,NOVREC |
---|
| 374 | DO 492 IJK=1,ILOT |
---|
| 375 | C(JA+J)=Z*((A(IA+I)+A(IC+I))+(A(IB+I)+A(ID+I))) |
---|
| 376 | C(JC+J)=Z*((A(IA+I)+A(IC+I))-(A(IB+I)+A(ID+I))) |
---|
| 377 | C(JB+J)=Z*(A(IA+I)-A(IC+I)) |
---|
| 378 | D(JB+J)=Z*(A(ID+I)-A(IB+I)) |
---|
| 379 | I=I+INC3 |
---|
| 380 | J=J+INC4 |
---|
| 381 | 492 CONTINUE |
---|
| 382 | IBASE=IBASE+INC1 |
---|
| 383 | JBASE=JBASE+INC2 |
---|
| 384 | 494 CONTINUE |
---|
| 385 | GO TO 900 |
---|
| 386 | |
---|
| 387 | ! CODING FOR FACTOR 5 |
---|
| 388 | ! ------------------- |
---|
| 389 | 500 CONTINUE |
---|
| 390 | IA=1 |
---|
| 391 | IB=IA+IINK |
---|
| 392 | IC=IB+IINK |
---|
| 393 | ID=IC+IINK |
---|
| 394 | IE=ID+IINK |
---|
| 395 | JA=1 |
---|
| 396 | JB=JA+(2*M-ILA)*INC2 |
---|
| 397 | JC=JB+2*M*INC2 |
---|
| 398 | JD=JC |
---|
| 399 | JE=JB |
---|
| 400 | |
---|
| 401 | IF (ILA.EQ.M) GO TO 590 |
---|
| 402 | |
---|
| 403 | DO 520 JL=1,ILA |
---|
| 404 | I=IBASE |
---|
| 405 | J=JBASE |
---|
| 406 | !DIR$ IVDEP |
---|
| 407 | !VOCL LOOP,NOVREC |
---|
| 408 | DO 510 IJK=1,ILOT |
---|
| 409 | A1=A(IB+I)+A(IE+I) |
---|
| 410 | A3=A(IB+I)-A(IE+I) |
---|
| 411 | A2=A(IC+I)+A(ID+I) |
---|
| 412 | A4=A(IC+I)-A(ID+I) |
---|
| 413 | A5=A(IA+I)-0.25*(A1+A2) |
---|
| 414 | A6=QRT5*(A1-A2) |
---|
| 415 | C(JA+J)=A(IA+I)+(A1+A2) |
---|
| 416 | C(JB+J)=A5+A6 |
---|
| 417 | C(JC+J)=A5-A6 |
---|
| 418 | D(JB+J)=-SIN72*A3-SIN36*A4 |
---|
| 419 | D(JC+J)=-SIN36*A3+SIN72*A4 |
---|
| 420 | I=I+INC3 |
---|
| 421 | J=J+INC4 |
---|
| 422 | 510 CONTINUE |
---|
| 423 | IBASE=IBASE+INC1 |
---|
| 424 | JBASE=JBASE+INC2 |
---|
| 425 | 520 CONTINUE |
---|
| 426 | JA=JA+JINK |
---|
| 427 | JINK=2*JINK |
---|
| 428 | JB=JB+JINK |
---|
| 429 | JC=JC+JINK |
---|
| 430 | JD=JD-JINK |
---|
| 431 | JE=JE-JINK |
---|
| 432 | IBASE=IBASE+IJUMP |
---|
| 433 | IJUMP=2*IJUMP+IINK |
---|
| 434 | IF (JB.EQ.JD) GO TO 560 |
---|
| 435 | DO 550 K=ILA,KSTOP,ILA |
---|
| 436 | KB=K+K |
---|
| 437 | KC=KB+KB |
---|
| 438 | KD=KC+KB |
---|
| 439 | KE=KD+KB |
---|
| 440 | C1=TRIGS(KB+1) |
---|
| 441 | S1=TRIGS(KB+2) |
---|
| 442 | C2=TRIGS(KC+1) |
---|
| 443 | S2=TRIGS(KC+2) |
---|
| 444 | C3=TRIGS(KD+1) |
---|
| 445 | S3=TRIGS(KD+2) |
---|
| 446 | C4=TRIGS(KE+1) |
---|
| 447 | S4=TRIGS(KE+2) |
---|
| 448 | JBASE=0 |
---|
| 449 | DO 540 JL=1,ILA |
---|
| 450 | I=IBASE |
---|
| 451 | J=JBASE |
---|
| 452 | !DIR$ IVDEP |
---|
| 453 | !VOCL LOOP,NOVREC |
---|
| 454 | DO 530 IJK=1,ILOT |
---|
| 455 | A1=(C1*A(IB+I)+S1*B(IB+I))+(C4*A(IE+I)+S4*B(IE+I)) |
---|
| 456 | A3=(C1*A(IB+I)+S1*B(IB+I))-(C4*A(IE+I)+S4*B(IE+I)) |
---|
| 457 | A2=(C2*A(IC+I)+S2*B(IC+I))+(C3*A(ID+I)+S3*B(ID+I)) |
---|
| 458 | A4=(C2*A(IC+I)+S2*B(IC+I))-(C3*A(ID+I)+S3*B(ID+I)) |
---|
| 459 | B1=(C1*B(IB+I)-S1*A(IB+I))+(C4*B(IE+I)-S4*A(IE+I)) |
---|
| 460 | B3=(C1*B(IB+I)-S1*A(IB+I))-(C4*B(IE+I)-S4*A(IE+I)) |
---|
| 461 | B2=(C2*B(IC+I)-S2*A(IC+I))+(C3*B(ID+I)-S3*A(ID+I)) |
---|
| 462 | B4=(C2*B(IC+I)-S2*A(IC+I))-(C3*B(ID+I)-S3*A(ID+I)) |
---|
| 463 | A5=A(IA+I)-0.25*(A1+A2) |
---|
| 464 | A6=QRT5*(A1-A2) |
---|
| 465 | B5=B(IA+I)-0.25*(B1+B2) |
---|
| 466 | B6=QRT5*(B1-B2) |
---|
| 467 | A10=A5+A6 |
---|
| 468 | A20=A5-A6 |
---|
| 469 | B10=B5+B6 |
---|
| 470 | B20=B5-B6 |
---|
| 471 | A11=SIN72*B3+SIN36*B4 |
---|
| 472 | A21=SIN36*B3-SIN72*B4 |
---|
| 473 | B11=SIN72*A3+SIN36*A4 |
---|
| 474 | B21=SIN36*A3-SIN72*A4 |
---|
| 475 | C(JA+J)=A(IA+I)+(A1+A2) |
---|
| 476 | C(JB+J)=A10+A11 |
---|
| 477 | C(JE+J)=A10-A11 |
---|
| 478 | C(JC+J)=A20+A21 |
---|
| 479 | C(JD+J)=A20-A21 |
---|
| 480 | D(JA+J)=B(IA+I)+(B1+B2) |
---|
| 481 | D(JB+J)=B10-B11 |
---|
| 482 | D(JE+J)=-(B10+B11) |
---|
| 483 | D(JC+J)=B20-B21 |
---|
| 484 | D(JD+J)=-(B20+B21) |
---|
| 485 | I=I+INC3 |
---|
| 486 | J=J+INC4 |
---|
| 487 | 530 CONTINUE |
---|
| 488 | IBASE=IBASE+INC1 |
---|
| 489 | JBASE=JBASE+INC2 |
---|
| 490 | 540 CONTINUE |
---|
| 491 | IBASE=IBASE+IJUMP |
---|
| 492 | JA=JA+JINK |
---|
| 493 | JB=JB+JINK |
---|
| 494 | JC=JC+JINK |
---|
| 495 | JD=JD-JINK |
---|
| 496 | JE=JE-JINK |
---|
| 497 | 550 CONTINUE |
---|
| 498 | IF (JB.GT.JD) GO TO 900 |
---|
| 499 | 560 CONTINUE |
---|
| 500 | JBASE=0 |
---|
| 501 | DO 580 JL=1,ILA |
---|
| 502 | I=IBASE |
---|
| 503 | J=JBASE |
---|
| 504 | !DIR$ IVDEP |
---|
| 505 | !VOCL LOOP,NOVREC |
---|
| 506 | DO 570 IJK=1,ILOT |
---|
| 507 | A1=A(IB+I)+A(IE+I) |
---|
| 508 | A3=A(IB+I)-A(IE+I) |
---|
| 509 | A2=A(IC+I)+A(ID+I) |
---|
| 510 | A4=A(IC+I)-A(ID+I) |
---|
| 511 | A5=A(IA+I)+0.25*(A3-A4) |
---|
| 512 | A6=QRT5*(A3+A4) |
---|
| 513 | C(JA+J)=A5+A6 |
---|
| 514 | C(JB+J)=A5-A6 |
---|
| 515 | C(JC+J)=A(IA+I)-(A3-A4) |
---|
| 516 | D(JA+J)=-SIN36*A1-SIN72*A2 |
---|
| 517 | D(JB+J)=-SIN72*A1+SIN36*A2 |
---|
| 518 | I=I+INC3 |
---|
| 519 | J=J+INC4 |
---|
| 520 | 570 CONTINUE |
---|
| 521 | IBASE=IBASE+INC1 |
---|
| 522 | JBASE=JBASE+INC2 |
---|
| 523 | 580 CONTINUE |
---|
| 524 | GO TO 900 |
---|
| 525 | |
---|
| 526 | 590 CONTINUE |
---|
| 527 | Z=1.0/FLOAT(N) |
---|
| 528 | ZQRT5=Z*QRT5 |
---|
| 529 | ZSIN36=Z*SIN36 |
---|
| 530 | ZSIN72=Z*SIN72 |
---|
| 531 | DO 594 JL=1,ILA |
---|
| 532 | I=IBASE |
---|
| 533 | J=JBASE |
---|
| 534 | !DIR$ IVDEP |
---|
| 535 | !VOCL LOOP,NOVREC |
---|
| 536 | DO 592 IJK=1,ILOT |
---|
| 537 | A1=A(IB+I)+A(IE+I) |
---|
| 538 | A3=A(IB+I)-A(IE+I) |
---|
| 539 | A2=A(IC+I)+A(ID+I) |
---|
| 540 | A4=A(IC+I)-A(ID+I) |
---|
| 541 | A5=Z*(A(IA+I)-0.25*(A1+A2)) |
---|
| 542 | A6=ZQRT5*(A1-A2) |
---|
| 543 | C(JA+J)=Z*(A(IA+I)+(A1+A2)) |
---|
| 544 | C(JB+J)=A5+A6 |
---|
| 545 | C(JC+J)=A5-A6 |
---|
| 546 | D(JB+J)=-ZSIN72*A3-ZSIN36*A4 |
---|
| 547 | D(JC+J)=-ZSIN36*A3+ZSIN72*A4 |
---|
| 548 | I=I+INC3 |
---|
| 549 | J=J+INC4 |
---|
| 550 | 592 CONTINUE |
---|
| 551 | IBASE=IBASE+INC1 |
---|
| 552 | JBASE=JBASE+INC2 |
---|
| 553 | 594 CONTINUE |
---|
| 554 | GO TO 900 |
---|
| 555 | |
---|
| 556 | ! CODING FOR FACTOR 6 |
---|
| 557 | ! ------------------- |
---|
| 558 | 600 CONTINUE |
---|
| 559 | IA=1 |
---|
| 560 | IB=IA+IINK |
---|
| 561 | IC=IB+IINK |
---|
| 562 | ID=IC+IINK |
---|
| 563 | IE=ID+IINK |
---|
| 564 | IF=IE+IINK |
---|
| 565 | JA=1 |
---|
| 566 | JB=JA+(2*M-ILA)*INC2 |
---|
| 567 | JC=JB+2*M*INC2 |
---|
| 568 | JD=JC+2*M*INC2 |
---|
| 569 | JE=JC |
---|
| 570 | JF=JB |
---|
| 571 | |
---|
| 572 | IF (ILA.EQ.M) GO TO 690 |
---|
| 573 | |
---|
| 574 | DO 620 JL=1,ILA |
---|
| 575 | I=IBASE |
---|
| 576 | J=JBASE |
---|
| 577 | !DIR$ IVDEP |
---|
| 578 | !VOCL LOOP,NOVREC |
---|
| 579 | DO 610 IJK=1,ILOT |
---|
| 580 | A11=(A(IC+I)+A(IF+I))+(A(IB+I)+A(IE+I)) |
---|
| 581 | C(JA+J)=(A(IA+I)+A(ID+I))+A11 |
---|
| 582 | C(JC+J)=(A(IA+I)+A(ID+I)-0.5*A11) |
---|
| 583 | D(JC+J)=SIN60*((A(IC+I)+A(IF+I))-(A(IB+I)+A(IE+I))) |
---|
| 584 | A11=(A(IC+I)-A(IF+I))+(A(IE+I)-A(IB+I)) |
---|
| 585 | C(JB+J)=(A(IA+I)-A(ID+I))-0.5*A11 |
---|
| 586 | D(JB+J)=SIN60*((A(IE+I)-A(IB+I))-(A(IC+I)-A(IF+I))) |
---|
| 587 | C(JD+J)=(A(IA+I)-A(ID+I))+A11 |
---|
| 588 | I=I+INC3 |
---|
| 589 | J=J+INC4 |
---|
| 590 | 610 CONTINUE |
---|
| 591 | IBASE=IBASE+INC1 |
---|
| 592 | JBASE=JBASE+INC2 |
---|
| 593 | 620 CONTINUE |
---|
| 594 | JA=JA+JINK |
---|
| 595 | JINK=2*JINK |
---|
| 596 | JB=JB+JINK |
---|
| 597 | JC=JC+JINK |
---|
| 598 | JD=JD-JINK |
---|
| 599 | JE=JE-JINK |
---|
| 600 | JF=JF-JINK |
---|
| 601 | IBASE=IBASE+IJUMP |
---|
| 602 | IJUMP=2*IJUMP+IINK |
---|
| 603 | IF (JC.EQ.JD) GO TO 660 |
---|
| 604 | DO 650 K=ILA,KSTOP,ILA |
---|
| 605 | KB=K+K |
---|
| 606 | KC=KB+KB |
---|
| 607 | KD=KC+KB |
---|
| 608 | KE=KD+KB |
---|
| 609 | KF=KE+KB |
---|
| 610 | C1=TRIGS(KB+1) |
---|
| 611 | S1=TRIGS(KB+2) |
---|
| 612 | C2=TRIGS(KC+1) |
---|
| 613 | S2=TRIGS(KC+2) |
---|
| 614 | C3=TRIGS(KD+1) |
---|
| 615 | S3=TRIGS(KD+2) |
---|
| 616 | C4=TRIGS(KE+1) |
---|
| 617 | S4=TRIGS(KE+2) |
---|
| 618 | C5=TRIGS(KF+1) |
---|
| 619 | S5=TRIGS(KF+2) |
---|
| 620 | JBASE=0 |
---|
| 621 | DO 640 JL=1,ILA |
---|
| 622 | I=IBASE |
---|
| 623 | J=JBASE |
---|
| 624 | !DIR$ IVDEP |
---|
| 625 | !VOCL LOOP,NOVREC |
---|
| 626 | DO 630 IJK=1,ILOT |
---|
| 627 | A1=C1*A(IB+I)+S1*B(IB+I) |
---|
| 628 | B1=C1*B(IB+I)-S1*A(IB+I) |
---|
| 629 | A2=C2*A(IC+I)+S2*B(IC+I) |
---|
| 630 | B2=C2*B(IC+I)-S2*A(IC+I) |
---|
| 631 | A3=C3*A(ID+I)+S3*B(ID+I) |
---|
| 632 | B3=C3*B(ID+I)-S3*A(ID+I) |
---|
| 633 | A4=C4*A(IE+I)+S4*B(IE+I) |
---|
| 634 | B4=C4*B(IE+I)-S4*A(IE+I) |
---|
| 635 | A5=C5*A(IF+I)+S5*B(IF+I) |
---|
| 636 | B5=C5*B(IF+I)-S5*A(IF+I) |
---|
| 637 | A11=(A2+A5)+(A1+A4) |
---|
| 638 | A20=(A(IA+I)+A3)-0.5*A11 |
---|
| 639 | A21=SIN60*((A2+A5)-(A1+A4)) |
---|
| 640 | B11=(B2+B5)+(B1+B4) |
---|
| 641 | B20=(B(IA+I)+B3)-0.5*B11 |
---|
| 642 | B21=SIN60*((B2+B5)-(B1+B4)) |
---|
| 643 | C(JA+J)=(A(IA+I)+A3)+A11 |
---|
| 644 | D(JA+J)=(B(IA+I)+B3)+B11 |
---|
| 645 | C(JC+J)=A20-B21 |
---|
| 646 | D(JC+J)=A21+B20 |
---|
| 647 | C(JE+J)=A20+B21 |
---|
| 648 | D(JE+J)=A21-B20 |
---|
| 649 | A11=(A2-A5)+(A4-A1) |
---|
| 650 | A20=(A(IA+I)-A3)-0.5*A11 |
---|
| 651 | A21=SIN60*((A4-A1)-(A2-A5)) |
---|
| 652 | B11=(B5-B2)-(B4-B1) |
---|
| 653 | B20=(B3-B(IA+I))-0.5*B11 |
---|
| 654 | B21=SIN60*((B5-B2)+(B4-B1)) |
---|
| 655 | C(JB+J)=A20-B21 |
---|
| 656 | D(JB+J)=A21-B20 |
---|
| 657 | C(JD+J)=A11+(A(IA+I)-A3) |
---|
| 658 | D(JD+J)=B11+(B3-B(IA+I)) |
---|
| 659 | C(JF+J)=A20+B21 |
---|
| 660 | D(JF+J)=A21+B20 |
---|
| 661 | I=I+INC3 |
---|
| 662 | J=J+INC4 |
---|
| 663 | 630 CONTINUE |
---|
| 664 | IBASE=IBASE+INC1 |
---|
| 665 | JBASE=JBASE+INC2 |
---|
| 666 | 640 CONTINUE |
---|
| 667 | IBASE=IBASE+IJUMP |
---|
| 668 | JA=JA+JINK |
---|
| 669 | JB=JB+JINK |
---|
| 670 | JC=JC+JINK |
---|
| 671 | JD=JD-JINK |
---|
| 672 | JE=JE-JINK |
---|
| 673 | JF=JF-JINK |
---|
| 674 | 650 CONTINUE |
---|
| 675 | IF (JC.GT.JD) GO TO 900 |
---|
| 676 | 660 CONTINUE |
---|
| 677 | JBASE=0 |
---|
| 678 | DO 680 JL=1,ILA |
---|
| 679 | I=IBASE |
---|
| 680 | J=JBASE |
---|
| 681 | !DIR$ IVDEP |
---|
| 682 | !VOCL LOOP,NOVREC |
---|
| 683 | DO 670 IJK=1,ILOT |
---|
| 684 | C(JA+J)=(A(IA+I)+0.5*(A(IC+I)-A(IE+I)))+ SIN60*(A(IB+I)-A(IF+I)) |
---|
| 685 | D(JA+J)=-(A(ID+I)+0.5*(A(IB+I)+A(IF+I)))-SIN60*(A(IC+I)+A(IE+I)) |
---|
| 686 | C(JB+J)=A(IA+I)-(A(IC+I)-A(IE+I)) |
---|
| 687 | D(JB+J)=A(ID+I)-(A(IB+I)+A(IF+I)) |
---|
| 688 | C(JC+J)=(A(IA+I)+0.5*(A(IC+I)-A(IE+I)))-SIN60*(A(IB+I)-A(IF+I)) |
---|
| 689 | D(JC+J)=-(A(ID+I)+0.5*(A(IB+I)+A(IF+I)))+SIN60*(A(IC+I)+A(IE+I)) |
---|
| 690 | I=I+INC3 |
---|
| 691 | J=J+INC4 |
---|
| 692 | 670 CONTINUE |
---|
| 693 | IBASE=IBASE+INC1 |
---|
| 694 | JBASE=JBASE+INC2 |
---|
| 695 | 680 CONTINUE |
---|
| 696 | GO TO 900 |
---|
| 697 | |
---|
| 698 | 690 CONTINUE |
---|
| 699 | Z=1.0/FLOAT(N) |
---|
| 700 | ZSIN60=Z*SIN60 |
---|
| 701 | DO 694 JL=1,ILA |
---|
| 702 | I=IBASE |
---|
| 703 | J=JBASE |
---|
| 704 | !DIR$ IVDEP |
---|
| 705 | !VOCL LOOP,NOVREC |
---|
| 706 | DO 692 IJK=1,ILOT |
---|
| 707 | A11=(A(IC+I)+A(IF+I))+(A(IB+I)+A(IE+I)) |
---|
| 708 | C(JA+J)=Z*((A(IA+I)+A(ID+I))+A11) |
---|
| 709 | C(JC+J)=Z*((A(IA+I)+A(ID+I))-0.5*A11) |
---|
| 710 | D(JC+J)=ZSIN60*((A(IC+I)+A(IF+I))-(A(IB+I)+A(IE+I))) |
---|
| 711 | A11=(A(IC+I)-A(IF+I))+(A(IE+I)-A(IB+I)) |
---|
| 712 | C(JB+J)=Z*((A(IA+I)-A(ID+I))-0.5*A11) |
---|
| 713 | D(JB+J)=ZSIN60*((A(IE+I)-A(IB+I))-(A(IC+I)-A(IF+I))) |
---|
| 714 | C(JD+J)=Z*((A(IA+I)-A(ID+I))+A11) |
---|
| 715 | I=I+INC3 |
---|
| 716 | J=J+INC4 |
---|
| 717 | 692 CONTINUE |
---|
| 718 | IBASE=IBASE+INC1 |
---|
| 719 | JBASE=JBASE+INC2 |
---|
| 720 | 694 CONTINUE |
---|
| 721 | GO TO 900 |
---|
| 722 | |
---|
| 723 | ! CODING FOR FACTOR 8 |
---|
| 724 | ! ------------------- |
---|
| 725 | 800 CONTINUE |
---|
| 726 | IBAD=3 |
---|
| 727 | IF (ILA.NE.M) GO TO 910 |
---|
| 728 | IA=1 |
---|
| 729 | IB=IA+IINK |
---|
| 730 | IC=IB+IINK |
---|
| 731 | ID=IC+IINK |
---|
| 732 | IE=ID+IINK |
---|
| 733 | IF=IE+IINK |
---|
| 734 | IG=IF+IINK |
---|
| 735 | IH=IG+IINK |
---|
| 736 | JA=1 |
---|
| 737 | JB=JA+ILA*INC2 |
---|
| 738 | JC=JB+2*M*INC2 |
---|
| 739 | JD=JC+2*M*INC2 |
---|
| 740 | JE=JD+2*M*INC2 |
---|
| 741 | Z=1.0/FLOAT(N) |
---|
| 742 | ZSIN45=Z*SQRT(0.5) |
---|
| 743 | |
---|
| 744 | DO 820 JL=1,ILA |
---|
| 745 | I=IBASE |
---|
| 746 | J=JBASE |
---|
| 747 | !DIR$ IVDEP |
---|
| 748 | !VOCL LOOP,NOVREC |
---|
| 749 | DO 810 IJK=1,ILOT |
---|
| 750 | C(JA+J)=Z*(((A(IA+I)+A(IE+I))+(A(IC+I)+A(IG+I)))+& |
---|
| 751 | ((A(ID+I)+A(IH+I))+(A(IB+I)+A(IF+I)))) |
---|
| 752 | C(JE+J)=Z*(((A(IA+I)+A(IE+I))+(A(IC+I)+A(IG+I)))-& |
---|
| 753 | ((A(ID+I)+A(IH+I))+(A(IB+I)+A(IF+I)))) |
---|
| 754 | C(JC+J)=Z*((A(IA+I)+A(IE+I))-(A(IC+I)+A(IG+I))) |
---|
| 755 | D(JC+J)=Z*((A(ID+I)+A(IH+I))-(A(IB+I)+A(IF+I))) |
---|
| 756 | C(JB+J)=Z*(A(IA+I)-A(IE+I))& |
---|
| 757 | +ZSIN45*((A(IH+I)-A(ID+I))-(A(IF+I)-A(IB+I))) |
---|
| 758 | C(JD+J)=Z*(A(IA+I)-A(IE+I))& |
---|
| 759 | -ZSIN45*((A(IH+I)-A(ID+I))-(A(IF+I)-A(IB+I))) |
---|
| 760 | D(JB+J)=ZSIN45*((A(IH+I)-A(ID+I))+(A(IF+I)-A(IB+I)))& |
---|
| 761 | +Z*(A(IG+I)-A(IC+I)) |
---|
| 762 | D(JD+J)=ZSIN45*((A(IH+I)-A(ID+I))+(A(IF+I)-A(IB+I)))& |
---|
| 763 | -Z*(A(IG+I)-A(IC+I)) |
---|
| 764 | I=I+INC3 |
---|
| 765 | J=J+INC4 |
---|
| 766 | 810 CONTINUE |
---|
| 767 | IBASE=IBASE+INC1 |
---|
| 768 | JBASE=JBASE+INC2 |
---|
| 769 | 820 CONTINUE |
---|
| 770 | |
---|
| 771 | ! RETURN |
---|
| 772 | ! ------ |
---|
| 773 | 900 CONTINUE |
---|
| 774 | IBAD=0 |
---|
| 775 | 910 CONTINUE |
---|
| 776 | IERR=IBAD |
---|
| 777 | RETURN |
---|
| 778 | END |
---|