[BACK]Return to dgesvdx.f CVS log [TXT] [DIR] Up to [local] / rpl / lapack / lapack

Diff for /rpl/lapack/lapack/dgesvdx.f between versions 1.3 and 1.8

version 1.3, 2016/08/27 15:34:22 version 1.8, 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 DGESVDX + dependencies   *> Download DGESVDX + dependencies
 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dgesvdx.f">   *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dgesvdx.f">
 *> [TGZ]</a>   *> [TGZ]</a>
 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dgesvdx.f">   *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dgesvdx.f">
 *> [ZIP]</a>   *> [ZIP]</a>
 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dgesvdx.f">   *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dgesvdx.f">
 *> [TXT]</a>  *> [TXT]</a>
 *> \endhtmlonly   *> \endhtmlonly
 *  *
 *  Definition:  *  Definition:
 *  ===========  *  ===========
 *  *
 *     SUBROUTINE DGESVDX( JOBU, JOBVT, RANGE, M, N, A, LDA, VL, VU,   *     SUBROUTINE DGESVDX( JOBU, JOBVT, RANGE, M, N, A, LDA, VL, VU,
 *    $                    IL, IU, NS, S, U, LDU, VT, LDVT, WORK,   *    $                    IL, IU, NS, S, U, LDU, VT, LDVT, WORK,
 *    $                    LWORK, IWORK, INFO )  *    $                    LWORK, IWORK, INFO )
 *        *
 *  *
 *     .. Scalar Arguments ..  *     .. Scalar Arguments ..
 *      CHARACTER          JOBU, JOBVT, RANGE  *      CHARACTER          JOBU, JOBVT, RANGE
Line 33 Line 33
 *     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 43 Line 43
 *>  DGESVDX computes the singular value decomposition (SVD) of a real  *>  DGESVDX computes the singular value decomposition (SVD) of a real
 *>  M-by-N matrix A, optionally computing the left and/or right singular  *>  M-by-N matrix A, optionally computing the left and/or right singular
 *>  vectors. The SVD is written  *>  vectors. The SVD is written
 *>   *>
 *>      A = U * SIGMA * transpose(V)  *>      A = U * SIGMA * transpose(V)
 *>   *>
 *>  where SIGMA is an M-by-N matrix which is zero except for its  *>  where SIGMA is an M-by-N matrix which is zero except for its
 *>  min(m,n) diagonal elements, U is an M-by-M orthogonal matrix, and  *>  min(m,n) diagonal elements, U is an M-by-M orthogonal matrix, and
 *>  V is an N-by-N orthogonal matrix.  The diagonal elements of SIGMA  *>  V is an N-by-N orthogonal matrix.  The diagonal elements of SIGMA
 *>  are the singular values of A; they are real and non-negative, and  *>  are the singular values of A; they are real and non-negative, and
 *>  are returned in descending order.  The first min(m,n) columns of  *>  are returned in descending order.  The first min(m,n) columns of
 *>  U and V are the left and right singular vectors of A.  *>  U and V are the left and right singular vectors of A.
 *>   *>
 *>  DGESVDX uses an eigenvalue problem for obtaining the SVD, which   *>  DGESVDX uses an eigenvalue problem for obtaining the SVD, which
 *>  allows for the computation of a subset of singular values and   *>  allows for the computation of a subset of singular values and
 *>  vectors. See DBDSVDX for details.  *>  vectors. See DBDSVDX for details.
 *>   *>
 *>  Note that the routine returns V**T, not V.  *>  Note that the routine returns V**T, not V.
 *> \endverbatim  *> \endverbatim
 *     *
 *  Arguments:  *  Arguments:
 *  ==========  *  ==========
 *  *
Line 68 Line 68
 *>          JOBU is CHARACTER*1  *>          JOBU is CHARACTER*1
 *>          Specifies options for computing all or part of the matrix U:  *>          Specifies options for computing all or part of the matrix U:
 *>          = 'V':  the first min(m,n) columns of U (the left singular  *>          = 'V':  the first min(m,n) columns of U (the left singular
 *>                  vectors) or as specified by RANGE are returned in   *>                  vectors) or as specified by RANGE are returned in
 *>                  the array U;  *>                  the array U;
 *>          = 'N':  no columns of U (no left singular vectors) are  *>          = 'N':  no columns of U (no left singular vectors) are
 *>                  computed.  *>                  computed.
