CTITLERFFTMLT -- ROUTINE TO PERFORM MULTIPLE, REAL FFT'S 00010000 C*********************************************************************** C COPYRIGHT ATLANTIC RICHFIELD COMPANY 1991 * C*********************************************************************** CA AUTHOR JJ CHEN 00020000 CA DESIGNER JJ CHEN 00030000 CA LANGUAGE FORTRAN 00040000 CA SYSTEM IBM (SEE CRAY MATH LIBRARY) 00050001 CA WRITTEN 02/10/89 00060000 C REVISED MM/DD/YY III ... 00070000 C REVISED 09/24/89 JJC REINITIALIZE THE SCALE FACTOR. 00071000 C REVISED 03/20/90 LWC CHANGE NAME FROM FFT991 TO RFFTMLT FOR 00071100 C UNICOS CONVERSION. 00071200 CA 00071300 CA THIS SUBROUTINE PERFORMS MULTIPLE REAL FFT'S USING THE 00071400 CA IBM ESSL SOFTWARE. 00071500 CA 00071600 CA CALL RFFTMLT(A, WORK. TRIGS, IFAX, INC, JUMP, N, M, ISIGN) 00071700 CA 00071800 CA IN/OUT A = ARRAY OF LENGTH M*(N+2) R4 00071900 CA INPUT WORK = WORK ARRAY OF LENGTH M*(N+1) R4 00072000 CA INPUT TRIGS = ARRAY SET UP BY FFTFAX R4 00072100 CA INPUT IFAX = ARRAY SET UP BY FFTFAX I4 00073000 CA INPUT INC = INCREMENT BETWEEN ELEMENTS OF A I4 00074000 CA INPUT JUMP = INCREMENT (IN WORD PAIRS) BETWEEN THE I4 00075000 CA FIRST ELEMENTS OF SUCCESSIVE VECTORS. 00076000 CA INPUT N = LENGTH OF EACH TRANSFORM I4 00077000 CA INPUT M = NUMBER OF TRANSFORMS TO BE DONE I4 00078000 CA INPUT ISIGN = FORWARD/INVERSE SWITCH I4 00079000 CA +1 = WAVE-NUMBER TO SPACE-TIME DOMAIN 00080000 CA -1 = SPACE-TIME TO WAVE-NUMBER DOMAIN 00090000 CA 00100000 SUBROUTINE RFFTMLT(A, WORK, TRIGS, IFAX, INC, JUMP, N, M, 00110000 * ISIGN) 00120000 C 00130000 IMPLICIT INTEGER (A-Z) 00140000 C 00150000 EXTERNAL ENOTRM 00160000 C 00170000 C ================================= 00180000 C INTEGER ARRAYS -- INPUT ARGUMENTS 00190000 C ================================= 00200000 C 00210000 INTEGER IFAX (1) 00220000 C 00230000 C ================================= 00240000 C COMPLEX ARRAYS -- INPUT ARGUMENTS 00250000 C ================================= 00260000 C 00270000 REAL A (1) 00280000 REAL TRIGS (1) 00290000 REAL WORK (1) 00300000 REAL B (1) 00301000 REAL BUF (1) 00302000 C 00303000 C ======================= 00304000 C INTEGER ARRAYS -- LOCAL 00305000 C ======================= 00306000 C 00307000 INTEGER INCPRE(160) 00308000 INTEGER IPTRM1(160) 00309000 INTEGER IPTRP1(160) 00310000 INTEGER LENGTH(160) 00320000 INTEGER LOTPRE(160) 00330000 C 00340000 C ============== 00350000 C REAL VARIABLES 00360000 C ============== 00370000 C 00380000 REAL SCALE 00390000 C 00400000 C =============== 00410000 C DATA STATEMENTS 00420000 C =============== 00430000 C 00440000 DATA INCPRE / 160*0 / 00450000 DATA IPTRM1 / 160*0 / 00460000 DATA IPTRP1 / 160*0 / 00470000 DATA LENGTH / 160*0 / 00480000 DATA LOTPRE / 160*0 / 00490000 DATA NAUX1 / 50000 / 00500000 DATA NAUX2 / 50000 / 00510000 C 00520000 C INITIALIZE THE ESSL ERROR OPTION TABLE 00530000 C 00540000 CALL EINFO(0) 00550000 C 00560000 C MAKE ERROR CODE ESV2015 A RECOVERABLE ERROR 00570000 C SUPPRESS PRINTING ALL ERROR MESSAGES FOR IT 00580000 C 00590000 CALL ERRSET(2015,0,-1,0,ENOTRM) 00600000 C 00610000 C =========== 00620000 C PERFORM FFT 00630000 C =========== 00640000 C 00650000 IX = -ISIGN 00660000 C 00670000 IF (IX .EQ. 1) THEN 00680000 SCALE = 2.0 / SQRT(N * 2.0) 00690000 ELSE 00700000 SCALE = 1.0 / SQRT(N * 2.0) 00710000 ENDIF 00720000 C 00730000 IAUX3 = 1 00740000 NAUX3 = 0 00750000 C 00751000 LENI = 2 * JUMP 00752000 CALL GETMN2 (B, LENI, I , LENO) 00753000 IF (LENI .NE. LENO) CALL XDUMPX 00754000 IAUX4 = I + 1 00755000 C 00756000 DO 10 I = 1, 160 00757000 IPTR = I 00758000 IF (LENGTH(IPTR) .EQ. 0) GO TO 20 00759000 IF (LENGTH(IPTR) .EQ. N) THEN 00760000 IF (INCPRE(IPTR) .EQ. INC) THEN 00770000 IF (LOTPRE(IPTR) .EQ. M) THEN 00780000 GO TO 30 00790000 ENDIF 00800000 ENDIF 00810000 ENDIF 00811000 10 CONTINUE 00812000 CALL XDUMPX 00813000 C 00814000 C INITIALIZE SIN/COS ARRAYS 00815000 C 00816000 20 CONTINUE 00817000 LENGTH(IPTR) = N 00818000 INCPRE(IPTR) = INC 00819000 LOTPRE(IPTR) = M 00820000 LENI = 4 * (NAUX1+NAUX2) 00830000 CALL GETMN2 (BUF, LENI, I, LENO) 00840000 IF (LENI .NE. LENO) CALL XDUMPX 00850000 IPTRP1(IPTR) = I + 1 00860000 IAUX1 = IPTRP1(IPTR) 00870000 IAUX2 = IAUX1 + 2*NAUX1 00880000 SCALE = 2.0 / SQRT(N * 2.0) 00890000 25 CALL SRCFT (1, A, JUMP, A, JUMP, N, 1, 1, SCALE, BUF(IAUX1), 00900000 * NAUX1, BUF(IAUX2), NAUX2, BUF(IAUX3), NAUX3, *40) 00910000 C 00920000 IPTRM1(IPTR) = IAUX2 + 2*NAUX2 00930000 IAUX1 = IPTRM1(IPTR) 00940000 IAUX2 = IAUX1 + 2*NAUX1 00950000 SCALE = 1.0 / SQRT(N * 2.0) 00960000 CALL SCRFT (1, A, JUMP, A, JUMP, N, 1, -1, SCALE, BUF(IAUX1), 00970000 * NAUX1, BUF(IAUX2), NAUX2, BUF(IAUX3), NAUX3) 00980000 C 00990000 C DO ACTUAL FFT 01000000 C 01010000 30 IF (IX .EQ. 1) THEN 01020000 IAUX1 = IPTRP1(IPTR) 01030000 IAUX2 = IAUX1 + 2*NAUX1 01040000 SCALE = 2.0 / SQRT(N * 2.0) 01050000 DO 33 JJ = 1, M 01060000 KK = (JJ-1) * JUMP + 1 01070000 CALL SCOPY(JUMP, A(KK), 1, B(IAUX4), 1) 01080000 CALL SRCFT (0, B(IAUX4), JUMP, B(IAUX4), JUMP, N, 1, IX, 01090000 * SCALE, BUF(IAUX1), NAUX1, BUF(IAUX2), NAUX2, 01100000 * BUF(IAUX3),NAUX3 ) 01110000 CALL SCOPY(JUMP, B(IAUX4), 1, A(KK), 1) 01120000 33 CONTINUE 01130000 ELSE 01140000 IAUX1 = IPTRM1(IPTR) 01150000 IAUX2 = IAUX1 + 2*NAUX1 01160000 SCALE = 1.0 / SQRT(N * 2.0) 01170000 DO 35 JJ = 1, M 01180000 KK = (JJ-1) * JUMP + 1 01190000 CALL SCOPY(JUMP, A(KK), 1, B(IAUX4), 1) 01200000 CALL SCRFT (0, B(IAUX4), JUMP, B(IAUX4), JUMP, N, 1, IX, 01201000 * SCALE, BUF(IAUX1), NAUX1, BUF(IAUX2), NAUX2, 01202000 * BUF(IAUX3),NAUX3 ) 01203000 CALL SCOPY(JUMP, B(IAUX4), 1, A(KK), 1) 01204000 35 CONTINUE 01205000 ENDIF 01206000 C 01207000 GO TO 8000 01208000 C 01209000 C CHECK THE RESULTING INPUT ARGUMENT VALUE 01210000 C IN NAUX AND TAKE THE DESIRED ACTION 01220000 C (RETURN THE NECESSARY SIZE FOR NAUX1, NAUX2) 01221000 C 01222000 40 CONTINUE 01223000 LENI = 4 * (NAUX1+NAUX2) 01224000 CALL GETMN2 (BUF, LENI, I, LENO) 01225000 IF (LENI .NE. LENO) CALL XDUMPX 01226000 IPTRP1(IPTR) = I + 1 01227000 IAUX1 = IPTRP1(IPTR) 01228000 IAUX2 = IAUX1 + 2 * NAUX1 01229000 GO TO 25 01230000 C 01240000 C ==== 01250000 C EXIT 01260000 C ==== 01270000 C 01280000 8000 RETURN 01281000 C 01282000 END 01283000