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