--- src/sgelsd.f	Thu Dec  5 20:18:10 2002
+++ debian/updates/src/sgelsd.f	Wed Dec  4 15:23:42 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 ..
       INTEGER            INFO, LDA, LDB, LWORK, M, N, NRHS, RANK
@@ -61,9 +62,10 @@
 *          The number of right hand sides, i.e., the number of columns
 *          of the matrices B and X. NRHS >= 0.
 *
-*  A       (input) REAL array, dimension (LDA,N)
+*  A       (input/output) REAL array, dimension (LDA,N)
 *          On entry, the M-by-N matrix A.
-*          On exit, A has been destroyed.
+*          On exit, the first min(m,n) rows of A are overwritten with
+*          its right singular vectors, stored rowwise.
 *
 *  LDA     (input) INTEGER
 *          The leading dimension of the array A.  LDA >= max(1,M).
@@ -95,24 +97,20 @@
 *          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
 *
 *  LWORK   (input) INTEGER
-*          The dimension of the array WORK. LWORK must be at least 1.
+*          The dimension of the array WORK. LWORK >= 1.
 *          The exact minimum amount of workspace needed depends on M,
-*          N and NRHS. As long as LWORK is at least
-*              12*N + 2*N*SMLSIZ + 8*N*NLVL + N*NRHS + (SMLSIZ+1)**2,
-*          if M is greater than or equal to N or
-*              12*M + 2*M*SMLSIZ + 8*M*NLVL + M*NRHS + (SMLSIZ+1)**2,
-*          if M is less than N, the code will execute correctly.
+*          N and NRHS.
+*          If M >= N, LWORK >=  11*N + 2*N*SMLSIZ + 8*N*NLVL + N*NRHS.
+*          If M < N, LWORK >=  11*M + 2*M*SMLSIZ + 8*M*NLVL + M*NRHS.
 *          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), and
-*             NLVL = MAX( 0, INT( LOG_2( MIN( M,N )/(SMLSIZ+1) ) ) + 1 )
+*              NLVL = INT( LOG_2( MIN( M,N )/(SMLSIZ+1) ) ) + 1
 *          For good performance, LWORK should generally be larger.
 *
-*          If LWORK = -1, then a workspace query is assumed; the routine
-*          only calculates the optimal size of the WORK array, returns
-*          this value as the first entry of the WORK array, and no error
-*          message related to LWORK is issued by XERBLA.
-*
+*          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 (LIWORK)
 *          LIWORK >= 3 * MINMN * NLVL + 11 * MINMN,
@@ -136,14 +134,15 @@
 *  =====================================================================
 *
 *     .. Parameters ..
+      INTEGER            LQUERV
+      PARAMETER          ( LQUERV = -1 )
       REAL               ZERO, ONE, TWO
       PARAMETER          ( ZERO = 0.0E0, ONE = 1.0E0, TWO = 2.0E0 )
 *     ..
 *     .. Local Scalars ..
-      LOGICAL            LQUERY
       INTEGER            IASCL, IBSCL, IE, IL, ITAU, ITAUP, ITAUQ,
      $                   LDWORK, MAXMN, MAXWRK, MINMN, MINWRK, MM,
-     $                   MNTHR, NLVL, NWORK, SMLSIZ, WLALSD
+     $                   MNTHR, NLVL, NWORK, SMLSIZ
       REAL               ANRM, BIGNUM, BNRM, EPS, SFMIN, SMLNUM
 *     ..
 *     .. External Subroutines ..
@@ -166,7 +165,6 @@
       MINMN = MIN( M, N )
       MAXMN = MAX( M, N )
       MNTHR = ILAENV( 6, 'SGELSD', ' ', M, N, NRHS, -1 )
-      LQUERY = ( LWORK.EQ.-1 )
       IF( M.LT.0 ) THEN
          INFO = -1
       ELSE IF( N.LT.0 ) THEN
@@ -190,8 +188,8 @@
 *
       MINWRK = 1
       MINMN = MAX( 1, MINMN )
