--- /dev/null
+*DECK RADFG
+ SUBROUTINE RADFG (IDO, IP, L1, IDL1, CC, C1, C2, CH, CH2, WA)
+C***BEGIN PROLOGUE RADFG
+C***SUBSIDIARY
+C***PURPOSE Calculate the fast Fourier transform of subvectors of
+C arbitrary length.
+C***LIBRARY SLATEC (FFTPACK)
+C***TYPE SINGLE PRECISION (RADFG-S)
+C***AUTHOR Swarztrauber, P. N., (NCAR)
+C***ROUTINES CALLED (NONE)
+C***REVISION HISTORY (YYMMDD)
+C 790601 DATE WRITTEN
+C 830401 Modified to use SLATEC library source file format.
+C 860115 Modified by Ron Boisvert to adhere to Fortran 77 by
+C (a) changing dummy array size declarations (1) to (*),
+C (b) changing references to intrinsic function FLOAT
+C to REAL, and
+C (c) changing definition of variable TPI by using
+C FORTRAN intrinsic function ATAN instead of a DATA
+C statement.
+C 881128 Modified by Dick Valent to meet prologue standards.
+C 890531 Changed all specific intrinsics to generic. (WRB)
+C 890831 Modified array declarations. (WRB)
+C 891214 Prologue converted to Version 4.0 format. (BAB)
+C 900402 Added TYPE section. (WRB)
+C***END PROLOGUE RADFG
+ DIMENSION CH(IDO,L1,*), CC(IDO,IP,*), C1(IDO,L1,*),
+ + C2(IDL1,*), CH2(IDL1,*), WA(*)
+C***FIRST EXECUTABLE STATEMENT RADFG
+ TPI = 8.*ATAN(1.)
+ ARG = TPI/IP
+ DCP = COS(ARG)
+ DSP = SIN(ARG)
+ IPPH = (IP+1)/2
+ IPP2 = IP+2
+ IDP2 = IDO+2
+ NBD = (IDO-1)/2
+ IF (IDO .EQ. 1) GO TO 119
+ DO 101 IK=1,IDL1
+ CH2(IK,1) = C2(IK,1)
+ 101 CONTINUE
+ DO 103 J=2,IP
+ DO 102 K=1,L1
+ CH(1,K,J) = C1(1,K,J)
+ 102 CONTINUE
+ 103 CONTINUE
+ IF (NBD .GT. L1) GO TO 107
+ IS = -IDO
+ DO 106 J=2,IP
+ IS = IS+IDO
+ IDIJ = IS
+ DO 105 I=3,IDO,2
+ IDIJ = IDIJ+2
+ DO 104 K=1,L1
+ CH(I-1,K,J) = WA(IDIJ-1)*C1(I-1,K,J)+WA(IDIJ)*C1(I,K,J)
+ CH(I,K,J) = WA(IDIJ-1)*C1(I,K,J)-WA(IDIJ)*C1(I-1,K,J)
+ 104 CONTINUE
+ 105 CONTINUE
+ 106 CONTINUE
+ GO TO 111
+ 107 IS = -IDO
+ DO 110 J=2,IP
+ IS = IS+IDO
+ DO 109 K=1,L1
+ IDIJ = IS
+CDIR$ IVDEP
+ DO 108 I=3,IDO,2
+ IDIJ = IDIJ+2
+ CH(I-1,K,J) = WA(IDIJ-1)*C1(I-1,K,J)+WA(IDIJ)*C1(I,K,J)
+ CH(I,K,J) = WA(IDIJ-1)*C1(I,K,J)-WA(IDIJ)*C1(I-1,K,J)
+ 108 CONTINUE
+ 109 CONTINUE
+ 110 CONTINUE
+ 111 IF (NBD .LT. L1) GO TO 115
+ DO 114 J=2,IPPH
+ JC = IPP2-J
+ DO 113 K=1,L1
+CDIR$ IVDEP
+ DO 112 I=3,IDO,2
+ C1(I-1,K,J) = CH(I-1,K,J)+CH(I-1,K,JC)
+ C1(I-1,K,JC) = CH(I,K,J)-CH(I,K,JC)
+ C1(I,K,J) = CH(I,K,J)+CH(I,K,JC)
+ C1(I,K,JC) = CH(I-1,K,JC)-CH(I-1,K,J)
+ 112 CONTINUE
+ 113 CONTINUE
+ 114 CONTINUE
+ GO TO 121
+ 115 DO 118 J=2,IPPH
+ JC = IPP2-J
+ DO 117 I=3,IDO,2
+ DO 116 K=1,L1
+ C1(I-1,K,J) = CH(I-1,K,J)+CH(I-1,K,JC)
+ C1(I-1,K,JC) = CH(I,K,J)-CH(I,K,JC)
+ C1(I,K,J) = CH(I,K,J)+CH(I,K,JC)
+ C1(I,K,JC) = CH(I-1,K,JC)-CH(I-1,K,J)
+ 116 CONTINUE
+ 117 CONTINUE
+ 118 CONTINUE
+ GO TO 121
+ 119 DO 120 IK=1,IDL1
+ C2(IK,1) = CH2(IK,1)
+ 120 CONTINUE
+ 121 DO 123 J=2,IPPH
+ JC = IPP2-J
+ DO 122 K=1,L1
+ C1(1,K,J) = CH(1,K,J)+CH(1,K,JC)
+ C1(1,K,JC) = CH(1,K,JC)-CH(1,K,J)
+ 122 CONTINUE
+ 123 CONTINUE
+C
+ AR1 = 1.
+ AI1 = 0.
+ DO 127 L=2,IPPH
+ LC = IPP2-L
+ AR1H = DCP*AR1-DSP*AI1
+ AI1 = DCP*AI1+DSP*AR1
+ AR1 = AR1H
+ DO 124 IK=1,IDL1
+ CH2(IK,L) = C2(IK,1)+AR1*C2(IK,2)
+ CH2(IK,LC) = AI1*C2(IK,IP)
+ 124 CONTINUE
+ DC2 = AR1
+ DS2 = AI1
+ AR2 = AR1
+ AI2 = AI1
+ DO 126 J=3,IPPH
+ JC = IPP2-J
+ AR2H = DC2*AR2-DS2*AI2
+ AI2 = DC2*AI2+DS2*AR2
+ AR2 = AR2H
+ DO 125 IK=1,IDL1
+ CH2(IK,L) = CH2(IK,L)+AR2*C2(IK,J)
+ CH2(IK,LC) = CH2(IK,LC)+AI2*C2(IK,JC)
+ 125 CONTINUE
+ 126 CONTINUE
+ 127 CONTINUE
+ DO 129 J=2,IPPH
+ DO 128 IK=1,IDL1
+ CH2(IK,1) = CH2(IK,1)+C2(IK,J)
+ 128 CONTINUE
+ 129 CONTINUE
+C
+ IF (IDO .LT. L1) GO TO 132
+ DO 131 K=1,L1
+ DO 130 I=1,IDO
+ CC(I,1,K) = CH(I,K,1)
+ 130 CONTINUE
+ 131 CONTINUE
+ GO TO 135
+ 132 DO 134 I=1,IDO
+ DO 133 K=1,L1
+ CC(I,1,K) = CH(I,K,1)
+ 133 CONTINUE
+ 134 CONTINUE
+ 135 DO 137 J=2,IPPH
+ JC = IPP2-J
+ J2 = J+J
+ DO 136 K=1,L1
+ CC(IDO,J2-2,K) = CH(1,K,J)
+ CC(1,J2-1,K) = CH(1,K,JC)
+ 136 CONTINUE
+ 137 CONTINUE
+ IF (IDO .EQ. 1) RETURN
+ IF (NBD .LT. L1) GO TO 141
+ DO 140 J=2,IPPH
+ JC = IPP2-J
+ J2 = J+J
+ DO 139 K=1,L1
+CDIR$ IVDEP
+ DO 138 I=3,IDO,2
+ IC = IDP2-I
+ CC(I-1,J2-1,K) = CH(I-1,K,J)+CH(I-1,K,JC)
+ CC(IC-1,J2-2,K) = CH(I-1,K,J)-CH(I-1,K,JC)
+ CC(I,J2-1,K) = CH(I,K,J)+CH(I,K,JC)
+ CC(IC,J2-2,K) = CH(I,K,JC)-CH(I,K,J)
+ 138 CONTINUE
+ 139 CONTINUE
+ 140 CONTINUE
+ RETURN
+ 141 DO 144 J=2,IPPH
+ JC = IPP2-J
+ J2 = J+J
+ DO 143 I=3,IDO,2
+ IC = IDP2-I
+ DO 142 K=1,L1
+ CC(I-1,J2-1,K) = CH(I-1,K,J)+CH(I-1,K,JC)
+ CC(IC-1,J2-2,K) = CH(I-1,K,J)-CH(I-1,K,JC)
+ CC(I,J2-1,K) = CH(I,K,J)+CH(I,K,JC)
+ CC(IC,J2-2,K) = CH(I,K,JC)-CH(I,K,J)
+ 142 CONTINUE
+ 143 CONTINUE
+ 144 CONTINUE
+ RETURN
+ END