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.

diff --git a/SRC/dgejsv.f b/SRC/dgejsv.f
index 213e85f..2dd98e1 100644 (file)
--- a/SRC/dgejsv.f
+++ b/SRC/dgejsv.f
@@ -489,7 +489,7 @@
       GOSCAL  = .TRUE.
       DO 1874 p = 1, N
          AAPP = ZERO
-         AAQQ = ZERO
+         AAQQ = ONE
          CALL DLASSQ( M, A(1,p), 1, AAPP, AAQQ )
          IF ( AAPP .GT. BIG ) THEN
             INFO = - 9
@@ -612,7 +612,7 @@
          IF ( L2TRAN ) THEN
             DO 1950 p = 1, M
                XSC   = ZERO
-               TEMP1 = ZERO
+               TEMP1 = ONE
                CALL DLASSQ( N, A(p,1), LDA, XSC, TEMP1 )
 *              DLASSQ gets both the ell_2 and the ell_infinity norm
 *              in one pass through the vector
@@ -643,7 +643,7 @@
       IF ( L2TRAN ) THEN
 *
          XSC   = ZERO
-         TEMP1 = ZERO
+         TEMP1 = ONE
          CALL DLASSQ( N, SVA, 1, XSC, TEMP1 )
          TEMP1 = ONE / TEMP1
 *
@@ -697,6 +697,7 @@
             KILL   = LSVEC
             LSVEC  = RSVEC
             RSVEC  = KILL
+            IF ( LSVEC ) N1 = N
 *
             ROWPIV = .TRUE.
          END IF
@@ -1475,10 +1476,10 @@
 *           Assemble the left singular vector matrix U (M x N).
 *
             IF ( N .LT. M ) THEN
-               CALL DLASET( 'A',  M-N, N, ZERO, ZERO, U(NR+1,1), LDU )
+               CALL DLASET( 'A',  M-N, N, ZERO, ZERO, U(N+1,1), LDU )
                IF ( N .LT. N1 ) THEN
                   CALL DLASET( 'A',N,  N1-N, ZERO, ZERO,  U(1,N+1),LDU )
