235 RECURSIVE SUBROUTINE dlaqz3( ILSCHUR, ILQ, ILZ, N, ILO, IHI, NW,
236 $ A, LDA, B, LDB, Q, LDQ, Z, LDZ, NS,
237 $ ND, ALPHAR, ALPHAI, BETA, QC, LDQC,
238 $ ZC, LDZC, WORK, LWORK, REC, INFO )
242 LOGICAL,
INTENT( IN ) :: ilschur, ilq, ilz
243 INTEGER,
INTENT( IN ) :: n, ilo, ihi, nw, lda, ldb, ldq, ldz,
244 $ ldqc, ldzc, lwork, rec
246 DOUBLE PRECISION,
INTENT( INOUT ) :: a( lda, * ), b( ldb, * ),
247 $ q( ldq, * ), z( ldz, * ), alphar( * ),
248 $ alphai( * ), beta( * )
249 INTEGER,
INTENT( OUT ) :: ns, nd, info
250 DOUBLE PRECISION :: qc( ldqc, * ), zc( ldzc, * ), work( * )
253 DOUBLE PRECISION :: zero, one, half
254 PARAMETER( zero = 0.0d0, one = 1.0d0, half = 0.5d0 )
258 INTEGER :: jw, kwtop, kwbot, istopm, istartm, k, k2, dtgexc_info,
259 $ ifst, ilst, lworkreq, qz_small_info
260 DOUBLE PRECISION :: s, smlnum, ulp, safmin, safmax, c1, s1, temp
265 DOUBLE PRECISION,
EXTERNAL ::
dlamch
270 jw = min( nw, ihi-ilo+1 )
272 IF ( kwtop .EQ. ilo )
THEN
275 s = a( kwtop, kwtop-1 )
281 CALL dtgexc( .true., .true., jw, a, lda, b, ldb, qc, ldqc, zc,
282 $ ldzc, ifst, ilst, work, -1, dtgexc_info )
283 lworkreq = int( work( 1 ) )
284 CALL dlaqz0(
'S',
'V',
'V', jw, 1, jw, a( kwtop, kwtop ), lda,
285 $ b( kwtop, kwtop ), ldb, alphar, alphai, beta, qc,
286 $ ldqc, zc, ldzc, work, -1, rec+1, qz_small_info )
287 lworkreq = max( lworkreq, int( work( 1 ) )+2*jw**2 )
288 lworkreq = max( lworkreq, n*nw, 2*nw**2+n )
289 IF ( lwork .EQ.-1 )
THEN
293 ELSE IF ( lwork .LT. lworkreq )
THEN
298 CALL xerbla(
'DLAQZ3', -info )
303 safmin =
dlamch(
'SAFE MINIMUM' )
305 ulp =
dlamch(
'PRECISION' )
306 smlnum = safmin*( dble( n )/ulp )
308 IF ( ihi .EQ. kwtop )
THEN
310 alphar( kwtop ) = a( kwtop, kwtop )
311 alphai( kwtop ) = zero
312 beta( kwtop ) = b( kwtop, kwtop )
315 IF ( abs( s ) .LE. max( smlnum, ulp*abs( a( kwtop,
319 IF ( kwtop .GT. ilo )
THEN
320 a( kwtop, kwtop-1 ) = zero
327 CALL dlacpy(
'ALL', jw, jw, a( kwtop, kwtop ), lda, work, jw )
328 CALL dlacpy(
'ALL', jw, jw, b( kwtop, kwtop ), ldb, work( jw**2+
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 ), lda )
344 CALL dlacpy(
'ALL', jw, jw, work( jw**2+1 ), jw, b( kwtop,
350 IF ( kwtop .EQ. ilo .OR. s .EQ. zero )
THEN
356 DO WHILE ( k .LE. jw )
358 IF ( kwbot-kwtop+1 .GE. 2 )
THEN
359 bulge = a( kwbot, kwbot-1 ) .NE. zero
364 temp = abs( a( kwbot, kwbot ) )+sqrt( abs( a( kwbot,
365 $ kwbot-1 ) ) )*sqrt( abs( a( kwbot-1, kwbot ) ) )
366 IF( temp .EQ. zero )
THEN
369 IF ( max( abs( s*qc( 1, kwbot-kwtop ) ), abs( s*qc( 1,
370 $ kwbot-kwtop+1 ) ) ) .LE. max( smlnum,
378 CALL dtgexc( .true., .true., jw, a( kwtop, kwtop ),
379 $ lda, b( kwtop, kwtop ), ldb, qc, ldqc,
380 $ zc, ldzc, ifst, ilst, work, lwork,
388 temp = abs( a( kwbot, kwbot ) )
389 IF( temp .EQ. zero )
THEN
392 IF ( ( abs( s*qc( 1, kwbot-kwtop+1 ) ) ) .LE. max( ulp*
393 $ temp, smlnum ) )
THEN
400 CALL dtgexc( .true., .true., jw, a( kwtop, kwtop ),
401 $ lda, b( kwtop, kwtop ), ldb, qc, ldqc,
402 $ zc, ldzc, ifst, ilst, work, lwork,
417 DO WHILE ( k .LE. ihi )
419 IF ( k .LT. ihi )
THEN
420 IF ( a( k+1, k ) .NE. zero )
THEN
426 CALL dlag2( a( k, k ), lda, b( k, k ), ldb, safmin,
427 $ beta( k ), beta( k+1 ), alphar( k ),
428 $ alphar( k+1 ), alphai( k ) )
429 alphai( k+1 ) = -alphai( k )
433 alphar( k ) = a( k, k )
435 beta( k ) = b( k, k )
440 IF ( kwtop .NE. ilo .AND. s .NE. zero )
THEN
442 a( kwtop:kwbot, kwtop-1 ) = a( kwtop, kwtop-1 )*qc( 1,
444 DO k = kwbot-1, kwtop, -1
445 CALL dlartg( a( k, kwtop-1 ), a( k+1, kwtop-1 ), c1, s1,
447 a( k, kwtop-1 ) = temp
448 a( k+1, kwtop-1 ) = zero
449 k2 = max( kwtop, k-1 )
450 CALL drot( ihi-k2+1, a( k, k2 ), lda, a( k+1, k2 ), lda, c1,
452 CALL drot( ihi-( k-1 )+1, b( k, k-1 ), ldb, b( k+1, k-1 ),
454 CALL drot( jw, qc( 1, k-kwtop+1 ), 1, qc( 1, k+1-kwtop+1 ),
462 DO WHILE ( k .GE. kwtop )
463 IF ( ( k .GE. kwtop+1 ) .AND. a( k+1, k-1 ) .NE. zero )
THEN
467 CALL dlaqz2( .true., .true., k2, kwtop, kwtop+jw-1,
468 $ kwbot, a, lda, b, ldb, jw, kwtop, qc,
469 $ ldqc, jw, kwtop, zc, ldzc )
479 CALL dlartg( b( k2+1, k2+1 ), b( k2+1, k2 ), c1, s1,
481 b( k2+1, k2+1 ) = temp
483 CALL drot( k2+2-istartm+1, a( istartm, k2+1 ), 1,
484 $ a( istartm, k2 ), 1, c1, s1 )
485 CALL drot( k2-istartm+1, b( istartm, k2+1 ), 1,
486 $ b( istartm, k2 ), 1, c1, s1 )
487 CALL drot( jw, zc( 1, k2+1-kwtop+1 ), 1, zc( 1,
488 $ k2-kwtop+1 ), 1, c1, s1 )
490 CALL dlartg( a( k2+1, k2 ), a( k2+2, k2 ), c1, s1,
494 CALL drot( istopm-k2, a( k2+1, k2+1 ), lda, a( k2+2,
495 $ k2+1 ), lda, c1, s1 )
496 CALL drot( istopm-k2, b( k2+1, k2+1 ), ldb, b( k2+2,
497 $ k2+1 ), ldb, c1, s1 )
498 CALL drot( jw, qc( 1, k2+1-kwtop+1 ), 1, qc( 1,
499 $ k2+2-kwtop+1 ), 1, c1, s1 )
504 CALL dlartg( b( kwbot, kwbot ), b( kwbot, kwbot-1 ), c1,
506 b( kwbot, kwbot ) = temp
507 b( kwbot, kwbot-1 ) = zero
508 CALL drot( kwbot-istartm, b( istartm, kwbot ), 1,
509 $ b( istartm, kwbot-1 ), 1, c1, s1 )
510 CALL drot( kwbot-istartm+1, a( istartm, kwbot ), 1,
511 $ a( istartm, kwbot-1 ), 1, c1, s1 )
512 CALL drot( jw, zc( 1, kwbot-kwtop+1 ), 1, zc( 1,
513 $ kwbot-1-kwtop+1 ), 1, c1, s1 )
530 IF ( istopm-ihi > 0 )
THEN
531 CALL dgemm(
'T',
'N', jw, istopm-ihi, jw, one, qc, ldqc,
532 $ a( kwtop, ihi+1 ), lda, zero, work, jw )
533 CALL dlacpy(
'ALL', jw, istopm-ihi, work, jw, a( kwtop,
535 CALL dgemm(
'T',
'N', jw, istopm-ihi, jw, one, qc, ldqc,
536 $ b( kwtop, ihi+1 ), ldb, zero, work, jw )
537 CALL dlacpy(
'ALL', jw, istopm-ihi, work, jw, b( kwtop,
541 CALL dgemm(
'N',
'N', n, jw, jw, one, q( 1, kwtop ), ldq, qc,
542 $ ldqc, zero, work, n )
543 CALL dlacpy(
'ALL', n, jw, work, n, q( 1, kwtop ), ldq )
546 IF ( kwtop-1-istartm+1 > 0 )
THEN
547 CALL dgemm(
'N',
'N', kwtop-istartm, jw, jw, one, a( istartm,
548 $ kwtop ), lda, zc, ldzc, zero, work,
550 CALL dlacpy(
'ALL', kwtop-istartm, jw, work, kwtop-istartm,
551 $ a( istartm, kwtop ), lda )
552 CALL dgemm(
'N',
'N', kwtop-istartm, jw, jw, one, b( istartm,
553 $ kwtop ), ldb, zc, ldzc, zero, work,
555 CALL dlacpy(
'ALL', kwtop-istartm, jw, work, kwtop-istartm,
556 $ b( istartm, kwtop ), ldb )
559 CALL dgemm(
'N',
'N', n, jw, jw, one, z( 1, kwtop ), ldz, zc,
560 $ ldzc, zero, work, n )
561 CALL dlacpy(
'ALL', n, jw, work, n, z( 1, kwtop ), ldz )
subroutine xerbla(srname, info)
subroutine dgemm(transa, transb, m, n, k, alpha, a, lda, b, ldb, beta, c, ldc)
DGEMM
subroutine dlacpy(uplo, m, n, a, lda, b, ldb)
DLACPY copies all or part of one two-dimensional array to another.
subroutine dlag2(a, lda, b, ldb, safmin, scale1, scale2, wr1, wr2, wi)
DLAG2 computes the eigenvalues of a 2-by-2 generalized eigenvalue problem, with scaling as necessary ...
double precision function dlamch(cmach)
DLAMCH
recursive subroutine dlaqz0(wants, wantq, wantz, n, ilo, ihi, a, lda, b, ldb, alphar, alphai, beta, q, ldq, z, ldz, work, lwork, rec, info)
DLAQZ0
subroutine dlaqz2(ilq, ilz, k, istartm, istopm, ihi, a, lda, b, ldb, nq, qstart, q, ldq, nz, zstart, z, ldz)
DLAQZ2
recursive subroutine dlaqz3(ilschur, ilq, ilz, n, ilo, ihi, nw, a, lda, b, ldb, q, ldq, z, ldz, ns, nd, alphar, alphai, beta, qc, ldqc, zc, ldzc, work, lwork, rec, info)
DLAQZ3
subroutine dlartg(f, g, c, s, r)
DLARTG generates a plane rotation with real cosine and real sine.
subroutine dlaset(uplo, m, n, alpha, beta, a, lda)
DLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values.
subroutine drot(n, dx, incx, dy, incy, c, s)
DROT
subroutine dtgexc(wantq, wantz, n, a, lda, b, ldb, q, ldq, z, ldz, ifst, ilst, work, lwork, info)
DTGEXC