Line data Source code
1 : MODULE m_qsf
2 :
3 : ! interface qsf
4 : ! module procedure qsfReal, qsfComplex
5 : ! end interface
6 :
7 : CONTAINS
8 :
9 38325 : SUBROUTINE qsf(h,y,z,ndim,isave)
10 : c ..................................................................
11 : c subroutine qsf
12 : c purpose
13 : c to integrate an equidistant table of function values.
14 : c description of parameters
15 : c h - the increment of arguement values.
16 : c y - the input vector of function values. y must be
17 : c dimensioned ndim.
18 : c z - if isave = 0,z is the final value of the integration.
19 : c if isave not 0,z is the resulting vector of integrals
20 : c for each y value. in this case it must be dimensioned
21 : c ndim.
22 : c ndim - the number of points in the table.
23 : c isave - determines whether a single value or multiple results
24 : c will be returned
25 : c isave = 0 - only the final value of the integra-
26 : c tion will be returned.
27 : c isave ne 0 - the results of integration at each
28 : c point in the table will be returned.
29 : c remarks
30 : c no action in case ndim less than 3.
31 : c subroutines and function subprograms required
32 : c none
33 : c method
34 : c integration is done by means of simpsons rule together with
35 : c newtons 3/8 rule or a combination of these two rules. trunca-
36 : c tion error is of order h**5 if ndim greater than 3. if ndim = 3
37 : c truncation error is of order h**4.
38 : c reference
39 : c (1) f.b.hildebrand, introduction to numerical analysis,
40 : c mcgraw-hill, new york/toronto/london, 1956, pp.71-76.
41 : c (2) r. zurmuehl, praktische mathematic fuer ingenieure und
42 : c physiker, springer, berlin/goettingen/heidelberg, 1963,
43 : c pp. 214-221.
44 : c ..................................................................
45 : C .. Scalar Arguments ..
46 :
47 : IMPLICIT NONE
48 : REAL h
49 : INTEGER isave,ndim
50 : C ..
51 : C .. Array Arguments ..
52 : REAL y(ndim),z(*)
53 : C ..
54 : C .. Local Scalars ..
55 : REAL aux,aux1,aux2,ht,sum1,sum2,val
56 : INTEGER i
57 : LOGICAL save
58 : C ..
59 38325 : save = isave .NE. 0
60 38325 : ht = h / 3.0
61 38325 : IF ( ndim == 5 ) THEN
62 : GOTO 140
63 38325 : ELSEIF ( ndim < 5 ) THEN
64 : GOTO 130
65 : ENDIF
66 : c ndim is greater than 5. preparations of integration loop
67 38325 : sum1 = y(2) + y(2)
68 38325 : sum1 = sum1 + sum1
69 38325 : sum1 = ht* (y(1)+sum1+y(3))
70 38325 : aux1 = y(4) + y(4)
71 38325 : aux1 = aux1 + aux1
72 38325 : aux1 = sum1 + ht* (y(3)+aux1+y(5))
73 38325 : aux2 = ht* (y(1)+3.875e0* (y(2)+y(5))+2.625e0* (y(3)+y(4))+y(6))
74 38325 : IF (.NOT.save) GO TO 30
75 38074 : 20 z(1) = 0.e0
76 38074 : sum2 = y(5) + y(5)
77 38074 : sum2 = sum2 + sum2
78 38074 : sum2 = aux2 - ht* (y(4)+sum2+y(6))
79 38074 : aux = y(3) + y(3)
80 38074 : aux = aux + aux
81 38074 : z(2) = sum2 - ht* (y(2)+aux+y(4))
82 38074 : z(3) = sum1
83 38325 : z(4) = sum2
84 : 30 CONTINUE
85 38325 : IF ( ndim <= 6 ) THEN
86 : GOTO 70
87 : ENDIF
88 : c integration loop
89 1866525 : 40 DO 60 i = 7,ndim,2
90 1828200 : sum1 = aux1
91 1828200 : sum2 = aux2
92 1828200 : aux1 = y(i-1) + y(i-1)
93 1828200 : aux1 = aux1 + aux1
94 1828200 : aux1 = sum1 + ht* (y(i-2)+aux1+y(i))
95 1828200 : IF (save) z(i-2) = sum1
96 1828200 : IF ( i >= ndim ) THEN
97 : GOTO 100
98 : ENDIF
99 1828200 : aux2 = y(i) + y(i)
100 1828200 : aux2 = aux2 + aux2
101 1828200 : aux2 = sum2 + ht* (y(i-1)+aux2+y(i+1))
102 1828200 : IF (save) z(i-1) = sum2
103 38325 : 60 CONTINUE
104 38325 : 70 IF (.NOT.save) GO TO 90
105 38074 : 80 z(ndim-1) = aux1
106 38074 : z(ndim) = aux2
107 38325 : RETURN
108 251 : 90 z(1) = aux2
109 251 : RETURN
110 0 : 100 IF (.NOT.save) GO TO 120
111 0 : 110 z(ndim-1) = sum2
112 0 : z(ndim) = aux1
113 0 : RETURN
114 0 : 120 z(1) = aux1
115 0 : RETURN
116 : c end of integration loop
117 : 130 CONTINUE
118 0 : IF ( ndim == 3 ) THEN
119 : GOTO 210
120 0 : ELSEIF ( ndim < 3 ) THEN
121 : GOTO 230
122 : ENDIF
123 : c ndim is equal to 4 or 5
124 0 : 140 sum2 = 1.125e0*ht* (y(1)+y(2)+y(2)+y(2)+y(3)+y(3)+y(3)+y(4))
125 0 : sum1 = y(2) + y(2)
126 0 : sum1 = sum1 + sum1
127 0 : sum1 = ht* (y(1)+sum1+y(3))
128 0 : IF (.NOT.save) GO TO 160
129 0 : 150 z(1) = 0.e0
130 0 : aux1 = y(3) + y(3)
131 0 : aux1 = aux1 + aux1
132 0 : z(2) = sum2 - ht* (y(2)+aux1+y(4))
133 : 160 CONTINUE
134 0 : IF ( ndim < 5 ) THEN
135 : GOTO 190
136 : ENDIF
137 0 : 170 aux1 = y(4) + y(4)
138 0 : aux1 = aux1 + aux1
139 0 : val = sum1 + ht* (y(3)+aux1+y(5))
140 0 : IF (save) GO TO 180
141 0 : z(1) = val
142 0 : RETURN
143 0 : 180 z(5) = val
144 0 : GO TO 200
145 0 : 190 IF (save) GO TO 200
146 0 : z(1) = sum2
147 0 : RETURN
148 0 : 200 z(3) = sum1
149 0 : z(4) = sum2
150 0 : RETURN
151 : c ndim is equal to 3
152 0 : 210 sum2 = y(2) + y(2)
153 0 : sum2 = sum2 + sum2
154 0 : val = ht* (y(1)+sum2+y(3))
155 0 : IF (save) GO TO 220
156 0 : z(1) = val
157 0 : RETURN
158 0 : 220 z(1) = 0.e0
159 0 : z(2) = ht* (1.25e0*y(1)+y(2)+y(2)-.25e0*y(3))
160 0 : z(3) = val
161 : 230 RETURN
162 : END SUBROUTINE qsf
163 :
164 0 : subroutine qsfComplex( h, y, z, ndim, isave )
165 :
166 : implicit none
167 :
168 : REAL h
169 : INTEGER isave,ndim
170 : COMPLEX y(ndim),z(*)
171 :
172 0 : real zr(ndim), zi(ndim)
173 :
174 0 : call qsf( h, real(y), zr, ndim, isave )
175 0 : call qsf( h, aimag(y), zi, ndim, isave )
176 0 : z(1:ndim) = cmplx( zr(1:ndim), zi(1:ndim) )
177 :
178 0 : end subroutine qsfComplex
179 :
180 : END MODULE m_qsf
|