-                  CALL DLASET( 'A',M-N,N1-N, ZERO, ONE,U(NR+1,N+1),LDU )
+                  CALL DLASET( 'A',M-N,N1-N, ZERO, ONE,U(N+1,N+1),LDU )
                END IF
             END IF
             CALL DORMQR( 'Left', 'No Tr', M, N1, N, A, LDA, WORK, U,
@@ -1582,11 +1583,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 DLASET( 'A',  M-N, N, ZERO, ZERO, U(NR+1,1), LDU )
-            IF ( N .LT. N1 ) THEN
-               CALL DLASET( 'A',N,  N1-N, ZERO, ZERO,  U(1,N+1),LDU )
-               CALL DLASET( 'A',M-N,N1-N, ZERO, ONE,U(NR+1,N+1),LDU )
+         IF ( NR .LT. M ) THEN
+            CALL DLASET( 'A',  M-NR, NR, ZERO, ZERO, U(NR+1,1), LDU )
+            IF ( NR .LT. N1 ) THEN
+               CALL DLASET( 'A',NR,  N1-NR, ZERO, ZERO,  U(1,NR+1),LDU )
+               CALL DLASET( 'A',M-NR,N1-NR, ZERO, ONE,U(NR+1,NR+1),LDU )
             END IF
          END IF
 *

diff --git a/SRC/dgesvj.f b/SRC/dgesvj.f
index b9d568e..f3264b1 100644 (file)
--- a/SRC/dgesvj.f
+++ b/SRC/dgesvj.f
@@ -93,7 +93,7 @@
 *  Arguments
 *  =========
 *
-*  JOBA    (input) CHARACTER*1
+*  JOBA    (input) CHARACTER* 1
 *          Specifies the structure of A.
 *          = 'L': The input matrix A is lower triangular;
 *          = 'U': The input matrix A is upper triangular;
@@ -222,7 +222,7 @@
 *                    orthogonal up to CTOL*EPS, EPS=DLAMCH('E').
 *                    It is required that CTOL >= ONE, i.e. it is not
 *                    allowed to force the routine to obtain orthogonality
-*                    below EPSILON.
+*                    below EPS.
 *          On exit :
 *          WORK(1) = SCALE is the scaling factor such that SCALE*SVA(1:N)
 *                    are the computed singular values of A.
@@ -262,9 +262,9 @@
 *     ..
 *     .. Local Scalars ..
       DOUBLE PRECISION   AAPP, AAPP0, AAPQ, AAQQ, APOAQ, AQOAP, BIG,
-     +                   BIGTHETA, CS, CTOL, EPSILON, LARGE, MXAAPQ,
+     +                   BIGTHETA, CS, CTOL, EPSLN, LARGE, MXAAPQ,
      +                   MXSINJ, ROOTBIG, ROOTEPS, ROOTSFMIN, ROOTTOL,
-     +                   SCALE, SFMIN, SMALL, SN, T, TEMP1, THETA,
+     +                   SKL, SFMIN, SMALL, SN, T, TEMP1, THETA,
      +                   THSIGN, TOL
       INTEGER            BLSKIP, EMPTSW, i, ibr, IERR, igl, IJBLSK, ir1,
      +                   ISWROT, jbc, jgl, KBL, LKAHEAD, MVL, N2, N34,
@@ -368,21 +368,21 @@
 *     ... and the machine dependent parameters are
 *[!]  (Make sure that DLAMCH() works properly on the target machine.)
 *
-      EPSILON = DLAMCH( 'Epsilon' )
-      ROOTEPS = DSQRT( EPSILON )
+      EPSLN = DLAMCH( 'Epsilon' )
+      ROOTEPS = DSQRT( EPSLN )
       SFMIN = DLAMCH( 'SafeMinimum' )
       ROOTSFMIN = DSQRT( SFMIN )
-      SMALL = SFMIN / EPSILON
+      SMALL = SFMIN / EPSLN
       BIG = DLAMCH( 'Overflow' )
 *     BIG         = ONE    / SFMIN
       ROOTBIG = ONE / ROOTSFMIN
       LARGE = BIG / DSQRT( DBLE( M*N ) )
       BIGTHETA = ONE / ROOTEPS
 *
-      TOL = CTOL*EPSILON
+      TOL = CTOL*EPSLN
       ROOTTOL = DSQRT( TOL )
 *
-      IF( DBLE( M )*EPSILON.GE.ONE ) THEN
+      IF( DBLE( M )*EPSLN.GE.ONE ) THEN
          INFO = -5
          CALL XERBLA( 'DGESVJ', -INFO )
          RETURN
@@ -407,7 +407,7 @@
 *     DSQRT(N)*max_i SVA(i) does not overflow. If INFinite entries
 *     in A are detected, the procedure returns with INFO=-6.
 *
-      SCALE = ONE / DSQRT( DBLE( M )*DBLE( N ) )
+      SKL= ONE / DSQRT( DBLE( M )*DBLE( N ) )
       NOSCALE = .TRUE.
       GOSCALE = .TRUE.
 *
@@ -415,7 +415,7 @@
 *        the input matrix is M-by-N lower triangular (trapezoidal)
          DO 1874 p = 1, N
             AAPP = ZERO
-            AAQQ = ZERO
+            AAQQ = ONE
             CALL DLASSQ( M-p+1, A( p, p ), 1, AAPP, AAQQ )
             IF( AAPP.GT.BIG ) THEN
                INFO = -6
@@ -427,11 +427,11 @@
                SVA( p ) = AAPP*AAQQ
             ELSE
                NOSCALE = .FALSE.
-               SVA( p ) = AAPP*( AAQQ*SCALE )
+               SVA( p ) = AAPP*( AAQQ*SKL)
                IF( GOSCALE ) THEN
                   GOSCALE = .FALSE.
                   DO 1873 q = 1, p - 1
-                     SVA( q ) = SVA( q )*SCALE
+                     SVA( q ) = SVA( q )*SKL
  1873             CONTINUE
                END IF
             END IF
@@ -440,7 +440,7 @@
 *        the input matrix is M-by-N upper triangular (trapezoidal)
          DO 2874 p = 1, N
             AAPP = ZERO
-            AAQQ = ZERO
+            AAQQ = ONE
             CALL DLASSQ( p, A( 1, p ), 1, AAPP, AAQQ )
             IF( AAPP.GT.BIG ) THEN
                INFO = -6
@@ -452,11 +452,11 @@
                SVA( p ) = AAPP*AAQQ
             ELSE
                NOSCALE = .FALSE.
-               SVA( p ) = AAPP*( AAQQ*SCALE )
+               SVA( p ) = AAPP*( AAQQ*SKL)
                IF( GOSCALE ) THEN
                   GOSCALE = .FALSE.
                   DO 2873 q = 1, p - 1
-                     SVA( q ) = SVA( q )*SCALE
+                     SVA( q ) = SVA( q )*SKL
  2873             CONTINUE
                END IF
             END IF
@@ -465,7 +465,7 @@
 *        the input matrix is M-by-N general dense
          DO 3874 p = 1, N
             AAPP = ZERO
-            AAQQ = ZERO
+            AAQQ = ONE
             CALL DLASSQ( M, A( 1, p ), 1, AAPP, AAQQ )
             IF( AAPP.GT.BIG ) THEN
                INFO = -6
@@ -477,18 +477,18 @@
                SVA( p ) = AAPP*AAQQ
             ELSE
                NOSCALE = .FALSE.
-               SVA( p ) = AAPP*( AAQQ*SCALE )
+               SVA( p ) = AAPP*( AAQQ*SKL)
                IF( GOSCALE ) THEN
                   GOSCALE = .FALSE.
                   DO 3873 q = 1, p - 1
-                     SVA( q ) = SVA( q )*SCALE
+                     SVA( q ) = SVA( q )*SKL
  3873             CONTINUE
                END IF
             END IF
  3874    CONTINUE
       END IF
 *
-      IF( NOSCALE )SCALE = ONE
+      IF( NOSCALE )SKL= ONE
 *
 *     Move the smaller part of the spectrum from the underflow threshold
 *(!)  Start by determining the position of the nonzero entries of the
@@ -517,9 +517,9 @@
 * #:) Quick return for one-column matrix
 *
       IF( N.EQ.1 ) THEN
-         IF( LSVEC )CALL DLASCL( 'G', 0, 0, SVA( 1 ), SCALE, M, 1,
+         IF( LSVEC )CALL DLASCL( 'G', 0, 0, SVA( 1 ), SKL, M, 1,
      +                           A( 1, 1 ), LDA, IERR )
-         WORK( 1 ) = ONE / SCALE
+         WORK( 1 ) = ONE / SKL
          IF( SVA( 1 ).GE.SFMIN ) THEN
             WORK( 2 ) = ONE
          ELSE
@@ -535,7 +535,7 @@
 *     Protect small singular values from underflow, and try to
 *     avoid underflows/overflows in computing Jacobi rotations.
 *
-      SN = DSQRT( SFMIN / EPSILON )
+      SN = DSQRT( SFMIN / EPSLN )
       TEMP1 = DSQRT( BIG / DBLE( N ) )
       IF( ( AAPP.LE.SN ) .OR. ( AAQQ.GE.TEMP1 ) .OR.
      +    ( ( SN.LE.AAQQ ) .AND. ( AAPP.LE.TEMP1 ) ) ) THEN
@@ -563,10 +563,10 @@
       IF( TEMP1.NE.ONE ) THEN
          CALL DLASCL( 'G', 0, 0, ONE, TEMP1, N, 1, SVA, N, IERR )
       END IF
-      SCALE = TEMP1*SCALE
-      IF( SCALE.NE.ONE ) THEN
-         CALL DLASCL( JOBA, 0, 0, ONE, SCALE, M, N, A, LDA, IERR )
-         SCALE = ONE / SCALE
+      SKL= TEMP1*SKL
+      IF( SKL.NE.ONE ) THEN
+         CALL DLASCL( JOBA, 0, 0, ONE, SKL, M, N, A, LDA, IERR )
+         SKL= ONE / SKL
       END IF
 *
 *     Row-cyclic Jacobi SVD algorithm with column pivoting
@@ -639,30 +639,30 @@
 *
             CALL DGSVJ0( JOBV, M-N34, N-N34, A( N34+1, N34+1 ), LDA,
      +                   WORK( N34+1 ), SVA( N34+1 ), MVL,
-     +                   V( N34*q+1, N34+1 ), LDV, EPSILON, SFMIN, TOL,
+     +                   V( N34*q+1, N34+1 ), LDV, EPSLN, SFMIN, TOL,
      +                   2, WORK( N+1 ), LWORK-N, IERR )
 *
             CALL DGSVJ0( JOBV, M-N2, N34-N2, A( N2+1, N2+1 ), LDA,
      +                   WORK( N2+1 ), SVA( N2+1 ), MVL,
-     +                   V( N2*q+1, N2+1 ), LDV, EPSILON, SFMIN, TOL, 2,
+     +                   V( N2*q+1, N2+1 ), LDV, EPSLN, SFMIN, TOL, 2,
      +                   WORK( N+1 ), LWORK-N, IERR )
 *
             CALL DGSVJ1( JOBV, M-N2, N-N2, N4, A( N2+1, N2+1 ), LDA,
      +                   WORK( N2+1 ), SVA( N2+1 ), MVL,
-     +                   V( N2*q+1, N2+1 ), LDV, EPSILON, SFMIN, TOL, 1,
+     +                   V( N2*q+1, N2+1 ), LDV, EPSLN, SFMIN, TOL, 1,
      +                   WORK( N+1 ), LWORK-N, IERR )
 *
             CALL DGSVJ0( JOBV, M-N4, N2-N4, A( N4+1, N4+1 ), LDA,
      +                   WORK( N4+1 ), SVA( N4+1 ), MVL,
-     +                   V( N4*q+1, N4+1 ), LDV, EPSILON, SFMIN, TOL, 1,
+     +                   V( N4*q+1, N4+1 ), LDV, EPSLN, SFMIN, TOL, 1,
      +                   WORK( N+1 ), LWORK-N, IERR )
 *
             CALL DGSVJ0( JOBV, M, N4, A, LDA, WORK, SVA, MVL, V, LDV,
-     +                   EPSILON, SFMIN, TOL, 1, WORK( N+1 ), LWORK-N,
+     +                   EPSLN, SFMIN, TOL, 1, WORK( N+1 ), LWORK-N,
      +                   IERR )
 *
             CALL DGSVJ1( JOBV, M, N2, N4, A, LDA, WORK, SVA, MVL, V,
-     +                   LDV, EPSILON, SFMIN, TOL, 1, WORK( N+1 ),
+     +                   LDV, EPSLN, SFMIN, TOL, 1, WORK( N+1 ),
      +                   LWORK-N, IERR )
 *
 *
@@ -670,21 +670,21 @@
 *
 *
             CALL DGSVJ0( JOBV, N4, N4, A, LDA, WORK, SVA, MVL, V, LDV,
-     +                   EPSILON, SFMIN, TOL, 2, WORK( N+1 ), LWORK-N,
+     +                   EPSLN, SFMIN, TOL, 2, WORK( N+1 ), LWORK-N,
      +                   IERR )
 *
             CALL DGSVJ0( JOBV, N2, N4, A( 1, N4+1 ), LDA, WORK( N4+1 ),
      +                   SVA( N4+1 ), MVL, V( N4*q+1, N4+1 ), LDV,
-     +                   EPSILON, SFMIN, TOL, 1, WORK( N+1 ), LWORK-N,
+     +                   EPSLN, SFMIN, TOL, 1, WORK( N+1 ), LWORK-N,
      +                   IERR )
 *
             CALL DGSVJ1( JOBV, N2, N2, N4, A, LDA, WORK, SVA, MVL, V,
-     +                   LDV, EPSILON, SFMIN, TOL, 1, WORK( N+1 ),
+     +                   LDV, EPSLN, SFMIN, TOL, 1, WORK( N+1 ),
      +                   LWORK-N, IERR )
 *
             CALL DGSVJ0( JOBV, N2+N4, N4, A( 1, N2+1 ), LDA,
      +                   WORK( N2+1 ), SVA( N2+1 ), MVL,
-     +                   V( N2*q+1, N2+1 ), LDV, EPSILON, SFMIN, TOL, 1,
+     +                   V( N2*q+1, N2+1 ), LDV, EPSLN, SFMIN, TOL, 1,
      +                   WORK( N+1 ), LWORK-N, IERR )
 
          END IF
@@ -753,7 +753,7 @@
                         SVA( p ) = DNRM2( M, A( 1, p ), 1 )*WORK( p )
                      ELSE
                         TEMP1 = ZERO
-                        AAPP = ZERO
+                        AAPP = ONE
                         CALL DLASSQ( M, A( 1, p ), 1, TEMP1, AAPP )
                         SVA( p ) = TEMP1*DSQRT( AAPP )*WORK( p )
                      END IF
@@ -841,8 +841,8 @@
      +                                              FASTR )
                                     SVA( q ) = AAQQ*DSQRT( DMAX1( ZERO,
      +                                         ONE+T*APOAQ*AAPQ ) )
-                                    AAPP = AAPP*DSQRT( ONE-T*AQOAP*
-     +                                     AAPQ )
+                                    AAPP = AAPP*DSQRT( DMAX1( ZERO,
+     +                                     ONE-T*AQOAP*AAPQ ) )
                                     MXSINJ = DMAX1( MXSINJ, DABS( T ) )
 *
                                  ELSE
@@ -989,7 +989,7 @@
      +                                         WORK( q )
                                  ELSE
                                     T = ZERO
-                                    AAQQ = ZERO
+                                    AAQQ = ONE
                                     CALL DLASSQ( M, A( 1, q ), 1, T,
      +                                           AAQQ )
                                     SVA( q ) = T*DSQRT( AAQQ )*WORK( q )
@@ -1002,7 +1002,7 @@
      +                                     WORK( p )
                                  ELSE
                                     T = ZERO
-                                    AAPP = ZERO
+                                    AAPP = ONE
                                     CALL DLASSQ( M, A( 1, p ), 1, T,
      +                                           AAPP )
                                     AAPP = T*DSQRT( AAPP )*WORK( p )
@@ -1164,8 +1164,8 @@
                                     MXSINJ = DMAX1( MXSINJ, DABS( SN ) )
                                     SVA( q ) = AAQQ*DSQRT( DMAX1( ZERO,
      +                                         ONE+T*APOAQ*AAPQ ) )
-                                    AAPP = AAPP*DSQRT( ONE-T*AQOAP*
-     +                                     AAPQ )
+                                    AAPP = AAPP*DSQRT( DMAX1( ZERO, 
+     +                                     ONE-T*AQOAP*AAPQ ) )
 *
                                     APOAQ = WORK( p ) / WORK( q )
                                     AQOAP = WORK( q ) / WORK( p )
@@ -1316,7 +1316,7 @@
      +                                         WORK( q )
                                  ELSE
                                     T = ZERO
-                                    AAQQ = ZERO
+                                    AAQQ = ONE
                                     CALL DLASSQ( M, A( 1, q ), 1, T,
      +                                           AAQQ )
                                     SVA( q ) = T*DSQRT( AAQQ )*WORK( q )
@@ -1329,7 +1329,7 @@
      +                                     WORK( p )
                                  ELSE
                                     T = ZERO
-                                    AAPP = ZERO
+                                    AAPP = ONE
                                     CALL DLASSQ( M, A( 1, p ), 1, T,
      +                                           AAPP )
                                     AAPP = T*DSQRT( AAPP )*WORK( p )
@@ -1395,7 +1395,7 @@
             SVA( N ) = DNRM2( M, A( 1, N ), 1 )*WORK( N )
          ELSE
             T = ZERO
-            AAPP = ZERO
+            AAPP = ONE
             CALL DLASSQ( M, A( 1, N ), 1, T, AAPP )
             SVA( N ) = T*DSQRT( AAPP )*WORK( N )
          END IF
@@ -1446,12 +1446,12 @@
          END IF
          IF( SVA( p ).NE.ZERO ) THEN
             N4 = N4 + 1
-            IF( SVA( p )*SCALE.GT.SFMIN )N2 = N2 + 1
+            IF( SVA( p )*SKL.GT.SFMIN )N2 = N2 + 1
          END IF
  5991 CONTINUE
       IF( SVA( N ).NE.ZERO ) THEN
          N4 = N4 + 1
-         IF( SVA( N )*SCALE.GT.SFMIN )N2 = N2 + 1
+         IF( SVA( N )*SKL.GT.SFMIN )N2 = N2 + 1
       END IF
 *
 *     Normalize the left singular vectors.
@@ -1478,17 +1478,17 @@
       END IF
 *
 *     Undo scaling, if necessary (and possible).
-      IF( ( ( SCALE.GT.ONE ) .AND. ( SVA( 1 ).LT.( BIG /
-     +    SCALE ) ) ) .OR. ( ( SCALE.LT.ONE ) .AND. ( SVA( N2 ).GT.
-     +    ( SFMIN / SCALE ) ) ) ) THEN
+      IF( ( ( SKL.GT.ONE ) .AND. ( SVA( 1 ).LT.( BIG /
+     +    SKL) ) ) .OR. ( ( SKL.LT.ONE ) .AND. ( SVA( N2 ).GT.
+     +    ( SFMIN / SKL) ) ) ) THEN
          DO 2400 p = 1, N
-            SVA( p ) = SCALE*SVA( p )
+            SVA( p ) = SKL*SVA( p )
  2400    CONTINUE
-         SCALE = ONE
+         SKL= ONE
       END IF
 *
-      WORK( 1 ) = SCALE
-*     The singular values of A are SCALE*SVA(1:N). If SCALE.NE.ONE
+      WORK( 1 ) = SKL
+*     The singular values of A are SKL*SVA(1:N). If SKL.NE.ONE
 *     then some of the singular values may overflow or underflow and
 *     the spectrum is given in this factored representation.
 *

diff --git a/SRC/dgsvj0.f b/SRC/dgsvj0.f
index 2df0181..0312690 100644 (file)
--- a/SRC/dgsvj0.f
+++ b/SRC/dgsvj0.f
@@ -183,9 +183,10 @@
          INFO = -3
       ELSE IF( LDA.LT.M ) THEN
          INFO = -5
-      ELSE IF( MV.LT.0 ) THEN
+      ELSE IF( ( RSVEC.OR.APPLV ) .AND. ( MV.LT.0 ) ) THEN
          INFO = -8
-      ELSE IF( LDV.LT.M ) THEN
+      ELSE IF( ( RSVEC.AND.( LDV.LT.N ) ).OR. 
+     &         ( APPLV.AND.( LDV.LT.MV ) ) ) THEN
          INFO = -10
       ELSE IF( TOL.LE.EPS ) THEN
          INFO = -13
@@ -307,7 +308,7 @@
                         SVA( p ) = DNRM2( M, A( 1, p ), 1 )*D( p )
                      ELSE
                         TEMP1 = ZERO
-                        AAPP = ZERO
+                        AAPP = ONE
                         CALL DLASSQ( M, A( 1, p ), 1, TEMP1, AAPP )
                         SVA( p ) = TEMP1*DSQRT( AAPP )*D( p )
                      END IF
@@ -391,8 +392,8 @@
      +                                              FASTR )
                                     SVA( q ) = AAQQ*DSQRT( DMAX1( ZERO,
      +                                         ONE+T*APOAQ*AAPQ ) )
