Diff for /rpl/lapack/lapack/dgesvd.f between versions 1.10 and 1.20

version 1.10, 2011/11/21 22:19:28 version 1.20, 2023/08/07 08:38:50
Line 2 Line 2
 *  *
 *  =========== DOCUMENTATION ===========  *  =========== DOCUMENTATION ===========
 *  *
 * Online html documentation available at   * Online html documentation available at
 *            http://www.netlib.org/lapack/explore-html/   *            http://www.netlib.org/lapack/explore-html/
 *  *
 *> \htmlonly  *> \htmlonly
 *> Download DGESVD + dependencies   *> Download DGESVD + dependencies
 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dgesvd.f">   *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dgesvd.f">
 *> [TGZ]</a>   *> [TGZ]</a>
 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dgesvd.f">   *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dgesvd.f">
 *> [ZIP]</a>   *> [ZIP]</a>
 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dgesvd.f">   *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dgesvd.f">
 *> [TXT]</a>  *> [TXT]</a>
 *> \endhtmlonly   *> \endhtmlonly
 *  *
 *  Definition:  *  Definition:
 *  ===========  *  ===========
 *  *
 *       SUBROUTINE DGESVD( JOBU, JOBVT, M, N, A, LDA, S, U, LDU, VT, LDVT,  *       SUBROUTINE DGESVD( JOBU, JOBVT, M, N, A, LDA, S, U, LDU, VT, LDVT,
 *                          WORK, LWORK, INFO )  *                          WORK, LWORK, INFO )
 *   *
 *       .. Scalar Arguments ..  *       .. Scalar Arguments ..
 *       CHARACTER          JOBU, JOBVT  *       CHARACTER          JOBU, JOBVT
 *       INTEGER            INFO, LDA, LDU, LDVT, LWORK, M, N  *       INTEGER            INFO, LDA, LDU, LDVT, LWORK, M, N
Line 29 Line 29
 *       DOUBLE PRECISION   A( LDA, * ), S( * ), U( LDU, * ),  *       DOUBLE PRECISION   A( LDA, * ), S( * ), U( LDU, * ),
 *      $                   VT( LDVT, * ), WORK( * )  *      $                   VT( LDVT, * ), WORK( * )
 *       ..  *       ..
 *    *
 *  *
 *> \par Purpose:  *> \par Purpose:
 *  =============  *  =============
Line 173 Line 173
 *>          LWORK is INTEGER  *>          LWORK is INTEGER
 *>          The dimension of the array WORK.  *>          The dimension of the array WORK.
 *>          LWORK >= MAX(1,5*MIN(M,N)) for the paths (see comments inside code):  *>          LWORK >= MAX(1,5*MIN(M,N)) for the paths (see comments inside code):
 *>             - PATH 1  (M much larger than N, JOBU='N')   *>             - PATH 1  (M much larger than N, JOBU='N')
 *>             - PATH 1t (N much larger than M, JOBVT='N')  *>             - PATH 1t (N much larger than M, JOBVT='N')
 *>          LWORK >= MAX(1,3*MIN(M,N)+MAX(M,N),5*MIN(M,N)) for the other paths  *>          LWORK >= MAX(1,3*MIN(M,N) + MAX(M,N),5*MIN(M,N)) for the other paths
 *>          For good performance, LWORK should generally be larger.  *>          For good performance, LWORK should generally be larger.
 *>  *>
 *>          If LWORK = -1, then a workspace query is assumed; the routine  *>          If LWORK = -1, then a workspace query is assumed; the routine
Line 198 Line 198
 *  Authors:  *  Authors:
 *  ========  *  ========
 *  *
 *> \author Univ. of Tennessee   *> \author Univ. of Tennessee
 *> \author Univ. of California Berkeley   *> \author Univ. of California Berkeley
 *> \author Univ. of Colorado Denver   *> \author Univ. of Colorado Denver
 *> \author NAG Ltd.   *> \author NAG Ltd.
 *  
 *> \date November 2011  
 *  *
 *> \ingroup doubleGEsing  *> \ingroup doubleGEsing
 *  *
Line 211 Line 209
       SUBROUTINE DGESVD( JOBU, JOBVT, M, N, A, LDA, S, U, LDU,        SUBROUTINE DGESVD( JOBU, JOBVT, M, N, A, LDA, S, U, LDU,
      $                   VT, LDVT, WORK, LWORK, INFO )       $                   VT, LDVT, WORK, LWORK, INFO )
 *  *
 *  -- LAPACK driver routine (version 3.4.0) --  *  -- LAPACK driver routine --
 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --  *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--  *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
 *     November 2011  
 *  *
 *     .. Scalar Arguments ..  *     .. Scalar Arguments ..
       CHARACTER          JOBU, JOBVT        CHARACTER          JOBU, JOBVT
Line 314 Line 311
             BDSPAC = 5*N              BDSPAC = 5*N
 *           Compute space needed for DGEQRF  *           Compute space needed for DGEQRF
             CALL DGEQRF( M, N, A, LDA, DUM(1), DUM(1), -1, IERR )              CALL DGEQRF( M, N, A, LDA, DUM(1), DUM(1), -1, IERR )
             LWORK_DGEQRF=DUM(1)              LWORK_DGEQRF = INT( DUM(1) )
 *           Compute space needed for DORGQR  *           Compute space needed for DORGQR
             CALL DORGQR( M, N, N, A, LDA, DUM(1), DUM(1), -1, IERR )              CALL DORGQR( M, N, N, A, LDA, DUM(1), DUM(1), -1, IERR )
             LWORK_DORGQR_N=DUM(1)              LWORK_DORGQR_N = INT( DUM(1) )
             CALL DORGQR( M, M, N, A, LDA, DUM(1), DUM(1), -1, IERR )              CALL DORGQR( M, M, N, A, LDA, DUM(1), DUM(1), -1, IERR )
             LWORK_DORGQR_M=DUM(1)              LWORK_DORGQR_M = INT( DUM(1) )
 *           Compute space needed for DGEBRD  *           Compute space needed for DGEBRD
             CALL DGEBRD( N, N, A, LDA, S, DUM(1), DUM(1),              CALL DGEBRD( N, N, A, LDA, S, DUM(1), DUM(1),
      $                   DUM(1), DUM(1), -1, IERR )       $                   DUM(1), DUM(1), -1, IERR )
             LWORK_DGEBRD=DUM(1)              LWORK_DGEBRD = INT( DUM(1) )
 *           Compute space needed for DORGBR P  *           Compute space needed for DORGBR P
             CALL DORGBR( 'P', N, N, N, A, LDA, DUM(1),              CALL DORGBR( 'P', N, N, N, A, LDA, DUM(1),
      $                   DUM(1), -1, IERR )       $                   DUM(1), -1, IERR )
             LWORK_DORGBR_P=DUM(1)              LWORK_DORGBR_P = INT( DUM(1) )
 *           Compute space needed for DORGBR Q  *           Compute space needed for DORGBR Q
             CALL DORGBR( 'Q', N, N, N, A, LDA, DUM(1),              CALL DORGBR( 'Q', N, N, N, A, LDA, DUM(1),
      $                   DUM(1), -1, IERR )       $                   DUM(1), -1, IERR )
             LWORK_DORGBR_Q=DUM(1)              LWORK_DORGBR_Q = INT( DUM(1) )
 *  *
             IF( M.GE.MNTHR ) THEN              IF( M.GE.MNTHR ) THEN
                IF( WNTUN ) THEN                 IF( WNTUN ) THEN
Line 339 Line 336
 *                 Path 1 (M much larger than N, JOBU='N')  *                 Path 1 (M much larger than N, JOBU='N')
 *  *
                   MAXWRK = N + LWORK_DGEQRF                    MAXWRK = N + LWORK_DGEQRF
                   MAXWRK = MAX( MAXWRK, 3*N+LWORK_DGEBRD )                    MAXWRK = MAX( MAXWRK, 3*N + LWORK_DGEBRD )
                   IF( WNTVO .OR. WNTVAS )                    IF( WNTVO .OR. WNTVAS )
      $               MAXWRK = MAX( MAXWRK, 3*N+LWORK_DORGBR_P )       $               MAXWRK = MAX( MAXWRK, 3*N + LWORK_DORGBR_P )
                   MAXWRK = MAX( MAXWRK, BDSPAC )                    MAXWRK = MAX( MAXWRK, BDSPAC )
                   MINWRK = MAX( 4*N, BDSPAC )                    MINWRK = MAX( 4*N, BDSPAC )
                ELSE IF( WNTUO .AND. WNTVN ) THEN                 ELSE IF( WNTUO .AND. WNTVN ) THEN
Line 349 Line 346
 *                 Path 2 (M much larger than N, JOBU='O', JOBVT='N')  *                 Path 2 (M much larger than N, JOBU='O', JOBVT='N')
 *  *
                   WRKBL = N + LWORK_DGEQRF                    WRKBL = N + LWORK_DGEQRF
                   WRKBL = MAX( WRKBL, N+LWORK_DORGQR_N )                    WRKBL = MAX( WRKBL, N + LWORK_DORGQR_N )
                   WRKBL = MAX( WRKBL, 3*N+LWORK_DGEBRD )                    WRKBL = MAX( WRKBL, 3*N + LWORK_DGEBRD )
                   WRKBL = MAX( WRKBL, 3*N+LWORK_DORGBR_Q )                    WRKBL = MAX( WRKBL, 3*N + LWORK_DORGBR_Q )
                   WRKBL = MAX( WRKBL, BDSPAC )                    WRKBL = MAX( WRKBL, BDSPAC )
                   MAXWRK = MAX( N*N+WRKBL, N*N+M*N+N )                    MAXWRK = MAX( N*N + WRKBL, N*N + M*N + N )
                   MINWRK = MAX( 3*N+M, BDSPAC )                    MINWRK = MAX( 3*N + M, BDSPAC )
                ELSE IF( WNTUO .AND. WNTVAS ) THEN                 ELSE IF( WNTUO .AND. WNTVAS ) THEN
 *  *
 *                 Path 3 (M much larger than N, JOBU='O', JOBVT='S' or  *                 Path 3 (M much larger than N, JOBU='O', JOBVT='S' or
 *                 'A')  *                 'A')
 *  *
                   WRKBL = N + LWORK_DGEQRF                    WRKBL = N + LWORK_DGEQRF
                   WRKBL = MAX( WRKBL, N+LWORK_DORGQR_N )                    WRKBL = MAX( WRKBL, N + LWORK_DORGQR_N )
                   WRKBL = MAX( WRKBL, 3*N+LWORK_DGEBRD )                    WRKBL = MAX( WRKBL, 3*N + LWORK_DGEBRD )
                   WRKBL = MAX( WRKBL, 3*N+LWORK_DORGBR_Q )                    WRKBL = MAX( WRKBL, 3*N + LWORK_DORGBR_Q )
                   WRKBL = MAX( WRKBL, 3*N+LWORK_DORGBR_P )                    WRKBL = MAX( WRKBL, 3*N + LWORK_DORGBR_P )
                   WRKBL = MAX( WRKBL, BDSPAC )                    WRKBL = MAX( WRKBL, BDSPAC )
                   MAXWRK = MAX( N*N+WRKBL, N*N+M*N+N )                    MAXWRK = MAX( N*N + WRKBL, N*N + M*N + N )
                   MINWRK = MAX( 3*N+M, BDSPAC )                    MINWRK = MAX( 3*N + M, BDSPAC )
                ELSE IF( WNTUS .AND. WNTVN ) THEN                 ELSE IF( WNTUS .AND. WNTVN ) THEN
 *  *
 *                 Path 4 (M much larger than N, JOBU='S', JOBVT='N')  *                 Path 4 (M much larger than N, JOBU='S', JOBVT='N')
 *  *
                   WRKBL = N + LWORK_DGEQRF                    WRKBL = N + LWORK_DGEQRF
                   WRKBL = MAX( WRKBL, N+LWORK_DORGQR_N )                    WRKBL = MAX( WRKBL, N + LWORK_DORGQR_N )
                   WRKBL = MAX( WRKBL, 3*N+LWORK_DGEBRD )                    WRKBL = MAX( WRKBL, 3*N + LWORK_DGEBRD )
                   WRKBL = MAX( WRKBL, 3*N+LWORK_DORGBR_Q )                    WRKBL = MAX( WRKBL, 3*N + LWORK_DORGBR_Q )
                   WRKBL = MAX( WRKBL, BDSPAC )                    WRKBL = MAX( WRKBL, BDSPAC )
                   MAXWRK = N*N + WRKBL                    MAXWRK = N*N + WRKBL
                   MINWRK = MAX( 3*N+M, BDSPAC )                    MINWRK = MAX( 3*N + M, BDSPAC )
                ELSE IF( WNTUS .AND. WNTVO ) THEN                 ELSE IF( WNTUS .AND. WNTVO ) THEN
 *  *
 *                 Path 5 (M much larger than N, JOBU='S', JOBVT='O')  *                 Path 5 (M much larger than N, JOBU='S', JOBVT='O')
 *  *
                   WRKBL = N + LWORK_DGEQRF                    WRKBL = N + LWORK_DGEQRF
                   WRKBL = MAX( WRKBL, N+LWORK_DORGQR_N )                    WRKBL = MAX( WRKBL, N + LWORK_DORGQR_N )
                   WRKBL = MAX( WRKBL, 3*N+LWORK_DGEBRD )                    WRKBL = MAX( WRKBL, 3*N + LWORK_DGEBRD )
                   WRKBL = MAX( WRKBL, 3*N+LWORK_DORGBR_Q )                    WRKBL = MAX( WRKBL, 3*N + LWORK_DORGBR_Q )
                   WRKBL = MAX( WRKBL, 3*N+LWORK_DORGBR_P )                    WRKBL = MAX( WRKBL, 3*N + LWORK_DORGBR_P )
                   WRKBL = MAX( WRKBL, BDSPAC )                    WRKBL = MAX( WRKBL, BDSPAC )
                   MAXWRK = 2*N*N + WRKBL                    MAXWRK = 2*N*N + WRKBL
                   MINWRK = MAX( 3*N+M, BDSPAC )                    MINWRK = MAX( 3*N + M, BDSPAC )
                ELSE IF( WNTUS .AND. WNTVAS ) THEN                 ELSE IF( WNTUS .AND. WNTVAS ) THEN
 *  *
 *                 Path 6 (M much larger than N, JOBU='S', JOBVT='S' or  *                 Path 6 (M much larger than N, JOBU='S', JOBVT='S' or
 *                 'A')  *                 'A')
 *  *
                   WRKBL = N + LWORK_DGEQRF                    WRKBL = N + LWORK_DGEQRF
                   WRKBL = MAX( WRKBL, N+LWORK_DORGQR_N )                    WRKBL = MAX( WRKBL, N + LWORK_DORGQR_N )
                   WRKBL = MAX( WRKBL, 3*N+LWORK_DGEBRD )                    WRKBL = MAX( WRKBL, 3*N + LWORK_DGEBRD )
                   WRKBL = MAX( WRKBL, 3*N+LWORK_DORGBR_Q )                    WRKBL = MAX( WRKBL, 3*N + LWORK_DORGBR_Q )
                   WRKBL = MAX( WRKBL, 3*N+LWORK_DORGBR_P )                    WRKBL = MAX( WRKBL, 3*N + LWORK_DORGBR_P )
                   WRKBL = MAX( WRKBL, BDSPAC )                    WRKBL = MAX( WRKBL, BDSPAC )
                   MAXWRK = N*N + WRKBL                    MAXWRK = N*N + WRKBL
                   MINWRK = MAX( 3*N+M, BDSPAC )                    MINWRK = MAX( 3*N + M, BDSPAC )
                ELSE IF( WNTUA .AND. WNTVN ) THEN                 ELSE IF( WNTUA .AND. WNTVN ) THEN
 *  *
 *                 Path 7 (M much larger than N, JOBU='A', JOBVT='N')  *                 Path 7 (M much larger than N, JOBU='A', JOBVT='N')
 *  *
                   WRKBL = N + LWORK_DGEQRF                    WRKBL = N + LWORK_DGEQRF
                   WRKBL = MAX( WRKBL, N+LWORK_DORGQR_M )                    WRKBL = MAX( WRKBL, N + LWORK_DORGQR_M )
                   WRKBL = MAX( WRKBL, 3*N+LWORK_DGEBRD )                    WRKBL = MAX( WRKBL, 3*N + LWORK_DGEBRD )
                   WRKBL = MAX( WRKBL, 3*N+LWORK_DORGBR_Q )                    WRKBL = MAX( WRKBL, 3*N + LWORK_DORGBR_Q )
                   WRKBL = MAX( WRKBL, BDSPAC )                    WRKBL = MAX( WRKBL, BDSPAC )
                   MAXWRK = N*N + WRKBL                    MAXWRK = N*N + WRKBL
                   MINWRK = MAX( 3*N+M, BDSPAC )                    MINWRK = MAX( 3*N + M, BDSPAC )
                ELSE IF( WNTUA .AND. WNTVO ) THEN                 ELSE IF( WNTUA .AND. WNTVO ) THEN
 *  *
 *                 Path 8 (M much larger than N, JOBU='A', JOBVT='O')  *                 Path 8 (M much larger than N, JOBU='A', JOBVT='O')
 *  *
                   WRKBL = N + LWORK_DGEQRF                    WRKBL = N + LWORK_DGEQRF
                   WRKBL = MAX( WRKBL, N+LWORK_DORGQR_M )                    WRKBL = MAX( WRKBL, N + LWORK_DORGQR_M )
                   WRKBL = MAX( WRKBL, 3*N+LWORK_DGEBRD )                    WRKBL = MAX( WRKBL, 3*N + LWORK_DGEBRD )
                   WRKBL = MAX( WRKBL, 3*N+LWORK_DORGBR_Q )                    WRKBL = MAX( WRKBL, 3*N + LWORK_DORGBR_Q )
                   WRKBL = MAX( WRKBL, 3*N+LWORK_DORGBR_P )                    WRKBL = MAX( WRKBL, 3*N + LWORK_DORGBR_P )
                   WRKBL = MAX( WRKBL, BDSPAC )                    WRKBL = MAX( WRKBL, BDSPAC )
                   MAXWRK = 2*N*N + WRKBL                    MAXWRK = 2*N*N + WRKBL
                   MINWRK = MAX( 3*N+M, BDSPAC )                    MINWRK = MAX( 3*N + M, BDSPAC )
                ELSE IF( WNTUA .AND. WNTVAS ) THEN                 ELSE IF( WNTUA .AND. WNTVAS ) THEN
 *  *
 *                 Path 9 (M much larger than N, JOBU='A', JOBVT='S' or  *                 Path 9 (M much larger than N, JOBU='A', JOBVT='S' or
 *                 'A')  *                 'A')
 *  *
                   WRKBL = N + LWORK_DGEQRF                    WRKBL = N + LWORK_DGEQRF
                   WRKBL = MAX( WRKBL, N+LWORK_DORGQR_M )                    WRKBL = MAX( WRKBL, N + LWORK_DORGQR_M )
                   WRKBL = MAX( WRKBL, 3*N+LWORK_DGEBRD )                    WRKBL = MAX( WRKBL, 3*N + LWORK_DGEBRD )
                   WRKBL = MAX( WRKBL, 3*N+LWORK_DORGBR_Q )                    WRKBL = MAX( WRKBL, 3*N + LWORK_DORGBR_Q )
                   WRKBL = MAX( WRKBL, 3*N+LWORK_DORGBR_P )                    WRKBL = MAX( WRKBL, 3*N + LWORK_DORGBR_P )
                   WRKBL = MAX( WRKBL, BDSPAC )                    WRKBL = MAX( WRKBL, BDSPAC )
                   MAXWRK = N*N + WRKBL                    MAXWRK = N*N + WRKBL
                   MINWRK = MAX( 3*N+M, BDSPAC )                    MINWRK = MAX( 3*N + M, BDSPAC )
                END IF                 END IF
             ELSE              ELSE
 *  *
