[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.1 and 1.9

version 1.1, 2015/11/26 11:44:21 version 1.9, 2023/08/07 08:39:19
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 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 )
 *        *
 *  *
 *     .. 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( * )
 *     ..  *     ..
 *  *
 *  *
 *  Purpose  *> \par Purpose:
 *  =======  *  =============
   *>
   *> \verbatim
   *>
   *>  ZGESVDX computes the singular value decomposition (SVD) of a complex
   *>  M-by-N matrix A, optionally computing the left and/or right singular
   *>  vectors. The SVD is written
   *>
   *>      A = U * SIGMA * transpose(V)
   *>
   *>  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 unitary matrix, and
   *>  V is an N-by-N unitary matrix.  The diagonal elements of SIGMA
   *>  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
   *>  U and V are the left and right singular vectors of A.
   *>
   *>  ZGESVDX uses an eigenvalue problem for obtaining the SVD, which
   *>  allows for the computation of a subset of singular values and
   *>  vectors. See DBDSVDX for details.
   *>
   *>  Note that the routine returns V**T, not V.
   *> \endverbatim
 *  *
 *  ZGESVDX computes the singular value decomposition (SVD) of a complex  
 *  M-by-N matrix A, optionally computing the left and/or right singular  
 *  vectors. The SVD is written  
 *   
 *       A = U * SIGMA * transpose(V)  
 *   
 *  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 unitary matrix, and  
 *  V is an N-by-N unitary matrix.  The diagonal elements of SIGMA  
 *  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  
 *  U and V are the left and right singular vectors of A.  
 *   
 *  ZGESVDX uses an eigenvalue problem for obtaining the SVD, which   
 *  allows for the computation of a subset of singular values and   
 *  vectors. See DBDSVDX for details.  
 *   
 *  Note that the routine returns V**T, not V.  
 *     
 *  Arguments:  *  Arguments:
 *  ==========  *  ==========
 *  *
Line 66 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 78 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 90 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 107 Line 110
 *>  *>
 *> \param[in,out] A  *> \param[in,out] A
 *> \verbatim  *> \verbatim
 *>          A is COMPLEX array, dimension (LDA,N)  *>          A is COMPLEX*16 array, dimension (LDA,N)
 *>          On entry, the M-by-N matrix A.  *>          On entry, the M-by-N matrix A.
 *>          On exit, the contents of A are destroyed.  *>          On exit, the contents of A are destroyed.
 *> \endverbatim  *> \endverbatim
Line 121 Line 124
 *> \param[in] VL  *> \param[in] VL
 *> \verbatim  *> \verbatim
 *>          VL is DOUBLE PRECISION  *>          VL is DOUBLE PRECISION
 *>          VL >=0.  *>          If RANGE='V', the lower bound of the interval to
   *>          be searched for singular values. VU > VL.
   *>          Not referenced if RANGE = 'A' or 'I'.
 *> \endverbatim  *> \endverbatim
 *>  *>
 *> \param[in] VU  *> \param[in] VU
 *> \verbatim  *> \verbatim
 *>          VU is DOUBLE PRECISION  *>          VU is DOUBLE PRECISION
 *>          If RANGE='V', the lower and upper bounds of the interval to  *>          If RANGE='V', the upper bound of the interval to
 *>          be searched for singular values. VU > VL.  *>          be searched for singular values. VU > VL.
 *>          Not referenced if RANGE = 'A' or 'I'.  *>          Not referenced if RANGE = 'A' or 'I'.
 *> \endverbatim  *> \endverbatim
Line 135 Line 140
 *> \param[in] IL  *> \param[in] IL
 *> \verbatim  *> \verbatim
 *>          IL is INTEGER  *>          IL is INTEGER
   *>          If RANGE='I', the index of the
   *>          smallest singular value to be returned.
   *>          1 <= IL <= IU <= min(M,N), if min(M,N) > 0.
   *>          Not referenced if RANGE = 'A' or 'V'.
 *> \endverbatim  *> \endverbatim
 *>  *>
 *> \param[in] IU  *> \param[in] IU
 *> \verbatim  *> \verbatim
 *>          IU is INTEGER  *>          IU is INTEGER
 *>          If RANGE='I', the indices (in ascending order) of the  *>          If RANGE='I', the index of the
 *>          smallest and largest singular values to be returned.  *>          largest singular value to be returned.
 *>          1 <= IL <= IU <= min(M,N), if min(M,N) > 0.  *>          1 <= IL <= IU <= min(M,N), if min(M,N) > 0.
 *>          Not referenced if RANGE = 'A' or 'V'.  *>          Not referenced if RANGE = 'A' or 'V'.
 *> \endverbatim  *> \endverbatim
