182
183
184
185
186
187
188 LOGICAL IEEE
189 INTEGER I0, ITER, N0, NDIV, NFAIL, PP
190 REAL DESIG, DMIN, DMIN1, DMIN2, DN, DN1, DN2, G,
191 $ QMAX, SIGMA, TAU
192
193
194 REAL Z( * )
195
196
197
198
199
200 REAL CBIAS
201 parameter( cbias = 1.50e0 )
202 REAL ZERO, QURTR, HALF, ONE, TWO, HUNDRD
203 parameter( zero = 0.0e0, qurtr = 0.250e0, half = 0.5e0,
204 $ one = 1.0e0, two = 2.0e0, hundrd = 100.0e0 )
205
206
207 INTEGER IPN4, J4, N0IN, NN, TTYPE
208 REAL EPS, S, T, TEMP, TOL, TOL2
209
210
212
213
214 REAL SLAMCH
215 LOGICAL SISNAN
217
218
219 INTRINSIC abs, max, min, sqrt
220
221
222
223 n0in = n0
224 eps =
slamch(
'Precision' )
225 tol = eps*hundrd
226 tol2 = tol**2
227
228
229
230 10 CONTINUE
231
232 IF( n0.LT.i0 )
233 $ RETURN
234 IF( n0.EQ.i0 )
235 $ GO TO 20
236 nn = 4*n0 + pp
237 IF( n0.EQ.( i0+1 ) )
238 $ GO TO 40
239
240
241
242 IF( z( nn-5 ).GT.tol2*( sigma+z( nn-3 ) ) .AND.
243 $ z( nn-2*pp-4 ).GT.tol2*z( nn-7 ) )
244 $ GO TO 30
245
246 20 CONTINUE
247
248 z( 4*n0-3 ) = z( 4*n0+pp-3 ) + sigma
249 n0 = n0 - 1
250 GO TO 10
251
252
253
254 30 CONTINUE
255
256 IF( z( nn-9 ).GT.tol2*sigma .AND.
257 $ z( nn-2*pp-8 ).GT.tol2*z( nn-11 ) )
258 $ GO TO 50
259
260 40 CONTINUE
261
262 IF( z( nn-3 ).GT.z( nn-7 ) ) THEN
263 s = z( nn-3 )
264 z( nn-3 ) = z( nn-7 )
265 z( nn-7 ) = s
266 END IF
267 t = half*( ( z( nn-7 )-z( nn-3 ) )+z( nn-5 ) )
268 IF( z( nn-5 ).GT.z( nn-3 )*tol2.AND.t.NE.zero ) THEN
269 s = z( nn-3 )*( z( nn-5 ) / t )
270 IF( s.LE.t ) THEN
271 s = z( nn-3 )*( z( nn-5 ) /
272 $ ( t*( one+sqrt( one+s / t ) ) ) )
273 ELSE
274 s = z( nn-3 )*( z( nn-5 ) / ( t+sqrt( t )*sqrt( t+s ) ) )
275 END IF
276 t = z( nn-7 ) + ( s+z( nn-5 ) )
277 z( nn-3 ) = z( nn-3 )*( z( nn-7 ) / t )
278 z( nn-7 ) = t
279 END IF
280 z( 4*n0-7 ) = z( nn-7 ) + sigma
281 z( 4*n0-3 ) = z( nn-3 ) + sigma
282 n0 = n0 - 2
283 GO TO 10
284
285 50 CONTINUE
286 IF( pp.EQ.2 )
287 $ pp = 0
288
289
290
291 IF( dmin.LE.zero .OR. n0.LT.n0in ) THEN
292 IF( cbias*z( 4*i0+pp-3 ).LT.z( 4*n0+pp-3 ) ) THEN
293 ipn4 = 4*( i0+n0 )
294 DO 60 j4 = 4*i0, 2*( i0+n0-1 ), 4
295 temp = z( j4-3 )
296 z( j4-3 ) = z( ipn4-j4-3 )
297 z( ipn4-j4-3 ) = temp
298 temp = z( j4-2 )
299 z( j4-2 ) = z( ipn4-j4-2 )
300 z( ipn4-j4-2 ) = temp
301 temp = z( j4-1 )
302 z( j4-1 ) = z( ipn4-j4-5 )
303 z( ipn4-j4-5 ) = temp
304 temp = z( j4 )
305 z( j4 ) = z( ipn4-j4-4 )
306 z( ipn4-j4-4 ) = temp
307 60 CONTINUE
308 IF( n0-i0.LE.4 ) THEN
309 z( 4*n0+pp-1 ) = z( 4*i0+pp-1 )
310 z( 4*n0-pp ) = z( 4*i0-pp )
311 END IF
312 dmin2 = min( dmin2, z( 4*n0+pp-1 ) )
313 z( 4*n0+pp-1 ) = min( z( 4*n0+pp-1 ), z( 4*i0+pp-1 ),
314 $ z( 4*i0+pp+3 ) )
315 z( 4*n0-pp ) = min( z( 4*n0-pp ), z( 4*i0-pp ),
316 $ z( 4*i0-pp+4 ) )
317 qmax = max( qmax, z( 4*i0+pp-3 ), z( 4*i0+pp+1 ) )
318 dmin = -zero
319 END IF
320 END IF
321
322
323
324 CALL slasq4( i0, n0, z, pp, n0in, dmin, dmin1, dmin2, dn, dn1,
325 $ dn2, tau, ttype, g )
326
327
328
329 70 CONTINUE
330
331 CALL slasq5( i0, n0, z, pp, tau, sigma, dmin, dmin1, dmin2, dn,
332 $ dn1, dn2, ieee, eps )
333
334 ndiv = ndiv + ( n0-i0+2 )
335 iter = iter + 1
336
337
338
339 IF( dmin.GE.zero .AND. dmin1.GE.zero ) THEN
340
341
342
343 GO TO 90
344
345 ELSE IF( dmin.LT.zero .AND. dmin1.GT.zero .AND.
346 $ z( 4*( n0-1 )-pp ).LT.tol*( sigma+dn1 ) .AND.
347 $ abs( dn ).LT.tol*sigma ) THEN
348
349
350
351 z( 4*( n0-1 )-pp+2 ) = zero
352 dmin = zero
353 GO TO 90
354 ELSE IF( dmin.LT.zero ) THEN
355
356
357
358 nfail = nfail + 1
359 IF( ttype.LT.-22 ) THEN
360
361
362
363 tau = zero
364 ELSE IF( dmin1.GT.zero ) THEN
365
366
367
368 tau = ( tau+dmin )*( one-two*eps )
369 ttype = ttype - 11
370 ELSE
371
372
373
374 tau = qurtr*tau
375 ttype = ttype - 12
376 END IF
377 GO TO 70
378 ELSE IF(
sisnan( dmin ) )
THEN
379
380
381
382 IF( tau.EQ.zero ) THEN
383 GO TO 80
384 ELSE
385 tau = zero
386 GO TO 70
387 END IF
388 ELSE
389
390
391
392 GO TO 80
393 END IF
394
395
396
397 80 CONTINUE
398 CALL slasq6( i0, n0, z, pp, dmin, dmin1, dmin2, dn, dn1, dn2 )
399 ndiv = ndiv + ( n0-i0+2 )
400 iter = iter + 1
401 tau = zero
402
403 90 CONTINUE
404 IF( tau.LT.sigma ) THEN
405 desig = desig + tau
406 t = sigma + desig
407 desig = desig - ( t-sigma )
408 ELSE
409 t = sigma + tau
410 desig = sigma - ( t-tau ) + desig
411 END IF
412 sigma = t
413
414 RETURN
415
416
417
logical function sisnan(sin)
SISNAN tests input for NaN.
real function slamch(cmach)
SLAMCH
subroutine slasq4(i0, n0, z, pp, n0in, dmin, dmin1, dmin2, dn, dn1, dn2, tau, ttype, g)
SLASQ4 computes an approximation to the smallest eigenvalue using values of d from the previous trans...
subroutine slasq5(i0, n0, z, pp, tau, sigma, dmin, dmin1, dmin2, dn, dnm1, dnm2, ieee, eps)
SLASQ5 computes one dqds transform in ping-pong form. Used by sbdsqr and sstegr.
subroutine slasq6(i0, n0, z, pp, dmin, dmin1, dmin2, dn, dnm1, dnm2)
SLASQ6 computes one dqd transform in ping-pong form. Used by sbdsqr and sstegr.