Bug in stgsna.f / dtgsna.f, which moves to incorrect SVD calculations. #103
Description
$ cat t.f
PROGRAM BUG_STGSNA
IMPLICIT NONE
INTEGER INFO, LWORK, M, MM, IWORK(20), i, j
REAL A(2,2), B(2,2), RWORK( 40 )
REAL VL(2,4), VR(2,4), S(2), DIF(2)
COMPLEX C( 2, 2 ), U( 2, 2 ), VT( 2, 2 )
COMPLEX WORK( 20 )
LOGICAL SELECT(2)
C(1,1) = (1.0,1.0)
C(2,1) = 1.0
C(1,2) = (-1.0,1.0)
C(2,2) = - 1.0
LWORK = 40
CALL CGESVD('A','A',2,2,C,2,S,U,2,VT,2,WORK,LWORK,RWORK,INFO)
SELECT(1) = .TRUE.
SELECT(2) = .TRUE.
A(1,1) = 1.0
A(1,2) = 1.0
A(2,1) = -A(1,2)
A(2,2) = A(1,1)
B(1,1) = 1.0
B(1,2) = 0.0
B(2,1) = 0.0
B(2,2) = 1.0
PRINT 10
PRINT 11,((A(i,j),j=1,2),i=1,2)
PRINT 12
PRINT 11,((B(i,j),j=1,2),i=1,2)
PRINT 13, S(2)
MM = 2
CALL STGEVC('B','S',SELECT,2,A,2,B,2,VL,2,VR,2,
& MM,M,RWORK,IWORK,INFO)
CALL STGSNA('B','S',SELECT,2,A,2,B,2,VL,2,
&VR,2,S,DIF,MM,M,RWORK,LWORK,IWORK,INFO)
PRINT 14, DIF(2)
10 FORMAT(1X,'INPUT matrix A')
11 FORMAT(1X,2F5.0)
12 FORMAT(1X,'INPUT matrix B')
13 FORMAT(1X,'Correct DIF(after CGESVD) ',F5.2)
14 FORMAT(1X,'DIF after STGSNA ',F5.2)
END
The result is:
INPUT matrix A
1. 1.
-1. 1.
INPUT matrix B
1. 0.
0. 1.
Correct DIF(after CGESVD) 0.87
DIF after STGSNA 0.62
Instead of the correct result:
INPUT matrix A
1. 1.
-1. 1.
INPUT matrix B
1. 0.
0. 1.
Correct DIF(after CGESVD) 0.87
DIF after STGSNA 0.87
Solution: correct in stgsna.f and dtgsna.f: (lines 637-639).
Instead of
637 ROOT1 = C1 + SQRT( C1*C1-4.0*C2 )
638 ROOT2 = C2 / ROOT1
639 ROOT1 = ROOT1 / TWO
Should be:
637 ROOT1 = C1 + SQRT( C1*C1-FOUR*C2 )
638 ROOT1 = ROOT1 / TWO
639 ROOT2 = C2 / ROOT1
Metadata
Metadata
Assignees
Labels
No labels