--- src/dgesdd.f	Thu Dec  5 20:18:05 2002
+++ debian/updates/src/dgesdd.f	Thu Dec  5 10:33:38 2002
@@ -4,7 +4,8 @@
 *  -- LAPACK driver routine (version 3.0) --
 *     Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
 *     Courant Institute, Argonne National Lab, and Rice University
-*     October 31, 1999
+*     June 30, 1999
+*     8-15-00:  Improve consistency of WS calculations (eca)
 *
 *     .. Scalar Arguments ..
       CHARACTER          JOBZ
@@ -116,16 +117,20 @@
 *  LWORK   (input) INTEGER
 *          The dimension of the array WORK. LWORK >= 1.
 *          If JOBZ = 'N',
-*            LWORK >= 3*min(M,N) + max(max(M,N),6*min(M,N)).
+*            LWORK >= max(14*min(M,N)+4, 10*min(M,N)+2+
+*                     SMLSIZ*(SMLSIZ+8)) + max(M,N)
+*          where SMLSIZ is returned by ILAENV and is equal to the
+*          maximum size of the subproblems at the bottom of the
+*          computation tree (usually about 25).
 *          If JOBZ = 'O',
-*            LWORK >= 3*min(M,N)*min(M,N) + 
-*                     max(max(M,N),5*min(M,N)*min(M,N)+4*min(M,N)).
+*            LWORK >= 5*min(M,N)*min(M,N) + max(M,N) + 9*min(M,N).
 *          If JOBZ = 'S' or 'A'
-*            LWORK >= 3*min(M,N)*min(M,N) +
-*                     max(max(M,N),4*min(M,N)*min(M,N)+4*min(M,N)).
+*            LWORK >= 4*min(M,N)*min(M,N) + max(M,N) + 9*min(M,N).
 *          For good performance, LWORK should generally be larger.
-*          If LWORK < 0 but other input arguments are legal, WORK(1)
-*          returns the optimal LWORK.
+*
+*          If LWORK = -1, a workspace query is assumed.  The optimal
+*          size for the WORK array is calculated and stored in WORK(1),
+*          and no other work except argument checking is performed.
 *
 *  IWORK   (workspace) INTEGER array, dimension (8*min(M,N))
 *
@@ -144,15 +149,17 @@
 *  =====================================================================
 *
 *     .. Parameters ..
+      INTEGER            LQUERV
+      PARAMETER          ( LQUERV = -1 )
       DOUBLE PRECISION   ZERO, ONE
-      PARAMETER          ( ZERO = 0.0D0, ONE = 1.0D0 )
+      PARAMETER          ( ZERO = 0.0D+0, ONE = 1.0D+0 )
 *     ..
 *     .. Local Scalars ..
-      LOGICAL            LQUERY, WNTQA, WNTQAS, WNTQN, WNTQO, WNTQS
-      INTEGER            BDSPAC, BLK, CHUNK, I, IE, IERR, IL,
+      LOGICAL            WNTQA, WNTQAS, WNTQN, WNTQO, WNTQS
+      INTEGER            BDSPAC, BDSPAN, BLK, CHUNK, I, IE, IERR, IL,
      $                   IR, ISCL, ITAU, ITAUP, ITAUQ, IU, IVT, LDWKVT,
      $                   LDWRKL, LDWRKR, LDWRKU, MAXWRK, MINMN, MINWRK,
-     $                   MNTHR, NWORK, WRKBL
+     $                   MNTHR, NWORK, SMLSIZ, WRKBL
       DOUBLE PRECISION   ANRM, BIGNUM, EPS, SMLNUM
 *     ..
 *     .. Local Arrays ..
@@ -168,7 +175,7 @@
       LOGICAL            LSAME
       INTEGER            ILAENV
       DOUBLE PRECISION   DLAMCH, DLANGE
-      EXTERNAL           DLAMCH, DLANGE, ILAENV, LSAME
+      EXTERNAL           LSAME, ILAENV, DLAMCH, DLANGE
 *     ..
 *     .. Intrinsic Functions ..
       INTRINSIC          DBLE, INT, MAX, MIN, SQRT