-                                    AAPP = AAPP*DSQRT( ONE-T*AQOAP*
-     +                                     AAPQ )
+                                    AAPP = AAPP*DSQRT( DMAX1( ZERO, 
+     +                                     ONE-T*AQOAP*AAPQ ) )
                                     MXSINJ = DMAX1( MXSINJ, DABS( T ) )
 *
                                  ELSE
@@ -535,7 +536,7 @@
      +                                         D( q )
                                  ELSE
                                     T = ZERO
-                                    AAQQ = ZERO
+                                    AAQQ = ONE
                                     CALL DLASSQ( M, A( 1, q ), 1, T,
      +                                           AAQQ )
                                     SVA( q ) = T*DSQRT( AAQQ )*D( q )
@@ -548,7 +549,7 @@
      +                                     D( p )
                                  ELSE
                                     T = ZERO
-                                    AAPP = ZERO
+                                    AAPP = ONE
                                     CALL DLASSQ( M, A( 1, p ), 1, T,
      +                                           AAPP )
                                     AAPP = T*DSQRT( AAPP )*D( p )
@@ -707,8 +708,8 @@
                                     MXSINJ = DMAX1( MXSINJ, DABS( SN ) )
                                     SVA( q ) = AAQQ*DSQRT( DMAX1( ZERO,
      +                                         ONE+T*APOAQ*AAPQ ) )