Line 80 Line 80
 *>           Specifies options for computing all or part of the matrix  *>           Specifies options for computing all or part of the matrix
 *>           V**T:  *>           V**T:
 *>           = 'V':  the first min(m,n) rows of V**T (the right singular  *>           = 'V':  the first min(m,n) rows of V**T (the right singular
 *>                   vectors) or as specified by RANGE are returned in   *>                   vectors) or as specified by RANGE are returned in
 *>                   the array VT;  *>                   the array VT;
 *>           = 'N':  no rows of V**T (no right singular vectors) are  *>           = 'N':  no rows of V**T (no right singular vectors) are
 *>                   computed.  *>                   computed.
Line 92 Line 92
 *>          = 'A': all singular values will be found.  *>          = 'A': all singular values will be found.
 *>          = 'V': all singular values in the half-open interval (VL,VU]  *>          = 'V': all singular values in the half-open interval (VL,VU]
 *>                 will be found.  *>                 will be found.
 *>          = 'I': the IL-th through IU-th singular values will be found.   *>          = 'I': the IL-th through IU-th singular values will be found.
 *> \endverbatim  *> \endverbatim
 *>  *>
 *> \param[in] M  *> \param[in] M
Line 157 Line 157
 *> \param[out] NS  *> \param[out] NS
 *> \verbatim  *> \verbatim
 *>          NS is INTEGER  *>          NS is INTEGER
 *>          The total number of singular values found,    *>          The total number of singular values found,
 *>          0 <= NS <= min(M,N).  *>          0 <= NS <= min(M,N).
 *>          If RANGE = 'A', NS = min(M,N); if RANGE = 'I', NS = IU-IL+1.  *>          If RANGE = 'A', NS = min(M,N); if RANGE = 'I', NS = IU-IL+1.
 *> \endverbatim  *> \endverbatim
Line 171 Line 171
 *> \param[out] U  *> \param[out] U
 *> \verbatim  *> \verbatim
 *>          U is DOUBLE PRECISION array, dimension (LDU,UCOL)  *>          U is DOUBLE PRECISION array, dimension (LDU,UCOL)
 *>          If JOBU = 'V', U contains columns of U (the left singular   *>          If JOBU = 'V', U contains columns of U (the left singular
 *>          vectors, stored columnwise) as specified by RANGE; if   *>          vectors, stored columnwise) as specified by RANGE; if
 *>          JOBU = 'N', U is not referenced.  *>          JOBU = 'N', U is not referenced.
 *>          Note: The user must ensure that UCOL >= NS; if RANGE = 'V',   *>          Note: The user must ensure that UCOL >= NS; if RANGE = 'V',
 *>          the exact value of NS is not known in advance and an upper  *>          the exact value of NS is not known in advance and an upper
 *>          bound must be used.  *>          bound must be used.
 *> \endverbatim  *> \endverbatim
Line 189 Line 189
 *> \param[out] VT  *> \param[out] VT
 *> \verbatim  *> \verbatim
 *>          VT is DOUBLE PRECISION array, dimension (LDVT,N)  *>          VT is DOUBLE PRECISION array, dimension (LDVT,N)
 *>          If JOBVT = 'V', VT contains the rows of V**T (the right singular   *>          If JOBVT = 'V', VT contains the rows of V**T (the right singular
 *>          vectors, stored rowwise) as specified by RANGE; if JOBVT = 'N',   *>          vectors, stored rowwise) as specified by RANGE; if JOBVT = 'N',
 *>          VT is not referenced.  *>          VT is not referenced.
 *>          Note: The user must ensure that LDVT >= NS; if RANGE = 'V',   *>          Note: The user must ensure that LDVT >= NS; if RANGE = 'V',
 *>          the exact value of NS is not known in advance and an upper   *>          the exact value of NS is not known in advance and an upper
 *>          bound must be used.  *>          bound must be used.
 *> \endverbatim  *> \endverbatim
 *>  *>
