Skip to content

Commit 128b0a9

Browse files
aravindh-krishnamoorthyaravindh-krishnamoorthy
committed
Issue #826: Initial version of the fix in c/d/s/z/stemr.f for 2x2 matrices.
1 parent 326b78e commit 128b0a9

File tree

4 files changed

+92
-25
lines changed

4 files changed

+92
-25
lines changed

‎SRC/cstemr.f‎

Lines changed: 28 additions & 6 deletions
Original file line numberDiff line numberDiff line change
@@ -329,7 +329,8 @@
329329
*> Jim Demmel, University of California, Berkeley, USA \n
330330
*> Inderjit Dhillon, University of Texas, Austin, USA \n
331331
*> Osni Marques, LBNL/NERSC, USA \n
332-
*> Christof Voemel, University of California, Berkeley, USA
332+
*> Christof Voemel, University of California, Berkeley, USA \n
333+
*> Aravindh Krishnamoorthy, FAU, Erlangen, Germany \n
333334
*
334335
* =====================================================================
335336
SUBROUTINE CSTEMR( JOBZ, RANGE, N, D, E, VL, VU, IL, IU,
@@ -361,7 +362,8 @@ SUBROUTINE CSTEMR( JOBZ, RANGE, N, D, E, VL, VU, IL, IU,
361362
$ MINRGP = 3.0E-3 )
362363
* ..
363364
* .. Local Scalars ..
364-
LOGICAL ALLEIG, INDEIG, LQUERY, VALEIG, WANTZ, ZQUERY
365+
LOGICAL ALLEIG, INDEIG, LQUERY, VALEIG, WANTZ, ZQUERY,
366+
$ LAESWAP
365367
INTEGER I, IBEGIN, IEND, IFIRST, IIL, IINDBL, IINDW,
366368
$ IINDWK, IINFO, IINSPL, IIU, ILAST, IN, INDD,
367369
$ INDE2, INDERR, INDGP, INDGRS, INDWRK, ITMP,
@@ -397,6 +399,7 @@ SUBROUTINE CSTEMR( JOBZ, RANGE, N, D, E, VL, VU, IL, IU,
397399
*
398400
LQUERY = ( ( LWORK.EQ.-1 ).OR.( LIWORK.EQ.-1 ) )
399401
ZQUERY = ( NZC.EQ.-1 )
402+
LAESWAP = .FALSE.
400403

401404
* SSTEMR needs WORK of size 6*N, IWORK of size 3*N.
402405
* In addition, SLARRE needs WORK of size 6*N, IWORK of size 5*N.
@@ -519,15 +522,29 @@ SUBROUTINE CSTEMR( JOBZ, RANGE, N, D, E, VL, VU, IL, IU,
519522
ELSE IF( WANTZ.AND.(.NOT.ZQUERY) ) THEN
520523
CALL SLAEV2( D(1), E(1), D(2), R1, R2, CS, SN )
521524
END IF
525+
* D/S/LAE2 and D/S/LAEV2 outputs satisfy |R1| >= |R2|. However,
526+
* the following code requires R1 >= R2. Hence, we correct
527+
* the order of R1, R2, CS, SN if R1 < R2 before further processing.
528+
IF( R1.LT.R2 ) THEN
529+
E(2) = R1
530+
R1 = R2
531+
R2 = E(2)
532+
LAESWAP = .TRUE.
533+
ENDIF
522534
IF( ALLEIG.OR.
523535
$ (VALEIG.AND.(R2.GT.WL).AND.
524536
$ (R2.LE.WU)).OR.
525537
$ (INDEIG.AND.(IIL.EQ.1)) ) THEN
526538
M = M+1
527539
W( M ) = R2
528540
IF( WANTZ.AND.(.NOT.ZQUERY) ) THEN
529-
Z( 1, M ) = -SN
530-
Z( 2, M ) = CS
541+
IF( LAESWAP ) THEN
542+
Z( 1, M ) = CS
543+
Z( 2, M ) = SN
544+
ELSE
545+
Z( 1, M ) = -SN
546+
Z( 2, M ) = CS
547+
ENDIF
531548
* Note: At most one of SN and CS can be zero.
532549
IF (SN.NE.ZERO) THEN
533550
IF (CS.NE.ZERO) THEN
@@ -550,8 +567,13 @@ SUBROUTINE CSTEMR( JOBZ, RANGE, N, D, E, VL, VU, IL, IU,
550567
M = M+1
551568
W( M ) = R1
552569
IF( WANTZ.AND.(.NOT.ZQUERY) ) THEN
553-
Z( 1, M ) = CS
554-
Z( 2, M ) = SN
570+
IF( LAESWAP ) THEN
571+
Z( 1, M ) = -SN
572+
Z( 2, M ) = CS
573+
ELSE
574+
Z( 1, M ) = CS
575+
Z( 2, M ) = SN
576+
ENDIF
555577
* Note: At most one of SN and CS can be zero.
556578
IF (SN.NE.ZERO) THEN
557579
IF (CS.NE.ZERO) THEN

‎SRC/dstemr.f‎

Lines changed: 9 additions & 8 deletions
Original file line numberDiff line numberDiff line change
@@ -312,7 +312,8 @@
312312
*> Jim Demmel, University of California, Berkeley, USA \n
313313
*> Inderjit Dhillon, University of Texas, Austin, USA \n
314314
*> Osni Marques, LBNL/NERSC, USA \n
315-
*> Christof Voemel, University of California, Berkeley, USA
315+
*> Christof Voemel, University of California, Berkeley, USA \n
316+
*> Aravindh Krishnamoorthy, FAU, Erlangen, Germany \n
316317
*
317318
* =====================================================================
318319
SUBROUTINE DSTEMR( JOBZ, RANGE, N, D, E, VL, VU, IL, IU,
@@ -345,7 +346,7 @@ SUBROUTINE DSTEMR( JOBZ, RANGE, N, D, E, VL, VU, IL, IU,
345346
* ..
346347
* .. Local Scalars ..
347348
LOGICAL ALLEIG, INDEIG, LQUERY, VALEIG, WANTZ, ZQUERY,
348-
$ DLAESWAP
349+
$ LAESWAP
349350
INTEGER I, IBEGIN, IEND, IFIRST, IIL, IINDBL, IINDW,
350351
$ IINDWK, IINFO, IINSPL, IIU, ILAST, IN, INDD,
351352
$ INDE2, INDERR, INDGP, INDGRS, INDWRK, ITMP,
@@ -381,7 +382,7 @@ SUBROUTINE DSTEMR( JOBZ, RANGE, N, D, E, VL, VU, IL, IU,
381382
*
382383
LQUERY = ( ( LWORK.EQ.-1 ).OR.( LIWORK.EQ.-1 ) )
383384
ZQUERY = ( NZC.EQ.-1 )
384-
DLAESWAP = .FALSE.
385+
LAESWAP = .FALSE.
385386

386387
* DSTEMR needs WORK of size 6*N, IWORK of size 3*N.
387388
* In addition, DLARRE needs WORK of size 6*N, IWORK of size 5*N.
@@ -504,14 +505,14 @@ SUBROUTINE DSTEMR( JOBZ, RANGE, N, D, E, VL, VU, IL, IU,
504505
ELSE IF( WANTZ.AND.(.NOT.ZQUERY) ) THEN
505506
CALL DLAEV2( D(1), E(1), D(2), R1, R2, CS, SN )
506507
END IF
507-
* DLAE2 and DLAEV2 outputs satisfy |R1| >= |R2|. However,
508-
* the following DSTEMR requires R1 >= R2. Hence, we correct
508+
* D/S/LAE2 and D/S/LAEV2 outputs satisfy |R1| >= |R2|. However,
509+
* the following code requires R1 >= R2. Hence, we correct
509510
* the order of R1, R2, CS, SN if R1 < R2 before further processing.
510511
IF( R1.LT.R2 ) THEN
511512
E(2) = R1
512513
R1 = R2
513514
R2 = E(2)
514-
DLAESWAP = .TRUE.
515+
LAESWAP = .TRUE.
515516
ENDIF
516517
IF( ALLEIG.OR.
517518
$ (VALEIG.AND.(R2.GT.WL).AND.
@@ -520,7 +521,7 @@ SUBROUTINE DSTEMR( JOBZ, RANGE, N, D, E, VL, VU, IL, IU,
520521
M = M+1
521522
W( M ) = R2
522523
IF( WANTZ.AND.(.NOT.ZQUERY) ) THEN
523-
IF( DLAESWAP ) THEN
524+
IF( LAESWAP ) THEN
524525
Z( 1, M ) = CS
525526
Z( 2, M ) = SN
526527
ELSE
@@ -549,7 +550,7 @@ SUBROUTINE DSTEMR( JOBZ, RANGE, N, D, E, VL, VU, IL, IU,
549550
M = M+1
550551
W( M ) = R1
551552
IF( WANTZ.AND.(.NOT.ZQUERY) ) THEN
552-
IF( DLAESWAP ) THEN
553+
IF( LAESWAP ) THEN
553554
Z( 1, M ) = -SN
554555
Z( 2, M ) = CS
555556
ELSE

‎SRC/sstemr.f‎

Lines changed: 28 additions & 6 deletions
Original file line numberDiff line numberDiff line change
@@ -312,7 +312,8 @@
312312
*> Jim Demmel, University of California, Berkeley, USA \n
313313
*> Inderjit Dhillon, University of Texas, Austin, USA \n
314314
*> Osni Marques, LBNL/NERSC, USA \n
315-
*> Christof Voemel, University of California, Berkeley, USA
315+
*> Christof Voemel, University of California, Berkeley, USA \n
316+
*> Aravindh Krishnamoorthy, FAU, Erlangen, Germany \n
316317
*
317318
* =====================================================================
318319
SUBROUTINE SSTEMR( JOBZ, RANGE, N, D, E, VL, VU, IL, IU,
@@ -344,7 +345,8 @@ SUBROUTINE SSTEMR( JOBZ, RANGE, N, D, E, VL, VU, IL, IU,
344345
$ MINRGP = 3.0E-3 )
345346
* ..
346347
* .. Local Scalars ..
347-
LOGICAL ALLEIG, INDEIG, LQUERY, VALEIG, WANTZ, ZQUERY
348+
LOGICAL ALLEIG, INDEIG, LQUERY, VALEIG, WANTZ, ZQUERY,
349+
$ LAESWAP
348350
INTEGER I, IBEGIN, IEND, IFIRST, IIL, IINDBL, IINDW,
349351
$ IINDWK, IINFO, IINSPL, IIU, ILAST, IN, INDD,
350352
$ INDE2, INDERR, INDGP, INDGRS, INDWRK, ITMP,
@@ -378,6 +380,7 @@ SUBROUTINE SSTEMR( JOBZ, RANGE, N, D, E, VL, VU, IL, IU,
378380
*
379381
LQUERY = ( ( LWORK.EQ.-1 ).OR.( LIWORK.EQ.-1 ) )
380382
ZQUERY = ( NZC.EQ.-1 )
383+
LAESWAP = .FALSE.
381384

382385
* SSTEMR needs WORK of size 6*N, IWORK of size 3*N.
383386
* In addition, SLARRE needs WORK of size 6*N, IWORK of size 5*N.
@@ -500,15 +503,29 @@ SUBROUTINE SSTEMR( JOBZ, RANGE, N, D, E, VL, VU, IL, IU,
500503
ELSE IF( WANTZ.AND.(.NOT.ZQUERY) ) THEN
501504
CALL SLAEV2( D(1), E(1), D(2), R1, R2, CS, SN )
502505
END IF
506+
* D/S/LAE2 and D/S/LAEV2 outputs satisfy |R1| >= |R2|. However,
507+
* the following code requires R1 >= R2. Hence, we correct
508+
* the order of R1, R2, CS, SN if R1 < R2 before further processing.
509+
IF( R1.LT.R2 ) THEN
510+
E(2) = R1
511+
R1 = R2
512+
R2 = E(2)
513+
LAESWAP = .TRUE.
514+
ENDIF
503515
IF( ALLEIG.OR.
504516
$ (VALEIG.AND.(R2.GT.WL).AND.
505517
$ (R2.LE.WU)).OR.
506518
$ (INDEIG.AND.(IIL.EQ.1)) ) THEN
507519
M = M+1
508520
W( M ) = R2
509521
IF( WANTZ.AND.(.NOT.ZQUERY) ) THEN
510-
Z( 1, M ) = -SN
511-
Z( 2, M ) = CS
522+
IF( LAESWAP ) THEN
523+
Z( 1, M ) = CS
524+
Z( 2, M ) = SN
525+
ELSE
526+
Z( 1, M ) = -SN
527+
Z( 2, M ) = CS
528+
ENDIF
512529
* Note: At most one of SN and CS can be zero.
513530
IF (SN.NE.ZERO) THEN
514531
IF (CS.NE.ZERO) THEN
@@ -531,8 +548,13 @@ SUBROUTINE SSTEMR( JOBZ, RANGE, N, D, E, VL, VU, IL, IU,
531548
M = M+1
532549
W( M ) = R1
533550
IF( WANTZ.AND.(.NOT.ZQUERY) ) THEN
534-
Z( 1, M ) = CS
535-
Z( 2, M ) = SN
551+
IF( LAESWAP ) THEN
552+
Z( 1, M ) = -SN
553+
Z( 2, M ) = CS
554+
ELSE
555+
Z( 1, M ) = CS
556+
Z( 2, M ) = SN
557+
ENDIF
536558
* Note: At most one of SN and CS can be zero.
537559
IF (SN.NE.ZERO) THEN
538560
IF (CS.NE.ZERO) THEN

‎SRC/zstemr.f‎

Lines changed: 27 additions & 5 deletions
Original file line numberDiff line numberDiff line change
@@ -330,6 +330,7 @@
330330
*> Inderjit Dhillon, University of Texas, Austin, USA \n
331331
*> Osni Marques, LBNL/NERSC, USA \n
332332
*> Christof Voemel, University of California, Berkeley, USA \n
333+
*> Aravindh Krishnamoorthy, FAU, Erlangen, Germany \n
333334
*
334335
* =====================================================================
335336
SUBROUTINE ZSTEMR( JOBZ, RANGE, N, D, E, VL, VU, IL, IU,
@@ -361,7 +362,8 @@ SUBROUTINE ZSTEMR( JOBZ, RANGE, N, D, E, VL, VU, IL, IU,
361362
$ MINRGP = 1.0D-3 )
362363
* ..
363364
* .. Local Scalars ..
364-
LOGICAL ALLEIG, INDEIG, LQUERY, VALEIG, WANTZ, ZQUERY
365+
LOGICAL ALLEIG, INDEIG, LQUERY, VALEIG, WANTZ, ZQUERY,
366+
$ LAESWAP
365367
INTEGER I, IBEGIN, IEND, IFIRST, IIL, IINDBL, IINDW,
366368
$ IINDWK, IINFO, IINSPL, IIU, ILAST, IN, INDD,
367369
$ INDE2, INDERR, INDGP, INDGRS, INDWRK, ITMP,
@@ -397,6 +399,7 @@ SUBROUTINE ZSTEMR( JOBZ, RANGE, N, D, E, VL, VU, IL, IU,
397399
*
398400
LQUERY = ( ( LWORK.EQ.-1 ).OR.( LIWORK.EQ.-1 ) )
399401
ZQUERY = ( NZC.EQ.-1 )
402+
LAESWAP = .FALSE.
400403

401404
* DSTEMR needs WORK of size 6*N, IWORK of size 3*N.
402405
* In addition, DLARRE needs WORK of size 6*N, IWORK of size 5*N.
@@ -519,15 +522,29 @@ SUBROUTINE ZSTEMR( JOBZ, RANGE, N, D, E, VL, VU, IL, IU,
519522
ELSE IF( WANTZ.AND.(.NOT.ZQUERY) ) THEN
520523
CALL DLAEV2( D(1), E(1), D(2), R1, R2, CS, SN )
521524
END IF
525+
* D/S/LAE2 and D/S/LAEV2 outputs satisfy |R1| >= |R2|. However,
526+
* the following code requires R1 >= R2. Hence, we correct
527+
* the order of R1, R2, CS, SN if R1 < R2 before further processing.
528+
IF( R1.LT.R2 ) THEN
529+
E(2) = R1
530+
R1 = R2
531+
R2 = E(2)
532+
LAESWAP = .TRUE.
533+
ENDIF
522534
IF( ALLEIG.OR.
523535
$ (VALEIG.AND.(R2.GT.WL).AND.
524536
$ (R2.LE.WU)).OR.
525537
$ (INDEIG.AND.(IIL.EQ.1)) ) THEN
526538
M = M+1
527539
W( M ) = R2
528540
IF( WANTZ.AND.(.NOT.ZQUERY) ) THEN
529-
Z( 1, M ) = -SN
530-
Z( 2, M ) = CS
541+
IF( LAESWAP ) THEN
542+
Z( 1, M ) = CS
543+
Z( 2, M ) = SN
544+
ELSE
545+
Z( 1, M ) = -SN
546+
Z( 2, M ) = CS
547+
ENDIF
531548
* Note: At most one of SN and CS can be zero.
532549
IF (SN.NE.ZERO) THEN
533550
IF (CS.NE.ZERO) THEN
@@ -550,8 +567,13 @@ SUBROUTINE ZSTEMR( JOBZ, RANGE, N, D, E, VL, VU, IL, IU,
550567
M = M+1
551568
W( M ) = R1
552569
IF( WANTZ.AND.(.NOT.ZQUERY) ) THEN
553-
Z( 1, M ) = CS
554-
Z( 2, M ) = SN
570+
IF( LAESWAP ) THEN
571+
Z( 1, M ) = -SN
572+
Z( 2, M ) = CS
573+
ELSE
574+
Z( 1, M ) = CS
575+
Z( 2, M ) = SN
576+
ENDIF
555577
* Note: At most one of SN and CS can be zero.
556578
IF (SN.NE.ZERO) THEN
557579
IF (CS.NE.ZERO) THEN

0 commit comments

Comments
 (0)