233      RECURSIVE SUBROUTINE dlaqz3( ILSCHUR, ILQ, ILZ, N, ILO, IHI,
 
  235     $                             A, LDA, B, LDB, Q, LDQ, Z, LDZ, NS,
 
  236     $                             ND, ALPHAR, ALPHAI, BETA, QC, LDQC,
 
  237     $                             ZC, LDZC, WORK, LWORK, REC, INFO )
 
  241      LOGICAL, 
INTENT( IN ) :: ilschur, ilq, ilz
 
  242      INTEGER, 
INTENT( IN ) :: n, ilo, ihi, nw, lda, ldb, ldq, ldz,
 
  243     $         ldqc, ldzc, lwork, rec
 
  245      DOUBLE PRECISION, 
INTENT( INOUT ) :: a( lda, * ), b( ldb, * ),
 
  246     $                  q( ldq, * ), z( ldz, * ), alphar( * ),
 
  247     $                  alphai( * ), beta( * )
 
  248      INTEGER, 
INTENT( OUT ) :: ns, nd, info
 
  249      DOUBLE PRECISION :: qc( ldqc, * ), zc( ldzc, * ), work( * )
 
  252      DOUBLE PRECISION :: zero, one, half
 
  253      parameter( zero = 0.0d0, one = 1.0d0, half = 0.5d0 )
 
  257      INTEGER :: jw, kwtop, kwbot, istopm, istartm, k, k2, dtgexc_info,
 
  258     $           ifst, ilst, lworkreq, qz_small_info
 
  259      DOUBLE PRECISION :: s, smlnum, ulp, safmin, safmax, c1, s1, temp
 
  264      DOUBLE PRECISION, 
EXTERNAL :: 
dlamch 
  269      jw = min( nw, ihi-ilo+1 )
 
  271      IF ( kwtop .EQ. ilo ) 
THEN 
  274         s = a( kwtop, kwtop-1 )
 
  280      CALL dtgexc( .true., .true., jw, a, lda, b, ldb, qc, ldqc, zc,
 
  281     $             ldzc, ifst, ilst, work, -1, dtgexc_info )
 
  282      lworkreq = int( work( 1 ) )
 
  283      CALL dlaqz0( 
'S', 
'V', 
'V', jw, 1, jw, a( kwtop, kwtop ), lda,
 
  284     $             b( kwtop, kwtop ), ldb, alphar, alphai, beta, qc,
 
  285     $             ldqc, zc, ldzc, work, -1, rec+1, qz_small_info )
 
  286      lworkreq = max( lworkreq, int( work( 1 ) )+2*jw**2 )
 
  287      lworkreq = max( lworkreq, n*nw, 2*nw**2+n )
 
  288      IF ( lwork .EQ.-1 ) 
THEN 
  292      ELSE IF ( lwork .LT. lworkreq ) 
THEN 
  297         CALL xerbla( 
'DLAQZ3', -info )
 
  302      safmin = 
dlamch( 
'SAFE MINIMUM' )
 
  304      ulp = 
dlamch( 
'PRECISION' )
 
  305      smlnum = safmin*( dble( n )/ulp )
 
  307      IF ( ihi .EQ. kwtop ) 
