floraachy
/
OpenBLAS
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

BUGFIX: removed fixes for bugs #148 and #149, because info for xerbla is wrong

This commit is contained in:
Werner Saar 2016-03-07 10:34:04 +01:00
parent 26b3b3a3e6
commit 10c2ebdfc5
4 changed files with 124 additions and 241 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

85
lapack-netlib/SRC/cgesvdx.f
View File

@ -170,7 +170,7 @@
*> 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 ILQFin advance and an upper
*> bound must be used. *> bound must be used.
*> \endverbatim *> \endverbatim
*> *>
@ -294,8 +294,8 @@
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, ITGKZ, IUTGK,
$ IUTGK, J, K, MAXWRK, MINMN, MINWRK, MNTHR $ J, K, MAXWRK, MINMN, MINWRK, MNTHR
REAL ABSTOL, ANRM, BIGNUM, EPS, SMLNUM REAL ABSTOL, ANRM, BIGNUM, EPS, SMLNUM
* .. * ..
* .. Local Arrays .. * .. Local Arrays ..
@ -367,14 +367,8 @@
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 ) THEN ELSE IF( WANTVT .AND. LDVT.LT.MINMN ) THEN
IF( INDS ) THEN INFO = -16
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
@ -396,24 +390,18 @@
* *
* Path 1 (M much larger than N) * Path 1 (M much larger than N)
* *
MINWRK = N*(N+5) MAXWRK = N + N*
MAXWRK = N + N*ILAENV(1,'CGEQRF',' ',M,N,-1,-1) $ ILAENV( 1, 'SGEQRF', ' ', M, N, -1, -1 )
MAXWRK = MAX(MAXWRK, MAXWRK = MAX( MAXWRK, N*N + N + 2*N*
$ N*N+2*N+2*N*ILAENV(1,'CGEBRD',' ',N,N,-1,-1)) $ ILAENV( 1, 'SGEBRD', ' ', N, N, -1, -1 ) )
IF (WANTU .OR. WANTVT) THEN MINWRK = N*(N+4)
MAXWRK = MAX(MAXWRK,
$ N*N+2*N+N*ILAENV(1,'CUNMQR','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)
* *
MINWRK = 3*N + M MAXWRK = 2*N + ( M+N )*
MAXWRK = 2*N + (M+N)*ILAENV(1,'CGEBRD',' ',M,N,-1,-1) $ ILAENV( 1, 'CGEBRD', ' ', M, N, -1, -1 )
IF (WANTU .OR. WANTVT) THEN MINWRK = 2*N + M
MAXWRK = MAX(MAXWRK,
$ 2*N+N*ILAENV(1,'CUNMQR','LN',N,N,N,-1))
END IF
END IF END IF
ELSE ELSE
MNTHR = ILAENV( 6, 'CGESVD', JOBU // JOBVT, M, N, 0, 0 ) MNTHR = ILAENV( 6, 'CGESVD', JOBU // JOBVT, M, N, 0, 0 )
@ -421,25 +409,18 @@
* *
* Path 1t (N much larger than M) * Path 1t (N much larger than M)
* *
MINWRK = M*(M+5) MAXWRK = M + M*
MAXWRK = M + M*ILAENV(1,'CGELQF',' ',M,N,-1,-1) $ ILAENV( 1, 'CGELQF', ' ', M, N, -1, -1 )
MAXWRK = MAX(MAXWRK, MAXWRK = MAX( MAXWRK, M*M + M + 2*M*
$ M*M+2*M+2*M*ILAENV(1,'CGEBRD',' ',M,M,-1,-1)) $ ILAENV( 1, 'CGEBRD', ' ', M, M, -1, -1 ) )
IF (WANTU .OR. WANTVT) THEN MINWRK = M*(M+4)
MAXWRK = MAX(MAXWRK,
$ M*M+2*M+M*ILAENV(1,'CUNMQR','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 )*
MINWRK = 3*M + N $ ILAENV( 1, 'CGEBRD', ' ', M, N, -1, -1 )
MAXWRK = 2*M + (M+N)*ILAENV(1,'CGEBRD',' ',M,N,-1,-1) MINWRK = 2*M + N
IF (WANTU .OR. WANTVT) THEN
MAXWRK = MAX(MAXWRK,
$ 2*M+M*ILAENV(1,'CUNMQR','LN',M,M,M,-1))
END IF
END IF END IF
END IF END IF
END IF END IF
@ -537,14 +518,14 @@
CALL CGEBRD( N, N, WORK( IQRF ), N, RWORK( ID ), CALL CGEBRD( 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) ITEMP = 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 SBDSVDX( 'U', JOBZ, RNGTGK, N, RWORK( ID ), CALL SBDSVDX( '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( ITEMP ),
$ IWORK, INFO) $ IWORK, INFO)
* *
* If needed, compute left singular vectors. * If needed, compute left singular vectors.
@ -558,7 +539,7 @@
END DO END DO
K = K + N K = K + N
END DO END DO
CALL CLASET( 'A', M-N, NS, CZERO, CZERO, U( N+1,1 ), LDU) CALL CLASET( 'A', M-N, N, CZERO, CZERO, U( N+1,1 ), LDU )
* *
* Call CUNMBR to compute QB*UB. * Call CUNMBR to compute QB*UB.
* (Workspace in WORK( ITEMP ): need N, prefer N*NB) * (Workspace in WORK( ITEMP ): need N, prefer N*NB)
@ -613,14 +594,14 @@
CALL CGEBRD( M, N, A, LDA, RWORK( ID ), RWORK( IE ), CALL CGEBRD( 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) ITEMP = 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 SBDSVDX( 'U', JOBZ, RNGTGK, N, RWORK( ID ), CALL SBDSVDX( '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( ITEMP ),
$ IWORK, INFO) $ IWORK, INFO)
* *
* If needed, compute left singular vectors. * If needed, compute left singular vectors.
@ -634,7 +615,7 @@
END DO END DO
K = K + N K = K + N
END DO END DO
CALL CLASET( 'A', M-N, NS, CZERO, CZERO, U( N+1,1 ), LDU) CALL CLASET( 'A', M-N, N, CZERO, CZERO, U( N+1,1 ), LDU )
* *
* Call CUNMBR to compute QB*UB. * Call CUNMBR to compute QB*UB.
* (Workspace in WORK( ITEMP ): need N, prefer N*NB) * (Workspace in WORK( ITEMP ): need N, prefer N*NB)
@ -700,14 +681,14 @@
CALL CGEBRD( M, M, WORK( ILQF ), M, RWORK( ID ), CALL CGEBRD( 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) ITEMP = 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 SBDSVDX( 'U', JOBZ, RNGTGK, M, RWORK( ID ), CALL SBDSVDX( '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( ITEMP ),
$ IWORK, INFO) $ IWORK, INFO)
* *
* If needed, compute left singular vectors. * If needed, compute left singular vectors.
@ -741,7 +722,7 @@
END DO END DO
K = K + M K = K + M
END DO END DO
CALL CLASET( 'A', NS, N-M, CZERO, CZERO, CALL CLASET( 'A', M, N-M, CZERO, CZERO,
$ VT( 1,M+1 ), LDVT ) $ VT( 1,M+1 ), LDVT )
* *
* Call CUNMBR to compute (VB**T)*(PB**T) * Call CUNMBR to compute (VB**T)*(PB**T)
@ -777,14 +758,14 @@
CALL CGEBRD( M, N, A, LDA, RWORK( ID ), RWORK( IE ), CALL CGEBRD( 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) ITEMP = 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 SBDSVDX( 'L', JOBZ, RNGTGK, M, RWORK( ID ), CALL SBDSVDX( '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( ITEMP ),
$ IWORK, INFO) $ IWORK, INFO)
* *
* If needed, compute left singular vectors. * If needed, compute left singular vectors.
@ -818,7 +799,7 @@
END DO END DO
K = K + M K = K + M
END DO END DO
CALL CLASET( 'A', NS, N-M, CZERO, CZERO, CALL CLASET( 'A', M, N-M, CZERO, CZERO,
$ VT( 1,M+1 ), LDVT ) $ VT( 1,M+1 ), LDVT )
* *
* Call CUNMBR to compute VB**T * PB**T * Call CUNMBR to compute VB**T * PB**T

74
lapack-netlib/SRC/dgesvdx.f
View File

@ -169,7 +169,7 @@
*> 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 ILQFin advance and an upper
*> bound must be used. *> bound must be used.
*> \endverbatim *> \endverbatim
*> *>
@ -357,14 +357,8 @@
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 ) THEN ELSE IF( WANTVT .AND. LDVT.LT.MINMN ) THEN
IF( INDS ) THEN INFO = -16
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
@ -386,34 +380,18 @@
* *
* Path 1 (M much larger than N) * Path 1 (M much larger than N)
* *
MAXWRK = N + MAXWRK = N*(N*2+16) +
$ 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*2+20) + 2*N*
$ ILAENV( 1, 'DGEBRD', ' ', N, N, -1, -1 ) ) $ ILAENV( 1, 'DGEBRD', ' ', N, N, -1, -1 ) )
IF (WANTU) THEN MINWRK = N*(N*2+21)
MAXWRK = MAX(MAXWRK,N*(N*3+6)+N*
$ ILAENV( 1, 'DORMQR', ' ', N, N, -1, -1 ) )
END IF
IF (WANTVT) THEN
MAXWRK = MAX(MAXWRK,N*(N*3+6)+N*
$ ILAENV( 1, 'DORMLQ', ' ', N, N, -1, -1 ) )
END IF
MINWRK = N*(N*3+20)
ELSE ELSE
* *
* Path 2 (M at least N, but not much larger) * Path 2 (M at least N, but not much larger)
* *
MAXWRK = 4*N + ( M+N )* MAXWRK = N*(N*2+19) + ( M+N )*
$ ILAENV( 1, 'DGEBRD', ' ', M, N, -1, -1 ) $ ILAENV( 1, 'DGEBRD', ' ', M, N, -1, -1 )
IF (WANTU) THEN MINWRK = N*(N*2+20) + M
MAXWRK = MAX(MAXWRK,N*(N*2+5)+N*
$ ILAENV( 1, 'DORMQR', ' ', N, N, -1, -1 ) )
END IF
IF (WANTVT) THEN
MAXWRK = MAX(MAXWRK,N*(N*2+5)+N*
$ ILAENV( 1, 'DORMLQ', ' ', N, N, -1, -1 ) )
END IF
MINWRK = MAX(N*(N*2+19),4*N+M)
END IF END IF
ELSE ELSE
MNTHR = ILAENV( 6, 'DGESVD', JOBU // JOBVT, M, N, 0, 0 ) MNTHR = ILAENV( 6, 'DGESVD', JOBU // JOBVT, M, N, 0, 0 )
@ -421,34 +399,18 @@
* *
* Path 1t (N much larger than M) * Path 1t (N much larger than M)
* *
MAXWRK = M + MAXWRK = M*(M*2+16) +
$ 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*2+20) + 2*M*
$ ILAENV( 1, 'DGEBRD', ' ', M, M, -1, -1 ) ) $ ILAENV( 1, 'DGEBRD', ' ', M, M, -1, -1 ) )
IF (WANTU) THEN MINWRK = M*(M*2+21)
MAXWRK = MAX(MAXWRK,M*(M*3+6)+M*
$ ILAENV( 1, 'DORMQR', ' ', M, M, -1, -1 ) )
END IF
IF (WANTVT) THEN
MAXWRK = MAX(MAXWRK,M*(M*3+6)+M*
$ ILAENV( 1, 'DORMLQ', ' ', M, M, -1, -1 ) )
END IF
MINWRK = M*(M*3+20)
ELSE ELSE
* *
* Path 2t (N at least M, but not much larger) * Path 2t (N greater than M, but not much larger)
* *
MAXWRK = 4*M + ( M+N )* MAXWRK = M*(M*2+19) + ( M+N )*
$ ILAENV( 1, 'DGEBRD', ' ', M, N, -1, -1 ) $ ILAENV( 1, 'DGEBRD', ' ', M, N, -1, -1 )
IF (WANTU) THEN MINWRK = M*(M*2+20) + N
MAXWRK = MAX(MAXWRK,M*(M*2+5)+M*
$ ILAENV( 1, 'DORMQR', ' ', M, M, -1, -1 ) )
END IF
IF (WANTVT) THEN
MAXWRK = MAX(MAXWRK,M*(M*2+5)+M*
$ ILAENV( 1, 'DORMLQ', ' ', M, M, -1, -1 ) )
END IF
MINWRK = MAX(M*(M*2+19),4*M+N)
END IF END IF
END IF END IF
END IF END IF
@ -560,7 +522,7 @@
CALL DCOPY( N, WORK( J ), 1, U( 1,I ), 1 ) CALL DCOPY( N, WORK( J ), 1, U( 1,I ), 1 )
J = J + N*2 J = J + N*2
END DO END DO
CALL DLASET( 'A', M-N, NS, ZERO, ZERO, U( N+1,1 ), LDU ) CALL DLASET( 'A', M-N, N, ZERO, ZERO, U( N+1,1 ), LDU )
* *
* 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)
@ -629,7 +591,7 @@
CALL DCOPY( N, WORK( J ), 1, U( 1,I ), 1 ) CALL DCOPY( N, WORK( J ), 1, U( 1,I ), 1 )
J = J + N*2 J = J + N*2
END DO END DO
CALL DLASET( 'A', M-N, NS, ZERO, ZERO, U( N+1,1 ), LDU ) CALL DLASET( 'A', M-N, N, ZERO, ZERO, U( N+1,1 ), LDU )
* *
* 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)
@ -725,7 +687,7 @@
CALL DCOPY( M, WORK( J ), 1, VT( I,1 ), LDVT ) CALL DCOPY( M, WORK( J ), 1, VT( I,1 ), LDVT )
J = J + M*2 J = J + M*2
END DO END DO
CALL DLASET( 'A', NS, N-M, ZERO, ZERO, VT( 1,M+1 ), LDVT) CALL DLASET( 'A', M, N-M, ZERO, ZERO, VT( 1,M+1 ), LDVT )
* *
* 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)
@ -794,7 +756,7 @@
CALL DCOPY( M, WORK( J ), 1, VT( I,1 ), LDVT ) CALL DCOPY( M, WORK( J ), 1, VT( I,1 ), LDVT )
J = J + M*2 J = J + M*2
END DO END DO
CALL DLASET( 'A', NS, N-M, ZERO, ZERO, VT( 1,M+1 ), LDVT) CALL DLASET( 'A', M, N-M, ZERO, ZERO, VT( 1,M+1 ), LDVT )
* *
* 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)