Line 447 Line 444
 *  *
                CALL DGEBRD( M, N, A, LDA, S, DUM(1), DUM(1),                 CALL DGEBRD( M, N, A, LDA, S, DUM(1), DUM(1),
      $                   DUM(1), DUM(1), -1, IERR )       $                   DUM(1), DUM(1), -1, IERR )
                LWORK_DGEBRD=DUM(1)                 LWORK_DGEBRD = INT( DUM(1) )
                MAXWRK = 3*N + LWORK_DGEBRD                 MAXWRK = 3*N + LWORK_DGEBRD
                IF( WNTUS .OR. WNTUO ) THEN                 IF( WNTUS .OR. WNTUO ) THEN
                   CALL DORGBR( 'Q', M, N, N, A, LDA, DUM(1),                    CALL DORGBR( 'Q', M, N, N, A, LDA, DUM(1),
      $                   DUM(1), -1, IERR )       $                   DUM(1), -1, IERR )
                   LWORK_DORGBR_Q=DUM(1)                    LWORK_DORGBR_Q = INT( DUM(1) )
                   MAXWRK = MAX( MAXWRK, 3*N+LWORK_DORGBR_Q )                    MAXWRK = MAX( MAXWRK, 3*N + LWORK_DORGBR_Q )
                END IF                 END IF
                IF( WNTUA ) THEN                 IF( WNTUA ) THEN
                   CALL DORGBR( 'Q', M, M, N, A, LDA, DUM(1),                    CALL DORGBR( 'Q', M, M, N, A, LDA, DUM(1),
      $                   DUM(1), -1, IERR )       $                   DUM(1), -1, IERR )
                   LWORK_DORGBR_Q=DUM(1)                    LWORK_DORGBR_Q = INT( DUM(1) )
                   MAXWRK = MAX( MAXWRK, 3*N+LWORK_DORGBR_Q )                    MAXWRK = MAX( MAXWRK, 3*N + LWORK_DORGBR_Q )
                END IF                 END IF
                IF( .NOT.WNTVN ) THEN                 IF( .NOT.WNTVN ) THEN
                  MAXWRK = MAX( MAXWRK, 3*N+LWORK_DORGBR_P )                   MAXWRK = MAX( MAXWRK, 3*N + LWORK_DORGBR_P )
                END IF                 END IF
                MAXWRK = MAX( MAXWRK, BDSPAC )                 MAXWRK = MAX( MAXWRK, BDSPAC )
                MINWRK = MAX( 3*N+M, BDSPAC )                 MINWRK = MAX( 3*N + M, BDSPAC )
             END IF              END IF
          ELSE IF( MINMN.GT.0 ) THEN           ELSE IF( MINMN.GT.0 ) THEN
 *  *
Line 475 Line 472
             BDSPAC = 5*M              BDSPAC = 5*M
 *           Compute space needed for DGELQF  *           Compute space needed for DGELQF
             CALL DGELQF( M, N, A, LDA, DUM(1), DUM(1), -1, IERR )              CALL DGELQF( M, N, A, LDA, DUM(1), DUM(1), -1, IERR )
             LWORK_DGELQF=DUM(1)              LWORK_DGELQF = INT( DUM(1) )
 *           Compute space needed for DORGLQ  *           Compute space needed for DORGLQ
             CALL DORGLQ( N, N, M, VT, LDVT, DUM(1), DUM(1), -1, IERR )              CALL DORGLQ( N, N, M, DUM(1), N, DUM(1), DUM(1), -1, IERR )
             LWORK_DORGLQ_N=DUM(1)              LWORK_DORGLQ_N = INT( DUM(1) )
             CALL DORGLQ( M, N, M, A, LDA, DUM(1), DUM(1), -1, IERR )              CALL DORGLQ( M, N, M, A, LDA, DUM(1), DUM(1), -1, IERR )
             LWORK_DORGLQ_M=DUM(1)              LWORK_DORGLQ_M = INT( DUM(1) )
 *           Compute space needed for DGEBRD  *           Compute space needed for DGEBRD
             CALL DGEBRD( M, M, A, LDA, S, DUM(1), DUM(1),              CALL DGEBRD( M, M, A, LDA, S, DUM(1), DUM(1),
      $                   DUM(1), DUM(1), -1, IERR )       $                   DUM(1), DUM(1), -1, IERR )
             LWORK_DGEBRD=DUM(1)              LWORK_DGEBRD = INT( DUM(1) )
 *            Compute space needed for DORGBR P  *            Compute space needed for DORGBR P
             CALL DORGBR( 'P', M, M, M, A, N, DUM(1),              CALL DORGBR( 'P', M, M, M, A, N, DUM(1),
      $                   DUM(1), -1, IERR )       $                   DUM(1), -1, IERR )
             LWORK_DORGBR_P=DUM(1)              LWORK_DORGBR_P = INT( DUM(1) )
 *           Compute space needed for DORGBR Q  *           Compute space needed for DORGBR Q
             CALL DORGBR( 'Q', M, M, M, A, N, DUM(1),              CALL DORGBR( 'Q', M, M, M, A, N, DUM(1),
      $                   DUM(1), -1, IERR )       $                   DUM(1), -1, IERR )
             LWORK_DORGBR_Q=DUM(1)              LWORK_DORGBR_Q = INT( DUM(1) )
             IF( N.GE.MNTHR ) THEN              IF( N.GE.MNTHR ) THEN
                IF( WNTVN ) THEN                 IF( WNTVN ) THEN
 *  *
 *                 Path 1t(N much larger than M, JOBVT='N')  *                 Path 1t(N much larger than M, JOBVT='N')
 *  *
                   MAXWRK = M + LWORK_DGELQF                    MAXWRK = M + LWORK_DGELQF
                   MAXWRK = MAX( MAXWRK, 3*M+LWORK_DGEBRD )                    MAXWRK = MAX( MAXWRK, 3*M + LWORK_DGEBRD )
                   IF( WNTUO .OR. WNTUAS )                    IF( WNTUO .OR. WNTUAS )
      $               MAXWRK = MAX( MAXWRK, 3*M+LWORK_DORGBR_Q )       $               MAXWRK = MAX( MAXWRK, 3*M + LWORK_DORGBR_Q )
                   MAXWRK = MAX( MAXWRK, BDSPAC )                    MAXWRK = MAX( MAXWRK, BDSPAC )
                   MINWRK = MAX( 4*M, BDSPAC )                    MINWRK = MAX( 4*M, BDSPAC )
                ELSE IF( WNTVO .AND. WNTUN ) THEN                 ELSE IF( WNTVO .AND. WNTUN ) THEN
Line 509 Line 506
 *                 Path 2t(N much larger than M, JOBU='N', JOBVT='O')  *                 Path 2t(N much larger than M, JOBU='N', JOBVT='O')
 *  *
                   WRKBL = M + LWORK_DGELQF                    WRKBL = M + LWORK_DGELQF
                   WRKBL = MAX( WRKBL, M+LWORK_DORGLQ_M )                    WRKBL = MAX( WRKBL, M + LWORK_DORGLQ_M )
                   WRKBL = MAX( WRKBL, 3*M+LWORK_DGEBRD )                    WRKBL = MAX( WRKBL, 3*M + LWORK_DGEBRD )
                   WRKBL = MAX( WRKBL, 3*M+LWORK_DORGBR_P )                    WRKBL = MAX( WRKBL, 3*M + LWORK_DORGBR_P )
                   WRKBL = MAX( WRKBL, BDSPAC )                    WRKBL = MAX( WRKBL, BDSPAC )
                   MAXWRK = MAX( M*M+WRKBL, M*M+M*N+M )                    MAXWRK = MAX( M*M + WRKBL, M*M + M*N + M )
                   MINWRK = MAX( 3*M+N, BDSPAC )                    MINWRK = MAX( 3*M + N, BDSPAC )
                ELSE IF( WNTVO .AND. WNTUAS ) THEN                 ELSE IF( WNTVO .AND. WNTUAS ) THEN
 *  *
 *                 Path 3t(N much larger than M, JOBU='S' or 'A',  *                 Path 3t(N much larger than M, JOBU='S' or 'A',
 *                 JOBVT='O')  *                 JOBVT='O')
 *  *
                   WRKBL = M + LWORK_DGELQF                    WRKBL = M + LWORK_DGELQF
                   WRKBL = MAX( WRKBL, M+LWORK_DORGLQ_M )                    WRKBL = MAX( WRKBL, M + LWORK_DORGLQ_M )
                   WRKBL = MAX( WRKBL, 3*M+LWORK_DGEBRD )                    WRKBL = MAX( WRKBL, 3*M + LWORK_DGEBRD )
                   WRKBL = MAX( WRKBL, 3*M+LWORK_DORGBR_P )                    WRKBL = MAX( WRKBL, 3*M + LWORK_DORGBR_P )
                   WRKBL = MAX( WRKBL, 3*M+LWORK_DORGBR_Q )                    WRKBL = MAX( WRKBL, 3*M + LWORK_DORGBR_Q )
                   WRKBL = MAX( WRKBL, BDSPAC )                    WRKBL = MAX( WRKBL, BDSPAC )
                   MAXWRK = MAX( M*M+WRKBL, M*M+M*N+M )                    MAXWRK = MAX( M*M + WRKBL, M*M + M*N + M )
                   MINWRK = MAX( 3*M+N, BDSPAC )                    MINWRK = MAX( 3*M + N, BDSPAC )
                ELSE IF( WNTVS .AND. WNTUN ) THEN                 ELSE IF( WNTVS .AND. WNTUN ) THEN
 *  *
 *                 Path 4t(N much larger than M, JOBU='N', JOBVT='S')  *                 Path 4t(N much larger than M, JOBU='N', JOBVT='S')
 *  *
                   WRKBL = M + LWORK_DGELQF                    WRKBL = M + LWORK_DGELQF
                   WRKBL = MAX( WRKBL, M+LWORK_DORGLQ_M )                    WRKBL = MAX( WRKBL, M + LWORK_DORGLQ_M )
                   WRKBL = MAX( WRKBL, 3*M+LWORK_DGEBRD )                    WRKBL = MAX( WRKBL, 3*M + LWORK_DGEBRD )
                   WRKBL = MAX( WRKBL, 3*M+LWORK_DORGBR_P )                    WRKBL = MAX( WRKBL, 3*M + LWORK_DORGBR_P )
                   WRKBL = MAX( WRKBL, BDSPAC )                    WRKBL = MAX( WRKBL, BDSPAC )
                   MAXWRK = M*M + WRKBL                    MAXWRK = M*M + WRKBL
                   MINWRK = MAX( 3*M+N, BDSPAC )                    MINWRK = MAX( 3*M + N, BDSPAC )
                ELSE IF( WNTVS .AND. WNTUO ) THEN                 ELSE IF( WNTVS .AND. WNTUO ) THEN
 *  *
 *                 Path 5t(N much larger than M, JOBU='O', JOBVT='S')  *                 Path 5t(N much larger than M, JOBU='O', JOBVT='S')
 *  *
                   WRKBL = M + LWORK_DGELQF                    WRKBL = M + LWORK_DGELQF
                   WRKBL = MAX( WRKBL, M+LWORK_DORGLQ_M )                    WRKBL = MAX( WRKBL, M + LWORK_DORGLQ_M )
                   WRKBL = MAX( WRKBL, 3*M+LWORK_DGEBRD )                    WRKBL = MAX( WRKBL, 3*M + LWORK_DGEBRD )
                   WRKBL = MAX( WRKBL, 3*M+LWORK_DORGBR_P )                    WRKBL = MAX( WRKBL, 3*M + LWORK_DORGBR_P )
                   WRKBL = MAX( WRKBL, 3*M+LWORK_DORGBR_Q )                    WRKBL = MAX( WRKBL, 3*M + LWORK_DORGBR_Q )
                   WRKBL = MAX( WRKBL, BDSPAC )                    WRKBL = MAX( WRKBL, BDSPAC )
                   MAXWRK = 2*M*M + WRKBL                    MAXWRK = 2*M*M + WRKBL
                   MINWRK = MAX( 3*M+N, BDSPAC )                    MINWRK = MAX( 3*M + N, BDSPAC )
                ELSE IF( WNTVS .AND. WNTUAS ) THEN                 ELSE IF( WNTVS .AND. WNTUAS ) THEN
 *  *
 *                 Path 6t(N much larger than M, JOBU='S' or 'A',  *                 Path 6t(N much larger than M, JOBU='S' or 'A',
 *                 JOBVT='S')  *                 JOBVT='S')
 *  *
                   WRKBL = M + LWORK_DGELQF                    WRKBL = M + LWORK_DGELQF
                   WRKBL = MAX( WRKBL, M+LWORK_DORGLQ_M )                    WRKBL = MAX( WRKBL, M + LWORK_DORGLQ_M )
                   WRKBL = MAX( WRKBL, 3*M+LWORK_DGEBRD )                    WRKBL = MAX( WRKBL, 3*M + LWORK_DGEBRD )
                   WRKBL = MAX( WRKBL, 3*M+LWORK_DORGBR_P )                    WRKBL = MAX( WRKBL, 3*M + LWORK_DORGBR_P )
                   WRKBL = MAX( WRKBL, 3*M+LWORK_DORGBR_Q )                    WRKBL = MAX( WRKBL, 3*M + LWORK_DORGBR_Q )
                   WRKBL = MAX( WRKBL, BDSPAC )                    WRKBL = MAX( WRKBL, BDSPAC )
                   MAXWRK = M*M + WRKBL                    MAXWRK = M*M + WRKBL
                   MINWRK = MAX( 3*M+N, BDSPAC )                    MINWRK = MAX( 3*M + N, BDSPAC )
                ELSE IF( WNTVA .AND. WNTUN ) THEN                 ELSE IF( WNTVA .AND. WNTUN ) THEN
 *  *
 *                 Path 7t(N much larger than M, JOBU='N', JOBVT='A')  *                 Path 7t(N much larger than M, JOBU='N', JOBVT='A')
 *  *
                   WRKBL = M + LWORK_DGELQF                    WRKBL = M + LWORK_DGELQF
                   WRKBL = MAX( WRKBL, M+LWORK_DORGLQ_N )                    WRKBL = MAX( WRKBL, M + LWORK_DORGLQ_N )
                   WRKBL = MAX( WRKBL, 3*M+LWORK_DGEBRD )                    WRKBL = MAX( WRKBL, 3*M + LWORK_DGEBRD )
                   WRKBL = MAX( WRKBL, 3*M+LWORK_DORGBR_P )                    WRKBL = MAX( WRKBL, 3*M + LWORK_DORGBR_P )
                   WRKBL = MAX( WRKBL, BDSPAC )                    WRKBL = MAX( WRKBL, BDSPAC )
                   MAXWRK = M*M + WRKBL                    MAXWRK = M*M + WRKBL
                   MINWRK = MAX( 3*M+N, BDSPAC )                    MINWRK = MAX( 3*M + N, BDSPAC )
                ELSE IF( WNTVA .AND. WNTUO ) THEN                 ELSE IF( WNTVA .AND. WNTUO ) THEN
 *  *
 *                 Path 8t(N much larger than M, JOBU='O', JOBVT='A')  *                 Path 8t(N much larger than M, JOBU='O', JOBVT='A')
 *  *
                   WRKBL = M + LWORK_DGELQF                    WRKBL = M + LWORK_DGELQF
                   WRKBL = MAX( WRKBL, M+LWORK_DORGLQ_N )                    WRKBL = MAX( WRKBL, M + LWORK_DORGLQ_N )
                   WRKBL = MAX( WRKBL, 3*M+LWORK_DGEBRD )                    WRKBL = MAX( WRKBL, 3*M + LWORK_DGEBRD )
                   WRKBL = MAX( WRKBL, 3*M+LWORK_DORGBR_P )                    WRKBL = MAX( WRKBL, 3*M + LWORK_DORGBR_P )
                   WRKBL = MAX( WRKBL, 3*M+LWORK_DORGBR_Q )                    WRKBL = MAX( WRKBL, 3*M + LWORK_DORGBR_Q )
                   WRKBL = MAX( WRKBL, BDSPAC )                    WRKBL = MAX( WRKBL, BDSPAC )
                   MAXWRK = 2*M*M + WRKBL                    MAXWRK = 2*M*M + WRKBL
                   MINWRK = MAX( 3*M+N, BDSPAC )                    MINWRK = MAX( 3*M + N, BDSPAC )
                ELSE IF( WNTVA .AND. WNTUAS ) THEN                 ELSE IF( WNTVA .AND. WNTUAS ) THEN
 *  *
 *                 Path 9t(N much larger than M, JOBU='S' or 'A',  *                 Path 9t(N much larger than M, JOBU='S' or 'A',
 *                 JOBVT='A')  *                 JOBVT='A')
 *  *
                   WRKBL = M + LWORK_DGELQF                    WRKBL = M + LWORK_DGELQF
                   WRKBL = MAX( WRKBL, M+LWORK_DORGLQ_N )                    WRKBL = MAX( WRKBL, M + LWORK_DORGLQ_N )
                   WRKBL = MAX( WRKBL, 3*M+LWORK_DGEBRD )                    WRKBL = MAX( WRKBL, 3*M + LWORK_DGEBRD )
                   WRKBL = MAX( WRKBL, 3*M+LWORK_DORGBR_P )                    WRKBL = MAX( WRKBL, 3*M + LWORK_DORGBR_P )
                   WRKBL = MAX( WRKBL, 3*M+LWORK_DORGBR_Q )                    WRKBL = MAX( WRKBL, 3*M + LWORK_DORGBR_Q )
                   WRKBL = MAX( WRKBL, BDSPAC )                    WRKBL = MAX( WRKBL, BDSPAC )
                   MAXWRK = M*M + WRKBL                    MAXWRK = M*M + WRKBL
                   MINWRK = MAX( 3*M+N, BDSPAC )                    MINWRK = MAX( 3*M + N, BDSPAC )
                END IF                 END IF
             ELSE              ELSE
 *  *
Line 607 Line 604
 *  *
                CALL DGEBRD( M, N, A, LDA, S, DUM(1), DUM(1),                 CALL DGEBRD( M, N, A, LDA, S, DUM(1), DUM(1),
      $                   DUM(1), DUM(1), -1, IERR )       $                   DUM(1), DUM(1), -1, IERR )
                LWORK_DGEBRD=DUM(1)                 LWORK_DGEBRD = INT( DUM(1) )
                MAXWRK = 3*M + LWORK_DGEBRD                 MAXWRK = 3*M + LWORK_DGEBRD
                IF( WNTVS .OR. WNTVO ) THEN                 IF( WNTVS .OR. WNTVO ) THEN
 *                Compute space needed for DORGBR P  *                Compute space needed for DORGBR P
                  CALL DORGBR( 'P', M, N, M, A, N, DUM(1),                   CALL DORGBR( 'P', M, N, M, A, N, DUM(1),
      $                   DUM(1), -1, IERR )       $                   DUM(1), -1, IERR )
                  LWORK_DORGBR_P=DUM(1)                   LWORK_DORGBR_P = INT( DUM(1) )
                  MAXWRK = MAX( MAXWRK, 3*M+LWORK_DORGBR_P )                   MAXWRK = MAX( MAXWRK, 3*M + LWORK_DORGBR_P )
                END IF                 END IF
                IF( WNTVA ) THEN                 IF( WNTVA ) THEN
                  CALL DORGBR( 'P', N, N, M, A, N, DUM(1),                   CALL DORGBR( 'P', N, N, M, A, N, DUM(1),
      $                   DUM(1), -1, IERR )       $                   DUM(1), -1, IERR )
                  LWORK_DORGBR_P=DUM(1)                   LWORK_DORGBR_P = INT( DUM(1) )
                  MAXWRK = MAX( MAXWRK, 3*M+LWORK_DORGBR_P )                   MAXWRK = MAX( MAXWRK, 3*M + LWORK_DORGBR_P )
                END IF                 END IF
                IF( .NOT.WNTUN ) THEN                 IF( .NOT.WNTUN ) THEN
                   MAXWRK = MAX( MAXWRK, 3*M+LWORK_DORGBR_Q )                    MAXWRK = MAX( MAXWRK, 3*M + LWORK_DORGBR_Q )
                END IF                 END IF
                MAXWRK = MAX( MAXWRK, BDSPAC )                 MAXWRK = MAX( MAXWRK, BDSPAC )
                MINWRK = MAX( 3*M+N, BDSPAC )                 MINWRK = MAX( 3*M + N, BDSPAC )
             END IF              END IF
          END IF           END IF
          MAXWRK = MAX( MAXWRK, MINWRK )           MAXWRK = MAX( MAXWRK, MINWRK )
