1***************************** -*- Mode: Fortran -*- ******************** 2* 3* cfft.f -- FFT routines. 4* 5* Copyright (C) 1994-95 K. Scott Hunziker. 6* Copyright (C) 1990-94 The Boeing Company. 7* 8* See the file COPYING for license, warranty, and permission details. 9* 10************************************************************************ 11* 12* $Id: cfft.f,v 1.2 1996/10/19 19:55:37 ksh Exp $ 13* 14 SUBROUTINE CFFTI(N,WSAVE) 15CSE 16 IMPLICIT REAL*8(A-H,O-Z) 17C***BEGIN PROLOGUE CFFTI 18C***DATE WRITTEN 790601 (YYMMDD) 19C***REVISION DATE 830401 (YYMMDD) 20C***CATEGORY NO. J1A2 21C***KEYWORDS FOURIER TRANSFORM 22C***AUTHOR SWARZTRAUBER, P. N., (NCAR) 23CPS 24C***PURPOSE Initialize for CFFTF and CFFTB. 25C***DESCRIPTION 26C 27C Subroutine CFFTI initializes the array WSAVE which is used in 28C both CFFTF and CFFTB. The prime factorization of N together with 29C a tabulation of the trigonometric functions are computed and 30C stored in WSAVE. 31CPE 32CAS 33C Input Parameter 34C 35C N the length of the sequence to be transformed 36C 37C Output Parameter 38C 39C WSAVE a work array which must be dimensioned at least 4*N+15. 40C The same work array can be used for both CFFTF and CFFTB 41C as long as N remains unchanged. Different WSAVE arrays 42C are required for different values of N. The contents of 43C WSAVE must not be changed between calls of CFFTF or CFFTB. 44CAE 45C***REFERENCES (NONE) 46C***ROUTINES CALLED CFFTI1 47C***END PROLOGUE CFFTI 48 DIMENSION WSAVE(1) 49C***FIRST EXECUTABLE STATEMENT CFFTI 50 IF (N .EQ. 1) RETURN 51 IW1 = N+N+1 52 IW2 = IW1+N+N 53 CALL CFFTI1 (N,WSAVE(IW1),WSAVE(IW2)) 54 RETURN 55 END 56* 57 SUBROUTINE CFFTI1(N,WA,IFAC) 58CSE 59 IMPLICIT REAL*8(A-H,O-Z) 60C***BEGIN PROLOGUE CFFTI1 61C***REFER TO CFFTI 62C***ROUTINES CALLED (NONE) 63C***END PROLOGUE CFFTI1 64 DIMENSION WA(1) ,IFAC(3) ,NTRYH(4) 65 DATA NTRYH(1),NTRYH(2),NTRYH(3),NTRYH(4)/3,4,2,5/ 66C***FIRST EXECUTABLE STATEMENT CFFTI1 67 NL = N 68 NF = 0 69 J = 0 70 101 J = J+1 71 IF (J-4) 102,102,103 72 102 NTRY = NTRYH(J) 73 GO TO 104 74 103 NTRY = NTRY+2 75 104 NQ = NL/NTRY 76 NR = NL-NTRY*NQ 77 IF (NR) 101,105,101 78 105 NF = NF+1 79 IFAC(NF+2) = NTRY 80 NL = NQ 81 IF (NTRY .NE. 2) GO TO 107 82 IF (NF .EQ. 1) GO TO 107 83 DO 106 I=2,NF 84 IB = NF-I+2 85 IFAC(IB+2) = IFAC(IB+1) 86 106 CONTINUE 87 IFAC(3) = 2 88 107 IF (NL .NE. 1) GO TO 104 89 IFAC(1) = N 90 IFAC(2) = NF 91 TPI = 6.28318530717959 92 ARGH = TPI/DBLE(N) 93 I = 2 94 L1 = 1 95 DO 110 K1=1,NF 96 IP = IFAC(K1+2) 97 LD = 0 98 L2 = L1*IP 99 IDO = N/L2 100 IDOT = IDO+IDO+2 101 IPM = IP-1 102 DO 109 J=1,IPM 103 I1 = I 104 WA(I-1) = 1. 105 WA(I) = 0. 106 LD = LD+L1 107 FI = 0. 108 ARGLD = DBLE(LD)*ARGH 109 DO 108 II=4,IDOT,2 110 I = I+2 111 FI = FI+1. 112 ARG = FI*ARGLD 113 WA(I-1) = COS(ARG) 114 WA(I) = SIN(ARG) 115 108 CONTINUE 116 IF (IP .LE. 5) GO TO 109 117 WA(I1-1) = WA(I-1) 118 WA(I1) = WA(I) 119 109 CONTINUE 120 L1 = L2 121 110 CONTINUE 122 RETURN 123 END 124* 125 SUBROUTINE CFFTF(N,C,WSAVE) 126CSE 127 IMPLICIT REAL*8(A-H,O-Z) 128C***BEGIN PROLOGUE CFFTF 129C***DATE WRITTEN 790601 (YYMMDD) 130C***REVISION DATE 800626 (YYMMDD) 131C***CATEGORY NO. J1A2 132C***KEYWORDS FOURIER TRANSFORM 133C***AUTHOR SWARZTRAUBER, P. N., (NCAR) 134CPS 135C***PURPOSE Forward transform of a complex, periodic sequence. 136C***DESCRIPTION 137C 138C Subroutine CFFTF computes the forward complex discrete Fourier 139C transform (the Fourier analysis). Equivalently, CFFTF computes 140C the Fourier coefficients of a complex periodic sequence. 141C The transform is defined below at output parameter C. 142C 143C The transform is not normalized. To obtain a normalized transform 144C the output must be divided by N. Otherwise a call of CFFTF 145C followed by a call of CFFTB will multiply the sequence by N. 146C 147C The array WSAVE which is used by subroutine CFFTF must be 148C initialized by calling subroutine CFFTI(N,WSAVE). 149CPE 150CAS 151C Input Parameters 152C 153C 154C N the length of the complex sequence C. The method is 155C more efficient when N is the product of small primes. 156C 157C C a complex array of length N which contains the sequence 158C 159C WSAVE a real work array which must be dimensioned at least 4*N+15 160C in the program that calls CFFTF. The WSAVE array must be 161C initialized by calling subroutine CFFTI(N,WSAVE), and a 162C different WSAVE array must be used for each different 163C value of N. This initialization does not have to be 164C repeated so long as N remains unchanged. Thus subsequent 165C transforms can be obtained faster than the first. 166C The same WSAVE array can be used by CFFTF and CFFTB. 167C 168C Output Parameters 169C 170C C for J=1,...,N 171C 172C C(J)=the sum from K=1,...,N of 173C 174C C(K)*EXP(-I*J*K*2*PI/N) 175C 176C where I=SQRT(-1) 177C 178C WSAVE contains initialization calculations which must not be 179C destroyed between calls of subroutine CFFTF or CFFTB 180CAE 181C***REFERENCES (NONE) 182C***ROUTINES CALLED CFFTF1 183C***END PROLOGUE CFFTF 184 DIMENSION C(1) ,WSAVE(1) 185C***FIRST EXECUTABLE STATEMENT CFFTF 186 IF (N .EQ. 1) RETURN 187 IW1 = N+N+1 188 IW2 = IW1+N+N 189 CALL CFFTF1 (N,C,WSAVE,WSAVE(IW1),WSAVE(IW2)) 190 RETURN 191 END 192* 193 SUBROUTINE CFFTF1(N,C,CH,WA,IFAC) 194CSE 195 IMPLICIT REAL*8(A-H,O-Z) 196C***BEGIN PROLOGUE CFFTF1 197C***REFER TO CFFTF 198C***ROUTINES CALLED PASSF,PASSF2,PASSF3,PASSF4,PASSF5 199C***END PROLOGUE CFFTF1 200 DIMENSION CH(1) ,C(1) ,WA(1) ,IFAC(2) 201C***FIRST EXECUTABLE STATEMENT CFFTF1 202 NF = IFAC(2) 203 NA = 0 204 L1 = 1 205 IW = 1 206 DO 116 K1=1,NF 207 IP = IFAC(K1+2) 208 L2 = IP*L1 209 IDO = N/L2 210 IDOT = IDO+IDO 211 IDL1 = IDOT*L1 212 IF (IP .NE. 4) GO TO 103 213 IX2 = IW+IDOT 214 IX3 = IX2+IDOT 215 IF (NA .NE. 0) GO TO 101 216 CALL PASSF4 (IDOT,L1,C,CH,WA(IW),WA(IX2),WA(IX3)) 217 GO TO 102 218 101 CALL PASSF4 (IDOT,L1,CH,C,WA(IW),WA(IX2),WA(IX3)) 219 102 NA = 1-NA 220 GO TO 115 221 103 IF (IP .NE. 2) GO TO 106 222 IF (NA .NE. 0) GO TO 104 223 CALL PASSF2 (IDOT,L1,C,CH,WA(IW)) 224 GO TO 105 225 104 CALL PASSF2 (IDOT,L1,CH,C,WA(IW)) 226 105 NA = 1-NA 227 GO TO 115 228 106 IF (IP .NE. 3) GO TO 109 229 IX2 = IW+IDOT 230 IF (NA .NE. 0) GO TO 107 231 CALL PASSF3 (IDOT,L1,C,CH,WA(IW),WA(IX2)) 232 GO TO 108 233 107 CALL PASSF3 (IDOT,L1,CH,C,WA(IW),WA(IX2)) 234 108 NA = 1-NA 235 GO TO 115 236 109 IF (IP .NE. 5) GO TO 112 237 IX2 = IW+IDOT 238 IX3 = IX2+IDOT 239 IX4 = IX3+IDOT 240 IF (NA .NE. 0) GO TO 110 241 CALL PASSF5 (IDOT,L1,C,CH,WA(IW),WA(IX2),WA(IX3),WA(IX4)) 242 GO TO 111 243 110 CALL PASSF5 (IDOT,L1,CH,C,WA(IW),WA(IX2),WA(IX3),WA(IX4)) 244 111 NA = 1-NA 245 GO TO 115 246 112 IF (NA .NE. 0) GO TO 113 247 CALL PASSF (NAC,IDOT,IP,L1,IDL1,C,C,C,CH,CH,WA(IW)) 248 GO TO 114 249 113 CALL PASSF (NAC,IDOT,IP,L1,IDL1,CH,CH,CH,C,C,WA(IW)) 250 114 IF (NAC .NE. 0) NA = 1-NA 251 115 L1 = L2 252 IW = IW+(IP-1)*IDOT 253 116 CONTINUE 254 IF (NA .EQ. 0) RETURN 255 N2 = N+N 256 DO 117 I=1,N2 257 C(I) = CH(I) 258 117 CONTINUE 259 RETURN 260 END 261* 262 SUBROUTINE PASSF(NAC,IDO,IP,L1,IDL1,CC,C1,C2,CH,CH2,WA) 263CSE 264 IMPLICIT REAL*8(A-H,O-Z) 265C***BEGIN PROLOGUE PASSF 266C***REFER TO CFFTF 267C***ROUTINES CALLED (NONE) 268C***END PROLOGUE PASSF 269 DIMENSION CH(IDO,L1,IP) ,CC(IDO,IP,L1) , 270 1 C1(IDO,L1,IP) ,WA(1) ,C2(IDL1,IP), 271 2 CH2(IDL1,IP) 272C***FIRST EXECUTABLE STATEMENT PASSF 273 IDOT = IDO/2 274 NT = IP*IDL1 275 IPP2 = IP+2 276 IPPH = (IP+1)/2 277 IDP = IP*IDO 278C 279 IF (IDO .LT. L1) GO TO 106 280 DO 103 J=2,IPPH 281 JC = IPP2-J 282 DO 102 K=1,L1 283CDIR$ IVDEP 284 DO 101 I=1,IDO 285 CH(I,K,J) = CC(I,J,K)+CC(I,JC,K) 286 CH(I,K,JC) = CC(I,J,K)-CC(I,JC,K) 287 101 CONTINUE 288 102 CONTINUE 289 103 CONTINUE 290 DO 105 K=1,L1 291CDIR$ IVDEP 292 DO 104 I=1,IDO 293 CH(I,K,1) = CC(I,1,K) 294 104 CONTINUE 295 105 CONTINUE 296 GO TO 112 297 106 DO 109 J=2,IPPH 298 JC = IPP2-J 299 DO 108 I=1,IDO 300CDIR$ IVDEP 301 DO 107 K=1,L1 302 CH(I,K,J) = CC(I,J,K)+CC(I,JC,K) 303 CH(I,K,JC) = CC(I,J,K)-CC(I,JC,K) 304 107 CONTINUE 305 108 CONTINUE 306 109 CONTINUE 307 DO 111 I=1,IDO 308CDIR$ IVDEP 309 DO 110 K=1,L1 310 CH(I,K,1) = CC(I,1,K) 311 110 CONTINUE 312 111 CONTINUE 313 112 IDL = 2-IDO 314 INC = 0 315 DO 116 L=2,IPPH 316 LC = IPP2-L 317 IDL = IDL+IDO 318CDIR$ IVDEP 319 DO 113 IK=1,IDL1 320 C2(IK,L) = CH2(IK,1)+WA(IDL-1)*CH2(IK,2) 321 C2(IK,LC) = -WA(IDL)*CH2(IK,IP) 322 113 CONTINUE 323 IDLJ = IDL 324 INC = INC+IDO 325 DO 115 J=3,IPPH 326 JC = IPP2-J 327 IDLJ = IDLJ+INC 328 IF (IDLJ .GT. IDP) IDLJ = IDLJ-IDP 329 WAR = WA(IDLJ-1) 330 WAI = WA(IDLJ) 331CDIR$ IVDEP 332 DO 114 IK=1,IDL1 333 C2(IK,L) = C2(IK,L)+WAR*CH2(IK,J) 334 C2(IK,LC) = C2(IK,LC)-WAI*CH2(IK,JC) 335 114 CONTINUE 336 115 CONTINUE 337 116 CONTINUE 338 DO 118 J=2,IPPH 339CDIR$ IVDEP 340 DO 117 IK=1,IDL1 341 CH2(IK,1) = CH2(IK,1)+CH2(IK,J) 342 117 CONTINUE 343 118 CONTINUE 344 DO 120 J=2,IPPH 345 JC = IPP2-J 346CDIR$ IVDEP 347 DO 119 IK=2,IDL1,2 348 CH2(IK-1,J) = C2(IK-1,J)-C2(IK,JC) 349 CH2(IK-1,JC) = C2(IK-1,J)+C2(IK,JC) 350 CH2(IK,J) = C2(IK,J)+C2(IK-1,JC) 351 CH2(IK,JC) = C2(IK,J)-C2(IK-1,JC) 352 119 CONTINUE 353 120 CONTINUE 354 NAC = 1 355 IF (IDO .EQ. 2) RETURN 356 NAC = 0 357CDIR$ IVDEP 358 DO 121 IK=1,IDL1 359 C2(IK,1) = CH2(IK,1) 360 121 CONTINUE 361 DO 123 J=2,IP 362CDIR$ IVDEP 363 DO 122 K=1,L1 364 C1(1,K,J) = CH(1,K,J) 365 C1(2,K,J) = CH(2,K,J) 366 122 CONTINUE 367 123 CONTINUE 368 IF (IDOT .GT. L1) GO TO 127 369 IDIJ = 0 370 DO 126 J=2,IP 371 IDIJ = IDIJ+2 372 DO 125 I=4,IDO,2 373 IDIJ = IDIJ+2 374CDIR$ IVDEP 375 DO 124 K=1,L1 376 C1(I-1,K,J) = WA(IDIJ-1)*CH(I-1,K,J)+WA(IDIJ)*CH(I,K,J) 377 C1(I,K,J) = WA(IDIJ-1)*CH(I,K,J)-WA(IDIJ)*CH(I-1,K,J) 378 124 CONTINUE 379 125 CONTINUE 380 126 CONTINUE 381 RETURN 382 127 IDJ = 2-IDO 383 DO 130 J=2,IP 384 IDJ = IDJ+IDO 385 DO 129 K=1,L1 386 IDIJ = IDJ 387CDIR$ IVDEP 388 DO 128 I=4,IDO,2 389 IDIJ = IDIJ+2 390 C1(I-1,K,J) = WA(IDIJ-1)*CH(I-1,K,J)+WA(IDIJ)*CH(I,K,J) 391 C1(I,K,J) = WA(IDIJ-1)*CH(I,K,J)-WA(IDIJ)*CH(I-1,K,J) 392 128 CONTINUE 393 129 CONTINUE 394 130 CONTINUE 395 RETURN 396 END 397* 398 SUBROUTINE PASSF2(IDO,L1,CC,CH,WA1) 399CSE 400 IMPLICIT REAL*8(A-H,O-Z) 401C***BEGIN PROLOGUE PASSF2 402C***REFER TO CFFTF 403C***ROUTINES CALLED (NONE) 404C***END PROLOGUE PASSF2 405 DIMENSION CC(IDO,2,L1) ,CH(IDO,L1,2) , 406 1 WA1(1) 407C***FIRST EXECUTABLE STATEMENT PASSF2 408 IF (IDO .GT. 2) GO TO 102 409 DO 101 K=1,L1 410 CH(1,K,1) = CC(1,1,K)+CC(1,2,K) 411 CH(1,K,2) = CC(1,1,K)-CC(1,2,K) 412 CH(2,K,1) = CC(2,1,K)+CC(2,2,K) 413 CH(2,K,2) = CC(2,1,K)-CC(2,2,K) 414 101 CONTINUE 415 RETURN 416 102 IF(IDO/2.LT.L1) GO TO 105 417 DO 104 K=1,L1 418CDIR$ IVDEP 419 DO 103 I=2,IDO,2 420 CH(I-1,K,1) = CC(I-1,1,K)+CC(I-1,2,K) 421 TR2 = CC(I-1,1,K)-CC(I-1,2,K) 422 CH(I,K,1) = CC(I,1,K)+CC(I,2,K) 423 TI2 = CC(I,1,K)-CC(I,2,K) 424 CH(I,K,2) = WA1(I-1)*TI2-WA1(I)*TR2 425 CH(I-1,K,2) = WA1(I-1)*TR2+WA1(I)*TI2 426 103 CONTINUE 427 104 CONTINUE 428 RETURN 429 105 DO 107 I=2,IDO,2 430CDIR$ IVDEP 431 DO 106 K=1,L1 432 CH(I-1,K,1) = CC(I-1,1,K)+CC(I-1,2,K) 433 TR2 = CC(I-1,1,K)-CC(I-1,2,K) 434 CH(I,K,1) = CC(I,1,K)+CC(I,2,K) 435 TI2 = CC(I,1,K)-CC(I,2,K) 436 CH(I,K,2) = WA1(I-1)*TI2-WA1(I)*TR2 437 CH(I-1,K,2) = WA1(I-1)*TR2+WA1(I)*TI2 438 106 CONTINUE 439 107 CONTINUE 440 RETURN 441 END 442* 443 SUBROUTINE PASSF3(IDO,L1,CC,CH,WA1,WA2) 444CSE 445 IMPLICIT REAL*8(A-H,O-Z) 446C***BEGIN PROLOGUE PASSF3 447C***REFER TO CFFTF 448C***ROUTINES CALLED (NONE) 449C***END PROLOGUE PASSF3 450 DIMENSION CC(IDO,3,L1) ,CH(IDO,L1,3) , 451 1 WA1(1) ,WA2(1) 452 DATA TAUR,TAUI /-.5,-.866025403784439/ 453C***FIRST EXECUTABLE STATEMENT PASSF3 454 IF (IDO .NE. 2) GO TO 102 455 DO 101 K=1,L1 456 TR2 = CC(1,2,K)+CC(1,3,K) 457 CR2 = CC(1,1,K)+TAUR*TR2 458 CH(1,K,1) = CC(1,1,K)+TR2 459 TI2 = CC(2,2,K)+CC(2,3,K) 460 CI2 = CC(2,1,K)+TAUR*TI2 461 CH(2,K,1) = CC(2,1,K)+TI2 462 CR3 = TAUI*(CC(1,2,K)-CC(1,3,K)) 463 CI3 = TAUI*(CC(2,2,K)-CC(2,3,K)) 464 CH(1,K,2) = CR2-CI3 465 CH(1,K,3) = CR2+CI3 466 CH(2,K,2) = CI2+CR3 467 CH(2,K,3) = CI2-CR3 468 101 CONTINUE 469 RETURN 470 102 IF(IDO/2.LT.L1) GO TO 105 471 DO 104 K=1,L1 472CDIR$ IVDEP 473 DO 103 I=2,IDO,2 474 TR2 = CC(I-1,2,K)+CC(I-1,3,K) 475 CR2 = CC(I-1,1,K)+TAUR*TR2 476 CH(I-1,K,1) = CC(I-1,1,K)+TR2 477 TI2 = CC(I,2,K)+CC(I,3,K) 478 CI2 = CC(I,1,K)+TAUR*TI2 479 CH(I,K,1) = CC(I,1,K)+TI2 480 CR3 = TAUI*(CC(I-1,2,K)-CC(I-1,3,K)) 481 CI3 = TAUI*(CC(I,2,K)-CC(I,3,K)) 482 DR2 = CR2-CI3 483 DR3 = CR2+CI3 484 DI2 = CI2+CR3 485 DI3 = CI2-CR3 486 CH(I,K,2) = WA1(I-1)*DI2-WA1(I)*DR2 487 CH(I-1,K,2) = WA1(I-1)*DR2+WA1(I)*DI2 488 CH(I,K,3) = WA2(I-1)*DI3-WA2(I)*DR3 489 CH(I-1,K,3) = WA2(I-1)*DR3+WA2(I)*DI3 490 103 CONTINUE 491 104 CONTINUE 492 RETURN 493 105 DO 107 I=2,IDO,2 494CDIR$ IVDEP 495 DO 106 K=1,L1 496 TR2 = CC(I-1,2,K)+CC(I-1,3,K) 497 CR2 = CC(I-1,1,K)+TAUR*TR2 498 CH(I-1,K,1) = CC(I-1,1,K)+TR2 499 TI2 = CC(I,2,K)+CC(I,3,K) 500 CI2 = CC(I,1,K)+TAUR*TI2 501 CH(I,K,1) = CC(I,1,K)+TI2 502 CR3 = TAUI*(CC(I-1,2,K)-CC(I-1,3,K)) 503 CI3 = TAUI*(CC(I,2,K)-CC(I,3,K)) 504 DR2 = CR2-CI3 505 DR3 = CR2+CI3 506 DI2 = CI2+CR3 507 DI3 = CI2-CR3 508 CH(I,K,2) = WA1(I-1)*DI2-WA1(I)*DR2 509 CH(I-1,K,2) = WA1(I-1)*DR2+WA1(I)*DI2 510 CH(I,K,3) = WA2(I-1)*DI3-WA2(I)*DR3 511 CH(I-1,K,3) = WA2(I-1)*DR3+WA2(I)*DI3 512 106 CONTINUE 513 107 CONTINUE 514 RETURN 515 END 516* 517 SUBROUTINE PASSF4(IDO,L1,CC,CH,WA1,WA2,WA3) 518CSE 519 IMPLICIT REAL*8(A-H,O-Z) 520C***BEGIN PROLOGUE PASSF4 521C***REFER TO CFFTF 522C***ROUTINES CALLED (NONE) 523C***END PROLOGUE PASSF4 524 DIMENSION CC(IDO,4,L1) ,CH(IDO,L1,4) , 525 1 WA1(1) ,WA2(1) ,WA3(1) 526C***FIRST EXECUTABLE STATEMENT PASSF4 527 IF (IDO .NE. 2) GO TO 102 528 DO 101 K=1,L1 529 TI1 = CC(2,1,K)-CC(2,3,K) 530 TI2 = CC(2,1,K)+CC(2,3,K) 531 TR4 = CC(2,2,K)-CC(2,4,K) 532 TI3 = CC(2,2,K)+CC(2,4,K) 533 TR1 = CC(1,1,K)-CC(1,3,K) 534 TR2 = CC(1,1,K)+CC(1,3,K) 535 TI4 = CC(1,4,K)-CC(1,2,K) 536 TR3 = CC(1,2,K)+CC(1,4,K) 537 CH(1,K,1) = TR2+TR3 538 CH(1,K,3) = TR2-TR3 539 CH(2,K,1) = TI2+TI3 540 CH(2,K,3) = TI2-TI3 541 CH(1,K,2) = TR1+TR4 542 CH(1,K,4) = TR1-TR4 543 CH(2,K,2) = TI1+TI4 544 CH(2,K,4) = TI1-TI4 545 101 CONTINUE 546 RETURN 547 102 IF(IDO/2.LT.L1) GO TO 105 548 DO 104 K=1,L1 549CDIR$ IVDEP 550 DO 103 I=2,IDO,2 551 TI1 = CC(I,1,K)-CC(I,3,K) 552 TI2 = CC(I,1,K)+CC(I,3,K) 553 TI3 = CC(I,2,K)+CC(I,4,K) 554 TR4 = CC(I,2,K)-CC(I,4,K) 555 TR1 = CC(I-1,1,K)-CC(I-1,3,K) 556 TR2 = CC(I-1,1,K)+CC(I-1,3,K) 557 TI4 = CC(I-1,4,K)-CC(I-1,2,K) 558 TR3 = CC(I-1,2,K)+CC(I-1,4,K) 559 CH(I-1,K,1) = TR2+TR3 560 CR3 = TR2-TR3 561 CH(I,K,1) = TI2+TI3 562 CI3 = TI2-TI3 563 CR2 = TR1+TR4 564 CR4 = TR1-TR4 565 CI2 = TI1+TI4 566 CI4 = TI1-TI4 567 CH(I-1,K,2) = WA1(I-1)*CR2+WA1(I)*CI2 568 CH(I,K,2) = WA1(I-1)*CI2-WA1(I)*CR2 569 CH(I-1,K,3) = WA2(I-1)*CR3+WA2(I)*CI3 570 CH(I,K,3) = WA2(I-1)*CI3-WA2(I)*CR3 571 CH(I-1,K,4) = WA3(I-1)*CR4+WA3(I)*CI4 572 CH(I,K,4) = WA3(I-1)*CI4-WA3(I)*CR4 573 103 CONTINUE 574 104 CONTINUE 575 RETURN 576 105 DO 107 I=2,IDO,2 577CDIR$ IVDEP 578 DO 106 K=1,L1 579 TI1 = CC(I,1,K)-CC(I,3,K) 580 TI2 = CC(I,1,K)+CC(I,3,K) 581 TI3 = CC(I,2,K)+CC(I,4,K) 582 TR4 = CC(I,2,K)-CC(I,4,K) 583 TR1 = CC(I-1,1,K)-CC(I-1,3,K) 584 TR2 = CC(I-1,1,K)+CC(I-1,3,K) 585 TI4 = CC(I-1,4,K)-CC(I-1,2,K) 586 TR3 = CC(I-1,2,K)+CC(I-1,4,K) 587 CH(I-1,K,1) = TR2+TR3 588 CR3 = TR2-TR3 589 CH(I,K,1) = TI2+TI3 590 CI3 = TI2-TI3 591 CR2 = TR1+TR4 592 CR4 = TR1-TR4 593 CI2 = TI1+TI4 594 CI4 = TI1-TI4 595 CH(I-1,K,2) = WA1(I-1)*CR2+WA1(I)*CI2 596 CH(I,K,2) = WA1(I-1)*CI2-WA1(I)*CR2 597 CH(I-1,K,3) = WA2(I-1)*CR3+WA2(I)*CI3 598 CH(I,K,3) = WA2(I-1)*CI3-WA2(I)*CR3 599 CH(I-1,K,4) = WA3(I-1)*CR4+WA3(I)*CI4 600 CH(I,K,4) = WA3(I-1)*CI4-WA3(I)*CR4 601 106 CONTINUE 602 107 CONTINUE 603 RETURN 604 END 605* 606 SUBROUTINE PASSF5(IDO,L1,CC,CH,WA1,WA2,WA3,WA4) 607CSE 608 IMPLICIT REAL*8(A-H,O-Z) 609C***BEGIN PROLOGUE PASSF5 610C***REFER TO CFFTF 611C***ROUTINES CALLED (NONE) 612C***END PROLOGUE PASSF5 613 DIMENSION CC(IDO,5,L1) ,CH(IDO,L1,5) , 614 1 WA1(1) ,WA2(1) ,WA3(1) ,WA4(1) 615 DATA TR11,TI11,TR12,TI12 /.309016994374947,-.951056516295154, 616 1-.809016994374947,-.587785252292473/ 617C***FIRST EXECUTABLE STATEMENT PASSF5 618 IF (IDO .NE. 2) GO TO 102 619 DO 101 K=1,L1 620 TI5 = CC(2,2,K)-CC(2,5,K) 621 TI2 = CC(2,2,K)+CC(2,5,K) 622 TI4 = CC(2,3,K)-CC(2,4,K) 623 TI3 = CC(2,3,K)+CC(2,4,K) 624 TR5 = CC(1,2,K)-CC(1,5,K) 625 TR2 = CC(1,2,K)+CC(1,5,K) 626 TR4 = CC(1,3,K)-CC(1,4,K) 627 TR3 = CC(1,3,K)+CC(1,4,K) 628 CH(1,K,1) = CC(1,1,K)+TR2+TR3 629 CH(2,K,1) = CC(2,1,K)+TI2+TI3 630 CR2 = CC(1,1,K)+TR11*TR2+TR12*TR3 631 CI2 = CC(2,1,K)+TR11*TI2+TR12*TI3 632 CR3 = CC(1,1,K)+TR12*TR2+TR11*TR3 633 CI3 = CC(2,1,K)+TR12*TI2+TR11*TI3 634 CR5 = TI11*TR5+TI12*TR4 635 CI5 = TI11*TI5+TI12*TI4 636 CR4 = TI12*TR5-TI11*TR4 637 CI4 = TI12*TI5-TI11*TI4 638 CH(1,K,2) = CR2-CI5 639 CH(1,K,5) = CR2+CI5 640 CH(2,K,2) = CI2+CR5 641 CH(2,K,3) = CI3+CR4 642 CH(1,K,3) = CR3-CI4 643 CH(1,K,4) = CR3+CI4 644 CH(2,K,4) = CI3-CR4 645 CH(2,K,5) = CI2-CR5 646 101 CONTINUE 647 RETURN 648 102 IF(IDO/2.LT.L1) GO TO 105 649 DO 104 K=1,L1 650CDIR$ IVDEP 651 DO 103 I=2,IDO,2 652 TI5 = CC(I,2,K)-CC(I,5,K) 653 TI2 = CC(I,2,K)+CC(I,5,K) 654 TI4 = CC(I,3,K)-CC(I,4,K) 655 TI3 = CC(I,3,K)+CC(I,4,K) 656 TR5 = CC(I-1,2,K)-CC(I-1,5,K) 657 TR2 = CC(I-1,2,K)+CC(I-1,5,K) 658 TR4 = CC(I-1,3,K)-CC(I-1,4,K) 659 TR3 = CC(I-1,3,K)+CC(I-1,4,K) 660 CH(I-1,K,1) = CC(I-1,1,K)+TR2+TR3 661 CH(I,K,1) = CC(I,1,K)+TI2+TI3 662 CR2 = CC(I-1,1,K)+TR11*TR2+TR12*TR3 663 CI2 = CC(I,1,K)+TR11*TI2+TR12*TI3 664 CR3 = CC(I-1,1,K)+TR12*TR2+TR11*TR3 665 CI3 = CC(I,1,K)+TR12*TI2+TR11*TI3 666 CR5 = TI11*TR5+TI12*TR4 667 CI5 = TI11*TI5+TI12*TI4 668 CR4 = TI12*TR5-TI11*TR4 669 CI4 = TI12*TI5-TI11*TI4 670 DR3 = CR3-CI4 671 DR4 = CR3+CI4 672 DI3 = CI3+CR4 673 DI4 = CI3-CR4 674 DR5 = CR2+CI5 675 DR2 = CR2-CI5 676 DI5 = CI2-CR5 677 DI2 = CI2+CR5 678 CH(I-1,K,2) = WA1(I-1)*DR2+WA1(I)*DI2 679 CH(I,K,2) = WA1(I-1)*DI2-WA1(I)*DR2 680 CH(I-1,K,3) = WA2(I-1)*DR3+WA2(I)*DI3 681 CH(I,K,3) = WA2(I-1)*DI3-WA2(I)*DR3 682 CH(I-1,K,4) = WA3(I-1)*DR4+WA3(I)*DI4 683 CH(I,K,4) = WA3(I-1)*DI4-WA3(I)*DR4 684 CH(I-1,K,5) = WA4(I-1)*DR5+WA4(I)*DI5 685 CH(I,K,5) = WA4(I-1)*DI5-WA4(I)*DR5 686 103 CONTINUE 687 104 CONTINUE 688 RETURN 689 105 DO 107 I=2,IDO,2 690CDIR$ IVDEP 691 DO 106 K=1,L1 692 TI5 = CC(I,2,K)-CC(I,5,K) 693 TI2 = CC(I,2,K)+CC(I,5,K) 694 TI4 = CC(I,3,K)-CC(I,4,K) 695 TI3 = CC(I,3,K)+CC(I,4,K) 696 TR5 = CC(I-1,2,K)-CC(I-1,5,K) 697 TR2 = CC(I-1,2,K)+CC(I-1,5,K) 698 TR4 = CC(I-1,3,K)-CC(I-1,4,K) 699 TR3 = CC(I-1,3,K)+CC(I-1,4,K) 700 CH(I-1,K,1) = CC(I-1,1,K)+TR2+TR3 701 CH(I,K,1) = CC(I,1,K)+TI2+TI3 702 CR2 = CC(I-1,1,K)+TR11*TR2+TR12*TR3 703 CI2 = CC(I,1,K)+TR11*TI2+TR12*TI3 704 CR3 = CC(I-1,1,K)+TR12*TR2+TR11*TR3 705 CI3 = CC(I,1,K)+TR12*TI2+TR11*TI3 706 CR5 = TI11*TR5+TI12*TR4 707 CI5 = TI11*TI5+TI12*TI4 708 CR4 = TI12*TR5-TI11*TR4 709 CI4 = TI12*TI5-TI11*TI4 710 DR3 = CR3-CI4 711 DR4 = CR3+CI4 712 DI3 = CI3+CR4 713 DI4 = CI3-CR4 714 DR5 = CR2+CI5 715 DR2 = CR2-CI5 716 DI5 = CI2-CR5 717 DI2 = CI2+CR5 718 CH(I-1,K,2) = WA1(I-1)*DR2+WA1(I)*DI2 719 CH(I,K,2) = WA1(I-1)*DI2-WA1(I)*DR2 720 CH(I-1,K,3) = WA2(I-1)*DR3+WA2(I)*DI3 721 CH(I,K,3) = WA2(I-1)*DI3-WA2(I)*DR3 722 CH(I-1,K,4) = WA3(I-1)*DR4+WA3(I)*DI4 723 CH(I,K,4) = WA3(I-1)*DI4-WA3(I)*DR4 724 CH(I-1,K,5) = WA4(I-1)*DR5+WA4(I)*DI5 725 CH(I,K,5) = WA4(I-1)*DI5-WA4(I)*DR5 726 106 CONTINUE 727 107 CONTINUE 728 RETURN 729 END 730* 731 SUBROUTINE CFFTB(N,C,WSAVE) 732CSE 733 IMPLICIT REAL*8(A-H,O-Z) 734C***BEGIN PROLOGUE CFFTB 735C***DATE WRITTEN 790601 (YYMMDD) 736C***REVISION DATE 830401 (YYMMDD) 737C***CATEGORY NO. J1A2 738C***KEYWORDS FOURIER TRANSFORM 739C***AUTHOR SWARZTRAUBER, P. N., (NCAR) 740CPS 741C***PURPOSE Unnormalized inverse of CFFTF. 742C***DESCRIPTION 743C 744C Subroutine CFFTB computes the backward complex discrete Fourier 745C transform (the Fourier synthesis). Equivalently, CFFTB computes 746C a complex periodic sequence from its Fourier coefficients. 747C The transform is defined below at output parameter C. 748C 749C A call of CFFTF followed by a call of CFFTB will multiply the 750C sequence by N. 751C 752C The array WSAVE which is used by subroutine CFFTB must be 753C initialized by calling subroutine CFFTI(N,WSAVE). 754CPE 755CAS 756C Input Parameters 757C 758C 759C N the length of the complex sequence C. The method is 760C more efficient when N is the product of small primes. 761C 762C C a complex array of length N which contains the sequence 763C 764C WSAVE a real work array which must be dimensioned at least 4*N+15 765C in the program that calls CFFTB. The WSAVE array must be 766C initialized by calling subroutine CFFTI(N,WSAVE), and a 767C different WSAVE array must be used for each different 768C value of N. This initialization does not have to be 769C repeated so long as N remains unchanged. Thus subsequent 770C transforms can be obtained faster than the first. 771C The same WSAVE array can be used by CFFTF and CFFTB. 772C 773C Output Parameters 774C 775C C For J=1,...,N 776C 777C C(J)=the sum from K=1,...,N of 778C 779C C(K)*EXP(I*J*K*2*PI/N) 780C 781C where I=SQRT(-1) 782C 783C WSAVE contains initialization calculations which must not be 784C destroyed between calls of subroutine CFFTF or CFFTB 785CAE 786C***REFERENCES (NONE) 787C***ROUTINES CALLED CFFTB1 788C***END PROLOGUE CFFTB 789 DIMENSION C(1) ,WSAVE(1) 790C***FIRST EXECUTABLE STATEMENT CFFTB 791 IF (N .EQ. 1) RETURN 792 IW1 = N+N+1 793 IW2 = IW1+N+N 794 CALL CFFTB1 (N,C,WSAVE,WSAVE(IW1),WSAVE(IW2)) 795 RETURN 796 END 797* 798 SUBROUTINE CFFTB1(N,C,CH,WA,IFAC) 799CSE 800 IMPLICIT REAL*8(A-H,O-Z) 801C***BEGIN PROLOGUE CFFTB1 802C***REFER TO CFFTB 803C***ROUTINES CALLED PASSB,PASSB2,PASSB3,PASSB4,PASSB5 804C***END PROLOGUE CFFTB1 805 DIMENSION CH(1) ,C(1) ,WA(1) ,IFAC(2) 806C***FIRST EXECUTABLE STATEMENT CFFTB1 807 NF = IFAC(2) 808 NA = 0 809 L1 = 1 810 IW = 1 811 DO 116 K1=1,NF 812 IP = IFAC(K1+2) 813 L2 = IP*L1 814 IDO = N/L2 815 IDOT = IDO+IDO 816 IDL1 = IDOT*L1 817 IF (IP .NE. 4) GO TO 103 818 IX2 = IW+IDOT 819 IX3 = IX2+IDOT 820 IF (NA .NE. 0) GO TO 101 821 CALL PASSB4 (IDOT,L1,C,CH,WA(IW),WA(IX2),WA(IX3)) 822 GO TO 102 823 101 CALL PASSB4 (IDOT,L1,CH,C,WA(IW),WA(IX2),WA(IX3)) 824 102 NA = 1-NA 825 GO TO 115 826 103 IF (IP .NE. 2) GO TO 106 827 IF (NA .NE. 0) GO TO 104 828 CALL PASSB2 (IDOT,L1,C,CH,WA(IW)) 829 GO TO 105 830 104 CALL PASSB2 (IDOT,L1,CH,C,WA(IW)) 831 105 NA = 1-NA 832 GO TO 115 833 106 IF (IP .NE. 3) GO TO 109 834 IX2 = IW+IDOT 835 IF (NA .NE. 0) GO TO 107 836 CALL PASSB3 (IDOT,L1,C,CH,WA(IW),WA(IX2)) 837 GO TO 108 838 107 CALL PASSB3 (IDOT,L1,CH,C,WA(IW),WA(IX2)) 839 108 NA = 1-NA 840 GO TO 115 841 109 IF (IP .NE. 5) GO TO 112 842 IX2 = IW+IDOT 843 IX3 = IX2+IDOT 844 IX4 = IX3+IDOT 845 IF (NA .NE. 0) GO TO 110 846 CALL PASSB5 (IDOT,L1,C,CH,WA(IW),WA(IX2),WA(IX3),WA(IX4)) 847 GO TO 111 848 110 CALL PASSB5 (IDOT,L1,CH,C,WA(IW),WA(IX2),WA(IX3),WA(IX4)) 849 111 NA = 1-NA 850 GO TO 115 851 112 IF (NA .NE. 0) GO TO 113 852 CALL PASSB (NAC,IDOT,IP,L1,IDL1,C,C,C,CH,CH,WA(IW)) 853 GO TO 114 854 113 CALL PASSB (NAC,IDOT,IP,L1,IDL1,CH,CH,CH,C,C,WA(IW)) 855 114 IF (NAC .NE. 0) NA = 1-NA 856 115 L1 = L2 857 IW = IW+(IP-1)*IDOT 858 116 CONTINUE 859 IF (NA .EQ. 0) RETURN 860 N2 = N+N 861 DO 117 I=1,N2 862 C(I) = CH(I) 863 117 CONTINUE 864 RETURN 865 END 866* 867 SUBROUTINE PASSB(NAC,IDO,IP,L1,IDL1,CC,C1,C2,CH,CH2,WA) 868CSE 869 IMPLICIT REAL*8(A-H,O-Z) 870C***BEGIN PROLOGUE PASSB 871C***REFER TO CFFTB 872C***ROUTINES CALLED (NONE) 873C***END PROLOGUE PASSB 874 DIMENSION CH(IDO,L1,IP) ,CC(IDO,IP,L1) , 875 1 C1(IDO,L1,IP) ,WA(1) ,C2(IDL1,IP), 876 2 CH2(IDL1,IP) 877C***FIRST EXECUTABLE STATEMENT PASSB 878 IDOT = IDO/2 879 NT = IP*IDL1 880 IPP2 = IP+2 881 IPPH = (IP+1)/2 882 IDP = IP*IDO 883C 884 IF (IDO .LT. L1) GO TO 106 885 DO 103 J=2,IPPH 886 JC = IPP2-J 887 DO 102 K=1,L1 888CDIR$ IVDEP 889 DO 101 I=1,IDO 890 CH(I,K,J) = CC(I,J,K)+CC(I,JC,K) 891 CH(I,K,JC) = CC(I,J,K)-CC(I,JC,K) 892 101 CONTINUE 893 102 CONTINUE 894 103 CONTINUE 895 DO 105 K=1,L1 896CDIR$ IVDEP 897 DO 104 I=1,IDO 898 CH(I,K,1) = CC(I,1,K) 899 104 CONTINUE 900 105 CONTINUE 901 GO TO 112 902 106 DO 109 J=2,IPPH 903 JC = IPP2-J 904 DO 108 I=1,IDO 905CDIR$ IVDEP 906 DO 107 K=1,L1 907 CH(I,K,J) = CC(I,J,K)+CC(I,JC,K) 908 CH(I,K,JC) = CC(I,J,K)-CC(I,JC,K) 909 107 CONTINUE 910 108 CONTINUE 911 109 CONTINUE 912 DO 111 I=1,IDO 913CDIR$ IVDEP 914 DO 110 K=1,L1 915 CH(I,K,1) = CC(I,1,K) 916 110 CONTINUE 917 111 CONTINUE 918 112 IDL = 2-IDO 919 INC = 0 920 DO 116 L=2,IPPH 921 LC = IPP2-L 922 IDL = IDL+IDO 923CDIR$ IVDEP 924 DO 113 IK=1,IDL1 925 C2(IK,L) = CH2(IK,1)+WA(IDL-1)*CH2(IK,2) 926 C2(IK,LC) = WA(IDL)*CH2(IK,IP) 927 113 CONTINUE 928 IDLJ = IDL 929 INC = INC+IDO 930 DO 115 J=3,IPPH 931 JC = IPP2-J 932 IDLJ = IDLJ+INC 933 IF (IDLJ .GT. IDP) IDLJ = IDLJ-IDP 934 WAR = WA(IDLJ-1) 935 WAI = WA(IDLJ) 936CDIR$ IVDEP 937 DO 114 IK=1,IDL1 938 C2(IK,L) = C2(IK,L)+WAR*CH2(IK,J) 939 C2(IK,LC) = C2(IK,LC)+WAI*CH2(IK,JC) 940 114 CONTINUE 941 115 CONTINUE 942 116 CONTINUE 943 DO 118 J=2,IPPH 944CDIR$ IVDEP 945 DO 117 IK=1,IDL1 946 CH2(IK,1) = CH2(IK,1)+CH2(IK,J) 947 117 CONTINUE 948 118 CONTINUE 949 DO 120 J=2,IPPH 950 JC = IPP2-J 951CDIR$ IVDEP 952 DO 119 IK=2,IDL1,2 953 CH2(IK-1,J) = C2(IK-1,J)-C2(IK,JC) 954 CH2(IK-1,JC) = C2(IK-1,J)+C2(IK,JC) 955 CH2(IK,J) = C2(IK,J)+C2(IK-1,JC) 956 CH2(IK,JC) = C2(IK,J)-C2(IK-1,JC) 957 119 CONTINUE 958 120 CONTINUE 959 NAC = 1 960 IF (IDO .EQ. 2) RETURN 961 NAC = 0 962 DO 121 IK=1,IDL1 963 C2(IK,1) = CH2(IK,1) 964 121 CONTINUE 965 DO 123 J=2,IP 966CDIR$ IVDEP 967 DO 122 K=1,L1 968 C1(1,K,J) = CH(1,K,J) 969 C1(2,K,J) = CH(2,K,J) 970 122 CONTINUE 971 123 CONTINUE 972 IF (IDOT .GT. L1) GO TO 127 973 IDIJ = 0 974 DO 126 J=2,IP 975 IDIJ = IDIJ+2 976 DO 125 I=4,IDO,2 977 IDIJ = IDIJ+2 978CDIR$ IVDEP 979 DO 124 K=1,L1 980 C1(I-1,K,J) = WA(IDIJ-1)*CH(I-1,K,J)-WA(IDIJ)*CH(I,K,J) 981 C1(I,K,J) = WA(IDIJ-1)*CH(I,K,J)+WA(IDIJ)*CH(I-1,K,J) 982 124 CONTINUE 983 125 CONTINUE 984 126 CONTINUE 985 RETURN 986 127 IDJ = 2-IDO 987 DO 130 J=2,IP 988 IDJ = IDJ+IDO 989 DO 129 K=1,L1 990 IDIJ = IDJ 991CDIR$ IVDEP 992 DO 128 I=4,IDO,2 993 IDIJ = IDIJ+2 994 C1(I-1,K,J) = WA(IDIJ-1)*CH(I-1,K,J)-WA(IDIJ)*CH(I,K,J) 995 C1(I,K,J) = WA(IDIJ-1)*CH(I,K,J)+WA(IDIJ)*CH(I-1,K,J) 996 128 CONTINUE 997 129 CONTINUE 998 130 CONTINUE 999 RETURN 1000 END 1001* 1002 SUBROUTINE PASSB2(IDO,L1,CC,CH,WA1) 1003CSE 1004 IMPLICIT REAL*8(A-H,O-Z) 1005C***BEGIN PROLOGUE PASSB2 1006C***REFER TO CFFTB 1007C***ROUTINES CALLED (NONE) 1008C***END PROLOGUE PASSB2 1009 DIMENSION CC(IDO,2,L1) ,CH(IDO,L1,2) , 1010 1 WA1(1) 1011C***FIRST EXECUTABLE STATEMENT PASSB2 1012 IF (IDO .GT. 2) GO TO 102 1013 DO 101 K=1,L1 1014 CH(1,K,1) = CC(1,1,K)+CC(1,2,K) 1015 CH(1,K,2) = CC(1,1,K)-CC(1,2,K) 1016 CH(2,K,1) = CC(2,1,K)+CC(2,2,K) 1017 CH(2,K,2) = CC(2,1,K)-CC(2,2,K) 1018 101 CONTINUE 1019 RETURN 1020 102 IF(IDO/2.LT.L1) GO TO 105 1021 DO 104 K=1,L1 1022CDIR$ IVDEP 1023 DO 103 I=2,IDO,2 1024 CH(I-1,K,1) = CC(I-1,1,K)+CC(I-1,2,K) 1025 TR2 = CC(I-1,1,K)-CC(I-1,2,K) 1026 CH(I,K,1) = CC(I,1,K)+CC(I,2,K) 1027 TI2 = CC(I,1,K)-CC(I,2,K) 1028 CH(I,K,2) = WA1(I-1)*TI2+WA1(I)*TR2 1029 CH(I-1,K,2) = WA1(I-1)*TR2-WA1(I)*TI2 1030 103 CONTINUE 1031 104 CONTINUE 1032 RETURN 1033 105 DO 107 I=2,IDO,2 1034CDIR$ IVDEP 1035 DO 106 K=1,L1 1036 CH(I-1,K,1) = CC(I-1,1,K)+CC(I-1,2,K) 1037 TR2 = CC(I-1,1,K)-CC(I-1,2,K) 1038 CH(I,K,1) = CC(I,1,K)+CC(I,2,K) 1039 TI2 = CC(I,1,K)-CC(I,2,K) 1040 CH(I,K,2) = WA1(I-1)*TI2+WA1(I)*TR2 1041 CH(I-1,K,2) = WA1(I-1)*TR2-WA1(I)*TI2 1042 106 CONTINUE 1043 107 CONTINUE 1044 RETURN 1045 END 1046* 1047 SUBROUTINE PASSB3(IDO,L1,CC,CH,WA1,WA2) 1048CSE 1049 IMPLICIT REAL*8(A-H,O-Z) 1050C***BEGIN PROLOGUE PASSB3 1051C***REFER TO CFFTB 1052C***ROUTINES CALLED (NONE) 1053C***END PROLOGUE PASSB3 1054 DIMENSION CC(IDO,3,L1) ,CH(IDO,L1,3) , 1055 1 WA1(1) ,WA2(1) 1056 DATA TAUR,TAUI /-.5,.866025403784439/ 1057C***FIRST EXECUTABLE STATEMENT PASSB3 1058 IF (IDO .NE. 2) GO TO 102 1059 DO 101 K=1,L1 1060 TR2 = CC(1,2,K)+CC(1,3,K) 1061 CR2 = CC(1,1,K)+TAUR*TR2 1062 CH(1,K,1) = CC(1,1,K)+TR2 1063 TI2 = CC(2,2,K)+CC(2,3,K) 1064 CI2 = CC(2,1,K)+TAUR*TI2 1065 CH(2,K,1) = CC(2,1,K)+TI2 1066 CR3 = TAUI*(CC(1,2,K)-CC(1,3,K)) 1067 CI3 = TAUI*(CC(2,2,K)-CC(2,3,K)) 1068 CH(1,K,2) = CR2-CI3 1069 CH(1,K,3) = CR2+CI3 1070 CH(2,K,2) = CI2+CR3 1071 CH(2,K,3) = CI2-CR3 1072 101 CONTINUE 1073 RETURN 1074 102 IF(IDO/2.LT.L1) GO TO 105 1075 DO 104 K=1,L1 1076CDIR$ IVDEP 1077 DO 103 I=2,IDO,2 1078 TR2 = CC(I-1,2,K)+CC(I-1,3,K) 1079 CR2 = CC(I-1,1,K)+TAUR*TR2 1080 CH(I-1,K,1) = CC(I-1,1,K)+TR2 1081 TI2 = CC(I,2,K)+CC(I,3,K) 1082 CI2 = CC(I,1,K)+TAUR*TI2 1083 CH(I,K,1) = CC(I,1,K)+TI2 1084 CR3 = TAUI*(CC(I-1,2,K)-CC(I-1,3,K)) 1085 CI3 = TAUI*(CC(I,2,K)-CC(I,3,K)) 1086 DR2 = CR2-CI3 1087 DR3 = CR2+CI3 1088 DI2 = CI2+CR3 1089 DI3 = CI2-CR3 1090 CH(I,K,2) = WA1(I-1)*DI2+WA1(I)*DR2 1091 CH(I-1,K,2) = WA1(I-1)*DR2-WA1(I)*DI2 1092 CH(I,K,3) = WA2(I-1)*DI3+WA2(I)*DR3 1093 CH(I-1,K,3) = WA2(I-1)*DR3-WA2(I)*DI3 1094 103 CONTINUE 1095 104 CONTINUE 1096 RETURN 1097 105 DO 107 I=2,IDO,2 1098CDIR$ IVDEP 1099 DO 106 K=1,L1 1100 TR2 = CC(I-1,2,K)+CC(I-1,3,K) 1101 CR2 = CC(I-1,1,K)+TAUR*TR2 1102 CH(I-1,K,1) = CC(I-1,1,K)+TR2 1103 TI2 = CC(I,2,K)+CC(I,3,K) 1104 CI2 = CC(I,1,K)+TAUR*TI2 1105 CH(I,K,1) = CC(I,1,K)+TI2 1106 CR3 = TAUI*(CC(I-1,2,K)-CC(I-1,3,K)) 1107 CI3 = TAUI*(CC(I,2,K)-CC(I,3,K)) 1108 DR2 = CR2-CI3 1109 DR3 = CR2+CI3 1110 DI2 = CI2+CR3 1111 DI3 = CI2-CR3 1112 CH(I,K,2) = WA1(I-1)*DI2+WA1(I)*DR2 1113 CH(I-1,K,2) = WA1(I-1)*DR2-WA1(I)*DI2 1114 CH(I,K,3) = WA2(I-1)*DI3+WA2(I)*DR3 1115 CH(I-1,K,3) = WA2(I-1)*DR3-WA2(I)*DI3 1116 106 CONTINUE 1117 107 CONTINUE 1118 RETURN 1119 END 1120* 1121 SUBROUTINE PASSB4(IDO,L1,CC,CH,WA1,WA2,WA3) 1122CSE 1123 IMPLICIT REAL*8(A-H,O-Z) 1124C***BEGIN PROLOGUE PASSB4 1125C***REFER TO CFFTB 1126C***ROUTINES CALLED (NONE) 1127C***END PROLOGUE PASSB4 1128 DIMENSION CC(IDO,4,L1) ,CH(IDO,L1,4) , 1129 1 WA1(1) ,WA2(1) ,WA3(1) 1130C***FIRST EXECUTABLE STATEMENT PASSB4 1131 IF (IDO .NE. 2) GO TO 102 1132 DO 101 K=1,L1 1133 TI1 = CC(2,1,K)-CC(2,3,K) 1134 TI2 = CC(2,1,K)+CC(2,3,K) 1135 TR4 = CC(2,4,K)-CC(2,2,K) 1136 TI3 = CC(2,2,K)+CC(2,4,K) 1137 TR1 = CC(1,1,K)-CC(1,3,K) 1138 TR2 = CC(1,1,K)+CC(1,3,K) 1139 TI4 = CC(1,2,K)-CC(1,4,K) 1140 TR3 = CC(1,2,K)+CC(1,4,K) 1141 CH(1,K,1) = TR2+TR3 1142 CH(1,K,3) = TR2-TR3 1143 CH(2,K,1) = TI2+TI3 1144 CH(2,K,3) = TI2-TI3 1145 CH(1,K,2) = TR1+TR4 1146 CH(1,K,4) = TR1-TR4 1147 CH(2,K,2) = TI1+TI4 1148 CH(2,K,4) = TI1-TI4 1149 101 CONTINUE 1150 RETURN 1151 102 IF(IDO/2.LT.L1) GO TO 105 1152 DO 104 K=1,L1 1153CDIR$ IVDEP 1154 DO 103 I=2,IDO,2 1155 TI1 = CC(I,1,K)-CC(I,3,K) 1156 TI2 = CC(I,1,K)+CC(I,3,K) 1157 TI3 = CC(I,2,K)+CC(I,4,K) 1158 TR4 = CC(I,4,K)-CC(I,2,K) 1159 TR1 = CC(I-1,1,K)-CC(I-1,3,K) 1160 TR2 = CC(I-1,1,K)+CC(I-1,3,K) 1161 TI4 = CC(I-1,2,K)-CC(I-1,4,K) 1162 TR3 = CC(I-1,2,K)+CC(I-1,4,K) 1163 CH(I-1,K,1) = TR2+TR3 1164 CR3 = TR2-TR3 1165 CH(I,K,1) = TI2+TI3 1166 CI3 = TI2-TI3 1167 CR2 = TR1+TR4 1168 CR4 = TR1-TR4 1169 CI2 = TI1+TI4 1170 CI4 = TI1-TI4 1171 CH(I-1,K,2) = WA1(I-1)*CR2-WA1(I)*CI2 1172 CH(I,K,2) = WA1(I-1)*CI2+WA1(I)*CR2 1173 CH(I-1,K,3) = WA2(I-1)*CR3-WA2(I)*CI3 1174 CH(I,K,3) = WA2(I-1)*CI3+WA2(I)*CR3 1175 CH(I-1,K,4) = WA3(I-1)*CR4-WA3(I)*CI4 1176 CH(I,K,4) = WA3(I-1)*CI4+WA3(I)*CR4 1177 103 CONTINUE 1178 104 CONTINUE 1179 RETURN 1180 105 DO 107 I=2,IDO,2 1181CDIR$ IVDEP 1182 DO 106 K=1,L1 1183 TI1 = CC(I,1,K)-CC(I,3,K) 1184 TI2 = CC(I,1,K)+CC(I,3,K) 1185 TI3 = CC(I,2,K)+CC(I,4,K) 1186 TR4 = CC(I,4,K)-CC(I,2,K) 1187 TR1 = CC(I-1,1,K)-CC(I-1,3,K) 1188 TR2 = CC(I-1,1,K)+CC(I-1,3,K) 1189 TI4 = CC(I-1,2,K)-CC(I-1,4,K) 1190 TR3 = CC(I-1,2,K)+CC(I-1,4,K) 1191 CH(I-1,K,1) = TR2+TR3 1192 CR3 = TR2-TR3 1193 CH(I,K,1) = TI2+TI3 1194 CI3 = TI2-TI3 1195 CR2 = TR1+TR4 1196 CR4 = TR1-TR4 1197 CI2 = TI1+TI4 1198 CI4 = TI1-TI4 1199 CH(I-1,K,2) = WA1(I-1)*CR2-WA1(I)*CI2 1200 CH(I,K,2) = WA1(I-1)*CI2+WA1(I)*CR2 1201 CH(I-1,K,3) = WA2(I-1)*CR3-WA2(I)*CI3 1202 CH(I,K,3) = WA2(I-1)*CI3+WA2(I)*CR3 1203 CH(I-1,K,4) = WA3(I-1)*CR4-WA3(I)*CI4 1204 CH(I,K,4) = WA3(I-1)*CI4+WA3(I)*CR4 1205 106 CONTINUE 1206 107 CONTINUE 1207 RETURN 1208 END 1209* 1210 SUBROUTINE PASSB5(IDO,L1,CC,CH,WA1,WA2,WA3,WA4) 1211CSE 1212 IMPLICIT REAL*8(A-H,O-Z) 1213C***BEGIN PROLOGUE PASSB5 1214C***REFER TO CFFTB 1215C***ROUTINES CALLED (NONE) 1216C***END PROLOGUE PASSB5 1217 DIMENSION CC(IDO,5,L1) ,CH(IDO,L1,5) , 1218 1 WA1(1) ,WA2(1) ,WA3(1) ,WA4(1) 1219 DATA TR11,TI11,TR12,TI12 /.309016994374947,.951056516295154, 1220 1-.809016994374947,.587785252292473/ 1221C***FIRST EXECUTABLE STATEMENT PASSB5 1222 IF (IDO .NE. 2) GO TO 102 1223 DO 101 K=1,L1 1224 TI5 = CC(2,2,K)-CC(2,5,K) 1225 TI2 = CC(2,2,K)+CC(2,5,K) 1226 TI4 = CC(2,3,K)-CC(2,4,K) 1227 TI3 = CC(2,3,K)+CC(2,4,K) 1228 TR5 = CC(1,2,K)-CC(1,5,K) 1229 TR2 = CC(1,2,K)+CC(1,5,K) 1230 TR4 = CC(1,3,K)-CC(1,4,K) 1231 TR3 = CC(1,3,K)+CC(1,4,K) 1232 CH(1,K,1) = CC(1,1,K)+TR2+TR3 1233 CH(2,K,1) = CC(2,1,K)+TI2+TI3 1234 CR2 = CC(1,1,K)+TR11*TR2+TR12*TR3 1235 CI2 = CC(2,1,K)+TR11*TI2+TR12*TI3 1236 CR3 = CC(1,1,K)+TR12*TR2+TR11*TR3 1237 CI3 = CC(2,1,K)+TR12*TI2+TR11*TI3 1238 CR5 = TI11*TR5+TI12*TR4 1239 CI5 = TI11*TI5+TI12*TI4 1240 CR4 = TI12*TR5-TI11*TR4 1241 CI4 = TI12*TI5-TI11*TI4 1242 CH(1,K,2) = CR2-CI5 1243 CH(1,K,5) = CR2+CI5 1244 CH(2,K,2) = CI2+CR5 1245 CH(2,K,3) = CI3+CR4 1246 CH(1,K,3) = CR3-CI4 1247 CH(1,K,4) = CR3+CI4 1248 CH(2,K,4) = CI3-CR4 1249 CH(2,K,5) = CI2-CR5 1250 101 CONTINUE 1251 RETURN 1252 102 IF(IDO/2.LT.L1) GO TO 105 1253 DO 104 K=1,L1 1254CDIR$ IVDEP 1255 DO 103 I=2,IDO,2 1256 TI5 = CC(I,2,K)-CC(I,5,K) 1257 TI2 = CC(I,2,K)+CC(I,5,K) 1258 TI4 = CC(I,3,K)-CC(I,4,K) 1259 TI3 = CC(I,3,K)+CC(I,4,K) 1260 TR5 = CC(I-1,2,K)-CC(I-1,5,K) 1261 TR2 = CC(I-1,2,K)+CC(I-1,5,K) 1262 TR4 = CC(I-1,3,K)-CC(I-1,4,K) 1263 TR3 = CC(I-1,3,K)+CC(I-1,4,K) 1264 CH(I-1,K,1) = CC(I-1,1,K)+TR2+TR3 1265 CH(I,K,1) = CC(I,1,K)+TI2+TI3 1266 CR2 = CC(I-1,1,K)+TR11*TR2+TR12*TR3 1267 CI2 = CC(I,1,K)+TR11*TI2+TR12*TI3 1268 CR3 = CC(I-1,1,K)+TR12*TR2+TR11*TR3 1269 CI3 = CC(I,1,K)+TR12*TI2+TR11*TI3 1270 CR5 = TI11*TR5+TI12*TR4 1271 CI5 = TI11*TI5+TI12*TI4 1272 CR4 = TI12*TR5-TI11*TR4 1273 CI4 = TI12*TI5-TI11*TI4 1274 DR3 = CR3-CI4 1275 DR4 = CR3+CI4 1276 DI3 = CI3+CR4 1277 DI4 = CI3-CR4 1278 DR5 = CR2+CI5 1279 DR2 = CR2-CI5 1280 DI5 = CI2-CR5 1281 DI2 = CI2+CR5 1282 CH(I-1,K,2) = WA1(I-1)*DR2-WA1(I)*DI2 1283 CH(I,K,2) = WA1(I-1)*DI2+WA1(I)*DR2 1284 CH(I-1,K,3) = WA2(I-1)*DR3-WA2(I)*DI3 1285 CH(I,K,3) = WA2(I-1)*DI3+WA2(I)*DR3 1286 CH(I-1,K,4) = WA3(I-1)*DR4-WA3(I)*DI4 1287 CH(I,K,4) = WA3(I-1)*DI4+WA3(I)*DR4 1288 CH(I-1,K,5) = WA4(I-1)*DR5-WA4(I)*DI5 1289 CH(I,K,5) = WA4(I-1)*DI5+WA4(I)*DR5 1290 103 CONTINUE 1291 104 CONTINUE 1292 RETURN 1293 105 DO 107 I=2,IDO,2 1294CDIR$ IVDEP 1295 DO 106 K=1,L1 1296 TI5 = CC(I,2,K)-CC(I,5,K) 1297 TI2 = CC(I,2,K)+CC(I,5,K) 1298 TI4 = CC(I,3,K)-CC(I,4,K) 1299 TI3 = CC(I,3,K)+CC(I,4,K) 1300 TR5 = CC(I-1,2,K)-CC(I-1,5,K) 1301 TR2 = CC(I-1,2,K)+CC(I-1,5,K) 1302 TR4 = CC(I-1,3,K)-CC(I-1,4,K) 1303 TR3 = CC(I-1,3,K)+CC(I-1,4,K) 1304 CH(I-1,K,1) = CC(I-1,1,K)+TR2+TR3 1305 CH(I,K,1) = CC(I,1,K)+TI2+TI3 1306 CR2 = CC(I-1,1,K)+TR11*TR2+TR12*TR3 1307 CI2 = CC(I,1,K)+TR11*TI2+TR12*TI3 1308 CR3 = CC(I-1,1,K)+TR12*TR2+TR11*TR3 1309 CI3 = CC(I,1,K)+TR12*TI2+TR11*TI3 1310 CR5 = TI11*TR5+TI12*TR4 1311 CI5 = TI11*TI5+TI12*TI4 1312 CR4 = TI12*TR5-TI11*TR4 1313 CI4 = TI12*TI5-TI11*TI4 1314 DR3 = CR3-CI4 1315 DR4 = CR3+CI4 1316 DI3 = CI3+CR4 1317 DI4 = CI3-CR4 1318 DR5 = CR2+CI5 1319 DR2 = CR2-CI5 1320 DI5 = CI2-CR5 1321 DI2 = CI2+CR5 1322 CH(I-1,K,2) = WA1(I-1)*DR2-WA1(I)*DI2 1323 CH(I,K,2) = WA1(I-1)*DI2+WA1(I)*DR2 1324 CH(I-1,K,3) = WA2(I-1)*DR3-WA2(I)*DI3 1325 CH(I,K,3) = WA2(I-1)*DI3+WA2(I)*DR3 1326 CH(I-1,K,4) = WA3(I-1)*DR4-WA3(I)*DI4 1327 CH(I,K,4) = WA3(I-1)*DI4+WA3(I)*DR4 1328 CH(I-1,K,5) = WA4(I-1)*DR5-WA4(I)*DI5 1329 CH(I,K,5) = WA4(I-1)*DI5+WA4(I)*DR5 1330 106 CONTINUE 1331 107 CONTINUE 1332 RETURN 1333 END 1334