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 CALL slabad( safmin, safmax )
304 ulp =
slamch(
'PRECISION' )
305 smlnum = safmin*( real( 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 slacpy(
'ALL', jw, jw, a( kwtop, kwtop ), lda, work, jw )
327 CALL slacpy(
'ALL', jw, jw, b( kwtop, kwtop ), ldb, work( jw**2+
331 CALL slaset(
'FULL', jw, jw, zero, one, qc, ldqc )
332 CALL slaset(
'FULL', jw, jw, zero, one, zc, ldzc )
333 CALL slaqz0(
'S',
'V',
'V', jw, 1, jw, a( kwtop, kwtop ), lda,
334 $ b( kwtop, kwtop ), ldb, alphar, alphai, beta, qc,
335 $ ldqc, zc, ldzc, work( 2*jw**2+1 ), lwork-2*jw**2,
336 $ rec+1, qz_small_info )
338 IF( qz_small_info .NE. 0 )
THEN
341 ns = jw-qz_small_info
342 CALL slacpy(
'ALL', jw, jw, work, jw, a( kwtop, kwtop ), lda )
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, c1,
451 CALL srot( ihi-( k-1 )+1, b( k, k-1 ), ldb, b( k+1, k-1 ),
453 CALL srot( jw, qc( 1, k-kwtop+1 ), 1, qc( 1, k+1-kwtop+1 ),
461 DO WHILE ( k .GE. kwtop )
462 IF ( ( k .GE. kwtop+1 ) .AND. a( k+1, k-1 ) .NE. zero )
THEN
466 CALL slaqz2( .true., .true., k2, kwtop, kwtop+jw-1,
467 $ kwbot, a, lda, b, ldb, jw, kwtop, qc,
468 $ ldqc, jw, kwtop, zc, ldzc )
478 CALL slartg( b( k2+1, k2+1 ), b( k2+1, k2 ), c1, s1,
480 b( k2+1, k2+1 ) = temp
482 CALL srot( k2+2-istartm+1, a( istartm, k2+1 ), 1,
483 $ a( istartm, k2 ), 1, c1, s1 )
484 CALL srot( k2-istartm+1, b( istartm, k2+1 ), 1,
485 $ b( istartm, k2 ), 1, c1, s1 )
486 CALL srot( jw, zc( 1, k2+1-kwtop+1 ), 1, zc( 1,
487 $ k2-kwtop+1 ), 1, c1, s1 )
489 CALL slartg( a( k2+1, k2 ), a( k2+2, k2 ), c1, s1,
493 CALL srot( istopm-k2, a( k2+1, k2+1 ), lda, a( k2+2,
494 $ k2+1 ), lda, c1, s1 )
495 CALL srot( istopm-k2, b( k2+1, k2+1 ), ldb, b( k2+2,
496 $ k2+1 ), ldb, c1, s1 )
497 CALL srot( jw, qc( 1, k2+1-kwtop+1 ), 1, qc( 1,
498 $ k2+2-kwtop+1 ), 1, c1, s1 )
503 CALL slartg( b( kwbot, kwbot ), b( kwbot, kwbot-1 ), c1,
505 b( kwbot, kwbot ) = temp
506 b( kwbot, kwbot-1 ) = zero
507 CALL srot( kwbot-istartm, b( istartm, kwbot ), 1,
508 $ b( istartm, kwbot-1 ), 1, c1, s1 )
509 CALL srot( kwbot-istartm+1, a( istartm, kwbot ), 1,
510 $ a( istartm, kwbot-1 ), 1, c1, s1 )
511 CALL srot( jw, zc( 1, kwbot-kwtop+1 ), 1, zc( 1,
512 $ kwbot-1-kwtop+1 ), 1, c1, s1 )
529 IF ( istopm-ihi > 0 )
THEN
530 CALL sgemm(
'T',
'N', jw, istopm-ihi, jw, one, qc, ldqc,
531 $ a( kwtop, ihi+1 ), lda, zero, work, jw )
532 CALL slacpy(
'ALL', jw, istopm-ihi, work, jw, a( kwtop,
534 CALL sgemm(
'T',
'N', jw, istopm-ihi, jw, one, qc, ldqc,
535 $ b( kwtop, ihi+1 ), ldb, zero, work, jw )
536 CALL slacpy(
'ALL', jw, istopm-ihi, work, jw, b( kwtop,
540 CALL sgemm(
'N',
'N', n, jw, jw, one, q( 1, kwtop ), ldq, qc,
541 $ ldqc, zero, work, n )
542 CALL slacpy(
'ALL', n, jw, work, n, q( 1, kwtop ), ldq )
545 IF ( kwtop-1-istartm+1 > 0 )
THEN
546 CALL sgemm(
'N',
'N', kwtop-istartm, jw, jw, one, a( istartm,
547 $ kwtop ), lda, zc, ldzc, zero, work,
549 CALL slacpy(
'ALL', kwtop-istartm, jw, work, kwtop-istartm,
550 $ a( istartm, kwtop ), lda )
551 CALL sgemm(
'N',
'N', kwtop-istartm, jw, jw, one, b( istartm,
552 $ kwtop ), ldb, zc, ldzc, zero, work,
554 CALL slacpy(
'ALL', kwtop-istartm, jw, work, kwtop-istartm,
555 $ b( istartm, kwtop ), ldb )
558 CALL sgemm(
'N',
'N', n, jw, jw, one, z( 1, kwtop ), ldz, zc,
559 $ ldzc, zero, work, n )
560 CALL slacpy(
'ALL', n, jw, work, n, z( 1, kwtop ), ldz )
subroutine slabad(SMALL, LARGE)
SLABAD
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 slacpy(UPLO, M, N, A, LDA, B, LDB)
SLACPY copies all or part of one two-dimensional array to another.
subroutine slartg(f, g, c, s, r)
SLARTG generates a plane rotation with real cosine and real sine.
subroutine xerbla(SRNAME, INFO)
XERBLA
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
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 stgexc(WANTQ, WANTZ, N, A, LDA, B, LDB, Q, LDQ, Z, LDZ, IFST, ILST, WORK, LWORK, INFO)
STGEXC
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 ...
subroutine srot(N, SX, INCX, SY, INCY, C, S)
SROT
subroutine sgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
SGEMM
real function slamch(CMACH)
SLAMCH