Line 149 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 163 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 ILQFin 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 181 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 206 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 228 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 247 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 November 2015  
 *  *
 *> \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.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 2015  
 *  *
 *     .. Scalar Arguments ..  *     .. Scalar Arguments ..
       CHARACTER          JOBU, JOBVT, RANGE        CHARACTER          JOBU, JOBVT, RANGE
Line 274 Line 280
 *     .. 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 291 Line 297
       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, ITEMPR, ITGKZ,
      $                   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
 *     ..  *     ..
 *     .. Local Arrays ..  *     .. Local Arrays ..
       DOUBLE PRECISION   DUM( 1 )        DOUBLE PRECISION   DUM( 1 )
 *     ..  *     ..
 *     .. External Subroutines ..  *     .. External Subroutines ..
       EXTERNAL           ZGEBRD, ZGELQF, ZGEQRF, ZLASCL, ZLASET,        EXTERNAL           ZGEBRD, ZGELQF, ZGEQRF, ZLASCL, ZLASET, ZLACPY,
      $                   DLASCL, XERBLA       $                   ZUNMLQ, ZUNMBR, ZUNMQR, DBDSVDX, DLASCL, XERBLA
 *     ..  *     ..
 *     .. External Functions ..  *     .. External Functions ..
       LOGICAL            LSAME        LOGICAL            LSAME
Line 364 Line 370
          IF( INFO.EQ.0 ) THEN           IF( INFO.EQ.0 ) THEN
             IF( WANTU .AND. LDU.LT.M ) THEN              IF( WANTU .AND. LDU.LT.M ) THEN
                INFO = -15                 INFO = -15
             ELSE IF( WANTVT .AND. LDVT.LT.MINMN ) THEN              ELSE IF( WANTVT ) THEN
                INFO = -16                 IF( INDS ) THEN
                      IF( LDVT.LT.IU-IL+1 ) THEN
                          INFO = -17
                      END IF
                  ELSE IF( LDVT.LT.MINMN ) THEN
                      INFO = -17
                  END IF
             END IF              END IF
          END IF           END IF
       END IF        END IF
Line 387 Line 399
 *  *
 *                 Path 1 (M much larger than N)  *                 Path 1 (M much larger than N)
 *  *
                   MAXWRK = N + N*                    MINWRK = N*(N+5)
      $                     ILAENV( 1, 'DGEQRF', ' ', M, N, -1, -1 )                    MAXWRK = N + N*ILAENV(1,'ZGEQRF',' ',M,N,-1,-1)
                   MAXWRK = MAX( MAXWRK, N*N + N + 2*N*                    MAXWRK = MAX(MAXWRK,
      $                     ILAENV( 1, 'DGEBRD', ' ', N, N, -1, -1 ) )       $                     N*N+2*N+2*N*ILAENV(1,'ZGEBRD',' ',N,N,-1,-1))
                   MINWRK = N*(N+4)                    IF (WANTU .OR. WANTVT) THEN
                        MAXWRK = MAX(MAXWRK,
        $                       N*N+2*N+N*ILAENV(1,'ZUNMQR','LN',N,N,N,-1))
                     END IF
                ELSE                 ELSE
 *  *
 *                 Path 2 (M at least N, but not much larger)  *                 Path 2 (M at least N, but not much larger)
 *  *
                   MAXWRK = 2*N + ( M+N )*                    MINWRK = 3*N + M
      $                     ILAENV( 1, 'ZGEBRD', ' ', M, N, -1, -1 )                    MAXWRK = 2*N + (M+N)*ILAENV(1,'ZGEBRD',' ',M,N,-1,-1)
                   MINWRK = 2*N + M                    IF (WANTU .OR. WANTVT) THEN
                        MAXWRK = MAX(MAXWRK,
        $                        2*N+N*ILAENV(1,'ZUNMQR','LN',N,N,N,-1))
                     END IF
                END IF                 END IF
             ELSE              ELSE
                MNTHR = ILAENV( 6, 'ZGESVD', JOBU // JOBVT, M, N, 0, 0 )                 MNTHR = ILAENV( 6, 'ZGESVD', JOBU // JOBVT, M, N, 0, 0 )
Line 406 Line 424
 *  *
 *                 Path 1t (N much larger than M)  *                 Path 1t (N much larger than M)
 *  *
                   MAXWRK = M + M*                    MINWRK = M*(M+5)
      $                     ILAENV( 1, 'ZGELQF', ' ', M, N, -1, -1 )                    MAXWRK = M + M*ILAENV(1,'ZGELQF',' ',M,N,-1,-1)
                   MAXWRK = MAX( MAXWRK, M*M + M + 2*M*                    MAXWRK = MAX(MAXWRK,
      $                     ILAENV( 1, 'ZGEBRD', ' ', M, M, -1, -1 ) )       $                     M*M+2*M+2*M*ILAENV(1,'ZGEBRD',' ',M,M,-1,-1))
                   MINWRK = M*(M+4)                    IF (WANTU .OR. WANTVT) THEN
                        MAXWRK = MAX(MAXWRK,
        $                       M*M+2*M+M*ILAENV(1,'ZUNMQR','LN',M,M,M,-1))
                     END IF
                ELSE                 ELSE
 *  *
 *                 Path 2t (N greater than M, but not much larger)  *                 Path 2t (N greater than M, but not much larger)
 *  *
                   MAXWRK = M*(M*2+19) + ( M+N )*  *
      $                     ILAENV( 1, 'ZGEBRD', ' ', M, N, -1, -1 )                    MINWRK = 3*M + N
                   MINWRK = 2*M + N                    MAXWRK = 2*M + (M+N)*ILAENV(1,'ZGEBRD',' ',M,N,-1,-1)
                     IF (WANTU .OR. WANTVT) THEN
                        MAXWRK = MAX(MAXWRK,
        $                        2*M+M*ILAENV(1,'ZUNMQR','LN',M,M,M,-1))
                     END IF
                END IF                 END IF
             END IF              END IF
          END IF           END IF
Line 444 Line 469
 *  *
 *     Set singular values indices accord to RANGE='A'.  *     Set singular values indices accord to RANGE='A'.
 *  *
       ALLS = LSAME( RANGE, 'A' )  
       INDS = LSAME( RANGE, 'I' )  
       IF( ALLS ) THEN        IF( ALLS ) THEN
          RNGTGK = 'I'           RNGTGK = 'I'
          ILTGK = 1           ILTGK = 1
Line 454 Line 477
          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 498 Line 521
             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 )
             ITEMP = 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( ITEMP ),       $                    RWORK( ITGKZ ), N*2, RWORK( ITEMPR ),
      $                    IWORK, INFO)       $                    IWORK, INFO)
 *  *
 *           If needed, compute left singular vectors.  *           If needed, compute left singular vectors.
