Included bug fix provided by Zlatco on Jacobi SVD

Email from Zlatco on August 24th 2010:
The problem that was reported (with zero matrix) is caused by bad initialization to xLASSQ.
It should be ZERO, ONE and not ZERO, ZERO. In fact, I had it ZERO, ONE throughout
the complete development of the code and decided to change it to ZERO, ZERO a the very
end to make it "more elegant". That was stupid, because xLASSQ does not touch those
variables in case of zero vector, leaving scaling at ZERO, and in the nonzero case the scaling
is between ONE and SQRT(N). So, in case of zero vector, division by a variable that
 is normally bigger than ONE causes division by zero.
I have corrected that and few other things, stress tested the code and it should be OK now.


README:
i)   In xgejsv.f and xgesvj.f input parameters SCALE and
     SUMSQ in xlassq.f are now initially set as SCALE = ZERO, SUMSQ=ONE.
     Setting them both to zero (without carefully reading xlassq.f) caused 
     problems with exactly zero columns.
ii)  There was a problem in the branch that computes only SIGMA and U of a
     rank deficient matrix. The computed numerical rank (NR) was incorrectly
     written as N in parameter lists of the corresponding calls.
iii) In xgsvj0.f, xgsvj1.f testing the input parameters is changed to prevent
     unnecessarily negative INFO in some situations.
iv) Minor changes, renaming some variables etc.

1 files changed, 11 insertions, 10 deletions

diff --git a/SRC/sgejsv.f b/SRC/sgejsv.f
index 840f2e1a..a88202d2 100644
--- a/SRC/sgejsv.f
+++ b/SRC/sgejsv.f
@@ -488,7 +488,7 @@
       GOSCAL  = .TRUE.
       DO 1874 p = 1, N
          AAPP = ZERO
-         AAQQ = ZERO
+         AAQQ = ONE
          CALL SLASSQ( M, A(1,p), 1, AAPP, AAQQ )
          IF ( AAPP .GT. BIG ) THEN
             INFO = - 9
@@ -611,7 +611,7 @@
          IF ( L2TRAN ) THEN
             DO 1950 p = 1, M
                XSC   = ZERO
-               TEMP1 = ZERO
+               TEMP1 = ONE
                CALL SLASSQ( N, A(p,1), LDA, XSC, TEMP1 )
 *              SLASSQ gets both the ell_2 and the ell_infinity norm
 *              in one pass through the vector
@@ -642,7 +642,7 @@
       IF ( L2TRAN ) THEN
 *
          XSC   = ZERO
-         TEMP1 = ZERO
+         TEMP1 = ONE
          CALL SLASSQ( N, SVA, 1, XSC, TEMP1 )
          TEMP1 = ONE / TEMP1
 *
@@ -696,6 +696,7 @@
             KILL   = LSVEC
             LSVEC  = RSVEC
             RSVEC  = KILL
+            IF ( LSVEC ) N1 = N 
 *
             ROWPIV = .TRUE.
          END IF
@@ -1473,10 +1474,10 @@
 *           Assemble the left singular vector matrix U (M x N).
 *
             IF ( N .LT. M ) THEN
-               CALL SLASET( 'A',  M-N, N, ZERO, ZERO, U(NR+1,1), LDU )
+               CALL SLASET( 'A',  M-N, N, ZERO, ZERO, U(N+1,1), LDU )
                IF ( N .LT. N1 ) THEN
                   CALL SLASET( 'A',N,  N1-N, ZERO, ZERO,  U(1,N+1),LDU )
-                  CALL SLASET( 'A',M-N,N1-N, ZERO, ONE,U(NR+1,N+1),LDU )
+                  CALL SLASET( 'A',M-N,N1-N, ZERO, ONE,U(N+1,N+1),LDU )
                END IF
             END IF
             CALL SORMQR( 'Left', 'No Tr', M, N1, N, A, LDA, WORK, U,
@@ -1580,11 +1581,11 @@
 *           At this moment, V contains the right singular vectors of A.
 *           Next, assemble the left singular vector matrix U (M x N).
 *
-         IF ( N .LT. M ) THEN
-            CALL SLASET( 'A',  M-N, N, ZERO, ZERO, U(NR+1,1), LDU )
-            IF ( N .LT. N1 ) THEN
-               CALL SLASET( 'A',N,  N1-N, ZERO, ZERO,  U(1,N+1),LDU )
-               CALL SLASET( 'A',M-N,N1-N, ZERO, ONE,U(NR+1,N+1),LDU )
+         IF ( NR .LT. M ) THEN
+            CALL SLASET( 'A',  M-NR, NR, ZERO, ZERO, U(NR+1,1), LDU )
+            IF ( NR .LT. N1 ) THEN
+               CALL SLASET( 'A',NR,  N1-NR, ZERO, ZERO,  U(1,NR+1),LDU )
+               CALL SLASET( 'A',M-NR,N1-NR, ZERO, ONE,U(NR+1,NR+1),LDU )
             END IF
          END IF
 *