117
118
119
120
121
122
123 INTEGER I0, N0, PP
124 DOUBLE PRECISION DMIN, DMIN1, DMIN2, DN, DNM1, DNM2
125
126
127 DOUBLE PRECISION Z( * )
128
129
130
131
132
133 DOUBLE PRECISION ZERO
134 parameter( zero = 0.0d0 )
135
136
137 INTEGER J4, J4P2
138 DOUBLE PRECISION D, EMIN, SAFMIN, TEMP
139
140
141 DOUBLE PRECISION DLAMCH
143
144
145 INTRINSIC min
146
147
148
149 IF( ( n0-i0-1 ).LE.0 )
150 $ RETURN
151
152 safmin =
dlamch(
'Safe minimum' )
153 j4 = 4*i0 + pp - 3
154 emin = z( j4+4 )
155 d = z( j4 )
156 dmin = d
157
158 IF( pp.EQ.0 ) THEN
159 DO 10 j4 = 4*i0, 4*( n0-3 ), 4
160 z( j4-2 ) = d + z( j4-1 )
161 IF( z( j4-2 ).EQ.zero ) THEN
162 z( j4 ) = zero
163 d = z( j4+1 )
164 dmin = d
165 emin = zero
166 ELSE IF( safmin*z( j4+1 ).LT.z( j4-2 ) .AND.
167 $ safmin*z( j4-2 ).LT.z( j4+1 ) ) THEN
168 temp = z( j4+1 ) / z( j4-2 )
169 z( j4 ) = z( j4-1 )*temp
170 d = d*temp
171 ELSE
172 z( j4 ) = z( j4+1 )*( z( j4-1 ) / z( j4-2 ) )
173 d = z( j4+1 )*( d / z( j4-2 ) )
174 END IF
175 dmin = min( dmin, d )
176 emin = min( emin, z( j4 ) )
177 10 CONTINUE
178 ELSE
179 DO 20 j4 = 4*i0, 4*( n0-3 ), 4
180 z( j4-3 ) = d + z( j4 )
181 IF( z( j4-3 ).EQ.zero ) THEN
182 z( j4-1 ) = zero
183 d = z( j4+2 )
184 dmin = d
185 emin = zero
186 ELSE IF( safmin*z( j4+2 ).LT.z( j4-3 ) .AND.
187 $ safmin*z( j4-3 ).LT.z( j4+2 ) ) THEN
188 temp = z( j4+2 ) / z( j4-3 )
189 z( j4-1 ) = z( j4 )*temp
190 d = d*temp
191 ELSE
192 z( j4-1 ) = z( j4+2 )*( z( j4 ) / z( j4-3 ) )
193 d = z( j4+2 )*( d / z( j4-3 ) )
194 END IF
195 dmin = min( dmin, d )
196 emin = min( emin, z( j4-1 ) )
197 20 CONTINUE
198 END IF
199
200
201
202 dnm2 = d
203 dmin2 = dmin
204 j4 = 4*( n0-2 ) - pp
205 j4p2 = j4 + 2*pp - 1
206 z( j4-2 ) = dnm2 + z( j4p2 )
207 IF( z( j4-2 ).EQ.zero ) THEN
208 z( j4 ) = zero
209 dnm1 = z( j4p2+2 )
210 dmin = dnm1
211 emin = zero
212 ELSE IF( safmin*z( j4p2+2 ).LT.z( j4-2 ) .AND.
213 $ safmin*z( j4-2 ).LT.z( j4p2+2 ) ) THEN
214 temp = z( j4p2+2 ) / z( j4-2 )
215 z( j4 ) = z( j4p2 )*temp
216 dnm1 = dnm2*temp
217 ELSE
218 z( j4 ) = z( j4p2+2 )*( z( j4p2 ) / z( j4-2 ) )
219 dnm1 = z( j4p2+2 )*( dnm2 / z( j4-2 ) )
220 END IF
221 dmin = min( dmin, dnm1 )
222
223 dmin1 = dmin
224 j4 = j4 + 4
225 j4p2 = j4 + 2*pp - 1
226 z( j4-2 ) = dnm1 + z( j4p2 )
227 IF( z( j4-2 ).EQ.zero ) THEN
228 z( j4 ) = zero
229 dn = z( j4p2+2 )
230 dmin = dn
231 emin = zero
232 ELSE IF( safmin*z( j4p2+2 ).LT.z( j4-2 ) .AND.
233 $ safmin*z( j4-2 ).LT.z( j4p2+2 ) ) THEN
234 temp = z( j4p2+2 ) / z( j4-2 )
235 z( j4 ) = z( j4p2 )*temp
236 dn = dnm1*temp
237 ELSE
238 z( j4 ) = z( j4p2+2 )*( z( j4p2 ) / z( j4-2 ) )
239 dn = z( j4p2+2 )*( dnm1 / z( j4-2 ) )
240 END IF
241 dmin = min( dmin, dn )
242
243 z( j4+2 ) = dn
244 z( 4*n0-pp ) = emin
245 RETURN
246
247
248
double precision function dlamch(cmach)
DLAMCH