@@ -187,7 +194,6 @@
       WNTQN = LSAME( JOBZ, 'N' )
       MINWRK = 1
       MAXWRK = 1
-      LQUERY = ( LWORK.EQ.-1 )
 *
       IF( .NOT.( WNTQA .OR. WNTQS .OR. WNTQO .OR. WNTQN ) ) THEN
          INFO = -1
@@ -206,6 +212,8 @@
          INFO = -10
       END IF
 *
+      SMLSIZ = ILAENV( 9, 'DGESDD', ' ', 0, 0, 0, 0 )
+*
 *     Compute workspace
 *      (Note: Comments in the code beginning "Workspace:" describe the
 *       minimal amount of workspace needed at that point in the code,
@@ -218,22 +226,19 @@
 *
 *           Compute space needed for DBDSDC
 *
-            IF( WNTQN ) THEN
-               BDSPAC = 7*N
-            ELSE
-               BDSPAC = 3*N*N + 4*N
-            END IF
+            BDSPAC = 3*N*N + 7*N
+            BDSPAN = MAX( 12*N+4, 8*N+2+SMLSIZ*( SMLSIZ+8 ) )
             IF( M.GE.MNTHR ) THEN
                IF( WNTQN ) THEN
 *
 *                 Path 1 (M much larger than N, JOBZ='N')
 *
-                  WRKBL = N + N*ILAENV( 1, 'DGEQRF', ' ', M, N, -1,
-     $                    -1 )
-                  WRKBL = MAX( WRKBL, 3*N+2*N*
-     $                    ILAENV( 1, 'DGEBRD', ' ', N, N, -1, -1 ) )
-                  MAXWRK = MAX( WRKBL, BDSPAC+N )
-                  MINWRK = BDSPAC + N
+                  MAXWRK = N + N*ILAENV( 1, 'DGEQRF', ' ', M, N, -1,
+     $                     -1 )
+                  MAXWRK = MAX( MAXWRK, 3*N+2*N*
+     $                     ILAENV( 1, 'DGEBRD', ' ', N, N, -1, -1 ) )
+                  MAXWRK = MAX( MAXWRK, BDSPAC )
+                  MINWRK = BDSPAC
                ELSE IF( WNTQO ) THEN
 *
 *                 Path 2 (M much larger than N, JOBZ='O')
@@ -247,9 +252,9 @@
      $                    ILAENV( 1, 'DORMBR', 'QLN', N, N, N, -1 ) )
                   WRKBL = MAX( WRKBL, 3*N+N*
      $                    ILAENV( 1, 'DORMBR', 'PRT', N, N, N, -1 ) )
-                  WRKBL = MAX( WRKBL, BDSPAC+3*N )
+                  WRKBL = MAX( WRKBL, BDSPAC+2*N )
                   MAXWRK = WRKBL + 2*N*N
-                  MINWRK = BDSPAC + 2*N*N + 3*N
+                  MINWRK = BDSPAC + 2*N*N + 2*N
                ELSE IF( WNTQS ) THEN
 *
 *                 Path 3 (M much larger than N, JOBZ='S')
@@ -263,9 +268,9 @@
      $                    ILAENV( 1, 'DORMBR', 'QLN', N, N, N, -1 ) )
                   WRKBL = MAX( WRKBL, 3*N+N*
      $                    ILAENV( 1, 'DORMBR', 'PRT', N, N, N, -1 ) )
-                  WRKBL = MAX( WRKBL, BDSPAC+3*N )
+                  WRKBL = MAX( WRKBL, BDSPAC+2*N )
                   MAXWRK = WRKBL + N*N
-                  MINWRK = BDSPAC + N*N + 3*N
+                  MINWRK = BDSPAC + N*N + 2*N
                ELSE IF( WNTQA ) THEN
 *
 *                 Path 4 (M much larger than N, JOBZ='A')
@@ -279,9 +284,9 @@
      $                    ILAENV( 1, 'DORMBR', 'QLN', N, N, N, -1 ) )
                   WRKBL = MAX( WRKBL, 3*N+N*
      $                    ILAENV( 1, 'DORMBR', 'PRT', N, N, N, -1 ) )
