234 RECURSIVE SUBROUTINE slaqz3( ILSCHUR, ILQ, ILZ, N, ILO, IHI, NW,
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 REAL,
INTENT( INOUT ) :: a( lda, * ), b( ldb, * ), q( ldq, * ),
246 $ z( ldz, * ), alphar( * ), alphai( * ), beta( * )
247 INTEGER,
INTENT( OUT ) :: ns, nd, info
248 REAL :: qc( ldqc, * ), zc( ldzc, * ), work( * )
251 REAL :: zero, one, half
252 PARAMETER( zero = 0.0, one = 1.0, half = 0.5 )
256 INTEGER :: jw, kwtop, kwbot, istopm, istartm, k, k2, stgexc_info,
257 $ ifst, ilst, lworkreq, qz_small_info
258 REAL :: s, smlnum, ulp, safmin, safmax, c1, s1, temp
268 jw = min( nw, ihi-ilo+1 )
270 IF ( kwtop .EQ. ilo )
THEN
273 s = a( kwtop, kwtop-1 )
279 CALL stgexc( .true., .true., jw, a, lda, b, ldb, qc, ldqc, zc,
280 $ ldzc, ifst, ilst, work, -1, stgexc_info )
281 lworkreq = int( work( 1 ) )
282 CALL slaqz0(
'S',
'V',
'V', jw, 1, jw, a( kwtop, kwtop ), lda,
283 $ b( kwtop, kwtop ), ldb, alphar, alphai, beta, qc,
284 $ ldqc, zc, ldzc, work, -1, rec+1, qz_small_info )
285 lworkreq = max( lworkreq, int( work( 1 ) )+2*jw**2 )
286 lworkreq = max( lworkreq, n*nw, 2*nw**2+n )
287 IF ( lwork .EQ.-1 )
THEN
291 ELSE IF ( lwork .LT. lworkreq )
THEN
296 CALL xerbla(
'SLAQZ3', -info )
301 safmin =
slamch(
'SAFE MINIMUM' )
303 ulp =
slamch(
'PRECISION' )
304 smlnum = safmin*( real( n )/ulp )
306 IF ( ihi .EQ. kwtop )
THEN
308 alphar( kwtop ) = a( kwtop, kwtop )
309 alphai( kwtop ) = zero
310 beta( kwtop ) = b( kwtop, kwtop )
313 IF ( abs( s ) .LE. max( smlnum, ulp*abs( a( kwtop,
317 IF ( kwtop .GT. ilo )
THEN
318 a( kwtop, kwtop-1 ) = zero
325 CALL slacpy(
'ALL', jw, jw, a( kwtop, kwtop ), lda, work, jw )
326 CALL slacpy(
'ALL', jw, jw, b( kwtop, kwtop ), ldb, work( jw**2+
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 ), lda )
342 CALL slacpy(
'ALL', jw, jw, work( jw**2+1 ), jw, b( kwtop,
348 IF ( kwtop .EQ. ilo .OR. s .EQ. zero )
THEN
354 DO WHILE ( k .LE. jw )
356 IF ( kwbot-kwtop+1 .GE. 2 )
THEN
357 bulge = a( kwbot, kwbot-1 ) .NE. zero
362 temp = abs( a( kwbot, kwbot ) )+sqrt( abs( a( kwbot,
363 $ kwbot-1 ) ) )*sqrt( abs( a( kwbot-1, kwbot ) ) )
364 IF( temp .EQ. zero )
THEN
367 IF ( max( abs( s*qc( 1, kwbot-kwtop ) ), abs( s*qc( 1,
368 $ kwbot-kwtop+1 ) ) ) .LE. max( smlnum,
376 CALL stgexc( .true., .true., jw, a( kwtop, kwtop ),
377 $ lda, b( kwtop, kwtop ), ldb, qc, ldqc,
378 $ zc, ldzc, ifst, ilst, work, lwork,
386 temp = abs( a( kwbot, kwbot ) )
387 IF( temp .EQ. zero )
THEN
390 IF ( ( abs( s*qc( 1, kwbot-kwtop+1 ) ) ) .LE. max( ulp*
391 $ temp, smlnum ) )
THEN
398 CALL stgexc( .true., .true., jw, a( kwtop, kwtop ),
399 $ lda, b( kwtop, kwtop ), ldb, qc, ldqc,
400 $ zc, ldzc, ifst, ilst, work, lwork,
415 DO WHILE ( k .LE. ihi )
417 IF ( k .LT. ihi )
THEN
418 IF ( a( k+1, k ) .NE. zero )
THEN
424 CALL slag2( a( k, k ), lda, b( k, k ), ldb, safmin,
425 $ beta( k ), beta( k+1 ), alphar( k ),
426 $ alphar( k+1 ), alphai( k ) )
427 alphai( k+1 ) = -alphai( k )
431 alphar( k ) = a( k, k )
433 beta( k ) = b( k, k )
438 IF ( kwtop .NE. ilo .AND. s .NE. zero )
THEN
440 a( kwtop:kwbot, kwtop-1 ) = a( kwtop, kwtop-1 )*qc( 1,
442 DO k = kwbot-1, kwtop, -1
443 CALL slartg( a( k, kwtop-1 ), a( k+1, kwtop-1 ), c1, s1,
445 a( k, kwtop-1 ) = temp
446 a( k+1, kwtop-1 ) = zero
447 k2 = max( kwtop, k-1 )
448 CALL srot( ihi-k2+1, a( k, k2 ), lda, a( k+1, k2 ), lda, c1,
450 CALL srot( ihi-( k-1 )+1, b( k, k-1 ), ldb, b( k+1, k-1 ),
452 CALL srot( jw, qc( 1, k-kwtop+1 ), 1, qc( 1, k+1-kwtop+1 ),
460 DO WHILE ( k .GE. kwtop )
461 IF ( ( k .GE. kwtop+1 ) .AND. a( k+1, k-1 ) .NE. zero )
THEN
465 CALL slaqz2( .true., .true., k2, kwtop, kwtop+jw-1,
466 $ kwbot, a, lda, b, ldb, jw, kwtop, qc,
467 $ ldqc, jw, kwtop, zc, ldzc )
477 CALL slartg( b( k2+1, k2+1 ), b( k2+1, k2 ), c1, s1,
479 b( k2+1, k2+1 ) = temp
481 CALL srot( k2+2-istartm+1, a( istartm, k2+1 ), 1,
482 $ a( istartm, k2 ), 1, c1, s1 )
483 CALL srot( k2-istartm+1, b( istartm, k2+1 ), 1,
484 $ b( istartm, k2 ), 1, c1, s1 )
485 CALL srot( jw, zc( 1, k2+1-kwtop+1 ), 1, zc( 1,
486 $ k2-kwtop+1 ), 1, c1, s1 )
488 CALL slartg( a( k2+1, k2 ), a( k2+2, k2 ), c1, s1,
492 CALL srot( istopm-k2, a( k2+1, k2+1 ), lda, a( k2+2,
493 $ k2+1 ), lda, c1, s1 )
494 CALL srot( istopm-k2, b( k2+1, k2+1 ), ldb, b( k2+2,
495 $ k2+1 ), ldb, c1, s1 )
496 CALL srot( jw, qc( 1, k2+1-kwtop+1 ), 1, qc( 1,
497 $ k2+2-kwtop+1 ), 1, c1, s1 )
502 CALL slartg( b( kwbot, kwbot ), b( kwbot, kwbot-1 ), c1,
504 b( kwbot, kwbot ) = temp
505 b( kwbot, kwbot-1 ) = zero
506 CALL srot( kwbot-istartm, b( istartm, kwbot ), 1,
507 $ b( istartm, kwbot-1 ), 1, c1, s1 )
508 CALL srot( kwbot-istartm+1, a( istartm, kwbot ), 1,
509 $ a( istartm, kwbot-1 ), 1, c1, s1 )
510 CALL srot( jw, zc( 1, kwbot-kwtop+1 ), 1, zc( 1,
511 $ kwbot-1-kwtop+1 ), 1, c1, s1 )
528 IF ( istopm-ihi > 0 )
THEN
529 CALL sgemm(
'T',
'N', jw, istopm-ihi, jw, one, qc, ldqc,
530 $ a( kwtop, ihi+1 ), lda, zero, work, jw )
531 CALL slacpy(
'ALL', jw, istopm-ihi, work, jw, a( kwtop,
533 CALL sgemm(
'T',
'N', jw, istopm-ihi, jw, one, qc, ldqc,
534 $ b( kwtop, ihi+1 ), ldb, zero, work, jw )
535 CALL slacpy(
'ALL', jw, istopm-ihi, work, jw, b( kwtop,
539 CALL sgemm(
'N',
'N', n, jw, jw, one, q( 1, kwtop ), ldq, qc,
540 $ ldqc, zero, work, n )
541 CALL slacpy(
'ALL', n, jw, work, n, q( 1, kwtop ), ldq )
544 IF ( kwtop-1-istartm+1 > 0 )
THEN
545 CALL sgemm(
'N',
'N', kwtop-istartm, jw, jw, one, a( istartm,
546 $ kwtop ), lda, zc, ldzc, zero, work,
548 CALL slacpy(
'ALL', kwtop-istartm, jw, work, kwtop-istartm,
549 $ a( istartm, kwtop ), lda )
550 CALL sgemm(
'N',
'N', kwtop-istartm, jw, jw, one, b( istartm,
551 $ kwtop ), ldb, zc, ldzc, zero, work,
553 CALL slacpy(
'ALL', kwtop-istartm, jw, work, kwtop-istartm,
554 $ b( istartm, kwtop ), ldb )
557 CALL sgemm(
'N',
'N', n, jw, jw, one, z( 1, kwtop ), ldz, zc,
558 $ ldzc, zero, work, n )
559 CALL slacpy(
'ALL', n, jw, work, n, z( 1, kwtop ), ldz )
subroutine xerbla(srname, info)
subroutine sgemm(transa, transb, m, n, k, alpha, a, lda, b, ldb, beta, c, ldc)
SGEMM
subroutine slacpy(uplo, m, n, a, lda, b, ldb)
SLACPY copies all or part of one two-dimensional array to another.
subroutine slag2(a, lda, b, ldb, safmin, scale1, scale2, wr1, wr2, wi)
SLAG2 computes the eigenvalues of a 2-by-2 generalized eigenvalue problem, with scaling as necessary ...
real function slamch(cmach)
SLAMCH
recursive subroutine slaqz0(wants, wantq, wantz, n, ilo, ihi, a, lda, b, ldb, alphar, alphai, beta, q, ldq, z, ldz, work, lwork, rec, info)
SLAQZ0
subroutine slaqz2(ilq, ilz, k, istartm, istopm, ihi, a, lda, b, ldb, nq, qstart, q, ldq, nz, zstart, z, ldz)
SLAQZ2
recursive subroutine slaqz3(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)
SLAQZ3
subroutine slartg(f, g, c, s, r)
SLARTG generates a plane rotation with real cosine and real sine.
subroutine slaset(uplo, m, n, alpha, beta, a, lda)
SLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values.
subroutine srot(n, sx, incx, sy, incy, c, s)
SROT
real function sroundup_lwork(lwork)
SROUNDUP_LWORK
subroutine stgexc(wantq, wantz, n, a, lda, b, ldb, q, ldq, z, ldz, ifst, ilst, work, lwork, info)
STGEXC