-                                    AAPP = AAPP*DSQRT( ONE-T*AQOAP*
-     +                                     AAPQ )
+                                    AAPP = AAPP*DSQRT( DMAX1( ZERO, 
+     +                                     ONE-T*AQOAP*AAPQ ) )
 *
                                     APOAQ = D( p ) / D( q )
                                     AQOAP = D( q ) / D( p )
@@ -856,7 +857,7 @@
      +                                         D( q )
                                  ELSE
                                     T = ZERO
-                                    AAQQ = ZERO
+                                    AAQQ = ONE
                                     CALL DLASSQ( M, A( 1, q ), 1, T,
      +                                           AAQQ )
                                     SVA( q ) = T*DSQRT( AAQQ )*D( q )
@@ -869,7 +870,7 @@
      +                                     D( p )
                                  ELSE
                                     T = ZERO
-                                    AAPP = ZERO
+                                    AAPP = ONE
                                     CALL DLASSQ( M, A( 1, p ), 1, T,
      +                                           AAPP )
                                     AAPP = T*DSQRT( AAPP )*D( p )
@@ -932,7 +933,7 @@
             SVA( N ) = DNRM2( M, A( 1, N ), 1 )*D( N )
          ELSE
             T = ZERO
-            AAPP = ZERO
+            AAPP = ONE
             CALL DLASSQ( M, A( 1, N ), 1, T, AAPP )
             SVA( N ) = T*DSQRT( AAPP )*D( N )
          END IF

diff --git a/SRC/dgsvj1.f b/SRC/dgsvj1.f
index b8794a2..05f163c 100644 (file)
--- a/SRC/dgsvj1.f
+++ b/SRC/dgsvj1.f
@@ -205,9 +205,10 @@
          INFO = -4
       ELSE IF( LDA.LT.M ) THEN
          INFO = -6
-      ELSE IF( MV.LT.0 ) THEN
+      ELSE IF( ( RSVEC.OR.APPLV ) .AND. ( MV.LT.0 ) ) THEN
          INFO = -9