Line 214 Line 214
 *> \verbatim  *> \verbatim
 *>          LWORK is INTEGER  *>          LWORK is INTEGER
 *>          The dimension of the array WORK.  *>          The dimension of the array WORK.
 *>          LWORK >= MAX(1,MIN(M,N)*(MIN(M,N)+4)) for the paths (see   *>          LWORK >= MAX(1,MIN(M,N)*(MIN(M,N)+4)) for the paths (see
 *>          comments inside the code):  *>          comments inside the code):
 *>             - PATH 1  (M much larger than N)   *>             - PATH 1  (M much larger than N)
 *>             - PATH 1t (N much larger than M)  *>             - PATH 1t (N much larger than M)
 *>          LWORK >= MAX(1,MIN(M,N)*2+MAX(M,N)) for the other paths.  *>          LWORK >= MAX(1,MIN(M,N)*2+MAX(M,N)) for the other paths.
 *>          For good performance, LWORK should generally be larger.  *>          For good performance, LWORK should generally be larger.
Line 230 Line 230
 *> \param[out] IWORK  *> \param[out] IWORK
 *> \verbatim  *> \verbatim
 *>          IWORK is INTEGER array, dimension (12*MIN(M,N))  *>          IWORK is INTEGER array, dimension (12*MIN(M,N))
 *>          If INFO = 0, the first NS elements of IWORK are zero. If INFO > 0,   *>          If INFO = 0, the first NS elements of IWORK are zero. If INFO > 0,
 *>          then IWORK contains the indices of the eigenvectors that failed   *>          then IWORK contains the indices of the eigenvectors that failed
 *>          to converge in DBDSVDX/DSTEVX.  *>          to converge in DBDSVDX/DSTEVX.
 *> \endverbatim  *> \endverbatim
 *>  *>
Line 249 Line 249
 *  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 June 2016  
 *  *
 *> \ingroup doubleGEsing  *> \ingroup doubleGEsing
 *  *
 *  =====================================================================  *  =====================================================================
       SUBROUTINE DGESVDX( JOBU, JOBVT, RANGE, M, N, A, LDA, VL, VU,         SUBROUTINE DGESVDX( JOBU, JOBVT, RANGE, M, N, A, LDA, VL, VU,
      $                    IL, IU, NS, S, U, LDU, VT, LDVT, WORK,        $                    IL, IU, NS, S, U, LDU, VT, LDVT, WORK,
      $                    LWORK, IWORK, INFO )       $                    LWORK, IWORK, INFO )
 *  *
 *  -- LAPACK driver routine (version 3.6.1) --  *  -- 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..--
 *     June 2016  
 *  *
 *     .. Scalar Arguments ..  *     .. Scalar Arguments ..
       CHARACTER          JOBU, JOBVT, RANGE        CHARACTER          JOBU, JOBVT, RANGE
Line 289 Line 286
       CHARACTER          JOBZ, RNGTGK        CHARACTER          JOBZ, RNGTGK
       LOGICAL            ALLS, INDS, LQUERY, VALS, WANTU, WANTVT        LOGICAL            ALLS, INDS, LQUERY, VALS, WANTU, WANTVT
       INTEGER            I, ID, IE, IERR, ILQF, ILTGK, IQRF, ISCL,        INTEGER            I, ID, IE, IERR, ILQF, ILTGK, IQRF, ISCL,
      $                   ITAU, ITAUP, ITAUQ, ITEMP, ITGKZ, IUTGK,        $                   ITAU, ITAUP, ITAUQ, ITEMP, ITGKZ, IUTGK,
      $                   J, MAXWRK, MINMN, MINWRK, MNTHR       $                   J, MAXWRK, MINMN, MINWRK, MNTHR
       DOUBLE PRECISION   ABSTOL, ANRM, BIGNUM, EPS, SMLNUM        DOUBLE PRECISION   ABSTOL, ANRM, BIGNUM, EPS, SMLNUM
 *     ..  *     ..
Line 299 Line 296
 *     .. External Subroutines ..  *     .. External Subroutines ..
       EXTERNAL           DBDSVDX, DGEBRD, DGELQF, DGEQRF, DLACPY,        EXTERNAL           DBDSVDX, DGEBRD, DGELQF, DGEQRF, DLACPY,
      $                   DLASCL, DLASET, DORMBR, DORMLQ, DORMQR,       $                   DLASCL, DLASET, DORMBR, DORMLQ, DORMQR,
      $                   DSCAL, XERBLA       $                   DCOPY, XERBLA
 *     ..  *     ..
 *     .. External Functions ..  *     .. External Functions ..
       LOGICAL            LSAME        LOGICAL            LSAME
       INTEGER            ILAENV        INTEGER            ILAENV
       DOUBLE PRECISION   DLAMCH, DLANGE, DNRM2        DOUBLE PRECISION   DLAMCH, DLANGE
       EXTERNAL           LSAME, ILAENV, DLAMCH, DLANGE, DNRM2        EXTERNAL           LSAME, ILAENV, DLAMCH, DLANGE
 *     ..  *     ..
 *     .. Intrinsic Functions ..  *     .. Intrinsic Functions ..
       INTRINSIC          MAX, MIN, SQRT        INTRINSIC          MAX, MIN, SQRT