THEN 
  309         alphar( kwtop ) = a( kwtop, kwtop )
 
  310         alphai( kwtop ) = zero
 
  311         beta( kwtop ) = b( kwtop, kwtop )
 
  314         IF ( abs( s ) .LE. max( smlnum, ulp*abs( a( kwtop,
 
  318            IF ( kwtop .GT. ilo ) 
THEN 
  319               a( kwtop, kwtop-1 ) = zero
 
  326      CALL dlacpy( 
'ALL', jw, jw, a( kwtop, kwtop ), lda, work, jw )
 
  327      CALL dlacpy( 
'ALL', jw, jw, b( kwtop, kwtop ), ldb,
 
  332      CALL dlaset( 
'FULL', jw, jw, zero, one, qc, ldqc )
 
  333      CALL dlaset( 
'FULL', jw, jw, zero, one, zc, ldzc )
 
  334      CALL dlaqz0( 
'S', 
'V', 
'V', jw, 1, jw, a( kwtop, kwtop ), lda,
 
  335     $             b( kwtop, kwtop ), ldb, alphar, alphai, beta, qc,
 
  336     $             ldqc, zc, ldzc, work( 2*jw**2+1 ), lwork-2*jw**2,
 
  337     $             rec+1, qz_small_info )
 
  339      IF( qz_small_info .NE. 0 ) 
THEN 
  342         ns = jw-qz_small_info
 
  343         CALL dlacpy( 
'ALL', jw, jw, work, jw, a( kwtop, kwtop ),
 
  345         CALL dlacpy( 
'ALL', jw, jw, work( jw**2+1 ), jw, b( kwtop,
 
  351      IF ( kwtop .EQ. ilo .OR. s .EQ. zero ) 
THEN 
  357         DO WHILE ( k .LE. jw )
 
  359            IF ( kwbot-kwtop+1 .GE. 2 ) 
THEN 
  360               bulge = a( kwbot, kwbot-1 ) .NE. zero
 
  365               temp = abs( a( kwbot, kwbot ) )+sqrt( abs( a( kwbot,
 
  366     $            kwbot-1 ) ) )*sqrt( abs( a( kwbot-1, kwbot ) ) )
 
  367               IF( temp .EQ. zero )
THEN 
  370               IF ( max( abs( s*qc( 1, kwbot-kwtop ) ), abs( s*qc( 1,
 
  371     $            kwbot-kwtop+1 ) ) ) .LE. max( smlnum,
 
  379                  CALL dtgexc( .true., .true., jw, a( kwtop, kwtop ),
 
  380     $                         lda, b( kwtop, kwtop ), ldb, qc, ldqc,
 
  381     $                         zc, ldzc, ifst, ilst, work, lwork,
 
  389               temp = abs( a( kwbot, kwbot ) )
 
  390               IF( temp .EQ. zero ) 
THEN 
  393               IF ( ( abs( s*qc( 1, kwbot-kwtop+1 ) ) ) .LE. max( ulp*
 
  394     $            temp, smlnum ) ) 
THEN 
  401                  CALL dtgexc( .true., .true., jw, a( kwtop, kwtop ),
 
  402     $                         lda, b( kwtop, kwtop ), ldb, qc, ldqc,
 
  403     $                         zc, ldzc, ifst, ilst, work, lwork,
 
  418      DO WHILE ( k .LE. ihi )
 
  420         IF ( k .LT. ihi ) 
THEN 
  421            IF ( a( k+1, k ) .NE. zero ) 
THEN 
  427            CALL dlag2( a( k, k ), lda, b( k, k ), ldb, safmin,
 
  428     $                  beta( k ), beta( k+1 ), alphar( k ),
 
  429     $                  alphar( k+1 ), alphai( k ) )
 
  430            alphai( k+1 ) = -alphai( k )
 
  434            alphar( k ) = a( k, k )
 
  436            beta( k ) = b( k, k )
 
  441      IF ( kwtop .NE. ilo .AND. s .NE. zero ) 
THEN 
  443         a( kwtop:kwbot, kwtop-1 ) = a( kwtop, kwtop-1 )*qc( 1,
 
  445         DO k = kwbot-1, kwtop, -1
 
  446            CALL dlartg( a( k, kwtop-1 ), a( k+1, kwtop-1 ), c1, s1,
 
  448            a( k, kwtop-1 ) = temp
 
  449            a( k+1, kwtop-1 ) = zero
 
  450            k2 = max( kwtop, k-1 )
 
  451            CALL drot( ihi-k2+1, a( k, k2 ), lda, a( k+1, k2 ), lda,
 
  454            CALL drot( ihi-( k-1 )+1, b( k, k-1 ), ldb, b( k+1,
 
  457            CALL drot( jw, qc( 1, k-kwtop+1 ), 1, qc( 1,
 
  466         DO WHILE ( k .GE. kwtop )
 
  467            IF ( ( k .GE. kwtop+1 ) .AND. a( k+1, k-1 ) .NE. zero ) 
THEN 
  471                  CALL dlaqz2( .true., .true., k2, kwtop, kwtop+jw-1,
 
  472     $                         kwbot, a, lda, b, ldb, jw, kwtop, qc,
 
  473     $                         ldqc, jw, kwtop, zc, ldzc )
 
  483                  CALL dlartg( b( k2+1, k2+1 ), b( k2+1, k2 ), c1,
 
  486                  b( k2+1, k2+1 ) = temp
 
  488                  CALL drot( k2+2-istartm+1, a( istartm, k2+1 ), 1,
 
  489     $                       a( istartm, k2 ), 1, c1, s1 )
 
  490                  CALL drot( k2-istartm+1, b( istartm, k2+1 ), 1,
 
  491     $                       b( istartm, k2 ), 1, c1, s1 )
 
  492                  CALL drot( jw, zc( 1, k2+1-kwtop+1 ), 1, zc( 1,
 
  493     $                       k2-kwtop+1 ), 1, c1, s1 )
 
  495                  CALL dlartg( a( k2+1, k2 ), a( k2+2, k2 ), c1, s1,
 
  499                  CALL drot( istopm-k2, a( k2+1, k2+1 ), lda,
 
  501     $                       k2+1 ), lda, c1, s1 )
 
  502                  CALL drot( istopm-k2, b( k2+1, k2+1 ), ldb,
 
  504     $                       k2+1 ), ldb, c1, s1 )
 
  505                  CALL drot( jw, qc( 1, k2+1-kwtop+1 ), 1, qc( 1,
 
  506     $                       k2+2-kwtop+1 ), 1, c1, s1 )
 
  511               CALL dlartg( b( kwbot, kwbot ), b( kwbot, kwbot-1 ),
 
  514               b( kwbot, kwbot ) = temp
 
  515               b( kwbot, kwbot-1 ) = zero
 
  516               CALL drot( kwbot-istartm, b( istartm, kwbot ), 1,
 
  517     $                    b( istartm, kwbot-1 ), 1, c1, s1 )
 
  518               CALL drot( kwbot-istartm+1, a( istartm, kwbot ), 1,
 
  519     $                    a( istartm, kwbot-1 ), 1, c1, s1 )
 
  520               CALL drot( jw, zc( 1, kwbot-kwtop+1 ), 1, zc( 1,
 
  521     $                    kwbot-1-kwtop+1 ), 1, c1, s1 )
 
  538      IF ( istopm-ihi > 0 ) 
THEN 
  539         CALL dgemm( 
'T', 
'N', jw, istopm-ihi, jw, one, qc, ldqc,
 
  540     $               a( kwtop, ihi+1 ), lda, zero, work, jw )
 
  541         CALL dlacpy( 
'ALL', jw, istopm-ihi, work, jw, a( kwtop,
 
  543         CALL dgemm( 
'T', 
'N', jw, istopm-ihi, jw, one, qc, ldqc,
 
  544     $               b( kwtop, ihi+1 ), ldb, zero, work, jw )
 
  545         CALL dlacpy( 
'ALL', jw, istopm-ihi, work, jw, b( kwtop,
 
  549         CALL dgemm( 
'N', 
'N', n, jw, jw, one, q( 1, kwtop ), ldq,
 
  551     $               ldqc, zero, work, n )
 
  552         CALL dlacpy( 
'ALL', n, jw, work, n, q( 1, kwtop ), ldq )
 
  555      IF ( kwtop-1-istartm+1 > 0 ) 
THEN 
  556         CALL dgemm( 
'N', 
'N', kwtop-istartm, jw, jw, one,
 
  558     $               kwtop ), lda, zc, ldzc, zero, work,
 
  560         CALL dlacpy( 
'ALL', kwtop-istartm, jw, work, kwtop-istartm,
 
  561     $                a( istartm, kwtop ), lda )
 
  562         CALL dgemm( 
'N', 
'N', kwtop-istartm, jw, jw, one,
 
  564     $               kwtop ), ldb, zc, ldzc, zero, work,
 
  566         CALL dlacpy( 
'ALL', kwtop-istartm, jw, work, kwtop-istartm,
 
  567     $                b( istartm, kwtop ), ldb )
 
  570         CALL dgemm( 
'N', 
'N', n, jw, jw, one, z( 1, kwtop ), ldz,
 
  572     $               ldzc, zero, work, n )
 
  573         CALL dlacpy( 
'ALL', n, jw, work, n, z( 1, kwtop ), ldz )