Bug in stgsna.f / dtgsna.f, which moves to incorrect SVD calculations.

[elivanova]

$ 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 

[elivanova]

Thank you very much. 22.12.2016 13:20, langou пишет:

…

Nevermind. I can edit your post for you. I just did it. Sorry for the noise. — You are receiving this because you authored the thread. Reply to this email directly, view it on GitHub <#103 (comment)>, or mute the thread <https://github.com/notifications/unsubscribe-auth/AT2ZSnqEHquyduczUEEXo9r_ui1PoSa9ks5rKk77gaJpZM4LT07X>.