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 )