273 SUBROUTINE slaqr2( WANTT, WANTZ, N, KTOP, KBOT, NW, H, LDH,
275 $ IHIZ, Z, LDZ, NS, ND, SR, SI, V, LDV, NH, T,
276 $ LDT, NV, WV, LDWV, WORK, LWORK )
283 INTEGER IHIZ, ILOZ, KBOT, KTOP, LDH, LDT, LDV, LDWV,
284 $ LDZ, LWORK, N, ND, NH, NS, NV, NW
288 REAL H( LDH, * ), SI( * ), SR( * ), T( LDT, * ),
289 $ V( LDV, * ), WORK( * ), WV( LDWV, * ),
296 PARAMETER ( ZERO = 0.0e0, one = 1.0e0 )
299 REAL AA, BB, CC, CS, DD, EVI, EVK, FOO, S,
300 $ SAFMAX, SAFMIN, SMLNUM, SN, TAU, ULP
301 INTEGER I, IFST, ILST, INFO, INFQR, J, JW, K, KCOL,
302 $ kend, kln, krow, kwtop, ltop, lwk1, lwk2,
304 LOGICAL BULGE, SORTED
307 REAL SLAMCH, SROUNDUP_LWORK
308 EXTERNAL SLAMCH, SROUNDUP_LWORK
317 INTRINSIC abs, int, max, min, real, sqrt
323 jw = min( nw, kbot-ktop+1 )
330 CALL sgehrd( jw, 1, jw-1, t, ldt, work, work, -1, info )
331 lwk1 = int( work( 1 ) )
335 CALL sormhr(
'R',
'N', jw, jw, 1, jw-1, t, ldt, work, v,
338 lwk2 = int( work( 1 ) )
342 lwkopt = jw + max( lwk1, lwk2 )
347 IF( lwork.EQ.-1 )
THEN
348 work( 1 ) = sroundup_lwork( lwkopt )
365 safmin = slamch(
'SAFE MINIMUM' )
366 safmax = one / safmin
367 ulp = slamch(
'PRECISION' )
368 smlnum = safmin*( real( n ) / ulp )
372 jw = min( nw, kbot-ktop+1 )
373 kwtop = kbot - jw + 1
374 IF( kwtop.EQ.ktop )
THEN
377 s = h( kwtop, kwtop-1 )
380 IF( kbot.EQ.kwtop )
THEN
384 sr( kwtop ) = h( kwtop, kwtop )
388 IF( abs( s ).LE.max( smlnum, ulp*abs( h( kwtop, kwtop ) ) ) )
393 $ h( kwtop, kwtop-1 ) = zero
405 CALL slacpy(
'U', jw, jw, h( kwtop, kwtop ), ldh, t, ldt )
406 CALL scopy( jw-1, h( kwtop+1, kwtop ), ldh+1, t( 2, 1 ),
409 CALL slaset(
'A', jw, jw, zero, one, v, ldv )
410 CALL slahqr( .true., .true., jw, 1, jw, t, ldt, sr( kwtop ),
411 $ si( kwtop ), 1, jw, v, ldv, infqr )
420 $ t( jw, jw-2 ) = zero
427 IF( ilst.LE.ns )
THEN
431 bulge = t( ns, ns-1 ).NE.zero
436 IF( .NOT.bulge )
THEN
440 foo = abs( t( ns, ns ) )
443 IF( abs( s*v( 1, ns ) ).LE.max( smlnum, ulp*foo ) )
THEN
454 CALL strexc(
'V', jw, t, ldt, v, ldv, ifst, ilst,
463 foo = abs( t( ns, ns ) ) + sqrt( abs( t( ns, ns-1 ) ) )*
464 $ sqrt( abs( t( ns-1, ns ) ) )
467 IF( max( abs( s*v( 1, ns ) ), abs( s*v( 1, ns-1 ) ) ).LE.
468 $ max( smlnum, ulp*foo ) )
THEN
480 CALL strexc(
'V', jw, t, ldt, v, ldv, ifst, ilst,
514 ELSE IF( t( i+1, i ).EQ.zero )
THEN
522 evi = abs( t( i, i ) )
524 evi = abs( t( i, i ) ) + sqrt( abs( t( i+1, i ) ) )*
525 $ sqrt( abs( t( i, i+1 ) ) )
529 evk = abs( t( k, k ) )
530 ELSE IF( t( k+1, k ).EQ.zero )
THEN
531 evk = abs( t( k, k ) )
533 evk = abs( t( k, k ) ) + sqrt( abs( t( k+1, k ) ) )*
534 $ sqrt( abs( t( k, k+1 ) ) )
537 IF( evi.GE.evk )
THEN
543 CALL strexc(
'V', jw, t, ldt, v, ldv, ifst, ilst,
554 ELSE IF( t( i+1, i ).EQ.zero )
THEN
569 IF( i.GE.infqr+1 )
THEN
570 IF( i.EQ.infqr+1 )
THEN
571 sr( kwtop+i-1 ) = t( i, i )
572 si( kwtop+i-1 ) = zero
574 ELSE IF( t( i, i-1 ).EQ.zero )
THEN
575 sr( kwtop+i-1 ) = t( i, i )
576 si( kwtop+i-1 ) = zero
583 CALL slanv2( aa, bb, cc, dd, sr( kwtop+i-2 ),
584 $ si( kwtop+i-2 ), sr( kwtop+i-1 ),
585 $ si( kwtop+i-1 ), cs, sn )
591 IF( ns.LT.jw .OR. s.EQ.zero )
THEN
592 IF( ns.GT.1 .AND. s.NE.zero )
THEN
596 CALL scopy( ns, v, ldv, work, 1 )
597 CALL slarfg( ns, work( 1 ), work( 2 ), 1, tau )
599 CALL slaset(
'L', jw-2, jw-2, zero, zero, t( 3, 1 ),
602 CALL slarf1f(
'L', ns, jw, work, 1, tau, t, ldt,
604 CALL slarf1f(
'R', ns, ns, work, 1, tau, t, ldt,
606 CALL slarf1f(
'R', jw, ns, work, 1, tau, v, ldv,
609 CALL sgehrd( jw, 1, ns, t, ldt, work, work( jw+1 ),
616 $ h( kwtop, kwtop-1 ) = s*v( 1, 1 )
617 CALL slacpy(
'U', jw, jw, t, ldt, h( kwtop, kwtop ), ldh )
618 CALL scopy( jw-1, t( 2, 1 ), ldt+1, h( kwtop+1, kwtop ),
624 IF( ns.GT.1 .AND. s.NE.zero )
625 $
CALL sormhr(
'R',
'N', jw, ns, 1, ns, t, ldt, work, v,
627 $ work( jw+1 ), lwork-jw, info )
636 DO 70 krow = ltop, kwtop - 1, nv
637 kln = min( nv, kwtop-krow )
638 CALL sgemm(
'N',
'N', kln, jw, jw, one, h( krow, kwtop ),
639 $ ldh, v, ldv, zero, wv, ldwv )
640 CALL slacpy(
'A', kln, jw, wv, ldwv, h( krow, kwtop ),
647 DO 80 kcol = kbot + 1, n, nh
648 kln = min( nh, n-kcol+1 )
649 CALL sgemm(
'C',
'N', jw, kln, jw, one, v, ldv,
650 $ h( kwtop, kcol ), ldh, zero, t, ldt )
651 CALL slacpy(
'A', jw, kln, t, ldt, h( kwtop, kcol ),
659 DO 90 krow = iloz, ihiz, nv
660 kln = min( nv, ihiz-krow+1 )
661 CALL sgemm(
'N',
'N', kln, jw, jw, one, z( krow,
663 $ ldz, v, ldv, zero, wv, ldwv )
664 CALL slacpy(
'A', kln, jw, wv, ldwv, z( krow, kwtop ),
684 work( 1 ) = sroundup_lwork( lwkopt )