74
lapack-netlib/SRC/sgesvdx.f
View File

@ -169,7 +169,7 @@
*> 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 ILQFin advance and an upper
*> bound must be used. *> bound must be used.
*> \endverbatim *> \endverbatim
*> *>
@ -357,14 +357,8 @@
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 ) THEN ELSE IF( WANTVT .AND. LDVT.LT.MINMN ) THEN
IF( INDS ) THEN INFO = -16
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
@ -386,34 +380,18 @@
* *
* Path 1 (M much larger than N) * Path 1 (M much larger than N)
* *
MAXWRK = N + MAXWRK = N*(N*2+16) +
$ N*ILAENV( 1, 'SGEQRF', ' ', M, N, -1, -1 ) $ N*ILAENV( 1, 'SGEQRF', ' ', M, N, -1, -1 )
MAXWRK = MAX( MAXWRK, N*(N+5) + 2*N* MAXWRK = MAX( MAXWRK, N*(N*2+20) + 2*N*
$ ILAENV( 1, 'SGEBRD', ' ', N, N, -1, -1 ) ) $ ILAENV( 1, 'SGEBRD', ' ', N, N, -1, -1 ) )
IF (WANTU) THEN MINWRK = N*(N*2+21)
MAXWRK = MAX(MAXWRK,N*(N*3+6)+N*
$ ILAENV( 1, 'SORMQR', ' ', N, N, -1, -1 ) )
END IF
IF (WANTVT) THEN
MAXWRK = MAX(MAXWRK,N*(N*3+6)+N*
$ ILAENV( 1, 'SORMLQ', ' ', N, N, -1, -1 ) )
END IF
MINWRK = N*(N*3+20)
ELSE ELSE
* *
* Path 2 (M at least N, but not much larger) * Path 2 (M at least N, but not much larger)
* *
MAXWRK = 4*N + ( M+N )* MAXWRK = N*(N*2+19) + ( M+N )*
$ ILAENV( 1, 'SGEBRD', ' ', M, N, -1, -1 ) $ ILAENV( 1, 'SGEBRD', ' ', M, N, -1, -1 )
IF (WANTU) THEN MINWRK = N*(N*2+20) + M
MAXWRK = MAX(MAXWRK,N*(N*2+5)+N*
$ ILAENV( 1, 'SORMQR', ' ', N, N, -1, -1 ) )
END IF
IF (WANTVT) THEN
MAXWRK = MAX(MAXWRK,N*(N*2+5)+N*
$ ILAENV( 1, 'SORMLQ', ' ', N, N, -1, -1 ) )
END IF
MINWRK = MAX(N*(N*2+19),4*N+M)
END IF END IF
ELSE ELSE
MNTHR = ILAENV( 6, 'SGESVD', JOBU // JOBVT, M, N, 0, 0 ) MNTHR = ILAENV( 6, 'SGESVD', JOBU // JOBVT, M, N, 0, 0 )
@ -421,34 +399,18 @@
* *
* Path 1t (N much larger than M) * Path 1t (N much larger than M)
* *
MAXWRK = M + MAXWRK = M*(M*2+16) +
$ M*ILAENV( 1, 'SGELQF', ' ', M, N, -1, -1 ) $ M*ILAENV( 1, 'SGELQF', ' ', M, N, -1, -1 )
MAXWRK = MAX( MAXWRK, M*(M+5) + 2*M* MAXWRK = MAX( MAXWRK, M*(M*2+20) + 2*M*
$ ILAENV( 1, 'SGEBRD', ' ', M, M, -1, -1 ) ) $ ILAENV( 1, 'SGEBRD', ' ', M, M, -1, -1 ) )
IF (WANTU) THEN MINWRK = M*(M*2+21)
MAXWRK = MAX(MAXWRK,M*(M*3+6)+M*
$ ILAENV( 1, 'SORMQR', ' ', M, M, -1, -1 ) )
END IF
IF (WANTVT) THEN
MAXWRK = MAX(MAXWRK,M*(M*3+6)+M*
$ ILAENV( 1, 'SORMLQ', ' ', M, M, -1, -1 ) )
END IF
MINWRK = M*(M*3+20)
ELSE ELSE
* *
* Path 2t (N at least M, but not much larger) * Path 2t (N greater than M, but not much larger)
* *
MAXWRK = 4*M + ( M+N )* MAXWRK = M*(M*2+19) + ( M+N )*
$ ILAENV( 1, 'SGEBRD', ' ', M, N, -1, -1 ) $ ILAENV( 1, 'SGEBRD', ' ', M, N, -1, -1 )
IF (WANTU) THEN MINWRK = M*(M*2+20) + N
MAXWRK = MAX(MAXWRK,M*(M*2+5)+M*
$ ILAENV( 1, 'SORMQR', ' ', M, M, -1, -1 ) )
END IF
IF (WANTVT) THEN
MAXWRK = MAX(MAXWRK,M*(M*2+5)+M*
$ ILAENV( 1, 'SORMLQ', ' ', M, M, -1, -1 ) )
END IF
MINWRK = MAX(M*(M*2+19),4*M+N)
END IF END IF
END IF END IF
END IF END IF
@ -560,7 +522,7 @@
CALL SCOPY( N, WORK( J ), 1, U( 1,I ), 1 ) CALL SCOPY( N, WORK( J ), 1, U( 1,I ), 1 )
J = J + N*2 J = J + N*2
END DO END DO
CALL SLASET( 'A', M-N, NS, ZERO, ZERO, U( N+1,1 ), LDU ) CALL SLASET( 'A', M-N, N, ZERO, ZERO, U( N+1,1 ), LDU )
* *
* Call SORMBR to compute QB*UB. * Call SORMBR to compute QB*UB.
* (Workspace in WORK( ITEMP ): need N, prefer N*NB) * (Workspace in WORK( ITEMP ): need N, prefer N*NB)
@ -629,7 +591,7 @@
CALL SCOPY( N, WORK( J ), 1, U( 1,I ), 1 ) CALL SCOPY( N, WORK( J ), 1, U( 1,I ), 1 )
J = J + N*2 J = J + N*2
END DO END DO
CALL SLASET( 'A', M-N, NS, ZERO, ZERO, U( N+1,1 ), LDU ) CALL SLASET( 'A', M-N, N, ZERO, ZERO, U( N+1,1 ), LDU )
* *
* Call SORMBR to compute QB*UB. * Call SORMBR to compute QB*UB.
* (Workspace in WORK( ITEMP ): need N, prefer N*NB) * (Workspace in WORK( ITEMP ): need N, prefer N*NB)
@ -725,7 +687,7 @@
CALL SCOPY( M, WORK( J ), 1, VT( I,1 ), LDVT ) CALL SCOPY( M, WORK( J ), 1, VT( I,1 ), LDVT )
J = J + M*2 J = J + M*2
END DO END DO
CALL SLASET( 'A', NS, N-M, ZERO, ZERO, VT( 1,M+1 ), LDVT) CALL SLASET( 'A', M, N-M, ZERO, ZERO, VT( 1,M+1 ), LDVT )
* *
* Call SORMBR to compute (VB**T)*(PB**T) * Call SORMBR to compute (VB**T)*(PB**T)
* (Workspace in WORK( ITEMP ): need M, prefer M*NB) * (Workspace in WORK( ITEMP ): need M, prefer M*NB)
@ -794,7 +756,7 @@
CALL SCOPY( M, WORK( J ), 1, VT( I,1 ), LDVT ) CALL SCOPY( M, WORK( J ), 1, VT( I,1 ), LDVT )
J = J + M*2 J = J + M*2
END DO END DO
CALL SLASET( 'A', NS, N-M, ZERO, ZERO, VT( 1,M+1 ), LDVT) CALL SLASET( 'A', M, N-M, ZERO, ZERO, VT( 1,M+1 ), LDVT )
* *
* Call SORMBR to compute VB**T * PB**T * Call SORMBR to compute VB**T * PB**T
* (Workspace in WORK( ITEMP ): need M, prefer M*NB) * (Workspace in WORK( ITEMP ): need M, prefer M*NB)

132
lapack-netlib/SRC/zgesvdx.f
View File

@ -36,30 +36,27 @@
* .. * ..
* *
* *
*> \par Purpose: * 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
*> ZGESVDX computes the singular value decomposition (SVD) of a complex * vectors. The SVD is written
*> M-by-N matrix A, optionally computing the left and/or right singular *
*> 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
*> * min(m,n) diagonal elements, U is an M-by-M unitary matrix, and
*> where SIGMA is an M-by-N matrix which is zero except for its * V is an N-by-N unitary matrix. The diagonal elements of SIGMA
*> min(m,n) diagonal elements, U is an M-by-M unitary matrix, and * are the singular values of A; they are real and non-negative, and
*> V is an N-by-N unitary matrix. The diagonal elements of SIGMA * are returned in descending order. The first min(m,n) columns of
*> are the singular values of A; they are real and non-negative, and * U and V are the left and right singular vectors of A.
*> 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
*> ZGESVDX uses an eigenvalue problem for obtaining the SVD, which * vectors. See DBDSVDX for details.
*> 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.
*>
*> Note that the routine returns V**T, not V.
*> \endverbatim
* *
* Arguments: * Arguments:
* ========== * ==========
@ -110,7 +107,7 @@
*> *>
*> \param[in,out] A *> \param[in,out] A
*> \verbatim *> \verbatim
*> A is COMPLEX*16 array, dimension (LDA,N) *> A is COMPLEX 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
@ -170,7 +167,7 @@
*> 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 ILQFin advance and an upper
*> bound must be used. *> bound must be used.
*> \endverbatim *> \endverbatim
*> *>
@ -294,8 +291,8 @@
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, ITGKZ, IUTGK,
$ IUTGK, J, K, MAXWRK, MINMN, MINWRK, MNTHR $ J, K, MAXWRK, MINMN, MINWRK, MNTHR
DOUBLE PRECISION ABSTOL, ANRM, BIGNUM, EPS, SMLNUM DOUBLE PRECISION ABSTOL, ANRM, BIGNUM, EPS, SMLNUM
* .. * ..
* .. Local Arrays .. * .. Local Arrays ..
@ -367,14 +364,8 @@
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 ) THEN ELSE IF( WANTVT .AND. LDVT.LT.MINMN ) THEN
IF( INDS ) THEN INFO = -16
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
@ -396,24 +387,18 @@
* *
* Path 1 (M much larger than N) * Path 1 (M much larger than N)
* *
MINWRK = N*(N+5) MAXWRK = N + N*
MAXWRK = N + N*ILAENV(1,'ZGEQRF',' ',M,N,-1,-1) $ ILAENV( 1, 'DGEQRF', ' ', M, N, -1, -1 )
MAXWRK = MAX(MAXWRK, MAXWRK = MAX( MAXWRK, N*N + N + 2*N*
$ N*N+2*N+2*N*ILAENV(1,'ZGEBRD',' ',N,N,-1,-1)) $ ILAENV( 1, 'DGEBRD', ' ', N, N, -1, -1 ) )
IF (WANTU .OR. WANTVT) THEN MINWRK = N*(N+4)
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)
* *
MINWRK = 3*N + M MAXWRK = 2*N + ( M+N )*
MAXWRK = 2*N + (M+N)*ILAENV(1,'ZGEBRD',' ',M,N,-1,-1) $ ILAENV( 1, 'ZGEBRD', ' ', M, N, -1, -1 )
IF (WANTU .OR. WANTVT) THEN MINWRK = 2*N + M
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 )
@ -421,25 +406,18 @@
* *
* Path 1t (N much larger than M) * Path 1t (N much larger than M)
* *
MINWRK = M*(M+5) MAXWRK = M + M*
MAXWRK = M + M*ILAENV(1,'ZGELQF',' ',M,N,-1,-1) $ ILAENV( 1, 'ZGELQF', ' ', M, N, -1, -1 )
MAXWRK = MAX(MAXWRK, MAXWRK = MAX( MAXWRK, M*M + M + 2*M*
$ M*M+2*M+2*M*ILAENV(1,'ZGEBRD',' ',M,M,-1,-1)) $ ILAENV( 1, 'ZGEBRD', ' ', M, M, -1, -1 ) )
IF (WANTU .OR. WANTVT) THEN MINWRK = M*(M+4)
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 )*
MINWRK = 3*M + N $ ILAENV( 1, 'ZGEBRD', ' ', M, N, -1, -1 )
MAXWRK = 2*M + (M+N)*ILAENV(1,'ZGEBRD',' ',M,N,-1,-1) MINWRK = 2*M + N
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
@ -537,14 +515,14 @@
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) ITEMP = 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( ITEMP ),
$ IWORK, INFO) $ IWORK, INFO)
* *
* If needed, compute left singular vectors. * If needed, compute left singular vectors.
@ -558,7 +536,7 @@
END DO END DO
K = K + N K = K + N
END DO END DO
CALL ZLASET( 'A', M-N, NS, CZERO, CZERO, U( N+1,1 ), LDU) CALL ZLASET( 'A', M-N, N, 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)
@ -613,14 +591,14 @@
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) ITEMP = 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( ITEMP ),
$ IWORK, INFO) $ IWORK, INFO)
* *
* If needed, compute left singular vectors. * If needed, compute left singular vectors.
@ -634,7 +612,7 @@
END DO END DO
K = K + N K = K + N
END DO END DO
CALL ZLASET( 'A', M-N, NS, CZERO, CZERO, U( N+1,1 ), LDU) CALL ZLASET( 'A', M-N, N, 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)
@ -700,14 +678,14 @@
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) ITEMP = 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( ITEMP ),
$ IWORK, INFO) $ IWORK, INFO)
* *
* If needed, compute left singular vectors. * If needed, compute left singular vectors.
@ -741,7 +719,7 @@
END DO END DO
K = K + M K = K + M
END DO END DO
CALL ZLASET( 'A', NS, N-M, CZERO, CZERO, CALL ZLASET( 'A', M, 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)
@ -777,14 +755,14 @@
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) ITEMP = 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( ITEMP ),
$ IWORK, INFO) $ IWORK, INFO)
* *
* If needed, compute left singular vectors. * If needed, compute left singular vectors.
@ -818,7 +796,7 @@
END DO END DO
K = K + M K = K + M
END DO END DO
CALL ZLASET( 'A', NS, N-M, CZERO, CZERO, CALL ZLASET( 'A', M, 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