-                  WRKBL = MAX( WRKBL, BDSPAC+3*N )
-                  MAXWRK = WRKBL + N*N
-                  MINWRK = BDSPAC + N*N + 3*N
+                  WRKBL = MAX( WRKBL, BDSPAC+2*N )
+                  MAXWRK = N*N + WRKBL
+                  MINWRK = BDSPAC + N*N + M + N
                END IF
             ELSE
 *
@@ -289,53 +294,49 @@
 *
                WRKBL = 3*N + ( M+N )*ILAENV( 1, 'DGEBRD', ' ', M, N, -1,
      $                 -1 )
-               IF( WNTQN ) THEN
-                  MAXWRK = MAX( WRKBL, BDSPAC+3*N )
-                  MINWRK = 3*N + MAX( M, BDSPAC )
-               ELSE IF( WNTQO ) THEN
+*               MAXWRK = MAX( WRKBL, BDSPAC+3*N )
+*               MINWRK = 3*N + MAX( M, BDSPAC )
+               IF( WNTQO ) THEN
                   WRKBL = MAX( WRKBL, 3*N+N*
      $                    ILAENV( 1, 'DORMBR', 'QLN', M, N, N, -1 ) )
                   WRKBL = MAX( WRKBL, 3*N+N*
      $                    ILAENV( 1, 'DORMBR', 'PRT', N, N, N, -1 ) )
-                  WRKBL = MAX( WRKBL, BDSPAC+3*N )
+                  WRKBL = MAX( WRKBL, BDSPAC+2*N+M )
                   MAXWRK = WRKBL + M*N
-                  MINWRK = 3*N + MAX( M, N*N+BDSPAC )
+                  MINWRK = BDSPAC + N*N + 2*N + M
                ELSE IF( WNTQS ) THEN
-                  WRKBL = MAX( WRKBL, 3*N+N*
-     $                    ILAENV( 1, 'DORMBR', 'QLN', M, N, N, -1 ) )
-                  WRKBL = MAX( WRKBL, 3*N+N*
-     $                    ILAENV( 1, 'DORMBR', 'PRT', N, N, N, -1 ) )
-                  MAXWRK = MAX( WRKBL, BDSPAC+3*N )
-                  MINWRK = 3*N + MAX( M, BDSPAC )
+                  WRKBL = MAX( WRKBL, 3*N+N*
+     $                     ILAENV( 1, 'DORMBR', 'QLN', M, N, N, -1 ) )
+                  WRKBL = MAX( WRKBL, 3*N+N*
+     $                     ILAENV( 1, 'DORMBR', 'PRT', N, N, N, -1 ) )
+                  MAXWRK = MAX( WRKBL, BDSPAC+2*N+M )
+                  MINWRK = BDSPAC + 2*N + M
                ELSE IF( WNTQA ) THEN
-                  WRKBL = MAX( WRKBL, 3*N+M*
-     $                    ILAENV( 1, 'DORMBR', 'QLN', M, M, N, -1 ) )
-                  WRKBL = MAX( WRKBL, 3*N+N*
-     $                    ILAENV( 1, 'DORMBR', 'PRT', N, N, N, -1 ) )
-                  MAXWRK = MAX( MAXWRK, BDSPAC+3*N )
-                  MINWRK = 3*N + MAX( M, BDSPAC )
+                  WRKBL = MAX( WRKBL, 3*N+M*
+     $                     ILAENV( 1, 'DORMBR', 'QLN', M, M, N, -1 ) )
+                  WRKBL = MAX( WRKBL, 3*N+N*
+     $                     ILAENV( 1, 'DORMBR', 'PRT', N, N, N, -1 ) )
+                  MAXWRK = MAX( WRKBL, BDSPAC+2*N+M )
+                  MINWRK = BDSPAC + 2*N + M
                END IF
             END IF
          ELSE
 *
 *           Compute space needed for DBDSDC
 *