Line 685 Line 682
                IWORK = ITAU + N                 IWORK = ITAU + N
 *  *
 *              Compute A=Q*R  *              Compute A=Q*R
 *              (Workspace: need 2*N, prefer N+N*NB)  *              (Workspace: need 2*N, prefer N + N*NB)
 *  *
                CALL DGEQRF( M, N, A, LDA, WORK( ITAU ), WORK( IWORK ),                 CALL DGEQRF( M, N, A, LDA, WORK( ITAU ), WORK( IWORK ),
      $                      LWORK-IWORK+1, IERR )       $                      LWORK-IWORK+1, IERR )
 *  *
 *              Zero out below R  *              Zero out below R
 *  *
                CALL DLASET( 'L', N-1, N-1, ZERO, ZERO, A( 2, 1 ), LDA )                 IF( N .GT. 1 ) THEN
                     CALL DLASET( 'L', N-1, N-1, ZERO, ZERO, A( 2, 1 ),
        $                         LDA )
                  END IF
                IE = 1                 IE = 1
                ITAUQ = IE + N                 ITAUQ = IE + N
                ITAUP = ITAUQ + N                 ITAUP = ITAUQ + N
                IWORK = ITAUP + N                 IWORK = ITAUP + N
 *  *
 *              Bidiagonalize R in A  *              Bidiagonalize R in A
 *              (Workspace: need 4*N, prefer 3*N+2*N*NB)  *              (Workspace: need 4*N, prefer 3*N + 2*N*NB)
 *  *
                CALL DGEBRD( N, N, A, LDA, S, WORK( IE ), WORK( ITAUQ ),                 CALL DGEBRD( N, N, A, LDA, S, WORK( IE ), WORK( ITAUQ ),
      $                      WORK( ITAUP ), WORK( IWORK ), LWORK-IWORK+1,       $                      WORK( ITAUP ), WORK( IWORK ), LWORK-IWORK+1,
Line 708 Line 708
                IF( WNTVO .OR. WNTVAS ) THEN                 IF( WNTVO .OR. WNTVAS ) THEN
 *  *
 *                 If right singular vectors desired, generate P'.  *                 If right singular vectors desired, generate P'.
 *                 (Workspace: need 4*N-1, prefer 3*N+(N-1)*NB)  *                 (Workspace: need 4*N-1, prefer 3*N + (N-1)*NB)
 *  *
                   CALL DORGBR( 'P', N, N, N, A, LDA, WORK( ITAUP ),                    CALL DORGBR( 'P', N, N, N, A, LDA, WORK( ITAUP ),
      $                         WORK( IWORK ), LWORK-IWORK+1, IERR )       $                         WORK( IWORK ), LWORK-IWORK+1, IERR )
Line 739 Line 739
 *                 Sufficient workspace for a fast algorithm  *                 Sufficient workspace for a fast algorithm
 *  *
                   IR = 1                    IR = 1
                   IF( LWORK.GE.MAX( WRKBL, LDA*N+N )+LDA*N ) THEN                    IF( LWORK.GE.MAX( WRKBL, LDA*N + N ) + LDA*N ) THEN
 *  *
 *                    WORK(IU) is LDA by N, WORK(IR) is LDA by N  *                    WORK(IU) is LDA by N, WORK(IR) is LDA by N
 *  *
                      LDWRKU = LDA                       LDWRKU = LDA
                      LDWRKR = LDA                       LDWRKR = LDA
                   ELSE IF( LWORK.GE.MAX( WRKBL, LDA*N+N )+N*N ) THEN                    ELSE IF( LWORK.GE.MAX( WRKBL, LDA*N + N ) + N*N ) THEN
 *  *
 *                    WORK(IU) is LDA by N, WORK(IR) is N by N  *                    WORK(IU) is LDA by N, WORK(IR) is N by N
 *  *
Line 762 Line 762
                   IWORK = ITAU + N                    IWORK = ITAU + N
 *  *
 *                 Compute A=Q*R  *                 Compute A=Q*R
 *                 (Workspace: need N*N+2*N, prefer N*N+N+N*NB)  *                 (Workspace: need N*N + 2*N, prefer N*N + N + N*NB)
 *  *
                   CALL DGEQRF( M, N, A, LDA, WORK( ITAU ),                    CALL DGEQRF( M, N, A, LDA, WORK( ITAU ),
      $                         WORK( IWORK ), LWORK-IWORK+1, IERR )       $                         WORK( IWORK ), LWORK-IWORK+1, IERR )
Line 774 Line 774
      $                         LDWRKR )       $                         LDWRKR )
 *  *
 *                 Generate Q in A  *                 Generate Q in A
 *                 (Workspace: need N*N+2*N, prefer N*N+N+N*NB)  *                 (Workspace: need N*N + 2*N, prefer N*N + N + N*NB)
 *  *
                   CALL DORGQR( M, N, N, A, LDA, WORK( ITAU ),                    CALL DORGQR( M, N, N, A, LDA, WORK( ITAU ),
      $                         WORK( IWORK ), LWORK-IWORK+1, IERR )       $                         WORK( IWORK ), LWORK-IWORK+1, IERR )
