232 RECURSIVE SUBROUTINE slaqz3( ILSCHUR, ILQ, ILZ, N, ILO, IHI,
234 $ A, LDA, B, LDB, Q, LDQ, Z, LDZ, NS,
235 $ ND, ALPHAR, ALPHAI, BETA, QC, LDQC,
236 $ ZC, LDZC, WORK, LWORK, REC, INFO )
240 LOGICAL,
INTENT( IN ) :: ilschur, ilq, ilz
241 INTEGER,
INTENT( IN ) :: n, ilo, ihi, nw, lda, ldb, ldq, ldz,
242 $ ldqc, ldzc, lwork, rec
244 REAL,
INTENT( INOUT ) :: a( lda, * ), b( ldb, * ), q( ldq, * ),
245 $ z( ldz, * ), alphar( * ), alphai( * ), beta( * )
246 INTEGER,
INTENT( OUT ) :: ns, nd, info
247 REAL :: qc( ldqc, * ), zc( ldzc, * ), work( * )
250 REAL :: zero, one, half
251 parameter( zero = 0.0, one = 1.0, half = 0.5 )
255 INTEGER :: jw, kwtop, kwbot, istopm, istartm, k, k2, stgexc_info,
256 $ ifst, ilst, lworkreq, qz_small_info
257 REAL :: s, smlnum, ulp, safmin, safmax, c1, s1, temp
267 jw = min( nw, ihi-ilo+1 )
269 IF ( kwtop .EQ. ilo )
THEN
272 s = a( kwtop, kwtop-1 )
278 CALL stgexc( .true., .true., jw, a, lda, b, ldb, qc, ldqc, zc,
279 $ ldzc, ifst, ilst, work, -1, stgexc_info )
280 lworkreq = int( work( 1 ) )
281 CALL slaqz0(
'S',
'V',
'V', jw, 1, jw, a( kwtop, kwtop ), lda,
282 $ b( kwtop, kwtop ), ldb, alphar, alphai, beta, qc,
283 $ ldqc, zc, ldzc, work, -1, rec+1, qz_small_info )
284 lworkreq = max( lworkreq, int( work( 1 ) )+2*jw**2 )
285 lworkreq = max( lworkreq, n*nw, 2*nw**2+n )
286 IF ( lwork .EQ.-1 )
THEN
290 ELSE IF ( lwork .LT. lworkreq )
THEN
295 CALL xerbla(
'SLAQZ3', -info )
300 safmin =
slamch(
'SAFE MINIMUM' )
302 ulp =
slamch(
'PRECISION' )
303 smlnum = safmin*( real( n )/ulp )
305 IF ( ihi .EQ. kwtop )
THEN
307 alphar( kwtop ) = a( kwtop, kwtop )
308 alphai( kwtop ) = zero
309 beta( kwtop ) = b( kwtop, kwtop )
312 IF ( abs( s ) .LE. max( smlnum, ulp*abs( a( kwtop,
316 IF ( kwtop .GT. ilo )
THEN
317 a( kwtop, kwtop-1 ) = zero
324 CALL slacpy(
'ALL', jw, jw, a( kwtop, kwtop ), lda, work, jw )
325 CALL slacpy(
'ALL', jw, jw, b( kwtop, kwtop ), ldb,
330 CALL slaset(
'FULL', jw, jw, zero, one, qc, ldqc )
331 CALL slaset(
'FULL', jw, jw, zero, one, zc, ldzc )
332 CALL slaqz0(
'S',
'V',
'V', jw, 1, jw, a( kwtop, kwtop ), lda,
333 $ b( kwtop, kwtop ), ldb, alphar, alphai, beta, qc,
334 $ ldqc, zc, ldzc, work( 2*jw**2+1 ), lwork-2*jw**2,
335 $ rec+1, qz_small_info )
337 IF( qz_small_info .NE. 0 )
THEN
340 ns = jw-qz_small_info
341 CALL slacpy(
'ALL', jw, jw, work, jw, a( kwtop, kwtop ),
343 CALL slacpy(
'ALL', jw, jw, work( jw**2+1 ), jw, b( kwtop,
349 IF ( kwtop .EQ. ilo .OR. s .EQ. zero )
THEN
355 DO WHILE ( k .LE. jw )
357 IF ( kwbot-kwtop+1 .GE. 2 )
THEN
358 bulge = a( kwbot, kwbot-1 ) .NE. zero
363 temp = abs( a( kwbot, kwbot ) )+sqrt( abs( a( kwbot,
364 $ kwbot-1 ) ) )*sqrt( abs( a( kwbot-1, kwbot ) ) )
365 IF( temp .EQ. zero )
THEN
368 IF ( max( abs( s*qc( 1, kwbot-kwtop ) ), abs( s*qc( 1,
369 $ kwbot-kwtop+1 ) ) ) .LE. max( smlnum,
377 CALL stgexc( .true., .true., jw, a( kwtop, kwtop ),
378 $ lda, b( kwtop, kwtop ), ldb, qc, ldqc,
379 $ zc, ldzc, ifst, ilst, work, lwork,
387 temp = abs( a( kwbot, kwbot ) )
388 IF( temp .EQ. zero )
THEN
391 IF ( ( abs( s*qc( 1, kwbot-kwtop+1 ) ) ) .LE. max( ulp*
392 $ temp, smlnum ) )
THEN
399 CALL stgexc( .true., .true., jw, a( kwtop, kwtop ),
400 $ lda, b( kwtop, kwtop ), ldb, qc, ldqc,
401 $ zc, ldzc, ifst, ilst, work, lwork,
416 DO WHILE ( k .LE. ihi )
418 IF ( k .LT. ihi )
THEN
419 IF ( a( k+1, k ) .NE. zero )
THEN
425 CALL slag2( a( k, k ), lda, b( k, k ), ldb, safmin,
426 $ beta( k ), beta( k+1 ), alphar( k ),
427 $ alphar( k+1 ), alphai( k ) )
428 alphai( k+1 ) = -alphai( k )
432 alphar( k ) = a( k, k )
434 beta( k ) = b( k, k )
439 IF ( kwtop .NE. ilo .AND. s .NE. zero )
THEN
441 a( kwtop:kwbot, kwtop-1 ) = a( kwtop, kwtop-1 )*qc( 1,
443 DO k = kwbot-1, kwtop, -1
444 CALL slartg( a( k, kwtop-1 ), a( k+1, kwtop-1 ), c1, s1,
446 a( k, kwtop-1 ) = temp
447 a( k+1, kwtop-1 ) = zero
448 k2 = max( kwtop, k-1 )
449 CALL srot( ihi-k2+1, a( k, k2 ), lda, a( k+1, k2 ), lda,
452 CALL srot( ihi-( k-1 )+1, b( k, k-1 ), ldb, b( k+1,
455 CALL srot( jw, qc( 1, k-kwtop+1 ), 1, qc( 1,
464 DO WHILE ( k .GE. kwtop )
465 IF ( ( k .GE. kwtop+1 ) .AND. a( k+1, k-1 ) .NE. zero )
THEN
469 CALL slaqz2( .true., .true., k2, kwtop, kwtop+jw-1,
470 $ kwbot, a, lda, b, ldb, jw, kwtop, qc,
471 $ ldqc, jw, kwtop, zc, ldzc )
481 CALL slartg( b( k2+1, k2+1 ), b( k2+1, k2 ), c1,
484 b( k2+1, k2+1 ) = temp
486 CALL srot( k2+2-istartm+1, a( istartm, k2+1 ), 1,
487 $ a( istartm, k2 ), 1, c1, s1 )
488 CALL srot( k2-istartm+1, b( istartm, k2+1 ), 1,
489 $ b( istartm, k2 ), 1, c1, s1 )
490 CALL srot( jw, zc( 1, k2+1-kwtop+1 ), 1, zc( 1,
491 $ k2-kwtop+1 ), 1, c1, s1 )
493 CALL slartg( a( k2+1, k2 ), a( k2+2, k2 ), c1, s1,
497 CALL srot( istopm-k2, a( k2+1, k2+1 ), lda,
499 $ k2+1 ), lda, c1, s1 )
500 CALL srot( istopm-k2, b( k2+1, k2+1 ), ldb,
502 $ k2+1 ), ldb, c1, s1 )
503 CALL srot( jw, qc( 1, k2+1-kwtop+1 ), 1, qc( 1,
504 $ k2+2-kwtop+1 ), 1, c1, s1 )
509 CALL slartg( b( kwbot, kwbot ), b( kwbot, kwbot-1 ),
512 b( kwbot, kwbot ) = temp
513 b( kwbot, kwbot-1 ) = zero
514 CALL srot( kwbot-istartm, b( istartm, kwbot ), 1,
515 $ b( istartm, kwbot-1 ), 1, c1, s1 )
516 CALL srot( kwbot-istartm+1, a( istartm, kwbot ), 1,
517 $ a( istartm, kwbot-1 ), 1, c1, s1 )
518 CALL srot( jw, zc( 1, kwbot-kwtop+1 ), 1, zc( 1,
519 $ kwbot-1-kwtop+1 ), 1, c1, s1 )
536 IF ( istopm-ihi > 0 )
THEN
537 CALL sgemm(
'T',
'N', jw, istopm-ihi, jw, one, qc, ldqc,
538 $ a( kwtop, ihi+1 ), lda, zero, work, jw )
539 CALL slacpy(
'ALL', jw, istopm-ihi, work, jw, a( kwtop,
541 CALL sgemm(
'T',
'N', jw, istopm-ihi, jw, one, qc, ldqc,
542 $ b( kwtop, ihi+1 ), ldb, zero, work, jw )
543 CALL slacpy(
'ALL', jw, istopm-ihi, work, jw, b( kwtop,
547 CALL sgemm(
'N',
'N', n, jw, jw, one, q( 1, kwtop ), ldq,
549 $ ldqc, zero, work, n )
550 CALL slacpy(
'ALL', n, jw, work, n, q( 1, kwtop ), ldq )
553 IF ( kwtop-1-istartm+1 > 0 )
THEN
554 CALL sgemm(
'N',
'N', kwtop-istartm, jw, jw, one,
556 $ kwtop ), lda, zc, ldzc, zero, work,
558 CALL slacpy(
'ALL', kwtop-istartm, jw, work, kwtop-istartm,
559 $ a( istartm, kwtop ), lda )
560 CALL sgemm(
'N',
'N', kwtop-istartm, jw, jw, one,
562 $ kwtop ), ldb, zc, ldzc, zero, work,
564 CALL slacpy(
'ALL', kwtop-istartm, jw, work, kwtop-istartm,
565 $ b( istartm, kwtop ), ldb )
568 CALL sgemm(
'N',
'N', n, jw, jw, one, z( 1, kwtop ), ldz,
570 $ ldzc, zero, work, n )
571 CALL slacpy(
'ALL', n, jw, work, n, z( 1, kwtop ), ldz )