-      ELSE IF( LDV.LT.M ) THEN
+      ELSE IF( ( RSVEC.AND.( LDV.LT.N ) ).OR. 
+     &         ( APPLV.AND.( LDV.LT.MV ) )  ) THEN
          INFO = -11
       ELSE IF( TOL.LE.EPS ) THEN
          INFO = -14
@@ -406,8 +407,8 @@
                                     MXSINJ = DMAX1( MXSINJ, DABS( SN ) )
                                     SVA( q ) = AAQQ*DSQRT( DMAX1( ZERO,
      +                                         ONE+T*APOAQ*AAPQ ) )
-                                    AAPP = AAPP*DSQRT( ONE-T*AQOAP*
-     +                                     AAPQ )
+                                    AAPP = AAPP*DSQRT( DMAX1( ZERO, 
+     +                                    ONE-T*AQOAP*AAPQ ) )
 
                                     APOAQ = D( p ) / D( q )
                                     AQOAP = D( q ) / D( p )
@@ -555,7 +556,7 @@
      +                                         D( q )
                                  ELSE
                                     T = ZERO
-                                    AAQQ = ZERO
+                                    AAQQ = ONE
                                     CALL DLASSQ( M, A( 1, q ), 1, T,
      +                                           AAQQ )
                                     SVA( q ) = T*DSQRT( AAQQ )*D( q )
@@ -568,7 +569,7 @@
      +                                     D( p )
                                  ELSE
                                     T = ZERO
-                                    AAPP = ZERO
+                                    AAPP = ONE
                                     CALL DLASSQ( M, A( 1, p ), 1, T,
      +                                           AAPP )
                                     AAPP = T*DSQRT( AAPP )*D( p )
@@ -635,7 +636,7 @@
             SVA( N ) = DNRM2( M, A( 1, N ), 1 )*D( N )
          ELSE
             T = ZERO
-            AAPP = ZERO
+            AAPP = ONE
             CALL DLASSQ( M, A( 1, N ), 1, T, AAPP )
             SVA( N ) = T*DSQRT( AAPP )*D( N )
          END IF

diff --git a/SRC/sgejsv.f b/SRC/sgejsv.f
index 840f2e1..a88202d 100644 (file)
--- 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
 *

diff --git a/SRC/sgesvj.f b/SRC/sgesvj.f
index f6939cc..e38f2f6 100644 (file)
--- a/SRC/sgesvj.f
+++ b/SRC/sgesvj.f
@@ -93,7 +93,7 @@
 *  Arguments
 *  =========
 *
-*  JOBA    (input) CHARACTER*1
+*  JOBA    (input) CHARACTER* 1
 *          Specifies the structure of A.
 *          = 'L': The input matrix A is lower triangular;
 *          = 'U': The input matrix A is upper triangular;
@@ -260,9 +260,9 @@
 *     ..
 *     .. Local Scalars ..
       REAL               AAPP, AAPP0, AAPQ, AAQQ, APOAQ, AQOAP, BIG,
-     +                   BIGTHETA, CS, CTOL, EPSILON, LARGE, MXAAPQ,
+     +                   BIGTHETA, CS, CTOL, EPSLN, LARGE, MXAAPQ,
      +                   MXSINJ, ROOTBIG, ROOTEPS, ROOTSFMIN, ROOTTOL,
-     +                   SCALE, SFMIN, SMALL, SN, T, TEMP1, THETA,
+     +                   SKL, SFMIN, SMALL, SN, T, TEMP1, THETA,
      +                   THSIGN, TOL
       INTEGER            BLSKIP, EMPTSW, i, ibr, IERR, igl, IJBLSK, ir1,
      +                   ISWROT, jbc, jgl, KBL, LKAHEAD, MVL, N2, N34,
@@ -364,20 +364,20 @@
 *     ... and the machine dependent parameters are
 *[!]  (Make sure that SLAMCH() works properly on the target machine.)
 *
-      EPSILON = SLAMCH( 'Epsilon' )
-      ROOTEPS = SQRT( EPSILON )
+      EPSLN = SLAMCH( 'Epsilon' )
+      ROOTEPS = SQRT( EPSLN )
       SFMIN = SLAMCH( 'SafeMinimum' )
       ROOTSFMIN = SQRT( SFMIN )
-      SMALL = SFMIN / EPSILON
+      SMALL = SFMIN / EPSLN
       BIG = SLAMCH( 'Overflow' )
       ROOTBIG = ONE / ROOTSFMIN
       LARGE = BIG / SQRT( FLOAT( M*N ) )
       BIGTHETA = ONE / ROOTEPS
 *
-      TOL = CTOL*EPSILON
+      TOL = CTOL*EPSLN
       ROOTTOL = SQRT( TOL )
 *
-      IF( FLOAT( M )*EPSILON.GE.ONE ) THEN
+      IF( FLOAT( M )*EPSLN.GE.ONE ) THEN
          INFO = -5
          CALL XERBLA( 'SGESVJ', -INFO )
          RETURN
@@ -402,7 +402,7 @@
 *     SQRT(N)*max_i SVA(i) does not overflow. If INFinite entries
 *     in A are detected, the procedure returns with INFO=-6.
 *
-      SCALE = ONE / SQRT( FLOAT( M )*FLOAT( N ) )
+      SKL = ONE / SQRT( FLOAT( M )*FLOAT( N ) )
       NOSCALE = .TRUE.
       GOSCALE = .TRUE.
 *
@@ -410,7 +410,7 @@
 *        the input matrix is M-by-N lower triangular (trapezoidal)
          DO 1874 p = 1, N
             AAPP = ZERO
-            AAQQ = ZERO
+            AAQQ = ONE
             CALL SLASSQ( M-p+1, A( p, p ), 1, AAPP, AAQQ )
             IF( AAPP.GT.BIG ) THEN
                INFO = -6
@@ -422,11 +422,11 @@
                SVA( p ) = AAPP*AAQQ
             ELSE
                NOSCALE = .FALSE.
-               SVA( p ) = AAPP*( AAQQ*SCALE )
+               SVA( p ) = AAPP*( AAQQ*SKL )
                IF( GOSCALE ) THEN
                   GOSCALE = .FALSE.
                   DO 1873 q = 1, p - 1
-                     SVA( q ) = SVA( q )*SCALE
+                     SVA( q ) = SVA( q )*SKL
  1873             CONTINUE
                END IF
             END IF
@@ -435,7 +435,7 @@
 *        the input matrix is M-by-N upper triangular (trapezoidal)
          DO 2874 p = 1, N
             AAPP = ZERO
-            AAQQ = ZERO
+            AAQQ = ONE
             CALL SLASSQ( p, A( 1, p ), 1, AAPP, AAQQ )
             IF( AAPP.GT.BIG ) THEN
                INFO = -6
@@ -447,11 +447,11 @@
                SVA( p ) = AAPP*AAQQ
             ELSE
                NOSCALE = .FALSE.