-            IF( WNTQN ) THEN
-               BDSPAC = 7*M
-            ELSE
-               BDSPAC = 3*M*M + 4*M
-            END IF
+            BDSPAC = 3*M*M + 7*M
+            BDSPAN = MAX( 12*M+4, 8*M+2+SMLSIZ*( SMLSIZ+8 ) )
             IF( N.GE.MNTHR ) THEN
                IF( WNTQN ) THEN
 *
 *                 Path 1t (N much larger than M, JOBZ='N')
 *
-                  WRKBL = M + M*ILAENV( 1, 'DGELQF', ' ', M, N, -1,
-     $                    -1 )
-                  WRKBL = MAX( WRKBL, 3*M+2*M*
-     $                    ILAENV( 1, 'DGEBRD', ' ', M, M, -1, -1 ) )
-                  MAXWRK = MAX( WRKBL, BDSPAC+M )
-                  MINWRK = BDSPAC + M
+                  MAXWRK = M + M*ILAENV( 1, 'DGELQF', ' ', M, N, -1,
+     $                     -1 )
+                  MAXWRK = MAX( MAXWRK, 3*M+2*M*
+     $                     ILAENV( 1, 'DGEBRD', ' ', M, M, -1, -1 ) )
+                  MAXWRK = MAX( MAXWRK, BDSPAC )
+                  MINWRK = BDSPAC
                ELSE IF( WNTQO ) THEN
 *
 *                 Path 2t (N much larger than M, JOBZ='O')
@@ -349,9 +350,9 @@
      $                    ILAENV( 1, 'DORMBR', 'QLN', M, M, M, -1 ) )
                   WRKBL = MAX( WRKBL, 3*M+M*
      $                    ILAENV( 1, 'DORMBR', 'PRT', M, M, M, -1 ) )
-                  WRKBL = MAX( WRKBL, BDSPAC+3*M )
+                  WRKBL = MAX( WRKBL, BDSPAC+2*M )
                   MAXWRK = WRKBL + 2*M*M
-                  MINWRK = BDSPAC + 2*M*M + 3*M
+                  MINWRK = BDSPAC + 2*M*M + 2*M
                ELSE IF( WNTQS ) THEN
 *
 *                 Path 3t (N much larger than M, JOBZ='S')
@@ -365,9 +366,9 @@
      $                    ILAENV( 1, 'DORMBR', 'QLN', M, M, M, -1 ) )
                   WRKBL = MAX( WRKBL, 3*M+M*
      $                    ILAENV( 1, 'DORMBR', 'PRT', M, M, M, -1 ) )
-                  WRKBL = MAX( WRKBL, BDSPAC+3*M )
+                  WRKBL = MAX( WRKBL, BDSPAC+2*M )
                   MAXWRK = WRKBL + M*M
-                  MINWRK = BDSPAC + M*M + 3*M
+                  MINWRK = BDSPAC + M*M + 2*M
                ELSE IF( WNTQA ) THEN
 *
 *                 Path 4t (N much larger than M, JOBZ='A')
@@ -381,9 +382,9 @@
      $                    ILAENV( 1, 'DORMBR', 'QLN', M, M, M, -1 ) )
                   WRKBL = MAX( WRKBL, 3*M+M*
      $                    ILAENV( 1, 'DORMBR', 'PRT', M, M, M, -1 ) )
-                  WRKBL = MAX( WRKBL, BDSPAC+3*M )
+                  WRKBL = MAX( WRKBL, BDSPAC+2*M )
                   MAXWRK = WRKBL + M*M
-                  MINWRK = BDSPAC + M*M + 3*M
+                  MINWRK = BDSPAC + M*M + M + N
                END IF
             ELSE
 *
@@ -391,52 +392,49 @@
 *
                WRKBL = 3*M + ( M+N )*ILAENV( 1, 'DGEBRD', ' ', M, N, -1,
      $                 -1 )
-               IF( WNTQN ) THEN
-                  MAXWRK = MAX( WRKBL, BDSPAC+3*M )
-                  MINWRK = 3*M + MAX( N, BDSPAC )
-               ELSE IF( WNTQO ) THEN
+*               MAXWRK = MAX( WRKBL, BDSPAC+3*M )
+*               MINWRK = 3*M + MAX( N, BDSPAC )
+               IF( WNTQO ) THEN
                   WRKBL = MAX( WRKBL, 3*M+M*
      $                    ILAENV( 1, 'DORMBR', 'QLN', M, M, N, -1 ) )
                   WRKBL = MAX( WRKBL, 3*M+M*
      $                    ILAENV( 1, 'DORMBR', 'PRT', M, N, M, -1 ) )