-      NLVL = MAX( INT( LOG( REAL( MINMN ) / REAL( SMLSIZ+1 ) ) / 
-     $       LOG( TWO ) ) + 1, 0 )
+      NLVL = INT( LOG( REAL( MINMN ) / REAL( SMLSIZ+1 ) ) / 
+     $       LOG( TWO ) ) + 1
 *
       IF( INFO.EQ.0 ) THEN
          MAXWRK = 0
@@ -216,12 +214,11 @@
      $               ILAENV( 1, 'SORMBR', 'QLT', MM, NRHS, N, -1 ) )
             MAXWRK = MAX( MAXWRK, 3*N+( N-1 )*
      $               ILAENV( 1, 'SORMBR', 'PLN', N, NRHS, N, -1 ) )
-            WLALSD = 9*N+2*N*SMLSIZ+8*N*NLVL+N*NRHS+(SMLSIZ+1)**2
-            MAXWRK = MAX( MAXWRK, 3*N+WLALSD )
-            MINWRK = MAX( 3*N+MM, 3*N+NRHS, 3*N+WLALSD )
+            MAXWRK = MAX( MAXWRK, 3*N+8*N+2*N*SMLSIZ+8*N*NLVL+N*NRHS )
+            MINWRK = MAX( 3*N+MM, 3*N+NRHS,
+     $               3*N+8*N+2*N*SMLSIZ+8*N*NLVL+N*NRHS )
          END IF
          IF( N.GT.M ) THEN
-            WLALSD = 9*M+2*M*SMLSIZ+8*M*NLVL+M*NRHS+(SMLSIZ+1)**2
             IF( N.GE.MNTHR ) THEN
 *
 *              Path 2a - underdetermined, with many more columns
@@ -241,7 +238,8 @@
                END IF
                MAXWRK = MAX( MAXWRK, M+NRHS*
      $                  ILAENV( 1, 'SORMLQ', 'LT', N, NRHS, M, -1 ) )
-               MAXWRK = MAX( MAXWRK, M*M+4*M+WLALSD )
+               MAXWRK = MAX( MAXWRK, M*M+4*M+8*M+2*M*SMLSIZ+8*M*NLVL+M*
+     $                  NRHS )
             ELSE
 *
 *              Path 2 - remaining underdetermined cases.
@@ -252,26 +250,25 @@
      $                  ILAENV( 1, 'SORMBR', 'QLT', M, NRHS, N, -1 ) )
                MAXWRK = MAX( MAXWRK, 3*M+M*
      $                  ILAENV( 1, 'SORMBR', 'PLN', N, NRHS, M, -1 ) )
-               MAXWRK = MAX( MAXWRK, 3*M+WLALSD )
+               MAXWRK = MAX( MAXWRK, 3*M+8*M+2*M*SMLSIZ+8*M*NLVL+M*
+     $                  NRHS )
             END IF
-            MINWRK = MAX( 3*M+NRHS, 3*M+M, 3*M+WLALSD )
+            MINWRK = MAX( 3*M+NRHS, 3*M+M,
+     $               3*M+8*M+2*M*SMLSIZ+8*M*NLVL+M*NRHS )
          END IF
          MINWRK = MIN( MINWRK, MAXWRK )
          WORK( 1 ) = MAXWRK
-         IF( LWORK.LT.MINWRK .AND. .NOT.LQUERY ) THEN
-            INFO = -12
-         END IF
+         IF( LWORK.LT.MINWRK .AND. LWORK.NE.LQUERV )
+     $      INFO = -12
       END IF
 *
+*     Quick returns
+*
       IF( INFO.NE.0 ) THEN
          CALL XERBLA( 'SGELSD', -INFO )
          RETURN
-      ELSE IF( LQUERY ) THEN
-         GO TO 10
       END IF
-*
-*     Quick return if possible.
-*
+      IF( LWORK.EQ.LQUERV ) RETURN
       IF( M.EQ.0 .OR. N.EQ.0 ) THEN
          RANK = 0
          RETURN
