300 RECURSIVE SUBROUTINE slaqz0( WANTS, WANTQ, WANTZ, N, ILO, IHI, A,
301 $ LDA, B, LDB, ALPHAR, ALPHAI, BETA,
302 $ Q, LDQ, Z, LDZ, WORK, LWORK, REC,
307 CHARACTER,
INTENT( IN ) :: wants, wantq, wantz
308 INTEGER,
INTENT( IN ) :: n, ilo, ihi, lda, ldb, ldq, ldz, lwork,
311 INTEGER,
INTENT( OUT ) :: info
313 REAL,
INTENT( INOUT ) :: a( lda, * ), b( ldb, * ), q( ldq, * ),
314 $ z( ldz, * ), alphar( * ), alphai( * ), beta( * ), work( * )
317 REAL :: zero, one, half
318 PARAMETER( zero = 0.0, one = 1.0, half = 0.5 )
321 REAL :: smlnum, ulp, eshift, safmin, safmax, c1, s1, temp, swap,
323 INTEGER :: istart, istop, iiter, maxit, istart2, k, ld, nshifts,
324 $ nblock, nw, nmin, nibble, n_undeflated, n_deflated,
325 $ ns, sweep_info, shiftpos, lworkreq, k2, istartm,
326 $ istopm, iwants, iwantq, iwantz, norm_info, aed_info,
327 $ nwr, nbr, nsr, itemp1, itemp2, rcost, i
328 LOGICAL :: ilschur, ilq, ilz
329 CHARACTER :: jbcmpz*3
335 LOGICAL,
EXTERNAL ::
lsame
336 INTEGER,
EXTERNAL ::
ilaenv
341 IF(
lsame( wants,
'E' ) )
THEN
344 ELSE IF(
lsame( wants,
'S' ) )
THEN
351 IF(
lsame( wantq,
'N' ) )
THEN
354 ELSE IF(
lsame( wantq,
'V' ) )
THEN
357 ELSE IF(
lsame( wantq,
'I' ) )
THEN
364 IF(
lsame( wantz,
'N' ) )
THEN
367 ELSE IF(
lsame( wantz,
'V' ) )
THEN
370 ELSE IF(
lsame( wantz,
'I' ) )
THEN
380 IF( iwants.EQ.0 )
THEN
382 ELSE IF( iwantq.EQ.0 )
THEN
384 ELSE IF( iwantz.EQ.0 )
THEN
386 ELSE IF( n.LT.0 )
THEN
388 ELSE IF( ilo.LT.1 )
THEN
390 ELSE IF( ihi.GT.n .OR. ihi.LT.ilo-1 )
THEN
392 ELSE IF( lda.LT.n )
THEN
394 ELSE IF( ldb.LT.n )
THEN
396 ELSE IF( ldq.LT.1 .OR. ( ilq .AND. ldq.LT.n ) )
THEN
398 ELSE IF( ldz.LT.1 .OR. ( ilz .AND. ldz.LT.n ) )
THEN
402 CALL xerbla(
'SLAQZ0', -info )
410 work( 1 ) = real( 1 )
417 jbcmpz( 1:1 ) = wants
418 jbcmpz( 2:2 ) = wantq
419 jbcmpz( 3:3 ) = wantz
421 nmin =
ilaenv( 12,
'SLAQZ0', jbcmpz, n, ilo, ihi, lwork )
423 nwr =
ilaenv( 13,
'SLAQZ0', jbcmpz, n, ilo, ihi, lwork )
425 nwr = min( ihi-ilo+1, ( n-1 ) / 3, nwr )
427 nibble =
ilaenv( 14,
'SLAQZ0', jbcmpz, n, ilo, ihi, lwork )
429 nsr =
ilaenv( 15,
'SLAQZ0', jbcmpz, n, ilo, ihi, lwork )
430 nsr = min( nsr, ( n+6 ) / 9, ihi-ilo )
431 nsr = max( 2, nsr-mod( nsr, 2 ) )
433 rcost =
ilaenv( 17,
'SLAQZ0', jbcmpz, n, ilo, ihi, lwork )
434 itemp1 = int( nsr/sqrt( 1+2*nsr/( real( rcost )/100*n ) ) )
435 itemp1 = ( ( itemp1-1 )/4 )*4+4
438 IF( n .LT. nmin .OR. rec .GE. 2 )
THEN
439 CALL shgeqz( wants, wantq, wantz, n, ilo, ihi, a, lda, b, ldb,
440 $ alphar, alphai, beta, q, ldq, z, ldz, work,
450 nw = max( nwr, nmin )
451 CALL slaqz3( ilschur, ilq, ilz, n, ilo, ihi, nw, a, lda, b, ldb,
452 $ q, ldq, z, ldz, n_undeflated, n_deflated, alphar,
453 $ alphai, beta, work, nw, work, nw, work, -1, rec,
455 itemp1 = int( work( 1 ) )
457 CALL slaqz4( ilschur, ilq, ilz, n, ilo, ihi, nsr, nbr, alphar,
458 $ alphai, beta, a, lda, b, ldb, q, ldq, z, ldz, work,
459 $ nbr, work, nbr, work, -1, sweep_info )
460 itemp2 = int( work( 1 ) )
462 lworkreq = max( itemp1+2*nw**2, itemp2+2*nbr**2 )
463 IF ( lwork .EQ.-1 )
THEN
464 work( 1 ) = real( lworkreq )
466 ELSE IF ( lwork .LT. lworkreq )
THEN
470 CALL xerbla(
'SLAQZ0', info )
476 IF( iwantq.EQ.3 )
CALL slaset(
'FULL', n, n, zero, one, q, ldq )
477 IF( iwantz.EQ.3 )
CALL slaset(
'FULL', n, n, zero, one, z, ldz )
480 safmin =
slamch(
'SAFE MINIMUM' )
482 CALL slabad( safmin, safmax )
483 ulp =
slamch(
'PRECISION' )
484 smlnum = safmin*( real( n )/ulp )
486 bnorm =
slanhs(
'F', ihi-ilo+1, b( ilo, ilo ), ldb, work )
487 btol = max( safmin, ulp*bnorm )
491 maxit = 3*( ihi-ilo+1 )
495 IF( iiter .GE. maxit )
THEN
499 IF ( istart+1 .GE. istop )
THEN
505 IF ( abs( a( istop-1, istop-2 ) ) .LE. max( smlnum,
506 $ ulp*( abs( a( istop-1, istop-1 ) )+abs( a( istop-2,
507 $ istop-2 ) ) ) ) )
THEN
508 a( istop-1, istop-2 ) = zero
512 ELSE IF ( abs( a( istop, istop-1 ) ) .LE. max( smlnum,
513 $ ulp*( abs( a( istop, istop ) )+abs( a( istop-1,
514 $ istop-1 ) ) ) ) )
THEN
515 a( istop, istop-1 ) = zero
521 IF ( abs( a( istart+2, istart+1 ) ) .LE. max( smlnum,
522 $ ulp*( abs( a( istart+1, istart+1 ) )+abs( a( istart+2,
523 $ istart+2 ) ) ) ) )
THEN
524 a( istart+2, istart+1 ) = zero
528 ELSE IF ( abs( a( istart+1, istart ) ) .LE. max( smlnum,
529 $ ulp*( abs( a( istart, istart ) )+abs( a( istart+1,
530 $ istart+1 ) ) ) ) )
THEN
531 a( istart+1, istart ) = zero
537 IF ( istart+1 .GE. istop )
THEN
543 DO k = istop, istart+1, -1
544 IF ( abs( a( k, k-1 ) ) .LE. max( smlnum, ulp*( abs( a( k,
545 $ k ) )+abs( a( k-1, k-1 ) ) ) ) )
THEN
564 DO WHILE ( k.GE.istart2 )
566 IF( abs( b( k, k ) ) .LT. btol )
THEN
570 DO k2 = k, istart2+1, -1
571 CALL slartg( b( k2-1, k2 ), b( k2-1, k2-1 ), c1, s1,
574 b( k2-1, k2-1 ) = zero
576 CALL srot( k2-2-istartm+1, b( istartm, k2 ), 1,
577 $ b( istartm, k2-1 ), 1, c1, s1 )
578 CALL srot( min( k2+1, istop )-istartm+1, a( istartm,
579 $ k2 ), 1, a( istartm, k2-1 ), 1, c1, s1 )
581 CALL srot( n, z( 1, k2 ), 1, z( 1, k2-1 ), 1, c1,
585 IF( k2.LT.istop )
THEN
586 CALL slartg( a( k2, k2-1 ), a( k2+1, k2-1 ), c1,
589 a( k2+1, k2-1 ) = zero
591 CALL srot( istopm-k2+1, a( k2, k2 ), lda, a( k2+1,
592 $ k2 ), lda, c1, s1 )
593 CALL srot( istopm-k2+1, b( k2, k2 ), ldb, b( k2+1,
594 $ k2 ), ldb, c1, s1 )
596 CALL srot( n, q( 1, k2 ), 1, q( 1, k2+1 ), 1,
603 IF( istart2.LT.istop )
THEN
604 CALL slartg( a( istart2, istart2 ), a( istart2+1,
605 $ istart2 ), c1, s1, temp )
606 a( istart2, istart2 ) = temp
607 a( istart2+1, istart2 ) = zero
609 CALL srot( istopm-( istart2+1 )+1, a( istart2,
610 $ istart2+1 ), lda, a( istart2+1,
611 $ istart2+1 ), lda, c1, s1 )
612 CALL srot( istopm-( istart2+1 )+1, b( istart2,
613 $ istart2+1 ), ldb, b( istart2+1,
614 $ istart2+1 ), ldb, c1, s1 )
616 CALL srot( n, q( 1, istart2 ), 1, q( 1,
617 $ istart2+1 ), 1, c1, s1 )
629 IF ( istart2 .GE. istop )
THEN
640 IF ( istop-istart2+1 .LT. nmin )
THEN
644 IF ( istop-istart+1 .LT. nmin )
THEN
655 CALL slaqz3( ilschur, ilq, ilz, n, istart2, istop, nw, a, lda,
656 $ b, ldb, q, ldq, z, ldz, n_undeflated, n_deflated,
657 $ alphar, alphai, beta, work, nw, work( nw**2+1 ),
658 $ nw, work( 2*nw**2+1 ), lwork-2*nw**2, rec,
661 IF ( n_deflated > 0 )
THEN
662 istop = istop-n_deflated
667 IF ( 100*n_deflated > nibble*( n_deflated+n_undeflated ) .OR.
668 $ istop-istart2+1 .LT. nmin )
THEN
676 ns = min( nshifts, istop-istart2 )
677 ns = min( ns, n_undeflated )
678 shiftpos = istop-n_deflated-n_undeflated+1
683 DO i = shiftpos, shiftpos+n_undeflated-1, 2
684 IF( alphai( i ).NE.-alphai( i+1 ) )
THEN
687 alphar( i ) = alphar( i+1 )
688 alphar( i+1 ) = alphar( i+2 )
692 alphai( i ) = alphai( i+1 )
693 alphai( i+1 ) = alphai( i+2 )
697 beta( i ) = beta( i+1 )
698 beta( i+1 ) = beta( i+2 )
703 IF ( mod( ld, 6 ) .EQ. 0 )
THEN
707 IF( ( real( maxit )*safmin )*abs( a( istop,
708 $ istop-1 ) ).LT.abs( a( istop-1, istop-1 ) ) )
THEN
709 eshift = a( istop, istop-1 )/b( istop-1, istop-1 )
711 eshift = eshift+one/( safmin*real( maxit ) )
713 alphar( shiftpos ) = one
714 alphar( shiftpos+1 ) = zero
715 alphai( shiftpos ) = zero
716 alphai( shiftpos+1 ) = zero
717 beta( shiftpos ) = eshift
718 beta( shiftpos+1 ) = eshift
725 CALL slaqz4( ilschur, ilq, ilz, n, istart2, istop, ns, nblock,
726 $ alphar( shiftpos ), alphai( shiftpos ),
727 $ beta( shiftpos ), a, lda, b, ldb, q, ldq, z, ldz,
728 $ work, nblock, work( nblock**2+1 ), nblock,
729 $ work( 2*nblock**2+1 ), lwork-2*nblock**2,
740 80
CALL shgeqz( wants, wantq, wantz, n, ilo, ihi, a, lda, b, ldb,
741 $ alphar, alphai, beta, q, ldq, z, ldz, work, lwork,
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 slartg(f, g, c, s, r)
SLARTG generates a plane rotation with real cosine and real sine.
integer function ilaenv(ISPEC, NAME, OPTS, N1, N2, N3, N4)
ILAENV
subroutine xerbla(SRNAME, INFO)
XERBLA
logical function lsame(CA, CB)
LSAME
subroutine slaqz4(ILSCHUR, ILQ, ILZ, N, ILO, IHI, NSHIFTS, NBLOCK_DESIRED, SR, SI, SS, A, LDA, B, LDB, Q, LDQ, Z, LDZ, QC, LDQC, ZC, LDZC, WORK, LWORK, INFO)
SLAQZ4
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 shgeqz(JOB, COMPQ, COMPZ, N, ILO, IHI, H, LDH, T, LDT, ALPHAR, ALPHAI, BETA, Q, LDQ, Z, LDZ, WORK, LWORK, INFO)
SHGEQZ
real function slanhs(NORM, N, A, LDA, WORK)
SLANHS returns the value of the 1-norm, Frobenius norm, infinity-norm, or the largest absolute value ...
subroutine srot(N, SX, INCX, SY, INCY, C, S)
SROT
real function slamch(CMACH)
SLAMCH