-                  WRKBL = MAX( WRKBL, BDSPAC+3*M )
+                  WRKBL = MAX( WRKBL, BDSPAC+2*M+N )
                   MAXWRK = WRKBL + M*N
-                  MINWRK = 3*M + MAX( N, M*M+BDSPAC )
+                  MINWRK = BDSPAC + M*M + 2*M + N
                ELSE IF( WNTQS ) THEN
-                  WRKBL = MAX( WRKBL, 3*M+M*
-     $                    ILAENV( 1, 'DORMBR', 'QLN', M, M, N, -1 ) )
-                  WRKBL = MAX( WRKBL, 3*M+M*
-     $                    ILAENV( 1, 'DORMBR', 'PRT', M, N, M, -1 ) )
-                  MAXWRK = MAX( WRKBL, BDSPAC+3*M )
-                  MINWRK = 3*M + MAX( N, BDSPAC )
+                  WRKBL = MAX( WRKBL, 3*M+M*
+     $                     ILAENV( 1, 'DORMBR', 'QLN', M, M, N, -1 ) )
+                  WRKBL = MAX( WRKBL, 3*M+M*
+     $                     ILAENV( 1, 'DORMBR', 'PRT', M, N, M, -1 ) )
+                  MAXWRK = MAX( WRKBL, BDSPAC+2*M+N )
+                  MINWRK = BDSPAC + 2*M + N
                ELSE IF( WNTQA ) THEN
-                  WRKBL = MAX( WRKBL, 3*M+M*
-     $                    ILAENV( 1, 'DORMBR', 'QLN', M, M, N, -1 ) )
-                  WRKBL = MAX( WRKBL, 3*M+M*
-     $                    ILAENV( 1, 'DORMBR', 'PRT', N, N, M, -1 ) )
-                  MAXWRK = MAX( WRKBL, BDSPAC+3*M )
-                  MINWRK = 3*M + MAX( N, BDSPAC )
+                  WRKBL = MAX( WRKBL, 3*M+M*
+     $                     ILAENV( 1, 'DORMBR', 'QLN', M, M, N, -1 ) )
+                  WRKBL = MAX( WRKBL, 3*M+N*
+     $                     ILAENV( 1, 'DORMBR', 'PRT', N, N, M, -1 ) )
+                  MAXWRK = MAX( WRKBL, BDSPAC+2*M+N )
+                  MINWRK = BDSPAC + 2*M + N
                END IF
             END IF
          END IF
+      END IF
+      IF( INFO.EQ.0 ) THEN
          WORK( 1 ) = MAXWRK
+         IF( LWORK.LT.MINWRK .AND. LWORK.NE.LQUERV )
+     $      INFO = -12
       END IF
 *
-      IF( LWORK.LT.MINWRK .AND. .NOT.LQUERY ) THEN
-         INFO = -12
-      END IF
+*     Quick returns
+*
       IF( INFO.NE.0 ) THEN
          CALL XERBLA( 'DGESDD', -INFO )
          RETURN
-      ELSE IF( LQUERY ) THEN
-         RETURN
       END IF
-*
-*     Quick return if possible
-*
+      IF( LWORK.EQ.LQUERV )
+     $   RETURN
       IF( M.EQ.0 .OR. N.EQ.0 ) THEN
-         IF( LWORK.GE.1 )
-     $      WORK( 1 ) = ONE
          RETURN
       END IF
 *
@@ -497,7 +495,7 @@
                NWORK = IE + N
 *
 *              Perform bidiagonal SVD, computing singular values only
-*              (Workspace: need N+BDSPAC)
+*              (Workspace: need BDSPAN)
 *
                CALL DBDSDC( 'U', 'N', N, S, WORK( IE ), DUM, 1, DUM, 1,
      $                      DUM, IDUM, WORK( NWORK ), IWORK, INFO )
