Skip to content

Commit 0d82f77

Browse files
thijssteelthijssteel
committed
improve exceptional shift in xlahqr
1 parent 2b9f0f5 commit 0d82f77

File tree

4 files changed

+66
-29
lines changed

4 files changed

+66
-29
lines changed

‎SRC/clahqr.f

Lines changed: 16 additions & 7 deletions
Original file line numberDiff line numberDiff line change
@@ -218,14 +218,16 @@ SUBROUTINE CLAHQR( WANTT, WANTZ, N, ILO, IHI, H, LDH, W, ILOZ,
218218
PARAMETER ( RZERO = 0.0e0, RONE = 1.0e0, HALF = 0.5e0 )
219219
REAL DAT1
220220
PARAMETER ( DAT1 = 3.0e0 / 4.0e0 )
221+
INTEGER KEXSH
222+
PARAMETER ( KEXSH = 6 )
221223
* ..
222224
* .. Local Scalars ..
223225
COMPLEX CDUM, H11, H11S, H22, SC, SUM, T, T1, TEMP, U,
224226
$ V2, X, Y
225227
REAL AA, AB, BA, BB, H10, H21, RTEMP, S, SAFMAX,
226228
$ SAFMIN, SMLNUM, SX, T2, TST, ULP
227229
INTEGER I, I1, I2, ITS, ITMAX, J, JHI, JLO, K, L, M,
228-
$ NH, NZ
230+
$ NH, NZ, KDEFL
229231
* ..
230232
* .. Local Arrays ..
231233
COMPLEX V( 2 )
@@ -315,6 +317,10 @@ SUBROUTINE CLAHQR( WANTT, WANTZ, N, ILO, IHI, H, LDH, W, ILOZ,
315317
*
316318
ITMAX = 30 * MAX( 10, NH )
317319
*
320+
* KDFL counts the number of iterations since a deflation
321+
*
322+
KDFL = -2
323+
*
318324
* The main loop begins here. I is the loop index and decreases from
319325
* IHI to ILO in steps of 1. Each iteration of the loop works
320326
* with the active submatrix in rows and columns L to I.
@@ -374,6 +380,7 @@ SUBROUTINE CLAHQR( WANTT, WANTZ, N, ILO, IHI, H, LDH, W, ILOZ,
374380
*
375381
IF( L.GE.I )
376382
$ GO TO 140
383+
KDEFL = KDEFL + 1
377384
*
378385
* Now the active submatrix is in rows and columns L to I. If
379386
* eigenvalues only are being computed, only the active submatrix
@@ -384,18 +391,18 @@ SUBROUTINE CLAHQR( WANTT, WANTZ, N, ILO, IHI, H, LDH, W, ILOZ,
384391
I2 = I
385392
END IF
386393
*
387-
IF( ITS.EQ.10 ) THEN
394+
IF( MOD(KDEFL,2*KEXSH).EQ.0 ) THEN
388395
*
389396
* Exceptional shift.
390397
*
391-
S = DAT1*ABS( REAL( H( L+1, L ) ) )
392-
T = S + H( L, L )
393-
ELSE IF( ITS.EQ.20 ) THEN
398+
S = DAT1*ABS( REAL( H( I, I-1 ) ) )
399+
T = S + H( I, I )
400+
ELSE IF( MOD(KDEFL,KEXSH).EQ.0 ) THEN
394401
*
395402
* Exceptional shift.
396403
*
397-
S = DAT1*ABS( REAL( H( I, I-1 ) ) )
398-
T = S + H( I, I )
404+
S = DAT1*ABS( REAL( H( L+1, L ) ) )
405+
T = S + H( L, L )
399406
ELSE
400407
*
401408
* Wilkinson's shift.
@@ -556,6 +563,8 @@ SUBROUTINE CLAHQR( WANTT, WANTZ, N, ILO, IHI, H, LDH, W, ILOZ,
556563
* H(I,I-1) is negligible: one eigenvalue has converged.
557564
*
558565
W( I ) = H( I, I )
566+
* reset deflation counter
567+
KDFL = 0
559568
*
560569
* return to start of the main loop with new value of I.
561570
*

‎SRC/dlahqr.f

Lines changed: 17 additions & 7 deletions
Original file line numberDiff line numberDiff line change
@@ -227,13 +227,16 @@ SUBROUTINE DLAHQR( WANTT, WANTZ, N, ILO, IHI, H, LDH, WR, WI,
227227
PARAMETER ( ZERO = 0.0d0, ONE = 1.0d0, TWO = 2.0d0 )
228228
DOUBLE PRECISION DAT1, DAT2
229229
PARAMETER ( DAT1 = 3.0d0 / 4.0d0, DAT2 = -0.4375d0 )
230+
INTEGER KEXSH
231+
PARAMETER ( KEXSH = 6 )
230232
* ..
231233
* .. Local Scalars ..
232234
DOUBLE PRECISION AA, AB, BA, BB, CS, DET, H11, H12, H21, H21S,
233235
$ H22, RT1I, RT1R, RT2I, RT2R, RTDISC, S, SAFMAX,
234236
$ SAFMIN, SMLNUM, SN, SUM, T1, T2, T3, TR, TST,
235237
$ ULP, V2, V3
236-
INTEGER I, I1, I2, ITS, ITMAX, J, K, L, M, NH, NR, NZ
238+
INTEGER I, I1, I2, ITS, ITMAX, J, K, L, M, NH, NR, NZ,
239+
$ KDEFL
237240
* ..
238241
* .. Local Arrays ..
239242
DOUBLE PRECISION V( 3 )
@@ -294,6 +297,10 @@ SUBROUTINE DLAHQR( WANTT, WANTZ, N, ILO, IHI, H, LDH, WR, WI,
294297
*
295298
ITMAX = 30 * MAX( 10, NH )
296299
*
300+
* KDFL counts the number of iterations since a deflation
301+
*
302+
KDFL = -2
303+
*
297304
* The main loop begins here. I is the loop index and decreases from
298305
* IHI to ILO in steps of 1 or 2. Each iteration of the loop works
299306
* with the active submatrix in rows and columns L to I.
@@ -353,6 +360,7 @@ SUBROUTINE DLAHQR( WANTT, WANTZ, N, ILO, IHI, H, LDH, WR, WI,
353360
*
354361
IF( L.GE.I-1 )
355362
$ GO TO 150
363+
KDEFL = KDEFL + 1
356364
*
357365
* Now the active submatrix is in rows and columns L to I. If
358366
* eigenvalues only are being computed, only the active submatrix
@@ -363,21 +371,21 @@ SUBROUTINE DLAHQR( WANTT, WANTZ, N, ILO, IHI, H, LDH, WR, WI,
363371
I2 = I
364372
END IF
365373
*
366-
IF( ITS.EQ.10 ) THEN
374+
IF( MOD(KDEFL,2*KEXSH).EQ.0 ) THEN
367375
*
368376
* Exceptional shift.
369377
*
370-
S = ABS( H( L+1, L ) ) + ABS( H( L+2, L+1 ) )
371-
H11 = DAT1*S + H( L, L )
378+
S = ABS( H( I, I-1 ) ) + ABS( H( I-1, I-2 ) )
379+
H11 = DAT1*S + H( I, I )
372380
H12 = DAT2*S
373381
H21 = S
374382
H22 = H11
375-
ELSE IF( ITS.EQ.20 ) THEN
383+
ELSE IF( MOD(KDEFL,KEXSH).EQ.0 ) THEN
376384
*
377385
* Exceptional shift.
378386
*
379-
S = ABS( H( I, I-1 ) ) + ABS( H( I-1, I-2 ) )
380-
H11 = DAT1*S + H( I, I )
387+
S = ABS( H( L+1, L ) ) + ABS( H( L+2, L+1 ) )
388+
H11 = DAT1*S + H( L, L )
381389
H12 = DAT2*S
382390
H21 = S
383391
H22 = H11
@@ -599,6 +607,8 @@ SUBROUTINE DLAHQR( WANTT, WANTZ, N, ILO, IHI, H, LDH, WR, WI,
599607
CALL DROT( NZ, Z( ILOZ, I-1 ), 1, Z( ILOZ, I ), 1, CS, SN )
600608
END IF
601609
END IF
610+
* reset deflation counter
611+
KDFL = 0
602612
*
603613
* return to start of the main loop with new value of I.
604614
*

‎SRC/slahqr.f

Lines changed: 17 additions & 8 deletions
Original file line numberDiff line numberDiff line change
@@ -227,13 +227,16 @@ SUBROUTINE SLAHQR( WANTT, WANTZ, N, ILO, IHI, H, LDH, WR, WI,
227227
PARAMETER ( ZERO = 0.0e0, ONE = 1.0e0, TWO = 2.0e0 )
228228
REAL DAT1, DAT2
229229
PARAMETER ( DAT1 = 3.0e0 / 4.0e0, DAT2 = -0.4375e0 )
230+
INTEGER KEXSH
231+
PARAMETER ( KEXSH = 6 )
230232
* ..
231233
* .. Local Scalars ..
232234
REAL AA, AB, BA, BB, CS, DET, H11, H12, H21, H21S,
233235
$ H22, RT1I, RT1R, RT2I, RT2R, RTDISC, S, SAFMAX,
234236
$ SAFMIN, SMLNUM, SN, SUM, T1, T2, T3, TR, TST,
235237
$ ULP, V2, V3
236-
INTEGER I, I1, I2, ITS, ITMAX, J, K, L, M, NH, NR, NZ
238+
INTEGER I, I1, I2, ITS, ITMAX, J, K, L, M, NH, NR, NZ,
239+
$ KDEFL
237240
* ..
238241
* .. Local Arrays ..
239242
REAL V( 3 )
@@ -294,6 +297,10 @@ SUBROUTINE SLAHQR( WANTT, WANTZ, N, ILO, IHI, H, LDH, WR, WI,
294297
*
295298
ITMAX = 30 * MAX( 10, NH )
296299
*
300+
* KDFL counts the number of iterations since a deflation
301+
*
302+
KDFL = -2
303+
*
297304
* The main loop begins here. I is the loop index and decreases from
298305
* IHI to ILO in steps of 1 or 2. Each iteration of the loop works
299306
* with the active submatrix in rows and columns L to I.
@@ -353,6 +360,7 @@ SUBROUTINE SLAHQR( WANTT, WANTZ, N, ILO, IHI, H, LDH, WR, WI,
353360
*
354361
IF( L.GE.I-1 )
355362
$ GO TO 150
363+
KDEFL = KDEFL + 1
356364
*
357365
* Now the active submatrix is in rows and columns L to I. If
358366
* eigenvalues only are being computed, only the active submatrix
@@ -363,25 +371,24 @@ SUBROUTINE SLAHQR( WANTT, WANTZ, N, ILO, IHI, H, LDH, WR, WI,
363371
I2 = I
364372
END IF
365373
*
366-
IF( ITS.EQ.10 ) THEN
374+
IF( MOD(KDEFL,2*KEXSH).EQ.0 ) THEN
367375
*
368376
* Exceptional shift.
369377
*
370-
S = ABS( H( L+1, L ) ) + ABS( H( L+2, L+1 ) )
371-
H11 = DAT1*S + H( L, L )
378+
S = ABS( H( I, I-1 ) ) + ABS( H( I-1, I-2 ) )
379+
H11 = DAT1*S + H( I, I )
372380
H12 = DAT2*S
373381
H21 = S
374382
H22 = H11
375-
ELSE IF( ITS.EQ.20 ) THEN
383+
ELSE IF( MOD(KDEFL,KEXSH).EQ.0 ) THEN
376384
*
377385
* Exceptional shift.
378386
*
379-
S = ABS( H( I, I-1 ) ) + ABS( H( I-1, I-2 ) )
380-
H11 = DAT1*S + H( I, I )
387+
S = ABS( H( L+1, L ) ) + ABS( H( L+2, L+1 ) )
388+
H11 = DAT1*S + H( L, L )
381389
H12 = DAT2*S
382390
H21 = S
383391
H22 = H11
384-
ELSE
385392
*
386393
* Prepare to use Francis' double shift
387394
* (i.e. 2nd degree generalized Rayleigh quotient)
@@ -599,6 +606,8 @@ SUBROUTINE SLAHQR( WANTT, WANTZ, N, ILO, IHI, H, LDH, WR, WI,
599606
CALL SROT( NZ, Z( ILOZ, I-1 ), 1, Z( ILOZ, I ), 1, CS, SN )
600607
END IF
601608
END IF
609+
* reset deflation counter
610+
KDFL = 0
602611
*
603612
* return to start of the main loop with new value of I.
604613
*

‎SRC/zlahqr.f

Lines changed: 16 additions & 7 deletions
Original file line numberDiff line numberDiff line change
@@ -218,14 +218,16 @@ SUBROUTINE ZLAHQR( WANTT, WANTZ, N, ILO, IHI, H, LDH, W, ILOZ,
218218
PARAMETER ( RZERO = 0.0d0, RONE = 1.0d0, HALF = 0.5d0 )
219219
DOUBLE PRECISION DAT1
220220
PARAMETER ( DAT1 = 3.0d0 / 4.0d0 )
221+
INTEGER KEXSH
222+
PARAMETER ( KEXSH = 6 )
221223
* ..
222224
* .. Local Scalars ..
223225
COMPLEX*16 CDUM, H11, H11S, H22, SC, SUM, T, T1, TEMP, U,
224226
$ V2, X, Y
225227
DOUBLE PRECISION AA, AB, BA, BB, H10, H21, RTEMP, S, SAFMAX,
226228
$ SAFMIN, SMLNUM, SX, T2, TST, ULP
227229
INTEGER I, I1, I2, ITS, ITMAX, J, JHI, JLO, K, L, M,
228-
$ NH, NZ
230+
$ NH, NZ, KDEFL
229231
* ..
230232
* .. Local Arrays ..
231233
COMPLEX*16 V( 2 )
@@ -315,6 +317,10 @@ SUBROUTINE ZLAHQR( WANTT, WANTZ, N, ILO, IHI, H, LDH, W, ILOZ,
315317
*
316318
ITMAX = 30 * MAX( 10, NH )
317319
*
320+
* KDFL counts the number of iterations since a deflation
321+
*
322+
KDFL = -2
323+
*
318324
* The main loop begins here. I is the loop index and decreases from
319325
* IHI to ILO in steps of 1. Each iteration of the loop works
320326
* with the active submatrix in rows and columns L to I.
@@ -374,6 +380,7 @@ SUBROUTINE ZLAHQR( WANTT, WANTZ, N, ILO, IHI, H, LDH, W, ILOZ,
374380
*
375381
IF( L.GE.I )
376382
$ GO TO 140
383+
KDEFL = KDEFL + 1
377384
*
378385
* Now the active submatrix is in rows and columns L to I. If
379386
* eigenvalues only are being computed, only the active submatrix
@@ -384,18 +391,18 @@ SUBROUTINE ZLAHQR( WANTT, WANTZ, N, ILO, IHI, H, LDH, W, ILOZ,
384391
I2 = I
385392
END IF
386393
*
387-
IF( ITS.EQ.10 ) THEN
394+
IF( MOD(KDEFL,2*KEXSH).EQ.0 ) THEN
388395
*
389396
* Exceptional shift.
390397
*
391-
S = DAT1*ABS( DBLE( H( L+1, L ) ) )
392-
T = S + H( L, L )
393-
ELSE IF( ITS.EQ.20 ) THEN
398+
S = DAT1*ABS( DBLE( H( I, I-1 ) ) )
399+
T = S + H( I, I )
400+
ELSE IF( MOD(KDEFL,KEXSH).EQ.0 ) THEN
394401
*
395402
* Exceptional shift.
396403
*
397-
S = DAT1*ABS( DBLE( H( I, I-1 ) ) )
398-
T = S + H( I, I )
404+
S = DAT1*ABS( DBLE( H( L+1, L ) ) )
405+
T = S + H( L, L )
399406
ELSE
400407
*
401408
* Wilkinson's shift.
@@ -557,6 +564,8 @@ SUBROUTINE ZLAHQR( WANTT, WANTZ, N, ILO, IHI, H, LDH, W, ILOZ,
557564
* H(I,I-1) is negligible: one eigenvalue has converged.
558565
*
559566
W( I ) = H( I, I )
567+
* reset deflation counter
568+
KDFL = 0
560569
*
561570
* return to start of the main loop with new value of I.
562571
*

0 commit comments

Comments
 (0)