-               SVA( p ) = AAPP*( AAQQ*SCALE )
+               SVA( p ) = AAPP*( AAQQ*SKL )
                IF( GOSCALE ) THEN
                   GOSCALE = .FALSE.
                   DO 2873 q = 1, p - 1
-                     SVA( q ) = SVA( q )*SCALE
+                     SVA( q ) = SVA( q )*SKL
  2873             CONTINUE
                END IF
             END IF
@@ -460,7 +460,7 @@
 *        the input matrix is M-by-N general dense
          DO 3874 p = 1, N
             AAPP = ZERO
-            AAQQ = ZERO
+            AAQQ = ONE
             CALL SLASSQ( M, A( 1, p ), 1, AAPP, AAQQ )
             IF( AAPP.GT.BIG ) THEN
                INFO = -6
@@ -472,18 +472,18 @@
                SVA( p ) = AAPP*AAQQ
             ELSE
                NOSCALE = .FALSE.
-               SVA( p ) = AAPP*( AAQQ*SCALE )
+               SVA( p ) = AAPP*( AAQQ*SKL )
                IF( GOSCALE ) THEN
                   GOSCALE = .FALSE.
                   DO 3873 q = 1, p - 1
-                     SVA( q ) = SVA( q )*SCALE
+                     SVA( q ) = SVA( q )*SKL
  3873             CONTINUE
                END IF
             END IF
  3874    CONTINUE
       END IF
 *
-      IF( NOSCALE )SCALE = ONE
+      IF( NOSCALE )SKL = ONE
 *
 *     Move the smaller part of the spectrum from the underflow threshold
 *(!)  Start by determining the position of the nonzero entries of the
@@ -512,9 +512,9 @@
 * #:) Quick return for one-column matrix
 *
       IF( N.EQ.1 ) THEN
-         IF( LSVEC )CALL SLASCL( 'G', 0, 0, SVA( 1 ), SCALE, M, 1,
+         IF( LSVEC )CALL SLASCL( 'G', 0, 0, SVA( 1 ), SKL, M, 1,
      +                           A( 1, 1 ), LDA, IERR )
-         WORK( 1 ) = ONE / SCALE
+         WORK( 1 ) = ONE / SKL
          IF( SVA( 1 ).GE.SFMIN ) THEN
             WORK( 2 ) = ONE
          ELSE
@@ -530,7 +530,7 @@
 *     Protect small singular values from underflow, and try to
 *     avoid underflows/overflows in computing Jacobi rotations.
 *
-      SN = SQRT( SFMIN / EPSILON )
+      SN = SQRT( SFMIN / EPSLN )
       TEMP1 = SQRT( BIG / FLOAT( N ) )
       IF( ( AAPP.LE.SN ) .OR. ( AAQQ.GE.TEMP1 ) .OR.
      +    ( ( SN.LE.AAQQ ) .AND. ( AAPP.LE.TEMP1 ) ) ) THEN
@@ -558,10 +558,10 @@
       IF( TEMP1.NE.ONE ) THEN
          CALL SLASCL( 'G', 0, 0, ONE, TEMP1, N, 1, SVA, N, IERR )
       END IF
-      SCALE = TEMP1*SCALE
-      IF( SCALE.NE.ONE ) THEN
-         CALL SLASCL( JOBA, 0, 0, ONE, SCALE, M, N, A, LDA, IERR )
-         SCALE = ONE / SCALE
+      SKL = TEMP1*SKL
+      IF( SKL.NE.ONE ) THEN
+         CALL SLASCL( JOBA, 0, 0, ONE, SKL, M, N, A, LDA, IERR )
+         SKL = ONE / SKL
       END IF
 *
 *     Row-cyclic Jacobi SVD algorithm with column pivoting
@@ -634,30 +634,30 @@
 *
             CALL SGSVJ0( JOBV, M-N34, N-N34, A( N34+1, N34+1 ), LDA,
      +                   WORK( N34+1 ), SVA( N34+1 ), MVL,
-     +                   V( N34*q+1, N34+1 ), LDV, EPSILON, SFMIN, TOL,
+     +                   V( N34*q+1, N34+1 ), LDV, EPSLN, SFMIN, TOL,
      +                   2, WORK( N+1 ), LWORK-N, IERR )
 *
             CALL SGSVJ0( JOBV, M-N2, N34-N2, A( N2+1, N2+1 ), LDA,
      +                   WORK( N2+1 ), SVA( N2+1 ), MVL,
-     +                   V( N2*q+1, N2+1 ), LDV, EPSILON, SFMIN, TOL, 2,
+     +                   V( N2*q+1, N2+1 ), LDV, EPSLN, SFMIN, TOL, 2,
      +                   WORK( N+1 ), LWORK-N, IERR )
 *
             CALL SGSVJ1( JOBV, M-N2, N-N2, N4, A( N2+1, N2+1 ), LDA,
      +                   WORK( N2+1 ), SVA( N2+1 ), MVL,
-     +                   V( N2*q+1, N2+1 ), LDV, EPSILON, SFMIN, TOL, 1,
+     +                   V( N2*q+1, N2+1 ), LDV, EPSLN, SFMIN, TOL, 1,
      +                   WORK( N+1 ), LWORK-N, IERR )
 *
             CALL SGSVJ0( JOBV, M-N4, N2-N4, A( N4+1, N4+1 ), LDA,
      +                   WORK( N4+1 ), SVA( N4+1 ), MVL,
-     +                   V( N4*q+1, N4+1 ), LDV, EPSILON, SFMIN, TOL, 1,
+     +                   V( N4*q+1, N4+1 ), LDV, EPSLN, SFMIN, TOL, 1,
      +                   WORK( N+1 ), LWORK-N, IERR )
 *
             CALL SGSVJ0( JOBV, M, N4, A, LDA, WORK, SVA, MVL, V, LDV,
-     +                   EPSILON, SFMIN, TOL, 1, WORK( N+1 ), LWORK-N,
+     +                   EPSLN, SFMIN, TOL, 1, WORK( N+1 ), LWORK-N,
      +                   IERR )
 *
             CALL SGSVJ1( JOBV, M, N2, N4, A, LDA, WORK, SVA, MVL, V,
-     +                   LDV, EPSILON, SFMIN, TOL, 1, WORK( N+1 ),
+     +                   LDV, EPSLN, SFMIN, TOL, 1, WORK( N+1 ),
      +                   LWORK-N, IERR )
 *
 *
@@ -665,21 +665,21 @@
 *
 *
             CALL SGSVJ0( JOBV, N4, N4, A, LDA, WORK, SVA, MVL, V, LDV,
-     +                   EPSILON, SFMIN, TOL, 2, WORK( N+1 ), LWORK-N,
+     +                   EPSLN, SFMIN, TOL, 2, WORK( N+1 ), LWORK-N,
      +                   IERR )
 *
             CALL SGSVJ0( JOBV, N2, N4, A( 1, N4+1 ), LDA, WORK( N4+1 ),
      +                   SVA( N4+1 ), MVL, V( N4*q+1, N4+1 ), LDV,
-     +                   EPSILON, SFMIN, TOL, 1, WORK( N+1 ), LWORK-N,
+     +                   EPSLN, SFMIN, TOL, 1, WORK( N+1 ), LWORK-N,
      +                   IERR )
 *
             CALL SGSVJ1( JOBV, N2, N2, N4, A, LDA, WORK, SVA, MVL, V,
-     +                   LDV, EPSILON, SFMIN, TOL, 1, WORK( N+1 ),
+     +                   LDV, EPSLN, SFMIN, TOL, 1, WORK( N+1 ),
      +                   LWORK-N, IERR )
 *
             CALL SGSVJ0( JOBV, N2+N4, N4, A( 1, N2+1 ), LDA,
      +                   WORK( N2+1 ), SVA( N2+1 ), MVL,
-     +                   V( N2*q+1, N2+1 ), LDV, EPSILON, SFMIN, TOL, 1,
+     +                   V( N2*q+1, N2+1 ), LDV, EPSLN, SFMIN, TOL, 1,
      +                   WORK( N+1 ), LWORK-N, IERR )
 
          END IF
@@ -747,7 +747,7 @@
                         SVA( p ) = SNRM2( M, A( 1, p ), 1 )*WORK( p )
                      ELSE
                         TEMP1 = ZERO
-                        AAPP = ZERO
+                        AAPP = ONE
                         CALL SLASSQ( M, A( 1, p ), 1, TEMP1, AAPP )
                         SVA( p ) = TEMP1*SQRT( AAPP )*WORK( p )
                      END IF
@@ -834,7 +834,8 @@
      +                                              FASTR )
                                     SVA( q ) = AAQQ*SQRT( AMAX1( ZERO,
      +                                         ONE+T*APOAQ*AAPQ ) )
