[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.2

version 1.1, 2015/11/26 11:44:21 version 1.2, 2016/08/27 15:27:12
Line 36 Line 36
 *     ..  *     ..
 *  *
 *  *
 *  Purpose  *> \par Purpose:
 *  =======  *  =============
 *  *>
 *  ZGESVDX computes the singular value decomposition (SVD) of a complex  *> \verbatim
 *  M-by-N matrix A, optionally computing the left and/or right singular  *>
 *  vectors. The SVD is written  *>  ZGESVDX computes the singular value decomposition (SVD) of a complex
 *   *>  M-by-N matrix A, optionally computing the left and/or right singular
 *       A = U * SIGMA * transpose(V)  *>  vectors. The SVD is written
 *   *>
 *  where SIGMA is an M-by-N matrix which is zero except for its  *>      A = U * SIGMA * transpose(V)
 *  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  *>  where SIGMA is an M-by-N matrix which is zero except for its
 *  are the singular values of A; they are real and non-negative, and  *>  min(m,n) diagonal elements, U is an M-by-M unitary matrix, and
 *  are returned in descending order.  The first min(m,n) columns of  *>  V is an N-by-N unitary matrix.  The diagonal elements of SIGMA
 *  U and V are the left and right singular vectors of A.  *>  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
 *  ZGESVDX uses an eigenvalue problem for obtaining the SVD, which   *>  U and V are the left and right singular vectors of A.
 *  allows for the computation of a subset of singular values and   *>
 *  vectors. See DBDSVDX for details.  *>  ZGESVDX uses an eigenvalue problem for obtaining the SVD, which
 *   *>  allows for the computation of a subset of singular values and
 *  Note that the routine returns V**T, not V.  *>  vectors. See DBDSVDX for details.
   *>
   *>  Note that the routine returns V**T, not V.
   *> \endverbatim
 *     *   
 *  Arguments:  *  Arguments:
 *  ==========  *  ==========
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 167 Line 176
 *>          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 252 Line 261
 *> \author Univ. of Colorado Denver   *> \author Univ. of Colorado Denver 
 *> \author NAG Ltd.   *> \author NAG Ltd. 
 *  *
 *> \date November 2015  *> \date June 2016
 *  *
 *> \ingroup complex16GEsing  *> \ingroup complex16GEsing
 *  *
Line 261 Line 270
      $                    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 (version 3.6.1) --
 *  -- 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  *     June 2016
 *  *
 *     .. Scalar Arguments ..  *     .. Scalar Arguments ..
       CHARACTER          JOBU, JOBVT, RANGE        CHARACTER          JOBU, JOBVT, RANGE
Line 291 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, 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 ..
Line 364 Line 373
          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 402
 *  *
 *                 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 427
 *  *
 *                 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 472
 *  *
 *     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 515 Line 541
             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 562
                   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)
Line 591 Line 617
             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 612 Line 638
                   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)
Line 678 Line 704
             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 719 Line 745
                   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)
Line 755 Line 781
             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.
Line 796 Line 822
                   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
Removed from v.1.1  
changed lines
  Added in v.1.2

CVSweb interface <joel.bertrand@systella.fr>