423 lines
14 KiB
Text
423 lines
14 KiB
Text
$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 ************
|