Line 784 Line 784
                   IWORK = ITAUP + N                    IWORK = ITAUP + N
 *  *
 *                 Bidiagonalize R in WORK(IR)  *                 Bidiagonalize R in WORK(IR)
 *                 (Workspace: need N*N+4*N, prefer N*N+3*N+2*N*NB)  *                 (Workspace: need N*N + 4*N, prefer N*N + 3*N + 2*N*NB)
 *  *
                   CALL DGEBRD( N, N, WORK( IR ), LDWRKR, S, WORK( IE ),                    CALL DGEBRD( N, N, WORK( IR ), LDWRKR, S, WORK( IE ),
      $                         WORK( ITAUQ ), WORK( ITAUP ),       $                         WORK( ITAUQ ), WORK( ITAUP ),
      $                         WORK( IWORK ), LWORK-IWORK+1, IERR )       $                         WORK( IWORK ), LWORK-IWORK+1, IERR )
 *  *
 *                 Generate left vectors bidiagonalizing R  *                 Generate left vectors bidiagonalizing R
 *                 (Workspace: need N*N+4*N, prefer N*N+3*N+N*NB)  *                 (Workspace: need N*N + 4*N, prefer N*N + 3*N + N*NB)
 *  *
                   CALL DORGBR( 'Q', N, N, N, WORK( IR ), LDWRKR,                    CALL DORGBR( 'Q', N, N, N, WORK( IR ), LDWRKR,
      $                         WORK( ITAUQ ), WORK( IWORK ),       $                         WORK( ITAUQ ), WORK( IWORK ),
Line 800 Line 800
 *  *
 *                 Perform bidiagonal QR iteration, computing left  *                 Perform bidiagonal QR iteration, computing left
 *                 singular vectors of R in WORK(IR)  *                 singular vectors of R in WORK(IR)
 *                 (Workspace: need N*N+BDSPAC)  *                 (Workspace: need N*N + BDSPAC)
 *  *
                   CALL DBDSQR( 'U', N, 0, N, 0, S, WORK( IE ), DUM, 1,                    CALL DBDSQR( 'U', N, 0, N, 0, S, WORK( IE ), DUM, 1,
      $                         WORK( IR ), LDWRKR, DUM, 1,       $                         WORK( IR ), LDWRKR, DUM, 1,
Line 809 Line 809
 *  *
 *                 Multiply Q in A by left singular vectors of R in  *                 Multiply Q in A by left singular vectors of R in
 *                 WORK(IR), storing result in WORK(IU) and copying to A  *                 WORK(IR), storing result in WORK(IU) and copying to A
 *                 (Workspace: need N*N+2*N, prefer N*N+M*N+N)  *                 (Workspace: need N*N + 2*N, prefer N*N + M*N + N)
 *  *
                   DO 10 I = 1, M, LDWRKU                    DO 10 I = 1, M, LDWRKU
                      CHUNK = MIN( M-I+1, LDWRKU )                       CHUNK = MIN( M-I+1, LDWRKU )
Line 830 Line 830
                   IWORK = ITAUP + N                    IWORK = ITAUP + N
 *  *
 *                 Bidiagonalize A  *                 Bidiagonalize A
 *                 (Workspace: need 3*N+M, prefer 3*N+(M+N)*NB)  *                 (Workspace: need 3*N + M, prefer 3*N + (M + N)*NB)
 *  *
                   CALL DGEBRD( M, N, A, LDA, S, WORK( IE ),                    CALL DGEBRD( M, N, A, LDA, S, WORK( IE ),
      $                         WORK( ITAUQ ), WORK( ITAUP ),       $                         WORK( ITAUQ ), WORK( ITAUP ),
      $                         WORK( IWORK ), LWORK-IWORK+1, IERR )       $                         WORK( IWORK ), LWORK-IWORK+1, IERR )
 *  *
 *                 Generate left vectors bidiagonalizing A  *                 Generate left vectors bidiagonalizing A
 *                 (Workspace: need 4*N, prefer 3*N+N*NB)  *                 (Workspace: need 4*N, prefer 3*N + N*NB)
 *  *
                   CALL DORGBR( 'Q', M, N, N, A, LDA, WORK( ITAUQ ),                    CALL DORGBR( 'Q', M, N, N, A, LDA, WORK( ITAUQ ),
      $                         WORK( IWORK ), LWORK-IWORK+1, IERR )       $                         WORK( IWORK ), LWORK-IWORK+1, IERR )
Line 863 Line 863
 *                 Sufficient workspace for a fast algorithm  *                 Sufficient workspace for a fast algorithm
 *  *
                   IR = 1                    IR = 1
                   IF( LWORK.GE.MAX( WRKBL, LDA*N+N )+LDA*N ) THEN                    IF( LWORK.GE.MAX( WRKBL, LDA*N + N ) + LDA*N ) THEN
 *  *
 *                    WORK(IU) is LDA by N and WORK(IR) is LDA by N  *                    WORK(IU) is LDA by N and WORK(IR) is LDA by N
 *  *
                      LDWRKU = LDA                       LDWRKU = LDA
                      LDWRKR = LDA                       LDWRKR = LDA
                   ELSE IF( LWORK.GE.MAX( WRKBL, LDA*N+N )+N*N ) THEN                    ELSE IF( LWORK.GE.MAX( WRKBL, LDA*N + N ) + N*N ) THEN
 *  *
 *                    WORK(IU) is LDA by N and WORK(IR) is N by N  *                    WORK(IU) is LDA by N and WORK(IR) is N by N
 *  *
Line 886 Line 886
                   IWORK = ITAU + N                    IWORK = ITAU + N
 *  *
 *                 Compute A=Q*R  *                 Compute A=Q*R
 *                 (Workspace: need N*N+2*N, prefer N*N+N+N*NB)  *                 (Workspace: need N*N + 2*N, prefer N*N + N + N*NB)
 *  *
                   CALL DGEQRF( M, N, A, LDA, WORK( ITAU ),                    CALL DGEQRF( M, N, A, LDA, WORK( ITAU ),
      $                         WORK( IWORK ), LWORK-IWORK+1, IERR )       $                         WORK( IWORK ), LWORK-IWORK+1, IERR )
Line 899 Line 899
      $                            VT( 2, 1 ), LDVT )       $                            VT( 2, 1 ), LDVT )
 *  *
 *                 Generate Q in A  *                 Generate Q in A
 *                 (Workspace: need N*N+2*N, prefer N*N+N+N*NB)  *                 (Workspace: need N*N + 2*N, prefer N*N + N + N*NB)
 *  *
                   CALL DORGQR( M, N, N, A, LDA, WORK( ITAU ),                    CALL DORGQR( M, N, N, A, LDA, WORK( ITAU ),
      $                         WORK( IWORK ), LWORK-IWORK+1, IERR )       $                         WORK( IWORK ), LWORK-IWORK+1, IERR )
Line 909 Line 909
                   IWORK = ITAUP + N                    IWORK = ITAUP + N
 *  *
 *                 Bidiagonalize R in VT, copying result to WORK(IR)  *                 Bidiagonalize R in VT, copying result to WORK(IR)
 *                 (Workspace: need N*N+4*N, prefer N*N+3*N+2*N*NB)  *                 (Workspace: need N*N + 4*N, prefer N*N + 3*N + 2*N*NB)
 *  *
                   CALL DGEBRD( N, N, VT, LDVT, S, WORK( IE ),                    CALL DGEBRD( N, N, VT, LDVT, S, WORK( IE ),
      $                         WORK( ITAUQ ), WORK( ITAUP ),       $                         WORK( ITAUQ ), WORK( ITAUP ),
Line 917 Line 917
                   CALL DLACPY( 'L', N, N, VT, LDVT, WORK( IR ), LDWRKR )                    CALL DLACPY( 'L', N, N, VT, LDVT, WORK( IR ), LDWRKR )
 *  *
 *                 Generate left vectors bidiagonalizing R in WORK(IR)  *                 Generate left vectors bidiagonalizing R in WORK(IR)
 *                 (Workspace: need N*N+4*N, prefer N*N+3*N+N*NB)  *                 (Workspace: need N*N + 4*N, prefer N*N + 3*N + N*NB)
 *  *
                   CALL DORGBR( 'Q', N, N, N, WORK( IR ), LDWRKR,                    CALL DORGBR( 'Q', N, N, N, WORK( IR ), LDWRKR,
      $                         WORK( ITAUQ ), WORK( IWORK ),       $                         WORK( ITAUQ ), WORK( IWORK ),
      $                         LWORK-IWORK+1, IERR )       $                         LWORK-IWORK+1, IERR )
 *  *
 *                 Generate right vectors bidiagonalizing R in VT  *                 Generate right vectors bidiagonalizing R in VT
 *                 (Workspace: need N*N+4*N-1, prefer N*N+3*N+(N-1)*NB)  *                 (Workspace: need N*N + 4*N-1, prefer N*N + 3*N + (N-1)*NB)
 *  *
                   CALL DORGBR( 'P', N, N, N, VT, LDVT, WORK( ITAUP ),                    CALL DORGBR( 'P', N, N, N, VT, LDVT, WORK( ITAUP ),
      $                         WORK( IWORK ), LWORK-IWORK+1, IERR )       $                         WORK( IWORK ), LWORK-IWORK+1, IERR )
Line 933 Line 933
 *                 Perform bidiagonal QR iteration, computing left  *                 Perform bidiagonal QR iteration, computing left
 *                 singular vectors of R in WORK(IR) and computing right  *                 singular vectors of R in WORK(IR) and computing right
 *                 singular vectors of R in VT  *                 singular vectors of R in VT
 *                 (Workspace: need N*N+BDSPAC)  *                 (Workspace: need N*N + BDSPAC)
 *  *
                   CALL DBDSQR( 'U', N, N, N, 0, S, WORK( IE ), VT, LDVT,                    CALL DBDSQR( 'U', N, N, N, 0, S, WORK( IE ), VT, LDVT,
      $                         WORK( IR ), LDWRKR, DUM, 1,       $                         WORK( IR ), LDWRKR, DUM, 1,
Line 942 Line 942
 *  *
 *                 Multiply Q in A by left singular vectors of R in  *                 Multiply Q in A by left singular vectors of R in
 *                 WORK(IR), storing result in WORK(IU) and copying to A  *                 WORK(IR), storing result in WORK(IU) and copying to A
 *                 (Workspace: need N*N+2*N, prefer N*N+M*N+N)  *                 (Workspace: need N*N + 2*N, prefer N*N + M*N + N)
 *  *
                   DO 20 I = 1, M, LDWRKU                    DO 20 I = 1, M, LDWRKU
                      CHUNK = MIN( M-I+1, LDWRKU )                       CHUNK = MIN( M-I+1, LDWRKU )
Line 961 Line 961
                   IWORK = ITAU + N                    IWORK = ITAU + N
 *  *
 *                 Compute A=Q*R  *                 Compute A=Q*R
 *                 (Workspace: need 2*N, prefer N+N*NB)  *                 (Workspace: need 2*N, prefer N + N*NB)
 *  *
                   CALL DGEQRF( M, N, A, LDA, WORK( ITAU ),                    CALL DGEQRF( M, N, A, LDA, WORK( ITAU ),
      $                         WORK( IWORK ), LWORK-IWORK+1, IERR )       $                         WORK( IWORK ), LWORK-IWORK+1, IERR )
Line 974 Line 974
      $                            VT( 2, 1 ), LDVT )       $                            VT( 2, 1 ), LDVT )
 *  *
 *                 Generate Q in A  *                 Generate Q in A
 *                 (Workspace: need 2*N, prefer N+N*NB)  *                 (Workspace: need 2*N, prefer N + N*NB)
 *  *
                   CALL DORGQR( M, N, N, A, LDA, WORK( ITAU ),                    CALL DORGQR( M, N, N, A, LDA, WORK( ITAU ),
      $                         WORK( IWORK ), LWORK-IWORK+1, IERR )       $                         WORK( IWORK ), LWORK-IWORK+1, IERR )
Line 984 Line 984
                   IWORK = ITAUP + N                    IWORK = ITAUP + N
 *  *
 *                 Bidiagonalize R in VT  *                 Bidiagonalize R in VT
 *                 (Workspace: need 4*N, prefer 3*N+2*N*NB)  *                 (Workspace: need 4*N, prefer 3*N + 2*N*NB)
 *  *
                   CALL DGEBRD( N, N, VT, LDVT, S, WORK( IE ),                    CALL DGEBRD( N, N, VT, LDVT, S, WORK( IE ),
      $                         WORK( ITAUQ ), WORK( ITAUP ),       $                         WORK( ITAUQ ), WORK( ITAUP ),
      $                         WORK( IWORK ), LWORK-IWORK+1, IERR )       $                         WORK( IWORK ), LWORK-IWORK+1, IERR )
 *  *
 *                 Multiply Q in A by left vectors bidiagonalizing R  *                 Multiply Q in A by left vectors bidiagonalizing R
 *                 (Workspace: need 3*N+M, prefer 3*N+M*NB)  *                 (Workspace: need 3*N + M, prefer 3*N + M*NB)
 *  *
                   CALL DORMBR( 'Q', 'R', 'N', M, N, N, VT, LDVT,                    CALL DORMBR( 'Q', 'R', 'N', M, N, N, VT, LDVT,
      $                         WORK( ITAUQ ), A, LDA, WORK( IWORK ),       $                         WORK( ITAUQ ), A, LDA, WORK( IWORK ),
      $                         LWORK-IWORK+1, IERR )       $                         LWORK-IWORK+1, IERR )
 *  *
 *                 Generate right vectors bidiagonalizing R in VT  *                 Generate right vectors bidiagonalizing R in VT
 *                 (Workspace: need 4*N-1, prefer 3*N+(N-1)*NB)  *                 (Workspace: need 4*N-1, prefer 3*N + (N-1)*NB)
 *  *
                   CALL DORGBR( 'P', N, N, N, VT, LDVT, WORK( ITAUP ),                    CALL DORGBR( 'P', N, N, N, VT, LDVT, WORK( ITAUP ),
      $                         WORK( IWORK ), LWORK-IWORK+1, IERR )       $                         WORK( IWORK ), LWORK-IWORK+1, IERR )
Line 1042 Line 1042
                      IWORK = ITAU + N                       IWORK = ITAU + N
 *  *
 *                    Compute A=Q*R  *                    Compute A=Q*R
 *                    (Workspace: need N*N+2*N, prefer N*N+N+N*NB)  *                    (Workspace: need N*N + 2*N, prefer N*N + N + N*NB)
 *  *
                      CALL DGEQRF( M, N, A, LDA, WORK( ITAU ),                       CALL DGEQRF( M, N, A, LDA, WORK( ITAU ),
      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )       $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
Line 1055 Line 1055
      $                            WORK( IR+1 ), LDWRKR )       $                            WORK( IR+1 ), LDWRKR )
 *  *
 *                    Generate Q in A  *                    Generate Q in A
 *                    (Workspace: need N*N+2*N, prefer N*N+N+N*NB)  *                    (Workspace: need N*N + 2*N, prefer N*N + N + N*NB)
 *  *
                      CALL DORGQR( M, N, N, A, LDA, WORK( ITAU ),                       CALL DORGQR( M, N, N, A, LDA, WORK( ITAU ),
      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )       $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
Line 1065 Line 1065
                      IWORK = ITAUP + N                       IWORK = ITAUP + N
 *  *
 *                    Bidiagonalize R in WORK(IR)  *                    Bidiagonalize R in WORK(IR)
 *                    (Workspace: need N*N+4*N, prefer N*N+3*N+2*N*NB)  *                    (Workspace: need N*N + 4*N, prefer N*N + 3*N + 2*N*NB)
 *  *
                      CALL DGEBRD( N, N, WORK( IR ), LDWRKR, S,                       CALL DGEBRD( N, N, WORK( IR ), LDWRKR, S,
      $                            WORK( IE ), WORK( ITAUQ ),       $                            WORK( IE ), WORK( ITAUQ ),
Line 1073 Line 1073
      $                            LWORK-IWORK+1, IERR )       $                            LWORK-IWORK+1, IERR )
 *  *
 *                    Generate left vectors bidiagonalizing R in WORK(IR)  *                    Generate left vectors bidiagonalizing R in WORK(IR)
 *                    (Workspace: need N*N+4*N, prefer N*N+3*N+N*NB)  *                    (Workspace: need N*N + 4*N, prefer N*N + 3*N + N*NB)
 *  *
                      CALL DORGBR( 'Q', N, N, N, WORK( IR ), LDWRKR,                       CALL DORGBR( 'Q', N, N, N, WORK( IR ), LDWRKR,
      $                            WORK( ITAUQ ), WORK( IWORK ),       $                            WORK( ITAUQ ), WORK( IWORK ),
Line 1082 Line 1082
 *  *
 *                    Perform bidiagonal QR iteration, computing left  *                    Perform bidiagonal QR iteration, computing left
 *                    singular vectors of R in WORK(IR)  *                    singular vectors of R in WORK(IR)
 *                    (Workspace: need N*N+BDSPAC)  *                    (Workspace: need N*N + BDSPAC)
 *  *
                      CALL DBDSQR( 'U', N, 0, N, 0, S, WORK( IE ), DUM,                       CALL DBDSQR( 'U', N, 0, N, 0, S, WORK( IE ), DUM,
      $                            1, WORK( IR ), LDWRKR, DUM, 1,       $                            1, WORK( IR ), LDWRKR, DUM, 1,
Line 1103 Line 1103
                      IWORK = ITAU + N                       IWORK = ITAU + N
 *  *
 *                    Compute A=Q*R, copying result to U  *                    Compute A=Q*R, copying result to U
 *                    (Workspace: need 2*N, prefer N+N*NB)  *                    (Workspace: need 2*N, prefer N + N*NB)
 *  *
                      CALL DGEQRF( M, N, A, LDA, WORK( ITAU ),                       CALL DGEQRF( M, N, A, LDA, WORK( ITAU ),
      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )       $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
                      CALL DLACPY( 'L', M, N, A, LDA, U, LDU )                       CALL DLACPY( 'L', M, N, A, LDA, U, LDU )
 *  *
 *                    Generate Q in U  *                    Generate Q in U
 *                    (Workspace: need 2*N, prefer N+N*NB)  *                    (Workspace: need 2*N, prefer N + N*NB)
 *  *
                      CALL DORGQR( M, N, N, U, LDU, WORK( ITAU ),                       CALL DORGQR( M, N, N, U, LDU, WORK( ITAU ),
      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )       $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
Line 1121 Line 1121
 *  *
 *                    Zero out below R in A  *                    Zero out below R in A
 *  *
                      CALL DLASET( 'L', N-1, N-1, ZERO, ZERO, A( 2, 1 ),                       IF( N .GT. 1 ) THEN
      $                            LDA )                          CALL DLASET( 'L', N-1, N-1, ZERO, ZERO,
        $                               A( 2, 1 ), LDA )
                        END IF
 *  *
 *                    Bidiagonalize R in A  *                    Bidiagonalize R in A
 *                    (Workspace: need 4*N, prefer 3*N+2*N*NB)  *                    (Workspace: need 4*N, prefer 3*N + 2*N*NB)
 *  *
                      CALL DGEBRD( N, N, A, LDA, S, WORK( IE ),                       CALL DGEBRD( N, N, A, LDA, S, WORK( IE ),
      $                            WORK( ITAUQ ), WORK( ITAUP ),       $                            WORK( ITAUQ ), WORK( ITAUP ),
      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )       $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
 *  *
 *                    Multiply Q in U by left vectors bidiagonalizing R  *                    Multiply Q in U by left vectors bidiagonalizing R
 *                    (Workspace: need 3*N+M, prefer 3*N+M*NB)  *                    (Workspace: need 3*N + M, prefer 3*N + M*NB)
 *  *
                      CALL DORMBR( 'Q', 'R', 'N', M, N, N, A, LDA,                       CALL DORMBR( 'Q', 'R', 'N', M, N, N, A, LDA,
      $                            WORK( ITAUQ ), U, LDU, WORK( IWORK ),       $                            WORK( ITAUQ ), U, LDU, WORK( IWORK ),
Line 1167 Line 1169
                         LDWRKU = LDA                          LDWRKU = LDA
                         IR = IU + LDWRKU*N                          IR = IU + LDWRKU*N
                         LDWRKR = LDA                          LDWRKR = LDA
                      ELSE IF( LWORK.GE.WRKBL+( LDA+N )*N ) THEN                       ELSE IF( LWORK.GE.WRKBL+( LDA + N )*N ) THEN
 *  *
 *                       WORK(IU) is LDA by N and WORK(IR) is N by N  *                       WORK(IU) is LDA by N and WORK(IR) is N by N
 *  *
Line 1186 Line 1188
                      IWORK = ITAU + N                       IWORK = ITAU + N
 *  *
 *                    Compute A=Q*R  *                    Compute A=Q*R
 *                    (Workspace: need 2*N*N+2*N, prefer 2*N*N+N+N*NB)  *                    (Workspace: need 2*N*N + 2*N, prefer 2*N*N + N + N*NB)
 *  *
                      CALL DGEQRF( M, N, A, LDA, WORK( ITAU ),                       CALL DGEQRF( M, N, A, LDA, WORK( ITAU ),
      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )       $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
Line 1199 Line 1201
      $                            WORK( IU+1 ), LDWRKU )       $                            WORK( IU+1 ), LDWRKU )
 *  *
 *                    Generate Q in A  *                    Generate Q in A
 *                    (Workspace: need 2*N*N+2*N, prefer 2*N*N+N+N*NB)  *                    (Workspace: need 2*N*N + 2*N, prefer 2*N*N + N + N*NB)
 *  *
                      CALL DORGQR( M, N, N, A, LDA, WORK( ITAU ),                       CALL DORGQR( M, N, N, A, LDA, WORK( ITAU ),
      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )       $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
Line 1210 Line 1212
 *  *
 *                    Bidiagonalize R in WORK(IU), copying result to  *                    Bidiagonalize R in WORK(IU), copying result to
 *                    WORK(IR)  *                    WORK(IR)
 *                    (Workspace: need 2*N*N+4*N,  *                    (Workspace: need 2*N*N + 4*N,
 *                                prefer 2*N*N+3*N+2*N*NB)  *                                prefer 2*N*N+3*N+2*N*NB)
 *  *
                      CALL DGEBRD( N, N, WORK( IU ), LDWRKU, S,                       CALL DGEBRD( N, N, WORK( IU ), LDWRKU, S,
Line 1221 Line 1223
      $                            WORK( IR ), LDWRKR )       $                            WORK( IR ), LDWRKR )
 *  *
 *                    Generate left bidiagonalizing vectors in WORK(IU)  *                    Generate left bidiagonalizing vectors in WORK(IU)
 *                    (Workspace: need 2*N*N+4*N, prefer 2*N*N+3*N+N*NB)  *                    (Workspace: need 2*N*N + 4*N, prefer 2*N*N + 3*N + N*NB)
 *  *
                      CALL DORGBR( 'Q', N, N, N, WORK( IU ), LDWRKU,                       CALL DORGBR( 'Q', N, N, N, WORK( IU ), LDWRKU,
      $                            WORK( ITAUQ ), WORK( IWORK ),       $                            WORK( ITAUQ ), WORK( IWORK ),
      $                            LWORK-IWORK+1, IERR )       $                            LWORK-IWORK+1, IERR )
 *  *
 *                    Generate right bidiagonalizing vectors in WORK(IR)  *                    Generate right bidiagonalizing vectors in WORK(IR)
 *                    (Workspace: need 2*N*N+4*N-1,  *                    (Workspace: need 2*N*N + 4*N-1,
 *                                prefer 2*N*N+3*N+(N-1)*NB)  *                                prefer 2*N*N+3*N+(N-1)*NB)
 *  *
                      CALL DORGBR( 'P', N, N, N, WORK( IR ), LDWRKR,                       CALL DORGBR( 'P', N, N, N, WORK( IR ), LDWRKR,
Line 1239 Line 1241
 *                    Perform bidiagonal QR iteration, computing left  *                    Perform bidiagonal QR iteration, computing left
 *                    singular vectors of R in WORK(IU) and computing  *                    singular vectors of R in WORK(IU) and computing
 *                    right singular vectors of R in WORK(IR)  *                    right singular vectors of R in WORK(IR)
 *                    (Workspace: need 2*N*N+BDSPAC)  *                    (Workspace: need 2*N*N + BDSPAC)
 *  *
                      CALL DBDSQR( 'U', N, N, N, 0, S, WORK( IE ),                       CALL DBDSQR( 'U', N, N, N, 0, S, WORK( IE ),
      $                            WORK( IR ), LDWRKR, WORK( IU ),       $                            WORK( IR ), LDWRKR, WORK( IU ),
Line 1266 Line 1268
                      IWORK = ITAU + N                       IWORK = ITAU + N
 *  *
 *                    Compute A=Q*R, copying result to U  *                    Compute A=Q*R, copying result to U
 *                    (Workspace: need 2*N, prefer N+N*NB)  *                    (Workspace: need 2*N, prefer N + N*NB)
 *  *
                      CALL DGEQRF( M, N, A, LDA, WORK( ITAU ),                       CALL DGEQRF( M, N, A, LDA, WORK( ITAU ),
      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )       $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
                      CALL DLACPY( 'L', M, N, A, LDA, U, LDU )                       CALL DLACPY( 'L', M, N, A, LDA, U, LDU )
 *  *
 *                    Generate Q in U  *                    Generate Q in U
 *                    (Workspace: need 2*N, prefer N+N*NB)  *                    (Workspace: need 2*N, prefer N + N*NB)
 *  *
                      CALL DORGQR( M, N, N, U, LDU, WORK( ITAU ),                       CALL DORGQR( M, N, N, U, LDU, WORK( ITAU ),
      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )       $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
Line 1284 Line 1286
 *  *
 *                    Zero out below R in A  *                    Zero out below R in A
 *  *
                      CALL DLASET( 'L', N-1, N-1, ZERO, ZERO, A( 2, 1 ),                       IF( N .GT. 1 ) THEN
      $                            LDA )                          CALL DLASET( 'L', N-1, N-1, ZERO, ZERO,
        $                               A( 2, 1 ), LDA )
                        END IF
 *  *
 *                    Bidiagonalize R in A  *                    Bidiagonalize R in A
 *                    (Workspace: need 4*N, prefer 3*N+2*N*NB)  *                    (Workspace: need 4*N, prefer 3*N + 2*N*NB)
 *  *
                      CALL DGEBRD( N, N, A, LDA, S, WORK( IE ),                       CALL DGEBRD( N, N, A, LDA, S, WORK( IE ),
      $                            WORK( ITAUQ ), WORK( ITAUP ),       $                            WORK( ITAUQ ), WORK( ITAUP ),
      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )       $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
 *  *
 *                    Multiply Q in U by left vectors bidiagonalizing R  *                    Multiply Q in U by left vectors bidiagonalizing R
 *                    (Workspace: need 3*N+M, prefer 3*N+M*NB)  *                    (Workspace: need 3*N + M, prefer 3*N + M*NB)
 *  *
                      CALL DORMBR( 'Q', 'R', 'N', M, N, N, A, LDA,                       CALL DORMBR( 'Q', 'R', 'N', M, N, N, A, LDA,
      $                            WORK( ITAUQ ), U, LDU, WORK( IWORK ),       $                            WORK( ITAUQ ), U, LDU, WORK( IWORK ),
      $                            LWORK-IWORK+1, IERR )       $                            LWORK-IWORK+1, IERR )
 *  *
 *                    Generate right vectors bidiagonalizing R in A  *                    Generate right vectors bidiagonalizing R in A
 *                    (Workspace: need 4*N-1, prefer 3*N+(N-1)*NB)  *                    (Workspace: need 4*N-1, prefer 3*N + (N-1)*NB)
 *  *
                      CALL DORGBR( 'P', N, N, N, A, LDA, WORK( ITAUP ),                       CALL DORGBR( 'P', N, N, N, A, LDA, WORK( ITAUP ),
      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )       $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
Line 1346 Line 1350
                      IWORK = ITAU + N                       IWORK = ITAU + N
 *  *
 *                    Compute A=Q*R  *                    Compute A=Q*R
 *                    (Workspace: need N*N+2*N, prefer N*N+N+N*NB)  *                    (Workspace: need N*N + 2*N, prefer N*N + N + N*NB)
 *  *
                      CALL DGEQRF( M, N, A, LDA, WORK( ITAU ),                       CALL DGEQRF( M, N, A, LDA, WORK( ITAU ),
      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )       $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
Line 1359 Line 1363
      $                            WORK( IU+1 ), LDWRKU )       $                            WORK( IU+1 ), LDWRKU )
 *  *
 *                    Generate Q in A  *                    Generate Q in A
 *                    (Workspace: need N*N+2*N, prefer N*N+N+N*NB)  *                    (Workspace: need N*N + 2*N, prefer N*N + N + N*NB)
 *  *
                      CALL DORGQR( M, N, N, A, LDA, WORK( ITAU ),                       CALL DORGQR( M, N, N, A, LDA, WORK( ITAU ),
      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )       $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
Line 1369 Line 1373
                      IWORK = ITAUP + N                       IWORK = ITAUP + N
 *  *
 *                    Bidiagonalize R in WORK(IU), copying result to VT  *                    Bidiagonalize R in WORK(IU), copying result to VT
 *                    (Workspace: need N*N+4*N, prefer N*N+3*N+2*N*NB)  *                    (Workspace: need N*N + 4*N, prefer N*N + 3*N + 2*N*NB)
 *  *
                      CALL DGEBRD( N, N, WORK( IU ), LDWRKU, S,                       CALL DGEBRD( N, N, WORK( IU ), LDWRKU, S,
      $                            WORK( IE ), WORK( ITAUQ ),       $                            WORK( IE ), WORK( ITAUQ ),
Line 1379 Line 1383
      $                            LDVT )       $                            LDVT )
 *  *
 *                    Generate left bidiagonalizing vectors in WORK(IU)  *                    Generate left bidiagonalizing vectors in WORK(IU)
 *                    (Workspace: need N*N+4*N, prefer N*N+3*N+N*NB)  *                    (Workspace: need N*N + 4*N, prefer N*N + 3*N + N*NB)
 *  *
                      CALL DORGBR( 'Q', N, N, N, WORK( IU ), LDWRKU,                       CALL DORGBR( 'Q', N, N, N, WORK( IU ), LDWRKU,
      $                            WORK( ITAUQ ), WORK( IWORK ),       $                            WORK( ITAUQ ), WORK( IWORK ),
      $                            LWORK-IWORK+1, IERR )       $                            LWORK-IWORK+1, IERR )
 *  *
 *                    Generate right bidiagonalizing vectors in VT  *                    Generate right bidiagonalizing vectors in VT
 *                    (Workspace: need N*N+4*N-1,  *                    (Workspace: need N*N + 4*N-1,
 *                                prefer N*N+3*N+(N-1)*NB)  *                                prefer N*N+3*N+(N-1)*NB)
 *  *
                      CALL DORGBR( 'P', N, N, N, VT, LDVT, WORK( ITAUP ),                       CALL DORGBR( 'P', N, N, N, VT, LDVT, WORK( ITAUP ),
Line 1396 Line 1400
 *                    Perform bidiagonal QR iteration, computing left  *                    Perform bidiagonal QR iteration, computing left
 *                    singular vectors of R in WORK(IU) and computing  *                    singular vectors of R in WORK(IU) and computing
 *                    right singular vectors of R in VT  *                    right singular vectors of R in VT
 *                    (Workspace: need N*N+BDSPAC)  *                    (Workspace: need N*N + BDSPAC)
 *  *
                      CALL DBDSQR( 'U', N, N, N, 0, S, WORK( IE ), VT,                       CALL DBDSQR( 'U', N, N, N, 0, S, WORK( IE ), VT,
      $                            LDVT, WORK( IU ), LDWRKU, DUM, 1,       $                            LDVT, WORK( IU ), LDWRKU, DUM, 1,
Line 1417 Line 1421
                      IWORK = ITAU + N                       IWORK = ITAU + N
 *  *
 *                    Compute A=Q*R, copying result to U  *                    Compute A=Q*R, copying result to U
 *                    (Workspace: need 2*N, prefer N+N*NB)  *                    (Workspace: need 2*N, prefer N + N*NB)
 *  *
                      CALL DGEQRF( M, N, A, LDA, WORK( ITAU ),                       CALL DGEQRF( M, N, A, LDA, WORK( ITAU ),
      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )       $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
                      CALL DLACPY( 'L', M, N, A, LDA, U, LDU )                       CALL DLACPY( 'L', M, N, A, LDA, U, LDU )
 *  *
 *                    Generate Q in U  *                    Generate Q in U
 *                    (Workspace: need 2*N, prefer N+N*NB)  *                    (Workspace: need 2*N, prefer N + N*NB)
 *  *
                      CALL DORGQR( M, N, N, U, LDU, WORK( ITAU ),                       CALL DORGQR( M, N, N, U, LDU, WORK( ITAU ),
      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )       $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
Line 1441 Line 1445
                      IWORK = ITAUP + N                       IWORK = ITAUP + N
 *  *
 *                    Bidiagonalize R in VT  *                    Bidiagonalize R in VT
 *                    (Workspace: need 4*N, prefer 3*N+2*N*NB)  *                    (Workspace: need 4*N, prefer 3*N + 2*N*NB)
 *  *
                      CALL DGEBRD( N, N, VT, LDVT, S, WORK( IE ),                       CALL DGEBRD( N, N, VT, LDVT, S, WORK( IE ),
      $                            WORK( ITAUQ ), WORK( ITAUP ),       $                            WORK( ITAUQ ), WORK( ITAUP ),
Line 1449 Line 1453
 *  *
 *                    Multiply Q in U by left bidiagonalizing vectors  *                    Multiply Q in U by left bidiagonalizing vectors
 *                    in VT  *                    in VT
 *                    (Workspace: need 3*N+M, prefer 3*N+M*NB)  *                    (Workspace: need 3*N + M, prefer 3*N + M*NB)
 *  *
                      CALL DORMBR( 'Q', 'R', 'N', M, N, N, VT, LDVT,                       CALL DORMBR( 'Q', 'R', 'N', M, N, N, VT, LDVT,
      $                            WORK( ITAUQ ), U, LDU, WORK( IWORK ),       $                            WORK( ITAUQ ), U, LDU, WORK( IWORK ),
      $                            LWORK-IWORK+1, IERR )       $                            LWORK-IWORK+1, IERR )
 *  *
 *                    Generate right bidiagonalizing vectors in VT  *                    Generate right bidiagonalizing vectors in VT
 *                    (Workspace: need 4*N-1, prefer 3*N+(N-1)*NB)  *                    (Workspace: need 4*N-1, prefer 3*N + (N-1)*NB)
 *  *
                      CALL DORGBR( 'P', N, N, N, VT, LDVT, WORK( ITAUP ),                       CALL DORGBR( 'P', N, N, N, VT, LDVT, WORK( ITAUP ),
      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )       $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
Line 1503 Line 1507
                      IWORK = ITAU + N                       IWORK = ITAU + N
 *  *
 *                    Compute A=Q*R, copying result to U  *                    Compute A=Q*R, copying result to U
 *                    (Workspace: need N*N+2*N, prefer N*N+N+N*NB)  *                    (Workspace: need N*N + 2*N, prefer N*N + N + N*NB)
 *  *
                      CALL DGEQRF( M, N, A, LDA, WORK( ITAU ),                       CALL DGEQRF( M, N, A, LDA, WORK( ITAU ),
      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )       $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
Line 1517 Line 1521
      $                            WORK( IR+1 ), LDWRKR )       $                            WORK( IR+1 ), LDWRKR )
 *  *
 *                    Generate Q in U  *                    Generate Q in U
 *                    (Workspace: need N*N+N+M, prefer N*N+N+M*NB)  *                    (Workspace: need N*N + N + M, prefer N*N + N + M*NB)
 *  *
                      CALL DORGQR( M, M, N, U, LDU, WORK( ITAU ),                       CALL DORGQR( M, M, N, U, LDU, WORK( ITAU ),
      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )       $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
Line 1527 Line 1531
                      IWORK = ITAUP + N                       IWORK = ITAUP + N
 *  *
 *                    Bidiagonalize R in WORK(IR)  *                    Bidiagonalize R in WORK(IR)
 *                    (Workspace: need N*N+4*N, prefer N*N+3*N+2*N*NB)  *                    (Workspace: need N*N + 4*N, prefer N*N + 3*N + 2*N*NB)
 *  *
                      CALL DGEBRD( N, N, WORK( IR ), LDWRKR, S,                       CALL DGEBRD( N, N, WORK( IR ), LDWRKR, S,
      $                            WORK( IE ), WORK( ITAUQ ),       $                            WORK( IE ), WORK( ITAUQ ),
Line 1535 Line 1539
      $                            LWORK-IWORK+1, IERR )       $                            LWORK-IWORK+1, IERR )
 *  *
 *                    Generate left bidiagonalizing vectors in WORK(IR)  *                    Generate left bidiagonalizing vectors in WORK(IR)
 *                    (Workspace: need N*N+4*N, prefer N*N+3*N+N*NB)  *                    (Workspace: need N*N + 4*N, prefer N*N + 3*N + N*NB)
 *  *
                      CALL DORGBR( 'Q', N, N, N, WORK( IR ), LDWRKR,                       CALL DORGBR( 'Q', N, N, N, WORK( IR ), LDWRKR,
      $                            WORK( ITAUQ ), WORK( IWORK ),       $                            WORK( ITAUQ ), WORK( IWORK ),
Line 1544 Line 1548
 *  *
 *                    Perform bidiagonal QR iteration, computing left  *                    Perform bidiagonal QR iteration, computing left
 *                    singular vectors of R in WORK(IR)  *                    singular vectors of R in WORK(IR)
 *                    (Workspace: need N*N+BDSPAC)  *                    (Workspace: need N*N + BDSPAC)
 *  *
                      CALL DBDSQR( 'U', N, 0, N, 0, S, WORK( IE ), DUM,                       CALL DBDSQR( 'U', N, 0, N, 0, S, WORK( IE ), DUM,
      $                            1, WORK( IR ), LDWRKR, DUM, 1,       $                            1, WORK( IR ), LDWRKR, DUM, 1,
Line 1569 Line 1573
                      IWORK = ITAU + N                       IWORK = ITAU + N
 *  *
 *                    Compute A=Q*R, copying result to U  *                    Compute A=Q*R, copying result to U
 *                    (Workspace: need 2*N, prefer N+N*NB)  *                    (Workspace: need 2*N, prefer N + N*NB)
 *  *
                      CALL DGEQRF( M, N, A, LDA, WORK( ITAU ),                       CALL DGEQRF( M, N, A, LDA, WORK( ITAU ),
      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )       $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
                      CALL DLACPY( 'L', M, N, A, LDA, U, LDU )                       CALL DLACPY( 'L', M, N, A, LDA, U, LDU )
 *  *
 *                    Generate Q in U  *                    Generate Q in U
 *                    (Workspace: need N+M, prefer N+M*NB)  *                    (Workspace: need N + M, prefer N + M*NB)
 *  *
                      CALL DORGQR( M, M, N, U, LDU, WORK( ITAU ),                       CALL DORGQR( M, M, N, U, LDU, WORK( ITAU ),
      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )       $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
Line 1587 Line 1591
 *  *
 *                    Zero out below R in A  *                    Zero out below R in A
 *  *
                      CALL DLASET( 'L', N-1, N-1, ZERO, ZERO, A( 2, 1 ),                       IF( N .GT. 1 ) THEN
      $                            LDA )                          CALL DLASET( 'L', N-1, N-1, ZERO, ZERO,
        $                                A( 2, 1 ), LDA )
                        END IF
 *  *
 *                    Bidiagonalize R in A  *                    Bidiagonalize R in A
 *                    (Workspace: need 4*N, prefer 3*N+2*N*NB)  *                    (Workspace: need 4*N, prefer 3*N + 2*N*NB)
 *  *
                      CALL DGEBRD( N, N, A, LDA, S, WORK( IE ),                       CALL DGEBRD( N, N, A, LDA, S, WORK( IE ),
      $                            WORK( ITAUQ ), WORK( ITAUP ),       $                            WORK( ITAUQ ), WORK( ITAUP ),
Line 1599 Line 1605
 *  *
 *                    Multiply Q in U by left bidiagonalizing vectors  *                    Multiply Q in U by left bidiagonalizing vectors
 *                    in A  *                    in A
 *                    (Workspace: need 3*N+M, prefer 3*N+M*NB)  *                    (Workspace: need 3*N + M, prefer 3*N + M*NB)
 *  *
                      CALL DORMBR( 'Q', 'R', 'N', M, N, N, A, LDA,                       CALL DORMBR( 'Q', 'R', 'N', M, N, N, A, LDA,
      $                            WORK( ITAUQ ), U, LDU, WORK( IWORK ),       $                            WORK( ITAUQ ), U, LDU, WORK( IWORK ),
Line 1634 Line 1640
                         LDWRKU = LDA                          LDWRKU = LDA
                         IR = IU + LDWRKU*N                          IR = IU + LDWRKU*N
                         LDWRKR = LDA                          LDWRKR = LDA
                      ELSE IF( LWORK.GE.WRKBL+( LDA+N )*N ) THEN                       ELSE IF( LWORK.GE.WRKBL+( LDA + N )*N ) THEN
 *  *
 *                       WORK(IU) is LDA by N and WORK(IR) is N by N  *                       WORK(IU) is LDA by N and WORK(IR) is N by N
 *  *
Line 1653 Line 1659
                      IWORK = ITAU + N                       IWORK = ITAU + N
 *  *
 *                    Compute A=Q*R, copying result to U  *                    Compute A=Q*R, copying result to U
 *                    (Workspace: need 2*N*N+2*N, prefer 2*N*N+N+N*NB)  *                    (Workspace: need 2*N*N + 2*N, prefer 2*N*N + N + N*NB)
 *  *
                      CALL DGEQRF( M, N, A, LDA, WORK( ITAU ),                       CALL DGEQRF( M, N, A, LDA, WORK( ITAU ),
      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )       $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
                      CALL DLACPY( 'L', M, N, A, LDA, U, LDU )                       CALL DLACPY( 'L', M, N, A, LDA, U, LDU )
 *  *
 *                    Generate Q in U  *                    Generate Q in U
 *                    (Workspace: need 2*N*N+N+M, prefer 2*N*N+N+M*NB)  *                    (Workspace: need 2*N*N + N + M, prefer 2*N*N + N + M*NB)
 *  *
                      CALL DORGQR( M, M, N, U, LDU, WORK( ITAU ),                       CALL DORGQR( M, M, N, U, LDU, WORK( ITAU ),
      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )       $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
Line 1678 Line 1684
 *  *
 *                    Bidiagonalize R in WORK(IU), copying result to  *                    Bidiagonalize R in WORK(IU), copying result to
 *                    WORK(IR)  *                    WORK(IR)
 *                    (Workspace: need 2*N*N+4*N,  *                    (Workspace: need 2*N*N + 4*N,
 *                                prefer 2*N*N+3*N+2*N*NB)  *                                prefer 2*N*N+3*N+2*N*NB)
 *  *
                      CALL DGEBRD( N, N, WORK( IU ), LDWRKU, S,                       CALL DGEBRD( N, N, WORK( IU ), LDWRKU, S,
Line 1689 Line 1695
      $                            WORK( IR ), LDWRKR )       $                            WORK( IR ), LDWRKR )
 *  *
 *                    Generate left bidiagonalizing vectors in WORK(IU)  *                    Generate left bidiagonalizing vectors in WORK(IU)
 *                    (Workspace: need 2*N*N+4*N, prefer 2*N*N+3*N+N*NB)  *                    (Workspace: need 2*N*N + 4*N, prefer 2*N*N + 3*N + N*NB)
 *  *
                      CALL DORGBR( 'Q', N, N, N, WORK( IU ), LDWRKU,                       CALL DORGBR( 'Q', N, N, N, WORK( IU ), LDWRKU,
      $                            WORK( ITAUQ ), WORK( IWORK ),       $                            WORK( ITAUQ ), WORK( IWORK ),
      $                            LWORK-IWORK+1, IERR )       $                            LWORK-IWORK+1, IERR )
 *  *
 *                    Generate right bidiagonalizing vectors in WORK(IR)  *                    Generate right bidiagonalizing vectors in WORK(IR)
 *                    (Workspace: need 2*N*N+4*N-1,  *                    (Workspace: need 2*N*N + 4*N-1,
 *                                prefer 2*N*N+3*N+(N-1)*NB)  *                                prefer 2*N*N+3*N+(N-1)*NB)
 *  *
                      CALL DORGBR( 'P', N, N, N, WORK( IR ), LDWRKR,                       CALL DORGBR( 'P', N, N, N, WORK( IR ), LDWRKR,
Line 1707 Line 1713
 *                    Perform bidiagonal QR iteration, computing left  *                    Perform bidiagonal QR iteration, computing left
 *                    singular vectors of R in WORK(IU) and computing  *                    singular vectors of R in WORK(IU) and computing
 *                    right singular vectors of R in WORK(IR)  *                    right singular vectors of R in WORK(IR)
 *                    (Workspace: need 2*N*N+BDSPAC)  *                    (Workspace: need 2*N*N + BDSPAC)
 *  *
                      CALL DBDSQR( 'U', N, N, N, 0, S, WORK( IE ),                       CALL DBDSQR( 'U', N, N, N, 0, S, WORK( IE ),
      $                            WORK( IR ), LDWRKR, WORK( IU ),       $                            WORK( IR ), LDWRKR, WORK( IU ),
Line 1737 Line 1743
                      IWORK = ITAU + N                       IWORK = ITAU + N
 *  *
 *                    Compute A=Q*R, copying result to U  *                    Compute A=Q*R, copying result to U
 *                    (Workspace: need 2*N, prefer N+N*NB)  *                    (Workspace: need 2*N, prefer N + N*NB)
 *  *
                      CALL DGEQRF( M, N, A, LDA, WORK( ITAU ),                       CALL DGEQRF( M, N, A, LDA, WORK( ITAU ),
      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )       $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
                      CALL DLACPY( 'L', M, N, A, LDA, U, LDU )                       CALL DLACPY( 'L', M, N, A, LDA, U, LDU )
 *  *
 *                    Generate Q in U  *                    Generate Q in U
 *                    (Workspace: need N+M, prefer N+M*NB)  *                    (Workspace: need N + M, prefer N + M*NB)
 *  *
                      CALL DORGQR( M, M, N, U, LDU, WORK( ITAU ),                       CALL DORGQR( M, M, N, U, LDU, WORK( ITAU ),
      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )       $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
Line 1755 Line 1761
 *  *
 *                    Zero out below R in A  *                    Zero out below R in A
 *  *
                      CALL DLASET( 'L', N-1, N-1, ZERO, ZERO, A( 2, 1 ),                       IF( N .GT. 1 ) THEN
      $                            LDA )                          CALL DLASET( 'L', N-1, N-1, ZERO, ZERO,
        $                                A( 2, 1 ), LDA )
                        END IF
 *  *
 *                    Bidiagonalize R in A  *                    Bidiagonalize R in A
 *                    (Workspace: need 4*N, prefer 3*N+2*N*NB)  *                    (Workspace: need 4*N, prefer 3*N + 2*N*NB)
 *  *
                      CALL DGEBRD( N, N, A, LDA, S, WORK( IE ),                       CALL DGEBRD( N, N, A, LDA, S, WORK( IE ),
      $                            WORK( ITAUQ ), WORK( ITAUP ),       $                            WORK( ITAUQ ), WORK( ITAUP ),
Line 1767 Line 1775
 *  *
 *                    Multiply Q in U by left bidiagonalizing vectors  *                    Multiply Q in U by left bidiagonalizing vectors
 *                    in A  *                    in A
 *                    (Workspace: need 3*N+M, prefer 3*N+M*NB)  *                    (Workspace: need 3*N + M, prefer 3*N + M*NB)
 *  *
                      CALL DORMBR( 'Q', 'R', 'N', M, N, N, A, LDA,                       CALL DORMBR( 'Q', 'R', 'N', M, N, N, A, LDA,
      $                            WORK( ITAUQ ), U, LDU, WORK( IWORK ),       $                            WORK( ITAUQ ), U, LDU, WORK( IWORK ),
      $                            LWORK-IWORK+1, IERR )       $                            LWORK-IWORK+1, IERR )
 *  *
 *                    Generate right bidiagonalizing vectors in A  *                    Generate right bidiagonalizing vectors in A
 *                    (Workspace: need 4*N-1, prefer 3*N+(N-1)*NB)  *                    (Workspace: need 4*N-1, prefer 3*N + (N-1)*NB)
 *  *
                      CALL DORGBR( 'P', N, N, N, A, LDA, WORK( ITAUP ),                       CALL DORGBR( 'P', N, N, N, A, LDA, WORK( ITAUP ),
      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )       $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
Line 1818 Line 1826
                      IWORK = ITAU + N                       IWORK = ITAU + N
 *  *
 *                    Compute A=Q*R, copying result to U  *                    Compute A=Q*R, copying result to U
 *                    (Workspace: need N*N+2*N, prefer N*N+N+N*NB)  *                    (Workspace: need N*N + 2*N, prefer N*N + N + N*NB)
 *  *
                      CALL DGEQRF( M, N, A, LDA, WORK( ITAU ),                       CALL DGEQRF( M, N, A, LDA, WORK( ITAU ),
      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )       $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
                      CALL DLACPY( 'L', M, N, A, LDA, U, LDU )                       CALL DLACPY( 'L', M, N, A, LDA, U, LDU )
 *  *
 *                    Generate Q in U  *                    Generate Q in U
 *                    (Workspace: need N*N+N+M, prefer N*N+N+M*NB)  *                    (Workspace: need N*N + N + M, prefer N*N + N + M*NB)
 *  *
                      CALL DORGQR( M, M, N, U, LDU, WORK( ITAU ),                       CALL DORGQR( M, M, N, U, LDU, WORK( ITAU ),
      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )       $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
Line 1842 Line 1850
                      IWORK = ITAUP + N                       IWORK = ITAUP + N
 *  *
 *                    Bidiagonalize R in WORK(IU), copying result to VT  *                    Bidiagonalize R in WORK(IU), copying result to VT
 *                    (Workspace: need N*N+4*N, prefer N*N+3*N+2*N*NB)  *                    (Workspace: need N*N + 4*N, prefer N*N + 3*N + 2*N*NB)
 *  *
                      CALL DGEBRD( N, N, WORK( IU ), LDWRKU, S,                       CALL DGEBRD( N, N, WORK( IU ), LDWRKU, S,
      $                            WORK( IE ), WORK( ITAUQ ),       $                            WORK( IE ), WORK( ITAUQ ),
Line 1852 Line 1860
      $                            LDVT )       $                            LDVT )
 *  *
 *                    Generate left bidiagonalizing vectors in WORK(IU)  *                    Generate left bidiagonalizing vectors in WORK(IU)
 *                    (Workspace: need N*N+4*N, prefer N*N+3*N+N*NB)  *                    (Workspace: need N*N + 4*N, prefer N*N + 3*N + N*NB)
 *  *
                      CALL DORGBR( 'Q', N, N, N, WORK( IU ), LDWRKU,                       CALL DORGBR( 'Q', N, N, N, WORK( IU ), LDWRKU,
      $                            WORK( ITAUQ ), WORK( IWORK ),       $                            WORK( ITAUQ ), WORK( IWORK ),
      $                            LWORK-IWORK+1, IERR )       $                            LWORK-IWORK+1, IERR )
 *  *
 *                    Generate right bidiagonalizing vectors in VT  *                    Generate right bidiagonalizing vectors in VT
 *                    (Workspace: need N*N+4*N-1,  *                    (Workspace: need N*N + 4*N-1,
 *                                prefer N*N+3*N+(N-1)*NB)  *                                prefer N*N+3*N+(N-1)*NB)
 *  *
                      CALL DORGBR( 'P', N, N, N, VT, LDVT, WORK( ITAUP ),                       CALL DORGBR( 'P', N, N, N, VT, LDVT, WORK( ITAUP ),
Line 1869 Line 1877
 *                    Perform bidiagonal QR iteration, computing left  *                    Perform bidiagonal QR iteration, computing left
 *                    singular vectors of R in WORK(IU) and computing  *                    singular vectors of R in WORK(IU) and computing
 *                    right singular vectors of R in VT  *                    right singular vectors of R in VT
 *                    (Workspace: need N*N+BDSPAC)  *                    (Workspace: need N*N + BDSPAC)
 *  *
                      CALL DBDSQR( 'U', N, N, N, 0, S, WORK( IE ), VT,                       CALL DBDSQR( 'U', N, N, N, 0, S, WORK( IE ), VT,
      $                            LDVT, WORK( IU ), LDWRKU, DUM, 1,       $                            LDVT, WORK( IU ), LDWRKU, DUM, 1,
Line 1894 Line 1902
                      IWORK = ITAU + N                       IWORK = ITAU + N
 *  *
 *                    Compute A=Q*R, copying result to U  *                    Compute A=Q*R, copying result to U
 *                    (Workspace: need 2*N, prefer N+N*NB)  *                    (Workspace: need 2*N, prefer N + N*NB)
 *  *
                      CALL DGEQRF( M, N, A, LDA, WORK( ITAU ),                       CALL DGEQRF( M, N, A, LDA, WORK( ITAU ),
      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )       $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
                      CALL DLACPY( 'L', M, N, A, LDA, U, LDU )                       CALL DLACPY( 'L', M, N, A, LDA, U, LDU )
 *  *
 *                    Generate Q in U  *                    Generate Q in U
 *                    (Workspace: need N+M, prefer N+M*NB)  *                    (Workspace: need N + M, prefer N + M*NB)
 *  *
                      CALL DORGQR( M, M, N, U, LDU, WORK( ITAU ),                       CALL DORGQR( M, M, N, U, LDU, WORK( ITAU ),
      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )       $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
Line 1918 Line 1926
                      IWORK = ITAUP + N                       IWORK = ITAUP + N
 *  *
 *                    Bidiagonalize R in VT  *                    Bidiagonalize R in VT
 *                    (Workspace: need 4*N, prefer 3*N+2*N*NB)  *                    (Workspace: need 4*N, prefer 3*N + 2*N*NB)
 *  *
                      CALL DGEBRD( N, N, VT, LDVT, S, WORK( IE ),                       CALL DGEBRD( N, N, VT, LDVT, S, WORK( IE ),
      $                            WORK( ITAUQ ), WORK( ITAUP ),       $                            WORK( ITAUQ ), WORK( ITAUP ),
Line 1926 Line 1934
 *  *
 *                    Multiply Q in U by left bidiagonalizing vectors  *                    Multiply Q in U by left bidiagonalizing vectors
 *                    in VT  *                    in VT
 *                    (Workspace: need 3*N+M, prefer 3*N+M*NB)  *                    (Workspace: need 3*N + M, prefer 3*N + M*NB)
 *  *
                      CALL DORMBR( 'Q', 'R', 'N', M, N, N, VT, LDVT,                       CALL DORMBR( 'Q', 'R', 'N', M, N, N, VT, LDVT,
      $                            WORK( ITAUQ ), U, LDU, WORK( IWORK ),       $                            WORK( ITAUQ ), U, LDU, WORK( IWORK ),
      $                            LWORK-IWORK+1, IERR )       $                            LWORK-IWORK+1, IERR )
 *  *
 *                    Generate right bidiagonalizing vectors in VT  *                    Generate right bidiagonalizing vectors in VT
 *                    (Workspace: need 4*N-1, prefer 3*N+(N-1)*NB)  *                    (Workspace: need 4*N-1, prefer 3*N + (N-1)*NB)
 *  *
                      CALL DORGBR( 'P', N, N, N, VT, LDVT, WORK( ITAUP ),                       CALL DORGBR( 'P', N, N, N, VT, LDVT, WORK( ITAUP ),
      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )       $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
Line 1967 Line 1975
             IWORK = ITAUP + N              IWORK = ITAUP + N
 *  *
 *           Bidiagonalize A  *           Bidiagonalize A
 *           (Workspace: need 3*N+M, prefer 3*N+(M+N)*NB)  *           (Workspace: need 3*N + M, prefer 3*N + (M + N)*NB)
 *  *
             CALL DGEBRD( M, N, A, LDA, S, WORK( IE ), WORK( ITAUQ ),              CALL DGEBRD( M, N, A, LDA, S, WORK( IE ), WORK( ITAUQ ),
      $                   WORK( ITAUP ), WORK( IWORK ), LWORK-IWORK+1,       $                   WORK( ITAUP ), WORK( IWORK ), LWORK-IWORK+1,
Line 1976 Line 1984
 *  *
 *              If left singular vectors desired in U, copy result to U  *              If left singular vectors desired in U, copy result to U
 *              and generate left bidiagonalizing vectors in U  *              and generate left bidiagonalizing vectors in U
 *              (Workspace: need 3*N+NCU, prefer 3*N+NCU*NB)  *              (Workspace: need 3*N + NCU, prefer 3*N + NCU*NB)
 *  *
                CALL DLACPY( 'L', M, N, A, LDA, U, LDU )                 CALL DLACPY( 'L', M, N, A, LDA, U, LDU )
                IF( WNTUS )                 IF( WNTUS )
Line 1990 Line 1998
 *  *
 *              If right singular vectors desired in VT, copy result to  *              If right singular vectors desired in VT, copy result to
 *              VT and generate right bidiagonalizing vectors in VT  *              VT and generate right bidiagonalizing vectors in VT
 *              (Workspace: need 4*N-1, prefer 3*N+(N-1)*NB)  *              (Workspace: need 4*N-1, prefer 3*N + (N-1)*NB)
 *  *
                CALL DLACPY( 'U', N, N, A, LDA, VT, LDVT )                 CALL DLACPY( 'U', N, N, A, LDA, VT, LDVT )
                CALL DORGBR( 'P', N, N, N, VT, LDVT, WORK( ITAUP ),                 CALL DORGBR( 'P', N, N, N, VT, LDVT, WORK( ITAUP ),
Line 2000 Line 2008
 *  *
 *              If left singular vectors desired in A, generate left  *              If left singular vectors desired in A, generate left
 *              bidiagonalizing vectors in A  *              bidiagonalizing vectors in A
 *              (Workspace: need 4*N, prefer 3*N+N*NB)  *              (Workspace: need 4*N, prefer 3*N + N*NB)
 *  *
                CALL DORGBR( 'Q', M, N, N, A, LDA, WORK( ITAUQ ),                 CALL DORGBR( 'Q', M, N, N, A, LDA, WORK( ITAUQ ),
      $                      WORK( IWORK ), LWORK-IWORK+1, IERR )       $                      WORK( IWORK ), LWORK-IWORK+1, IERR )
Line 2009 Line 2017
 *  *
 *              If right singular vectors desired in A, generate right  *              If right singular vectors desired in A, generate right
 *              bidiagonalizing vectors in A  *              bidiagonalizing vectors in A
 *              (Workspace: need 4*N-1, prefer 3*N+(N-1)*NB)  *              (Workspace: need 4*N-1, prefer 3*N + (N-1)*NB)
 *  *
                CALL DORGBR( 'P', N, N, N, A, LDA, WORK( ITAUP ),                 CALL DORGBR( 'P', N, N, N, A, LDA, WORK( ITAUP ),
      $                      WORK( IWORK ), LWORK-IWORK+1, IERR )       $                      WORK( IWORK ), LWORK-IWORK+1, IERR )
Line 2071 Line 2079
                IWORK = ITAU + M                 IWORK = ITAU + M
 *  *
 *              Compute A=L*Q  *              Compute A=L*Q
 *              (Workspace: need 2*M, prefer M+M*NB)  *              (Workspace: need 2*M, prefer M + M*NB)
 *  *
                CALL DGELQF( M, N, A, LDA, WORK( ITAU ), WORK( IWORK ),                 CALL DGELQF( M, N, A, LDA, WORK( ITAU ), WORK( IWORK ),
      $                      LWORK-IWORK+1, IERR )       $                      LWORK-IWORK+1, IERR )
Line 2085 Line 2093
                IWORK = ITAUP + M                 IWORK = ITAUP + M
 *  *
 *              Bidiagonalize L in A  *              Bidiagonalize L in A
 *              (Workspace: need 4*M, prefer 3*M+2*M*NB)  *              (Workspace: need 4*M, prefer 3*M + 2*M*NB)
 *  *
                CALL DGEBRD( M, M, A, LDA, S, WORK( IE ), WORK( ITAUQ ),                 CALL DGEBRD( M, M, A, LDA, S, WORK( IE ), WORK( ITAUQ ),
      $                      WORK( ITAUP ), WORK( IWORK ), LWORK-IWORK+1,       $                      WORK( ITAUP ), WORK( IWORK ), LWORK-IWORK+1,
Line 2093 Line 2101
                IF( WNTUO .OR. WNTUAS ) THEN                 IF( WNTUO .OR. WNTUAS ) THEN
 *  *
 *                 If left singular vectors desired, generate Q  *                 If left singular vectors desired, generate Q
 *                 (Workspace: need 4*M, prefer 3*M+M*NB)  *                 (Workspace: need 4*M, prefer 3*M + M*NB)
 *  *
                   CALL DORGBR( 'Q', M, M, M, A, LDA, WORK( ITAUQ ),                    CALL DORGBR( 'Q', M, M, M, A, LDA, WORK( ITAUQ ),
      $                         WORK( IWORK ), LWORK-IWORK+1, IERR )       $                         WORK( IWORK ), LWORK-IWORK+1, IERR )
Line 2126 Line 2134
 *                 Sufficient workspace for a fast algorithm  *                 Sufficient workspace for a fast algorithm
 *  *
                   IR = 1                    IR = 1
                   IF( LWORK.GE.MAX( WRKBL, LDA*N+M )+LDA*M ) THEN                    IF( LWORK.GE.MAX( WRKBL, LDA*N + M ) + LDA*M ) THEN
 *  *
 *                    WORK(IU) is LDA by N and WORK(IR) is LDA by M  *                    WORK(IU) is LDA by N and WORK(IR) is LDA by M
 *  *
                      LDWRKU = LDA                       LDWRKU = LDA
                      CHUNK = N                       CHUNK = N
                      LDWRKR = LDA                       LDWRKR = LDA
                   ELSE IF( LWORK.GE.MAX( WRKBL, LDA*N+M )+M*M ) THEN                    ELSE IF( LWORK.GE.MAX( WRKBL, LDA*N + M ) + M*M ) THEN
 *  *
 *                    WORK(IU) is LDA by N and WORK(IR) is M by M  *                    WORK(IU) is LDA by N and WORK(IR) is M by M
 *  *
Line 2152 Line 2160
                   IWORK = ITAU + M                    IWORK = ITAU + M
 *  *
 *                 Compute A=L*Q  *                 Compute A=L*Q
 *                 (Workspace: need M*M+2*M, prefer M*M+M+M*NB)  *                 (Workspace: need M*M + 2*M, prefer M*M + M + M*NB)
 *  *
                   CALL DGELQF( M, N, A, LDA, WORK( ITAU ),                    CALL DGELQF( M, N, A, LDA, WORK( ITAU ),
      $                         WORK( IWORK ), LWORK-IWORK+1, IERR )       $                         WORK( IWORK ), LWORK-IWORK+1, IERR )
Line 2164 Line 2172
      $                         WORK( IR+LDWRKR ), LDWRKR )       $                         WORK( IR+LDWRKR ), LDWRKR )
 *  *
 *                 Generate Q in A  *                 Generate Q in A
 *                 (Workspace: need M*M+2*M, prefer M*M+M+M*NB)  *                 (Workspace: need M*M + 2*M, prefer M*M + M + M*NB)
 *  *
                   CALL DORGLQ( M, N, M, A, LDA, WORK( ITAU ),                    CALL DORGLQ( M, N, M, A, LDA, WORK( ITAU ),
      $                         WORK( IWORK ), LWORK-IWORK+1, IERR )       $                         WORK( IWORK ), LWORK-IWORK+1, IERR )
Line 2174 Line 2182
                   IWORK = ITAUP + M                    IWORK = ITAUP + M
 *  *
 *                 Bidiagonalize L in WORK(IR)  *                 Bidiagonalize L in WORK(IR)
 *                 (Workspace: need M*M+4*M, prefer M*M+3*M+2*M*NB)  *                 (Workspace: need M*M + 4*M, prefer M*M + 3*M + 2*M*NB)
 *  *
                   CALL DGEBRD( M, M, WORK( IR ), LDWRKR, S, WORK( IE ),                    CALL DGEBRD( M, M, WORK( IR ), LDWRKR, S, WORK( IE ),
      $                         WORK( ITAUQ ), WORK( ITAUP ),       $                         WORK( ITAUQ ), WORK( ITAUP ),
      $                         WORK( IWORK ), LWORK-IWORK+1, IERR )       $                         WORK( IWORK ), LWORK-IWORK+1, IERR )
 *  *
 *                 Generate right vectors bidiagonalizing L  *                 Generate right vectors bidiagonalizing L
 *                 (Workspace: need M*M+4*M-1, prefer M*M+3*M+(M-1)*NB)  *                 (Workspace: need M*M + 4*M-1, prefer M*M + 3*M + (M-1)*NB)
 *  *
                   CALL DORGBR( 'P', M, M, M, WORK( IR ), LDWRKR,                    CALL DORGBR( 'P', M, M, M, WORK( IR ), LDWRKR,
      $                         WORK( ITAUP ), WORK( IWORK ),       $                         WORK( ITAUP ), WORK( IWORK ),
Line 2190 Line 2198
 *  *
 *                 Perform bidiagonal QR iteration, computing right  *                 Perform bidiagonal QR iteration, computing right
 *                 singular vectors of L in WORK(IR)  *                 singular vectors of L in WORK(IR)
 *                 (Workspace: need M*M+BDSPAC)  *                 (Workspace: need M*M + BDSPAC)
 *  *
                   CALL DBDSQR( 'U', M, M, 0, 0, S, WORK( IE ),                    CALL DBDSQR( 'U', M, M, 0, 0, S, WORK( IE ),
      $                         WORK( IR ), LDWRKR, DUM, 1, DUM, 1,       $                         WORK( IR ), LDWRKR, DUM, 1, DUM, 1,
Line 2199 Line 2207
 *  *
 *                 Multiply right singular vectors of L in WORK(IR) by Q  *                 Multiply right singular vectors of L in WORK(IR) by Q
 *                 in A, storing result in WORK(IU) and copying to A  *                 in A, storing result in WORK(IU) and copying to A
 *                 (Workspace: need M*M+2*M, prefer M*M+M*N+M)  *                 (Workspace: need M*M + 2*M, prefer M*M + M*N + M)
 *  *
                   DO 30 I = 1, N, CHUNK                    DO 30 I = 1, N, CHUNK
                      BLK = MIN( N-I+1, CHUNK )                       BLK = MIN( N-I+1, CHUNK )
Line 2220 Line 2228
                   IWORK = ITAUP + M                    IWORK = ITAUP + M
 *  *
 *                 Bidiagonalize A  *                 Bidiagonalize A
 *                 (Workspace: need 3*M+N, prefer 3*M+(M+N)*NB)  *                 (Workspace: need 3*M + N, prefer 3*M + (M + N)*NB)
 *  *
                   CALL DGEBRD( M, N, A, LDA, S, WORK( IE ),                    CALL DGEBRD( M, N, A, LDA, S, WORK( IE ),
      $                         WORK( ITAUQ ), WORK( ITAUP ),       $                         WORK( ITAUQ ), WORK( ITAUP ),
      $                         WORK( IWORK ), LWORK-IWORK+1, IERR )       $                         WORK( IWORK ), LWORK-IWORK+1, IERR )
 *  *
 *                 Generate right vectors bidiagonalizing A  *                 Generate right vectors bidiagonalizing A
 *                 (Workspace: need 4*M, prefer 3*M+M*NB)  *                 (Workspace: need 4*M, prefer 3*M + M*NB)
 *  *
                   CALL DORGBR( 'P', M, N, M, A, LDA, WORK( ITAUP ),                    CALL DORGBR( 'P', M, N, M, A, LDA, WORK( ITAUP ),
      $                         WORK( IWORK ), LWORK-IWORK+1, IERR )       $                         WORK( IWORK ), LWORK-IWORK+1, IERR )
Line 2253 Line 2261
 *                 Sufficient workspace for a fast algorithm  *                 Sufficient workspace for a fast algorithm
 *  *
                   IR = 1                    IR = 1
                   IF( LWORK.GE.MAX( WRKBL, LDA*N+M )+LDA*M ) THEN                    IF( LWORK.GE.MAX( WRKBL, LDA*N + M ) + LDA*M ) THEN
 *  *
 *                    WORK(IU) is LDA by N and WORK(IR) is LDA by M  *                    WORK(IU) is LDA by N and WORK(IR) is LDA by M
 *  *
                      LDWRKU = LDA                       LDWRKU = LDA
                      CHUNK = N                       CHUNK = N
                      LDWRKR = LDA                       LDWRKR = LDA
                   ELSE IF( LWORK.GE.MAX( WRKBL, LDA*N+M )+M*M ) THEN                    ELSE IF( LWORK.GE.MAX( WRKBL, LDA*N + M ) + M*M ) THEN
 *  *
 *                    WORK(IU) is LDA by N and WORK(IR) is M by M  *                    WORK(IU) is LDA by N and WORK(IR) is M by M
 *  *
Line 2279 Line 2287
                   IWORK = ITAU + M                    IWORK = ITAU + M
 *  *
 *                 Compute A=L*Q  *                 Compute A=L*Q
 *                 (Workspace: need M*M+2*M, prefer M*M+M+M*NB)  *                 (Workspace: need M*M + 2*M, prefer M*M + M + M*NB)
 *  *
                   CALL DGELQF( M, N, A, LDA, WORK( ITAU ),                    CALL DGELQF( M, N, A, LDA, WORK( ITAU ),
      $                         WORK( IWORK ), LWORK-IWORK+1, IERR )       $                         WORK( IWORK ), LWORK-IWORK+1, IERR )
Line 2291 Line 2299
      $                         LDU )       $                         LDU )
 *  *
 *                 Generate Q in A  *                 Generate Q in A
 *                 (Workspace: need M*M+2*M, prefer M*M+M+M*NB)  *                 (Workspace: need M*M + 2*M, prefer M*M + M + M*NB)
 *  *
                   CALL DORGLQ( M, N, M, A, LDA, WORK( ITAU ),                    CALL DORGLQ( M, N, M, A, LDA, WORK( ITAU ),
      $                         WORK( IWORK ), LWORK-IWORK+1, IERR )       $                         WORK( IWORK ), LWORK-IWORK+1, IERR )
Line 2301 Line 2309
                   IWORK = ITAUP + M                    IWORK = ITAUP + M
 *  *
 *                 Bidiagonalize L in U, copying result to WORK(IR)  *                 Bidiagonalize L in U, copying result to WORK(IR)
 *                 (Workspace: need M*M+4*M, prefer M*M+3*M+2*M*NB)  *                 (Workspace: need M*M + 4*M, prefer M*M + 3*M + 2*M*NB)
 *  *
                   CALL DGEBRD( M, M, U, LDU, S, WORK( IE ),                    CALL DGEBRD( M, M, U, LDU, S, WORK( IE ),
      $                         WORK( ITAUQ ), WORK( ITAUP ),       $                         WORK( ITAUQ ), WORK( ITAUP ),
Line 2309 Line 2317
                   CALL DLACPY( 'U', M, M, U, LDU, WORK( IR ), LDWRKR )                    CALL DLACPY( 'U', M, M, U, LDU, WORK( IR ), LDWRKR )
 *  *
 *                 Generate right vectors bidiagonalizing L in WORK(IR)  *                 Generate right vectors bidiagonalizing L in WORK(IR)
 *                 (Workspace: need M*M+4*M-1, prefer M*M+3*M+(M-1)*NB)  *                 (Workspace: need M*M + 4*M-1, prefer M*M + 3*M + (M-1)*NB)
 *  *
                   CALL DORGBR( 'P', M, M, M, WORK( IR ), LDWRKR,                    CALL DORGBR( 'P', M, M, M, WORK( IR ), LDWRKR,
      $                         WORK( ITAUP ), WORK( IWORK ),       $                         WORK( ITAUP ), WORK( IWORK ),
      $                         LWORK-IWORK+1, IERR )       $                         LWORK-IWORK+1, IERR )
 *  *
 *                 Generate left vectors bidiagonalizing L in U  *                 Generate left vectors bidiagonalizing L in U
 *                 (Workspace: need M*M+4*M, prefer M*M+3*M+M*NB)  *                 (Workspace: need M*M + 4*M, prefer M*M + 3*M + M*NB)
 *  *
                   CALL DORGBR( 'Q', M, M, M, U, LDU, WORK( ITAUQ ),                    CALL DORGBR( 'Q', M, M, M, U, LDU, WORK( ITAUQ ),
      $                         WORK( IWORK ), LWORK-IWORK+1, IERR )       $                         WORK( IWORK ), LWORK-IWORK+1, IERR )
Line 2325 Line 2333
 *                 Perform bidiagonal QR iteration, computing left  *                 Perform bidiagonal QR iteration, computing left
 *                 singular vectors of L in U, and computing right  *                 singular vectors of L in U, and computing right
 *                 singular vectors of L in WORK(IR)  *                 singular vectors of L in WORK(IR)
 *                 (Workspace: need M*M+BDSPAC)  *                 (Workspace: need M*M + BDSPAC)
 *  *
                   CALL DBDSQR( 'U', M, M, M, 0, S, WORK( IE ),                    CALL DBDSQR( 'U', M, M, M, 0, S, WORK( IE ),
      $                         WORK( IR ), LDWRKR, U, LDU, DUM, 1,       $                         WORK( IR ), LDWRKR, U, LDU, DUM, 1,
Line 2334 Line 2342
 *  *
 *                 Multiply right singular vectors of L in WORK(IR) by Q  *                 Multiply right singular vectors of L in WORK(IR) by Q
 *                 in A, storing result in WORK(IU) and copying to A  *                 in A, storing result in WORK(IU) and copying to A
 *                 (Workspace: need M*M+2*M, prefer M*M+M*N+M))  *                 (Workspace: need M*M + 2*M, prefer M*M + M*N + M))
 *  *
                   DO 40 I = 1, N, CHUNK                    DO 40 I = 1, N, CHUNK
                      BLK = MIN( N-I+1, CHUNK )                       BLK = MIN( N-I+1, CHUNK )
Line 2353 Line 2361
                   IWORK = ITAU + M                    IWORK = ITAU + M
 *  *
 *                 Compute A=L*Q  *                 Compute A=L*Q
 *                 (Workspace: need 2*M, prefer M+M*NB)  *                 (Workspace: need 2*M, prefer M + M*NB)
 *  *
                   CALL DGELQF( M, N, A, LDA, WORK( ITAU ),                    CALL DGELQF( M, N, A, LDA, WORK( ITAU ),
      $                         WORK( IWORK ), LWORK-IWORK+1, IERR )       $                         WORK( IWORK ), LWORK-IWORK+1, IERR )
Line 2365 Line 2373
      $                         LDU )       $                         LDU )
 *  *
 *                 Generate Q in A  *                 Generate Q in A
 *                 (Workspace: need 2*M, prefer M+M*NB)  *                 (Workspace: need 2*M, prefer M + M*NB)
 *  *
                   CALL DORGLQ( M, N, M, A, LDA, WORK( ITAU ),                    CALL DORGLQ( M, N, M, A, LDA, WORK( ITAU ),
      $                         WORK( IWORK ), LWORK-IWORK+1, IERR )       $                         WORK( IWORK ), LWORK-IWORK+1, IERR )
Line 2375 Line 2383
                   IWORK = ITAUP + M                    IWORK = ITAUP + M
 *  *
 *                 Bidiagonalize L in U  *                 Bidiagonalize L in U
 *                 (Workspace: need 4*M, prefer 3*M+2*M*NB)  *                 (Workspace: need 4*M, prefer 3*M + 2*M*NB)
 *  *
                   CALL DGEBRD( M, M, U, LDU, S, WORK( IE ),                    CALL DGEBRD( M, M, U, LDU, S, WORK( IE ),
      $                         WORK( ITAUQ ), WORK( ITAUP ),       $                         WORK( ITAUQ ), WORK( ITAUP ),
      $                         WORK( IWORK ), LWORK-IWORK+1, IERR )       $                         WORK( IWORK ), LWORK-IWORK+1, IERR )
 *  *
 *                 Multiply right vectors bidiagonalizing L by Q in A  *                 Multiply right vectors bidiagonalizing L by Q in A
 *                 (Workspace: need 3*M+N, prefer 3*M+N*NB)  *                 (Workspace: need 3*M + N, prefer 3*M + N*NB)
 *  *
                   CALL DORMBR( 'P', 'L', 'T', M, N, M, U, LDU,                    CALL DORMBR( 'P', 'L', 'T', M, N, M, U, LDU,
      $                         WORK( ITAUP ), A, LDA, WORK( IWORK ),       $                         WORK( ITAUP ), A, LDA, WORK( IWORK ),
      $                         LWORK-IWORK+1, IERR )       $                         LWORK-IWORK+1, IERR )
 *  *
 *                 Generate left vectors bidiagonalizing L in U  *                 Generate left vectors bidiagonalizing L in U
 *                 (Workspace: need 4*M, prefer 3*M+M*NB)  *                 (Workspace: need 4*M, prefer 3*M + M*NB)
 *  *
                   CALL DORGBR( 'Q', M, M, M, U, LDU, WORK( ITAUQ ),                    CALL DORGBR( 'Q', M, M, M, U, LDU, WORK( ITAUQ ),
      $                         WORK( IWORK ), LWORK-IWORK+1, IERR )       $                         WORK( IWORK ), LWORK-IWORK+1, IERR )
Line 2433 Line 2441
                      IWORK = ITAU + M                       IWORK = ITAU + M
 *  *
 *                    Compute A=L*Q  *                    Compute A=L*Q
 *                    (Workspace: need M*M+2*M, prefer M*M+M+M*NB)  *                    (Workspace: need M*M + 2*M, prefer M*M + M + M*NB)
 *  *
                      CALL DGELQF( M, N, A, LDA, WORK( ITAU ),                       CALL DGELQF( M, N, A, LDA, WORK( ITAU ),
      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )       $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
Line 2446 Line 2454
      $                            WORK( IR+LDWRKR ), LDWRKR )       $                            WORK( IR+LDWRKR ), LDWRKR )
 *  *
 *                    Generate Q in A  *                    Generate Q in A
 *                    (Workspace: need M*M+2*M, prefer M*M+M+M*NB)  *                    (Workspace: need M*M + 2*M, prefer M*M + M + M*NB)
 *  *
                      CALL DORGLQ( M, N, M, A, LDA, WORK( ITAU ),                       CALL DORGLQ( M, N, M, A, LDA, WORK( ITAU ),
      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )       $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
Line 2456 Line 2464
                      IWORK = ITAUP + M                       IWORK = ITAUP + M
 *  *
 *                    Bidiagonalize L in WORK(IR)  *                    Bidiagonalize L in WORK(IR)
 *                    (Workspace: need M*M+4*M, prefer M*M+3*M+2*M*NB)  *                    (Workspace: need M*M + 4*M, prefer M*M + 3*M + 2*M*NB)
 *  *
                      CALL DGEBRD( M, M, WORK( IR ), LDWRKR, S,                       CALL DGEBRD( M, M, WORK( IR ), LDWRKR, S,
      $                            WORK( IE ), WORK( ITAUQ ),       $                            WORK( IE ), WORK( ITAUQ ),
Line 2465 Line 2473
 *  *
 *                    Generate right vectors bidiagonalizing L in  *                    Generate right vectors bidiagonalizing L in
 *                    WORK(IR)  *                    WORK(IR)
 *                    (Workspace: need M*M+4*M, prefer M*M+3*M+(M-1)*NB)  *                    (Workspace: need M*M + 4*M, prefer M*M + 3*M + (M-1)*NB)
 *  *
                      CALL DORGBR( 'P', M, M, M, WORK( IR ), LDWRKR,                       CALL DORGBR( 'P', M, M, M, WORK( IR ), LDWRKR,
      $                            WORK( ITAUP ), WORK( IWORK ),       $                            WORK( ITAUP ), WORK( IWORK ),
Line 2474 Line 2482
 *  *
 *                    Perform bidiagonal QR iteration, computing right  *                    Perform bidiagonal QR iteration, computing right
 *                    singular vectors of L in WORK(IR)  *                    singular vectors of L in WORK(IR)
 *                    (Workspace: need M*M+BDSPAC)  *                    (Workspace: need M*M + BDSPAC)
 *  *
                      CALL DBDSQR( 'U', M, M, 0, 0, S, WORK( IE ),                       CALL DBDSQR( 'U', M, M, 0, 0, S, WORK( IE ),
      $                            WORK( IR ), LDWRKR, DUM, 1, DUM, 1,       $                            WORK( IR ), LDWRKR, DUM, 1, DUM, 1,
Line 2495 Line 2503
                      IWORK = ITAU + M                       IWORK = ITAU + M
 *  *
 *                    Compute A=L*Q  *                    Compute A=L*Q
 *                    (Workspace: need 2*M, prefer M+M*NB)  *                    (Workspace: need 2*M, prefer M + M*NB)
 *  *
                      CALL DGELQF( M, N, A, LDA, WORK( ITAU ),                       CALL DGELQF( M, N, A, LDA, WORK( ITAU ),
      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )       $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
Line 2505 Line 2513
                      CALL DLACPY( 'U', M, N, A, LDA, VT, LDVT )                       CALL DLACPY( 'U', M, N, A, LDA, VT, LDVT )
 *  *
 *                    Generate Q in VT  *                    Generate Q in VT
 *                    (Workspace: need 2*M, prefer M+M*NB)  *                    (Workspace: need 2*M, prefer M + M*NB)
 *  *
                      CALL DORGLQ( M, N, M, VT, LDVT, WORK( ITAU ),                       CALL DORGLQ( M, N, M, VT, LDVT, WORK( ITAU ),
      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )       $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
Line 2520 Line 2528
      $                            LDA )       $                            LDA )
 *  *
 *                    Bidiagonalize L in A  *                    Bidiagonalize L in A
 *                    (Workspace: need 4*M, prefer 3*M+2*M*NB)  *                    (Workspace: need 4*M, prefer 3*M + 2*M*NB)
 *  *
                      CALL DGEBRD( M, M, A, LDA, S, WORK( IE ),                       CALL DGEBRD( M, M, A, LDA, S, WORK( IE ),
      $                            WORK( ITAUQ ), WORK( ITAUP ),       $                            WORK( ITAUQ ), WORK( ITAUP ),
      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )       $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
 *  *
 *                    Multiply right vectors bidiagonalizing L by Q in VT  *                    Multiply right vectors bidiagonalizing L by Q in VT
 *                    (Workspace: need 3*M+N, prefer 3*M+N*NB)  *                    (Workspace: need 3*M + N, prefer 3*M + N*NB)
 *  *
                      CALL DORMBR( 'P', 'L', 'T', M, N, M, A, LDA,                       CALL DORMBR( 'P', 'L', 'T', M, N, M, A, LDA,
      $                            WORK( ITAUP ), VT, LDVT,       $                            WORK( ITAUP ), VT, LDVT,
Line 2562 Line 2570
                         LDWRKU = LDA                          LDWRKU = LDA
                         IR = IU + LDWRKU*M                          IR = IU + LDWRKU*M
                         LDWRKR = LDA                          LDWRKR = LDA
                      ELSE IF( LWORK.GE.WRKBL+( LDA+M )*M ) THEN                       ELSE IF( LWORK.GE.WRKBL+( LDA + M )*M ) THEN
 *  *
 *                       WORK(IU) is LDA by M and WORK(IR) is M by M  *                       WORK(IU) is LDA by M and WORK(IR) is M by M
 *  *
Line 2581 Line 2589
                      IWORK = ITAU + M                       IWORK = ITAU + M
 *  *
 *                    Compute A=L*Q  *                    Compute A=L*Q
 *                    (Workspace: need 2*M*M+2*M, prefer 2*M*M+M+M*NB)  *                    (Workspace: need 2*M*M + 2*M, prefer 2*M*M + M + M*NB)
 *  *
                      CALL DGELQF( M, N, A, LDA, WORK( ITAU ),                       CALL DGELQF( M, N, A, LDA, WORK( ITAU ),
      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )       $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
Line 2594 Line 2602
      $                            WORK( IU+LDWRKU ), LDWRKU )       $                            WORK( IU+LDWRKU ), LDWRKU )
 *  *
 *                    Generate Q in A  *                    Generate Q in A
 *                    (Workspace: need 2*M*M+2*M, prefer 2*M*M+M+M*NB)  *                    (Workspace: need 2*M*M + 2*M, prefer 2*M*M + M + M*NB)
 *  *
                      CALL DORGLQ( M, N, M, A, LDA, WORK( ITAU ),                       CALL DORGLQ( M, N, M, A, LDA, WORK( ITAU ),
      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )       $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
Line 2605 Line 2613
 *  *
 *                    Bidiagonalize L in WORK(IU), copying result to  *                    Bidiagonalize L in WORK(IU), copying result to
 *                    WORK(IR)  *                    WORK(IR)
 *                    (Workspace: need 2*M*M+4*M,  *                    (Workspace: need 2*M*M + 4*M,
 *                                prefer 2*M*M+3*M+2*M*NB)  *                                prefer 2*M*M+3*M+2*M*NB)
 *  *
                      CALL DGEBRD( M, M, WORK( IU ), LDWRKU, S,                       CALL DGEBRD( M, M, WORK( IU ), LDWRKU, S,
Line 2616 Line 2624
      $                            WORK( IR ), LDWRKR )       $                            WORK( IR ), LDWRKR )
 *  *
 *                    Generate right bidiagonalizing vectors in WORK(IU)  *                    Generate right bidiagonalizing vectors in WORK(IU)
 *                    (Workspace: need 2*M*M+4*M-1,  *                    (Workspace: need 2*M*M + 4*M-1,
 *                                prefer 2*M*M+3*M+(M-1)*NB)  *                                prefer 2*M*M+3*M+(M-1)*NB)
 *  *
                      CALL DORGBR( 'P', M, M, M, WORK( IU ), LDWRKU,                       CALL DORGBR( 'P', M, M, M, WORK( IU ), LDWRKU,
Line 2624 Line 2632
      $                            LWORK-IWORK+1, IERR )       $                            LWORK-IWORK+1, IERR )
 *  *
 *                    Generate left bidiagonalizing vectors in WORK(IR)  *                    Generate left bidiagonalizing vectors in WORK(IR)
 *                    (Workspace: need 2*M*M+4*M, prefer 2*M*M+3*M+M*NB)  *                    (Workspace: need 2*M*M + 4*M, prefer 2*M*M + 3*M + M*NB)
 *  *
                      CALL DORGBR( 'Q', M, M, M, WORK( IR ), LDWRKR,                       CALL DORGBR( 'Q', M, M, M, WORK( IR ), LDWRKR,
      $                            WORK( ITAUQ ), WORK( IWORK ),       $                            WORK( ITAUQ ), WORK( IWORK ),
Line 2634 Line 2642
 *                    Perform bidiagonal QR iteration, computing left  *                    Perform bidiagonal QR iteration, computing left
 *                    singular vectors of L in WORK(IR) and computing  *                    singular vectors of L in WORK(IR) and computing
 *                    right singular vectors of L in WORK(IU)  *                    right singular vectors of L in WORK(IU)
 *                    (Workspace: need 2*M*M+BDSPAC)  *                    (Workspace: need 2*M*M + BDSPAC)
 *  *
                      CALL DBDSQR( 'U', M, M, M, 0, S, WORK( IE ),                       CALL DBDSQR( 'U', M, M, M, 0, S, WORK( IE ),
      $                            WORK( IU ), LDWRKU, WORK( IR ),       $                            WORK( IU ), LDWRKU, WORK( IR ),
Line 2661 Line 2669
                      IWORK = ITAU + M                       IWORK = ITAU + M
 *  *
 *                    Compute A=L*Q, copying result to VT  *                    Compute A=L*Q, copying result to VT
 *                    (Workspace: need 2*M, prefer M+M*NB)  *                    (Workspace: need 2*M, prefer M + M*NB)
 *  *
                      CALL DGELQF( M, N, A, LDA, WORK( ITAU ),                       CALL DGELQF( M, N, A, LDA, WORK( ITAU ),
      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )       $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
                      CALL DLACPY( 'U', M, N, A, LDA, VT, LDVT )                       CALL DLACPY( 'U', M, N, A, LDA, VT, LDVT )
 *  *
 *                    Generate Q in VT  *                    Generate Q in VT
 *                    (Workspace: need 2*M, prefer M+M*NB)  *                    (Workspace: need 2*M, prefer M + M*NB)
 *  *
                      CALL DORGLQ( M, N, M, VT, LDVT, WORK( ITAU ),                       CALL DORGLQ( M, N, M, VT, LDVT, WORK( ITAU ),
      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )       $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
Line 2683 Line 2691
      $                            LDA )       $                            LDA )
 *  *
 *                    Bidiagonalize L in A  *                    Bidiagonalize L in A
 *                    (Workspace: need 4*M, prefer 3*M+2*M*NB)  *                    (Workspace: need 4*M, prefer 3*M + 2*M*NB)
 *  *
                      CALL DGEBRD( M, M, A, LDA, S, WORK( IE ),                       CALL DGEBRD( M, M, A, LDA, S, WORK( IE ),
      $                            WORK( ITAUQ ), WORK( ITAUP ),       $                            WORK( ITAUQ ), WORK( ITAUP ),
      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )       $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
 *  *
 *                    Multiply right vectors bidiagonalizing L by Q in VT  *                    Multiply right vectors bidiagonalizing L by Q in VT
 *                    (Workspace: need 3*M+N, prefer 3*M+N*NB)  *                    (Workspace: need 3*M + N, prefer 3*M + N*NB)
 *  *
                      CALL DORMBR( 'P', 'L', 'T', M, N, M, A, LDA,                       CALL DORMBR( 'P', 'L', 'T', M, N, M, A, LDA,
      $                            WORK( ITAUP ), VT, LDVT,       $                            WORK( ITAUP ), VT, LDVT,
      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )       $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
 *  *
 *                    Generate left bidiagonalizing vectors of L in A  *                    Generate left bidiagonalizing vectors of L in A
 *                    (Workspace: need 4*M, prefer 3*M+M*NB)  *                    (Workspace: need 4*M, prefer 3*M + M*NB)
 *  *
                      CALL DORGBR( 'Q', M, M, M, A, LDA, WORK( ITAUQ ),                       CALL DORGBR( 'Q', M, M, M, A, LDA, WORK( ITAUQ ),
      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )       $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
Line 2741 Line 2749
                      IWORK = ITAU + M                       IWORK = ITAU + M
 *  *
 *                    Compute A=L*Q  *                    Compute A=L*Q
 *                    (Workspace: need M*M+2*M, prefer M*M+M+M*NB)  *                    (Workspace: need M*M + 2*M, prefer M*M + M + M*NB)
 *  *
                      CALL DGELQF( M, N, A, LDA, WORK( ITAU ),                       CALL DGELQF( M, N, A, LDA, WORK( ITAU ),
      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )       $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
Line 2754 Line 2762
      $                            WORK( IU+LDWRKU ), LDWRKU )       $                            WORK( IU+LDWRKU ), LDWRKU )
 *  *
 *                    Generate Q in A  *                    Generate Q in A
 *                    (Workspace: need M*M+2*M, prefer M*M+M+M*NB)  *                    (Workspace: need M*M + 2*M, prefer M*M + M + M*NB)
 *  *
                      CALL DORGLQ( M, N, M, A, LDA, WORK( ITAU ),                       CALL DORGLQ( M, N, M, A, LDA, WORK( ITAU ),
      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )       $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
Line 2764 Line 2772
                      IWORK = ITAUP + M                       IWORK = ITAUP + M
 *  *
 *                    Bidiagonalize L in WORK(IU), copying result to U  *                    Bidiagonalize L in WORK(IU), copying result to U
 *                    (Workspace: need M*M+4*M, prefer M*M+3*M+2*M*NB)  *                    (Workspace: need M*M + 4*M, prefer M*M + 3*M + 2*M*NB)
 *  *
                      CALL DGEBRD( M, M, WORK( IU ), LDWRKU, S,                       CALL DGEBRD( M, M, WORK( IU ), LDWRKU, S,
      $                            WORK( IE ), WORK( ITAUQ ),       $                            WORK( IE ), WORK( ITAUQ ),
Line 2774 Line 2782
      $                            LDU )       $                            LDU )
 *  *
 *                    Generate right bidiagonalizing vectors in WORK(IU)  *                    Generate right bidiagonalizing vectors in WORK(IU)
 *                    (Workspace: need M*M+4*M-1,  *                    (Workspace: need M*M + 4*M-1,
 *                                prefer M*M+3*M+(M-1)*NB)  *                                prefer M*M+3*M+(M-1)*NB)
 *  *
                      CALL DORGBR( 'P', M, M, M, WORK( IU ), LDWRKU,                       CALL DORGBR( 'P', M, M, M, WORK( IU ), LDWRKU,
Line 2782 Line 2790
      $                            LWORK-IWORK+1, IERR )       $                            LWORK-IWORK+1, IERR )
 *  *
 *                    Generate left bidiagonalizing vectors in U  *                    Generate left bidiagonalizing vectors in U
 *                    (Workspace: need M*M+4*M, prefer M*M+3*M+M*NB)  *                    (Workspace: need M*M + 4*M, prefer M*M + 3*M + M*NB)
 *  *
                      CALL DORGBR( 'Q', M, M, M, U, LDU, WORK( ITAUQ ),                       CALL DORGBR( 'Q', M, M, M, U, LDU, WORK( ITAUQ ),
      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )       $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
Line 2791 Line 2799
 *                    Perform bidiagonal QR iteration, computing left  *                    Perform bidiagonal QR iteration, computing left
 *                    singular vectors of L in U and computing right  *                    singular vectors of L in U and computing right
 *                    singular vectors of L in WORK(IU)  *                    singular vectors of L in WORK(IU)
 *                    (Workspace: need M*M+BDSPAC)  *                    (Workspace: need M*M + BDSPAC)
 *  *
                      CALL DBDSQR( 'U', M, M, M, 0, S, WORK( IE ),                       CALL DBDSQR( 'U', M, M, M, 0, S, WORK( IE ),
      $                            WORK( IU ), LDWRKU, U, LDU, DUM, 1,       $                            WORK( IU ), LDWRKU, U, LDU, DUM, 1,
Line 2812 Line 2820
                      IWORK = ITAU + M                       IWORK = ITAU + M
 *  *
 *                    Compute A=L*Q, copying result to VT  *                    Compute A=L*Q, copying result to VT
 *                    (Workspace: need 2*M, prefer M+M*NB)  *                    (Workspace: need 2*M, prefer M + M*NB)
 *  *
                      CALL DGELQF( M, N, A, LDA, WORK( ITAU ),                       CALL DGELQF( M, N, A, LDA, WORK( ITAU ),
      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )       $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
                      CALL DLACPY( 'U', M, N, A, LDA, VT, LDVT )                       CALL DLACPY( 'U', M, N, A, LDA, VT, LDVT )
 *  *
 *                    Generate Q in VT  *                    Generate Q in VT
 *                    (Workspace: need 2*M, prefer M+M*NB)  *                    (Workspace: need 2*M, prefer M + M*NB)
 *  *
                      CALL DORGLQ( M, N, M, VT, LDVT, WORK( ITAU ),                       CALL DORGLQ( M, N, M, VT, LDVT, WORK( ITAU ),
      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )       $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
Line 2835 Line 2843
                      IWORK = ITAUP + M                       IWORK = ITAUP + M
 *  *
 *                    Bidiagonalize L in U  *                    Bidiagonalize L in U
 *                    (Workspace: need 4*M, prefer 3*M+2*M*NB)  *                    (Workspace: need 4*M, prefer 3*M + 2*M*NB)
 *  *
                      CALL DGEBRD( M, M, U, LDU, S, WORK( IE ),                       CALL DGEBRD( M, M, U, LDU, S, WORK( IE ),
      $                            WORK( ITAUQ ), WORK( ITAUP ),       $                            WORK( ITAUQ ), WORK( ITAUP ),
Line 2843 Line 2851
 *  *
 *                    Multiply right bidiagonalizing vectors in U by Q  *                    Multiply right bidiagonalizing vectors in U by Q
 *                    in VT  *                    in VT
 *                    (Workspace: need 3*M+N, prefer 3*M+N*NB)  *                    (Workspace: need 3*M + N, prefer 3*M + N*NB)
 *  *
                      CALL DORMBR( 'P', 'L', 'T', M, N, M, U, LDU,                       CALL DORMBR( 'P', 'L', 'T', M, N, M, U, LDU,
      $                            WORK( ITAUP ), VT, LDVT,       $                            WORK( ITAUP ), VT, LDVT,
      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )       $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
 *  *
 *                    Generate left bidiagonalizing vectors in U  *                    Generate left bidiagonalizing vectors in U
 *                    (Workspace: need 4*M, prefer 3*M+M*NB)  *                    (Workspace: need 4*M, prefer 3*M + M*NB)
 *  *
                      CALL DORGBR( 'Q', M, M, M, U, LDU, WORK( ITAUQ ),                       CALL DORGBR( 'Q', M, M, M, U, LDU, WORK( ITAUQ ),
      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )       $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
Line 2877 Line 2885
 *                 N right singular vectors to be computed in VT and  *                 N right singular vectors to be computed in VT and
 *                 no left singular vectors to be computed  *                 no left singular vectors to be computed
 *  *
                   IF( LWORK.GE.M*M+MAX( N+M, 4*M, BDSPAC ) ) THEN                    IF( LWORK.GE.M*M+MAX( N + M, 4*M, BDSPAC ) ) THEN
 *  *
 *                    Sufficient workspace for a fast algorithm  *                    Sufficient workspace for a fast algorithm
 *  *
Line 2897 Line 2905
                      IWORK = ITAU + M                       IWORK = ITAU + M
 *  *
 *                    Compute A=L*Q, copying result to VT  *                    Compute A=L*Q, copying result to VT
 *                    (Workspace: need M*M+2*M, prefer M*M+M+M*NB)  *                    (Workspace: need M*M + 2*M, prefer M*M + M + M*NB)
 *  *
                      CALL DGELQF( M, N, A, LDA, WORK( ITAU ),                       CALL DGELQF( M, N, A, LDA, WORK( ITAU ),
      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )       $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
Line 2911 Line 2919
      $                            WORK( IR+LDWRKR ), LDWRKR )       $                            WORK( IR+LDWRKR ), LDWRKR )
 *  *
 *                    Generate Q in VT  *                    Generate Q in VT
 *                    (Workspace: need M*M+M+N, prefer M*M+M+N*NB)  *                    (Workspace: need M*M + M + N, prefer M*M + M + N*NB)
 *  *
                      CALL DORGLQ( N, N, M, VT, LDVT, WORK( ITAU ),                       CALL DORGLQ( N, N, M, VT, LDVT, WORK( ITAU ),
      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )       $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
Line 2921 Line 2929
                      IWORK = ITAUP + M                       IWORK = ITAUP + M
 *  *
 *                    Bidiagonalize L in WORK(IR)  *                    Bidiagonalize L in WORK(IR)
 *                    (Workspace: need M*M+4*M, prefer M*M+3*M+2*M*NB)  *                    (Workspace: need M*M + 4*M, prefer M*M + 3*M + 2*M*NB)
 *  *
                      CALL DGEBRD( M, M, WORK( IR ), LDWRKR, S,                       CALL DGEBRD( M, M, WORK( IR ), LDWRKR, S,
      $                            WORK( IE ), WORK( ITAUQ ),       $                            WORK( IE ), WORK( ITAUQ ),
Line 2929 Line 2937
      $                            LWORK-IWORK+1, IERR )       $                            LWORK-IWORK+1, IERR )
 *  *
 *                    Generate right bidiagonalizing vectors in WORK(IR)  *                    Generate right bidiagonalizing vectors in WORK(IR)
 *                    (Workspace: need M*M+4*M-1,  *                    (Workspace: need M*M + 4*M-1,
 *                                prefer M*M+3*M+(M-1)*NB)  *                                prefer M*M+3*M+(M-1)*NB)
 *  *
                      CALL DORGBR( 'P', M, M, M, WORK( IR ), LDWRKR,                       CALL DORGBR( 'P', M, M, M, WORK( IR ), LDWRKR,
Line 2939 Line 2947
 *  *
 *                    Perform bidiagonal QR iteration, computing right  *                    Perform bidiagonal QR iteration, computing right
 *                    singular vectors of L in WORK(IR)  *                    singular vectors of L in WORK(IR)
 *                    (Workspace: need M*M+BDSPAC)  *                    (Workspace: need M*M + BDSPAC)
 *  *
                      CALL DBDSQR( 'U', M, M, 0, 0, S, WORK( IE ),                       CALL DBDSQR( 'U', M, M, 0, 0, S, WORK( IE ),
      $                            WORK( IR ), LDWRKR, DUM, 1, DUM, 1,       $                            WORK( IR ), LDWRKR, DUM, 1, DUM, 1,
Line 2964 Line 2972
                      IWORK = ITAU + M                       IWORK = ITAU + M
 *  *
 *                    Compute A=L*Q, copying result to VT  *                    Compute A=L*Q, copying result to VT
 *                    (Workspace: need 2*M, prefer M+M*NB)  *                    (Workspace: need 2*M, prefer M + M*NB)
 *  *
                      CALL DGELQF( M, N, A, LDA, WORK( ITAU ),                       CALL DGELQF( M, N, A, LDA, WORK( ITAU ),
      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )       $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
                      CALL DLACPY( 'U', M, N, A, LDA, VT, LDVT )                       CALL DLACPY( 'U', M, N, A, LDA, VT, LDVT )
 *  *
 *                    Generate Q in VT  *                    Generate Q in VT
 *                    (Workspace: need M+N, prefer M+N*NB)  *                    (Workspace: need M + N, prefer M + N*NB)
 *  *
                      CALL DORGLQ( N, N, M, VT, LDVT, WORK( ITAU ),                       CALL DORGLQ( N, N, M, VT, LDVT, WORK( ITAU ),
      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )       $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
Line 2986 Line 2994
      $                            LDA )       $                            LDA )
 *  *
 *                    Bidiagonalize L in A  *                    Bidiagonalize L in A
 *                    (Workspace: need 4*M, prefer 3*M+2*M*NB)  *                    (Workspace: need 4*M, prefer 3*M + 2*M*NB)
 *  *
                      CALL DGEBRD( M, M, A, LDA, S, WORK( IE ),                       CALL DGEBRD( M, M, A, LDA, S, WORK( IE ),
      $                            WORK( ITAUQ ), WORK( ITAUP ),       $                            WORK( ITAUQ ), WORK( ITAUP ),
Line 2994 Line 3002
 *  *
 *                    Multiply right bidiagonalizing vectors in A by Q  *                    Multiply right bidiagonalizing vectors in A by Q
 *                    in VT  *                    in VT
 *                    (Workspace: need 3*M+N, prefer 3*M+N*NB)  *                    (Workspace: need 3*M + N, prefer 3*M + N*NB)
 *  *
                      CALL DORMBR( 'P', 'L', 'T', M, N, M, A, LDA,                       CALL DORMBR( 'P', 'L', 'T', M, N, M, A, LDA,
      $                            WORK( ITAUP ), VT, LDVT,       $                            WORK( ITAUP ), VT, LDVT,
Line 3017 Line 3025
 *                 N right singular vectors to be computed in VT and  *                 N right singular vectors to be computed in VT and
 *                 M left singular vectors to be overwritten on A  *                 M left singular vectors to be overwritten on A
 *  *
                   IF( LWORK.GE.2*M*M+MAX( N+M, 4*M, BDSPAC ) ) THEN                    IF( LWORK.GE.2*M*M+MAX( N + M, 4*M, BDSPAC ) ) THEN
 *  *
 *                    Sufficient workspace for a fast algorithm  *                    Sufficient workspace for a fast algorithm
 *  *
Line 3029 Line 3037
                         LDWRKU = LDA                          LDWRKU = LDA
                         IR = IU + LDWRKU*M                          IR = IU + LDWRKU*M
                         LDWRKR = LDA                          LDWRKR = LDA
                      ELSE IF( LWORK.GE.WRKBL+( LDA+M )*M ) THEN                       ELSE IF( LWORK.GE.WRKBL+( LDA + M )*M ) THEN
 *  *
 *                       WORK(IU) is LDA by M and WORK(IR) is M by M  *                       WORK(IU) is LDA by M and WORK(IR) is M by M
 *  *
Line 3048 Line 3056
                      IWORK = ITAU + M                       IWORK = ITAU + M
 *  *
 *                    Compute A=L*Q, copying result to VT  *                    Compute A=L*Q, copying result to VT
 *                    (Workspace: need 2*M*M+2*M, prefer 2*M*M+M+M*NB)  *                    (Workspace: need 2*M*M + 2*M, prefer 2*M*M + M + M*NB)
 *  *
                      CALL DGELQF( M, N, A, LDA, WORK( ITAU ),                       CALL DGELQF( M, N, A, LDA, WORK( ITAU ),
      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )       $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
                      CALL DLACPY( 'U', M, N, A, LDA, VT, LDVT )                       CALL DLACPY( 'U', M, N, A, LDA, VT, LDVT )
 *  *
 *                    Generate Q in VT  *                    Generate Q in VT
 *                    (Workspace: need 2*M*M+M+N, prefer 2*M*M+M+N*NB)  *                    (Workspace: need 2*M*M + M + N, prefer 2*M*M + M + N*NB)
 *  *
                      CALL DORGLQ( N, N, M, VT, LDVT, WORK( ITAU ),                       CALL DORGLQ( N, N, M, VT, LDVT, WORK( ITAU ),
      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )       $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
Line 3073 Line 3081
 *  *
 *                    Bidiagonalize L in WORK(IU), copying result to  *                    Bidiagonalize L in WORK(IU), copying result to
 *                    WORK(IR)  *                    WORK(IR)
 *                    (Workspace: need 2*M*M+4*M,  *                    (Workspace: need 2*M*M + 4*M,
 *                                prefer 2*M*M+3*M+2*M*NB)  *                                prefer 2*M*M+3*M+2*M*NB)
 *  *
                      CALL DGEBRD( M, M, WORK( IU ), LDWRKU, S,                       CALL DGEBRD( M, M, WORK( IU ), LDWRKU, S,
Line 3084 Line 3092
      $                            WORK( IR ), LDWRKR )       $                            WORK( IR ), LDWRKR )
 *  *
 *                    Generate right bidiagonalizing vectors in WORK(IU)  *                    Generate right bidiagonalizing vectors in WORK(IU)
 *                    (Workspace: need 2*M*M+4*M-1,  *                    (Workspace: need 2*M*M + 4*M-1,
 *                                prefer 2*M*M+3*M+(M-1)*NB)  *                                prefer 2*M*M+3*M+(M-1)*NB)
 *  *
                      CALL DORGBR( 'P', M, M, M, WORK( IU ), LDWRKU,                       CALL DORGBR( 'P', M, M, M, WORK( IU ), LDWRKU,
Line 3092 Line 3100
      $                            LWORK-IWORK+1, IERR )       $                            LWORK-IWORK+1, IERR )
 *  *
 *                    Generate left bidiagonalizing vectors in WORK(IR)  *                    Generate left bidiagonalizing vectors in WORK(IR)
 *                    (Workspace: need 2*M*M+4*M, prefer 2*M*M+3*M+M*NB)  *                    (Workspace: need 2*M*M + 4*M, prefer 2*M*M + 3*M + M*NB)
 *  *
                      CALL DORGBR( 'Q', M, M, M, WORK( IR ), LDWRKR,                       CALL DORGBR( 'Q', M, M, M, WORK( IR ), LDWRKR,
      $                            WORK( ITAUQ ), WORK( IWORK ),       $                            WORK( ITAUQ ), WORK( IWORK ),
Line 3102 Line 3110
 *                    Perform bidiagonal QR iteration, computing left  *                    Perform bidiagonal QR iteration, computing left
 *                    singular vectors of L in WORK(IR) and computing  *                    singular vectors of L in WORK(IR) and computing
 *                    right singular vectors of L in WORK(IU)  *                    right singular vectors of L in WORK(IU)
 *                    (Workspace: need 2*M*M+BDSPAC)  *                    (Workspace: need 2*M*M + BDSPAC)
 *  *
                      CALL DBDSQR( 'U', M, M, M, 0, S, WORK( IE ),                       CALL DBDSQR( 'U', M, M, M, 0, S, WORK( IE ),
      $                            WORK( IU ), LDWRKU, WORK( IR ),       $                            WORK( IU ), LDWRKU, WORK( IR ),
Line 3132 Line 3140
                      IWORK = ITAU + M                       IWORK = ITAU + M
 *  *
 *                    Compute A=L*Q, copying result to VT  *                    Compute A=L*Q, copying result to VT
 *                    (Workspace: need 2*M, prefer M+M*NB)  *                    (Workspace: need 2*M, prefer M + M*NB)
 *  *
                      CALL DGELQF( M, N, A, LDA, WORK( ITAU ),                       CALL DGELQF( M, N, A, LDA, WORK( ITAU ),
      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )       $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
                      CALL DLACPY( 'U', M, N, A, LDA, VT, LDVT )                       CALL DLACPY( 'U', M, N, A, LDA, VT, LDVT )
 *  *
 *                    Generate Q in VT  *                    Generate Q in VT
 *                    (Workspace: need M+N, prefer M+N*NB)  *                    (Workspace: need M + N, prefer M + N*NB)
 *  *
                      CALL DORGLQ( N, N, M, VT, LDVT, WORK( ITAU ),                       CALL DORGLQ( N, N, M, VT, LDVT, WORK( ITAU ),
      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )       $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
Line 3154 Line 3162
      $                            LDA )       $                            LDA )
 *  *
 *                    Bidiagonalize L in A  *                    Bidiagonalize L in A
 *                    (Workspace: need 4*M, prefer 3*M+2*M*NB)  *                    (Workspace: need 4*M, prefer 3*M + 2*M*NB)
 *  *
                      CALL DGEBRD( M, M, A, LDA, S, WORK( IE ),                       CALL DGEBRD( M, M, A, LDA, S, WORK( IE ),
      $                            WORK( ITAUQ ), WORK( ITAUP ),       $                            WORK( ITAUQ ), WORK( ITAUP ),
Line 3162 Line 3170
 *  *
 *                    Multiply right bidiagonalizing vectors in A by Q  *                    Multiply right bidiagonalizing vectors in A by Q
 *                    in VT  *                    in VT
 *                    (Workspace: need 3*M+N, prefer 3*M+N*NB)  *                    (Workspace: need 3*M + N, prefer 3*M + N*NB)
 *  *
                      CALL DORMBR( 'P', 'L', 'T', M, N, M, A, LDA,                       CALL DORMBR( 'P', 'L', 'T', M, N, M, A, LDA,
      $                            WORK( ITAUP ), VT, LDVT,       $                            WORK( ITAUP ), VT, LDVT,
      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )       $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
 *  *
 *                    Generate left bidiagonalizing vectors in A  *                    Generate left bidiagonalizing vectors in A
 *                    (Workspace: need 4*M, prefer 3*M+M*NB)  *                    (Workspace: need 4*M, prefer 3*M + M*NB)
 *  *
                      CALL DORGBR( 'Q', M, M, M, A, LDA, WORK( ITAUQ ),                       CALL DORGBR( 'Q', M, M, M, A, LDA, WORK( ITAUQ ),
      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )       $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
Line 3193 Line 3201
 *                 N right singular vectors to be computed in VT and  *                 N right singular vectors to be computed in VT and
 *                 M left singular vectors to be computed in U  *                 M left singular vectors to be computed in U
 *  *
                   IF( LWORK.GE.M*M+MAX( N+M, 4*M, BDSPAC ) ) THEN                    IF( LWORK.GE.M*M+MAX( N + M, 4*M, BDSPAC ) ) THEN
 *  *
 *                    Sufficient workspace for a fast algorithm  *                    Sufficient workspace for a fast algorithm
 *  *
Line 3213 Line 3221
                      IWORK = ITAU + M                       IWORK = ITAU + M
 *  *
 *                    Compute A=L*Q, copying result to VT  *                    Compute A=L*Q, copying result to VT
 *                    (Workspace: need M*M+2*M, prefer M*M+M+M*NB)  *                    (Workspace: need M*M + 2*M, prefer M*M + M + M*NB)
 *  *
                      CALL DGELQF( M, N, A, LDA, WORK( ITAU ),                       CALL DGELQF( M, N, A, LDA, WORK( ITAU ),
      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )       $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
                      CALL DLACPY( 'U', M, N, A, LDA, VT, LDVT )                       CALL DLACPY( 'U', M, N, A, LDA, VT, LDVT )
 *  *
 *                    Generate Q in VT  *                    Generate Q in VT
 *                    (Workspace: need M*M+M+N, prefer M*M+M+N*NB)  *                    (Workspace: need M*M + M + N, prefer M*M + M + N*NB)
 *  *
                      CALL DORGLQ( N, N, M, VT, LDVT, WORK( ITAU ),                       CALL DORGLQ( N, N, M, VT, LDVT, WORK( ITAU ),
      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )       $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
Line 3237 Line 3245
                      IWORK = ITAUP + M                       IWORK = ITAUP + M
 *  *
 *                    Bidiagonalize L in WORK(IU), copying result to U  *                    Bidiagonalize L in WORK(IU), copying result to U
 *                    (Workspace: need M*M+4*M, prefer M*M+3*M+2*M*NB)  *                    (Workspace: need M*M + 4*M, prefer M*M + 3*M + 2*M*NB)
 *  *
                      CALL DGEBRD( M, M, WORK( IU ), LDWRKU, S,                       CALL DGEBRD( M, M, WORK( IU ), LDWRKU, S,
      $                            WORK( IE ), WORK( ITAUQ ),       $                            WORK( IE ), WORK( ITAUQ ),
Line 3247 Line 3255
      $                            LDU )       $                            LDU )
 *  *
 *                    Generate right bidiagonalizing vectors in WORK(IU)  *                    Generate right bidiagonalizing vectors in WORK(IU)
 *                    (Workspace: need M*M+4*M, prefer M*M+3*M+(M-1)*NB)  *                    (Workspace: need M*M + 4*M, prefer M*M + 3*M + (M-1)*NB)
 *  *
                      CALL DORGBR( 'P', M, M, M, WORK( IU ), LDWRKU,                       CALL DORGBR( 'P', M, M, M, WORK( IU ), LDWRKU,
      $                            WORK( ITAUP ), WORK( IWORK ),       $                            WORK( ITAUP ), WORK( IWORK ),
      $                            LWORK-IWORK+1, IERR )       $                            LWORK-IWORK+1, IERR )
 *  *
 *                    Generate left bidiagonalizing vectors in U  *                    Generate left bidiagonalizing vectors in U
 *                    (Workspace: need M*M+4*M, prefer M*M+3*M+M*NB)  *                    (Workspace: need M*M + 4*M, prefer M*M + 3*M + M*NB)
 *  *
                      CALL DORGBR( 'Q', M, M, M, U, LDU, WORK( ITAUQ ),                       CALL DORGBR( 'Q', M, M, M, U, LDU, WORK( ITAUQ ),
      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )       $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
Line 3263 Line 3271
 *                    Perform bidiagonal QR iteration, computing left  *                    Perform bidiagonal QR iteration, computing left
 *                    singular vectors of L in U and computing right  *                    singular vectors of L in U and computing right
 *                    singular vectors of L in WORK(IU)  *                    singular vectors of L in WORK(IU)
 *                    (Workspace: need M*M+BDSPAC)  *                    (Workspace: need M*M + BDSPAC)
 *  *
                      CALL DBDSQR( 'U', M, M, M, 0, S, WORK( IE ),                       CALL DBDSQR( 'U', M, M, M, 0, S, WORK( IE ),
      $                            WORK( IU ), LDWRKU, U, LDU, DUM, 1,       $                            WORK( IU ), LDWRKU, U, LDU, DUM, 1,
Line 3288 Line 3296
                      IWORK = ITAU + M                       IWORK = ITAU + M
 *  *
 *                    Compute A=L*Q, copying result to VT  *                    Compute A=L*Q, copying result to VT
 *                    (Workspace: need 2*M, prefer M+M*NB)  *                    (Workspace: need 2*M, prefer M + M*NB)
 *  *
                      CALL DGELQF( M, N, A, LDA, WORK( ITAU ),                       CALL DGELQF( M, N, A, LDA, WORK( ITAU ),
      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )       $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
                      CALL DLACPY( 'U', M, N, A, LDA, VT, LDVT )                       CALL DLACPY( 'U', M, N, A, LDA, VT, LDVT )
 *  *
 *                    Generate Q in VT  *                    Generate Q in VT
 *                    (Workspace: need M+N, prefer M+N*NB)  *                    (Workspace: need M + N, prefer M + N*NB)
 *  *
                      CALL DORGLQ( N, N, M, VT, LDVT, WORK( ITAU ),                       CALL DORGLQ( N, N, M, VT, LDVT, WORK( ITAU ),
      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )       $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
Line 3311 Line 3319
                      IWORK = ITAUP + M                       IWORK = ITAUP + M
 *  *
 *                    Bidiagonalize L in U  *                    Bidiagonalize L in U
 *                    (Workspace: need 4*M, prefer 3*M+2*M*NB)  *                    (Workspace: need 4*M, prefer 3*M + 2*M*NB)
 *  *
                      CALL DGEBRD( M, M, U, LDU, S, WORK( IE ),                       CALL DGEBRD( M, M, U, LDU, S, WORK( IE ),
      $                            WORK( ITAUQ ), WORK( ITAUP ),       $                            WORK( ITAUQ ), WORK( ITAUP ),
Line 3319 Line 3327
 *  *
 *                    Multiply right bidiagonalizing vectors in U by Q  *                    Multiply right bidiagonalizing vectors in U by Q
 *                    in VT  *                    in VT
 *                    (Workspace: need 3*M+N, prefer 3*M+N*NB)  *                    (Workspace: need 3*M + N, prefer 3*M + N*NB)
 *  *
                      CALL DORMBR( 'P', 'L', 'T', M, N, M, U, LDU,                       CALL DORMBR( 'P', 'L', 'T', M, N, M, U, LDU,
      $                            WORK( ITAUP ), VT, LDVT,       $                            WORK( ITAUP ), VT, LDVT,
      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )       $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
 *  *
 *                    Generate left bidiagonalizing vectors in U  *                    Generate left bidiagonalizing vectors in U
 *                    (Workspace: need 4*M, prefer 3*M+M*NB)  *                    (Workspace: need 4*M, prefer 3*M + M*NB)
 *  *
                      CALL DORGBR( 'Q', M, M, M, U, LDU, WORK( ITAUQ ),                       CALL DORGBR( 'Q', M, M, M, U, LDU, WORK( ITAUQ ),
      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )       $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
Line 3360 Line 3368
             IWORK = ITAUP + M              IWORK = ITAUP + M
 *  *
 *           Bidiagonalize A  *           Bidiagonalize A
 *           (Workspace: need 3*M+N, prefer 3*M+(M+N)*NB)  *           (Workspace: need 3*M + N, prefer 3*M + (M + N)*NB)
 *  *
             CALL DGEBRD( M, N, A, LDA, S, WORK( IE ), WORK( ITAUQ ),              CALL DGEBRD( M, N, A, LDA, S, WORK( IE ), WORK( ITAUQ ),
      $                   WORK( ITAUP ), WORK( IWORK ), LWORK-IWORK+1,       $                   WORK( ITAUP ), WORK( IWORK ), LWORK-IWORK+1,
Line 3369 Line 3377
 *  *
 *              If left singular vectors desired in U, copy result to U  *              If left singular vectors desired in U, copy result to U
 *              and generate left bidiagonalizing vectors in U  *              and generate left bidiagonalizing vectors in U
 *              (Workspace: need 4*M-1, prefer 3*M+(M-1)*NB)  *              (Workspace: need 4*M-1, prefer 3*M + (M-1)*NB)
 *  *
                CALL DLACPY( 'L', M, M, A, LDA, U, LDU )                 CALL DLACPY( 'L', M, M, A, LDA, U, LDU )
                CALL DORGBR( 'Q', M, M, N, U, LDU, WORK( ITAUQ ),                 CALL DORGBR( 'Q', M, M, N, U, LDU, WORK( ITAUQ ),
Line 3379 Line 3387
 *  *
 *              If right singular vectors desired in VT, copy result to  *              If right singular vectors desired in VT, copy result to
 *              VT and generate right bidiagonalizing vectors in VT  *              VT and generate right bidiagonalizing vectors in VT
 *              (Workspace: need 3*M+NRVT, prefer 3*M+NRVT*NB)  *              (Workspace: need 3*M + NRVT, prefer 3*M + NRVT*NB)
 *  *
                CALL DLACPY( 'U', M, N, A, LDA, VT, LDVT )                 CALL DLACPY( 'U', M, N, A, LDA, VT, LDVT )
                IF( WNTVA )                 IF( WNTVA )
Line 3393 Line 3401
 *  *
 *              If left singular vectors desired in A, generate left  *              If left singular vectors desired in A, generate left
 *              bidiagonalizing vectors in A  *              bidiagonalizing vectors in A
 *              (Workspace: need 4*M-1, prefer 3*M+(M-1)*NB)  *              (Workspace: need 4*M-1, prefer 3*M + (M-1)*NB)
 *  *
                CALL DORGBR( 'Q', M, M, N, A, LDA, WORK( ITAUQ ),                 CALL DORGBR( 'Q', M, M, N, A, LDA, WORK( ITAUQ ),
      $                      WORK( IWORK ), LWORK-IWORK+1, IERR )       $                      WORK( IWORK ), LWORK-IWORK+1, IERR )
Line 3402 Line 3410
 *  *
 *              If right singular vectors desired in A, generate right  *              If right singular vectors desired in A, generate right
 *              bidiagonalizing vectors in A  *              bidiagonalizing vectors in A
 *              (Workspace: need 4*M, prefer 3*M+M*NB)  *              (Workspace: need 4*M, prefer 3*M + M*NB)
 *  *
                CALL DORGBR( 'P', M, N, M, A, LDA, WORK( ITAUP ),                 CALL DORGBR( 'P', M, N, M, A, LDA, WORK( ITAUP ),
      $                      WORK( IWORK ), LWORK-IWORK+1, IERR )       $                      WORK( IWORK ), LWORK-IWORK+1, IERR )

Removed from v.1.10  
changed lines
  Added in v.1.20


CVSweb interface <joel.bertrand@systella.fr>