Line 536 Line 559
                   END DO                    END DO
                   K = K + N                    K = K + N
                END DO                 END DO
                CALL ZLASET( 'A', M-N, N, CZERO, CZERO, U( N+1,1 ), LDU )                 CALL ZLASET( 'A', M-N, NS, CZERO, CZERO, U( N+1,1 ), LDU)
 *  *
 *              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 568 Line 591
 *              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 584 Line 607
 *  *
             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 )
             ITEMP = 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( ITEMP ),        $                    RWORK( ITGKZ ), N*2, RWORK( ITEMPR ),
      $                    IWORK, INFO)       $                    IWORK, INFO)
 *  *
 *           If needed, compute left singular vectors.  *           If needed, compute left singular vectors.
Line 606 Line 629
             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
                   K = K + N                    K = K + N
                END DO                 END DO
                CALL ZLASET( 'A', M-N, N, CZERO, CZERO, U( N+1,1 ), LDU )                 CALL ZLASET( 'A', M-N, NS, CZERO, CZERO, U( N+1,1 ), LDU)
 *  *
 *              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 637 Line 660
 *              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 650 Line 673
          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 665 Line 688
 *           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 673 Line 696
             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 )
             ITEMP = 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( ITEMP ),        $                    RWORK( ITGKZ ), M*2, RWORK( ITEMPR ),
      $                    IWORK, INFO)       $                    IWORK, INFO)
 *  *
 *           If needed, compute left singular vectors.  *           If needed, compute left singular vectors.
Line 703 Line 726
 *              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 719 Line 742
                   END DO                    END DO
                   K = K + M                    K = K + M
                END DO                 END DO
                CALL ZLASET( 'A', M, N-M, CZERO, CZERO,                  CALL ZLASET( 'A', NS, N-M, CZERO, CZERO,
      $                      VT( 1,M+1 ), LDVT )       $                      VT( 1,M+1 ), LDVT )
 *  *
 *              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 )
             ITEMP = 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( ITEMP ),        $                    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 780 Line 803
 *              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 796 Line 819
                   END DO                    END DO
                   K = K + M                    K = K + M
                END DO                 END DO
                CALL ZLASET( 'A', M, N-M, CZERO, CZERO,                  CALL ZLASET( 'A', NS, N-M, CZERO, CZERO,
      $                      VT( 1,M+1 ), LDVT )       $                      VT( 1,M+1 ), LDVT )
 *  *
 *              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.1  
changed lines
  Added in v.1.9

CVSweb interface <joel.bertrand@systella.fr>