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

Diff for /rpl/lapack/lapack/zgesvdx.f between versions 1.3 and 1.4

version 1.3, 2016/08/27 15:34:47 version 1.4, 2017/06/17 10:54:11
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 ZGESVDX + dependencies   *> Download ZGESVDX + dependencies
 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zgesvdx.f">   *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zgesvdx.f">
 *> [TGZ]</a>   *> [TGZ]</a>
 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zgesvdx.f">   *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zgesvdx.f">
 *> [ZIP]</a>   *> [ZIP]</a>
 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zgesvdx.f">   *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zgesvdx.f">
 *> [TXT]</a>  *> [TXT]</a>
 *> \endhtmlonly   *> \endhtmlonly
 *  *
 *  Definition:  *  Definition:
 *  ===========  *  ===========
 *  *
 *     SUBROUTINE CGESVDX( JOBU, JOBVT, RANGE, M, N, A, LDA, VL, VU,   *     SUBROUTINE CGESVDX( 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, RWORK, IWORK, INFO )  *    $                    LWORK, RWORK, IWORK, INFO )
 *        *
 *  *
 *     .. Scalar Arguments ..  *     .. Scalar Arguments ..
 *      CHARACTER          JOBU, JOBVT, RANGE  *      CHARACTER          JOBU, JOBVT, RANGE
Line 31 Line 31
 *     .. Array Arguments ..  *     .. Array Arguments ..
 *      INTEGER            IWORK( * )  *      INTEGER            IWORK( * )
 *      DOUBLE PRECISION   S( * ), RWORK( * )  *      DOUBLE PRECISION   S( * ), RWORK( * )
 *      COMPLEX*16         A( LDA, * ), U( LDU, * ), VT( LDVT, * ),   *      COMPLEX*16         A( LDA, * ), U( LDU, * ), VT( LDVT, * ),
 *     $                   WORK( * )  *     $                   WORK( * )
 *     ..  *     ..
 *  *
Line 60 Line 60
 *>  *>
 *>  Note that the routine returns V**T, not V.  *>  Note that the routine returns V**T, not V.
 *> \endverbatim  *> \endverbatim
 *     *
 *  Arguments:  *  Arguments:
 *  ==========  *  ==========
 *  *
Line 69 Line 69
 *>          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 81 Line 81
 *>           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 93 Line 93
 *>          = '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 158 Line 158
 *> \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 172 Line 172
 *> \param[out] U  *> \param[out] U
 *> \verbatim  *> \verbatim
 *>          U is COMPLEX*16 array, dimension (LDU,UCOL)  *>          U is COMPLEX*16 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 190 Line 190
 *> \param[out] VT  *> \param[out] VT
 *> \verbatim  *> \verbatim
 *>          VT is COMPLEX*16 array, dimension (LDVT,N)  *>          VT is COMPLEX*16 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 215 Line 215
 *> \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 237 Line 237
 *> \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 256 Line 256
 *  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  *> \date June 2016
 *  *
 *> \ingroup complex16GEsing  *> \ingroup complex16GEsing
 *  *
 *  =====================================================================  *  =====================================================================
       SUBROUTINE ZGESVDX( JOBU, JOBVT, RANGE, M, N, A, LDA, VL, VU,         SUBROUTINE ZGESVDX( 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, RWORK, IWORK, INFO )       $                    LWORK, RWORK, IWORK, INFO )
 *  *
 *  -- LAPACK driver routine (version 3.6.1) --  *  -- LAPACK driver routine (version 3.7.0) --
 *  -- 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  *     June 2016
Line 283 Line 283
 *     .. Array Arguments ..  *     .. Array Arguments ..
       INTEGER            IWORK( * )        INTEGER            IWORK( * )
       DOUBLE PRECISION   S( * ), RWORK( * )        DOUBLE PRECISION   S( * ), RWORK( * )
       COMPLEX*16         A( LDA, * ), U( LDU, * ), VT( LDVT, * ),         COMPLEX*16         A( LDA, * ), U( LDU, * ), VT( LDVT, * ),
      $                   WORK( * )       $                   WORK( * )
 *     ..  *     ..
 *  *
Line 300 Line 300
       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, ITEMPR, ITGKZ,        $                   ITAU, ITAUP, ITAUQ, ITEMP, ITEMPR, ITGKZ,
      $                   IUTGK, J, K, MAXWRK, MINMN, MINWRK, MNTHR       $                   IUTGK, J, K, MAXWRK, MINMN, MINWRK, MNTHR
       DOUBLE PRECISION   ABSTOL, ANRM, BIGNUM, EPS, SMLNUM        DOUBLE PRECISION   ABSTOL, ANRM, BIGNUM, EPS, SMLNUM
 *     ..  *     ..
