 |
PUMA
219
Portable University Model of the Atmosphere
|
|
PUMA
219
Portable University Model of the Atmosphere
|
00001 ! ============= 00002 ! MODULE FFTMOD 00003 ! ============= 00004 00005 ! alternate module with factors 2, 3, 4, 5, 6, 8 00006 ! this is a FORTRAN-90 version of the old FFT991 package 00007 ! use this for resolutions not supported by the <fftmod>, 00008 ! e.g. T63 (NLAT=96), T106 (NLAT=160) 00009 00010 module fftmod 00011 parameter(NRES = 8) 00012 integer :: nallowed(NRES) = (/ 64, 96, 128, 192, 256, 320, 384, 512 /) 00013 ! T21 - N64 : 8-4-2 00014 ! T31 - N96 : 8-4-3 00015 ! T42 - N128 : 8-4-4 00016 ! T63 - N192 : 8-6-4 00017 ! T85 - N256 : 8-4-4-2 00018 ! T106 - N320 : 8-5-4-2 00019 ! T127 - N384 : 8-4-4-3 00020 ! T170 - N512 : 8-4-4-4 00021 integer :: lastn = 0 00022 integer :: ifax(10) 00023 real,allocatable :: trigs(:) 00024 end module fftmod 00025 00026 00027 ! ================ 00028 ! SUBROUTINE GP2FC 00029 ! ================ 00030 00031 subroutine gp2fc(a,n,lot) 00032 use fftmod 00033 real :: a(n,lot) 00034 real :: worka(n+2,lot) 00035 real :: workb(n+2,lot) 00036 00037 if (n /= lastn) then 00038 if (allocated(trigs)) deallocate(trigs) 00039 allocate(trigs(n)) 00040 lastn = n 00041 call fftini(n) 00042 endif 00043 00044 worka(1:n,:) = a(:,:) 00045 worka(n+1:n+2,:) = 0.0 00046 00047 call fft991(worka,workb,trigs,ifax,1,n+2,n,lot,-1) 00048 00049 a(:,:) = worka(1:n,:) 00050 00051 return 00052 end subroutine gp2fc 00053 00054 00055 ! ================ 00056 ! SUBROUTINE FC2GP 00057 ! ================ 00058 00059 subroutine fc2gp(a,n,lot) 00060 use fftmod 00061 real :: a(n,lot) 00062 real :: worka(n+2,lot) 00063 real :: workb(n+2,lot) 00064 00065 if (n /= lastn) then 00066 if (allocated(trigs)) deallocate(trigs) 00067 allocate(trigs(n)) 00068 lastn = n 00069 call fftini(n) 00070 endif 00071 00072 worka(1:n,:) = a(:,:) 00073 worka(n+1:n+2,:) = 0.0 00074 00075 call fft991(worka,workb,trigs,ifax,1,n+2,n,lot,1) 00076 00077 a(:,:) = worka(1:n,:) 00078 00079 end subroutine fc2gp 00080 00081 00082 ! ================= 00083 ! SUBROUTINE FFTINI 00084 ! ================= 00085 00086 subroutine fftini(n) 00087 use fftmod 00088 logical labort 00089 00090 ! check for allowed values of n 00091 00092 labort = .true. 00093 do j = 1 , NRES 00094 if (n == nallowed(j)) labort = .false. 00095 enddo 00096 00097 if (labort) then 00098 write (*,*) '*** FFT does not support n = ',n,' ***' 00099 write (*,*) 'Following resolutions may be used:' 00100 write (*,*) '----------------------------------' 00101 do j = 1 , NRES 00102 write (*,1000) nallowed(j), nallowed(j)/2, nallowed(j)/3 00103 enddo 00104 stop 00105 endif 00106 1000 format(' NLON=',I5,' NLAT=',I5,' NTRU=',I5) 00107 00108 call set99(trigs,ifax,n) ! Factorization of n and sine/cosine values 00109 00110 return 00111 end subroutine fftini 00112 00113 00114 00115 00116 SUBROUTINE FFT991(A,WORK,TRIGS,IFAX,INC,JUMP,N,LOT,ISIGN) 00117 ! SUBROUTINE FFT991 - MULTIPLE FAST REAL PERIODIC TRANSFORM 00118 ! SUPERSEDES PREVIOUS ROUTINE 'FFT991' 00119 ! 00120 ! REAL TRANSFORM OF LENGTH N PERFORMED BY REMOVING REDUNDANT 00121 ! OPERATIONS FROM COMPLEX TRANSFORM OF LENGTH N 00122 ! 00123 ! A IS THE ARRAY CONTAINING INPUT & OUTPUT DATA 00124 ! WORK IS AN AREA OF SIZE (N+1)*MIN(LOT,64) 00125 ! TRIGS IS A PREVIOUSLY PREPARED LIST OF TRIG FUNCTION VALUES 00126 ! IFAX IS A PREVIOUSLY PREPARED LIST OF FACTORS OF N 00127 ! INC IS THE INCREMENT WITHIN EACH DATA 'VECTOR' 00128 ! (E.G. INC=1 FOR CONSECUTIVELY STORED DATA) 00129 ! JUMP IS THE INCREMENT BETWEEN THE START OF EACH DATA VECTOR 00130 ! N IS THE LENGTH OF THE DATA VECTORS 00131 ! LOT IS THE NUMBER OF DATA VECTORS 00132 ! ISIGN = +1 FOR TRANSFORM FROM SPECTRAL TO GRIDPOINT 00133 ! = -1 FOR TRANSFORM FROM GRIDPOINT TO SPECTRAL 00134 ! 00135 ! ORDERING OF COEFFICIENTS: 00136 ! A(0),B(0),A(1),B(1),A(2),B(2),...,A(N/2),B(N/2) 00137 ! WHERE B(0)=B(N/2)=0; (N+2) LOCATIONS REQUIRED 00138 ! 00139 ! ORDERING OF DATA: 00140 ! X(0),X(1),X(2),...,X(N-1), 0 , 0 ; (N+2) LOCATIONS REQUIRED 00141 ! 00142 ! VECTORIZATION IS ACHIEVED ON CRAY BY DOING THE TRANSFORMS 00143 ! IN PARALLEL 00144 ! 00145 ! N MUST BE COMPOSED OF FACTORS 2,3 & 5 BUT DOES NOT HAVE TO BE EVEN 00146 ! 00147 ! DEFINITION OF TRANSFORMS: 00148 ! ------------------------- 00149 ! 00150 ! ISIGN=+1: X(J)=SUM(K=0,...,N-1)(C(K)*EXP(2*I*J*K*PI/N)) 00151 ! WHERE C(K)=A(K)+I*B(K) AND C(N-K)=A(K)-I*B(K) 00152 ! 00153 ! ISIGN=-1: A(K)=(1/N)*SUM(J=0,...,N-1)(X(J)*COS(2*J*K*PI/N)) 00154 ! B(K)=-(1/N)*SUM(J=0,...,N-1)(X(J)*SIN(2*J*K*PI/N)) 00155 ! 00156 DIMENSION A(*),WORK(*),TRIGS(N),IFAX(10) 00157 ! 00158 NFAX=IFAX(1) 00159 NX=N+1 00160 IF (MOD(N,2).EQ.1) NX=N 00161 NBLOX=1+(LOT-1)/64 00162 NVEX=LOT-(NBLOX-1)*64 00163 IF (ISIGN.EQ.-1) GO TO 300 00164 ! 00165 ! ISIGN=+1, SPECTRAL TO GRIDPOINT TRANSFORM 00166 ! ----------------------------------------- 00167 100 CONTINUE 00168 ISTART=1 00169 DO 220 NB=1,NBLOX 00170 IA=ISTART 00171 I=ISTART 00172 DO 110 J=1,NVEX 00173 A(I+INC)=0.5*A(I) 00174 I=I+JUMP 00175 110 CONTINUE 00176 IF (MOD(N,2).EQ.1) GO TO 130 00177 I=ISTART+N*INC 00178 DO 120 J=1,NVEX 00179 A(I)=0.5*A(I) 00180 I=I+JUMP 00181 120 CONTINUE 00182 130 CONTINUE 00183 IA=ISTART+INC 00184 LA=1 00185 IGO=+1 00186 ! 00187 DO 160 K=1,NFAX 00188 IFAC=IFAX(K+1) 00189 IERR=-1 00190 IF (IGO.EQ.-1) GO TO 140 00191 CALL RPASSM(A(IA),A(IA+LA*INC),WORK(1),WORK(IFAC*LA+1),TRIGS, & 00192 INC,1,JUMP,NX,NVEX,N,IFAC,LA,IERR) 00193 GO TO 150 00194 140 CONTINUE 00195 CALL RPASSM(WORK(1),WORK(LA+1),A(IA),A(IA+IFAC*LA*INC),TRIGS, & 00196 1,INC,NX,JUMP,NVEX,N,IFAC,LA,IERR) 00197 150 CONTINUE 00198 IF (IERR.NE.0) GO TO 500 00199 LA=IFAC*LA 00200 IGO=-IGO 00201 IA=ISTART 00202 160 CONTINUE 00203 ! 00204 ! IF NECESSARY, COPY RESULTS BACK TO A 00205 ! ------------------------------------ 00206 IF (MOD(NFAX,2).EQ.0) GO TO 190 00207 IBASE=1 00208 JBASE=IA 00209 DO 180 JJ=1,NVEX 00210 I=IBASE 00211 J=JBASE 00212 DO 170 II=1,N 00213 A(J)=WORK(I) 00214 I=I+1 00215 J=J+INC 00216 170 CONTINUE 00217 IBASE=IBASE+NX 00218 JBASE=JBASE+JUMP 00219 180 CONTINUE 00220 190 CONTINUE 00221 ! 00222 ! FILL IN ZEROS AT END 00223 ! -------------------- 00224 IX=ISTART+N*INC 00225 DO 210 J=1,NVEX 00226 A(IX)=0.0 00227 A(IX+INC)=0.0 00228 IX=IX+JUMP 00229 210 CONTINUE 00230 ! 00231 ISTART=ISTART+NVEX*JUMP 00232 NVEX=64 00233 220 CONTINUE 00234 RETURN 00235 ! 00236 ! ISIGN=-1, GRIDPOINT TO SPECTRAL TRANSFORM 00237 ! ----------------------------------------- 00238 300 CONTINUE 00239 ISTART=1 00240 DO 410 NB=1,NBLOX 00241 IA=ISTART 00242 LA=N 00243 IGO=+1 00244 ! 00245 DO 340 K=1,NFAX 00246 IFAC=IFAX(NFAX+2-K) 00247 LA=LA/IFAC 00248 IERR=-1 00249 IF (IGO.EQ.-1) GO TO 320 00250 CALL QPASSM(A(IA),A(IA+IFAC*LA*INC),WORK(1),WORK(LA+1),TRIGS, & 00251 INC,1,JUMP,NX,NVEX,N,IFAC,LA,IERR) 00252 GO TO 330 00253 320 CONTINUE 00254 CALL QPASSM(WORK(1),WORK(IFAC*LA+1),A(IA),A(IA+LA*INC),TRIGS, & 00255 1,INC,NX,JUMP,NVEX,N,IFAC,LA,IERR) 00256 330 CONTINUE 00257 IF (IERR.NE.0) GO TO 500 00258 IGO=-IGO 00259 IA=ISTART+INC 00260 340 CONTINUE 00261 ! 00262 ! IF NECESSARY, COPY RESULTS BACK TO A 00263 ! ------------------------------------ 00264 IF (MOD(NFAX,2).EQ.0) GO TO 370 00265 IBASE=1 00266 JBASE=IA 00267 DO 360 JJ=1,NVEX 00268 I=IBASE 00269 J=JBASE 00270 DO 350 II=1,N 00271 A(J)=WORK(I) 00272 I=I+1 00273 J=J+INC 00274 350 CONTINUE 00275 IBASE=IBASE+NX 00276 JBASE=JBASE+JUMP 00277 360 CONTINUE 00278 370 CONTINUE 00279 ! 00280 ! SHIFT A(0) & FILL IN ZERO IMAG PARTS 00281 ! ------------------------------------ 00282 IX=ISTART 00283 DO 380 J=1,NVEX 00284 A(IX)=A(IX+INC) 00285 A(IX+INC)=0.0 00286 IX=IX+JUMP 00287 380 CONTINUE 00288 IF (MOD(N,2).EQ.1) GO TO 400 00289 IZ=ISTART+(N+1)*INC 00290 DO 390 J=1,NVEX 00291 A(IZ)=0.0 00292 IZ=IZ+JUMP 00293 390 CONTINUE 00294 400 CONTINUE 00295 ! 00296 ISTART=ISTART+NVEX*JUMP 00297 NVEX=64 00298 410 CONTINUE 00299 RETURN 00300 ! 00301 ! ERROR MESSAGES 00302 ! -------------- 00303 500 CONTINUE 00304 GO TO (510,530,550) IERR 00305 510 CONTINUE 00306 WRITE(6,520) NVEX 00307 520 FORMAT(16H1VECTOR LENGTH =,I4,17H, GREATER THAN 64) 00308 GO TO 570 00309 530 CONTINUE 00310 WRITE(6,540) IFAC 00311 540 FORMAT( 9H1FACTOR =,I3,17H, NOT CATERED FOR) 00312 GO TO 570 00313 550 CONTINUE 00314 WRITE(6,560) IFAC 00315 560 FORMAT(9H1FACTOR =,I3,31H, ONLY CATERED FOR IF LA*IFAC=N) 00316 570 CONTINUE 00317 RETURN 00318 END 00319 00320 00321 00322 SUBROUTINE SET99(TRIGS,IFAX,N) 00323 ! SUBROUTINE SET99 - COMPUTES FACTORS OF N & TRIGONOMETRIC 00324 ! FUNCTINS REQUIRED BY FFT99 & FFT991 00325 ! 00326 DIMENSION TRIGS(N),IFAX(10),JFAX(10),LFAX(7) 00327 DATA LFAX/6,8,5,4,3,2,1/ 00328 IXXX=1 00329 ! 00330 DEL=4.0*ASIN(1.0)/FLOAT(N) 00331 NIL=0 00332 NHL=(N/2)-1 00333 DO 10 K=NIL,NHL 00334 ANGLE=FLOAT(K)*DEL 00335 TRIGS(2*K+1)=COS(ANGLE) 00336 TRIGS(2*K+2)=SIN(ANGLE) 00337 10 CONTINUE 00338 ! 00339 ! FIND FACTORS OF N (8,6,5,4,3,2; ONLY ONE 8 ALLOWED) 00340 ! LOOK FOR SIXES FIRST, STORE FACTORS IN DESCENDING ORDER 00341 NU=N 00342 IFAC=6 00343 K=0 00344 L=1 00345 20 CONTINUE 00346 IF (MOD(NU,IFAC).NE.0) GO TO 30 00347 K=K+1 00348 JFAX(K)=IFAC 00349 IF (IFAC.NE.8) GO TO 25 00350 IF (K.EQ.1) GO TO 25 00351 JFAX(1)=8 00352 JFAX(K)=6 00353 25 CONTINUE 00354 NU=NU/IFAC 00355 IF (NU.EQ.1) GO TO 50 00356 IF (IFAC.NE.8) GO TO 20 00357 30 CONTINUE 00358 L=L+1 00359 IFAC=LFAX(L) 00360 IF (IFAC.GT.1) GO TO 20 00361 ! 00362 WRITE(6,40) N 00363 40 FORMAT(4H1N =,I4,27H - CONTAINS ILLEGAL FACTORS) 00364 RETURN 00365 ! 00366 ! NOW REVERSE ORDER OF FACTORS 00367 50 CONTINUE 00368 NFAX=K 00369 IFAX(1)=NFAX 00370 DO 60 I=1,NFAX 00371 IFAX(NFAX+2-I)=JFAX(I) 00372 60 CONTINUE 00373 IFAX(10)=N 00374 RETURN 00375 END 00376 00377 00378 00379 00380 SUBROUTINE QPASSM(A,B,C,D,TRIGS,INC1,INC2,INC3,INC4,LOT,N,IFAC,LA,IERR) 00381 ! SUBROUTINE QPASSM - PERFORMS ONE PASS THROUGH DATA AS PART 00382 ! OF MULTIPLE REAL FFT (FOURIER ANALYSIS) ROUTINE 00383 ! 00384 ! A IS FIRST REAL INPUT VECTOR 00385 ! EQUIVALENCE B(1) WITH A(IFAC*LA*INC1+1) 00386 ! C IS FIRST REAL OUTPUT VECTOR 00387 ! EQUIVALENCE D(1) WITH C(LA*INC2+1) 00388 ! TRIGS IS A PRECALCULATED LIST OF SINES & COSINES 00389 ! INC1 IS THE ADDRESSING INCREMENT FOR A 00390 ! INC2 IS THE ADDRESSING INCREMENT FOR C 00391 ! INC3 IS THE INCREMENT BETWEEN INPUT VECTORS A 00392 ! INC4 IS THE INCREMENT BETWEEN OUTPUT VECTORS C 00393 ! LOT IS THE NUMBER OF VECTORS 00394 ! N IS THE LENGTH OF THE VECTORS 00395 ! IFAC IS THE CURRENT FACTOR OF N 00396 ! LA = N/(PRODUCT OF FACTORS USED SO FAR) 00397 ! IERR IS AN ERROR INDICATOR: 00398 ! 0 - PASS COMPLETED WITHOUT ERROR 00399 ! 1 - LOT GREATER THAN 64 00400 ! 2 - IFAC NOT CATERED FOR 00401 ! 3 - IFAC ONLY CATERED FOR IF LA=N/IFAC 00402 ! 00403 !----------------------------------------------------------------------- 00404 ! 00405 DIMENSION A(*),B(*),C(*),D(*),TRIGS(N) 00406 ! 00407 parameter(SIN36 = 0.587785252292473) 00408 parameter(SIN72 = 0.951056516295154) 00409 parameter(QRT5 = 0.559016994374947) 00410 parameter(SIN60 = 0.866025403784437) 00411 ! 00412 M=N/IFAC 00413 IINK=LA*INC1 00414 JINK=LA*INC2 00415 IJUMP=(IFAC-1)*IINK 00416 KSTOP=(N-IFAC)/(2*IFAC) 00417 ! 00418 IBAD=1 00419 IF (LOT.GT.64) GO TO 910 00420 IBASE=0 00421 JBASE=0 00422 IGO=IFAC-1 00423 IF (IGO.EQ.7) IGO=6 00424 IBAD=2 00425 IF (IGO.LT.1.OR.IGO.GT.6) GO TO 910 00426 GO TO (200,300,400,500,600,800),IGO 00427 ! 00428 ! CODING FOR FACTOR 2 00429 ! ------------------- 00430 200 CONTINUE 00431 IA=1 00432 IB=IA+IINK 00433 JA=1 00434 JB=JA+(2*M-LA)*INC2 00435 ! 00436 IF (LA.EQ.M) GO TO 290 00437 ! 00438 DO 220 L=1,LA 00439 I=IBASE 00440 J=JBASE 00441 !DIR$ IVDEP 00442 DO 210 IJK=1,LOT 00443 C(JA+J)=A(IA+I)+A(IB+I) 00444 C(JB+J)=A(IA+I)-A(IB+I) 00445 I=I+INC3 00446 J=J+INC4 00447 210 CONTINUE 00448 IBASE=IBASE+INC1 00449 JBASE=JBASE+INC2 00450 220 CONTINUE 00451 JA=JA+JINK 00452 JINK=2*JINK 00453 JB=JB-JINK 00454 IBASE=IBASE+IJUMP 00455 IJUMP=2*IJUMP+IINK 00456 IF (JA.EQ.JB) GO TO 260 00457 DO 250 K=LA,KSTOP,LA 00458 KB=K+K 00459 C1=TRIGS(KB+1) 00460 S1=TRIGS(KB+2) 00461 JBASE=0 00462 DO 240 L=1,LA 00463 I=IBASE 00464 J=JBASE 00465 !DIR$ IVDEP 00466 DO 230 IJK=1,LOT 00467 C(JA+J)=A(IA+I)+(C1*A(IB+I)+S1*B(IB+I)) 00468 C(JB+J)=A(IA+I)-(C1*A(IB+I)+S1*B(IB+I)) 00469 D(JA+J)=(C1*B(IB+I)-S1*A(IB+I))+B(IA+I) 00470 D(JB+J)=(C1*B(IB+I)-S1*A(IB+I))-B(IA+I) 00471 I=I+INC3 00472 J=J+INC4 00473 230 CONTINUE 00474 IBASE=IBASE+INC1 00475 JBASE=JBASE+INC2 00476 240 CONTINUE 00477 IBASE=IBASE+IJUMP 00478 JA=JA+JINK 00479 JB=JB-JINK 00480 250 CONTINUE 00481 IF (JA.GT.JB) GO TO 900 00482 260 CONTINUE 00483 JBASE=0 00484 DO 280 L=1,LA 00485 I=IBASE 00486 J=JBASE 00487 !DIR$ IVDEP 00488 DO 270 IJK=1,LOT 00489 C(JA+J)=A(IA+I) 00490 D(JA+J)=-A(IB+I) 00491 I=I+INC3 00492 J=J+INC4 00493 270 CONTINUE 00494 IBASE=IBASE+INC1 00495 JBASE=JBASE+INC2 00496 280 CONTINUE 00497 GO TO 900 00498 ! 00499 290 CONTINUE 00500 Z=1.0/FLOAT(N) 00501 DO 294 L=1,LA 00502 I=IBASE 00503 J=JBASE 00504 !DIR$ IVDEP 00505 DO 292 IJK=1,LOT 00506 C(JA+J)=Z*(A(IA+I)+A(IB+I)) 00507 C(JB+J)=Z*(A(IA+I)-A(IB+I)) 00508 I=I+INC3 00509 J=J+INC4 00510 292 CONTINUE 00511 IBASE=IBASE+INC1 00512 JBASE=JBASE+INC2 00513 294 CONTINUE 00514 GO TO 900 00515 ! 00516 ! CODING FOR FACTOR 3 00517 ! ------------------- 00518 300 CONTINUE 00519 IA=1 00520 IB=IA+IINK 00521 IC=IB+IINK 00522 JA=1 00523 JB=JA+(2*M-LA)*INC2 00524 JC=JB 00525 ! 00526 IF (LA.EQ.M) GO TO 390 00527 ! 00528 DO 320 L=1,LA 00529 I=IBASE 00530 J=JBASE 00531 !DIR$ IVDEP 00532 DO 310 IJK=1,LOT 00533 C(JA+J)=A(IA+I)+(A(IB+I)+A(IC+I)) 00534 C(JB+J)=A(IA+I)-0.5*(A(IB+I)+A(IC+I)) 00535 D(JB+J)=SIN60*(A(IC+I)-A(IB+I)) 00536 I=I+INC3 00537 J=J+INC4 00538 310 CONTINUE 00539 IBASE=IBASE+INC1 00540 JBASE=JBASE+INC2 00541 320 CONTINUE 00542 JA=JA+JINK 00543 JINK=2*JINK 00544 JB=JB+JINK 00545 JC=JC-JINK 00546 IBASE=IBASE+IJUMP 00547 IJUMP=2*IJUMP+IINK 00548 IF (JA.EQ.JC) GO TO 360 00549 DO 350 K=LA,KSTOP,LA 00550 KB=K+K 00551 KC=KB+KB 00552 C1=TRIGS(KB+1) 00553 S1=TRIGS(KB+2) 00554 C2=TRIGS(KC+1) 00555 S2=TRIGS(KC+2) 00556 JBASE=0 00557 DO 340 L=1,LA 00558 I=IBASE 00559 J=JBASE 00560 !DIR$ IVDEP 00561 DO 330 IJK=1,LOT 00562 A1=(C1*A(IB+I)+S1*B(IB+I))+(C2*A(IC+I)+S2*B(IC+I)) 00563 B1=(C1*B(IB+I)-S1*A(IB+I))+(C2*B(IC+I)-S2*A(IC+I)) 00564 A2=A(IA+I)-0.5*A1 00565 B2=B(IA+I)-0.5*B1 00566 A3=SIN60*((C1*A(IB+I)+S1*B(IB+I))-(C2*A(IC+I)+S2*B(IC+I))) 00567 B3=SIN60*((C1*B(IB+I)-S1*A(IB+I))-(C2*B(IC+I)-S2*A(IC+I))) 00568 C(JA+J)=A(IA+I)+A1 00569 D(JA+J)=B(IA+I)+B1 00570 C(JB+J)=A2+B3 00571 D(JB+J)=B2-A3 00572 C(JC+J)=A2-B3 00573 D(JC+J)=-(B2+A3) 00574 I=I+INC3 00575 J=J+INC4 00576 330 CONTINUE 00577 IBASE=IBASE+INC1 00578 JBASE=JBASE+INC2 00579 340 CONTINUE 00580 IBASE=IBASE+IJUMP 00581 JA=JA+JINK 00582 JB=JB+JINK 00583 JC=JC-JINK 00584 350 CONTINUE 00585 IF (JA.GT.JC) GO TO 900 00586 360 CONTINUE 00587 JBASE=0 00588 DO 380 L=1,LA 00589 I=IBASE 00590 J=JBASE 00591 !DIR$ IVDEP 00592 DO 370 IJK=1,LOT 00593 C(JA+J)=A(IA+I)+0.5*(A(IB+I)-A(IC+I)) 00594 D(JA+J)=-SIN60*(A(IB+I)+A(IC+I)) 00595 C(JB+J)=A(IA+I)-(A(IB+I)-A(IC+I)) 00596 I=I+INC3 00597 J=J+INC4 00598 370 CONTINUE 00599 IBASE=IBASE+INC1 00600 JBASE=JBASE+INC2 00601 380 CONTINUE 00602 GO TO 900 00603 ! 00604 390 CONTINUE 00605 Z=1.0/FLOAT(N) 00606 ZSIN60=Z*SIN60 00607 DO 394 L=1,LA 00608 I=IBASE 00609 J=JBASE 00610 !DIR$ IVDEP 00611 DO 392 IJK=1,LOT 00612 C(JA+J)=Z*(A(IA+I)+(A(IB+I)+A(IC+I))) 00613 C(JB+J)=Z*(A(IA+I)-0.5*(A(IB+I)+A(IC+I))) 00614 D(JB+J)=ZSIN60*(A(IC+I)-A(IB+I)) 00615 I=I+INC3 00616 J=J+INC4 00617 392 CONTINUE 00618 IBASE=IBASE+INC1 00619 JBASE=JBASE+INC2 00620 394 CONTINUE 00621 GO TO 900 00622 ! 00623 ! CODING FOR FACTOR 4 00624 ! ------------------- 00625 400 CONTINUE 00626 IA=1 00627 IB=IA+IINK 00628 IC=IB+IINK 00629 ID=IC+IINK 00630 JA=1 00631 JB=JA+(2*M-LA)*INC2 00632 JC=JB+2*M*INC2 00633 JD=JB 00634 ! 00635 IF (LA.EQ.M) GO TO 490 00636 ! 00637 DO 420 L=1,LA 00638 I=IBASE 00639 J=JBASE 00640 !DIR$ IVDEP 00641 DO 410 IJK=1,LOT 00642 C(JA+J)=(A(IA+I)+A(IC+I))+(A(IB+I)+A(ID+I)) 00643 C(JC+J)=(A(IA+I)+A(IC+I))-(A(IB+I)+A(ID+I)) 00644 C(JB+J)=A(IA+I)-A(IC+I) 00645 D(JB+J)=A(ID+I)-A(IB+I) 00646 I=I+INC3 00647 J=J+INC4 00648 410 CONTINUE 00649 IBASE=IBASE+INC1 00650 JBASE=JBASE+INC2 00651 420 CONTINUE 00652 JA=JA+JINK 00653 JINK=2*JINK 00654 JB=JB+JINK 00655 JC=JC-JINK 00656 JD=JD-JINK 00657 IBASE=IBASE+IJUMP 00658 IJUMP=2*IJUMP+IINK 00659 IF (JB.EQ.JC) GO TO 460 00660 DO 450 K=LA,KSTOP,LA 00661 KB=K+K 00662 KC=KB+KB 00663 KD=KC+KB 00664 C1=TRIGS(KB+1) 00665 S1=TRIGS(KB+2) 00666 C2=TRIGS(KC+1) 00667 S2=TRIGS(KC+2) 00668 C3=TRIGS(KD+1) 00669 S3=TRIGS(KD+2) 00670 JBASE=0 00671 DO 440 L=1,LA 00672 I=IBASE 00673 J=JBASE 00674 !DIR$ IVDEP 00675 DO 430 IJK=1,LOT 00676 A0=A(IA+I)+(C2*A(IC+I)+S2*B(IC+I)) 00677 A2=A(IA+I)-(C2*A(IC+I)+S2*B(IC+I)) 00678 A1=(C1*A(IB+I)+S1*B(IB+I))+(C3*A(ID+I)+S3*B(ID+I)) 00679 A3=(C1*A(IB+I)+S1*B(IB+I))-(C3*A(ID+I)+S3*B(ID+I)) 00680 B0=B(IA+I)+(C2*B(IC+I)-S2*A(IC+I)) 00681 B2=B(IA+I)-(C2*B(IC+I)-S2*A(IC+I)) 00682 B1=(C1*B(IB+I)-S1*A(IB+I))+(C3*B(ID+I)-S3*A(ID+I)) 00683 B3=(C1*B(IB+I)-S1*A(IB+I))-(C3*B(ID+I)-S3*A(ID+I)) 00684 C(JA+J)=A0+A1 00685 C(JC+J)=A0-A1 00686 D(JA+J)=B0+B1 00687 D(JC+J)=B1-B0 00688 C(JB+J)=A2+B3 00689 C(JD+J)=A2-B3 00690 D(JB+J)=B2-A3 00691 D(JD+J)=-(B2+A3) 00692 I=I+INC3 00693 J=J+INC4 00694 430 CONTINUE 00695 IBASE=IBASE+INC1 00696 JBASE=JBASE+INC2 00697 440 CONTINUE 00698 IBASE=IBASE+IJUMP 00699 JA=JA+JINK 00700 JB=JB+JINK 00701 JC=JC-JINK 00702 JD=JD-JINK 00703 450 CONTINUE 00704 IF (JB.GT.JC) GO TO 900 00705 460 CONTINUE 00706 SIN45=SQRT(0.5) 00707 JBASE=0 00708 DO 480 L=1,LA 00709 I=IBASE 00710 J=JBASE 00711 !DIR$ IVDEP 00712 DO 470 IJK=1,LOT 00713 C(JA+J)=A(IA+I)+SIN45*(A(IB+I)-A(ID+I)) 00714 C(JB+J)=A(IA+I)-SIN45*(A(IB+I)-A(ID+I)) 00715 D(JA+J)=-A(IC+I)-SIN45*(A(IB+I)+A(ID+I)) 00716 D(JB+J)=A(IC+I)-SIN45*(A(IB+I)+A(ID+I)) 00717 I=I+INC3 00718 J=J+INC4 00719 470 CONTINUE 00720 IBASE=IBASE+INC1 00721 JBASE=JBASE+INC2 00722 480 CONTINUE 00723 GO TO 900 00724 ! 00725 490 CONTINUE 00726 Z=1.0/FLOAT(N) 00727 DO 494 L=1,LA 00728 I=IBASE 00729 J=JBASE 00730 !DIR$ IVDEP 00731 DO 492 IJK=1,LOT 00732 C(JA+J)=Z*((A(IA+I)+A(IC+I))+(A(IB+I)+A(ID+I))) 00733 C(JC+J)=Z*((A(IA+I)+A(IC+I))-(A(IB+I)+A(ID+I))) 00734 C(JB+J)=Z*(A(IA+I)-A(IC+I)) 00735 D(JB+J)=Z*(A(ID+I)-A(IB+I)) 00736 I=I+INC3 00737 J=J+INC4 00738 492 CONTINUE 00739 IBASE=IBASE+INC1 00740 JBASE=JBASE+INC2 00741 494 CONTINUE 00742 GO TO 900 00743 ! 00744 ! CODING FOR FACTOR 5 00745 ! ------------------- 00746 500 CONTINUE 00747 IA=1 00748 IB=IA+IINK 00749 IC=IB+IINK 00750 ID=IC+IINK 00751 IE=ID+IINK 00752 JA=1 00753 JB=JA+(2*M-LA)*INC2 00754 JC=JB+2*M*INC2 00755 JD=JC 00756 JE=JB 00757 ! 00758 IF (LA.EQ.M) GO TO 590 00759 ! 00760 DO 520 L=1,LA 00761 I=IBASE 00762 J=JBASE 00763 !DIR$ IVDEP 00764 DO 510 IJK=1,LOT 00765 A1=A(IB+I)+A(IE+I) 00766 A3=A(IB+I)-A(IE+I) 00767 A2=A(IC+I)+A(ID+I) 00768 A4=A(IC+I)-A(ID+I) 00769 A5=A(IA+I)-0.25*(A1+A2) 00770 A6=QRT5*(A1-A2) 00771 C(JA+J)=A(IA+I)+(A1+A2) 00772 C(JB+J)=A5+A6 00773 C(JC+J)=A5-A6 00774 D(JB+J)=-SIN72*A3-SIN36*A4 00775 D(JC+J)=-SIN36*A3+SIN72*A4 00776 I=I+INC3 00777 J=J+INC4 00778 510 CONTINUE 00779 IBASE=IBASE+INC1 00780 JBASE=JBASE+INC2 00781 520 CONTINUE 00782 JA=JA+JINK 00783 JINK=2*JINK 00784 JB=JB+JINK 00785 JC=JC+JINK 00786 JD=JD-JINK 00787 JE=JE-JINK 00788 IBASE=IBASE+IJUMP 00789 IJUMP=2*IJUMP+IINK 00790 IF (JB.EQ.JD) GO TO 560 00791 DO 550 K=LA,KSTOP,LA 00792 KB=K+K 00793 KC=KB+KB 00794 KD=KC+KB 00795 KE=KD+KB 00796 C1=TRIGS(KB+1) 00797 S1=TRIGS(KB+2) 00798 C2=TRIGS(KC+1) 00799 S2=TRIGS(KC+2) 00800 C3=TRIGS(KD+1) 00801 S3=TRIGS(KD+2) 00802 C4=TRIGS(KE+1) 00803 S4=TRIGS(KE+2) 00804 JBASE=0 00805 DO 540 L=1,LA 00806 I=IBASE 00807 J=JBASE 00808 !DIR$ IVDEP 00809 DO 530 IJK=1,LOT 00810 A1=(C1*A(IB+I)+S1*B(IB+I))+(C4*A(IE+I)+S4*B(IE+I)) 00811 A3=(C1*A(IB+I)+S1*B(IB+I))-(C4*A(IE+I)+S4*B(IE+I)) 00812 A2=(C2*A(IC+I)+S2*B(IC+I))+(C3*A(ID+I)+S3*B(ID+I)) 00813 A4=(C2*A(IC+I)+S2*B(IC+I))-(C3*A(ID+I)+S3*B(ID+I)) 00814 B1=(C1*B(IB+I)-S1*A(IB+I))+(C4*B(IE+I)-S4*A(IE+I)) 00815 B3=(C1*B(IB+I)-S1*A(IB+I))-(C4*B(IE+I)-S4*A(IE+I)) 00816 B2=(C2*B(IC+I)-S2*A(IC+I))+(C3*B(ID+I)-S3*A(ID+I)) 00817 B4=(C2*B(IC+I)-S2*A(IC+I))-(C3*B(ID+I)-S3*A(ID+I)) 00818 A5=A(IA+I)-0.25*(A1+A2) 00819 A6=QRT5*(A1-A2) 00820 B5=B(IA+I)-0.25*(B1+B2) 00821 B6=QRT5*(B1-B2) 00822 A10=A5+A6 00823 A20=A5-A6 00824 B10=B5+B6 00825 B20=B5-B6 00826 A11=SIN72*B3+SIN36*B4 00827 A21=SIN36*B3-SIN72*B4 00828 B11=SIN72*A3+SIN36*A4 00829 B21=SIN36*A3-SIN72*A4 00830 C(JA+J)=A(IA+I)+(A1+A2) 00831 C(JB+J)=A10+A11 00832 C(JE+J)=A10-A11 00833 C(JC+J)=A20+A21 00834 C(JD+J)=A20-A21 00835 D(JA+J)=B(IA+I)+(B1+B2) 00836 D(JB+J)=B10-B11 00837 D(JE+J)=-(B10+B11) 00838 D(JC+J)=B20-B21 00839 D(JD+J)=-(B20+B21) 00840 I=I+INC3 00841 J=J+INC4 00842 530 CONTINUE 00843 IBASE=IBASE+INC1 00844 JBASE=JBASE+INC2 00845 540 CONTINUE 00846 IBASE=IBASE+IJUMP 00847 JA=JA+JINK 00848 JB=JB+JINK 00849 JC=JC+JINK 00850 JD=JD-JINK 00851 JE=JE-JINK 00852 550 CONTINUE 00853 IF (JB.GT.JD) GO TO 900 00854 560 CONTINUE 00855 JBASE=0 00856 DO 580 L=1,LA 00857 I=IBASE 00858 J=JBASE 00859 !DIR$ IVDEP 00860 DO 570 IJK=1,LOT 00861 A1=A(IB+I)+A(IE+I) 00862 A3=A(IB+I)-A(IE+I) 00863 A2=A(IC+I)+A(ID+I) 00864 A4=A(IC+I)-A(ID+I) 00865 A5=A(IA+I)+0.25*(A3-A4) 00866 A6=QRT5*(A3+A4) 00867 C(JA+J)=A5+A6 00868 C(JB+J)=A5-A6 00869 C(JC+J)=A(IA+I)-(A3-A4) 00870 D(JA+J)=-SIN36*A1-SIN72*A2 00871 D(JB+J)=-SIN72*A1+SIN36*A2 00872 I=I+INC3 00873 J=J+INC4 00874 570 CONTINUE 00875 IBASE=IBASE+INC1 00876 JBASE=JBASE+INC2 00877 580 CONTINUE 00878 GO TO 900 00879 ! 00880 590 CONTINUE 00881 Z=1.0/FLOAT(N) 00882 ZQRT5=Z*QRT5 00883 ZSIN36=Z*SIN36 00884 ZSIN72=Z*SIN72 00885 DO 594 L=1,LA 00886 I=IBASE 00887 J=JBASE 00888 !DIR$ IVDEP 00889 DO 592 IJK=1,LOT 00890 A1=A(IB+I)+A(IE+I) 00891 A3=A(IB+I)-A(IE+I) 00892 A2=A(IC+I)+A(ID+I) 00893 A4=A(IC+I)-A(ID+I) 00894 A5=Z*(A(IA+I)-0.25*(A1+A2)) 00895 A6=ZQRT5*(A1-A2) 00896 C(JA+J)=Z*(A(IA+I)+(A1+A2)) 00897 C(JB+J)=A5+A6 00898 C(JC+J)=A5-A6 00899 D(JB+J)=-ZSIN72*A3-ZSIN36*A4 00900 D(JC+J)=-ZSIN36*A3+ZSIN72*A4 00901 I=I+INC3 00902 J=J+INC4 00903 592 CONTINUE 00904 IBASE=IBASE+INC1 00905 JBASE=JBASE+INC2 00906 594 CONTINUE 00907 GO TO 900 00908 ! 00909 ! CODING FOR FACTOR 6 00910 ! ------------------- 00911 600 CONTINUE 00912 IA=1 00913 IB=IA+IINK 00914 IC=IB+IINK 00915 ID=IC+IINK 00916 IE=ID+IINK 00917 IF=IE+IINK 00918 JA=1 00919 JB=JA+(2*M-LA)*INC2 00920 JC=JB+2*M*INC2 00921 JD=JC+2*M*INC2 00922 JE=JC 00923 JF=JB 00924 ! 00925 IF (LA.EQ.M) GO TO 690 00926 ! 00927 DO 620 L=1,LA 00928 I=IBASE 00929 J=JBASE 00930 !DIR$ IVDEP 00931 DO 610 IJK=1,LOT 00932 A11=(A(IC+I)+A(IF+I))+(A(IB+I)+A(IE+I)) 00933 C(JA+J)=(A(IA+I)+A(ID+I))+A11 00934 C(JC+J)=(A(IA+I)+A(ID+I)-0.5*A11) 00935 D(JC+J)=SIN60*((A(IC+I)+A(IF+I))-(A(IB+I)+A(IE+I))) 00936 A11=(A(IC+I)-A(IF+I))+(A(IE+I)-A(IB+I)) 00937 C(JB+J)=(A(IA+I)-A(ID+I))-0.5*A11 00938 D(JB+J)=SIN60*((A(IE+I)-A(IB+I))-(A(IC+I)-A(IF+I))) 00939 C(JD+J)=(A(IA+I)-A(ID+I))+A11 00940 I=I+INC3 00941 J=J+INC4 00942 610 CONTINUE 00943 IBASE=IBASE+INC1 00944 JBASE=JBASE+INC2 00945 620 CONTINUE 00946 JA=JA+JINK 00947 JINK=2*JINK 00948 JB=JB+JINK 00949 JC=JC+JINK 00950 JD=JD-JINK 00951 JE=JE-JINK 00952 JF=JF-JINK 00953 IBASE=IBASE+IJUMP 00954 IJUMP=2*IJUMP+IINK 00955 IF (JC.EQ.JD) GO TO 660 00956 DO 650 K=LA,KSTOP,LA 00957 KB=K+K 00958 KC=KB+KB 00959 KD=KC+KB 00960 KE=KD+KB 00961 KF=KE+KB 00962 C1=TRIGS(KB+1) 00963 S1=TRIGS(KB+2) 00964 C2=TRIGS(KC+1) 00965 S2=TRIGS(KC+2) 00966 C3=TRIGS(KD+1) 00967 S3=TRIGS(KD+2) 00968 C4=TRIGS(KE+1) 00969 S4=TRIGS(KE+2) 00970 C5=TRIGS(KF+1) 00971 S5=TRIGS(KF+2) 00972 JBASE=0 00973 DO 640 L=1,LA 00974 I=IBASE 00975 J=JBASE 00976 !DIR$ IVDEP 00977 DO 630 IJK=1,LOT 00978 A1=C1*A(IB+I)+S1*B(IB+I) 00979 B1=C1*B(IB+I)-S1*A(IB+I) 00980 A2=C2*A(IC+I)+S2*B(IC+I) 00981 B2=C2*B(IC+I)-S2*A(IC+I) 00982 A3=C3*A(ID+I)+S3*B(ID+I) 00983 B3=C3*B(ID+I)-S3*A(ID+I) 00984 A4=C4*A(IE+I)+S4*B(IE+I) 00985 B4=C4*B(IE+I)-S4*A(IE+I) 00986 A5=C5*A(IF+I)+S5*B(IF+I) 00987 B5=C5*B(IF+I)-S5*A(IF+I) 00988 A11=(A2+A5)+(A1+A4) 00989 A20=(A(IA+I)+A3)-0.5*A11 00990 A21=SIN60*((A2+A5)-(A1+A4)) 00991 B11=(B2+B5)+(B1+B4) 00992 B20=(B(IA+I)+B3)-0.5*B11 00993 B21=SIN60*((B2+B5)-(B1+B4)) 00994 C(JA+J)=(A(IA+I)+A3)+A11 00995 D(JA+J)=(B(IA+I)+B3)+B11 00996 C(JC+J)=A20-B21 00997 D(JC+J)=A21+B20 00998 C(JE+J)=A20+B21 00999 D(JE+J)=A21-B20 01000 A11=(A2-A5)+(A4-A1) 01001 A20=(A(IA+I)-A3)-0.5*A11 01002 A21=SIN60*((A4-A1)-(A2-A5)) 01003 B11=(B5-B2)-(B4-B1) 01004 B20=(B3-B(IA+I))-0.5*B11 01005 B21=SIN60*((B5-B2)+(B4-B1)) 01006 C(JB+J)=A20-B21 01007 D(JB+J)=A21-B20 01008 C(JD+J)=A11+(A(IA+I)-A3) 01009 D(JD+J)=B11+(B3-B(IA+I)) 01010 C(JF+J)=A20+B21 01011 D(JF+J)=A21+B20 01012 I=I+INC3 01013 J=J+INC4 01014 630 CONTINUE 01015 IBASE=IBASE+INC1 01016 JBASE=JBASE+INC2 01017 640 CONTINUE 01018 IBASE=IBASE+IJUMP 01019 JA=JA+JINK 01020 JB=JB+JINK 01021 JC=JC+JINK 01022 JD=JD-JINK 01023 JE=JE-JINK 01024 JF=JF-JINK 01025 650 CONTINUE 01026 IF (JC.GT.JD) GO TO 900 01027 660 CONTINUE 01028 JBASE=0 01029 DO 680 L=1,LA 01030 I=IBASE 01031 J=JBASE 01032 !DIR$ IVDEP 01033 DO 670 IJK=1,LOT 01034 C(JA+J)=(A(IA+I)+0.5*(A(IC+I)-A(IE+I)))+ SIN60*(A(IB+I)-A(IF+I)) 01035 D(JA+J)=-(A(ID+I)+0.5*(A(IB+I)+A(IF+I)))-SIN60*(A(IC+I)+A(IE+I)) 01036 C(JB+J)=A(IA+I)-(A(IC+I)-A(IE+I)) 01037 D(JB+J)=A(ID+I)-(A(IB+I)+A(IF+I)) 01038 C(JC+J)=(A(IA+I)+0.5*(A(IC+I)-A(IE+I)))-SIN60*(A(IB+I)-A(IF+I)) 01039 D(JC+J)=-(A(ID+I)+0.5*(A(IB+I)+A(IF+I)))+SIN60*(A(IC+I)+A(IE+I)) 01040 I=I+INC3 01041 J=J+INC4 01042 670 CONTINUE 01043 IBASE=IBASE+INC1 01044 JBASE=JBASE+INC2 01045 680 CONTINUE 01046 GO TO 900 01047 ! 01048 690 CONTINUE 01049 Z=1.0/FLOAT(N) 01050 ZSIN60=Z*SIN60 01051 DO 694 L=1,LA 01052 I=IBASE 01053 J=JBASE 01054 !DIR$ IVDEP 01055 DO 692 IJK=1,LOT 01056 A11=(A(IC+I)+A(IF+I))+(A(IB+I)+A(IE+I)) 01057 C(JA+J)=Z*((A(IA+I)+A(ID+I))+A11) 01058 C(JC+J)=Z*((A(IA+I)+A(ID+I))-0.5*A11) 01059 D(JC+J)=ZSIN60*((A(IC+I)+A(IF+I))-(A(IB+I)+A(IE+I))) 01060 A11=(A(IC+I)-A(IF+I))+(A(IE+I)-A(IB+I)) 01061 C(JB+J)=Z*((A(IA+I)-A(ID+I))-0.5*A11) 01062 D(JB+J)=ZSIN60*((A(IE+I)-A(IB+I))-(A(IC+I)-A(IF+I))) 01063 C(JD+J)=Z*((A(IA+I)-A(ID+I))+A11) 01064 I=I+INC3 01065 J=J+INC4 01066 692 CONTINUE 01067 IBASE=IBASE+INC1 01068 JBASE=JBASE+INC2 01069 694 CONTINUE 01070 GO TO 900 01071 ! 01072 ! CODING FOR FACTOR 8 01073 ! ------------------- 01074 800 CONTINUE 01075 IBAD=3 01076 IF (LA.NE.M) GO TO 910 01077 IA=1 01078 IB=IA+IINK 01079 IC=IB+IINK 01080 ID=IC+IINK 01081 IE=ID+IINK 01082 IF=IE+IINK 01083 IG=IF+IINK 01084 IH=IG+IINK 01085 JA=1 01086 JB=JA+LA*INC2 01087 JC=JB+2*M*INC2 01088 JD=JC+2*M*INC2 01089 JE=JD+2*M*INC2 01090 Z=1.0/FLOAT(N) 01091 ZSIN45=Z*SQRT(0.5) 01092 ! 01093 DO 820 L=1,LA 01094 I=IBASE 01095 J=JBASE 01096 !DIR$ IVDEP 01097 DO 810 IJK=1,LOT 01098 C(JA+J)=Z*(((A(IA+I)+A(IE+I))+(A(IC+I)+A(IG+I)))+ & 01099 ((A(ID+I)+A(IH+I))+(A(IB+I)+A(IF+I)))) 01100 C(JE+J)=Z*(((A(IA+I)+A(IE+I))+(A(IC+I)+A(IG+I)))- & 01101 ((A(ID+I)+A(IH+I))+(A(IB+I)+A(IF+I)))) 01102 C(JC+J)=Z*((A(IA+I)+A(IE+I))-(A(IC+I)+A(IG+I))) 01103 D(JC+J)=Z*((A(ID+I)+A(IH+I))-(A(IB+I)+A(IF+I))) 01104 C(JB+J)=Z*(A(IA+I)-A(IE+I)) & 01105 +ZSIN45*((A(IH+I)-A(ID+I))-(A(IF+I)-A(IB+I))) 01106 C(JD+J)=Z*(A(IA+I)-A(IE+I)) & 01107 -ZSIN45*((A(IH+I)-A(ID+I))-(A(IF+I)-A(IB+I))) 01108 D(JB+J)=ZSIN45*((A(IH+I)-A(ID+I))+(A(IF+I)-A(IB+I))) & 01109 +Z*(A(IG+I)-A(IC+I)) 01110 D(JD+J)=ZSIN45*((A(IH+I)-A(ID+I))+(A(IF+I)-A(IB+I))) & 01111 -Z*(A(IG+I)-A(IC+I)) 01112 I=I+INC3 01113 J=J+INC4 01114 810 CONTINUE 01115 IBASE=IBASE+INC1 01116 JBASE=JBASE+INC2 01117 820 CONTINUE 01118 ! 01119 ! RETURN 01120 ! ------ 01121 900 CONTINUE 01122 IBAD=0 01123 910 CONTINUE 01124 IERR=IBAD 01125 RETURN 01126 END 01127 01128 01129 01130 SUBROUTINE RPASSM(A,B,C,D,TRIGS,INC1,INC2,INC3,INC4,LOT,N,IFAC,LA,IERR) 01131 ! SUBROUTINE 'RPASSM' - PERFORMS ONE PASS THROUGH DATA AS PART 01132 ! OF MULTIPLE REAL FFT (FOURIER SYNTHESIS) ROUTINE 01133 ! 01134 ! A IS FIRST REAL INPUT VECTOR 01135 ! EQUIVALENCE B(1) WITH A (LA*INC1+1) 01136 ! C IS FIRST REAL OUTPUT VECTOR 01137 ! EQUIVALENCE D(1) WITH C(IFAC*LA*INC2+1) 01138 ! TRIGS IS A PRECALCULATED LIST OF SINES & COSINES 01139 ! INC1 IS THE ADDRESSING INCREMENT FOR A 01140 ! INC2 IS THE ADDRESSING INCREMENT FOR C 01141 ! INC3 IS THE INCREMENT BETWEEN INPUT VECTORS A 01142 ! INC4 IS THE INCREMENT BETWEEN OUTPUT VECTORS C 01143 ! LOT IS THE NUMBER OF VECTORS 01144 ! N IS THE LENGTH OF THE VECTORS 01145 ! IFAC IS THE CURRENT FACTOR OF N 01146 ! LA IS THE PRODUCT OF PREVIOUS FACTORS 01147 ! IERR IS AN ERROR INDICATOR: 01148 ! 0 - PASS COMPLETED WITHOUT ERROR 01149 ! 1 - LOT GREATER THAN 64 01150 ! 2 - IFAC NOT CATERED FOR 01151 ! 3 - IFAC ONLY CATERED FOR IF LA=N/IFAC 01152 ! 01153 !----------------------------------------------------------------------- 01154 ! 01155 DIMENSION A(*),B(*),C(*),D(*),TRIGS(N) 01156 ! 01157 DIMENSION A10(64),A11(64),A20(64),A21(64),B10(64),B11(64),B20(64),B21(64) 01158 parameter(SIN36 = 0.587785252292473) 01159 parameter(SIN72 = 0.951056516295154) 01160 parameter(QRT5 = 0.559016994374947) 01161 parameter(SIN60 = 0.866025403784437) 01162 ! 01163 M=N/IFAC 01164 IINK=LA*INC1 01165 JINK=LA*INC2 01166 JUMP=(IFAC-1)*JINK 01167 KSTOP=(N-IFAC)/(2*IFAC) 01168 ! 01169 IBAD=1 01170 IF (LOT.GT.64) GO TO 910 01171 IBASE=0 01172 JBASE=0 01173 IGO=IFAC-1 01174 IF (IGO.EQ.7) IGO=6 01175 IBAD=2 01176 IF (IGO.LT.1.OR.IGO.GT.6) GO TO 910 01177 GO TO (200,300,400,500,600,800),IGO 01178 ! 01179 ! CODING FOR FACTOR 2 01180 ! ------------------- 01181 200 CONTINUE 01182 IA=1 01183 IB=IA+(2*M-LA)*INC1 01184 JA=1 01185 JB=JA+JINK 01186 ! 01187 IF (LA.EQ.M) GO TO 290 01188 ! 01189 DO 220 L=1,LA 01190 I=IBASE 01191 J=JBASE 01192 !DIR$ IVDEP 01193 DO 210 IJK=1,LOT 01194 C(JA+J)=A(IA+I)+A(IB+I) 01195 C(JB+J)=A(IA+I)-A(IB+I) 01196 I=I+INC3 01197 J=J+INC4 01198 210 CONTINUE 01199 IBASE=IBASE+INC1 01200 JBASE=JBASE+INC2 01201 220 CONTINUE 01202 IA=IA+IINK 01203 IINK=2*IINK 01204 IB=IB-IINK 01205 IBASE=0 01206 JBASE=JBASE+JUMP 01207 JUMP=2*JUMP+JINK 01208 IF (IA.EQ.IB) GO TO 260 01209 DO 250 K=LA,KSTOP,LA 01210 KB=K+K 01211 C1=TRIGS(KB+1) 01212 S1=TRIGS(KB+2) 01213 IBASE=0 01214 DO 240 L=1,LA 01215 I=IBASE 01216 J=JBASE 01217 !DIR$ IVDEP 01218 DO 230 IJK=1,LOT 01219 C(JA+J)=A(IA+I)+A(IB+I) 01220 D(JA+J)=B(IA+I)-B(IB+I) 01221 C(JB+J)=C1*(A(IA+I)-A(IB+I))-S1*(B(IA+I)+B(IB+I)) 01222 D(JB+J)=S1*(A(IA+I)-A(IB+I))+C1*(B(IA+I)+B(IB+I)) 01223 I=I+INC3 01224 J=J+INC4 01225 230 CONTINUE 01226 IBASE=IBASE+INC1 01227 JBASE=JBASE+INC2 01228 240 CONTINUE 01229 IA=IA+IINK 01230 IB=IB-IINK 01231 JBASE=JBASE+JUMP 01232 250 CONTINUE 01233 IF (IA.GT.IB) GO TO 900 01234 260 CONTINUE 01235 IBASE=0 01236 DO 280 L=1,LA 01237 I=IBASE 01238 J=JBASE 01239 !DIR$ IVDEP 01240 DO 270 IJK=1,LOT 01241 C(JA+J)=A(IA+I) 01242 C(JB+J)=-B(IA+I) 01243 I=I+INC3 01244 J=J+INC4 01245 270 CONTINUE 01246 IBASE=IBASE+INC1 01247 JBASE=JBASE+INC2 01248 280 CONTINUE 01249 GO TO 900 01250 ! 01251 290 CONTINUE 01252 DO 294 L=1,LA 01253 I=IBASE 01254 J=JBASE 01255 !DIR$ IVDEP 01256 DO 292 IJK=1,LOT 01257 C(JA+J)=2.0*(A(IA+I)+A(IB+I)) 01258 C(JB+J)=2.0*(A(IA+I)-A(IB+I)) 01259 I=I+INC3 01260 J=J+INC4 01261 292 CONTINUE 01262 IBASE=IBASE+INC1 01263 JBASE=JBASE+INC2 01264 294 CONTINUE 01265 GO TO 900 01266 ! 01267 ! CODING FOR FACTOR 3 01268 ! ------------------- 01269 300 CONTINUE 01270 IA=1 01271 IB=IA+(2*M-LA)*INC1 01272 IC=IB 01273 JA=1 01274 JB=JA+JINK 01275 JC=JB+JINK 01276 ! 01277 IF (LA.EQ.M) GO TO 390 01278 ! 01279 DO 320 L=1,LA 01280 I=IBASE 01281 J=JBASE 01282 !DIR$ IVDEP 01283 DO 310 IJK=1,LOT 01284 C(JA+J)=A(IA+I)+A(IB+I) 01285 C(JB+J)=(A(IA+I)-0.5*A(IB+I))-(SIN60*(B(IB+I))) 01286 C(JC+J)=(A(IA+I)-0.5*A(IB+I))+(SIN60*(B(IB+I))) 01287 I=I+INC3 01288 J=J+INC4 01289 310 CONTINUE 01290 IBASE=IBASE+INC1 01291 JBASE=JBASE+INC2 01292 320 CONTINUE 01293 IA=IA+IINK 01294 IINK=2*IINK 01295 IB=IB+IINK 01296 IC=IC-IINK 01297 JBASE=JBASE+JUMP 01298 JUMP=2*JUMP+JINK 01299 IF (IA.EQ.IC) GO TO 360 01300 DO 350 K=LA,KSTOP,LA 01301 KB=K+K 01302 KC=KB+KB 01303 C1=TRIGS(KB+1) 01304 S1=TRIGS(KB+2) 01305 C2=TRIGS(KC+1) 01306 S2=TRIGS(KC+2) 01307 IBASE=0 01308 DO 340 L=1,LA 01309 I=IBASE 01310 J=JBASE 01311 !DIR$ IVDEP 01312 DO 330 IJK=1,LOT 01313 C(JA+J)=A(IA+I)+(A(IB+I)+A(IC+I)) 01314 D(JA+J)=B(IA+I)+(B(IB+I)-B(IC+I)) 01315 C(JB+J)=C1*((A(IA+I)-0.5*(A(IB+I)+A(IC+I)))-(SIN60*(B(IB+ & 01316 I)+B(IC+I)))) -S1*((B(IA+I)-0.5*(B(IB+I)-B(IC+I)))+(SIN60 & 01317 *(A(IB+I)-A(IC+I)))) 01318 D(JB+J)=S1*((A(IA+I)-0.5*(A(IB+I)+A(IC+I)))-(SIN60*(B(IB+ & 01319 I)+B(IC+I)))) +C1*((B(IA+I)-0.5*(B(IB+I)-B(IC+I)))+(SIN60 & 01320 *(A(IB+I)-A(IC+I)))) 01321 C(JC+J)=C2*((A(IA+I)-0.5*(A(IB+I)+A(IC+I)))+(SIN60*(B(IB+ & 01322 I)+B(IC+I)))) -S2*((B(IA+I)-0.5*(B(IB+I)-B(IC+I)))-(SIN60 & 01323 *(A(IB+I)-A(IC+I)))) 01324 D(JC+J)=S2*((A(IA+I)-0.5*(A(IB+I)+A(IC+I)))+(SIN60*(B(IB+ & 01325 I)+B(IC+I)))) +C2*((B(IA+I)-0.5*(B(IB+I)-B(IC+I)))-(SIN60 & 01326 *(A(IB+I)-A(IC+I)))) 01327 I=I+INC3 01328 J=J+INC4 01329 330 CONTINUE 01330 IBASE=IBASE+INC1 01331 JBASE=JBASE+INC2 01332 340 CONTINUE 01333 IA=IA+IINK 01334 IB=IB+IINK 01335 IC=IC-IINK 01336 JBASE=JBASE+JUMP 01337 350 CONTINUE 01338 IF (IA.GT.IC) GO TO 900 01339 360 CONTINUE 01340 IBASE=0 01341 DO 380 L=1,LA 01342 I=IBASE 01343 J=JBASE 01344 !DIR$ IVDEP 01345 DO 370 IJK=1,LOT 01346 C(JA+J)=A(IA+I)+A(IB+I) 01347 C(JB+J)=(0.5*A(IA+I)-A(IB+I))-(SIN60*B(IA+I)) 01348 C(JC+J)=-(0.5*A(IA+I)-A(IB+I))-(SIN60*B(IA+I)) 01349 I=I+INC3 01350 J=J+INC4 01351 370 CONTINUE 01352 IBASE=IBASE+INC1 01353 JBASE=JBASE+INC2 01354 380 CONTINUE 01355 GO TO 900 01356 ! 01357 390 CONTINUE 01358 SSIN60=2.0*SIN60 01359 DO 394 L=1,LA 01360 I=IBASE 01361 J=JBASE 01362 !DIR$ IVDEP 01363 DO 392 IJK=1,LOT 01364 C(JA+J)=2.0*(A(IA+I)+A(IB+I)) 01365 C(JB+J)=(2.0*A(IA+I)-A(IB+I))-(SSIN60*B(IB+I)) 01366 C(JC+J)=(2.0*A(IA+I)-A(IB+I))+(SSIN60*B(IB+I)) 01367 I=I+INC3 01368 J=J+INC4 01369 392 CONTINUE 01370 IBASE=IBASE+INC1 01371 JBASE=JBASE+INC2 01372 394 CONTINUE 01373 GO TO 900 01374 ! 01375 ! CODING FOR FACTOR 4 01376 ! ------------------- 01377 400 CONTINUE 01378 IA=1 01379 IB=IA+(2*M-LA)*INC1 01380 IC=IB+2*M*INC1 01381 ID=IB 01382 JA=1 01383 JB=JA+JINK 01384 JC=JB+JINK 01385 JD=JC+JINK 01386 ! 01387 IF (LA.EQ.M) GO TO 490 01388 ! 01389 DO 420 L=1,LA 01390 I=IBASE 01391 J=JBASE 01392 !DIR$ IVDEP 01393 DO 410 IJK=1,LOT 01394 C(JA+J)=(A(IA+I)+A(IC+I))+A(IB+I) 01395 C(JB+J)=(A(IA+I)-A(IC+I))-B(IB+I) 01396 C(JC+J)=(A(IA+I)+A(IC+I))-A(IB+I) 01397 C(JD+J)=(A(IA+I)-A(IC+I))+B(IB+I) 01398 I=I+INC3 01399 J=J+INC4 01400 410 CONTINUE 01401 IBASE=IBASE+INC1 01402 JBASE=JBASE+INC2 01403 420 CONTINUE 01404 IA=IA+IINK 01405 IINK=2*IINK 01406 IB=IB+IINK 01407 IC=IC-IINK 01408 ID=ID-IINK 01409 JBASE=JBASE+JUMP 01410 JUMP=2*JUMP+JINK 01411 IF (IB.EQ.IC) GO TO 460 01412 DO 450 K=LA,KSTOP,LA 01413 KB=K+K 01414 KC=KB+KB 01415 KD=KC+KB 01416 C1=TRIGS(KB+1) 01417 S1=TRIGS(KB+2) 01418 C2=TRIGS(KC+1) 01419 S2=TRIGS(KC+2) 01420 C3=TRIGS(KD+1) 01421 S3=TRIGS(KD+2) 01422 IBASE=0 01423 DO 440 L=1,LA 01424 I=IBASE 01425 J=JBASE 01426 !DIR$ IVDEP 01427 DO 430 IJK=1,LOT 01428 C(JA+J)=(A(IA+I)+A(IC+I))+(A(IB+I)+A(ID+I)) 01429 D(JA+J)=(B(IA+I)-B(IC+I))+(B(IB+I)-B(ID+I)) 01430 C(JC+J)=C2*((A(IA+I)+A(IC+I))-(A(IB+I)+A(ID+I))) & 01431 -S2*((B(IA+I)-B(IC+I))-(B(IB+I)-B(ID+I))) 01432 D(JC+J)=S2*((A(IA+I)+A(IC+I))-(A(IB+I)+A(ID+I))) & 01433 +C2*((B(IA+I)-B(IC+I))-(B(IB+I)-B(ID+I))) 01434 C(JB+J)=C1*((A(IA+I)-A(IC+I))-(B(IB+I)+B(ID+I))) & 01435 -S1*((B(IA+I)+B(IC+I))+(A(IB+I)-A(ID+I))) 01436 D(JB+J)=S1*((A(IA+I)-A(IC+I))-(B(IB+I)+B(ID+I))) & 01437 +C1*((B(IA+I)+B(IC+I))+(A(IB+I)-A(ID+I))) 01438 C(JD+J)=C3*((A(IA+I)-A(IC+I))+(B(IB+I)+B(ID+I))) & 01439 -S3*((B(IA+I)+B(IC+I))-(A(IB+I)-A(ID+I))) 01440 D(JD+J)=S3*((A(IA+I)-A(IC+I))+(B(IB+I)+B(ID+I))) & 01441 +C3*((B(IA+I)+B(IC+I))-(A(IB+I)-A(ID+I))) 01442 I=I+INC3 01443 J=J+INC4 01444 430 CONTINUE 01445 IBASE=IBASE+INC1 01446 JBASE=JBASE+INC2 01447 440 CONTINUE 01448 IA=IA+IINK 01449 IB=IB+IINK 01450 IC=IC-IINK 01451 ID=ID-IINK 01452 JBASE=JBASE+JUMP 01453 450 CONTINUE 01454 IF (IB.GT.IC) GO TO 900 01455 460 CONTINUE 01456 IBASE=0 01457 SIN45=SQRT(0.5) 01458 DO 480 L=1,LA 01459 I=IBASE 01460 J=JBASE 01461 !DIR$ IVDEP 01462 DO 470 IJK=1,LOT 01463 C(JA+J)=A(IA+I)+A(IB+I) 01464 C(JB+J)=SIN45*((A(IA+I)-A(IB+I))-(B(IA+I)+B(IB+I))) 01465 C(JC+J)=B(IB+I)-B(IA+I) 01466 C(JD+J)=-SIN45*((A(IA+I)-A(IB+I))+(B(IA+I)+B(IB+I))) 01467 I=I+INC3 01468 J=J+INC4 01469 470 CONTINUE 01470 IBASE=IBASE+INC1 01471 JBASE=JBASE+INC2 01472 480 CONTINUE 01473 GO TO 900 01474 ! 01475 490 CONTINUE 01476 DO 494 L=1,LA 01477 I=IBASE 01478 J=JBASE 01479 !DIR$ IVDEP 01480 DO 492 IJK=1,LOT 01481 C(JA+J)=2.0*((A(IA+I)+A(IC+I))+A(IB+I)) 01482 C(JB+J)=2.0*((A(IA+I)-A(IC+I))-B(IB+I)) 01483 C(JC+J)=2.0*((A(IA+I)+A(IC+I))-A(IB+I)) 01484 C(JD+J)=2.0*((A(IA+I)-A(IC+I))+B(IB+I)) 01485 I=I+INC3 01486 J=J+INC4 01487 492 CONTINUE 01488 IBASE=IBASE+INC1 01489 JBASE=JBASE+INC2 01490 494 CONTINUE 01491 GO TO 900 01492 ! 01493 ! CODING FOR FACTOR 5 01494 ! ------------------- 01495 500 CONTINUE 01496 IA=1 01497 IB=IA+(2*M-LA)*INC1 01498 IC=IB+2*M*INC1 01499 ID=IC 01500 IE=IB 01501 JA=1 01502 JB=JA+JINK 01503 JC=JB+JINK 01504 JD=JC+JINK 01505 JE=JD+JINK 01506 ! 01507 IF (LA.EQ.M) GO TO 590 01508 ! 01509 DO 520 L=1,LA 01510 I=IBASE 01511 J=JBASE 01512 !DIR$ IVDEP 01513 DO 510 IJK=1,LOT 01514 C(JA+J)=A(IA+I)+(A(IB+I)+A(IC+I)) 01515 C(JB+J)=((A(IA+I)-0.25*(A(IB+I)+A(IC+I)))+QRT5*(A(IB+I)-A(IC+I & 01516 ))) -(SIN72*B(IB+I)+SIN36*B(IC+I)) 01517 C(JC+J)=((A(IA+I)-0.25*(A(IB+I)+A(IC+I)))-QRT5*(A(IB+I)-A(IC+I & 01518 ))) -(SIN36*B(IB+I)-SIN72*B(IC+I)) 01519 C(JD+J)=((A(IA+I)-0.25*(A(IB+I)+A(IC+I)))-QRT5*(A(IB+I)-A(IC+I & 01520 ))) +(SIN36*B(IB+I)-SIN72*B(IC+I)) 01521 C(JE+J)=((A(IA+I)-0.25*(A(IB+I)+A(IC+I)))+QRT5*(A(IB+I)-A(IC+I & 01522 ))) +(SIN72*B(IB+I)+SIN36*B(IC+I)) 01523 I=I+INC3 01524 J=J+INC4 01525 510 CONTINUE 01526 IBASE=IBASE+INC1 01527 JBASE=JBASE+INC2 01528 520 CONTINUE 01529 IA=IA+IINK 01530 IINK=2*IINK 01531 IB=IB+IINK 01532 IC=IC+IINK 01533 ID=ID-IINK 01534 IE=IE-IINK 01535 JBASE=JBASE+JUMP 01536 JUMP=2*JUMP+JINK 01537 IF (IB.EQ.ID) GO TO 560 01538 DO 550 K=LA,KSTOP,LA 01539 KB=K+K 01540 KC=KB+KB 01541 KD=KC+KB 01542 KE=KD+KB 01543 C1=TRIGS(KB+1) 01544 S1=TRIGS(KB+2) 01545 C2=TRIGS(KC+1) 01546 S2=TRIGS(KC+2) 01547 C3=TRIGS(KD+1) 01548 S3=TRIGS(KD+2) 01549 C4=TRIGS(KE+1) 01550 S4=TRIGS(KE+2) 01551 IBASE=0 01552 DO 540 L=1,LA 01553 I=IBASE 01554 J=JBASE 01555 !DIR$ IVDEP 01556 DO 530 IJK=1,LOT 01557 ! 01558 A10(IJK)=(A(IA+I)-0.25*((A(IB+I)+A(IE+I))+(A(IC+I)+A(ID+I))) & 01559 ) +QRT5*((A(IB+I)+A(IE+I))-(A(IC+I)+A(ID+I))) 01560 A20(IJK)=(A(IA+I)-0.25*((A(IB+I)+A(IE+I))+(A(IC+I)+A(ID+I))) & 01561 ) -QRT5*((A(IB+I)+A(IE+I))-(A(IC+I)+A(ID+I))) 01562 B10(IJK)=(B(IA+I)-0.25*((B(IB+I)-B(IE+I))+(B(IC+I)-B(ID+I))) & 01563 ) +QRT5*((B(IB+I)-B(IE+I))-(B(IC+I)-B(ID+I))) 01564 B20(IJK)=(B(IA+I)-0.25*((B(IB+I)-B(IE+I))+(B(IC+I)-B(ID+I))) & 01565 ) -QRT5*((B(IB+I)-B(IE+I))-(B(IC+I)-B(ID+I))) 01566 A11(IJK)=SIN72*(B(IB+I)+B(IE+I))+SIN36*(B(IC+I)+B(ID+I)) 01567 A21(IJK)=SIN36*(B(IB+I)+B(IE+I))-SIN72*(B(IC+I)+B(ID+I)) 01568 B11(IJK)=SIN72*(A(IB+I)-A(IE+I))+SIN36*(A(IC+I)-A(ID+I)) 01569 B21(IJK)=SIN36*(A(IB+I)-A(IE+I))-SIN72*(A(IC+I)-A(ID+I)) 01570 ! 01571 C(JA+J)=A(IA+I)+((A(IB+I)+A(IE+I))+(A(IC+I)+A(ID+I))) 01572 D(JA+J)=B(IA+I)+((B(IB+I)-B(IE+I))+(B(IC+I)-B(ID+I))) 01573 C(JB+J)=C1*(A10(IJK)-A11(IJK))-S1*(B10(IJK)+B11(IJK)) 01574 D(JB+J)=S1*(A10(IJK)-A11(IJK))+C1*(B10(IJK)+B11(IJK)) 01575 C(JE+J)=C4*(A10(IJK)+A11(IJK))-S4*(B10(IJK)-B11(IJK)) 01576 D(JE+J)=S4*(A10(IJK)+A11(IJK))+C4*(B10(IJK)-B11(IJK)) 01577 C(JC+J)=C2*(A20(IJK)-A21(IJK))-S2*(B20(IJK)+B21(IJK)) 01578 D(JC+J)=S2*(A20(IJK)-A21(IJK))+C2*(B20(IJK)+B21(IJK)) 01579 C(JD+J)=C3*(A20(IJK)+A21(IJK))-S3*(B20(IJK)-B21(IJK)) 01580 D(JD+J)=S3*(A20(IJK)+A21(IJK))+C3*(B20(IJK)-B21(IJK)) 01581 ! 01582 I=I+INC3 01583 J=J+INC4 01584 530 CONTINUE 01585 IBASE=IBASE+INC1 01586 JBASE=JBASE+INC2 01587 540 CONTINUE 01588 IA=IA+IINK 01589 IB=IB+IINK 01590 IC=IC+IINK 01591 ID=ID-IINK 01592 IE=IE-IINK 01593 JBASE=JBASE+JUMP 01594 550 CONTINUE 01595 IF (IB.GT.ID) GO TO 900 01596 560 CONTINUE 01597 IBASE=0 01598 DO 580 L=1,LA 01599 I=IBASE 01600 J=JBASE 01601 !DIR$ IVDEP 01602 DO 570 IJK=1,LOT 01603 C(JA+J)=(A(IA+I)+A(IB+I))+A(IC+I) 01604 C(JB+J)=(QRT5*(A(IA+I)-A(IB+I))+(0.25*(A(IA+I)+A(IB+I))-A(IC+I & 01605 ))) -(SIN36*B(IA+I)+SIN72*B(IB+I)) 01606 C(JE+J)=-(QRT5*(A(IA+I)-A(IB+I))+(0.25*(A(IA+I)+A(IB+I))-A(IC+ & 01607 I))) -(SIN36*B(IA+I)+SIN72*B(IB+I)) 01608 C(JC+J)=(QRT5*(A(IA+I)-A(IB+I))-(0.25*(A(IA+I)+A(IB+I))-A(IC+I & 01609 ))) -(SIN72*B(IA+I)-SIN36*B(IB+I)) 01610 C(JD+J)=-(QRT5*(A(IA+I)-A(IB+I))-(0.25*(A(IA+I)+A(IB+I))-A(IC+ & 01611 I))) -(SIN72*B(IA+I)-SIN36*B(IB+I)) 01612 I=I+INC3 01613 J=J+INC4 01614 570 CONTINUE 01615 IBASE=IBASE+INC1 01616 JBASE=JBASE+INC2 01617 580 CONTINUE 01618 GO TO 900 01619 ! 01620 590 CONTINUE 01621 QQRT5=2.0*QRT5 01622 SSIN36=2.0*SIN36 01623 SSIN72=2.0*SIN72 01624 DO 594 L=1,LA 01625 I=IBASE 01626 J=JBASE 01627 !DIR$ IVDEP 01628 DO 592 IJK=1,LOT 01629 C(JA+J)=2.0*(A(IA+I)+(A(IB+I)+A(IC+I))) 01630 C(JB+J)=(2.0*(A(IA+I)-0.25*(A(IB+I)+A(IC+I))) & 01631 +QQRT5*(A(IB+I)-A(IC+I)))-(SSIN72*B(IB+I)+SSIN36*B(IC+I)) 01632 C(JC+J)=(2.0*(A(IA+I)-0.25*(A(IB+I)+A(IC+I))) & 01633 -QQRT5*(A(IB+I)-A(IC+I)))-(SSIN36*B(IB+I)-SSIN72*B(IC+I)) 01634 C(JD+J)=(2.0*(A(IA+I)-0.25*(A(IB+I)+A(IC+I))) & 01635 -QQRT5*(A(IB+I)-A(IC+I)))+(SSIN36*B(IB+I)-SSIN72*B(IC+I)) 01636 C(JE+J)=(2.0*(A(IA+I)-0.25*(A(IB+I)+A(IC+I))) & 01637 +QQRT5*(A(IB+I)-A(IC+I)))+(SSIN72*B(IB+I)+SSIN36*B(IC+I)) 01638 I=I+INC3 01639 J=J+INC4 01640 592 CONTINUE 01641 IBASE=IBASE+INC1 01642 JBASE=JBASE+INC2 01643 594 CONTINUE 01644 GO TO 900 01645 ! 01646 ! CODING FOR FACTOR 6 01647 ! ------------------- 01648 600 CONTINUE 01649 IA=1 01650 IB=IA+(2*M-LA)*INC1 01651 IC=IB+2*M*INC1 01652 ID=IC+2*M*INC1 01653 IE=IC 01654 IF=IB 01655 JA=1 01656 JB=JA+JINK 01657 JC=JB+JINK 01658 JD=JC+JINK 01659 JE=JD+JINK 01660 JF=JE+JINK 01661 ! 01662 IF (LA.EQ.M) GO TO 690 01663 ! 01664 DO 620 L=1,LA 01665 I=IBASE 01666 J=JBASE 01667 !DIR$ IVDEP 01668 DO 610 IJK=1,LOT 01669 C(JA+J)=(A(IA+I)+A(ID+I))+(A(IB+I)+A(IC+I)) 01670 C(JD+J)=(A(IA+I)-A(ID+I))-(A(IB+I)-A(IC+I)) 01671 C(JB+J)=((A(IA+I)-A(ID+I))+0.5*(A(IB+I)-A(IC+I))) & 01672 -(SIN60*(B(IB+I)+B(IC+I))) 01673 C(JF+J)=((A(IA+I)-A(ID+I))+0.5*(A(IB+I)-A(IC+I))) & 01674 +(SIN60*(B(IB+I)+B(IC+I))) 01675 C(JC+J)=((A(IA+I)+A(ID+I))-0.5*(A(IB+I)+A(IC+I))) & 01676 -(SIN60*(B(IB+I)-B(IC+I))) 01677 C(JE+J)=((A(IA+I)+A(ID+I))-0.5*(A(IB+I)+A(IC+I))) & 01678 +(SIN60*(B(IB+I)-B(IC+I))) 01679 I=I+INC3 01680 J=J+INC4 01681 610 CONTINUE 01682 IBASE=IBASE+INC1 01683 JBASE=JBASE+INC2 01684 620 CONTINUE 01685 IA=IA+IINK 01686 IINK=2*IINK 01687 IB=IB+IINK 01688 IC=IC+IINK 01689 ID=ID-IINK 01690 IE=IE-IINK 01691 IF=IF-IINK 01692 JBASE=JBASE+JUMP 01693 JUMP=2*JUMP+JINK 01694 IF (IC.EQ.ID) GO TO 660 01695 DO 650 K=LA,KSTOP,LA 01696 KB=K+K 01697 KC=KB+KB 01698 KD=KC+KB 01699 KE=KD+KB 01700 KF=KE+KB 01701 C1=TRIGS(KB+1) 01702 S1=TRIGS(KB+2) 01703 C2=TRIGS(KC+1) 01704 S2=TRIGS(KC+2) 01705 C3=TRIGS(KD+1) 01706 S3=TRIGS(KD+2) 01707 C4=TRIGS(KE+1) 01708 S4=TRIGS(KE+2) 01709 C5=TRIGS(KF+1) 01710 S5=TRIGS(KF+2) 01711 IBASE=0 01712 DO 640 L=1,LA 01713 I=IBASE 01714 J=JBASE 01715 !DIR$ IVDEP 01716 DO 630 IJK=1,LOT 01717 ! 01718 A11(IJK)= (A(IE+I)+A(IB+I))+(A(IC+I)+A(IF+I)) 01719 A20(IJK)=(A(IA+I)+A(ID+I))-0.5*A11(IJK) 01720 A21(IJK)=SIN60*((A(IE+I)+A(IB+I))-(A(IC+I)+A(IF+I))) 01721 B11(IJK)= (B(IB+I)-B(IE+I))+(B(IC+I)-B(IF+I)) 01722 B20(IJK)=(B(IA+I)-B(ID+I))-0.5*B11(IJK) 01723 B21(IJK)=SIN60*((B(IB+I)-B(IE+I))-(B(IC+I)-B(IF+I))) 01724 ! 01725 C(JA+J)=(A(IA+I)+A(ID+I))+A11(IJK) 01726 D(JA+J)=(B(IA+I)-B(ID+I))+B11(IJK) 01727 C(JC+J)=C2*(A20(IJK)-B21(IJK))-S2*(B20(IJK)+A21(IJK)) 01728 D(JC+J)=S2*(A20(IJK)-B21(IJK))+C2*(B20(IJK)+A21(IJK)) 01729 C(JE+J)=C4*(A20(IJK)+B21(IJK))-S4*(B20(IJK)-A21(IJK)) 01730 D(JE+J)=S4*(A20(IJK)+B21(IJK))+C4*(B20(IJK)-A21(IJK)) 01731 ! 01732 A11(IJK)=(A(IE+I)-A(IB+I))+(A(IC+I)-A(IF+I)) 01733 B11(IJK)=(B(IE+I)+B(IB+I))-(B(IC+I)+B(IF+I)) 01734 A20(IJK)=(A(IA+I)-A(ID+I))-0.5*A11(IJK) 01735 A21(IJK)=SIN60*((A(IE+I)-A(IB+I))-(A(IC+I)-A(IF+I))) 01736 B20(IJK)=(B(IA+I)+B(ID+I))+0.5*B11(IJK) 01737 B21(IJK)=SIN60*((B(IE+I)+B(IB+I))+(B(IC+I)+B(IF+I))) 01738 ! 01739 C(JD+J)=C3*((A(IA+I)-A(ID+I))+A11(IJK))-S3*((B(IA+I)+B(ID+I & 01740 ))-B11(IJK)) 01741 D(JD+J)=S3*((A(IA+I)-A(ID+I))+A11(IJK))+C3*((B(IA+I)+B(ID+I & 01742 ))-B11(IJK)) 01743 C(JB+J)=C1*(A20(IJK)-B21(IJK))-S1*(B20(IJK)-A21(IJK)) 01744 D(JB+J)=S1*(A20(IJK)-B21(IJK))+C1*(B20(IJK)-A21(IJK)) 01745 C(JF+J)=C5*(A20(IJK)+B21(IJK))-S5*(B20(IJK)+A21(IJK)) 01746 D(JF+J)=S5*(A20(IJK)+B21(IJK))+C5*(B20(IJK)+A21(IJK)) 01747 ! 01748 I=I+INC3 01749 J=J+INC4 01750 630 CONTINUE 01751 IBASE=IBASE+INC1 01752 JBASE=JBASE+INC2 01753 640 CONTINUE 01754 IA=IA+IINK 01755 IB=IB+IINK 01756 IC=IC+IINK 01757 ID=ID-IINK 01758 IE=IE-IINK 01759 IF=IF-IINK 01760 JBASE=JBASE+JUMP 01761 650 CONTINUE 01762 IF (IC.GT.ID) GO TO 900 01763 660 CONTINUE 01764 IBASE=0 01765 DO 680 L=1,LA 01766 I=IBASE 01767 J=JBASE 01768 !DIR$ IVDEP 01769 DO 670 IJK=1,LOT 01770 C(JA+J)=A(IB+I)+(A(IA+I)+A(IC+I)) 01771 C(JD+J)=B(IB+I)-(B(IA+I)+B(IC+I)) 01772 C(JB+J)=(SIN60*(A(IA+I)-A(IC+I)))-(0.5*(B(IA+I)+B(IC+I))+B(IB+I)) 01773 C(JF+J)=-(SIN60*(A(IA+I)-A(IC+I)))-(0.5*(B(IA+I)+B(IC+I))+B(IB+I)) 01774 C(JC+J)=SIN60*(B(IC+I)-B(IA+I))+(0.5*(A(IA+I)+A(IC+I))-A(IB+I)) 01775 C(JE+J)=SIN60*(B(IC+I)-B(IA+I))-(0.5*(A(IA+I)+A(IC+I))-A(IB+I)) 01776 I=I+INC3 01777 J=J+INC4 01778 670 CONTINUE 01779 IBASE=IBASE+INC1 01780 JBASE=JBASE+INC2 01781 680 CONTINUE 01782 GO TO 900 01783 ! 01784 690 CONTINUE 01785 SSIN60=2.0*SIN60 01786 DO 694 L=1,LA 01787 I=IBASE 01788 J=JBASE 01789 !DIR$ IVDEP 01790 DO 692 IJK=1,LOT 01791 C(JA+J)=(2.0*(A(IA+I)+A(ID+I)))+(2.0*(A(IB+I)+A(IC+I))) 01792 C(JD+J)=(2.0*(A(IA+I)-A(ID+I)))-(2.0*(A(IB+I)-A(IC+I))) 01793 C(JB+J)=(2.0*(A(IA+I)-A(ID+I))+(A(IB+I)-A(IC+I))) & 01794 -(SSIN60*(B(IB+I)+B(IC+I))) 01795 C(JF+J)=(2.0*(A(IA+I)-A(ID+I))+(A(IB+I)-A(IC+I))) & 01796 +(SSIN60*(B(IB+I)+B(IC+I))) 01797 C(JC+J)=(2.0*(A(IA+I)+A(ID+I))-(A(IB+I)+A(IC+I))) & 01798 -(SSIN60*(B(IB+I)-B(IC+I))) 01799 C(JE+J)=(2.0*(A(IA+I)+A(ID+I))-(A(IB+I)+A(IC+I))) & 01800 +(SSIN60*(B(IB+I)-B(IC+I))) 01801 I=I+INC3 01802 J=J+INC4 01803 692 CONTINUE 01804 IBASE=IBASE+INC1 01805 JBASE=JBASE+INC2 01806 694 CONTINUE 01807 GO TO 900 01808 ! 01809 ! CODING FOR FACTOR 8 01810 ! ------------------- 01811 800 CONTINUE 01812 IBAD=3 01813 IF (LA.NE.M) GO TO 910 01814 IA=1 01815 IB=IA+LA*INC1 01816 IC=IB+2*LA*INC1 01817 ID=IC+2*LA*INC1 01818 IE=ID+2*LA*INC1 01819 JA=1 01820 JB=JA+JINK 01821 JC=JB+JINK 01822 JD=JC+JINK 01823 JE=JD+JINK 01824 JF=JE+JINK 01825 JG=JF+JINK 01826 JH=JG+JINK 01827 SSIN45=SQRT(2.0) 01828 ! 01829 DO 820 L=1,LA 01830 I=IBASE 01831 J=JBASE 01832 !DIR$ IVDEP 01833 DO 810 IJK=1,LOT 01834 C(JA+J)=2.0*(((A(IA+I)+A(IE+I))+A(IC+I))+(A(IB+I)+A(ID+I))) 01835 C(JE+J)=2.0*(((A(IA+I)+A(IE+I))+A(IC+I))-(A(IB+I)+A(ID+I))) 01836 C(JC+J)=2.0*(((A(IA+I)+A(IE+I))-A(IC+I))-(B(IB+I)-B(ID+I))) 01837 C(JG+J)=2.0*(((A(IA+I)+A(IE+I))-A(IC+I))+(B(IB+I)-B(ID+I))) 01838 C(JB+J)=2.0*((A(IA+I)-A(IE+I))-B(IC+I)) & 01839 +SSIN45*((A(IB+I)-A(ID+I))-(B(IB+I)+B(ID+I))) 01840 C(JF+J)=2.0*((A(IA+I)-A(IE+I))-B(IC+I)) & 01841 -SSIN45*((A(IB+I)-A(ID+I))-(B(IB+I)+B(ID+I))) 01842 C(JD+J)=2.0*((A(IA+I)-A(IE+I))+B(IC+I)) & 01843 -SSIN45*((A(IB+I)-A(ID+I))+(B(IB+I)+B(ID+I))) 01844 C(JH+J)=2.0*((A(IA+I)-A(IE+I))+B(IC+I)) & 01845 +SSIN45*((A(IB+I)-A(ID+I))+(B(IB+I)+B(ID+I))) 01846 I=I+INC3 01847 J=J+INC4 01848 810 CONTINUE 01849 IBASE=IBASE+INC1 01850 JBASE=JBASE+INC2 01851 820 CONTINUE 01852 ! 01853 ! RETURN 01854 ! ------ 01855 900 CONTINUE 01856 IBAD=0 01857 910 CONTINUE 01858 IERR=IBAD 01859 RETURN 01860 END 01861 01862
