Skip to content

Commit 643f6db

Browse files
weslleyspereiraweslleyspereira
authored
Merge pull request #478 from thijssteel/exceptional-shifts-in-xlahqr
improve exceptional shift in xlahqr
2 parents 3b4c73e + 99fefbe commit 643f6db

4 files changed

Lines changed: 70 additions & 29 deletions

File tree

‎SRC/clahqr.f‎

Lines changed: 17 additions & 7 deletions
Original file line numberDiff line numberDiff line change
@@ -194,6 +194,7 @@
194194
* =====================================================================
195195
SUBROUTINE CLAHQR( WANTT, WANTZ, N, ILO, IHI, H, LDH, W, ILOZ,
196196
$ IHIZ, Z, LDZ, INFO )
197+
IMPLICIT NONE
197198
*
198199
* -- LAPACK auxiliary routine (version 3.7.0) --
199200
* -- LAPACK is a software package provided by Univ. of Tennessee, --
@@ -218,14 +219,16 @@ SUBROUTINE CLAHQR( WANTT, WANTZ, N, ILO, IHI, H, LDH, W, ILOZ,
218219
PARAMETER ( RZERO = 0.0e0, RONE = 1.0e0, HALF = 0.5e0 )
219220
REAL DAT1
220221
PARAMETER ( DAT1 = 3.0e0 / 4.0e0 )
222+
INTEGER KEXSH
223+
PARAMETER ( KEXSH = 6 )
221224
* ..
222225
* .. Local Scalars ..
223226
COMPLEX CDUM, H11, H11S, H22, SC, SUM, T, T1, TEMP, U,
224227
$ V2, X, Y
225228
REAL AA, AB, BA, BB, H10, H21, RTEMP, S, SAFMAX,
226229
$ SAFMIN, SMLNUM, SX, T2, TST, ULP
227230
INTEGER I, I1, I2, ITS, ITMAX, J, JHI, JLO, K, L, M,
228-
$ NH, NZ
231+
$ NH, NZ, KDEFL
229232
* ..
230233
* .. Local Arrays ..
231234
COMPLEX V( 2 )
@@ -315,6 +318,10 @@ SUBROUTINE CLAHQR( WANTT, WANTZ, N, ILO, IHI, H, LDH, W, ILOZ,
315318
*
316319
ITMAX = 30 * MAX( 10, NH )
317320
*
321+
* KDEFL counts the number of iterations since a deflation
322+
*
323+
KDEFL = -2
324+
*
318325
* The main loop begins here. I is the loop index and decreases from
319326
* IHI to ILO in steps of 1. Each iteration of the loop works
320327
* with the active submatrix in rows and columns L to I.
@@ -374,6 +381,7 @@ SUBROUTINE CLAHQR( WANTT, WANTZ, N, ILO, IHI, H, LDH, W, ILOZ,
374381
*
375382
IF( L.GE.I )
376383
$ GO TO 140
384+
KDEFL = KDEFL + 1
377385
*
378386
* Now the active submatrix is in rows and columns L to I. If
379387
* eigenvalues only are being computed, only the active submatrix
@@ -384,18 +392,18 @@ SUBROUTINE CLAHQR( WANTT, WANTZ, N, ILO, IHI, H, LDH, W, ILOZ,
384392
I2 = I
385393
END IF
386394
*
387-
IF( ITS.EQ.10 ) THEN
395+
IF( MOD(KDEFL,2*KEXSH).EQ.0 ) THEN
388396
*
389397
* Exceptional shift.
390398
*
391-
S = DAT1*ABS( REAL( H( L+1, L ) ) )
392-
T = S + H( L, L )
393-
ELSE IF( ITS.EQ.20 ) THEN
399+
S = DAT1*ABS( REAL( H( I, I-1 ) ) )
400+
T = S + H( I, I )
401+
ELSE IF( MOD(KDEFL,KEXSH).EQ.0 ) THEN
394402
*
395403
* Exceptional shift.
396404
*
397-
S = DAT1*ABS( REAL( H( I, I-1 ) ) )
398-
T = S + H( I, I )
405+
S = DAT1*ABS( REAL( H( L+1, L ) ) )
406+
T = S + H( L, L )
399407
ELSE
400408
*
401409
* Wilkinson's shift.
@@ -556,6 +564,8 @@ SUBROUTINE CLAHQR( WANTT, WANTZ, N, ILO, IHI, H, LDH, W, ILOZ,
556564
* H(I,I-1) is negligible: one eigenvalue has converged.
557565
*
558566
W( I ) = H( I, I )
567+
* reset deflation counter
568+
KDEFL = 0
559569
*
560570
* return to start of the main loop with new value of I.
561571
*

‎SRC/dlahqr.f‎

Lines changed: 18 additions & 7 deletions
Original file line numberDiff line numberDiff line change
@@ -206,6 +206,7 @@
206206
* =====================================================================
207207
SUBROUTINE DLAHQR( WANTT, WANTZ, N, ILO, IHI, H, LDH, WR, WI,
208208
$ ILOZ, IHIZ, Z, LDZ, INFO )
209+
IMPLICIT NONE
209210
*
210211
* -- LAPACK auxiliary routine (version 3.7.0) --
211212
* -- LAPACK is a software package provided by Univ. of Tennessee, --
@@ -227,13 +228,16 @@ SUBROUTINE DLAHQR( WANTT, WANTZ, N, ILO, IHI, H, LDH, WR, WI,
227228
PARAMETER ( ZERO = 0.0d0, ONE = 1.0d0, TWO = 2.0d0 )
228229
DOUBLE PRECISION DAT1, DAT2
229230
PARAMETER ( DAT1 = 3.0d0 / 4.0d0, DAT2 = -0.4375d0 )
231+
INTEGER KEXSH
232+
PARAMETER ( KEXSH = 6 )
230233
* ..
231234
* .. Local Scalars ..
232235
DOUBLE PRECISION AA, AB, BA, BB, CS, DET, H11, H12, H21, H21S,
233236
$ H22, RT1I, RT1R, RT2I, RT2R, RTDISC, S, SAFMAX,
234237
$ SAFMIN, SMLNUM, SN, SUM, T1, T2, T3, TR, TST,
235238
$ ULP, V2, V3
236-
INTEGER I, I1, I2, ITS, ITMAX, J, K, L, M, NH, NR, NZ
239+
INTEGER I, I1, I2, ITS, ITMAX, J, K, L, M, NH, NR, NZ,
240+
$ KDEFL
237241
* ..
238242
* .. Local Arrays ..
239243
DOUBLE PRECISION V( 3 )
@@ -294,6 +298,10 @@ SUBROUTINE DLAHQR( WANTT, WANTZ, N, ILO, IHI, H, LDH, WR, WI,
294298
*
295299
ITMAX = 30 * MAX( 10, NH )
296300
*
301+
* KDEFL counts the number of iterations since a deflation
302+
*
303+
KDEFL = -2
304+
*
297305
* The main loop begins here. I is the loop index and decreases from
298306
* IHI to ILO in steps of 1 or 2. Each iteration of the loop works
299307
* with the active submatrix in rows and columns L to I.
@@ -353,6 +361,7 @@ SUBROUTINE DLAHQR( WANTT, WANTZ, N, ILO, IHI, H, LDH, WR, WI,
353361
*
354362
IF( L.GE.I-1 )
355363
$ GO TO 150
364+
KDEFL = KDEFL + 1
356365
*
357366
* Now the active submatrix is in rows and columns L to I. If
358367
* eigenvalues only are being computed, only the active submatrix
@@ -363,21 +372,21 @@ SUBROUTINE DLAHQR( WANTT, WANTZ, N, ILO, IHI, H, LDH, WR, WI,
363372
I2 = I
364373
END IF
365374
*
366-
IF( ITS.EQ.10 ) THEN
375+
IF( MOD(KDEFL,2*KEXSH).EQ.0 ) THEN
367376
*
368377
* Exceptional shift.
369378
*
370-
S = ABS( H( L+1, L ) ) + ABS( H( L+2, L+1 ) )
371-
H11 = DAT1*S + H( L, L )
379+
S = ABS( H( I, I-1 ) ) + ABS( H( I-1, I-2 ) )
380+
H11 = DAT1*S + H( I, I )
372381
H12 = DAT2*S
373382
H21 = S
374383
H22 = H11
375-
ELSE IF( ITS.EQ.20 ) THEN
384+
ELSE IF( MOD(KDEFL,KEXSH).EQ.0 ) THEN
376385
*
377386
* Exceptional shift.
378387
*
379-
S = ABS( H( I, I-1 ) ) + ABS( H( I-1, I-2 ) )
380-
H11 = DAT1*S + H( I, I )
388+
S = ABS( H( L+1, L ) ) + ABS( H( L+2, L+1 ) )
389+
H11 = DAT1*S + H( L, L )
381390
H12 = DAT2*S
382391
H21 = S
383392
H22 = H11
@@ -599,6 +608,8 @@ SUBROUTINE DLAHQR( WANTT, WANTZ, N, ILO, IHI, H, LDH, WR, WI,
599608
CALL DROT( NZ, Z( ILOZ, I-1 ), 1, Z( ILOZ, I ), 1, CS, SN )
600609
END IF
601610
END IF
611+
* reset deflation counter
612+
KDEFL = 0
602613
*
603614
* return to start of the main loop with new value of I.
604615
*

‎SRC/slahqr.f‎

Lines changed: 18 additions & 8 deletions
Original file line numberDiff line numberDiff line change
@@ -206,6 +206,7 @@
206206
* =====================================================================
207207
SUBROUTINE SLAHQR( WANTT, WANTZ, N, ILO, IHI, H, LDH, WR, WI,
208208
$ ILOZ, IHIZ, Z, LDZ, INFO )
209+
IMPLICIT NONE
209210
*
210211
* -- LAPACK auxiliary routine (version 3.7.0) --
211212
* -- LAPACK is a software package provided by Univ. of Tennessee, --
@@ -227,13 +228,16 @@ SUBROUTINE SLAHQR( WANTT, WANTZ, N, ILO, IHI, H, LDH, WR, WI,
227228
PARAMETER ( ZERO = 0.0e0, ONE = 1.0e0, TWO = 2.0e0 )
228229
REAL DAT1, DAT2
229230
PARAMETER ( DAT1 = 3.0e0 / 4.0e0, DAT2 = -0.4375e0 )
231+
INTEGER KEXSH
232+
PARAMETER ( KEXSH = 6 )
230233
* ..
231234
* .. Local Scalars ..
232235
REAL AA, AB, BA, BB, CS, DET, H11, H12, H21, H21S,
233236
$ H22, RT1I, RT1R, RT2I, RT2R, RTDISC, S, SAFMAX,
234237
$ SAFMIN, SMLNUM, SN, SUM, T1, T2, T3, TR, TST,
235238
$ ULP, V2, V3
236-
INTEGER I, I1, I2, ITS, ITMAX, J, K, L, M, NH, NR, NZ
239+
INTEGER I, I1, I2, ITS, ITMAX, J, K, L, M, NH, NR, NZ,
240+
$ KDEFL
237241
* ..
238242
* .. Local Arrays ..
239243
REAL V( 3 )
@@ -294,6 +298,10 @@ SUBROUTINE SLAHQR( WANTT, WANTZ, N, ILO, IHI, H, LDH, WR, WI,
294298
*
295299
ITMAX = 30 * MAX( 10, NH )
296300
*
301+
* KDEFL counts the number of iterations since a deflation
302+
*
303+
KDEFL = -2
304+
*
297305
* The main loop begins here. I is the loop index and decreases from
298306
* IHI to ILO in steps of 1 or 2. Each iteration of the loop works
299307
* with the active submatrix in rows and columns L to I.
@@ -353,6 +361,7 @@ SUBROUTINE SLAHQR( WANTT, WANTZ, N, ILO, IHI, H, LDH, WR, WI,
353361
*
354362
IF( L.GE.I-1 )
355363
$ GO TO 150
364+
KDEFL = KDEFL + 1
356365
*
357366
* Now the active submatrix is in rows and columns L to I. If
358367
* eigenvalues only are being computed, only the active submatrix
@@ -363,25 +372,24 @@ SUBROUTINE SLAHQR( WANTT, WANTZ, N, ILO, IHI, H, LDH, WR, WI,
363372
I2 = I
364373
END IF
365374
*
366-
IF( ITS.EQ.10 ) THEN
375+
IF( MOD(KDEFL,2*KEXSH).EQ.0 ) THEN
367376
*
368377
* Exceptional shift.
369378
*
370-
S = ABS( H( L+1, L ) ) + ABS( H( L+2, L+1 ) )
371-
H11 = DAT1*S + H( L, L )
379+
S = ABS( H( I, I-1 ) ) + ABS( H( I-1, I-2 ) )
380+
H11 = DAT1*S + H( I, I )
372381
H12 = DAT2*S
373382
H21 = S
374383
H22 = H11
375-
ELSE IF( ITS.EQ.20 ) THEN
384+
ELSE IF( MOD(KDEFL,KEXSH).EQ.0 ) THEN
376385
*
377386
* Exceptional shift.
378387
*
379-
S = ABS( H( I, I-1 ) ) + ABS( H( I-1, I-2 ) )
380-
H11 = DAT1*S + H( I, I )
388+
S = ABS( H( L+1, L ) ) + ABS( H( L+2, L+1 ) )
389+
H11 = DAT1*S + H( L, L )
381390
H12 = DAT2*S
382391
H21 = S
383392
H22 = H11
384-
ELSE
385393
*
386394
* Prepare to use Francis' double shift
387395
* (i.e. 2nd degree generalized Rayleigh quotient)
@@ -599,6 +607,8 @@ SUBROUTINE SLAHQR( WANTT, WANTZ, N, ILO, IHI, H, LDH, WR, WI,
599607
CALL SROT( NZ, Z( ILOZ, I-1 ), 1, Z( ILOZ, I ), 1, CS, SN )
600608
END IF
601609
END IF
610+
* reset deflation counter
611+
KDEFL = 0
602612
*
603613
* return to start of the main loop with new value of I.
604614
*

‎SRC/zlahqr.f‎

Lines changed: 17 additions & 7 deletions
Original file line numberDiff line numberDiff line change
@@ -194,6 +194,7 @@
194194
* =====================================================================
195195
SUBROUTINE ZLAHQR( WANTT, WANTZ, N, ILO, IHI, H, LDH, W, ILOZ,
196196
$ IHIZ, Z, LDZ, INFO )
197+
IMPLICIT NONE
197198
*
198199
* -- LAPACK auxiliary routine (version 3.7.0) --
199200
* -- LAPACK is a software package provided by Univ. of Tennessee, --
@@ -218,14 +219,16 @@ SUBROUTINE ZLAHQR( WANTT, WANTZ, N, ILO, IHI, H, LDH, W, ILOZ,
218219
PARAMETER ( RZERO = 0.0d0, RONE = 1.0d0, HALF = 0.5d0 )
219220
DOUBLE PRECISION DAT1
220221
PARAMETER ( DAT1 = 3.0d0 / 4.0d0 )
222+
INTEGER KEXSH
223+
PARAMETER ( KEXSH = 6 )
221224
* ..
222225
* .. Local Scalars ..
223226
COMPLEX*16 CDUM, H11, H11S, H22, SC, SUM, T, T1, TEMP, U,
224227
$ V2, X, Y
225228
DOUBLE PRECISION AA, AB, BA, BB, H10, H21, RTEMP, S, SAFMAX,
226229
$ SAFMIN, SMLNUM, SX, T2, TST, ULP
227230
INTEGER I, I1, I2, ITS, ITMAX, J, JHI, JLO, K, L, M,
228-
$ NH, NZ
231+
$ NH, NZ, KDEFL
229232
* ..
230233
* .. Local Arrays ..
231234
COMPLEX*16 V( 2 )
@@ -315,6 +318,10 @@ SUBROUTINE ZLAHQR( WANTT, WANTZ, N, ILO, IHI, H, LDH, W, ILOZ,
315318
*
316319
ITMAX = 30 * MAX( 10, NH )
317320
*
321+
* KDEFL counts the number of iterations since a deflation
322+
*
323+
KDEFL = -2
324+
*
318325
* The main loop begins here. I is the loop index and decreases from
319326
* IHI to ILO in steps of 1. Each iteration of the loop works
320327
* with the active submatrix in rows and columns L to I.
@@ -374,6 +381,7 @@ SUBROUTINE ZLAHQR( WANTT, WANTZ, N, ILO, IHI, H, LDH, W, ILOZ,
374381
*
375382
IF( L.GE.I )
376383
$ GO TO 140
384+
KDEFL = KDEFL + 1
377385
*
378386
* Now the active submatrix is in rows and columns L to I. If
379387
* eigenvalues only are being computed, only the active submatrix
@@ -384,18 +392,18 @@ SUBROUTINE ZLAHQR( WANTT, WANTZ, N, ILO, IHI, H, LDH, W, ILOZ,
384392
I2 = I
385393
END IF
386394
*
387-
IF( ITS.EQ.10 ) THEN
395+
IF( MOD(KDEFL,2*KEXSH).EQ.0 ) THEN
388396
*
389397
* Exceptional shift.
390398
*
391-
S = DAT1*ABS( DBLE( H( L+1, L ) ) )
392-
T = S + H( L, L )
393-
ELSE IF( ITS.EQ.20 ) THEN
399+
S = DAT1*ABS( DBLE( H( I, I-1 ) ) )
400+
T = S + H( I, I )
401+
ELSE IF( MOD(KDEFL,KEXSH).EQ.0 ) THEN
394402
*
395403
* Exceptional shift.
396404
*
397-
S = DAT1*ABS( DBLE( H( I, I-1 ) ) )
398-
T = S + H( I, I )
405+
S = DAT1*ABS( DBLE( H( L+1, L ) ) )
406+
T = S + H( L, L )
399407
ELSE
400408
*
401409
* Wilkinson's shift.
@@ -557,6 +565,8 @@ SUBROUTINE ZLAHQR( WANTT, WANTZ, N, ILO, IHI, H, LDH, W, ILOZ,
557565
* H(I,I-1) is negligible: one eigenvalue has converged.
558566
*
559567
W( I ) = H( I, I )
568+
* reset deflation counter
569+
KDEFL = 0
560570
*
561571
* return to start of the main loop with new value of I.
562572
*

0 commit comments

Comments
 (0)