$NetBSD: patch-ag,v 1.5 2008/06/21 20:00:22 joerg Exp $ Reorder functions in use order to fix compilation with f2c. --- f_source/sciport/fftd.f.orig 2008-06-21 13:43:11.000000000 +0000 +++ f_source/sciport/fftd.f @@ -1,3 +1,192 @@ +C----------------------------------------------- ************ +C CABLE2 +C ************ + SUBROUTINE SPDCABLE2(NN,WORK) +C + DOUBLE PRECISION WORK(2,NN),P2,TWOPI + DATA TWOPI /6.2831853071795864769252867665590057683943388D0/ +C + N = 2 * NN + P2 = TWOPI/N + DO 10 I=1, NN + WORK(1,I) = DCOS(P2 * (I-1)) + WORK(2,I) = DSIN(P2 * (I-1)) + 10 CONTINUE + RETURN + END + +C----------------------------------------------- ************ +C RCONV2 +C ************ + SUBROUTINE SPDRCONV2(N,CY,C,CH) +C + DOUBLE COMPLEX CY(1) + DOUBLE PRECISION CH(N/2,2),P(2,1),C(2,N/2) + DOUBLE PRECISION X,Y,Z,Z1 +C + N2 = N/2 + P(1,1) = (C(1,1) + C(2,1)) * 2 + P(2,1) = (C(1,1) - C(2,1)) * 2 + CY(1) = DCMPLX(P(1,1),0.0D0) + CY(N2+1) = DCMPLX(P(2,1),0.0D0) + K = N2 + DO 10 I=2, N2 + X = C(1,I)+C(1,K) + Y = C(2,I)+C(2,K) + Z = C(1,I)-C(1,K) + Z1= C(2,I)-C(2,K) + P(1,1) = X + CH(I,1) * Y - CH(I,2) * Z + P(2,1) = Z1 - CH(I,2) * Y - CH(I,1) * Z + CY(I) = DCMPLX(P(1,1),P(2,1)) + K = K - 1 + 10 CONTINUE + RETURN + END + +C----------------------------------------------- ************ +C RCONV1 +C ************ + SUBROUTINE SPDRCONV1(N,CY,C,CH) +C + DOUBLE COMPLEX CY(1) + DOUBLE PRECISION CH(N/2,2),P(2,1),C(2,N/2) + DOUBLE PRECISION X,Y,Z,Z1 +C + N2 = N/2 + P(1,1) = (C(1,1) + C(2,1)) * 2 + P(2,1) = (C(1,1) - C(2,1)) * 2 + CY(1) = DCMPLX(P(1,1),0.0D0) + CY(N2+1) = DCMPLX(P(2,1),0.0D0) + K = N2 + DO 10 I=2, N2 + X = C(1,I)+C(1,K) + Y = C(2,I)+C(2,K) + Z = C(1,I)-C(1,K) + Z1= C(2,I)-C(2,K) + P(1,1) = X + CH(I,1) * Y + CH(I,2) * Z + P(2,1) = Z1 + CH(I,2) * Y - CH(I,1) * Z + CY(I) = DCMPLX(P(1,1),P(2,1)) + K = K - 1 + 10 CONTINUE + RETURN + END + +C----------------------------------------------- ************ +C RTOCK3 +C ************ + SUBROUTINE SPDRTOCK3(LS,NS,C,CH,CH2) +C + DOUBLE COMPLEX WYK,C(NS,LS,2),CH(NS,2,LS) + DOUBLE PRECISION CH2(2,NS,LS,2) +C + IF (LS .GT. NS) GOTO 30 + DO 600 I=1, LS + DO 600 J=1, NS + WYK =CONJG(DCMPLX(CH2(1,1,I,1),CH2(1,1,I,2))) * CH(J,2,I) + C(J,I,1) = CH(J,1,I) + WYK + C(J,I,2) = CH(J,1,I) - WYK + 600 CONTINUE + RETURN + 30 CONTINUE + DO 800 J=1, NS + DO 800 I=1, LS + WYK =CONJG(DCMPLX(CH2(1,1,I,1),CH2(1,1,I,2))) * CH(J,2,I) + C(J,I,1) = CH(J,1,I) + WYK + C(J,I,2) = CH(J,1,I) - WYK + 800 CONTINUE + RETURN + END + +C----------------------------------------------- ************ +C RABLE1 +C ************ + SUBROUTINE SPDRABLE1(NN,WORK) +C + DOUBLE PRECISION WORK(NN,2),P2,TWOPI + DATA TWOPI /6.2831853071795864769252867665590057683943388D0/ +C + N = 2 * NN + P2 = TWOPI/N + DO 10 I=1, NN + WORK(I,1) = DCOS(P2 * (I-1)) + WORK(I,2) = DSIN(P2 * (I-1)) + 10 CONTINUE + RETURN + END + +C----------------------------------------------- ************ +C RTOCK2 +C ************ + SUBROUTINE SPDRTOCK2(LS,NS,C,CH,CH2) +C + DOUBLE COMPLEX WYK,C(NS,LS,2),CH(NS,2,LS) + DOUBLE PRECISION CH2(2,NS,LS,2) +C + IF (LS .GT. NS) GOTO 20 + DO 200 I=1, LS + DO 200 J=1, NS + WYK = DCMPLX(CH2(1,1,I,1),CH2(1,1,I,2)) * CH(J,2,I) + C(J,I,1) = CH(J,1,I) + WYK + C(J,I,2) = CH(J,1,I) - WYK + 200 CONTINUE + RETURN + 20 CONTINUE + DO 400 J=1, NS + DO 400 I=1, LS + WYK = DCMPLX(CH2(1,1,I,1),CH2(1,1,I,2)) * CH(J,2,I) + C(J,I,1) = CH(J,1,I) + WYK + C(J,I,2) = CH(J,1,I) - WYK + 400 CONTINUE + RETURN + END + +C----------------------------------------------- ************ +C CRBLE1 +C ************ + SUBROUTINE SPDCRBLE1(NN,WORK) +C + DOUBLE PRECISION WORK(NN,2),P2,TWOPI + DATA TWOPI /6.2831853071795864769252867665590057683943388D0/ +C + N = 2 * NN + P2 = TWOPI/N + DO 10 I=1, NN + WORK(I,1) = DCOS(P2 * (I-1)) + WORK(I,2) = DSIN(P2 * (I-1)) + 10 CONTINUE + RETURN + END + +C----------------------------------------------- ************ +C CRFORM +C ************ + SUBROUTINE SPDCRFORM(IX,NS,NDIV2,CX,C,CH2) +C + DOUBLE COMPLEX CX(1),WYK1,C(NS,2),WYK + DOUBLE PRECISION CH2(NDIV2,2) +C + IF (IX .GT. 0) GOTO 50 + K = NS + 1 + DO 10 I=1, NS + WYK = CONJG(CX(NDIV2-I+2)) + C(I,1)= CX(I)+WYK + (CX(I) - WYK) * DCMPLX(CH2(I,2),CH2(I,1)) + WYK1 = CONJG(CX(NDIV2-K+2)) + C(I,2)= CX(K)+WYK1+ (CX(K) -WYK1) * DCMPLX(CH2(K,2),CH2(K,1)) + K = K + 1 + 10 CONTINUE + RETURN + 50 CONTINUE + K = NS + 1 + DO 20 I=1, NS + WYK = CONJG(CX(NDIV2-I+2)) + C(I,1)= CX(I)+WYK + (CX(I) - WYK) * DCMPLX(-CH2(I,2),CH2(I,1)) + WYK1 = CONJG(CX(NDIV2-K+2)) + C(I,2)= CX(K)+WYK1 +(CX(K) -WYK1) * DCMPLX(-CH2(K,2),CH2(K,1)) + K = K + 1 + 20 CONTINUE + RETURN + END + C------------------------------------------------------------- ************ C CRFFT2 C ************ @@ -62,36 +251,6 @@ C END C----------------------------------------------- ************ -C CRFORM -C ************ - SUBROUTINE SPDCRFORM(IX,NS,NDIV2,CX,C,CH2) -C - DOUBLE COMPLEX CX(1),WYK1,C(NS,2),WYK - DOUBLE PRECISION CH2(NDIV2,2) -C - IF (IX .GT. 0) GOTO 50 - K = NS + 1 - DO 10 I=1, NS - WYK = CONJG(CX(NDIV2-I+2)) - C(I,1)= CX(I)+WYK + (CX(I) - WYK) * DCMPLX(CH2(I,2),CH2(I,1)) - WYK1 = CONJG(CX(NDIV2-K+2)) - C(I,2)= CX(K)+WYK1+ (CX(K) -WYK1) * DCMPLX(CH2(K,2),CH2(K,1)) - K = K + 1 - 10 CONTINUE - RETURN - 50 CONTINUE - K = NS + 1 - DO 20 I=1, NS - WYK = CONJG(CX(NDIV2-I+2)) - C(I,1)= CX(I)+WYK + (CX(I) - WYK) * DCMPLX(-CH2(I,2),CH2(I,1)) - WYK1 = CONJG(CX(NDIV2-K+2)) - C(I,2)= CX(K)+WYK1 +(CX(K) -WYK1) * DCMPLX(-CH2(K,2),CH2(K,1)) - K = K + 1 - 20 CONTINUE - RETURN - END - -C----------------------------------------------- ************ C CROCK1 C ************ SUBROUTINE SPDCROCK1(NS,C,CH) @@ -157,23 +316,6 @@ C RETURN END -C----------------------------------------------- ************ -C CRBLE1 -C ************ - SUBROUTINE SPDCRBLE1(NN,WORK) -C - DOUBLE PRECISION WORK(NN,2),P2,TWOPI - DATA TWOPI /6.2831853071795864769252867665590057683943388D0/ -C - N = 2 * NN - P2 = TWOPI/N - DO 10 I=1, NN - WORK(I,1) = DCOS(P2 * (I-1)) - WORK(I,2) = DSIN(P2 * (I-1)) - 10 CONTINUE - RETURN - END - C------------------------------------------------------------- ************ C RCFFT2 C ************ @@ -236,32 +378,6 @@ C END C----------------------------------------------- ************ -C RTOCK2 -C ************ - SUBROUTINE SPDRTOCK2(LS,NS,C,CH,CH2) -C - DOUBLE COMPLEX WYK,C(NS,LS,2),CH(NS,2,LS) - DOUBLE PRECISION CH2(2,NS,LS,2) -C - IF (LS .GT. NS) GOTO 20 - DO 200 I=1, LS - DO 200 J=1, NS - WYK = DCMPLX(CH2(1,1,I,1),CH2(1,1,I,2)) * CH(J,2,I) - C(J,I,1) = CH(J,1,I) + WYK - C(J,I,2) = CH(J,1,I) - WYK - 200 CONTINUE - RETURN - 20 CONTINUE - DO 400 J=1, NS - DO 400 I=1, LS - WYK = DCMPLX(CH2(1,1,I,1),CH2(1,1,I,2)) * CH(J,2,I) - C(J,I,1) = CH(J,1,I) + WYK - C(J,I,2) = CH(J,1,I) - WYK - 400 CONTINUE - RETURN - END - -C----------------------------------------------- ************ C RTOCK1 C ************ SUBROUTINE SPDRTOCK1(NS,C,CH) @@ -275,105 +391,6 @@ C RETURN END -C----------------------------------------------- ************ -C RTOCK3 -C ************ - SUBROUTINE SPDRTOCK3(LS,NS,C,CH,CH2) -C - DOUBLE COMPLEX WYK,C(NS,LS,2),CH(NS,2,LS) - DOUBLE PRECISION CH2(2,NS,LS,2) -C - IF (LS .GT. NS) GOTO 30 - DO 600 I=1, LS - DO 600 J=1, NS - WYK =CONJG(DCMPLX(CH2(1,1,I,1),CH2(1,1,I,2))) * CH(J,2,I) - C(J,I,1) = CH(J,1,I) + WYK - C(J,I,2) = CH(J,1,I) - WYK - 600 CONTINUE - RETURN - 30 CONTINUE - DO 800 J=1, NS - DO 800 I=1, LS - WYK =CONJG(DCMPLX(CH2(1,1,I,1),CH2(1,1,I,2))) * CH(J,2,I) - C(J,I,1) = CH(J,1,I) + WYK - C(J,I,2) = CH(J,1,I) - WYK - 800 CONTINUE - RETURN - END - -C----------------------------------------------- ************ -C RABLE1 -C ************ - SUBROUTINE SPDRABLE1(NN,WORK) -C - DOUBLE PRECISION WORK(NN,2),P2,TWOPI - DATA TWOPI /6.2831853071795864769252867665590057683943388D0/ -C - N = 2 * NN - P2 = TWOPI/N - DO 10 I=1, NN - WORK(I,1) = DCOS(P2 * (I-1)) - WORK(I,2) = DSIN(P2 * (I-1)) - 10 CONTINUE - RETURN - END - -C----------------------------------------------- ************ -C RCONV1 -C ************ - SUBROUTINE SPDRCONV1(N,CY,C,CH) -C - DOUBLE COMPLEX CY(1) - DOUBLE PRECISION CH(N/2,2),P(2,1),C(2,N/2) - DOUBLE PRECISION X,Y,Z,Z1 -C - N2 = N/2 - P(1,1) = (C(1,1) + C(2,1)) * 2 - P(2,1) = (C(1,1) - C(2,1)) * 2 - CY(1) = DCMPLX(P(1,1),0.0D0) - CY(N2+1) = DCMPLX(P(2,1),0.0D0) - K = N2 - DO 10 I=2, N2 - X = C(1,I)+C(1,K) - Y = C(2,I)+C(2,K) - Z = C(1,I)-C(1,K) - Z1= C(2,I)-C(2,K) - P(1,1) = X + CH(I,1) * Y + CH(I,2) * Z - P(2,1) = Z1 + CH(I,2) * Y - CH(I,1) * Z - CY(I) = DCMPLX(P(1,1),P(2,1)) - K = K - 1 - 10 CONTINUE - RETURN - END - -C----------------------------------------------- ************ -C RCONV2 -C ************ - SUBROUTINE SPDRCONV2(N,CY,C,CH) -C - DOUBLE COMPLEX CY(1) - DOUBLE PRECISION CH(N/2,2),P(2,1),C(2,N/2) - DOUBLE PRECISION X,Y,Z,Z1 -C - N2 = N/2 - P(1,1) = (C(1,1) + C(2,1)) * 2 - P(2,1) = (C(1,1) - C(2,1)) * 2 - CY(1) = DCMPLX(P(1,1),0.0D0) - CY(N2+1) = DCMPLX(P(2,1),0.0D0) - K = N2 - DO 10 I=2, N2 - X = C(1,I)+C(1,K) - Y = C(2,I)+C(2,K) - Z = C(1,I)-C(1,K) - Z1= C(2,I)-C(2,K) - P(1,1) = X + CH(I,1) * Y - CH(I,2) * Z - P(2,1) = Z1 - CH(I,2) * Y - CH(I,1) * Z - CY(I) = DCMPLX(P(1,1),P(2,1)) - K = K - 1 - 10 CONTINUE - RETURN - END - C------------------------------------------------------------- ************ C CFFT2 C ************ @@ -516,23 +533,6 @@ C RETURN END -C----------------------------------------------- ************ -C CABLE2 -C ************ - SUBROUTINE SPDCABLE2(NN,WORK) -C - DOUBLE PRECISION WORK(2,NN),P2,TWOPI - DATA TWOPI /6.2831853071795864769252867665590057683943388D0/ -C - N = 2 * NN - P2 = TWOPI/N - DO 10 I=1, NN - WORK(1,I) = DCOS(P2 * (I-1)) - WORK(2,I) = DSIN(P2 * (I-1)) - 10 CONTINUE - RETURN - END - C------------------------------------------------------------- ************ C ABORT C ************