-                                    AAPP = AAPP*SQRT( ONE-T*AQOAP*AAPQ )
+                                    AAPP = AAPP*SQRT( AMAX1( ZERO, 
+     +                                         ONE-T*AQOAP*AAPQ ) )
                                     MXSINJ = AMAX1( MXSINJ, ABS( T ) )
 *
                                  ELSE
@@ -981,7 +982,7 @@
      +                                         WORK( q )
                                  ELSE
                                     T = ZERO
-                                    AAQQ = ZERO
+                                    AAQQ = ONE
                                     CALL SLASSQ( M, A( 1, q ), 1, T,
      +                                           AAQQ )
                                     SVA( q ) = T*SQRT( AAQQ )*WORK( q )
@@ -994,7 +995,7 @@
      +                                     WORK( p )
                                  ELSE
                                     T = ZERO
-                                    AAPP = ZERO
+                                    AAPP = ONE
                                     CALL SLASSQ( M, A( 1, p ), 1, T,
      +                                           AAPP )
                                     AAPP = T*SQRT( AAPP )*WORK( p )
@@ -1155,7 +1156,8 @@
                                     MXSINJ = AMAX1( MXSINJ, ABS( SN ) )
                                     SVA( q ) = AAQQ*SQRT( AMAX1( ZERO,
      +                                         ONE+T*APOAQ*AAPQ ) )
-                                    AAPP = AAPP*SQRT( ONE-T*AQOAP*AAPQ )
+                                    AAPP = AAPP*SQRT( AMAX1( ZERO,  
+     +                                         ONE-T*AQOAP*AAPQ ) )
 *
                                     APOAQ = WORK( p ) / WORK( q )
                                     AQOAP = WORK( q ) / WORK( p )
@@ -1306,7 +1308,7 @@
      +                                         WORK( q )
                                  ELSE
                                     T = ZERO
-                                    AAQQ = ZERO
+                                    AAQQ = ONE
                                     CALL SLASSQ( M, A( 1, q ), 1, T,
      +                                           AAQQ )
                                     SVA( q ) = T*SQRT( AAQQ )*WORK( q )
@@ -1319,7 +1321,7 @@
      +                                     WORK( p )
                                  ELSE
                                     T = ZERO
-                                    AAPP = ZERO
+                                    AAPP = ONE
                                     CALL SLASSQ( M, A( 1, p ), 1, T,
      +                                           AAPP )
                                     AAPP = T*SQRT( AAPP )*WORK( p )
@@ -1385,7 +1387,7 @@
             SVA( N ) = SNRM2( M, A( 1, N ), 1 )*WORK( N )
          ELSE
             T = ZERO
-            AAPP = ZERO
+            AAPP = ONE
             CALL SLASSQ( M, A( 1, N ), 1, T, AAPP )
             SVA( N ) = T*SQRT( AAPP )*WORK( N )
          END IF
@@ -1436,12 +1438,12 @@
          END IF
          IF( SVA( p ).NE.ZERO ) THEN
             N4 = N4 + 1
-            IF( SVA( p )*SCALE.GT.SFMIN )N2 = N2 + 1
+            IF( SVA( p )*SKL.GT.SFMIN )N2 = N2 + 1
          END IF
  5991 CONTINUE
       IF( SVA( N ).NE.ZERO ) THEN
          N4 = N4 + 1
-         IF( SVA( N )*SCALE.GT.SFMIN )N2 = N2 + 1
+         IF( SVA( N )*SKL.GT.SFMIN )N2 = N2 + 1
       END IF
 *
 *     Normalize the left singular vectors.
@@ -1468,17 +1470,17 @@
       END IF
 *
 *     Undo scaling, if necessary (and possible).
-      IF( ( ( SCALE.GT.ONE ) .AND. ( SVA( 1 ).LT.( BIG /
-     +    SCALE ) ) ) .OR. ( ( SCALE.LT.ONE ) .AND. ( SVA( N2 ).GT.
-     +    ( SFMIN / SCALE ) ) ) ) THEN
+      IF( ( ( SKL.GT.ONE ) .AND. ( SVA( 1 ).LT.( BIG /
+     +    SKL ) ) ) .OR. ( ( SKL.LT.ONE ) .AND. ( SVA( N2 ).GT.
+     +    ( SFMIN / SKL ) ) ) ) THEN
          DO 2400 p = 1, N
-            SVA( p ) = SCALE*SVA( p )
+            SVA( p ) = SKL*SVA( p )
  2400    CONTINUE
-         SCALE = ONE
+         SKL = ONE
       END IF
 *
-      WORK( 1 ) = SCALE
-*     The singular values of A are SCALE*SVA(1:N). If SCALE.NE.ONE
+      WORK( 1 ) = SKL
+*     The singular values of A are SKL*SVA(1:N). If SKL.NE.ONE
 *     then some of the singular values may overflow or underflow and
 *     the spectrum is given in this factored representation.
 *