@@ -512,10 +510,10 @@
 *
 *              WORK(IR) is LDWRKR by N
 *
-               IF( LWORK.GE.LDA*N+N*N+3*N+BDSPAC ) THEN
+               IF( LWORK.GE.LDA*N+4*N*N+9*N ) THEN
                   LDWRKR = LDA
                ELSE
-                  LDWRKR = ( LWORK-N*N-3*N-BDSPAC ) / N
+                  LDWRKR = ( LWORK-4*N*N-9*N ) / N
                END IF
                ITAU = IR + LDWRKR*N
                NWORK = ITAU + N
@@ -557,7 +555,7 @@
 *              Perform bidiagonal SVD, computing left singular vectors
 *              of bidiagonal matrix in WORK(IU) and computing right
 *              singular vectors of bidiagonal matrix in VT
-*              (Workspace: need N+N*N+BDSPAC)
+*              (Workspace: need 2*N*N+BDSPAC)
 *
                CALL DBDSDC( 'U', 'I', N, S, WORK( IE ), WORK( IU ), N,
      $                      VT, LDVT, DUM, IDUM, WORK( NWORK ), IWORK,
@@ -633,7 +631,7 @@
 *              Perform bidiagonal SVD, computing left singular vectors
 *              of bidiagoal matrix in U and computing right singular
 *              vectors of bidiagonal matrix in VT
-*              (Workspace: need N+BDSPAC)
+*              (Workspace: need N*N+BDSPAC)
 *
                CALL DBDSDC( 'U', 'I', N, S, WORK( IE ), U, LDU, VT,
      $                      LDVT, DUM, IDUM, WORK( NWORK ), IWORK,
@@ -681,7 +679,7 @@
                CALL DLACPY( 'L', M, N, A, LDA, U, LDU )
 *
 *              Generate Q in U
-*              (Workspace: need N*N+2*N, prefer N*N+N+N*NB)
+*              (Workspace: need N*N+N+M, prefer N*N+N+M*NB)
                CALL DORGQR( M, M, N, U, LDU, WORK( ITAU ),
      $                      WORK( NWORK ), LWORK-NWORK+1, IERR )
 *
@@ -703,7 +701,7 @@
 *              Perform bidiagonal SVD, computing left singular vectors
 *              of bidiagonal matrix in WORK(IU) and computing right
 *              singular vectors of bidiagonal matrix in VT
-*              (Workspace: need N+N*N+BDSPAC)
+*              (Workspace: need N*N+BDSPAC)
 *
                CALL DBDSDC( 'U', 'I', N, S, WORK( IE ), WORK( IU ), N,
      $                      VT, LDVT, DUM, IDUM, WORK( NWORK ), IWORK,
@@ -754,13 +752,13 @@
             IF( WNTQN ) THEN
 *
 *              Perform bidiagonal SVD, only computing singular values
-*              (Workspace: need N+BDSPAC)
+*              (Workspace: need BDSPAN)
 *
                CALL DBDSDC( 'U', 'N', N, S, WORK( IE ), DUM, 1, DUM, 1,
      $                      DUM, IDUM, WORK( NWORK ), IWORK, INFO )
             ELSE IF( WNTQO ) THEN
                IU = NWORK
-               IF( LWORK.GE.M*N+3*N+BDSPAC ) THEN
+               IF( LWORK.GE.M*N+3*N*N+9*N ) THEN
 *
 *                 WORK( IU ) is M by N
 *
@@ -785,7 +783,7 @@
 *              Perform bidiagonal SVD, computing left singular vectors
 *              of bidiagonal matrix in WORK(IU) and computing right
 *              singular vectors of bidiagonal matrix in VT
-*              (Workspace: need N+N*N+BDSPAC)
+*              (Workspace: need N*N+BDSPAC)
 *
                CALL DBDSDC( 'U', 'I', N, S, WORK( IE ), WORK( IU ),
      $                      LDWRKU, VT, LDVT, DUM, IDUM, WORK( NWORK ),
@@ -798,7 +796,7 @@
      $                      WORK( ITAUP ), VT, LDVT, WORK( NWORK ),
      $                      LWORK-NWORK+1, IERR )
 *
-               IF( LWORK.GE.M*N+3*N+BDSPAC ) THEN
+               IF( LWORK.GE.M*N+3*N*N+9*N ) THEN
 *
 *                 Overwrite WORK(IU) by left singular vectors of A
 *                 (Workspace: need N*N+2*N, prefer N*N+N+N*NB)
@@ -838,7 +836,7 @@
 *              Perform bidiagonal SVD, computing left singular vectors
 *              of bidiagonal matrix in U and computing right singular
 *              vectors of bidiagonal matrix in VT
-*              (Workspace: need N+BDSPAC)
+*              (Workspace: need BDSPAC)
 *
                CALL DLASET( 'F', M, N, ZERO, ZERO, U, LDU )
                CALL DBDSDC( 'U', 'I', N, S, WORK( IE ), U, LDU, VT,
@@ -855,12 +853,12 @@
                CALL DORMBR( 'P', 'R', 'T', N, N, N, A, LDA,
      $                      WORK( ITAUP ), VT, LDVT, WORK( NWORK ),
      $                      LWORK-NWORK+1, IERR )
-            ELSE IF( WNTQA ) THEN
+            ELSE
 *
 *              Perform bidiagonal SVD, computing left singular vectors
 *              of bidiagonal matrix in U and computing right singular
 *              vectors of bidiagonal matrix in VT
-*              (Workspace: need N+BDSPAC)
+*              (Workspace: need BDSPAC)
 *
                CALL DLASET( 'F', M, M, ZERO, ZERO, U, LDU )
                CALL DBDSDC( 'U', 'I', N, S, WORK( IE ), U, LDU, VT,
@@ -925,7 +923,7 @@
                NWORK = IE + M
 *
 *              Perform bidiagonal SVD, computing singular values only
-*              (Workspace: need M+BDSPAC)
+*              (Workspace: need BDSPAN)
 *
                CALL DBDSDC( 'U', 'N', M, S, WORK( IE ), DUM, 1, DUM, 1,
      $                      DUM, IDUM, WORK( NWORK ), IWORK, INFO )
@@ -941,7 +939,7 @@
 *              IVT is M by M
 *
                IL = IVT + M*M
-               IF( LWORK.GE.M*N+M*M+3*M+BDSPAC ) THEN
+               IF( LWORK.GE.M*N+4*M*M+9*M ) THEN
 *
 *                 WORK(IL) is M by N
 *
@@ -986,7 +984,7 @@
 *              Perform bidiagonal SVD, computing left singular vectors
 *              of bidiagonal matrix in U, and computing right singular
 *              vectors of bidiagonal matrix in WORK(IVT)
-*              (Workspace: need M+M*M+BDSPAC)
+*              (Workspace: need 2*M*M+BDSPAC)
 *
                CALL DBDSDC( 'U', 'I', M, S, WORK( IE ), U, LDU,
      $                      WORK( IVT ), M, DUM, IDUM, WORK( NWORK ),
@@ -1061,7 +1059,7 @@
 *              Perform bidiagonal SVD, computing left singular vectors
 *              of bidiagonal matrix in U and computing right singular
 *              vectors of bidiagonal matrix in VT
-*              (Workspace: need M+BDSPAC)
+*              (Workspace: need M*M+BDSPAC)
 *
                CALL DBDSDC( 'U', 'I', M, S, WORK( IE ), U, LDU, VT,
      $                      LDVT, DUM, IDUM, WORK( NWORK ), IWORK,
@@ -1108,7 +1106,7 @@
                CALL DLACPY( 'U', M, N, A, LDA, VT, LDVT )
 *
 *              Generate Q in VT
-*              (Workspace: need M*M+2*M, prefer M*M+M+M*NB)
+*              (Workspace: need M*M+M+N, prefer M*M+M+N*NB)
 *
                CALL DORGLQ( N, N, M, VT, LDVT, WORK( ITAU ),
      $                      WORK( NWORK ), LWORK-NWORK+1, IERR )
@@ -1131,7 +1129,7 @@
 *              Perform bidiagonal SVD, computing left singular vectors
 *              of bidiagonal matrix in U and computing right singular
 *              vectors of bidiagonal matrix in WORK(IVT)
-*              (Workspace: need M+M*M+BDSPAC)
+*              (Workspace: need M*M+BDSPAC)
 *
                CALL DBDSDC( 'U', 'I', M, S, WORK( IE ), U, LDU,
      $                      WORK( IVT ), LDWKVT, DUM, IDUM,
@@ -1182,14 +1180,14 @@
             IF( WNTQN ) THEN
 *
 *              Perform bidiagonal SVD, only computing singular values
-*              (Workspace: need M+BDSPAC)
+*              (Workspace: need BDSPAN)
 *
                CALL DBDSDC( 'L', 'N', M, S, WORK( IE ), DUM, 1, DUM, 1,
      $                      DUM, IDUM, WORK( NWORK ), IWORK, INFO )
             ELSE IF( WNTQO ) THEN
                LDWKVT = M
                IVT = NWORK
-               IF( LWORK.GE.M*N+3*M+BDSPAC ) THEN
+               IF( LWORK.GE.M*N+3*M*M+9*M ) THEN
 *
 *                 WORK( IVT ) is M by N
 *
@@ -1224,7 +1222,7 @@
      $                      WORK( ITAUQ ), U, LDU, WORK( NWORK ),
      $                      LWORK-NWORK+1, IERR )
 *
-               IF( LWORK.GE.M*N+3*M+BDSPAC ) THEN
+               IF( LWORK.GE.M*N+3*M*M+9*M ) THEN
 *
 *                 Overwrite WORK(IVT) by left singular vectors of A
 *                 (Workspace: need M*M+2*M, prefer M*M+M+M*NB)
@@ -1263,7 +1261,7 @@
 *              Perform bidiagonal SVD, computing left singular vectors
 *              of bidiagonal matrix in U and computing right singular
 *              vectors of bidiagonal matrix in VT
-*              (Workspace: need M+BDSPAC)
+*              (Workspace: need BDSPAC)
 *
                CALL DLASET( 'F', M, N, ZERO, ZERO, VT, LDVT )
                CALL DBDSDC( 'L', 'I', M, S, WORK( IE ), U, LDU, VT,
@@ -1280,12 +1278,12 @@
                CALL DORMBR( 'P', 'R', 'T', M, N, M, A, LDA,
      $                      WORK( ITAUP ), VT, LDVT, WORK( NWORK ),
      $                      LWORK-NWORK+1, IERR )
-            ELSE IF( WNTQA ) THEN
+            ELSE
 *
 *              Perform bidiagonal SVD, computing left singular vectors
 *              of bidiagonal matrix in U and computing right singular
 *              vectors of bidiagonal matrix in VT
-*              (Workspace: need M+BDSPAC)
+*              (Workspace: need BDSPAC)
 *
                CALL DLASET( 'F', N, N, ZERO, ZERO, VT, LDVT )
                CALL DBDSDC( 'L', 'I', M, S, WORK( IE ), U, LDU, VT,
@@ -1319,9 +1317,15 @@
          IF( ANRM.GT.BIGNUM )
      $      CALL DLASCL( 'G', 0, 0, BIGNUM, ANRM, MINMN, 1, S, MINMN,
      $                   IERR )
+         IF( INFO.NE.0 .AND. ANRM.GT.BIGNUM )
+     $      CALL DLASCL( 'G', 0, 0, BIGNUM, ANRM, MINMN-1, 1, WORK( 2 ),
+     $                   MINMN, IERR )
          IF( ANRM.LT.SMLNUM )
      $      CALL DLASCL( 'G', 0, 0, SMLNUM, ANRM, MINMN, 1, S, MINMN,
      $                   IERR )
+         IF( INFO.NE.0 .AND. ANRM.LT.SMLNUM )
+     $      CALL DLASCL( 'G', 0, 0, SMLNUM, ANRM, MINMN-1, 1, WORK( 2 ),
+     $                   MINMN, IERR )
       END IF
 *
 *     Return optimal workspace in WORK(1)
