Skip to content

Commit 34b804e

Browse files
langoulangou
authored
Merge pull request #698 from thijssteel/qz-revert-infinite-eigenvalue-criterion
Use normwise criterion for inf eigenvalues in QZ
2 parents 64023f5 + e703366 commit 34b804e

File tree

8 files changed

+34
-68
lines changed

8 files changed

+34
-68
lines changed

‎SRC/chgeqz.f

Lines changed: 2 additions & 7 deletions
Original file line numberDiff line numberDiff line change
@@ -523,9 +523,7 @@ SUBROUTINE CHGEQZ( JOB, COMPQ, COMPZ, N, ILO, IHI, H, LDH, T, LDT,
523523
END IF
524524
END IF
525525
*
526-
IF( ABS( T( ILAST, ILAST ) ).LE.MAX( SAFMIN, ULP*(
527-
$ ABS( T( ILAST - 1, ILAST ) ) + ABS( T( ILAST-1, ILAST-1 )
528-
$ ) ) ) ) THEN
526+
IF( ABS( T( ILAST, ILAST ) ).LE.BTOL ) THEN
529527
T( ILAST, ILAST ) = CZERO
530528
GO TO 50
531529
END IF
@@ -551,10 +549,7 @@ SUBROUTINE CHGEQZ( JOB, COMPQ, COMPZ, N, ILO, IHI, H, LDH, T, LDT,
551549
*
552550
* Test 2: for T(j,j)=0
553551
*
554-
TEMP = ABS ( T( J, J + 1 ) )
555-
IF ( J .GT. ILO )
556-
$ TEMP = TEMP + ABS ( T( J - 1, J ) )
557-
IF( ABS( T( J, J ) ).LT.MAX( SAFMIN,ULP*TEMP ) ) THEN
552+
IF( ABS( T( J, J ) ).LT.BTOL ) THEN
558553
T( J, J ) = CZERO
559554
*
560555
* Test 1a: Check for 2 consecutive small subdiagonals in A

‎SRC/claqz0.f

Lines changed: 6 additions & 10 deletions
Original file line numberDiff line numberDiff line change
@@ -299,7 +299,7 @@ RECURSIVE SUBROUTINE CLAQZ0( WANTS, WANTQ, WANTZ, N, ILO, IHI, A,
299299
PARAMETER( ZERO = 0.0, ONE = 1.0, HALF = 0.5 )
300300

301301
* Local scalars
302-
REAL :: SMLNUM, ULP, SAFMIN, SAFMAX, C1, TEMPR
302+
REAL :: SMLNUM, ULP, SAFMIN, SAFMAX, C1, TEMPR, BNORM, BTOL
303303
COMPLEX :: ESHIFT, S1, TEMP
304304
INTEGER :: ISTART, ISTOP, IITER, MAXIT, ISTART2, K, LD, NSHIFTS,
305305
$ NBLOCK, NW, NMIN, NIBBLE, N_UNDEFLATED, N_DEFLATED,
@@ -312,7 +312,7 @@ RECURSIVE SUBROUTINE CLAQZ0( WANTS, WANTQ, WANTZ, N, ILO, IHI, A,
312312
* External Functions
313313
EXTERNAL :: XERBLA, CHGEQZ, CLAQZ2, CLAQZ3, CLASET, SLABAD,
314314
$ CLARTG, CROT
315-
REAL, EXTERNAL :: SLAMCH
315+
REAL, EXTERNAL :: SLAMCH, CLANHS
316316
LOGICAL, EXTERNAL :: LSAME
317317
INTEGER, EXTERNAL :: ILAENV
318318

@@ -466,6 +466,9 @@ RECURSIVE SUBROUTINE CLAQZ0( WANTS, WANTQ, WANTZ, N, ILO, IHI, A,
466466
ULP = SLAMCH( 'PRECISION' )
467467
SMLNUM = SAFMIN*( REAL( N )/ULP )
468468

469+
BNORM = CLANHS( 'F', IHI-ILO+1, B( ILO, ILO ), LDB, RWORK )
470+
BTOL = MAX( SAFMIN, ULP*BNORM )
471+
469472
ISTART = ILO
470473
ISTOP = IHI
471474
MAXIT = 30*( IHI-ILO+1 )
@@ -528,15 +531,8 @@ RECURSIVE SUBROUTINE CLAQZ0( WANTS, WANTQ, WANTZ, N, ILO, IHI, A,
528531
* slow down the method when many infinite eigenvalues are present
529532
K = ISTOP
530533
DO WHILE ( K.GE.ISTART2 )
531-
TEMPR = ZERO
532-
IF( K .LT. ISTOP ) THEN
533-
TEMPR = TEMPR+ABS( B( K, K+1 ) )
534-
END IF
535-
IF( K .GT. ISTART2 ) THEN
536-
TEMPR = TEMPR+ABS( B( K-1, K ) )
537-
END IF
538534

539-
IF( ABS( B( K, K ) ) .LT. MAX( SMLNUM, ULP*TEMPR ) ) THEN
535+
IF( ABS( B( K, K ) ) .LT. BTOL ) THEN
540536
* A diagonal element of B is negligable, move it
541537
* to the top and deflate it
542538

‎SRC/dhgeqz.f

Lines changed: 2 additions & 7 deletions
Original file line numberDiff line numberDiff line change
@@ -536,9 +536,7 @@ SUBROUTINE DHGEQZ( JOB, COMPQ, COMPZ, N, ILO, IHI, H, LDH, T, LDT,
536536
END IF
537537
END IF
538538
*
539-
IF( ABS( T( ILAST, ILAST ) ).LE.MAX( SAFMIN, ULP*(
540-
$ ABS( T( ILAST - 1, ILAST ) ) + ABS( T( ILAST-1, ILAST-1 )
541-
$ ) ) ) ) THEN
539+
IF( ABS( T( ILAST, ILAST ) ).LE.BTOL ) THEN
542540
T( ILAST, ILAST ) = ZERO
543541
GO TO 70
544542
END IF
@@ -564,10 +562,7 @@ SUBROUTINE DHGEQZ( JOB, COMPQ, COMPZ, N, ILO, IHI, H, LDH, T, LDT,
564562
*
565563
* Test 2: for T(j,j)=0
566564
*
567-
TEMP = ABS ( T( J, J + 1 ) )
568-
IF ( J .GT. ILO )
569-
$ TEMP = TEMP + ABS ( T( J - 1, J ) )
570-
IF( ABS( T( J, J ) ).LT.MAX( SAFMIN,ULP*TEMP ) ) THEN
565+
IF( ABS( T( J, J ) ).LT.BTOL ) THEN
571566
T( J, J ) = ZERO
572567
*
573568
* Test 1a: Check for 2 consecutive small subdiagonals in A

‎SRC/dlaqz0.f

Lines changed: 6 additions & 10 deletions
Original file line numberDiff line numberDiff line change
@@ -322,7 +322,7 @@ RECURSIVE SUBROUTINE DLAQZ0( WANTS, WANTQ, WANTZ, N, ILO, IHI, A,
322322

323323
* Local scalars
324324
DOUBLE PRECISION :: SMLNUM, ULP, ESHIFT, SAFMIN, SAFMAX, C1, S1,
325-
$ TEMP, SWAP
325+
$ TEMP, SWAP, BNORM, BTOL
326326
INTEGER :: ISTART, ISTOP, IITER, MAXIT, ISTART2, K, LD, NSHIFTS,
327327
$ NBLOCK, NW, NMIN, NIBBLE, N_UNDEFLATED, N_DEFLATED,
328328
$ NS, SWEEP_INFO, SHIFTPOS, LWORKREQ, K2, ISTARTM,
@@ -334,7 +334,7 @@ RECURSIVE SUBROUTINE DLAQZ0( WANTS, WANTQ, WANTZ, N, ILO, IHI, A,
334334
* External Functions
335335
EXTERNAL :: XERBLA, DHGEQZ, DLASET, DLAQZ3, DLAQZ4, DLABAD,
336336
$ DLARTG, DROT
337-
DOUBLE PRECISION, EXTERNAL :: DLAMCH
337+
DOUBLE PRECISION, EXTERNAL :: DLAMCH, DLANHS
338338
LOGICAL, EXTERNAL :: LSAME
339339
INTEGER, EXTERNAL :: ILAENV
340340

@@ -486,6 +486,9 @@ RECURSIVE SUBROUTINE DLAQZ0( WANTS, WANTQ, WANTZ, N, ILO, IHI, A,
486486
ULP = DLAMCH( 'PRECISION' )
487487
SMLNUM = SAFMIN*( DBLE( N )/ULP )
488488

489+
BNORM = DLANHS( 'F', IHI-ILO+1, B( ILO, ILO ), LDB, WORK )
490+
BTOL = MAX( SAFMIN, ULP*BNORM )
491+
489492
ISTART = ILO
490493
ISTOP = IHI
491494
MAXIT = 3*( IHI-ILO+1 )
@@ -562,15 +565,8 @@ RECURSIVE SUBROUTINE DLAQZ0( WANTS, WANTQ, WANTZ, N, ILO, IHI, A,
562565
* slow down the method when many infinite eigenvalues are present
563566
K = ISTOP
564567
DO WHILE ( K.GE.ISTART2 )
565-
TEMP = ZERO
566-
IF( K .LT. ISTOP ) THEN
567-
TEMP = TEMP+ABS( B( K, K+1 ) )
568-
END IF
569-
IF( K .GT. ISTART2 ) THEN
570-
TEMP = TEMP+ABS( B( K-1, K ) )
571-
END IF
572568

573-
IF( ABS( B( K, K ) ) .LT. MAX( SMLNUM, ULP*TEMP ) ) THEN
569+
IF( ABS( B( K, K ) ) .LT. BTOL ) THEN
574570
* A diagonal element of B is negligable, move it
575571
* to the top and deflate it
576572

‎SRC/shgeqz.f

Lines changed: 2 additions & 7 deletions
Original file line numberDiff line numberDiff line change
@@ -536,9 +536,7 @@ SUBROUTINE SHGEQZ( JOB, COMPQ, COMPZ, N, ILO, IHI, H, LDH, T, LDT,
536536
END IF
537537
END IF
538538
*
539-
IF( ABS( T( ILAST, ILAST ) ).LE.MAX( SAFMIN, ULP*(
540-
$ ABS( T( ILAST - 1, ILAST ) ) + ABS( T( ILAST-1, ILAST-1 )
541-
$ ) ) ) ) THEN
539+
IF( ABS( T( ILAST, ILAST ) ).LE.BTOL ) THEN
542540
T( ILAST, ILAST ) = ZERO
543541
GO TO 70
544542
END IF
@@ -564,10 +562,7 @@ SUBROUTINE SHGEQZ( JOB, COMPQ, COMPZ, N, ILO, IHI, H, LDH, T, LDT,
564562
*
565563
* Test 2: for T(j,j)=0
566564
*
567-
TEMP = ABS ( T( J, J + 1 ) )
568-
IF ( J .GT. ILO )
569-
$ TEMP = TEMP + ABS ( T( J - 1, J ) )
570-
IF( ABS( T( J, J ) ).LT.MAX( SAFMIN,ULP*TEMP ) ) THEN
565+
IF( ABS( T( J, J ) ).LT.BTOL ) THEN
571566
T( J, J ) = ZERO
572567
*
573568
* Test 1a: Check for 2 consecutive small subdiagonals in A

‎SRC/slaqz0.f

Lines changed: 7 additions & 10 deletions
Original file line numberDiff line numberDiff line change
@@ -318,7 +318,8 @@ RECURSIVE SUBROUTINE SLAQZ0( WANTS, WANTQ, WANTZ, N, ILO, IHI, A,
318318
PARAMETER( ZERO = 0.0, ONE = 1.0, HALF = 0.5 )
319319

320320
* Local scalars
321-
REAL :: SMLNUM, ULP, ESHIFT, SAFMIN, SAFMAX, C1, S1, TEMP, SWAP
321+
REAL :: SMLNUM, ULP, ESHIFT, SAFMIN, SAFMAX, C1, S1, TEMP, SWAP,
322+
$ BNORM, BTOL
322323
INTEGER :: ISTART, ISTOP, IITER, MAXIT, ISTART2, K, LD, NSHIFTS,
323324
$ NBLOCK, NW, NMIN, NIBBLE, N_UNDEFLATED, N_DEFLATED,
324325
$ NS, SWEEP_INFO, SHIFTPOS, LWORKREQ, K2, ISTARTM,
@@ -330,7 +331,7 @@ RECURSIVE SUBROUTINE SLAQZ0( WANTS, WANTQ, WANTZ, N, ILO, IHI, A,
330331
* External Functions
331332
EXTERNAL :: XERBLA, SHGEQZ, SLAQZ3, SLAQZ4, SLASET, SLABAD,
332333
$ SLARTG, SROT
333-
REAL, EXTERNAL :: SLAMCH
334+
REAL, EXTERNAL :: SLAMCH, SLANHS
334335
LOGICAL, EXTERNAL :: LSAME
335336
INTEGER, EXTERNAL :: ILAENV
336337

@@ -482,6 +483,9 @@ RECURSIVE SUBROUTINE SLAQZ0( WANTS, WANTQ, WANTZ, N, ILO, IHI, A,
482483
ULP = SLAMCH( 'PRECISION' )
483484
SMLNUM = SAFMIN*( REAL( N )/ULP )
484485

486+
BNORM = SLANHS( 'F', IHI-ILO+1, B( ILO, ILO ), LDB, WORK )
487+
BTOL = MAX( SAFMIN, ULP*BNORM )
488+
485489
ISTART = ILO
486490
ISTOP = IHI
487491
MAXIT = 3*( IHI-ILO+1 )
@@ -558,15 +562,8 @@ RECURSIVE SUBROUTINE SLAQZ0( WANTS, WANTQ, WANTZ, N, ILO, IHI, A,
558562
* slow down the method when many infinite eigenvalues are present
559563
K = ISTOP
560564
DO WHILE ( K.GE.ISTART2 )
561-
TEMP = ZERO
562-
IF( K .LT. ISTOP ) THEN
563-
TEMP = TEMP+ABS( B( K, K+1 ) )
564-
END IF
565-
IF( K .GT. ISTART2 ) THEN
566-
TEMP = TEMP+ABS( B( K-1, K ) )
567-
END IF
568565

569-
IF( ABS( B( K, K ) ) .LT. MAX( SMLNUM, ULP*TEMP ) ) THEN
566+
IF( ABS( B( K, K ) ) .LT. BTOL ) THEN
570567
* A diagonal element of B is negligable, move it
571568
* to the top and deflate it
572569

‎SRC/zhgeqz.f

Lines changed: 2 additions & 7 deletions
Original file line numberDiff line numberDiff line change
@@ -524,9 +524,7 @@ SUBROUTINE ZHGEQZ( JOB, COMPQ, COMPZ, N, ILO, IHI, H, LDH, T, LDT,
524524
END IF
525525
END IF
526526
*
527-
IF( ABS( T( ILAST, ILAST ) ).LE.MAX( SAFMIN, ULP*(
528-
$ ABS( T( ILAST - 1, ILAST ) ) + ABS( T( ILAST-1, ILAST-1 )
529-
$ ) ) ) ) THEN
527+
IF( ABS( T( ILAST, ILAST ) ).LE.BTOL ) THEN
530528
T( ILAST, ILAST ) = CZERO
531529
GO TO 50
532530
END IF
@@ -552,10 +550,7 @@ SUBROUTINE ZHGEQZ( JOB, COMPQ, COMPZ, N, ILO, IHI, H, LDH, T, LDT,
552550
*
553551
* Test 2: for T(j,j)=0
554552
*
555-
TEMP = ABS ( T( J, J + 1 ) )
556-
IF ( J .GT. ILO )
557-
$ TEMP = TEMP + ABS ( T( J - 1, J ) )
558-
IF( ABS( T( J, J ) ).LT.MAX( SAFMIN,ULP*TEMP ) ) THEN
553+
IF( ABS( T( J, J ) ).LT.BTOL ) THEN
559554
T( J, J ) = CZERO
560555
*
561556
* Test 1a: Check for 2 consecutive small subdiagonals in A

‎SRC/zlaqz0.f

Lines changed: 7 additions & 10 deletions
Original file line numberDiff line numberDiff line change
@@ -300,7 +300,8 @@ RECURSIVE SUBROUTINE ZLAQZ0( WANTS, WANTQ, WANTZ, N, ILO, IHI, A,
300300
PARAMETER( ZERO = 0.0D0, ONE = 1.0D0, HALF = 0.5D0 )
301301

302302
* Local scalars
303-
DOUBLE PRECISION :: SMLNUM, ULP, SAFMIN, SAFMAX, C1, TEMPR
303+
DOUBLE PRECISION :: SMLNUM, ULP, SAFMIN, SAFMAX, C1, TEMPR,
304+
$ BNORM, BTOL
304305
COMPLEX*16 :: ESHIFT, S1, TEMP
305306
INTEGER :: ISTART, ISTOP, IITER, MAXIT, ISTART2, K, LD, NSHIFTS,
306307
$ NBLOCK, NW, NMIN, NIBBLE, N_UNDEFLATED, N_DEFLATED,
@@ -313,7 +314,7 @@ RECURSIVE SUBROUTINE ZLAQZ0( WANTS, WANTQ, WANTZ, N, ILO, IHI, A,
313314
* External Functions
314315
EXTERNAL :: XERBLA, ZHGEQZ, ZLAQZ2, ZLAQZ3, ZLASET, DLABAD,
315316
$ ZLARTG, ZROT
316-
DOUBLE PRECISION, EXTERNAL :: DLAMCH
317+
DOUBLE PRECISION, EXTERNAL :: DLAMCH, ZLANHS
317318
LOGICAL, EXTERNAL :: LSAME
318319
INTEGER, EXTERNAL :: ILAENV
319320

@@ -467,6 +468,9 @@ RECURSIVE SUBROUTINE ZLAQZ0( WANTS, WANTQ, WANTZ, N, ILO, IHI, A,
467468
ULP = DLAMCH( 'PRECISION' )
468469
SMLNUM = SAFMIN*( DBLE( N )/ULP )
469470

471+
BNORM = ZLANHS( 'F', IHI-ILO+1, B( ILO, ILO ), LDB, RWORK )
472+
BTOL = MAX( SAFMIN, ULP*BNORM )
473+
470474
ISTART = ILO
471475
ISTOP = IHI
472476
MAXIT = 30*( IHI-ILO+1 )
@@ -529,15 +533,8 @@ RECURSIVE SUBROUTINE ZLAQZ0( WANTS, WANTQ, WANTZ, N, ILO, IHI, A,
529533
* slow down the method when many infinite eigenvalues are present
530534
K = ISTOP
531535
DO WHILE ( K.GE.ISTART2 )
532-
TEMPR = ZERO
533-
IF( K .LT. ISTOP ) THEN
534-
TEMPR = TEMPR+ABS( B( K, K+1 ) )
535-
END IF
536-
IF( K .GT. ISTART2 ) THEN
537-
TEMPR = TEMPR+ABS( B( K-1, K ) )
538-
END IF
539536

540-
IF( ABS( B( K, K ) ) .LT. MAX( SMLNUM, ULP*TEMPR ) ) THEN
537+
IF( ABS( B( K, K ) ) .LT. BTOL ) THEN
541538
* A diagonal element of B is negligable, move it
542539
* to the top and deflate it
543540

0 commit comments

Comments
 (0)