Line 392 Line 389
 *  *
 *                 Path 1 (M much larger than N)  *                 Path 1 (M much larger than N)
 *  *
                   MAXWRK = N +                     MAXWRK = N +
      $                     N*ILAENV( 1, 'DGEQRF', ' ', M, N, -1, -1 )       $                     N*ILAENV( 1, 'DGEQRF', ' ', M, N, -1, -1 )
                   MAXWRK = MAX( MAXWRK, N*(N+5) + 2*N*                    MAXWRK = MAX( MAXWRK, N*(N+5) + 2*N*
      $                     ILAENV( 1, 'DGEBRD', ' ', N, N, -1, -1 ) )       $                     ILAENV( 1, 'DGEBRD', ' ', N, N, -1, -1 ) )
Line 427 Line 424
 *  *
 *                 Path 1t (N much larger than M)  *                 Path 1t (N much larger than M)
 *  *
                   MAXWRK = M +                     MAXWRK = M +
      $                     M*ILAENV( 1, 'DGELQF', ' ', M, N, -1, -1 )       $                     M*ILAENV( 1, 'DGELQF', ' ', M, N, -1, -1 )
                   MAXWRK = MAX( MAXWRK, M*(M+5) + 2*M*                    MAXWRK = MAX( MAXWRK, M*(M+5) + 2*M*
      $                     ILAENV( 1, 'DGEBRD', ' ', M, M, -1, -1 ) )       $                     ILAENV( 1, 'DGEBRD', ' ', M, M, -1, -1 ) )
Line 489 Line 486
          RNGTGK = 'I'           RNGTGK = 'I'
          ILTGK = IL           ILTGK = IL
          IUTGK = IU           IUTGK = IU
       ELSE              ELSE
          RNGTGK = 'V'           RNGTGK = 'V'
          ILTGK = 0           ILTGK = 0
          IUTGK = 0           IUTGK = 0
Line 533 Line 530
             ITEMP = ITAU + N              ITEMP = ITAU + N
             CALL DGEQRF( M, N, A, LDA, WORK( ITAU ), WORK( ITEMP ),              CALL DGEQRF( M, N, A, LDA, WORK( ITAU ), WORK( ITEMP ),
      $                   LWORK-ITEMP+1, INFO )       $                   LWORK-ITEMP+1, INFO )
 *    *
 *           Copy R into WORK and bidiagonalize it:  *           Copy R into WORK and bidiagonalize it:
 *           (Workspace: need N*N+5*N, prefer N*N+4*N+2*N*NB)  *           (Workspace: need N*N+5*N, prefer N*N+4*N+2*N*NB)
 *  *
Line 542 Line 539
             IE = ID + N              IE = ID + N
             ITAUQ = IE + N              ITAUQ = IE + N
             ITAUP = ITAUQ + N              ITAUP = ITAUQ + N
             ITEMP = ITAUP + N                           ITEMP = ITAUP + N
             CALL DLACPY( 'U', N, N, A, LDA, WORK( IQRF ), N )              CALL DLACPY( 'U', N, N, A, LDA, WORK( IQRF ), N )
             CALL DLASET( 'L', N-1, N-1, ZERO, ZERO, WORK( IQRF+1 ), N )              CALL DLASET( 'L', N-1, N-1, ZERO, ZERO, WORK( IQRF+1 ), N )
             CALL DGEBRD( N, N, WORK( IQRF ), N, WORK( ID ), WORK( IE ),               CALL DGEBRD( N, N, WORK( IQRF ), N, WORK( ID ), WORK( IE ),
      $                   WORK( ITAUQ ), WORK( ITAUP ), WORK( ITEMP ),       $                   WORK( ITAUQ ), WORK( ITAUP ), WORK( ITEMP ),
      $                   LWORK-ITEMP+1, INFO )       $                   LWORK-ITEMP+1, INFO )
 *  *
 *           Solve eigenvalue problem TGK*Z=Z*S.  *           Solve eigenvalue problem TGK*Z=Z*S.
 *           (Workspace: need 14*N + 2*N*(N+1))            *           (Workspace: need 14*N + 2*N*(N+1))
 *              *
             ITGKZ = ITEMP              ITGKZ = ITEMP
             ITEMP = ITGKZ + N*(N*2+1)              ITEMP = ITGKZ + N*(N*2+1)
             CALL DBDSVDX( 'U', JOBZ, RNGTGK, N, WORK( ID ), WORK( IE ),               CALL DBDSVDX( 'U', JOBZ, RNGTGK, N, WORK( ID ), WORK( IE ),
      $                    VL, VU, ILTGK, IUTGK, NS, S, WORK( ITGKZ ),       $                    VL, VU, ILTGK, IUTGK, NS, S, WORK( ITGKZ ),
      $                    N*2, WORK( ITEMP ), IWORK, INFO)       $                    N*2, WORK( ITEMP ), IWORK, INFO)
 *  *
Line 571 Line 568
 *              Call DORMBR to compute QB*UB.  *              Call DORMBR to compute QB*UB.
 *              (Workspace in WORK( ITEMP ): need N, prefer N*NB)  *              (Workspace in WORK( ITEMP ): need N, prefer N*NB)
 *  *
                CALL DORMBR( 'Q', 'L', 'N', N, NS, N, WORK( IQRF ), N,                  CALL DORMBR( 'Q', 'L', 'N', N, NS, N, WORK( IQRF ), N,
      $                      WORK( ITAUQ ), U, LDU, WORK( ITEMP ),        $                      WORK( ITAUQ ), U, LDU, WORK( ITEMP ),
      $                      LWORK-ITEMP+1, INFO )       $                      LWORK-ITEMP+1, INFO )
 *  *
 *              Call DORMQR to compute Q*(QB*UB).  *              Call DORMQR to compute Q*(QB*UB).
 *              (Workspace in WORK( ITEMP ): need N, prefer N*NB)  *              (Workspace in WORK( ITEMP ): need N, prefer N*NB)
 *  *
                CALL DORMQR( 'L', 'N', M, NS, N, A, LDA,                  CALL DORMQR( 'L', 'N', M, NS, N, A, LDA,
      $                      WORK( ITAU ), U, LDU, WORK( ITEMP ),       $                      WORK( ITAU ), U, LDU, WORK( ITEMP ),
      $                      LWORK-ITEMP+1, INFO )       $                      LWORK-ITEMP+1, INFO )
             END IF                END IF
 *        *
 *           If needed, compute right singular vectors.  *           If needed, compute right singular vectors.
 *  *
             IF( WANTVT) THEN              IF( WANTVT) THEN
Line 595 Line 592
 *              Call DORMBR to compute VB**T * PB**T  *              Call DORMBR to compute VB**T * PB**T
 *              (Workspace in WORK( ITEMP ): need N, prefer N*NB)  *              (Workspace in WORK( ITEMP ): need N, prefer N*NB)
 *  *
                CALL DORMBR( 'P', 'R', 'T', NS, N, N, WORK( IQRF ), N,                  CALL DORMBR( 'P', 'R', 'T', NS, N, N, WORK( IQRF ), N,
      $                      WORK( ITAUP ), VT, LDVT, WORK( ITEMP ),       $                      WORK( ITAUP ), VT, LDVT, WORK( ITEMP ),
      $                      LWORK-ITEMP+1, INFO )       $                      LWORK-ITEMP+1, INFO )
             END IF              END IF
Line 613 Line 610
             IE = ID + N              IE = ID + N
             ITAUQ = IE + N              ITAUQ = IE + N
             ITAUP = ITAUQ + N              ITAUP = ITAUQ + N
             ITEMP = ITAUP + N                        ITEMP = ITAUP + N
             CALL DGEBRD( M, N, A, LDA, WORK( ID ), WORK( IE ),               CALL DGEBRD( M, N, A, LDA, WORK( ID ), WORK( IE ),
      $                   WORK( ITAUQ ), WORK( ITAUP ), WORK( ITEMP ),       $                   WORK( ITAUQ ), WORK( ITAUP ), WORK( ITEMP ),
      $                   LWORK-ITEMP+1, INFO )       $                   LWORK-ITEMP+1, INFO )
 *  *
 *           Solve eigenvalue problem TGK*Z=Z*S.  *           Solve eigenvalue problem TGK*Z=Z*S.
 *           (Workspace: need 14*N + 2*N*(N+1))            *           (Workspace: need 14*N + 2*N*(N+1))
 *             *
             ITGKZ = ITEMP              ITGKZ = ITEMP
             ITEMP = ITGKZ + N*(N*2+1)              ITEMP = ITGKZ + N*(N*2+1)
             CALL DBDSVDX( 'U', JOBZ, RNGTGK, N, WORK( ID ), WORK( IE ),               CALL DBDSVDX( 'U', JOBZ, RNGTGK, N, WORK( ID ), WORK( IE ),
      $                    VL, VU, ILTGK, IUTGK, NS, S, WORK( ITGKZ ),       $                    VL, VU, ILTGK, IUTGK, NS, S, WORK( ITGKZ ),
      $                    N*2, WORK( ITEMP ), IWORK, INFO)       $                    N*2, WORK( ITEMP ), IWORK, INFO)
 *  *
Line 639 Line 636
 *  *
 *              Call DORMBR to compute QB*UB.  *              Call DORMBR to compute QB*UB.
 *              (Workspace in WORK( ITEMP ): need N, prefer N*NB)  *              (Workspace in WORK( ITEMP ): need N, prefer N*NB)
 *     *
                CALL DORMBR( 'Q', 'L', 'N', M, NS, N, A, LDA,                  CALL DORMBR( 'Q', 'L', 'N', M, NS, N, A, LDA,
      $                      WORK( ITAUQ ), U, LDU, WORK( ITEMP ),        $                      WORK( ITAUQ ), U, LDU, WORK( ITEMP ),
      $                      LWORK-ITEMP+1, IERR )       $                      LWORK-ITEMP+1, IERR )
             END IF                END IF
 *        *
 *           If needed, compute right singular vectors.  *           If needed, compute right singular vectors.
 *  *
             IF( WANTVT) THEN              IF( WANTVT) THEN
Line 657 Line 654
 *              Call DORMBR to compute VB**T * PB**T  *              Call DORMBR to compute VB**T * PB**T
 *              (Workspace in WORK( ITEMP ): need N, prefer N*NB)  *              (Workspace in WORK( ITEMP ): need N, prefer N*NB)
 *  *
                CALL DORMBR( 'P', 'R', 'T', NS, N, N, A, LDA,                  CALL DORMBR( 'P', 'R', 'T', NS, N, N, A, LDA,
      $                      WORK( ITAUP ), VT, LDVT, WORK( ITEMP ),       $                      WORK( ITAUP ), VT, LDVT, WORK( ITEMP ),
      $                      LWORK-ITEMP+1, IERR )       $                      LWORK-ITEMP+1, IERR )
             END IF              END IF
          END IF                                        END IF
       ELSE        ELSE
 *  *
 *        A has more columns than rows. If A has sufficiently more  *        A has more columns than rows. If A has sufficiently more
Line 670 Line 667
          IF( N.GE.MNTHR ) THEN           IF( N.GE.MNTHR ) THEN
 *  *
 *           Path 1t (N much larger than M):  *           Path 1t (N much larger than M):
 *           A = L * Q = ( QB * B * PB**T ) * Q   *           A = L * Q = ( QB * B * PB**T ) * Q
 *                     = ( QB * ( UB * S * VB**T ) * PB**T ) * Q  *                     = ( QB * ( UB * S * VB**T ) * PB**T ) * Q
 *           U = QB * UB ; V**T = VB**T * PB**T * Q  *           U = QB * UB ; V**T = VB**T * PB**T * Q
 *  *
Line 693 Line 690
             ITEMP = ITAUP + M              ITEMP = ITAUP + M
             CALL DLACPY( 'L', M, M, A, LDA, WORK( ILQF ), M )              CALL DLACPY( 'L', M, M, A, LDA, WORK( ILQF ), M )
             CALL DLASET( 'U', M-1, M-1, ZERO, ZERO, WORK( ILQF+M ), M )              CALL DLASET( 'U', M-1, M-1, ZERO, ZERO, WORK( ILQF+M ), M )
             CALL DGEBRD( M, M, WORK( ILQF ), M, WORK( ID ), WORK( IE ),               CALL DGEBRD( M, M, WORK( ILQF ), M, WORK( ID ), WORK( IE ),
      $                   WORK( ITAUQ ), WORK( ITAUP ), WORK( ITEMP ),       $                   WORK( ITAUQ ), WORK( ITAUP ), WORK( ITEMP ),
      $                   LWORK-ITEMP+1, INFO )       $                   LWORK-ITEMP+1, INFO )
 *  *
 *           Solve eigenvalue problem TGK*Z=Z*S.  *           Solve eigenvalue problem TGK*Z=Z*S.
 *           (Workspace: need 2*M*M+14*M)            *           (Workspace: need 2*M*M+14*M)
 *  *
             ITGKZ = ITEMP              ITGKZ = ITEMP
             ITEMP = ITGKZ + M*(M*2+1)              ITEMP = ITGKZ + M*(M*2+1)
             CALL DBDSVDX( 'U', JOBZ, RNGTGK, M, WORK( ID ), WORK( IE ),               CALL DBDSVDX( 'U', JOBZ, RNGTGK, M, WORK( ID ), WORK( IE ),
      $                    VL, VU, ILTGK, IUTGK, NS, S, WORK( ITGKZ ),       $                    VL, VU, ILTGK, IUTGK, NS, S, WORK( ITGKZ ),
      $                    M*2, WORK( ITEMP ), IWORK, INFO)       $                    M*2, WORK( ITEMP ), IWORK, INFO)
 *  *
Line 718 Line 715
 *              Call DORMBR to compute QB*UB.  *              Call DORMBR to compute QB*UB.
 *              (Workspace in WORK( ITEMP ): need M, prefer M*NB)  *              (Workspace in WORK( ITEMP ): need M, prefer M*NB)
 *  *
                CALL DORMBR( 'Q', 'L', 'N', M, NS, M, WORK( ILQF ), M,                  CALL DORMBR( 'Q', 'L', 'N', M, NS, M, WORK( ILQF ), M,
      $                      WORK( ITAUQ ), U, LDU, WORK( ITEMP ),        $                      WORK( ITAUQ ), U, LDU, WORK( ITEMP ),
      $                      LWORK-ITEMP+1, INFO )       $                      LWORK-ITEMP+1, INFO )
             END IF                END IF
 *        *
 *           If needed, compute right singular vectors.  *           If needed, compute right singular vectors.
 *  *
             IF( WANTVT) THEN              IF( WANTVT) THEN
Line 736 Line 733
 *              Call DORMBR to compute (VB**T)*(PB**T)  *              Call DORMBR to compute (VB**T)*(PB**T)
 *              (Workspace in WORK( ITEMP ): need M, prefer M*NB)  *              (Workspace in WORK( ITEMP ): need M, prefer M*NB)
 *  *
                CALL DORMBR( 'P', 'R', 'T', NS, M, M, WORK( ILQF ), M,                  CALL DORMBR( 'P', 'R', 'T', NS, M, M, WORK( ILQF ), M,
      $                      WORK( ITAUP ), VT, LDVT, WORK( ITEMP ),       $                      WORK( ITAUP ), VT, LDVT, WORK( ITEMP ),
      $                      LWORK-ITEMP+1, INFO )       $                      LWORK-ITEMP+1, INFO )
 *  *
 *              Call DORMLQ to compute ((VB**T)*(PB**T))*Q.  *              Call DORMLQ to compute ((VB**T)*(PB**T))*Q.
 *              (Workspace in WORK( ITEMP ): need M, prefer M*NB)  *              (Workspace in WORK( ITEMP ): need M, prefer M*NB)
 *  *
                CALL DORMLQ( 'R', 'N', NS, N, M, A, LDA,                  CALL DORMLQ( 'R', 'N', NS, N, M, A, LDA,
      $                      WORK( ITAU ), VT, LDVT, WORK( ITEMP ),       $                      WORK( ITAU ), VT, LDVT, WORK( ITEMP ),
      $                      LWORK-ITEMP+1, INFO )       $                      LWORK-ITEMP+1, INFO )
             END IF                END IF
          ELSE           ELSE
 *  *
 *           Path 2t (N greater than M, but not much larger)  *           Path 2t (N greater than M, but not much larger)
 *           Reduce to bidiagonal form without LQ decomposition  *           Reduce to bidiagonal form without LQ decomposition
 *           A = QB * B * PB**T = QB * ( UB * S * VB**T ) * PB**T  *           A = QB * B * PB**T = QB * ( UB * S * VB**T ) * PB**T
 *           U = QB * UB; V**T = VB**T * PB**T             *           U = QB * UB; V**T = VB**T * PB**T
 *  *
 *           Bidiagonalize A  *           Bidiagonalize A
 *           (Workspace: need 4*M+N, prefer 4*M+(M+N)*NB)  *           (Workspace: need 4*M+N, prefer 4*M+(M+N)*NB)
Line 762 Line 759
             ITAUQ = IE + M              ITAUQ = IE + M
             ITAUP = ITAUQ + M              ITAUP = ITAUQ + M
             ITEMP = ITAUP + M              ITEMP = ITAUP + M
             CALL DGEBRD( M, N, A, LDA, WORK( ID ), WORK( IE ),               CALL DGEBRD( M, N, A, LDA, WORK( ID ), WORK( IE ),
      $                   WORK( ITAUQ ), WORK( ITAUP ), WORK( ITEMP ),       $                   WORK( ITAUQ ), WORK( ITAUP ), WORK( ITEMP ),
      $                   LWORK-ITEMP+1, INFO )       $                   LWORK-ITEMP+1, INFO )
 *  *
 *           Solve eigenvalue problem TGK*Z=Z*S.  *           Solve eigenvalue problem TGK*Z=Z*S.
 *           (Workspace: need 2*M*M+14*M)            *           (Workspace: need 2*M*M+14*M)
 *  *
             ITGKZ = ITEMP              ITGKZ = ITEMP
             ITEMP = ITGKZ + M*(M*2+1)              ITEMP = ITGKZ + M*(M*2+1)
             CALL DBDSVDX( 'L', JOBZ, RNGTGK, M, WORK( ID ), WORK( IE ),               CALL DBDSVDX( 'L', JOBZ, RNGTGK, M, WORK( ID ), WORK( IE ),
      $                    VL, VU, ILTGK, IUTGK, NS, S, WORK( ITGKZ ),       $                    VL, VU, ILTGK, IUTGK, NS, S, WORK( ITGKZ ),
      $                    M*2, WORK( ITEMP ), IWORK, INFO)       $                    M*2, WORK( ITEMP ), IWORK, INFO)
 *   *
 *           If needed, compute left singular vectors.  *           If needed, compute left singular vectors.
 *  *
             IF( WANTU ) THEN              IF( WANTU ) THEN
Line 787 Line 784
 *              Call DORMBR to compute QB*UB.  *              Call DORMBR to compute QB*UB.
 *              (Workspace in WORK( ITEMP ): need M, prefer M*NB)  *              (Workspace in WORK( ITEMP ): need M, prefer M*NB)
 *  *
                CALL DORMBR( 'Q', 'L', 'N', M, NS, N, A, LDA,                  CALL DORMBR( 'Q', 'L', 'N', M, NS, N, A, LDA,
      $                      WORK( ITAUQ ), U, LDU, WORK( ITEMP ),        $                      WORK( ITAUQ ), U, LDU, WORK( ITEMP ),
      $                      LWORK-ITEMP+1, INFO )       $                      LWORK-ITEMP+1, INFO )
             END IF                END IF
 *        *
 *           If needed, compute right singular vectors.  *           If needed, compute right singular vectors.
 *  *
             IF( WANTVT) THEN              IF( WANTVT) THEN
Line 805 Line 802
 *              Call DORMBR to compute VB**T * PB**T  *              Call DORMBR to compute VB**T * PB**T
 *              (Workspace in WORK( ITEMP ): need M, prefer M*NB)  *              (Workspace in WORK( ITEMP ): need M, prefer M*NB)
 *  *
                CALL DORMBR( 'P', 'R', 'T', NS, N, M, A, LDA,                  CALL DORMBR( 'P', 'R', 'T', NS, N, M, A, LDA,
      $                      WORK( ITAUP ), VT, LDVT, WORK( ITEMP ),       $                      WORK( ITAUP ), VT, LDVT, WORK( ITEMP ),
      $                      LWORK-ITEMP+1, INFO )       $                      LWORK-ITEMP+1, INFO )
             END IF               END IF
          END IF           END IF
       END IF        END IF
 *  *
Removed from v.1.3  
changed lines
  Added in v.1.8

CVSweb interface <joel.bertrand@systella.fr>