diff --git a/SRC/sgsvj0.f b/SRC/sgsvj0.f
index 1adcf29..44833fb 100644 (file)
--- a/SRC/sgsvj0.f
+++ b/SRC/sgsvj0.f
@@ -186,9 +186,10 @@
          INFO = -3
       ELSE IF( LDA.LT.M ) THEN
          INFO = -5
-      ELSE IF( MV.LT.0 ) THEN
+      ELSE IF( ( RSVEC.OR.APPLV ) .AND. ( MV.LT.0 ) ) THEN
          INFO = -8
-      ELSE IF( LDV.LT.M ) THEN
+      ELSE IF( ( RSVEC.AND.( LDV.LT.N ) ).OR. 
+     &         ( APPLV.AND.( LDV.LT.MV ) ) ) THEN
          INFO = -10
       ELSE IF( TOL.LE.EPS ) THEN
          INFO = -13
@@ -310,7 +311,7 @@
                         SVA( p ) = SNRM2( M, A( 1, p ), 1 )*D( p )
                      ELSE
                         TEMP1 = ZERO
-                        AAPP = ZERO
+                        AAPP = ONE
                         CALL SLASSQ( M, A( 1, p ), 1, TEMP1, AAPP )
                         SVA( p ) = TEMP1*SQRT( AAPP )*D( p )
                      END IF
@@ -393,7 +394,8 @@
      +                                              FASTR )
                                     SVA( q ) = AAQQ*SQRT( AMAX1( ZERO,
      +                                         ONE+T*APOAQ*AAPQ ) )
-                                    AAPP = AAPP*SQRT( ONE-T*AQOAP*AAPQ )
+                                    AAPP = AAPP*SQRT( AMAX1( ZERO, 
+     +                                         ONE-T*AQOAP*AAPQ ) )
                                     MXSINJ = AMAX1( MXSINJ, ABS( T ) )
 *
                                  ELSE
@@ -536,7 +538,7 @@
      +                                         D( q )
                                  ELSE
                                     T = ZERO
-                                    AAQQ = ZERO
+                                    AAQQ = ONE
                                     CALL SLASSQ( M, A( 1, q ), 1, T,
      +                                           AAQQ )
                                     SVA( q ) = T*SQRT( AAQQ )*D( q )
@@ -549,7 +551,7 @@
      +                                     D( p )
                                  ELSE
                                     T = ZERO
-                                    AAPP = ZERO
+                                    AAPP = ONE
                                     CALL SLASSQ( M, A( 1, p ), 1, T,
      +                                           AAPP )
                                     AAPP = T*SQRT( AAPP )*D( p )
@@ -707,7 +709,8 @@
                                     MXSINJ = AMAX1( MXSINJ, ABS( SN ) )
                                     SVA( q ) = AAQQ*SQRT( AMAX1( ZERO,
      +                                         ONE+T*APOAQ*AAPQ ) )
-                                    AAPP = AAPP*SQRT( ONE-T*AQOAP*AAPQ )
+                                    AAPP = AAPP*SQRT( AMAX1( ZERO, 
+     +                                         ONE-T*AQOAP*AAPQ ) )
 *
                                     APOAQ = D( p ) / D( q )
                                     AQOAP = D( q ) / D( p )
@@ -855,7 +858,7 @@
      +                                         D( q )
                                  ELSE
                                     T = ZERO
-                                    AAQQ = ZERO
+                                    AAQQ = ONE
                                     CALL SLASSQ( M, A( 1, q ), 1, T,
      +                                           AAQQ )
                                     SVA( q ) = T*SQRT( AAQQ )*D( q )
@@ -868,7 +871,7 @@
      +                                     D( p )
                                  ELSE
                                     T = ZERO
-                                    AAPP = ZERO
+                                    AAPP = ONE
                                     CALL SLASSQ( M, A( 1, p ), 1, T,
      +                                           AAPP )
                                     AAPP = T*SQRT( AAPP )*D( p )
@@ -931,7 +934,7 @@
             SVA( N ) = SNRM2( M, A( 1, N ), 1 )*D( N )
          ELSE
             T = ZERO
-            AAPP = ZERO
+            AAPP = ONE
             CALL SLASSQ( M, A( 1, N ), 1, T, AAPP )
             SVA( N ) = T*SQRT( AAPP )*D( N )
          END IF

diff --git a/SRC/sgsvj1.f b/SRC/sgsvj1.f
index e8707ff..0021d57 100644 (file)
--- a/SRC/sgsvj1.f
+++ b/SRC/sgsvj1.f
@@ -205,9 +205,10 @@
          INFO = -4
       ELSE IF( LDA.LT.M ) THEN
          INFO = -6
-      ELSE IF( MV.LT.0 ) THEN
+      ELSE IF( ( RSVEC.OR.APPLV ) .AND. ( MV.LT.0 ) ) THEN
          INFO = -9
-      ELSE IF( LDV.LT.M ) THEN
+      ELSE IF( ( RSVEC.AND.( LDV.LT.N ) ).OR. 
+     &         ( APPLV.AND.( LDV.LT.MV ) )  ) THEN
          INFO = -11
       ELSE IF( TOL.LE.EPS ) THEN
          INFO = -14
@@ -405,7 +406,8 @@
                                     MXSINJ = AMAX1( MXSINJ, ABS( SN ) )
                                     SVA( q ) = AAQQ*SQRT( AMAX1( ZERO,
      +                                         ONE+T*APOAQ*AAPQ ) )
-                                    AAPP = AAPP*SQRT( ONE-T*AQOAP*AAPQ )
+                                    AAPP = AAPP*SQRT( AMAX1( ZERO, 
+     +                                         ONE-T*AQOAP*AAPQ ) )
 
                                     APOAQ = D( p ) / D( q )
                                     AQOAP = D( q ) / D( p )
@@ -553,7 +555,7 @@
      +                                         D( q )
                                  ELSE
                                     T = ZERO
-                                    AAQQ = ZERO
+                                    AAQQ = ONE
                                     CALL SLASSQ( M, A( 1, q ), 1, T,
      +                                           AAQQ )
                                     SVA( q ) = T*SQRT( AAQQ )*D( q )
@@ -566,7 +568,7 @@
      +                                     D( p )
                                  ELSE
                                     T = ZERO
-                                    AAPP = ZERO
+                                    AAPP = ONE
                                     CALL SLASSQ( M, A( 1, p ), 1, T,
      +                                           AAPP )
                                     AAPP = T*SQRT( AAPP )*D( p )
@@ -633,7 +635,7 @@
             SVA( N ) = SNRM2( M, A( 1, N ), 1 )*D( N )
          ELSE
             T = ZERO
-            AAPP = ZERO
+            AAPP = ONE
             CALL SLASSQ( M, A( 1, N ), 1, T, AAPP )
             SVA( N ) = T*SQRT( AAPP )*D( N )
          END IF