Line 480 Line 480
          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 524 Line 524
             ITEMP = ITAU + N              ITEMP = ITAU + N
             CALL ZGEQRF( M, N, A, LDA, WORK( ITAU ), WORK( ITEMP ),              CALL ZGEQRF( 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+3*N, prefer N*N+N+2*N*NB)  *           (Workspace: need N*N+3*N, prefer N*N+N+2*N*NB)
 *     *
             IQRF = ITEMP              IQRF = ITEMP
             ITAUQ = ITEMP + N*N              ITAUQ = ITEMP + N*N
             ITAUP = ITAUQ + N              ITAUP = ITAUQ + N
             ITEMP = ITAUP + N              ITEMP = ITAUP + N
             ID = 1                 ID = 1
             IE = ID + N              IE = ID + N
             ITGKZ = IE + N              ITGKZ = IE + N
             CALL ZLACPY( 'U', N, N, A, LDA, WORK( IQRF ), N )              CALL ZLACPY( 'U', N, N, A, LDA, WORK( IQRF ), N )
             CALL ZLASET( 'L', N-1, N-1, CZERO, CZERO,              CALL ZLASET( 'L', N-1, N-1, CZERO, CZERO,
      $                   WORK( IQRF+1 ), N )       $                   WORK( IQRF+1 ), N )
             CALL ZGEBRD( N, N, WORK( IQRF ), N, RWORK( ID ),               CALL ZGEBRD( N, N, WORK( IQRF ), N, RWORK( ID ),
      $                   RWORK( IE ), WORK( ITAUQ ), WORK( ITAUP ),       $                   RWORK( IE ), WORK( ITAUQ ), WORK( ITAUP ),
      $                   WORK( ITEMP ), LWORK-ITEMP+1, INFO )       $                   WORK( ITEMP ), LWORK-ITEMP+1, INFO )
             ITEMPR = ITGKZ + N*(N*2+1)              ITEMPR = ITGKZ + N*(N*2+1)
 *  *
 *           Solve eigenvalue problem TGK*Z=Z*S.  *           Solve eigenvalue problem TGK*Z=Z*S.
 *           (Workspace: need 2*N*N+14*N)            *           (Workspace: need 2*N*N+14*N)
 *                     *
             CALL DBDSVDX( 'U', JOBZ, RNGTGK, N, RWORK( ID ),              CALL DBDSVDX( 'U', JOBZ, RNGTGK, N, RWORK( ID ),
      $                    RWORK( IE ), VL, VU, ILTGK, IUTGK, NS, S,       $                    RWORK( IE ), VL, VU, ILTGK, IUTGK, NS, S,
      $                    RWORK( ITGKZ ), N*2, RWORK( ITEMPR ),       $                    RWORK( ITGKZ ), N*2, RWORK( ITEMPR ),
Line 567 Line 567
 *              Call ZUNMBR to compute QB*UB.  *              Call ZUNMBR to compute QB*UB.
 *              (Workspace in WORK( ITEMP ): need N, prefer N*NB)  *              (Workspace in WORK( ITEMP ): need N, prefer N*NB)
 *  *
                CALL ZUNMBR( 'Q', 'L', 'N', N, NS, N, WORK( IQRF ), N,                  CALL ZUNMBR( '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 ZUNMQR to compute Q*(QB*UB).  *              Call ZUNMQR to compute Q*(QB*UB).
 *              (Workspace in WORK( ITEMP ): need N, prefer N*NB)  *              (Workspace in WORK( ITEMP ): need N, prefer N*NB)
 *  *
                CALL ZUNMQR( 'L', 'N', M, NS, N, A, LDA,                  CALL ZUNMQR( '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 594 Line 594
 *              Call ZUNMBR to compute VB**T * PB**T  *              Call ZUNMBR 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 ZUNMBR( 'P', 'R', 'C', NS, N, N, WORK( IQRF ), N,                  CALL ZUNMBR( 'P', 'R', 'C', 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 610 Line 610
 *  *
             ITAUQ = 1              ITAUQ = 1
             ITAUP = ITAUQ + N              ITAUP = ITAUQ + N
             ITEMP = ITAUP + N                        ITEMP = ITAUP + N
             ID = 1              ID = 1
             IE = ID + N              IE = ID + N
             ITGKZ = IE + N              ITGKZ = IE + N
             CALL ZGEBRD( M, N, A, LDA, RWORK( ID ), RWORK( IE ),               CALL ZGEBRD( M, N, A, LDA, RWORK( ID ), RWORK( IE ),
      $                   WORK( ITAUQ ), WORK( ITAUP ), WORK( ITEMP ),       $                   WORK( ITAUQ ), WORK( ITAUP ), WORK( ITEMP ),
      $                   LWORK-ITEMP+1, INFO )       $                   LWORK-ITEMP+1, INFO )
             ITEMPR = ITGKZ + N*(N*2+1)              ITEMPR = ITGKZ + N*(N*2+1)
 *  *
 *           Solve eigenvalue problem TGK*Z=Z*S.  *           Solve eigenvalue problem TGK*Z=Z*S.
 *           (Workspace: need 2*N*N+14*N)            *           (Workspace: need 2*N*N+14*N)
 *                       *
             CALL DBDSVDX( 'U', JOBZ, RNGTGK, N, RWORK( ID ),              CALL DBDSVDX( 'U', JOBZ, RNGTGK, N, RWORK( ID ),
      $                    RWORK( IE ), VL, VU, ILTGK, IUTGK, NS, S,        $                    RWORK( IE ), VL, VU, ILTGK, IUTGK, NS, S,
      $                    RWORK( ITGKZ ), N*2, RWORK( ITEMPR ),        $                    RWORK( ITGKZ ), N*2, RWORK( ITEMPR ),
      $                    IWORK, INFO)       $                    IWORK, INFO)
 *  *
 *           If needed, compute left singular vectors.  *           If needed, compute left singular vectors.
Line 632 Line 632
             IF( WANTU ) THEN              IF( WANTU ) THEN
                K = ITGKZ                 K = ITGKZ
                DO I = 1, NS                 DO I = 1, NS
                   DO J = 1, N                                  DO J = 1, N
                      U( J, I ) = DCMPLX( RWORK( K ), ZERO )                       U( J, I ) = DCMPLX( RWORK( K ), ZERO )
                      K = K + 1                       K = K + 1
                   END DO                    END DO
Line 642 Line 642
 *  *
 *              Call ZUNMBR to compute QB*UB.  *              Call ZUNMBR to compute QB*UB.
 *              (Workspace in WORK( ITEMP ): need N, prefer N*NB)  *              (Workspace in WORK( ITEMP ): need N, prefer N*NB)
 *     *
                CALL ZUNMBR( 'Q', 'L', 'N', M, NS, N, A, LDA,                  CALL ZUNMBR( '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 663 Line 663
 *              Call ZUNMBR to compute VB**T * PB**T  *              Call ZUNMBR 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 ZUNMBR( 'P', 'R', 'C', NS, N, N, A, LDA,                  CALL ZUNMBR( 'P', 'R', 'C', 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 676 Line 676
          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 691 Line 691
 *           Copy L into WORK and bidiagonalize it:  *           Copy L into WORK and bidiagonalize it:
 *           (Workspace in WORK( ITEMP ): need M*M+3*M, prefer M*M+M+2*M*NB)  *           (Workspace in WORK( ITEMP ): need M*M+3*M, prefer M*M+M+2*M*NB)
 *  *
             ILQF = ITEMP                     ILQF = ITEMP
             ITAUQ = ILQF + M*M              ITAUQ = ILQF + M*M
             ITAUP = ITAUQ + M              ITAUP = ITAUQ + M
             ITEMP = ITAUP + M              ITEMP = ITAUP + M
Line 699 Line 699
             IE = ID + M              IE = ID + M
             ITGKZ = IE + M              ITGKZ = IE + M
             CALL ZLACPY( 'L', M, M, A, LDA, WORK( ILQF ), M )              CALL ZLACPY( 'L', M, M, A, LDA, WORK( ILQF ), M )
             CALL ZLASET( 'U', M-1, M-1, CZERO, CZERO,               CALL ZLASET( 'U', M-1, M-1, CZERO, CZERO,
      $                   WORK( ILQF+M ), M )       $                   WORK( ILQF+M ), M )
             CALL ZGEBRD( M, M, WORK( ILQF ), M, RWORK( ID ),              CALL ZGEBRD( M, M, WORK( ILQF ), M, RWORK( ID ),
      $                   RWORK( IE ), WORK( ITAUQ ), WORK( ITAUP ),        $                   RWORK( IE ), WORK( ITAUQ ), WORK( ITAUP ),
      $                   WORK( ITEMP ), LWORK-ITEMP+1, INFO )       $                   WORK( ITEMP ), LWORK-ITEMP+1, INFO )
             ITEMPR = ITGKZ + M*(M*2+1)              ITEMPR = ITGKZ + M*(M*2+1)
 *  *
 *           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)
 *  *
             CALL DBDSVDX( 'U', JOBZ, RNGTGK, M, RWORK( ID ),              CALL DBDSVDX( 'U', JOBZ, RNGTGK, M, RWORK( ID ),
      $                    RWORK( IE ), VL, VU, ILTGK, IUTGK, NS, S,        $                    RWORK( IE ), VL, VU, ILTGK, IUTGK, NS, S,
      $                    RWORK( ITGKZ ), M*2, RWORK( ITEMPR ),        $                    RWORK( ITGKZ ), M*2, RWORK( ITEMPR ),
      $                    IWORK, INFO)       $                    IWORK, INFO)
 *  *
 *           If needed, compute left singular vectors.  *           If needed, compute left singular vectors.
Line 729 Line 729
 *              Call ZUNMBR to compute QB*UB.  *              Call ZUNMBR to compute QB*UB.
 *              (Workspace in WORK( ITEMP ): need M, prefer M*NB)  *              (Workspace in WORK( ITEMP ): need M, prefer M*NB)
 *  *
                CALL ZUNMBR( 'Q', 'L', 'N', M, NS, M, WORK( ILQF ), M,                  CALL ZUNMBR( '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 751 Line 751
 *              Call ZUNMBR to compute (VB**T)*(PB**T)  *              Call ZUNMBR 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 ZUNMBR( 'P', 'R', 'C', NS, M, M, WORK( ILQF ), M,                  CALL ZUNMBR( 'P', 'R', 'C', 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 ZUNMLQ to compute ((VB**T)*(PB**T))*Q.  *              Call ZUNMLQ 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 ZUNMLQ( 'R', 'N', NS, N, M, A, LDA,                  CALL ZUNMLQ( '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 2*M+N, prefer 2*M+(M+N)*NB)  *           (Workspace: need 2*M+N, prefer 2*M+(M+N)*NB)
 *              *
             ITAUQ = 1              ITAUQ = 1
             ITAUP = ITAUQ + M              ITAUP = ITAUQ + M
             ITEMP = ITAUP + M              ITEMP = ITAUP + M
             ID = 1              ID = 1
             IE = ID + M              IE = ID + M
             ITGKZ = IE + M              ITGKZ = IE + M
             CALL ZGEBRD( M, N, A, LDA, RWORK( ID ), RWORK( IE ),               CALL ZGEBRD( M, N, A, LDA, RWORK( ID ), RWORK( IE ),
      $                   WORK( ITAUQ ), WORK( ITAUP ), WORK( ITEMP ),       $                   WORK( ITAUQ ), WORK( ITAUP ), WORK( ITEMP ),
      $                   LWORK-ITEMP+1, INFO )       $                   LWORK-ITEMP+1, INFO )
             ITEMPR = ITGKZ + M*(M*2+1)              ITEMPR = ITGKZ + M*(M*2+1)
 *  *
 *           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)
 *            *
             CALL DBDSVDX( 'L', JOBZ, RNGTGK, M, RWORK( ID ),               CALL DBDSVDX( 'L', JOBZ, RNGTGK, M, RWORK( ID ),
      $                    RWORK( IE ), VL, VU, ILTGK, IUTGK, NS, S,        $                    RWORK( IE ), VL, VU, ILTGK, IUTGK, NS, S,
      $                    RWORK( ITGKZ ), M*2, RWORK( ITEMPR ),        $                    RWORK( ITGKZ ), M*2, RWORK( ITEMPR ),
      $                    IWORK, INFO)       $                    IWORK, INFO)
 *   *
 *           If needed, compute left singular vectors.  *           If needed, compute left singular vectors.
 *  *
             IF( WANTU ) THEN              IF( WANTU ) THEN
Line 806 Line 806
 *              Call ZUNMBR to compute QB*UB.  *              Call ZUNMBR to compute QB*UB.
 *              (Workspace in WORK( ITEMP ): need M, prefer M*NB)  *              (Workspace in WORK( ITEMP ): need M, prefer M*NB)
 *  *
                CALL ZUNMBR( 'Q', 'L', 'N', M, NS, N, A, LDA,                  CALL ZUNMBR( '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 828 Line 828
 *              Call ZUNMBR to compute VB**T * PB**T  *              Call ZUNMBR 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 ZUNMBR( 'P', 'R', 'C', NS, N, M, A, LDA,                  CALL ZUNMBR( 'P', 'R', 'C', 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.4

CVSweb interface <joel.bertrand@systella.fr>