Commit c240a96f authored by Raymond Toy's avatar Raymond Toy
Browse files

Import most fixes from FFTPACK51

Most if the fixes are due to computing the lengths (and indices) for
LENSAV.  Apparently the lengths should be log2(N) instead of log(N).

But the fix for c1fm1f.f fixes a genuine bug in the complex fft which
caused wrong results for non-power of 2 FFT sizes.

Also fixed a couple of typos and incorrect return ier codes.
parent f8319349
...@@ -44,7 +44,7 @@ c write(*,*) wa(iw),wa(iw+1) ...@@ -44,7 +44,7 @@ c write(*,*) wa(iw),wa(iw+1)
65 CALL C1F5KF (IDO,L1,NA,CH,2,C,INC2,WA(IW)) 65 CALL C1F5KF (IDO,L1,NA,CH,2,C,INC2,WA(IW))
GO TO 120 GO TO 120
56 CALL C1FGKF (IDO,IP,L1,LID,NA,C,C,INC2,CH,CH, 56 CALL C1FGKF (IDO,IP,L1,LID,NA,C,C,INC2,CH,CH,
1 1,WA(IW)) 1 2,WA(IW))
GO TO 120 GO TO 120
66 CALL C1FGKF (IDO,IP,L1,LID,NA,CH,CH,2,C,C, 66 CALL C1FGKF (IDO,IP,L1,LID,NA,CH,CH,2,C,C,
1 INC2,WA(IW)) 1 INC2,WA(IW))
......
...@@ -18,13 +18,13 @@ C ...@@ -18,13 +18,13 @@ C
C C
IF (LENC .LT. INC*(N-1) + 1) THEN IF (LENC .LT. INC*(N-1) + 1) THEN
IER = 1 IER = 1
CALL XERFFT ('CFFT1B ', 6) CALL XERFFT ('CFFT1B ', 4)
ELSEIF (LENSAV .LT. 2*N + INT(LOG(REAL(N))) + 4) THEN ELSEIF (LENSAV .LT. 2*N + INT(LOG(REAL(N))/LOG(2.)) + 4) THEN
IER = 2 IER = 2
CALL XERFFT ('CFFT1B ', 8) CALL XERFFT ('CFFT1B ', 6)
ELSEIF (LENWRK .LT. 2*N) THEN ELSEIF (LENWRK .LT. 2*N) THEN
IER = 3 IER = 3
CALL XERFFT ('CFFT1B ', 10) CALL XERFFT ('CFFT1B ', 8)
ENDIF ENDIF
C C
IF (N .EQ. 1) RETURN IF (N .EQ. 1) RETURN
......
...@@ -18,13 +18,13 @@ C ...@@ -18,13 +18,13 @@ C
C C
IF (LENC .LT. INC*(N-1) + 1) THEN IF (LENC .LT. INC*(N-1) + 1) THEN
IER = 1 IER = 1
CALL XERFFT ('CFFT1F ', 6) CALL XERFFT ('CFFT1F ', 4)
ELSEIF (LENSAV .LT. 2*N + INT(LOG(REAL(N))) + 4) THEN ELSEIF (LENSAV .LT. 2*N + INT(LOG(REAL(N))/LOG(2.)) + 4) THEN
IER = 2 IER = 2
CALL XERFFT ('CFFT1F ', 8) CALL XERFFT ('CFFT1F ', 6)
ELSEIF (LENWRK .LT. 2*N) THEN ELSEIF (LENWRK .LT. 2*N) THEN
IER = 3 IER = 3
CALL XERFFT ('CFFTMF ', 10) CALL XERFFT ('CFFT1F ', 8)
ENDIF ENDIF
C C
IF (N .EQ. 1) RETURN IF (N .EQ. 1) RETURN
......
...@@ -14,7 +14,7 @@ CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC ...@@ -14,7 +14,7 @@ CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C C
IER = 0 IER = 0
C C
IF (LENSAV .LT. 2*N + INT(LOG(REAL(N))) + 4) THEN IF (LENSAV .LT. 2*N + INT(LOG(REAL(N))/LOG(2.)) + 4) THEN
IER = 2 IER = 2
CALL XERFFT ('CFFTMI ', 3) CALL XERFFT ('CFFTMI ', 3)
ENDIF ENDIF
......
...@@ -22,8 +22,8 @@ C ...@@ -22,8 +22,8 @@ C
IER = 5 IER = 5
CALL XERFFT ('CFFT2B', -2) CALL XERFFT ('CFFT2B', -2)
GO TO 100 GO TO 100
ELSEIF (LENSAV .LT. 2*L + INT(LOG(REAL(L))) + ELSEIF (LENSAV .LT. 2*L + INT(LOG(REAL(L))/LOG(2.)) +
1 2*M + INT(LOG(REAL(M))) +8) THEN 1 2*M + INT(LOG(REAL(M))/LOG(2.)) +8) THEN
IER = 2 IER = 2
CALL XERFFT ('CFFT2B', 6) CALL XERFFT ('CFFT2B', 6)
GO TO 100 GO TO 100
...@@ -34,9 +34,9 @@ C ...@@ -34,9 +34,9 @@ C
ENDIF ENDIF
C C
C Transform X lines of C array C Transform X lines of C array
IW = 2*L+INT(LOG(REAL(L))*LOG(2.)) + 3 IW = 2*L+INT(LOG(REAL(L))/LOG(2.)) + 3
CALL CFFTMB(L, 1, M, LDIM, C, (L-1) + LDIM*(M-1) +1, CALL CFFTMB(L, 1, M, LDIM, C, (L-1) + LDIM*(M-1) +1,
1 WSAVE(IW), 2*M + INT(LOG(REAL(M))) + 4, 1 WSAVE(IW), 2*M + INT(LOG(REAL(M))/LOG(2.)) + 4,
2 WORK, 2*L*M, IER1) 2 WORK, 2*L*M, IER1)
IF (IER1 .NE. 0) THEN IF (IER1 .NE. 0) THEN
IER = 20 IER = 20
...@@ -47,7 +47,7 @@ C ...@@ -47,7 +47,7 @@ C
C Transform Y lines of C array C Transform Y lines of C array
IW = 1 IW = 1
CALL CFFTMB (M, LDIM, L, 1, C, (M-1)*LDIM + L, CALL CFFTMB (M, LDIM, L, 1, C, (M-1)*LDIM + L,
1 WSAVE(IW), 2*L + INT(LOG(REAL(L))) + 4, 1 WSAVE(IW), 2*L + INT(LOG(REAL(L))/LOG(2.)) + 4,
2 WORK, 2*M*L, IER1) 2 WORK, 2*M*L, IER1)
IF (IER1 .NE. 0) THEN IF (IER1 .NE. 0) THEN
IER = 20 IER = 20
......
...@@ -22,8 +22,8 @@ C ...@@ -22,8 +22,8 @@ C
IER = 5 IER = 5
CALL XERFFT ('CFFT2F', -2) CALL XERFFT ('CFFT2F', -2)
GO TO 100 GO TO 100
ELSEIF (LENSAV .LT. 2*L + INT(LOG(REAL(L))) + ELSEIF (LENSAV .LT. 2*L + INT(LOG(REAL(L))/LOG(2.)) +
1 2*M + INT(LOG(REAL(M))) +8) THEN 1 2*M + INT(LOG(REAL(M))/LOG(2.)) +8) THEN
IER = 2 IER = 2
CALL XERFFT ('CFFT2F', 6) CALL XERFFT ('CFFT2F', 6)
GO TO 100 GO TO 100
...@@ -34,9 +34,9 @@ C ...@@ -34,9 +34,9 @@ C
ENDIF ENDIF
C C
C Transform X lines of C array C Transform X lines of C array
IW = 2*L+INT(LOG(REAL(L))*LOG(2.)) + 3 IW = 2*L+INT(LOG(REAL(L))/LOG(2.)) + 3
CALL CFFTMF(L, 1, M, LDIM, C, (L-1) + LDIM*(M-1) +1, CALL CFFTMF(L, 1, M, LDIM, C, (L-1) + LDIM*(M-1) +1,
1 WSAVE(IW), 2*M + INT(LOG(REAL(M))) + 4, 1 WSAVE(IW), 2*M + INT(LOG(REAL(M))/LOG(2.)) + 4,
2 WORK, 2*L*M, IER1) 2 WORK, 2*L*M, IER1)
IF (IER1 .NE. 0) THEN IF (IER1 .NE. 0) THEN
IER = 20 IER = 20
...@@ -47,7 +47,7 @@ C ...@@ -47,7 +47,7 @@ C
C Transform Y lines of C array C Transform Y lines of C array
IW = 1 IW = 1
CALL CFFTMF (M, LDIM, L, 1, C, (M-1)*LDIM + L, CALL CFFTMF (M, LDIM, L, 1, C, (M-1)*LDIM + L,
1 WSAVE(IW), 2*L + INT(LOG(REAL(L))) + 4, 1 WSAVE(IW), 2*L + INT(LOG(REAL(L))/LOG(2.)) + 4,
2 WORK, 2*M*L, IER1) 2 WORK, 2*M*L, IER1)
IF (IER1 .NE. 0) THEN IF (IER1 .NE. 0) THEN
IER = 20 IER = 20
......
...@@ -16,21 +16,22 @@ C Initialize error return ...@@ -16,21 +16,22 @@ C Initialize error return
C C
IER = 0 IER = 0
C C
IF (LENSAV .LT. 2*L + INT(LOG(REAL(L))) + IF (LENSAV .LT. 2*L + INT(LOG(REAL(L))/LOG(2.)) +
1 2*M + INT(LOG(REAL(M))) +8) THEN 1 2*M + INT(LOG(REAL(M))/LOG(2.)) +8) THEN
IER = 2 IER = 2
CALL XERFFT ('CFFT2I', 4) CALL XERFFT ('CFFT2I', 4)
GO TO 100 GO TO 100
ENDIF ENDIF
C C
CALL CFFTMI (L, WSAVE(1), 2*L + INT(LOG(REAL(L))) + 4, IER1) CALL CFFTMI (L, WSAVE(1), 2*L + INT(LOG(REAL(L))/LOG(2.)) + 4,
1 IER1)
IF (IER1 .NE. 0) THEN IF (IER1 .NE. 0) THEN
IER = 20 IER = 20
CALL XERFFT ('CFFT2I',-5) CALL XERFFT ('CFFT2I',-5)
GO TO 100 GO TO 100
ENDIF ENDIF
CALL CFFTMI (M, WSAVE(2*L+INT(LOG(REAL(L))*LOG(2.)) + 3), CALL CFFTMI (M, WSAVE(2*L+INT(LOG(REAL(L))/LOG(2.)) + 3),
1 2*M + INT(LOG(REAL(M))) + 4, IER1) 1 2*M + INT(LOG(REAL(M))/LOG(2.)) + 4, IER1)
IF (IER1 .NE. 0) THEN IF (IER1 .NE. 0) THEN
IER = 20 IER = 20
CALL XERFFT ('CFFT2I',-5) CALL XERFFT ('CFFT2I',-5)
......
...@@ -20,7 +20,7 @@ C ...@@ -20,7 +20,7 @@ C
IF (LENC .LT. (LOT-1)*JUMP + INC*(N-1) + 1) THEN IF (LENC .LT. (LOT-1)*JUMP + INC*(N-1) + 1) THEN
IER = 1 IER = 1
CALL XERFFT ('CFFTMB ', 6) CALL XERFFT ('CFFTMB ', 6)
ELSEIF (LENSAV .LT. 2*N + INT(LOG(REAL(N))) + 4) THEN ELSEIF (LENSAV .LT. 2*N + INT(LOG(REAL(N))/LOG(2.)) + 4) THEN
IER = 2 IER = 2
CALL XERFFT ('CFFTMB ', 8) CALL XERFFT ('CFFTMB ', 8)
ELSEIF (LENWRK .LT. 2*LOT*N) THEN ELSEIF (LENWRK .LT. 2*LOT*N) THEN
......
...@@ -20,7 +20,7 @@ C ...@@ -20,7 +20,7 @@ C
IF (LENC .LT. (LOT-1)*JUMP + INC*(N-1) + 1) THEN IF (LENC .LT. (LOT-1)*JUMP + INC*(N-1) + 1) THEN
IER = 1 IER = 1
CALL XERFFT ('CFFTMF ', 6) CALL XERFFT ('CFFTMF ', 6)
ELSEIF (LENSAV .LT. 2*N + INT(LOG(REAL(N))) + 4) THEN ELSEIF (LENSAV .LT. 2*N + INT(LOG(REAL(N))/LOG(2.)) + 4) THEN
IER = 2 IER = 2
CALL XERFFT ('CFFTMF ', 8) CALL XERFFT ('CFFTMF ', 8)
ELSEIF (LENWRK .LT. 2*LOT*N) THEN ELSEIF (LENWRK .LT. 2*LOT*N) THEN
......
...@@ -14,7 +14,7 @@ CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC ...@@ -14,7 +14,7 @@ CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C C
IER = 0 IER = 0
C C
IF (LENSAV .LT. 2*N + INT(LOG(REAL(N))) + 4) THEN IF (LENSAV .LT. 2*N + INT(LOG(REAL(N))/LOG(2.)) + 4) THEN
IER = 2 IER = 2
CALL XERFFT ('CFFTMI ', 3) CALL XERFFT ('CFFTMI ', 3)
ENDIF ENDIF
......
...@@ -19,7 +19,7 @@ C ...@@ -19,7 +19,7 @@ C
IER = 1 IER = 1
CALL XERFFT ('COSQ1B', 6) CALL XERFFT ('COSQ1B', 6)
GO TO 300 GO TO 300
ELSEIF (LENSAV .LT. 2*N + INT(LOG(REAL(N))) +4) THEN ELSEIF (LENSAV .LT. 2*N + INT(LOG(REAL(N))/LOG(2.)) +4) THEN
IER = 2 IER = 2
CALL XERFFT ('COSQ1B', 8) CALL XERFFT ('COSQ1B', 8)
GO TO 300 GO TO 300
......
...@@ -19,7 +19,7 @@ C ...@@ -19,7 +19,7 @@ C
IER = 1 IER = 1
CALL XERFFT ('COSQ1F', 6) CALL XERFFT ('COSQ1F', 6)
GO TO 300 GO TO 300
ELSEIF (LENSAV .LT. 2*N + INT(LOG(REAL(N))) +4) THEN ELSEIF (LENSAV .LT. 2*N + INT(LOG(REAL(N))/LOG(2.)) +4) THEN
IER = 2 IER = 2
CALL XERFFT ('COSQ1F', 8) CALL XERFFT ('COSQ1F', 8)
GO TO 300 GO TO 300
......
...@@ -13,7 +13,7 @@ CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC ...@@ -13,7 +13,7 @@ CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
REAL WSAVE(LENSAV) REAL WSAVE(LENSAV)
C C
IER = 0 IER = 0
IF (LENSAV .LT. 2*N + INT(LOG(REAL(N))) +4) THEN IF (LENSAV .LT. 2*N + INT(LOG(REAL(N))/LOG(2.)) +4) THEN
IER = 2 IER = 2
CALL XERFFT ('COSQ1I', 3) CALL XERFFT ('COSQ1I', 3)
GO TO 300 GO TO 300
...@@ -26,7 +26,7 @@ C ...@@ -26,7 +26,7 @@ C
FK = FK+1. FK = FK+1.
WSAVE(K) = COS(FK*DT) WSAVE(K) = COS(FK*DT)
101 CONTINUE 101 CONTINUE
LNSV = N + INT(LOG(REAL(N))) +4 LNSV = N + INT(LOG(REAL(N))/LOG(2.)) +4
CALL RFFT1I (N, WSAVE(N+1), LNSV, IER1) CALL RFFT1I (N, WSAVE(N+1), LNSV, IER1)
IF (IER1 .NE. 0) THEN IF (IER1 .NE. 0) THEN
IER = 20 IER = 20
......
...@@ -23,7 +23,7 @@ CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC ...@@ -23,7 +23,7 @@ CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
IF (MODN .NE. 0) GO TO 302 IF (MODN .NE. 0) GO TO 302
X(1,N) = .5*X(1,N) X(1,N) = .5*X(1,N)
302 LENX = INC*(N-1) + 1 302 LENX = INC*(N-1) + 1
LNSV = N + INT(LOG(REAL(N))) + 4 LNSV = N + INT(LOG(REAL(N))/LOG(2.)) + 4
LNWK = N LNWK = N
C C
CALL RFFT1B(N,INC,X,LENX,WSAVE(N+1),LNSV,WORK,LNWK,IER1) CALL RFFT1B(N,INC,X,LENX,WSAVE(N+1),LNSV,WORK,LNWK,IER1)
......
...@@ -29,7 +29,7 @@ CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC ...@@ -29,7 +29,7 @@ CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
IF (MODN .NE. 0) GO TO 303 IF (MODN .NE. 0) GO TO 303
X(1,NS2+1) = WSAVE(NS2)*WORK(NS2+1) X(1,NS2+1) = WSAVE(NS2)*WORK(NS2+1)
303 LENX = INC*(N-1) + 1 303 LENX = INC*(N-1) + 1
LNSV = N + INT(LOG(REAL(N))) + 4 LNSV = N + INT(LOG(REAL(N))/LOG(2.)) + 4
LNWK = N LNWK = N
C C
CALL RFFT1F(N,INC,X,LENX,WSAVE(N+1),LNSV,WORK,LNWK,IER1) CALL RFFT1F(N,INC,X,LENX,WSAVE(N+1),LNSV,WORK,LNWK,IER1)
......
...@@ -20,7 +20,7 @@ C ...@@ -20,7 +20,7 @@ C
IER = 1 IER = 1
CALL XERFFT ('COSQMB', 6) CALL XERFFT ('COSQMB', 6)
GO TO 300 GO TO 300
ELSEIF (LENSAV .LT. 2*N + INT(LOG(REAL(N))) +4) THEN ELSEIF (LENSAV .LT. 2*N + INT(LOG(REAL(N))/LOG(2.)) +4) THEN
IER = 2 IER = 2
CALL XERFFT ('COSQMB', 8) CALL XERFFT ('COSQMB', 8)
GO TO 300 GO TO 300
......
...@@ -20,7 +20,7 @@ C ...@@ -20,7 +20,7 @@ C
IER = 1 IER = 1
CALL XERFFT ('COSQMF', 6) CALL XERFFT ('COSQMF', 6)
GO TO 300 GO TO 300
ELSEIF (LENSAV .LT. 2*N + INT(LOG(REAL(N))) +4) THEN ELSEIF (LENSAV .LT. 2*N + INT(LOG(REAL(N))/LOG(2.)) +4) THEN
IER = 2 IER = 2
CALL XERFFT ('COSQMF', 8) CALL XERFFT ('COSQMF', 8)
GO TO 300 GO TO 300
......
...@@ -14,7 +14,7 @@ CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC ...@@ -14,7 +14,7 @@ CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C C
IER = 0 IER = 0
C C
IF (LENSAV .LT. 2*N + INT(LOG(REAL(N))) +4) THEN IF (LENSAV .LT. 2*N + INT(LOG(REAL(N))/LOG(2.)) +4) THEN
IER = 2 IER = 2
CALL XERFFT ('COSQMI', 3) CALL XERFFT ('COSQMI', 3)
GO TO 300 GO TO 300
...@@ -27,7 +27,7 @@ C ...@@ -27,7 +27,7 @@ C
FK = FK+1. FK = FK+1.
WSAVE(K) = COS(FK*DT) WSAVE(K) = COS(FK*DT)
101 CONTINUE 101 CONTINUE
LNSV = N + INT(LOG(REAL(N))) +4 LNSV = N + INT(LOG(REAL(N))/LOG(2.)) +4
CALL RFFTMI (N, WSAVE(N+1), LNSV, IER1) CALL RFFTMI (N, WSAVE(N+1), LNSV, IER1)
IF (IER1 .NE. 0) THEN IF (IER1 .NE. 0) THEN
IER = 20 IER = 20
......
...@@ -18,7 +18,7 @@ C ...@@ -18,7 +18,7 @@ C
IER = 1 IER = 1
CALL XERFFT ('COST1B', 6) CALL XERFFT ('COST1B', 6)
GO TO 100 GO TO 100
ELSEIF (LENSAV .LT. 2*N + INT(LOG(REAL(N))) +4) THEN ELSEIF (LENSAV .LT. 2*N + INT(LOG(REAL(N))/LOG(2.)) +4) THEN
IER = 2 IER = 2
CALL XERFFT ('COST1B', 8) CALL XERFFT ('COST1B', 8)
GO TO 100 GO TO 100
...@@ -27,6 +27,8 @@ C ...@@ -27,6 +27,8 @@ C
CALL XERFFT ('COST1B', 10) CALL XERFFT ('COST1B', 10)
GO TO 100 GO TO 100
ENDIF ENDIF
C
IF (N .EQ. 1) RETURN
C C
CALL COSTB1 (N,INC,X,WSAVE,WORK,IER1) CALL COSTB1 (N,INC,X,WSAVE,WORK,IER1)
IF (IER1 .NE. 0) THEN IF (IER1 .NE. 0) THEN
......
...@@ -18,7 +18,7 @@ C ...@@ -18,7 +18,7 @@ C
IER = 1 IER = 1
CALL XERFFT ('COST1F', 6) CALL XERFFT ('COST1F', 6)
GO TO 100 GO TO 100
ELSEIF (LENSAV .LT. 2*N + INT(LOG(REAL(N))) +4) THEN ELSEIF (LENSAV .LT. 2*N + INT(LOG(REAL(N))/LOG(2.)) +4) THEN
IER = 2 IER = 2
CALL XERFFT ('COST1F', 8) CALL XERFFT ('COST1F', 8)
GO TO 100 GO TO 100
...@@ -27,6 +27,8 @@ C ...@@ -27,6 +27,8 @@ C
CALL XERFFT ('COST1F', 10) CALL XERFFT ('COST1F', 10)
GO TO 100 GO TO 100
ENDIF ENDIF
C
IF (N .EQ. 1) RETURN
C C
CALL COSTF1(N,INC,X,WSAVE,WORK,IER1) CALL COSTF1(N,INC,X,WSAVE,WORK,IER1)
IF (IER1 .NE. 0) THEN IF (IER1 .NE. 0) THEN
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment