version 1.11, 2012/08/22 09:48:30
|
version 1.18, 2017/06/17 11:06:43
|
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 ZGESDD + dependencies |
*> Download ZGESDD + dependencies |
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zgesdd.f"> |
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zgesdd.f"> |
*> [TGZ]</a> |
*> [TGZ]</a> |
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zgesdd.f"> |
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zgesdd.f"> |
*> [ZIP]</a> |
*> [ZIP]</a> |
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zgesdd.f"> |
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zgesdd.f"> |
*> [TXT]</a> |
*> [TXT]</a> |
*> \endhtmlonly |
*> \endhtmlonly |
* |
* |
* Definition: |
* Definition: |
* =========== |
* =========== |
* |
* |
* SUBROUTINE ZGESDD( JOBZ, M, N, A, LDA, S, U, LDU, VT, LDVT, WORK, |
* SUBROUTINE ZGESDD( JOBZ, M, N, A, LDA, S, U, LDU, VT, LDVT, |
* LWORK, RWORK, IWORK, INFO ) |
* WORK, LWORK, RWORK, IWORK, INFO ) |
* |
* |
* .. Scalar Arguments .. |
* .. Scalar Arguments .. |
* CHARACTER JOBZ |
* CHARACTER JOBZ |
* INTEGER INFO, LDA, LDU, LDVT, LWORK, M, N |
* INTEGER INFO, LDA, LDU, LDVT, LWORK, M, N |
Line 31
|
Line 31
|
* COMPLEX*16 A( LDA, * ), U( LDU, * ), VT( LDVT, * ), |
* COMPLEX*16 A( LDA, * ), U( LDU, * ), VT( LDVT, * ), |
* $ WORK( * ) |
* $ WORK( * ) |
* .. |
* .. |
* |
* |
* |
* |
*> \par Purpose: |
*> \par Purpose: |
* ============= |
* ============= |
Line 135
|
Line 135
|
*> \param[in] LDU |
*> \param[in] LDU |
*> \verbatim |
*> \verbatim |
*> LDU is INTEGER |
*> LDU is INTEGER |
*> The leading dimension of the array U. LDU >= 1; if |
*> The leading dimension of the array U. LDU >= 1; |
*> JOBZ = 'S' or 'A' or JOBZ = 'O' and M < N, LDU >= M. |
*> if JOBZ = 'S' or 'A' or JOBZ = 'O' and M < N, LDU >= M. |
*> \endverbatim |
*> \endverbatim |
*> |
*> |
*> \param[out] VT |
*> \param[out] VT |
Line 152
|
Line 152
|
*> \param[in] LDVT |
*> \param[in] LDVT |
*> \verbatim |
*> \verbatim |
*> LDVT is INTEGER |
*> LDVT is INTEGER |
*> The leading dimension of the array VT. LDVT >= 1; if |
*> The leading dimension of the array VT. LDVT >= 1; |
*> JOBZ = 'A' or JOBZ = 'O' and M >= N, LDVT >= N; |
*> if JOBZ = 'A' or JOBZ = 'O' and M >= N, LDVT >= N; |
*> if JOBZ = 'S', LDVT >= min(M,N). |
*> if JOBZ = 'S', LDVT >= min(M,N). |
*> \endverbatim |
*> \endverbatim |
*> |
*> |
Line 167
|
Line 167
|
*> \verbatim |
*> \verbatim |
*> LWORK is INTEGER |
*> LWORK is INTEGER |
*> The dimension of the array WORK. LWORK >= 1. |
*> The dimension of the array WORK. LWORK >= 1. |
*> if JOBZ = 'N', LWORK >= 2*min(M,N)+max(M,N). |
|
*> if JOBZ = 'O', |
|
*> LWORK >= 2*min(M,N)*min(M,N)+2*min(M,N)+max(M,N). |
|
*> if JOBZ = 'S' or 'A', |
|
*> LWORK >= min(M,N)*min(M,N)+2*min(M,N)+max(M,N). |
|
*> For good performance, LWORK should generally be larger. |
|
*> |
|
*> If LWORK = -1, a workspace query is assumed. The optimal |
*> If LWORK = -1, a workspace query is assumed. The optimal |
*> size for the WORK array is calculated and stored in WORK(1), |
*> size for the WORK array is calculated and stored in WORK(1), |
*> and no other work except argument checking is performed. |
*> and no other work except argument checking is performed. |
|
*> |
|
*> Let mx = max(M,N) and mn = min(M,N). |
|
*> If JOBZ = 'N', LWORK >= 2*mn + mx. |
|
*> If JOBZ = 'O', LWORK >= 2*mn*mn + 2*mn + mx. |
|
*> If JOBZ = 'S', LWORK >= mn*mn + 3*mn. |
|
*> If JOBZ = 'A', LWORK >= mn*mn + 2*mn + mx. |
|
*> These are not tight minimums in all cases; see comments inside code. |
|
*> For good performance, LWORK should generally be larger; |
|
*> a query is recommended. |
*> \endverbatim |
*> \endverbatim |
*> |
*> |
*> \param[out] RWORK |
*> \param[out] RWORK |
*> \verbatim |
*> \verbatim |
*> RWORK is DOUBLE PRECISION array, dimension (MAX(1,LRWORK)) |
*> RWORK is DOUBLE PRECISION array, dimension (MAX(1,LRWORK)) |
*> If JOBZ = 'N', LRWORK >= 5*min(M,N). |
*> Let mx = max(M,N) and mn = min(M,N). |
*> Otherwise, |
*> If JOBZ = 'N', LRWORK >= 5*mn (LAPACK <= 3.6 needs 7*mn); |
*> LRWORK >= min(M,N)*max(5*min(M,N)+7,2*max(M,N)+2*min(M,N)+1) |
*> else if mx >> mn, LRWORK >= 5*mn*mn + 5*mn; |
|
*> else LRWORK >= max( 5*mn*mn + 5*mn, |
|
*> 2*mx*mn + 2*mn*mn + mn ). |
*> \endverbatim |
*> \endverbatim |
*> |
*> |
*> \param[out] IWORK |
*> \param[out] IWORK |
Line 203
|
Line 207
|
* 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 2011 |
*> \date June 2016 |
* |
* |
*> \ingroup complex16GEsing |
*> \ingroup complex16GEsing |
* |
* |
Line 219
|
Line 223
|
*> California at Berkeley, USA |
*> California at Berkeley, USA |
*> |
*> |
* ===================================================================== |
* ===================================================================== |
SUBROUTINE ZGESDD( JOBZ, M, N, A, LDA, S, U, LDU, VT, LDVT, WORK, |
SUBROUTINE ZGESDD( JOBZ, M, N, A, LDA, S, U, LDU, VT, LDVT, |
$ LWORK, RWORK, IWORK, INFO ) |
$ WORK, LWORK, RWORK, IWORK, INFO ) |
|
implicit none |
* |
* |
* -- LAPACK driver routine (version 3.4.0) -- |
* -- 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..-- |
* November 2011 |
* June 2016 |
* |
* |
* .. Scalar Arguments .. |
* .. Scalar Arguments .. |
CHARACTER JOBZ |
CHARACTER JOBZ |
Line 241
|
Line 246
|
* ===================================================================== |
* ===================================================================== |
* |
* |
* .. Parameters .. |
* .. Parameters .. |
INTEGER LQUERV |
|
PARAMETER ( LQUERV = -1 ) |
|
COMPLEX*16 CZERO, CONE |
COMPLEX*16 CZERO, CONE |
PARAMETER ( CZERO = ( 0.0D+0, 0.0D+0 ), |
PARAMETER ( CZERO = ( 0.0D+0, 0.0D+0 ), |
$ CONE = ( 1.0D+0, 0.0D+0 ) ) |
$ CONE = ( 1.0D+0, 0.0D+0 ) ) |
Line 250
|
Line 253
|
PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 ) |
PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 ) |
* .. |
* .. |
* .. Local Scalars .. |
* .. Local Scalars .. |
LOGICAL WNTQA, WNTQAS, WNTQN, WNTQO, WNTQS |
LOGICAL LQUERY, WNTQA, WNTQAS, WNTQN, WNTQO, WNTQS |
INTEGER BLK, CHUNK, I, IE, IERR, IL, IR, IRU, IRVT, |
INTEGER BLK, CHUNK, I, IE, IERR, IL, IR, IRU, IRVT, |
$ ISCL, ITAU, ITAUP, ITAUQ, IU, IVT, LDWKVT, |
$ ISCL, ITAU, ITAUP, ITAUQ, IU, IVT, LDWKVT, |
$ LDWRKL, LDWRKR, LDWRKU, MAXWRK, MINMN, MINWRK, |
$ LDWRKL, LDWRKR, LDWRKU, MAXWRK, MINMN, MINWRK, |
$ MNTHR1, MNTHR2, NRWORK, NWORK, WRKBL |
$ MNTHR1, MNTHR2, NRWORK, NWORK, WRKBL |
|
INTEGER LWORK_ZGEBRD_MN, LWORK_ZGEBRD_MM, |
|
$ LWORK_ZGEBRD_NN, LWORK_ZGELQF_MN, |
|
$ LWORK_ZGEQRF_MN, |
|
$ LWORK_ZUNGBR_P_MN, LWORK_ZUNGBR_P_NN, |
|
$ LWORK_ZUNGBR_Q_MN, LWORK_ZUNGBR_Q_MM, |
|
$ LWORK_ZUNGLQ_MN, LWORK_ZUNGLQ_NN, |
|
$ LWORK_ZUNGQR_MM, LWORK_ZUNGQR_MN, |
|
$ LWORK_ZUNMBR_PRC_MM, LWORK_ZUNMBR_QLN_MM, |
|
$ LWORK_ZUNMBR_PRC_MN, LWORK_ZUNMBR_QLN_MN, |
|
$ LWORK_ZUNMBR_PRC_NN, LWORK_ZUNMBR_QLN_NN |
DOUBLE PRECISION ANRM, BIGNUM, EPS, SMLNUM |
DOUBLE PRECISION ANRM, BIGNUM, EPS, SMLNUM |
* .. |
* .. |
* .. Local Arrays .. |
* .. Local Arrays .. |
INTEGER IDUM( 1 ) |
INTEGER IDUM( 1 ) |
DOUBLE PRECISION DUM( 1 ) |
DOUBLE PRECISION DUM( 1 ) |
|
COMPLEX*16 CDUM( 1 ) |
* .. |
* .. |
* .. External Subroutines .. |
* .. External Subroutines .. |
EXTERNAL DBDSDC, DLASCL, XERBLA, ZGEBRD, ZGELQF, ZGEMM, |
EXTERNAL DBDSDC, DLASCL, XERBLA, ZGEBRD, ZGELQF, ZGEMM, |
Line 268
|
Line 282
|
* .. |
* .. |
* .. External Functions .. |
* .. External Functions .. |
LOGICAL LSAME |
LOGICAL LSAME |
INTEGER ILAENV |
|
DOUBLE PRECISION DLAMCH, ZLANGE |
DOUBLE PRECISION DLAMCH, ZLANGE |
EXTERNAL LSAME, ILAENV, DLAMCH, ZLANGE |
EXTERNAL LSAME, DLAMCH, ZLANGE |
* .. |
* .. |
* .. Intrinsic Functions .. |
* .. Intrinsic Functions .. |
INTRINSIC INT, MAX, MIN, SQRT |
INTRINSIC INT, MAX, MIN, SQRT |
Line 279
|
Line 292
|
* |
* |
* Test the input arguments |
* Test the input arguments |
* |
* |
INFO = 0 |
INFO = 0 |
MINMN = MIN( M, N ) |
MINMN = MIN( M, N ) |
MNTHR1 = INT( MINMN*17.0D0 / 9.0D0 ) |
MNTHR1 = INT( MINMN*17.0D0 / 9.0D0 ) |
MNTHR2 = INT( MINMN*5.0D0 / 3.0D0 ) |
MNTHR2 = INT( MINMN*5.0D0 / 3.0D0 ) |
WNTQA = LSAME( JOBZ, 'A' ) |
WNTQA = LSAME( JOBZ, 'A' ) |
WNTQS = LSAME( JOBZ, 'S' ) |
WNTQS = LSAME( JOBZ, 'S' ) |
WNTQAS = WNTQA .OR. WNTQS |
WNTQAS = WNTQA .OR. WNTQS |
WNTQO = LSAME( JOBZ, 'O' ) |
WNTQO = LSAME( JOBZ, 'O' ) |
WNTQN = LSAME( JOBZ, 'N' ) |
WNTQN = LSAME( JOBZ, 'N' ) |
|
LQUERY = ( LWORK.EQ.-1 ) |
MINWRK = 1 |
MINWRK = 1 |
MAXWRK = 1 |
MAXWRK = 1 |
* |
* |
Line 309
|
Line 323
|
END IF |
END IF |
* |
* |
* Compute workspace |
* Compute workspace |
* (Note: Comments in the code beginning "Workspace:" describe the |
* Note: Comments in the code beginning "Workspace:" describe the |
* minimal amount of workspace needed at that point in the code, |
* minimal amount of workspace allocated at that point in the code, |
* as well as the preferred amount for good performance. |
* as well as the preferred amount for good performance. |
* CWorkspace refers to complex workspace, and RWorkspace to |
* CWorkspace refers to complex workspace, and RWorkspace to |
* real workspace. NB refers to the optimal block size for the |
* real workspace. NB refers to the optimal block size for the |
* immediately following subroutine, as returned by ILAENV.) |
* immediately following subroutine, as returned by ILAENV.) |
* |
* |
IF( INFO.EQ.0 .AND. M.GT.0 .AND. N.GT.0 ) THEN |
IF( INFO.EQ.0 ) THEN |
IF( M.GE.N ) THEN |
MINWRK = 1 |
|
MAXWRK = 1 |
|
IF( M.GE.N .AND. MINMN.GT.0 ) THEN |
* |
* |
* There is no complex work space needed for bidiagonal SVD |
* There is no complex work space needed for bidiagonal SVD |
* The real work space needed for bidiagonal SVD is BDSPAC |
* The real work space needed for bidiagonal SVD (dbdsdc) is |
* for computing singular values and singular vectors; BDSPAN |
* BDSPAC = 3*N*N + 4*N for singular values and vectors; |
* for computing singular values only. |
* BDSPAC = 4*N for singular values only; |
* BDSPAC = 5*N*N + 7*N |
* not including e, RU, and RVT matrices. |
* BDSPAN = MAX(7*N+4, 3*N+2+SMLSIZ*(SMLSIZ+8)) |
* |
|
* Compute space preferred for each routine |
|
CALL ZGEBRD( M, N, CDUM(1), M, DUM(1), DUM(1), CDUM(1), |
|
$ CDUM(1), CDUM(1), -1, IERR ) |
|
LWORK_ZGEBRD_MN = INT( CDUM(1) ) |
|
* |
|
CALL ZGEBRD( N, N, CDUM(1), N, DUM(1), DUM(1), CDUM(1), |
|
$ CDUM(1), CDUM(1), -1, IERR ) |
|
LWORK_ZGEBRD_NN = INT( CDUM(1) ) |
|
* |
|
CALL ZGEQRF( M, N, CDUM(1), M, CDUM(1), CDUM(1), -1, IERR ) |
|
LWORK_ZGEQRF_MN = INT( CDUM(1) ) |
|
* |
|
CALL ZUNGBR( 'P', N, N, N, CDUM(1), N, CDUM(1), CDUM(1), |
|
$ -1, IERR ) |
|
LWORK_ZUNGBR_P_NN = INT( CDUM(1) ) |
|
* |
|
CALL ZUNGBR( 'Q', M, M, N, CDUM(1), M, CDUM(1), CDUM(1), |
|
$ -1, IERR ) |
|
LWORK_ZUNGBR_Q_MM = INT( CDUM(1) ) |
|
* |
|
CALL ZUNGBR( 'Q', M, N, N, CDUM(1), M, CDUM(1), CDUM(1), |
|
$ -1, IERR ) |
|
LWORK_ZUNGBR_Q_MN = INT( CDUM(1) ) |
|
* |
|
CALL ZUNGQR( M, M, N, CDUM(1), M, CDUM(1), CDUM(1), |
|
$ -1, IERR ) |
|
LWORK_ZUNGQR_MM = INT( CDUM(1) ) |
|
* |
|
CALL ZUNGQR( M, N, N, CDUM(1), M, CDUM(1), CDUM(1), |
|
$ -1, IERR ) |
|
LWORK_ZUNGQR_MN = INT( CDUM(1) ) |
|
* |
|
CALL ZUNMBR( 'P', 'R', 'C', N, N, N, CDUM(1), N, CDUM(1), |
|
$ CDUM(1), N, CDUM(1), -1, IERR ) |
|
LWORK_ZUNMBR_PRC_NN = INT( CDUM(1) ) |
|
* |
|
CALL ZUNMBR( 'Q', 'L', 'N', M, M, N, CDUM(1), M, CDUM(1), |
|
$ CDUM(1), M, CDUM(1), -1, IERR ) |
|
LWORK_ZUNMBR_QLN_MM = INT( CDUM(1) ) |
|
* |
|
CALL ZUNMBR( 'Q', 'L', 'N', M, N, N, CDUM(1), M, CDUM(1), |
|
$ CDUM(1), M, CDUM(1), -1, IERR ) |
|
LWORK_ZUNMBR_QLN_MN = INT( CDUM(1) ) |
|
* |
|
CALL ZUNMBR( 'Q', 'L', 'N', N, N, N, CDUM(1), N, CDUM(1), |
|
$ CDUM(1), N, CDUM(1), -1, IERR ) |
|
LWORK_ZUNMBR_QLN_NN = INT( CDUM(1) ) |
* |
* |
IF( M.GE.MNTHR1 ) THEN |
IF( M.GE.MNTHR1 ) THEN |
IF( WNTQN ) THEN |
IF( WNTQN ) THEN |
* |
* |
* Path 1 (M much larger than N, JOBZ='N') |
* Path 1 (M >> N, JOBZ='N') |
* |
* |
MAXWRK = N + N*ILAENV( 1, 'ZGEQRF', ' ', M, N, -1, |
MAXWRK = N + LWORK_ZGEQRF_MN |
$ -1 ) |
MAXWRK = MAX( MAXWRK, 2*N + LWORK_ZGEBRD_NN ) |
MAXWRK = MAX( MAXWRK, 2*N+2*N* |
|
$ ILAENV( 1, 'ZGEBRD', ' ', N, N, -1, -1 ) ) |
|
MINWRK = 3*N |
MINWRK = 3*N |
ELSE IF( WNTQO ) THEN |
ELSE IF( WNTQO ) THEN |
* |
* |
* Path 2 (M much larger than N, JOBZ='O') |
* Path 2 (M >> N, JOBZ='O') |
* |
* |
WRKBL = N + N*ILAENV( 1, 'ZGEQRF', ' ', M, N, -1, -1 ) |
WRKBL = N + LWORK_ZGEQRF_MN |
WRKBL = MAX( WRKBL, N+N*ILAENV( 1, 'ZUNGQR', ' ', M, |
WRKBL = MAX( WRKBL, N + LWORK_ZUNGQR_MN ) |
$ N, N, -1 ) ) |
WRKBL = MAX( WRKBL, 2*N + LWORK_ZGEBRD_NN ) |
WRKBL = MAX( WRKBL, 2*N+2*N* |
WRKBL = MAX( WRKBL, 2*N + LWORK_ZUNMBR_QLN_NN ) |
$ ILAENV( 1, 'ZGEBRD', ' ', N, N, -1, -1 ) ) |
WRKBL = MAX( WRKBL, 2*N + LWORK_ZUNMBR_PRC_NN ) |
WRKBL = MAX( WRKBL, 2*N+N* |
|
$ ILAENV( 1, 'ZUNMBR', 'QLN', N, N, N, -1 ) ) |
|
WRKBL = MAX( WRKBL, 2*N+N* |
|
$ ILAENV( 1, 'ZUNMBR', 'PRC', N, N, N, -1 ) ) |
|
MAXWRK = M*N + N*N + WRKBL |
MAXWRK = M*N + N*N + WRKBL |
MINWRK = 2*N*N + 3*N |
MINWRK = 2*N*N + 3*N |
ELSE IF( WNTQS ) THEN |
ELSE IF( WNTQS ) THEN |
* |
* |
* Path 3 (M much larger than N, JOBZ='S') |
* Path 3 (M >> N, JOBZ='S') |
* |
* |
WRKBL = N + N*ILAENV( 1, 'ZGEQRF', ' ', M, N, -1, -1 ) |
WRKBL = N + LWORK_ZGEQRF_MN |
WRKBL = MAX( WRKBL, N+N*ILAENV( 1, 'ZUNGQR', ' ', M, |
WRKBL = MAX( WRKBL, N + LWORK_ZUNGQR_MN ) |
$ N, N, -1 ) ) |
WRKBL = MAX( WRKBL, 2*N + LWORK_ZGEBRD_NN ) |
WRKBL = MAX( WRKBL, 2*N+2*N* |
WRKBL = MAX( WRKBL, 2*N + LWORK_ZUNMBR_QLN_NN ) |
$ ILAENV( 1, 'ZGEBRD', ' ', N, N, -1, -1 ) ) |
WRKBL = MAX( WRKBL, 2*N + LWORK_ZUNMBR_PRC_NN ) |
WRKBL = MAX( WRKBL, 2*N+N* |
|
$ ILAENV( 1, 'ZUNMBR', 'QLN', N, N, N, -1 ) ) |
|
WRKBL = MAX( WRKBL, 2*N+N* |
|
$ ILAENV( 1, 'ZUNMBR', 'PRC', N, N, N, -1 ) ) |
|
MAXWRK = N*N + WRKBL |
MAXWRK = N*N + WRKBL |
MINWRK = N*N + 3*N |
MINWRK = N*N + 3*N |
ELSE IF( WNTQA ) THEN |
ELSE IF( WNTQA ) THEN |
* |
* |
* Path 4 (M much larger than N, JOBZ='A') |
* Path 4 (M >> N, JOBZ='A') |
* |
* |
WRKBL = N + N*ILAENV( 1, 'ZGEQRF', ' ', M, N, -1, -1 ) |
WRKBL = N + LWORK_ZGEQRF_MN |
WRKBL = MAX( WRKBL, N+M*ILAENV( 1, 'ZUNGQR', ' ', M, |
WRKBL = MAX( WRKBL, N + LWORK_ZUNGQR_MM ) |
$ M, N, -1 ) ) |
WRKBL = MAX( WRKBL, 2*N + LWORK_ZGEBRD_NN ) |
WRKBL = MAX( WRKBL, 2*N+2*N* |
WRKBL = MAX( WRKBL, 2*N + LWORK_ZUNMBR_QLN_NN ) |
$ ILAENV( 1, 'ZGEBRD', ' ', N, N, -1, -1 ) ) |
WRKBL = MAX( WRKBL, 2*N + LWORK_ZUNMBR_PRC_NN ) |
WRKBL = MAX( WRKBL, 2*N+N* |
|
$ ILAENV( 1, 'ZUNMBR', 'QLN', N, N, N, -1 ) ) |
|
WRKBL = MAX( WRKBL, 2*N+N* |
|
$ ILAENV( 1, 'ZUNMBR', 'PRC', N, N, N, -1 ) ) |
|
MAXWRK = N*N + WRKBL |
MAXWRK = N*N + WRKBL |
MINWRK = N*N + 2*N + M |
MINWRK = N*N + MAX( 3*N, N + M ) |
END IF |
END IF |
ELSE IF( M.GE.MNTHR2 ) THEN |
ELSE IF( M.GE.MNTHR2 ) THEN |
* |
* |
* Path 5 (M much larger than N, but not as much as MNTHR1) |
* Path 5 (M >> N, but not as much as MNTHR1) |
* |
* |
MAXWRK = 2*N + ( M+N )*ILAENV( 1, 'ZGEBRD', ' ', M, N, |
MAXWRK = 2*N + LWORK_ZGEBRD_MN |
$ -1, -1 ) |
|
MINWRK = 2*N + M |
MINWRK = 2*N + M |
IF( WNTQO ) THEN |
IF( WNTQO ) THEN |
MAXWRK = MAX( MAXWRK, 2*N+N* |
* Path 5o (M >> N, JOBZ='O') |
$ ILAENV( 1, 'ZUNGBR', 'P', N, N, N, -1 ) ) |
MAXWRK = MAX( MAXWRK, 2*N + LWORK_ZUNGBR_P_NN ) |
MAXWRK = MAX( MAXWRK, 2*N+N* |
MAXWRK = MAX( MAXWRK, 2*N + LWORK_ZUNGBR_Q_MN ) |
$ ILAENV( 1, 'ZUNGBR', 'Q', M, N, N, -1 ) ) |
|
MAXWRK = MAXWRK + M*N |
MAXWRK = MAXWRK + M*N |
MINWRK = MINWRK + N*N |
MINWRK = MINWRK + N*N |
ELSE IF( WNTQS ) THEN |
ELSE IF( WNTQS ) THEN |
MAXWRK = MAX( MAXWRK, 2*N+N* |
* Path 5s (M >> N, JOBZ='S') |
$ ILAENV( 1, 'ZUNGBR', 'P', N, N, N, -1 ) ) |
MAXWRK = MAX( MAXWRK, 2*N + LWORK_ZUNGBR_P_NN ) |
MAXWRK = MAX( MAXWRK, 2*N+N* |
MAXWRK = MAX( MAXWRK, 2*N + LWORK_ZUNGBR_Q_MN ) |
$ ILAENV( 1, 'ZUNGBR', 'Q', M, N, N, -1 ) ) |
|
ELSE IF( WNTQA ) THEN |
ELSE IF( WNTQA ) THEN |
MAXWRK = MAX( MAXWRK, 2*N+N* |
* Path 5a (M >> N, JOBZ='A') |
$ ILAENV( 1, 'ZUNGBR', 'P', N, N, N, -1 ) ) |
MAXWRK = MAX( MAXWRK, 2*N + LWORK_ZUNGBR_P_NN ) |
MAXWRK = MAX( MAXWRK, 2*N+M* |
MAXWRK = MAX( MAXWRK, 2*N + LWORK_ZUNGBR_Q_MM ) |
$ ILAENV( 1, 'ZUNGBR', 'Q', M, M, N, -1 ) ) |
|
END IF |
END IF |
ELSE |
ELSE |
* |
* |
* Path 6 (M at least N, but not much larger) |
* Path 6 (M >= N, but not much larger) |
* |
* |
MAXWRK = 2*N + ( M+N )*ILAENV( 1, 'ZGEBRD', ' ', M, N, |
MAXWRK = 2*N + LWORK_ZGEBRD_MN |
$ -1, -1 ) |
|
MINWRK = 2*N + M |
MINWRK = 2*N + M |
IF( WNTQO ) THEN |
IF( WNTQO ) THEN |
MAXWRK = MAX( MAXWRK, 2*N+N* |
* Path 6o (M >= N, JOBZ='O') |
$ ILAENV( 1, 'ZUNMBR', 'PRC', N, N, N, -1 ) ) |
MAXWRK = MAX( MAXWRK, 2*N + LWORK_ZUNMBR_PRC_NN ) |
MAXWRK = MAX( MAXWRK, 2*N+N* |
MAXWRK = MAX( MAXWRK, 2*N + LWORK_ZUNMBR_QLN_MN ) |
$ ILAENV( 1, 'ZUNMBR', 'QLN', M, N, N, -1 ) ) |
|
MAXWRK = MAXWRK + M*N |
MAXWRK = MAXWRK + M*N |
MINWRK = MINWRK + N*N |
MINWRK = MINWRK + N*N |
ELSE IF( WNTQS ) THEN |
ELSE IF( WNTQS ) THEN |
MAXWRK = MAX( MAXWRK, 2*N+N* |
* Path 6s (M >= N, JOBZ='S') |
$ ILAENV( 1, 'ZUNMBR', 'PRC', N, N, N, -1 ) ) |
MAXWRK = MAX( MAXWRK, 2*N + LWORK_ZUNMBR_QLN_MN ) |
MAXWRK = MAX( MAXWRK, 2*N+N* |
MAXWRK = MAX( MAXWRK, 2*N + LWORK_ZUNMBR_PRC_NN ) |
$ ILAENV( 1, 'ZUNMBR', 'QLN', M, N, N, -1 ) ) |
|
ELSE IF( WNTQA ) THEN |
ELSE IF( WNTQA ) THEN |
MAXWRK = MAX( MAXWRK, 2*N+N* |
* Path 6a (M >= N, JOBZ='A') |
$ ILAENV( 1, 'ZUNGBR', 'PRC', N, N, N, -1 ) ) |
MAXWRK = MAX( MAXWRK, 2*N + LWORK_ZUNMBR_QLN_MM ) |
MAXWRK = MAX( MAXWRK, 2*N+M* |
MAXWRK = MAX( MAXWRK, 2*N + LWORK_ZUNMBR_PRC_NN ) |
$ ILAENV( 1, 'ZUNGBR', 'QLN', M, M, N, -1 ) ) |
|
END IF |
END IF |
END IF |
END IF |
ELSE |
ELSE IF( MINMN.GT.0 ) THEN |
* |
* |
* There is no complex work space needed for bidiagonal SVD |
* There is no complex work space needed for bidiagonal SVD |
* The real work space needed for bidiagonal SVD is BDSPAC |
* The real work space needed for bidiagonal SVD (dbdsdc) is |
* for computing singular values and singular vectors; BDSPAN |
* BDSPAC = 3*M*M + 4*M for singular values and vectors; |
* for computing singular values only. |
* BDSPAC = 4*M for singular values only; |
* BDSPAC = 5*M*M + 7*M |
* not including e, RU, and RVT matrices. |
* BDSPAN = MAX(7*M+4, 3*M+2+SMLSIZ*(SMLSIZ+8)) |
* |
|
* Compute space preferred for each routine |
|
CALL ZGEBRD( M, N, CDUM(1), M, DUM(1), DUM(1), CDUM(1), |
|
$ CDUM(1), CDUM(1), -1, IERR ) |
|
LWORK_ZGEBRD_MN = INT( CDUM(1) ) |
|
* |
|
CALL ZGEBRD( M, M, CDUM(1), M, DUM(1), DUM(1), CDUM(1), |
|
$ CDUM(1), CDUM(1), -1, IERR ) |
|
LWORK_ZGEBRD_MM = INT( CDUM(1) ) |
|
* |
|
CALL ZGELQF( M, N, CDUM(1), M, CDUM(1), CDUM(1), -1, IERR ) |
|
LWORK_ZGELQF_MN = INT( CDUM(1) ) |
|
* |
|
CALL ZUNGBR( 'P', M, N, M, CDUM(1), M, CDUM(1), CDUM(1), |
|
$ -1, IERR ) |
|
LWORK_ZUNGBR_P_MN = INT( CDUM(1) ) |
|
* |
|
CALL ZUNGBR( 'P', N, N, M, CDUM(1), N, CDUM(1), CDUM(1), |
|
$ -1, IERR ) |
|
LWORK_ZUNGBR_P_NN = INT( CDUM(1) ) |
|
* |
|
CALL ZUNGBR( 'Q', M, M, N, CDUM(1), M, CDUM(1), CDUM(1), |
|
$ -1, IERR ) |
|
LWORK_ZUNGBR_Q_MM = INT( CDUM(1) ) |
|
* |
|
CALL ZUNGLQ( M, N, M, CDUM(1), M, CDUM(1), CDUM(1), |
|
$ -1, IERR ) |
|
LWORK_ZUNGLQ_MN = INT( CDUM(1) ) |
|
* |
|
CALL ZUNGLQ( N, N, M, CDUM(1), N, CDUM(1), CDUM(1), |
|
$ -1, IERR ) |
|
LWORK_ZUNGLQ_NN = INT( CDUM(1) ) |
|
* |
|
CALL ZUNMBR( 'P', 'R', 'C', M, M, M, CDUM(1), M, CDUM(1), |
|
$ CDUM(1), M, CDUM(1), -1, IERR ) |
|
LWORK_ZUNMBR_PRC_MM = INT( CDUM(1) ) |
|
* |
|
CALL ZUNMBR( 'P', 'R', 'C', M, N, M, CDUM(1), M, CDUM(1), |
|
$ CDUM(1), M, CDUM(1), -1, IERR ) |
|
LWORK_ZUNMBR_PRC_MN = INT( CDUM(1) ) |
|
* |
|
CALL ZUNMBR( 'P', 'R', 'C', N, N, M, CDUM(1), N, CDUM(1), |
|
$ CDUM(1), N, CDUM(1), -1, IERR ) |
|
LWORK_ZUNMBR_PRC_NN = INT( CDUM(1) ) |
|
* |
|
CALL ZUNMBR( 'Q', 'L', 'N', M, M, M, CDUM(1), M, CDUM(1), |
|
$ CDUM(1), M, CDUM(1), -1, IERR ) |
|
LWORK_ZUNMBR_QLN_MM = INT( CDUM(1) ) |
* |
* |
IF( N.GE.MNTHR1 ) THEN |
IF( N.GE.MNTHR1 ) THEN |
IF( WNTQN ) THEN |
IF( WNTQN ) THEN |
* |
* |
* Path 1t (N much larger than M, JOBZ='N') |
* Path 1t (N >> M, JOBZ='N') |
* |
* |
MAXWRK = M + M*ILAENV( 1, 'ZGELQF', ' ', M, N, -1, |
MAXWRK = M + LWORK_ZGELQF_MN |
$ -1 ) |
MAXWRK = MAX( MAXWRK, 2*M + LWORK_ZGEBRD_MM ) |
MAXWRK = MAX( MAXWRK, 2*M+2*M* |
|
$ ILAENV( 1, 'ZGEBRD', ' ', M, M, -1, -1 ) ) |
|
MINWRK = 3*M |
MINWRK = 3*M |
ELSE IF( WNTQO ) THEN |
ELSE IF( WNTQO ) THEN |
* |
* |
* Path 2t (N much larger than M, JOBZ='O') |
* Path 2t (N >> M, JOBZ='O') |
* |
* |
WRKBL = M + M*ILAENV( 1, 'ZGELQF', ' ', M, N, -1, -1 ) |
WRKBL = M + LWORK_ZGELQF_MN |
WRKBL = MAX( WRKBL, M+M*ILAENV( 1, 'ZUNGLQ', ' ', M, |
WRKBL = MAX( WRKBL, M + LWORK_ZUNGLQ_MN ) |
$ N, M, -1 ) ) |
WRKBL = MAX( WRKBL, 2*M + LWORK_ZGEBRD_MM ) |
WRKBL = MAX( WRKBL, 2*M+2*M* |
WRKBL = MAX( WRKBL, 2*M + LWORK_ZUNMBR_QLN_MM ) |
$ ILAENV( 1, 'ZGEBRD', ' ', M, M, -1, -1 ) ) |
WRKBL = MAX( WRKBL, 2*M + LWORK_ZUNMBR_PRC_MM ) |
WRKBL = MAX( WRKBL, 2*M+M* |
|
$ ILAENV( 1, 'ZUNMBR', 'PRC', M, M, M, -1 ) ) |
|
WRKBL = MAX( WRKBL, 2*M+M* |
|
$ ILAENV( 1, 'ZUNMBR', 'QLN', M, M, M, -1 ) ) |
|
MAXWRK = M*N + M*M + WRKBL |
MAXWRK = M*N + M*M + WRKBL |
MINWRK = 2*M*M + 3*M |
MINWRK = 2*M*M + 3*M |
ELSE IF( WNTQS ) THEN |
ELSE IF( WNTQS ) THEN |
* |
* |
* Path 3t (N much larger than M, JOBZ='S') |
* Path 3t (N >> M, JOBZ='S') |
* |
* |
WRKBL = M + M*ILAENV( 1, 'ZGELQF', ' ', M, N, -1, -1 ) |
WRKBL = M + LWORK_ZGELQF_MN |
WRKBL = MAX( WRKBL, M+M*ILAENV( 1, 'ZUNGLQ', ' ', M, |
WRKBL = MAX( WRKBL, M + LWORK_ZUNGLQ_MN ) |
$ N, M, -1 ) ) |
WRKBL = MAX( WRKBL, 2*M + LWORK_ZGEBRD_MM ) |
WRKBL = MAX( WRKBL, 2*M+2*M* |
WRKBL = MAX( WRKBL, 2*M + LWORK_ZUNMBR_QLN_MM ) |
$ ILAENV( 1, 'ZGEBRD', ' ', M, M, -1, -1 ) ) |
WRKBL = MAX( WRKBL, 2*M + LWORK_ZUNMBR_PRC_MM ) |
WRKBL = MAX( WRKBL, 2*M+M* |
|
$ ILAENV( 1, 'ZUNMBR', 'PRC', M, M, M, -1 ) ) |
|
WRKBL = MAX( WRKBL, 2*M+M* |
|
$ ILAENV( 1, 'ZUNMBR', 'QLN', M, M, M, -1 ) ) |
|
MAXWRK = M*M + WRKBL |
MAXWRK = M*M + WRKBL |
MINWRK = M*M + 3*M |
MINWRK = M*M + 3*M |
ELSE IF( WNTQA ) THEN |
ELSE IF( WNTQA ) THEN |
* |
* |
* Path 4t (N much larger than M, JOBZ='A') |
* Path 4t (N >> M, JOBZ='A') |
* |
* |
WRKBL = M + M*ILAENV( 1, 'ZGELQF', ' ', M, N, -1, -1 ) |
WRKBL = M + LWORK_ZGELQF_MN |
WRKBL = MAX( WRKBL, M+N*ILAENV( 1, 'ZUNGLQ', ' ', N, |
WRKBL = MAX( WRKBL, M + LWORK_ZUNGLQ_NN ) |
$ N, M, -1 ) ) |
WRKBL = MAX( WRKBL, 2*M + LWORK_ZGEBRD_MM ) |
WRKBL = MAX( WRKBL, 2*M+2*M* |
WRKBL = MAX( WRKBL, 2*M + LWORK_ZUNMBR_QLN_MM ) |
$ ILAENV( 1, 'ZGEBRD', ' ', M, M, -1, -1 ) ) |
WRKBL = MAX( WRKBL, 2*M + LWORK_ZUNMBR_PRC_MM ) |
WRKBL = MAX( WRKBL, 2*M+M* |
|
$ ILAENV( 1, 'ZUNMBR', 'PRC', M, M, M, -1 ) ) |
|
WRKBL = MAX( WRKBL, 2*M+M* |
|
$ ILAENV( 1, 'ZUNMBR', 'QLN', M, M, M, -1 ) ) |
|
MAXWRK = M*M + WRKBL |
MAXWRK = M*M + WRKBL |
MINWRK = M*M + 2*M + N |
MINWRK = M*M + MAX( 3*M, M + N ) |
END IF |
END IF |
ELSE IF( N.GE.MNTHR2 ) THEN |
ELSE IF( N.GE.MNTHR2 ) THEN |
* |
* |
* Path 5t (N much larger than M, but not as much as MNTHR1) |
* Path 5t (N >> M, but not as much as MNTHR1) |
* |
* |
MAXWRK = 2*M + ( M+N )*ILAENV( 1, 'ZGEBRD', ' ', M, N, |
MAXWRK = 2*M + LWORK_ZGEBRD_MN |
$ -1, -1 ) |
|
MINWRK = 2*M + N |
MINWRK = 2*M + N |
IF( WNTQO ) THEN |
IF( WNTQO ) THEN |
MAXWRK = MAX( MAXWRK, 2*M+M* |
* Path 5to (N >> M, JOBZ='O') |
$ ILAENV( 1, 'ZUNGBR', 'P', M, N, M, -1 ) ) |
MAXWRK = MAX( MAXWRK, 2*M + LWORK_ZUNGBR_Q_MM ) |
MAXWRK = MAX( MAXWRK, 2*M+M* |
MAXWRK = MAX( MAXWRK, 2*M + LWORK_ZUNGBR_P_MN ) |
$ ILAENV( 1, 'ZUNGBR', 'Q', M, M, N, -1 ) ) |
|
MAXWRK = MAXWRK + M*N |
MAXWRK = MAXWRK + M*N |
MINWRK = MINWRK + M*M |
MINWRK = MINWRK + M*M |
ELSE IF( WNTQS ) THEN |
ELSE IF( WNTQS ) THEN |
MAXWRK = MAX( MAXWRK, 2*M+M* |
* Path 5ts (N >> M, JOBZ='S') |
$ ILAENV( 1, 'ZUNGBR', 'P', M, N, M, -1 ) ) |
MAXWRK = MAX( MAXWRK, 2*M + LWORK_ZUNGBR_Q_MM ) |
MAXWRK = MAX( MAXWRK, 2*M+M* |
MAXWRK = MAX( MAXWRK, 2*M + LWORK_ZUNGBR_P_MN ) |
$ ILAENV( 1, 'ZUNGBR', 'Q', M, M, N, -1 ) ) |
|
ELSE IF( WNTQA ) THEN |
ELSE IF( WNTQA ) THEN |
MAXWRK = MAX( MAXWRK, 2*M+N* |
* Path 5ta (N >> M, JOBZ='A') |
$ ILAENV( 1, 'ZUNGBR', 'P', N, N, M, -1 ) ) |
MAXWRK = MAX( MAXWRK, 2*M + LWORK_ZUNGBR_Q_MM ) |
MAXWRK = MAX( MAXWRK, 2*M+M* |
MAXWRK = MAX( MAXWRK, 2*M + LWORK_ZUNGBR_P_NN ) |
$ ILAENV( 1, 'ZUNGBR', 'Q', M, M, N, -1 ) ) |
|
END IF |
END IF |
ELSE |
ELSE |
* |
* |
* Path 6t (N greater than M, but not much larger) |
* Path 6t (N > M, but not much larger) |
* |
* |
MAXWRK = 2*M + ( M+N )*ILAENV( 1, 'ZGEBRD', ' ', M, N, |
MAXWRK = 2*M + LWORK_ZGEBRD_MN |
$ -1, -1 ) |
|
MINWRK = 2*M + N |
MINWRK = 2*M + N |
IF( WNTQO ) THEN |
IF( WNTQO ) THEN |
MAXWRK = MAX( MAXWRK, 2*M+M* |
* Path 6to (N > M, JOBZ='O') |
$ ILAENV( 1, 'ZUNMBR', 'PRC', M, N, M, -1 ) ) |
MAXWRK = MAX( MAXWRK, 2*M + LWORK_ZUNMBR_QLN_MM ) |
MAXWRK = MAX( MAXWRK, 2*M+M* |
MAXWRK = MAX( MAXWRK, 2*M + LWORK_ZUNMBR_PRC_MN ) |
$ ILAENV( 1, 'ZUNMBR', 'QLN', M, M, N, -1 ) ) |
|
MAXWRK = MAXWRK + M*N |
MAXWRK = MAXWRK + M*N |
MINWRK = MINWRK + M*M |
MINWRK = MINWRK + M*M |
ELSE IF( WNTQS ) THEN |
ELSE IF( WNTQS ) THEN |
MAXWRK = MAX( MAXWRK, 2*M+M* |
* Path 6ts (N > M, JOBZ='S') |
$ ILAENV( 1, 'ZUNGBR', 'PRC', M, N, M, -1 ) ) |
MAXWRK = MAX( MAXWRK, 2*M + LWORK_ZUNMBR_QLN_MM ) |
MAXWRK = MAX( MAXWRK, 2*M+M* |
MAXWRK = MAX( MAXWRK, 2*M + LWORK_ZUNMBR_PRC_MN ) |
$ ILAENV( 1, 'ZUNGBR', 'QLN', M, M, N, -1 ) ) |
|
ELSE IF( WNTQA ) THEN |
ELSE IF( WNTQA ) THEN |
MAXWRK = MAX( MAXWRK, 2*M+N* |
* Path 6ta (N > M, JOBZ='A') |
$ ILAENV( 1, 'ZUNGBR', 'PRC', N, N, M, -1 ) ) |
MAXWRK = MAX( MAXWRK, 2*M + LWORK_ZUNMBR_QLN_MM ) |
MAXWRK = MAX( MAXWRK, 2*M+M* |
MAXWRK = MAX( MAXWRK, 2*M + LWORK_ZUNMBR_PRC_NN ) |
$ ILAENV( 1, 'ZUNGBR', 'QLN', M, M, N, -1 ) ) |
|
END IF |
END IF |
END IF |
END IF |
END IF |
END IF |
Line 554
|
Line 620
|
END IF |
END IF |
IF( INFO.EQ.0 ) THEN |
IF( INFO.EQ.0 ) THEN |
WORK( 1 ) = MAXWRK |
WORK( 1 ) = MAXWRK |
IF( LWORK.LT.MINWRK .AND. LWORK.NE.LQUERV ) |
IF( LWORK.LT.MINWRK .AND. .NOT. LQUERY ) THEN |
$ INFO = -13 |
INFO = -12 |
|
END IF |
END IF |
END IF |
* |
* |
* Quick returns |
|
* |
|
IF( INFO.NE.0 ) THEN |
IF( INFO.NE.0 ) THEN |
CALL XERBLA( 'ZGESDD', -INFO ) |
CALL XERBLA( 'ZGESDD', -INFO ) |
RETURN |
RETURN |
|
ELSE IF( LQUERY ) THEN |
|
RETURN |
END IF |
END IF |
IF( LWORK.EQ.LQUERV ) |
* |
$ RETURN |
* Quick return if possible |
|
* |
IF( M.EQ.0 .OR. N.EQ.0 ) THEN |
IF( M.EQ.0 .OR. N.EQ.0 ) THEN |
RETURN |
RETURN |
END IF |
END IF |
Line 598
|
Line 666
|
* |
* |
IF( WNTQN ) THEN |
IF( WNTQN ) THEN |
* |
* |
* Path 1 (M much larger than N, JOBZ='N') |
* Path 1 (M >> N, JOBZ='N') |
* No singular vectors to be computed |
* No singular vectors to be computed |
* |
* |
ITAU = 1 |
ITAU = 1 |
NWORK = ITAU + N |
NWORK = ITAU + N |
* |
* |
* Compute A=Q*R |
* Compute A=Q*R |
* (CWorkspace: need 2*N, prefer N+N*NB) |
* CWorkspace: need N [tau] + N [work] |
* (RWorkspace: need 0) |
* CWorkspace: prefer N [tau] + N*NB [work] |
|
* RWorkspace: need 0 |
* |
* |
CALL ZGEQRF( M, N, A, LDA, WORK( ITAU ), WORK( NWORK ), |
CALL ZGEQRF( M, N, A, LDA, WORK( ITAU ), WORK( NWORK ), |
$ LWORK-NWORK+1, IERR ) |
$ LWORK-NWORK+1, IERR ) |
Line 621
|
Line 690
|
NWORK = ITAUP + N |
NWORK = ITAUP + N |
* |
* |
* Bidiagonalize R in A |
* Bidiagonalize R in A |
* (CWorkspace: need 3*N, prefer 2*N+2*N*NB) |
* CWorkspace: need 2*N [tauq, taup] + N [work] |
* (RWorkspace: need N) |
* CWorkspace: prefer 2*N [tauq, taup] + 2*N*NB [work] |
|
* RWorkspace: need N [e] |
* |
* |
CALL ZGEBRD( N, N, A, LDA, S, RWORK( IE ), WORK( ITAUQ ), |
CALL ZGEBRD( N, N, A, LDA, S, RWORK( IE ), WORK( ITAUQ ), |
$ WORK( ITAUP ), WORK( NWORK ), LWORK-NWORK+1, |
$ WORK( ITAUP ), WORK( NWORK ), LWORK-NWORK+1, |
Line 630
|
Line 700
|
NRWORK = IE + N |
NRWORK = IE + N |
* |
* |
* Perform bidiagonal SVD, compute singular values only |
* Perform bidiagonal SVD, compute singular values only |
* (CWorkspace: 0) |
* CWorkspace: need 0 |
* (RWorkspace: need BDSPAN) |
* RWorkspace: need N [e] + BDSPAC |
* |
* |
CALL DBDSDC( 'U', 'N', N, S, RWORK( IE ), DUM, 1, DUM, 1, |
CALL DBDSDC( 'U', 'N', N, S, RWORK( IE ), DUM,1,DUM,1, |
$ DUM, IDUM, RWORK( NRWORK ), IWORK, INFO ) |
$ DUM, IDUM, RWORK( NRWORK ), IWORK, INFO ) |
* |
* |
ELSE IF( WNTQO ) THEN |
ELSE IF( WNTQO ) THEN |
* |
* |
* Path 2 (M much larger than N, JOBZ='O') |
* Path 2 (M >> N, JOBZ='O') |
* N left singular vectors to be overwritten on A and |
* N left singular vectors to be overwritten on A and |
* N right singular vectors to be computed in VT |
* N right singular vectors to be computed in VT |
* |
* |
Line 648
|
Line 718
|
* |
* |
LDWRKU = N |
LDWRKU = N |
IR = IU + LDWRKU*N |
IR = IU + LDWRKU*N |
IF( LWORK.GE.M*N+N*N+3*N ) THEN |
IF( LWORK .GE. M*N + N*N + 3*N ) THEN |
* |
* |
* WORK(IR) is M by N |
* WORK(IR) is M by N |
* |
* |
LDWRKR = M |
LDWRKR = M |
ELSE |
ELSE |
LDWRKR = ( LWORK-N*N-3*N ) / N |
LDWRKR = ( LWORK - N*N - 3*N ) / N |
END IF |
END IF |
ITAU = IR + LDWRKR*N |
ITAU = IR + LDWRKR*N |
NWORK = ITAU + N |
NWORK = ITAU + N |
* |
* |
* Compute A=Q*R |
* Compute A=Q*R |
* (CWorkspace: need N*N+2*N, prefer M*N+N+N*NB) |
* CWorkspace: need N*N [U] + N*N [R] + N [tau] + N [work] |
* (RWorkspace: 0) |
* CWorkspace: prefer N*N [U] + N*N [R] + N [tau] + N*NB [work] |
|
* RWorkspace: need 0 |
* |
* |
CALL ZGEQRF( M, N, A, LDA, WORK( ITAU ), WORK( NWORK ), |
CALL ZGEQRF( M, N, A, LDA, WORK( ITAU ), WORK( NWORK ), |
$ LWORK-NWORK+1, IERR ) |
$ LWORK-NWORK+1, IERR ) |
Line 673
|
Line 744
|
$ LDWRKR ) |
$ LDWRKR ) |
* |
* |
* Generate Q in A |
* Generate Q in A |
* (CWorkspace: need 2*N, prefer N+N*NB) |
* CWorkspace: need N*N [U] + N*N [R] + N [tau] + N [work] |
* (RWorkspace: 0) |
* CWorkspace: prefer N*N [U] + N*N [R] + N [tau] + N*NB [work] |
|
* RWorkspace: need 0 |
* |
* |
CALL ZUNGQR( M, N, N, A, LDA, WORK( ITAU ), |
CALL ZUNGQR( M, N, N, A, LDA, WORK( ITAU ), |
$ WORK( NWORK ), LWORK-NWORK+1, IERR ) |
$ WORK( NWORK ), LWORK-NWORK+1, IERR ) |
Line 684
|
Line 756
|
NWORK = ITAUP + N |
NWORK = ITAUP + N |
* |
* |
* Bidiagonalize R in WORK(IR) |
* Bidiagonalize R in WORK(IR) |
* (CWorkspace: need N*N+3*N, prefer M*N+2*N+2*N*NB) |
* CWorkspace: need N*N [U] + N*N [R] + 2*N [tauq, taup] + N [work] |
* (RWorkspace: need N) |
* CWorkspace: prefer N*N [U] + N*N [R] + 2*N [tauq, taup] + 2*N*NB [work] |
|
* RWorkspace: need N [e] |
* |
* |
CALL ZGEBRD( N, N, WORK( IR ), LDWRKR, S, RWORK( IE ), |
CALL ZGEBRD( N, N, WORK( IR ), LDWRKR, S, RWORK( IE ), |
$ WORK( ITAUQ ), WORK( ITAUP ), WORK( NWORK ), |
$ WORK( ITAUQ ), WORK( ITAUP ), WORK( NWORK ), |
Line 694
|
Line 767
|
* Perform bidiagonal SVD, computing left singular vectors |
* Perform bidiagonal SVD, computing left singular vectors |
* of R in WORK(IRU) and computing right singular vectors |
* of R in WORK(IRU) and computing right singular vectors |
* of R in WORK(IRVT) |
* of R in WORK(IRVT) |
* (CWorkspace: need 0) |
* CWorkspace: need 0 |
* (RWorkspace: need BDSPAC) |
* RWorkspace: need N [e] + N*N [RU] + N*N [RVT] + BDSPAC |
* |
* |
IRU = IE + N |
IRU = IE + N |
IRVT = IRU + N*N |
IRVT = IRU + N*N |
Line 706
|
Line 779
|
* |
* |
* Copy real matrix RWORK(IRU) to complex matrix WORK(IU) |
* Copy real matrix RWORK(IRU) to complex matrix WORK(IU) |
* Overwrite WORK(IU) by the left singular vectors of R |
* Overwrite WORK(IU) by the left singular vectors of R |
* (CWorkspace: need 2*N*N+3*N, prefer M*N+N*N+2*N+N*NB) |
* CWorkspace: need N*N [U] + N*N [R] + 2*N [tauq, taup] + N [work] |
* (RWorkspace: 0) |
* CWorkspace: prefer N*N [U] + N*N [R] + 2*N [tauq, taup] + N*NB [work] |
|
* RWorkspace: need 0 |
* |
* |
CALL ZLACP2( 'F', N, N, RWORK( IRU ), N, WORK( IU ), |
CALL ZLACP2( 'F', N, N, RWORK( IRU ), N, WORK( IU ), |
$ LDWRKU ) |
$ LDWRKU ) |
Line 717
|
Line 791
|
* |
* |
* Copy real matrix RWORK(IRVT) to complex matrix VT |
* Copy real matrix RWORK(IRVT) to complex matrix VT |
* Overwrite VT by the right singular vectors of R |
* Overwrite VT by the right singular vectors of R |
* (CWorkspace: need N*N+3*N, prefer M*N+2*N+N*NB) |
* CWorkspace: need N*N [U] + N*N [R] + 2*N [tauq, taup] + N [work] |
* (RWorkspace: 0) |
* CWorkspace: prefer N*N [U] + N*N [R] + 2*N [tauq, taup] + N*NB [work] |
|
* RWorkspace: need 0 |
* |
* |
CALL ZLACP2( 'F', N, N, RWORK( IRVT ), N, VT, LDVT ) |
CALL ZLACP2( 'F', N, N, RWORK( IRVT ), N, VT, LDVT ) |
CALL ZUNMBR( 'P', 'R', 'C', N, N, N, WORK( IR ), LDWRKR, |
CALL ZUNMBR( 'P', 'R', 'C', N, N, N, WORK( IR ), LDWRKR, |
Line 727
|
Line 802
|
* |
* |
* Multiply Q in A by left singular vectors of R in |
* Multiply Q in A by left singular vectors of R in |
* WORK(IU), storing result in WORK(IR) and copying to A |
* WORK(IU), storing result in WORK(IR) and copying to A |
* (CWorkspace: need 2*N*N, prefer N*N+M*N) |
* CWorkspace: need N*N [U] + N*N [R] |
* (RWorkspace: 0) |
* CWorkspace: prefer N*N [U] + M*N [R] |
|
* RWorkspace: need 0 |
* |
* |
DO 10 I = 1, M, LDWRKR |
DO 10 I = 1, M, LDWRKR |
CHUNK = MIN( M-I+1, LDWRKR ) |
CHUNK = MIN( M-I+1, LDWRKR ) |
Line 741
|
Line 817
|
* |
* |
ELSE IF( WNTQS ) THEN |
ELSE IF( WNTQS ) THEN |
* |
* |
* Path 3 (M much larger than N, JOBZ='S') |
* Path 3 (M >> N, JOBZ='S') |
* N left singular vectors to be computed in U and |
* N left singular vectors to be computed in U and |
* N right singular vectors to be computed in VT |
* N right singular vectors to be computed in VT |
* |
* |
Line 754
|
Line 830
|
NWORK = ITAU + N |
NWORK = ITAU + N |
* |
* |
* Compute A=Q*R |
* Compute A=Q*R |
* (CWorkspace: need N*N+2*N, prefer N*N+N+N*NB) |
* CWorkspace: need N*N [R] + N [tau] + N [work] |
* (RWorkspace: 0) |
* CWorkspace: prefer N*N [R] + N [tau] + N*NB [work] |
|
* RWorkspace: need 0 |
* |
* |
CALL ZGEQRF( M, N, A, LDA, WORK( ITAU ), WORK( NWORK ), |
CALL ZGEQRF( M, N, A, LDA, WORK( ITAU ), WORK( NWORK ), |
$ LWORK-NWORK+1, IERR ) |
$ LWORK-NWORK+1, IERR ) |
Line 767
|
Line 844
|
$ LDWRKR ) |
$ LDWRKR ) |
* |
* |
* Generate Q in A |
* Generate Q in A |
* (CWorkspace: need 2*N, prefer N+N*NB) |
* CWorkspace: need N*N [R] + N [tau] + N [work] |
* (RWorkspace: 0) |
* CWorkspace: prefer N*N [R] + N [tau] + N*NB [work] |
|
* RWorkspace: need 0 |
* |
* |
CALL ZUNGQR( M, N, N, A, LDA, WORK( ITAU ), |
CALL ZUNGQR( M, N, N, A, LDA, WORK( ITAU ), |
$ WORK( NWORK ), LWORK-NWORK+1, IERR ) |
$ WORK( NWORK ), LWORK-NWORK+1, IERR ) |
Line 778
|
Line 856
|
NWORK = ITAUP + N |
NWORK = ITAUP + N |
* |
* |
* Bidiagonalize R in WORK(IR) |
* Bidiagonalize R in WORK(IR) |
* (CWorkspace: need N*N+3*N, prefer N*N+2*N+2*N*NB) |
* CWorkspace: need N*N [R] + 2*N [tauq, taup] + N [work] |
* (RWorkspace: need N) |
* CWorkspace: prefer N*N [R] + 2*N [tauq, taup] + 2*N*NB [work] |
|
* RWorkspace: need N [e] |
* |
* |
CALL ZGEBRD( N, N, WORK( IR ), LDWRKR, S, RWORK( IE ), |
CALL ZGEBRD( N, N, WORK( IR ), LDWRKR, S, RWORK( IE ), |
$ WORK( ITAUQ ), WORK( ITAUP ), WORK( NWORK ), |
$ WORK( ITAUQ ), WORK( ITAUP ), WORK( NWORK ), |
Line 788
|
Line 867
|
* Perform bidiagonal SVD, computing left singular vectors |
* Perform bidiagonal SVD, computing left singular vectors |
* of bidiagonal matrix in RWORK(IRU) and computing right |
* of bidiagonal matrix in RWORK(IRU) and computing right |
* singular vectors of bidiagonal matrix in RWORK(IRVT) |
* singular vectors of bidiagonal matrix in RWORK(IRVT) |
* (CWorkspace: need 0) |
* CWorkspace: need 0 |
* (RWorkspace: need BDSPAC) |
* RWorkspace: need N [e] + N*N [RU] + N*N [RVT] + BDSPAC |
* |
* |
IRU = IE + N |
IRU = IE + N |
IRVT = IRU + N*N |
IRVT = IRU + N*N |
Line 800
|
Line 879
|
* |
* |
* Copy real matrix RWORK(IRU) to complex matrix U |
* Copy real matrix RWORK(IRU) to complex matrix U |
* Overwrite U by left singular vectors of R |
* Overwrite U by left singular vectors of R |
* (CWorkspace: need N*N+3*N, prefer N*N+2*N+N*NB) |
* CWorkspace: need N*N [R] + 2*N [tauq, taup] + N [work] |
* (RWorkspace: 0) |
* CWorkspace: prefer N*N [R] + 2*N [tauq, taup] + N*NB [work] |
|
* RWorkspace: need 0 |
* |
* |
CALL ZLACP2( 'F', N, N, RWORK( IRU ), N, U, LDU ) |
CALL ZLACP2( 'F', N, N, RWORK( IRU ), N, U, LDU ) |
CALL ZUNMBR( 'Q', 'L', 'N', N, N, N, WORK( IR ), LDWRKR, |
CALL ZUNMBR( 'Q', 'L', 'N', N, N, N, WORK( IR ), LDWRKR, |
Line 810
|
Line 890
|
* |
* |
* Copy real matrix RWORK(IRVT) to complex matrix VT |
* Copy real matrix RWORK(IRVT) to complex matrix VT |
* Overwrite VT by right singular vectors of R |
* Overwrite VT by right singular vectors of R |
* (CWorkspace: need N*N+3*N, prefer N*N+2*N+N*NB) |
* CWorkspace: need N*N [R] + 2*N [tauq, taup] + N [work] |
* (RWorkspace: 0) |
* CWorkspace: prefer N*N [R] + 2*N [tauq, taup] + N*NB [work] |
|
* RWorkspace: need 0 |
* |
* |
CALL ZLACP2( 'F', N, N, RWORK( IRVT ), N, VT, LDVT ) |
CALL ZLACP2( 'F', N, N, RWORK( IRVT ), N, VT, LDVT ) |
CALL ZUNMBR( 'P', 'R', 'C', N, N, N, WORK( IR ), LDWRKR, |
CALL ZUNMBR( 'P', 'R', 'C', N, N, N, WORK( IR ), LDWRKR, |
Line 820
|
Line 901
|
* |
* |
* Multiply Q in A by left singular vectors of R in |
* Multiply Q in A by left singular vectors of R in |
* WORK(IR), storing result in U |
* WORK(IR), storing result in U |
* (CWorkspace: need N*N) |
* CWorkspace: need N*N [R] |
* (RWorkspace: 0) |
* RWorkspace: need 0 |
* |
* |
CALL ZLACPY( 'F', N, N, U, LDU, WORK( IR ), LDWRKR ) |
CALL ZLACPY( 'F', N, N, U, LDU, WORK( IR ), LDWRKR ) |
CALL ZGEMM( 'N', 'N', M, N, N, CONE, A, LDA, WORK( IR ), |
CALL ZGEMM( 'N', 'N', M, N, N, CONE, A, LDA, WORK( IR ), |
Line 829
|
Line 910
|
* |
* |
ELSE IF( WNTQA ) THEN |
ELSE IF( WNTQA ) THEN |
* |
* |
* Path 4 (M much larger than N, JOBZ='A') |
* Path 4 (M >> N, JOBZ='A') |
* M left singular vectors to be computed in U and |
* M left singular vectors to be computed in U and |
* N right singular vectors to be computed in VT |
* N right singular vectors to be computed in VT |
* |
* |
Line 842
|
Line 923
|
NWORK = ITAU + N |
NWORK = ITAU + N |
* |
* |
* Compute A=Q*R, copying result to U |
* Compute A=Q*R, copying result to U |
* (CWorkspace: need 2*N, prefer N+N*NB) |
* CWorkspace: need N*N [U] + N [tau] + N [work] |
* (RWorkspace: 0) |
* CWorkspace: prefer N*N [U] + N [tau] + N*NB [work] |
|
* RWorkspace: need 0 |
* |
* |
CALL ZGEQRF( M, N, A, LDA, WORK( ITAU ), WORK( NWORK ), |
CALL ZGEQRF( M, N, A, LDA, WORK( ITAU ), WORK( NWORK ), |
$ LWORK-NWORK+1, IERR ) |
$ LWORK-NWORK+1, IERR ) |
CALL ZLACPY( 'L', M, N, A, LDA, U, LDU ) |
CALL ZLACPY( 'L', M, N, A, LDA, U, LDU ) |
* |
* |
* Generate Q in U |
* Generate Q in U |
* (CWorkspace: need N+M, prefer N+M*NB) |
* CWorkspace: need N*N [U] + N [tau] + M [work] |
* (RWorkspace: 0) |
* CWorkspace: prefer N*N [U] + N [tau] + M*NB [work] |
|
* RWorkspace: need 0 |
* |
* |
CALL ZUNGQR( M, M, N, U, LDU, WORK( ITAU ), |
CALL ZUNGQR( M, M, N, U, LDU, WORK( ITAU ), |
$ WORK( NWORK ), LWORK-NWORK+1, IERR ) |
$ WORK( NWORK ), LWORK-NWORK+1, IERR ) |
Line 866
|
Line 949
|
NWORK = ITAUP + N |
NWORK = ITAUP + N |
* |
* |
* Bidiagonalize R in A |
* Bidiagonalize R in A |
* (CWorkspace: need 3*N, prefer 2*N+2*N*NB) |
* CWorkspace: need N*N [U] + 2*N [tauq, taup] + N [work] |
* (RWorkspace: need N) |
* CWorkspace: prefer N*N [U] + 2*N [tauq, taup] + 2*N*NB [work] |
|
* RWorkspace: need N [e] |
* |
* |
CALL ZGEBRD( N, N, A, LDA, S, RWORK( IE ), WORK( ITAUQ ), |
CALL ZGEBRD( N, N, A, LDA, S, RWORK( IE ), WORK( ITAUQ ), |
$ WORK( ITAUP ), WORK( NWORK ), LWORK-NWORK+1, |
$ WORK( ITAUP ), WORK( NWORK ), LWORK-NWORK+1, |
Line 879
|
Line 963
|
* Perform bidiagonal SVD, computing left singular vectors |
* Perform bidiagonal SVD, computing left singular vectors |
* of bidiagonal matrix in RWORK(IRU) and computing right |
* of bidiagonal matrix in RWORK(IRU) and computing right |
* singular vectors of bidiagonal matrix in RWORK(IRVT) |
* singular vectors of bidiagonal matrix in RWORK(IRVT) |
* (CWorkspace: need 0) |
* CWorkspace: need 0 |
* (RWorkspace: need BDSPAC) |
* RWorkspace: need N [e] + N*N [RU] + N*N [RVT] + BDSPAC |
* |
* |
CALL DBDSDC( 'U', 'I', N, S, RWORK( IE ), RWORK( IRU ), |
CALL DBDSDC( 'U', 'I', N, S, RWORK( IE ), RWORK( IRU ), |
$ N, RWORK( IRVT ), N, DUM, IDUM, |
$ N, RWORK( IRVT ), N, DUM, IDUM, |
Line 888
|
Line 972
|
* |
* |
* Copy real matrix RWORK(IRU) to complex matrix WORK(IU) |
* Copy real matrix RWORK(IRU) to complex matrix WORK(IU) |
* Overwrite WORK(IU) by left singular vectors of R |
* Overwrite WORK(IU) by left singular vectors of R |
* (CWorkspace: need N*N+3*N, prefer N*N+2*N+N*NB) |
* CWorkspace: need N*N [U] + 2*N [tauq, taup] + N [work] |
* (RWorkspace: 0) |
* CWorkspace: prefer N*N [U] + 2*N [tauq, taup] + N*NB [work] |
|
* RWorkspace: need 0 |
* |
* |
CALL ZLACP2( 'F', N, N, RWORK( IRU ), N, WORK( IU ), |
CALL ZLACP2( 'F', N, N, RWORK( IRU ), N, WORK( IU ), |
$ LDWRKU ) |
$ LDWRKU ) |
Line 899
|
Line 984
|
* |
* |
* Copy real matrix RWORK(IRVT) to complex matrix VT |
* Copy real matrix RWORK(IRVT) to complex matrix VT |
* Overwrite VT by right singular vectors of R |
* Overwrite VT by right singular vectors of R |
* (CWorkspace: need 3*N, prefer 2*N+N*NB) |
* CWorkspace: need N*N [U] + 2*N [tauq, taup] + N [work] |
* (RWorkspace: 0) |
* CWorkspace: prefer N*N [U] + 2*N [tauq, taup] + N*NB [work] |
|
* RWorkspace: need 0 |
* |
* |
CALL ZLACP2( 'F', N, N, RWORK( IRVT ), N, VT, LDVT ) |
CALL ZLACP2( 'F', N, N, RWORK( IRVT ), N, VT, LDVT ) |
CALL ZUNMBR( 'P', 'R', 'C', N, N, N, A, LDA, |
CALL ZUNMBR( 'P', 'R', 'C', N, N, N, A, LDA, |
Line 909
|
Line 995
|
* |
* |
* Multiply Q in U by left singular vectors of R in |
* Multiply Q in U by left singular vectors of R in |
* WORK(IU), storing result in A |
* WORK(IU), storing result in A |
* (CWorkspace: need N*N) |
* CWorkspace: need N*N [U] |
* (RWorkspace: 0) |
* RWorkspace: need 0 |
* |
* |
CALL ZGEMM( 'N', 'N', M, N, N, CONE, U, LDU, WORK( IU ), |
CALL ZGEMM( 'N', 'N', M, N, N, CONE, U, LDU, WORK( IU ), |
$ LDWRKU, CZERO, A, LDA ) |
$ LDWRKU, CZERO, A, LDA ) |
Line 925
|
Line 1011
|
* |
* |
* MNTHR2 <= M < MNTHR1 |
* MNTHR2 <= M < MNTHR1 |
* |
* |
* Path 5 (M much larger than N, but not as much as MNTHR1) |
* Path 5 (M >> N, but not as much as MNTHR1) |
* Reduce to bidiagonal form without QR decomposition, use |
* Reduce to bidiagonal form without QR decomposition, use |
* ZUNGBR and matrix multiplication to compute singular vectors |
* ZUNGBR and matrix multiplication to compute singular vectors |
* |
* |
Line 936
|
Line 1022
|
NWORK = ITAUP + N |
NWORK = ITAUP + N |
* |
* |
* Bidiagonalize A |
* Bidiagonalize A |
* (CWorkspace: need 2*N+M, prefer 2*N+(M+N)*NB) |
* CWorkspace: need 2*N [tauq, taup] + M [work] |
* (RWorkspace: need N) |
* CWorkspace: prefer 2*N [tauq, taup] + (M+N)*NB [work] |
|
* RWorkspace: need N [e] |
* |
* |
CALL ZGEBRD( M, N, A, LDA, S, RWORK( IE ), WORK( ITAUQ ), |
CALL ZGEBRD( M, N, A, LDA, S, RWORK( IE ), WORK( ITAUQ ), |
$ WORK( ITAUP ), WORK( NWORK ), LWORK-NWORK+1, |
$ WORK( ITAUP ), WORK( NWORK ), LWORK-NWORK+1, |
$ IERR ) |
$ IERR ) |
IF( WNTQN ) THEN |
IF( WNTQN ) THEN |
* |
* |
|
* Path 5n (M >> N, JOBZ='N') |
* Compute singular values only |
* Compute singular values only |
* (Cworkspace: 0) |
* CWorkspace: need 0 |
* (Rworkspace: need BDSPAN) |
* RWorkspace: need N [e] + BDSPAC |
* |
* |
CALL DBDSDC( 'U', 'N', N, S, RWORK( IE ), DUM, 1, DUM, 1, |
CALL DBDSDC( 'U', 'N', N, S, RWORK( IE ), DUM, 1,DUM,1, |
$ DUM, IDUM, RWORK( NRWORK ), IWORK, INFO ) |
$ DUM, IDUM, RWORK( NRWORK ), IWORK, INFO ) |
ELSE IF( WNTQO ) THEN |
ELSE IF( WNTQO ) THEN |
IU = NWORK |
IU = NWORK |
Line 956
|
Line 1044
|
IRVT = IRU + N*N |
IRVT = IRU + N*N |
NRWORK = IRVT + N*N |
NRWORK = IRVT + N*N |
* |
* |
|
* Path 5o (M >> N, JOBZ='O') |
* Copy A to VT, generate P**H |
* Copy A to VT, generate P**H |
* (Cworkspace: need 2*N, prefer N+N*NB) |
* CWorkspace: need 2*N [tauq, taup] + N [work] |
* (Rworkspace: 0) |
* CWorkspace: prefer 2*N [tauq, taup] + N*NB [work] |
|
* RWorkspace: need 0 |
* |
* |
CALL ZLACPY( 'U', N, N, A, LDA, VT, LDVT ) |
CALL ZLACPY( 'U', N, N, A, LDA, VT, LDVT ) |
CALL ZUNGBR( 'P', N, N, N, VT, LDVT, WORK( ITAUP ), |
CALL ZUNGBR( 'P', N, N, N, VT, LDVT, WORK( ITAUP ), |
$ WORK( NWORK ), LWORK-NWORK+1, IERR ) |
$ WORK( NWORK ), LWORK-NWORK+1, IERR ) |
* |
* |
* Generate Q in A |
* Generate Q in A |
* (CWorkspace: need 2*N, prefer N+N*NB) |
* CWorkspace: need 2*N [tauq, taup] + N [work] |
* (RWorkspace: 0) |
* CWorkspace: prefer 2*N [tauq, taup] + N*NB [work] |
|
* RWorkspace: need 0 |
* |
* |
CALL ZUNGBR( 'Q', M, N, N, A, LDA, WORK( ITAUQ ), |
CALL ZUNGBR( 'Q', M, N, N, A, LDA, WORK( ITAUQ ), |
$ WORK( NWORK ), LWORK-NWORK+1, IERR ) |
$ WORK( NWORK ), LWORK-NWORK+1, IERR ) |
* |
* |
IF( LWORK.GE.M*N+3*N ) THEN |
IF( LWORK .GE. M*N + 3*N ) THEN |
* |
* |
* WORK( IU ) is M by N |
* WORK( IU ) is M by N |
* |
* |
Line 980
|
Line 1071
|
* |
* |
* WORK(IU) is LDWRKU by N |
* WORK(IU) is LDWRKU by N |
* |
* |
LDWRKU = ( LWORK-3*N ) / N |
LDWRKU = ( LWORK - 3*N ) / N |
END IF |
END IF |
NWORK = IU + LDWRKU*N |
NWORK = IU + LDWRKU*N |
* |
* |
* Perform bidiagonal SVD, computing left singular vectors |
* Perform bidiagonal SVD, computing left singular vectors |
* of bidiagonal matrix in RWORK(IRU) and computing right |
* of bidiagonal matrix in RWORK(IRU) and computing right |
* singular vectors of bidiagonal matrix in RWORK(IRVT) |
* singular vectors of bidiagonal matrix in RWORK(IRVT) |
* (CWorkspace: need 0) |
* CWorkspace: need 0 |
* (RWorkspace: need BDSPAC) |
* RWorkspace: need N [e] + N*N [RU] + N*N [RVT] + BDSPAC |
* |
* |
CALL DBDSDC( 'U', 'I', N, S, RWORK( IE ), RWORK( IRU ), |
CALL DBDSDC( 'U', 'I', N, S, RWORK( IE ), RWORK( IRU ), |
$ N, RWORK( IRVT ), N, DUM, IDUM, |
$ N, RWORK( IRVT ), N, DUM, IDUM, |
Line 996
|
Line 1087
|
* |
* |
* Multiply real matrix RWORK(IRVT) by P**H in VT, |
* Multiply real matrix RWORK(IRVT) by P**H in VT, |
* storing the result in WORK(IU), copying to VT |
* storing the result in WORK(IU), copying to VT |
* (Cworkspace: need 0) |
* CWorkspace: need 2*N [tauq, taup] + N*N [U] |
* (Rworkspace: need 3*N*N) |
* RWorkspace: need N [e] + N*N [RU] + N*N [RVT] + 2*N*N [rwork] |
* |
* |
CALL ZLARCM( N, N, RWORK( IRVT ), N, VT, LDVT, |
CALL ZLARCM( N, N, RWORK( IRVT ), N, VT, LDVT, |
$ WORK( IU ), LDWRKU, RWORK( NRWORK ) ) |
$ WORK( IU ), LDWRKU, RWORK( NRWORK ) ) |
Line 1005
|
Line 1096
|
* |
* |
* Multiply Q in A by real matrix RWORK(IRU), storing the |
* Multiply Q in A by real matrix RWORK(IRU), storing the |
* result in WORK(IU), copying to A |
* result in WORK(IU), copying to A |
* (CWorkspace: need N*N, prefer M*N) |
* CWorkspace: need 2*N [tauq, taup] + N*N [U] |
* (Rworkspace: need 3*N*N, prefer N*N+2*M*N) |
* CWorkspace: prefer 2*N [tauq, taup] + M*N [U] |
|
* RWorkspace: need N [e] + N*N [RU] + 2*N*N [rwork] |
|
* RWorkspace: prefer N [e] + N*N [RU] + 2*M*N [rwork] < N + 5*N*N since M < 2*N here |
* |
* |
NRWORK = IRVT |
NRWORK = IRVT |
DO 20 I = 1, M, LDWRKU |
DO 20 I = 1, M, LDWRKU |
Line 1019
|
Line 1112
|
* |
* |
ELSE IF( WNTQS ) THEN |
ELSE IF( WNTQS ) THEN |
* |
* |
|
* Path 5s (M >> N, JOBZ='S') |
* Copy A to VT, generate P**H |
* Copy A to VT, generate P**H |
* (Cworkspace: need 2*N, prefer N+N*NB) |
* CWorkspace: need 2*N [tauq, taup] + N [work] |
* (Rworkspace: 0) |
* CWorkspace: prefer 2*N [tauq, taup] + N*NB [work] |
|
* RWorkspace: need 0 |
* |
* |
CALL ZLACPY( 'U', N, N, A, LDA, VT, LDVT ) |
CALL ZLACPY( 'U', N, N, A, LDA, VT, LDVT ) |
CALL ZUNGBR( 'P', N, N, N, VT, LDVT, WORK( ITAUP ), |
CALL ZUNGBR( 'P', N, N, N, VT, LDVT, WORK( ITAUP ), |
$ WORK( NWORK ), LWORK-NWORK+1, IERR ) |
$ WORK( NWORK ), LWORK-NWORK+1, IERR ) |
* |
* |
* Copy A to U, generate Q |
* Copy A to U, generate Q |
* (Cworkspace: need 2*N, prefer N+N*NB) |
* CWorkspace: need 2*N [tauq, taup] + N [work] |
* (Rworkspace: 0) |
* CWorkspace: prefer 2*N [tauq, taup] + N*NB [work] |
|
* RWorkspace: need 0 |
* |
* |
CALL ZLACPY( 'L', M, N, A, LDA, U, LDU ) |
CALL ZLACPY( 'L', M, N, A, LDA, U, LDU ) |
CALL ZUNGBR( 'Q', M, N, N, U, LDU, WORK( ITAUQ ), |
CALL ZUNGBR( 'Q', M, N, N, U, LDU, WORK( ITAUQ ), |
Line 1038
|
Line 1134
|
* Perform bidiagonal SVD, computing left singular vectors |
* Perform bidiagonal SVD, computing left singular vectors |
* of bidiagonal matrix in RWORK(IRU) and computing right |
* of bidiagonal matrix in RWORK(IRU) and computing right |
* singular vectors of bidiagonal matrix in RWORK(IRVT) |
* singular vectors of bidiagonal matrix in RWORK(IRVT) |
* (CWorkspace: need 0) |
* CWorkspace: need 0 |
* (RWorkspace: need BDSPAC) |
* RWorkspace: need N [e] + N*N [RU] + N*N [RVT] + BDSPAC |
* |
* |
IRU = NRWORK |
IRU = NRWORK |
IRVT = IRU + N*N |
IRVT = IRU + N*N |
Line 1050
|
Line 1146
|
* |
* |
* Multiply real matrix RWORK(IRVT) by P**H in VT, |
* Multiply real matrix RWORK(IRVT) by P**H in VT, |
* storing the result in A, copying to VT |
* storing the result in A, copying to VT |
* (Cworkspace: need 0) |
* CWorkspace: need 0 |
* (Rworkspace: need 3*N*N) |
* RWorkspace: need N [e] + N*N [RU] + N*N [RVT] + 2*N*N [rwork] |
* |
* |
CALL ZLARCM( N, N, RWORK( IRVT ), N, VT, LDVT, A, LDA, |
CALL ZLARCM( N, N, RWORK( IRVT ), N, VT, LDVT, A, LDA, |
$ RWORK( NRWORK ) ) |
$ RWORK( NRWORK ) ) |
Line 1059
|
Line 1155
|
* |
* |
* Multiply Q in U by real matrix RWORK(IRU), storing the |
* Multiply Q in U by real matrix RWORK(IRU), storing the |
* result in A, copying to U |
* result in A, copying to U |
* (CWorkspace: need 0) |
* CWorkspace: need 0 |
* (Rworkspace: need N*N+2*M*N) |
* RWorkspace: need N [e] + N*N [RU] + 2*M*N [rwork] < N + 5*N*N since M < 2*N here |
* |
* |
NRWORK = IRVT |
NRWORK = IRVT |
CALL ZLACRM( M, N, U, LDU, RWORK( IRU ), N, A, LDA, |
CALL ZLACRM( M, N, U, LDU, RWORK( IRU ), N, A, LDA, |
Line 1068
|
Line 1164
|
CALL ZLACPY( 'F', M, N, A, LDA, U, LDU ) |
CALL ZLACPY( 'F', M, N, A, LDA, U, LDU ) |
ELSE |
ELSE |
* |
* |
|
* Path 5a (M >> N, JOBZ='A') |
* Copy A to VT, generate P**H |
* Copy A to VT, generate P**H |
* (Cworkspace: need 2*N, prefer N+N*NB) |
* CWorkspace: need 2*N [tauq, taup] + N [work] |
* (Rworkspace: 0) |
* CWorkspace: prefer 2*N [tauq, taup] + N*NB [work] |
|
* RWorkspace: need 0 |
* |
* |
CALL ZLACPY( 'U', N, N, A, LDA, VT, LDVT ) |
CALL ZLACPY( 'U', N, N, A, LDA, VT, LDVT ) |
CALL ZUNGBR( 'P', N, N, N, VT, LDVT, WORK( ITAUP ), |
CALL ZUNGBR( 'P', N, N, N, VT, LDVT, WORK( ITAUP ), |
$ WORK( NWORK ), LWORK-NWORK+1, IERR ) |
$ WORK( NWORK ), LWORK-NWORK+1, IERR ) |
* |
* |
* Copy A to U, generate Q |
* Copy A to U, generate Q |
* (Cworkspace: need 2*N, prefer N+N*NB) |
* CWorkspace: need 2*N [tauq, taup] + M [work] |
* (Rworkspace: 0) |
* CWorkspace: prefer 2*N [tauq, taup] + M*NB [work] |
|
* RWorkspace: need 0 |
* |
* |
CALL ZLACPY( 'L', M, N, A, LDA, U, LDU ) |
CALL ZLACPY( 'L', M, N, A, LDA, U, LDU ) |
CALL ZUNGBR( 'Q', M, M, N, U, LDU, WORK( ITAUQ ), |
CALL ZUNGBR( 'Q', M, M, N, U, LDU, WORK( ITAUQ ), |
Line 1087
|
Line 1186
|
* Perform bidiagonal SVD, computing left singular vectors |
* Perform bidiagonal SVD, computing left singular vectors |
* of bidiagonal matrix in RWORK(IRU) and computing right |
* of bidiagonal matrix in RWORK(IRU) and computing right |
* singular vectors of bidiagonal matrix in RWORK(IRVT) |
* singular vectors of bidiagonal matrix in RWORK(IRVT) |
* (CWorkspace: need 0) |
* CWorkspace: need 0 |
* (RWorkspace: need BDSPAC) |
* RWorkspace: need N [e] + N*N [RU] + N*N [RVT] + BDSPAC |
* |
* |
IRU = NRWORK |
IRU = NRWORK |
IRVT = IRU + N*N |
IRVT = IRU + N*N |
Line 1099
|
Line 1198
|
* |
* |
* Multiply real matrix RWORK(IRVT) by P**H in VT, |
* Multiply real matrix RWORK(IRVT) by P**H in VT, |
* storing the result in A, copying to VT |
* storing the result in A, copying to VT |
* (Cworkspace: need 0) |
* CWorkspace: need 0 |
* (Rworkspace: need 3*N*N) |
* RWorkspace: need N [e] + N*N [RU] + N*N [RVT] + 2*N*N [rwork] |
* |
* |
CALL ZLARCM( N, N, RWORK( IRVT ), N, VT, LDVT, A, LDA, |
CALL ZLARCM( N, N, RWORK( IRVT ), N, VT, LDVT, A, LDA, |
$ RWORK( NRWORK ) ) |
$ RWORK( NRWORK ) ) |
Line 1108
|
Line 1207
|
* |
* |
* Multiply Q in U by real matrix RWORK(IRU), storing the |
* Multiply Q in U by real matrix RWORK(IRU), storing the |
* result in A, copying to U |
* result in A, copying to U |
* (CWorkspace: 0) |
* CWorkspace: need 0 |
* (Rworkspace: need 3*N*N) |
* RWorkspace: need N [e] + N*N [RU] + 2*M*N [rwork] < N + 5*N*N since M < 2*N here |
* |
* |
NRWORK = IRVT |
NRWORK = IRVT |
CALL ZLACRM( M, N, U, LDU, RWORK( IRU ), N, A, LDA, |
CALL ZLACRM( M, N, U, LDU, RWORK( IRU ), N, A, LDA, |
Line 1121
|
Line 1220
|
* |
* |
* M .LT. MNTHR2 |
* M .LT. MNTHR2 |
* |
* |
* Path 6 (M at least N, but not much larger) |
* Path 6 (M >= N, but not much larger) |
* Reduce to bidiagonal form without QR decomposition |
* Reduce to bidiagonal form without QR decomposition |
* Use ZUNMBR to compute singular vectors |
* Use ZUNMBR to compute singular vectors |
* |
* |
Line 1132
|
Line 1231
|
NWORK = ITAUP + N |
NWORK = ITAUP + N |
* |
* |
* Bidiagonalize A |
* Bidiagonalize A |
* (CWorkspace: need 2*N+M, prefer 2*N+(M+N)*NB) |
* CWorkspace: need 2*N [tauq, taup] + M [work] |
* (RWorkspace: need N) |
* CWorkspace: prefer 2*N [tauq, taup] + (M+N)*NB [work] |
|
* RWorkspace: need N [e] |
* |
* |
CALL ZGEBRD( M, N, A, LDA, S, RWORK( IE ), WORK( ITAUQ ), |
CALL ZGEBRD( M, N, A, LDA, S, RWORK( IE ), WORK( ITAUQ ), |
$ WORK( ITAUP ), WORK( NWORK ), LWORK-NWORK+1, |
$ WORK( ITAUP ), WORK( NWORK ), LWORK-NWORK+1, |
$ IERR ) |
$ IERR ) |
IF( WNTQN ) THEN |
IF( WNTQN ) THEN |
* |
* |
|
* Path 6n (M >= N, JOBZ='N') |
* Compute singular values only |
* Compute singular values only |
* (Cworkspace: 0) |
* CWorkspace: need 0 |
* (Rworkspace: need BDSPAN) |
* RWorkspace: need N [e] + BDSPAC |
* |
* |
CALL DBDSDC( 'U', 'N', N, S, RWORK( IE ), DUM, 1, DUM, 1, |
CALL DBDSDC( 'U', 'N', N, S, RWORK( IE ), DUM,1,DUM,1, |
$ DUM, IDUM, RWORK( NRWORK ), IWORK, INFO ) |
$ DUM, IDUM, RWORK( NRWORK ), IWORK, INFO ) |
ELSE IF( WNTQO ) THEN |
ELSE IF( WNTQO ) THEN |
IU = NWORK |
IU = NWORK |
IRU = NRWORK |
IRU = NRWORK |
IRVT = IRU + N*N |
IRVT = IRU + N*N |
NRWORK = IRVT + N*N |
NRWORK = IRVT + N*N |
IF( LWORK.GE.M*N+3*N ) THEN |
IF( LWORK .GE. M*N + 3*N ) THEN |
* |
* |
* WORK( IU ) is M by N |
* WORK( IU ) is M by N |
* |
* |
Line 1160
|
Line 1261
|
* |
* |
* WORK( IU ) is LDWRKU by N |
* WORK( IU ) is LDWRKU by N |
* |
* |
LDWRKU = ( LWORK-3*N ) / N |
LDWRKU = ( LWORK - 3*N ) / N |
END IF |
END IF |
NWORK = IU + LDWRKU*N |
NWORK = IU + LDWRKU*N |
* |
* |
|
* Path 6o (M >= N, JOBZ='O') |
* Perform bidiagonal SVD, computing left singular vectors |
* Perform bidiagonal SVD, computing left singular vectors |
* of bidiagonal matrix in RWORK(IRU) and computing right |
* of bidiagonal matrix in RWORK(IRU) and computing right |
* singular vectors of bidiagonal matrix in RWORK(IRVT) |
* singular vectors of bidiagonal matrix in RWORK(IRVT) |
* (CWorkspace: need 0) |
* CWorkspace: need 0 |
* (RWorkspace: need BDSPAC) |
* RWorkspace: need N [e] + N*N [RU] + N*N [RVT] + BDSPAC |
* |
* |
CALL DBDSDC( 'U', 'I', N, S, RWORK( IE ), RWORK( IRU ), |
CALL DBDSDC( 'U', 'I', N, S, RWORK( IE ), RWORK( IRU ), |
$ N, RWORK( IRVT ), N, DUM, IDUM, |
$ N, RWORK( IRVT ), N, DUM, IDUM, |
Line 1176
|
Line 1278
|
* |
* |
* Copy real matrix RWORK(IRVT) to complex matrix VT |
* Copy real matrix RWORK(IRVT) to complex matrix VT |
* Overwrite VT by right singular vectors of A |
* Overwrite VT by right singular vectors of A |
* (Cworkspace: need 2*N, prefer N+N*NB) |
* CWorkspace: need 2*N [tauq, taup] + N*N [U] + N [work] |
* (Rworkspace: need 0) |
* CWorkspace: prefer 2*N [tauq, taup] + N*N [U] + N*NB [work] |
|
* RWorkspace: need N [e] + N*N [RU] + N*N [RVT] |
* |
* |
CALL ZLACP2( 'F', N, N, RWORK( IRVT ), N, VT, LDVT ) |
CALL ZLACP2( 'F', N, N, RWORK( IRVT ), N, VT, LDVT ) |
CALL ZUNMBR( 'P', 'R', 'C', N, N, N, A, LDA, |
CALL ZUNMBR( 'P', 'R', 'C', N, N, N, A, LDA, |
$ WORK( ITAUP ), VT, LDVT, WORK( NWORK ), |
$ WORK( ITAUP ), VT, LDVT, WORK( NWORK ), |
$ LWORK-NWORK+1, IERR ) |
$ LWORK-NWORK+1, IERR ) |
* |
* |
IF( LWORK.GE.M*N+3*N ) THEN |
IF( LWORK .GE. M*N + 3*N ) THEN |
* |
* |
* Copy real matrix RWORK(IRU) to complex matrix WORK(IU) |
* Path 6o-fast |
* Overwrite WORK(IU) by left singular vectors of A, copying |
* Copy real matrix RWORK(IRU) to complex matrix WORK(IU) |
* to A |
* Overwrite WORK(IU) by left singular vectors of A, copying |
* (Cworkspace: need M*N+2*N, prefer M*N+N+N*NB) |
* to A |
* (Rworkspace: need 0) |
* CWorkspace: need 2*N [tauq, taup] + M*N [U] + N [work] |
|
* CWorkspace: prefer 2*N [tauq, taup] + M*N [U] + N*NB [work] |
|
* RWorkspace: need N [e] + N*N [RU] |
* |
* |
CALL ZLASET( 'F', M, N, CZERO, CZERO, WORK( IU ), |
CALL ZLASET( 'F', M, N, CZERO, CZERO, WORK( IU ), |
$ LDWRKU ) |
$ LDWRKU ) |
Line 1202
|
Line 1307
|
CALL ZLACPY( 'F', M, N, WORK( IU ), LDWRKU, A, LDA ) |
CALL ZLACPY( 'F', M, N, WORK( IU ), LDWRKU, A, LDA ) |
ELSE |
ELSE |
* |
* |
|
* Path 6o-slow |
* Generate Q in A |
* Generate Q in A |
* (Cworkspace: need 2*N, prefer N+N*NB) |
* CWorkspace: need 2*N [tauq, taup] + N*N [U] + N [work] |
* (Rworkspace: need 0) |
* CWorkspace: prefer 2*N [tauq, taup] + N*N [U] + N*NB [work] |
|
* RWorkspace: need 0 |
* |
* |
CALL ZUNGBR( 'Q', M, N, N, A, LDA, WORK( ITAUQ ), |
CALL ZUNGBR( 'Q', M, N, N, A, LDA, WORK( ITAUQ ), |
$ WORK( NWORK ), LWORK-NWORK+1, IERR ) |
$ WORK( NWORK ), LWORK-NWORK+1, IERR ) |
* |
* |
* Multiply Q in A by real matrix RWORK(IRU), storing the |
* Multiply Q in A by real matrix RWORK(IRU), storing the |
* result in WORK(IU), copying to A |
* result in WORK(IU), copying to A |
* (CWorkspace: need N*N, prefer M*N) |
* CWorkspace: need 2*N [tauq, taup] + N*N [U] |
* (Rworkspace: need 3*N*N, prefer N*N+2*M*N) |
* CWorkspace: prefer 2*N [tauq, taup] + M*N [U] |
|
* RWorkspace: need N [e] + N*N [RU] + 2*N*N [rwork] |
|
* RWorkspace: prefer N [e] + N*N [RU] + 2*M*N [rwork] < N + 5*N*N since M < 2*N here |
* |
* |
NRWORK = IRVT |
NRWORK = IRVT |
DO 30 I = 1, M, LDWRKU |
DO 30 I = 1, M, LDWRKU |
Line 1227
|
Line 1336
|
* |
* |
ELSE IF( WNTQS ) THEN |
ELSE IF( WNTQS ) THEN |
* |
* |
|
* Path 6s (M >= N, JOBZ='S') |
* Perform bidiagonal SVD, computing left singular vectors |
* Perform bidiagonal SVD, computing left singular vectors |
* of bidiagonal matrix in RWORK(IRU) and computing right |
* of bidiagonal matrix in RWORK(IRU) and computing right |
* singular vectors of bidiagonal matrix in RWORK(IRVT) |
* singular vectors of bidiagonal matrix in RWORK(IRVT) |
* (CWorkspace: need 0) |
* CWorkspace: need 0 |
* (RWorkspace: need BDSPAC) |
* RWorkspace: need N [e] + N*N [RU] + N*N [RVT] + BDSPAC |
* |
* |
IRU = NRWORK |
IRU = NRWORK |
IRVT = IRU + N*N |
IRVT = IRU + N*N |
Line 1242
|
Line 1352
|
* |
* |
* Copy real matrix RWORK(IRU) to complex matrix U |
* Copy real matrix RWORK(IRU) to complex matrix U |
* Overwrite U by left singular vectors of A |
* Overwrite U by left singular vectors of A |
* (CWorkspace: need 3*N, prefer 2*N+N*NB) |
* CWorkspace: need 2*N [tauq, taup] + N [work] |
* (RWorkspace: 0) |
* CWorkspace: prefer 2*N [tauq, taup] + N*NB [work] |
|
* RWorkspace: need N [e] + N*N [RU] + N*N [RVT] |
* |
* |
CALL ZLASET( 'F', M, N, CZERO, CZERO, U, LDU ) |
CALL ZLASET( 'F', M, N, CZERO, CZERO, U, LDU ) |
CALL ZLACP2( 'F', N, N, RWORK( IRU ), N, U, LDU ) |
CALL ZLACP2( 'F', N, N, RWORK( IRU ), N, U, LDU ) |
Line 1253
|
Line 1364
|
* |
* |
* Copy real matrix RWORK(IRVT) to complex matrix VT |
* Copy real matrix RWORK(IRVT) to complex matrix VT |
* Overwrite VT by right singular vectors of A |
* Overwrite VT by right singular vectors of A |
* (CWorkspace: need 3*N, prefer 2*N+N*NB) |
* CWorkspace: need 2*N [tauq, taup] + N [work] |
* (RWorkspace: 0) |
* CWorkspace: prefer 2*N [tauq, taup] + N*NB [work] |
|
* RWorkspace: need N [e] + N*N [RU] + N*N [RVT] |
* |
* |
CALL ZLACP2( 'F', N, N, RWORK( IRVT ), N, VT, LDVT ) |
CALL ZLACP2( 'F', N, N, RWORK( IRVT ), N, VT, LDVT ) |
CALL ZUNMBR( 'P', 'R', 'C', N, N, N, A, LDA, |
CALL ZUNMBR( 'P', 'R', 'C', N, N, N, A, LDA, |
Line 1262
|
Line 1374
|
$ LWORK-NWORK+1, IERR ) |
$ LWORK-NWORK+1, IERR ) |
ELSE |
ELSE |
* |
* |
|
* Path 6a (M >= N, JOBZ='A') |
* Perform bidiagonal SVD, computing left singular vectors |
* Perform bidiagonal SVD, computing left singular vectors |
* of bidiagonal matrix in RWORK(IRU) and computing right |
* of bidiagonal matrix in RWORK(IRU) and computing right |
* singular vectors of bidiagonal matrix in RWORK(IRVT) |
* singular vectors of bidiagonal matrix in RWORK(IRVT) |
* (CWorkspace: need 0) |
* CWorkspace: need 0 |
* (RWorkspace: need BDSPAC) |
* RWorkspace: need N [e] + N*N [RU] + N*N [RVT] + BDSPAC |
* |
* |
IRU = NRWORK |
IRU = NRWORK |
IRVT = IRU + N*N |
IRVT = IRU + N*N |
Line 1285
|
Line 1398
|
* |
* |
* Copy real matrix RWORK(IRU) to complex matrix U |
* Copy real matrix RWORK(IRU) to complex matrix U |
* Overwrite U by left singular vectors of A |
* Overwrite U by left singular vectors of A |
* (CWorkspace: need 2*N+M, prefer 2*N+M*NB) |
* CWorkspace: need 2*N [tauq, taup] + M [work] |
* (RWorkspace: 0) |
* CWorkspace: prefer 2*N [tauq, taup] + M*NB [work] |
|
* RWorkspace: need N [e] + N*N [RU] + N*N [RVT] |
* |
* |
CALL ZLACP2( 'F', N, N, RWORK( IRU ), N, U, LDU ) |
CALL ZLACP2( 'F', N, N, RWORK( IRU ), N, U, LDU ) |
CALL ZUNMBR( 'Q', 'L', 'N', M, M, N, A, LDA, |
CALL ZUNMBR( 'Q', 'L', 'N', M, M, N, A, LDA, |
Line 1295
|
Line 1409
|
* |
* |
* Copy real matrix RWORK(IRVT) to complex matrix VT |
* Copy real matrix RWORK(IRVT) to complex matrix VT |
* Overwrite VT by right singular vectors of A |
* Overwrite VT by right singular vectors of A |
* (CWorkspace: need 3*N, prefer 2*N+N*NB) |
* CWorkspace: need 2*N [tauq, taup] + N [work] |
* (RWorkspace: 0) |
* CWorkspace: prefer 2*N [tauq, taup] + N*NB [work] |
|
* RWorkspace: need N [e] + N*N [RU] + N*N [RVT] |
* |
* |
CALL ZLACP2( 'F', N, N, RWORK( IRVT ), N, VT, LDVT ) |
CALL ZLACP2( 'F', N, N, RWORK( IRVT ), N, VT, LDVT ) |
CALL ZUNMBR( 'P', 'R', 'C', N, N, N, A, LDA, |
CALL ZUNMBR( 'P', 'R', 'C', N, N, N, A, LDA, |
Line 1316
|
Line 1431
|
* |
* |
IF( WNTQN ) THEN |
IF( WNTQN ) THEN |
* |
* |
* Path 1t (N much larger than M, JOBZ='N') |
* Path 1t (N >> M, JOBZ='N') |
* No singular vectors to be computed |
* No singular vectors to be computed |
* |
* |
ITAU = 1 |
ITAU = 1 |
NWORK = ITAU + M |
NWORK = ITAU + M |
* |
* |
* Compute A=L*Q |
* Compute A=L*Q |
* (CWorkspace: need 2*M, prefer M+M*NB) |
* CWorkspace: need M [tau] + M [work] |
* (RWorkspace: 0) |
* CWorkspace: prefer M [tau] + M*NB [work] |
|
* RWorkspace: need 0 |
* |
* |
CALL ZGELQF( M, N, A, LDA, WORK( ITAU ), WORK( NWORK ), |
CALL ZGELQF( M, N, A, LDA, WORK( ITAU ), WORK( NWORK ), |
$ LWORK-NWORK+1, IERR ) |
$ LWORK-NWORK+1, IERR ) |
Line 1339
|
Line 1455
|
NWORK = ITAUP + M |
NWORK = ITAUP + M |
* |
* |
* Bidiagonalize L in A |
* Bidiagonalize L in A |
* (CWorkspace: need 3*M, prefer 2*M+2*M*NB) |
* CWorkspace: need 2*M [tauq, taup] + M [work] |
* (RWorkspace: need M) |
* CWorkspace: prefer 2*M [tauq, taup] + 2*M*NB [work] |
|
* RWorkspace: need M [e] |
* |
* |
CALL ZGEBRD( M, M, A, LDA, S, RWORK( IE ), WORK( ITAUQ ), |
CALL ZGEBRD( M, M, A, LDA, S, RWORK( IE ), WORK( ITAUQ ), |
$ WORK( ITAUP ), WORK( NWORK ), LWORK-NWORK+1, |
$ WORK( ITAUP ), WORK( NWORK ), LWORK-NWORK+1, |
Line 1348
|
Line 1465
|
NRWORK = IE + M |
NRWORK = IE + M |
* |
* |
* Perform bidiagonal SVD, compute singular values only |
* Perform bidiagonal SVD, compute singular values only |
* (CWorkspace: 0) |
* CWorkspace: need 0 |
* (RWorkspace: need BDSPAN) |
* RWorkspace: need M [e] + BDSPAC |
* |
* |
CALL DBDSDC( 'U', 'N', M, S, RWORK( IE ), DUM, 1, DUM, 1, |
CALL DBDSDC( 'U', 'N', M, S, RWORK( IE ), DUM,1,DUM,1, |
$ DUM, IDUM, RWORK( NRWORK ), IWORK, INFO ) |
$ DUM, IDUM, RWORK( NRWORK ), IWORK, INFO ) |
* |
* |
ELSE IF( WNTQO ) THEN |
ELSE IF( WNTQO ) THEN |
* |
* |
* Path 2t (N much larger than M, JOBZ='O') |
* Path 2t (N >> M, JOBZ='O') |
* M right singular vectors to be overwritten on A and |
* M right singular vectors to be overwritten on A and |
* M left singular vectors to be computed in U |
* M left singular vectors to be computed in U |
* |
* |
Line 1366
|
Line 1483
|
* WORK(IVT) is M by M |
* WORK(IVT) is M by M |
* |
* |
IL = IVT + LDWKVT*M |
IL = IVT + LDWKVT*M |
IF( LWORK.GE.M*N+M*M+3*M ) THEN |
IF( LWORK .GE. M*N + M*M + 3*M ) THEN |
* |
* |
* WORK(IL) M by N |
* WORK(IL) M by N |
* |
* |
Line 1377
|
Line 1494
|
* WORK(IL) is M by CHUNK |
* WORK(IL) is M by CHUNK |
* |
* |
LDWRKL = M |
LDWRKL = M |
CHUNK = ( LWORK-M*M-3*M ) / M |
CHUNK = ( LWORK - M*M - 3*M ) / M |
END IF |
END IF |
ITAU = IL + LDWRKL*CHUNK |
ITAU = IL + LDWRKL*CHUNK |
NWORK = ITAU + M |
NWORK = ITAU + M |
* |
* |
* Compute A=L*Q |
* Compute A=L*Q |
* (CWorkspace: need 2*M, prefer M+M*NB) |
* CWorkspace: need M*M [VT] + M*M [L] + M [tau] + M [work] |
* (RWorkspace: 0) |
* CWorkspace: prefer M*M [VT] + M*M [L] + M [tau] + M*NB [work] |
|
* RWorkspace: need 0 |
* |
* |
CALL ZGELQF( M, N, A, LDA, WORK( ITAU ), WORK( NWORK ), |
CALL ZGELQF( M, N, A, LDA, WORK( ITAU ), WORK( NWORK ), |
$ LWORK-NWORK+1, IERR ) |
$ LWORK-NWORK+1, IERR ) |
Line 1396
|
Line 1514
|
$ WORK( IL+LDWRKL ), LDWRKL ) |
$ WORK( IL+LDWRKL ), LDWRKL ) |
* |
* |
* Generate Q in A |
* Generate Q in A |
* (CWorkspace: need M*M+2*M, prefer M*M+M+M*NB) |
* CWorkspace: need M*M [VT] + M*M [L] + M [tau] + M [work] |
* (RWorkspace: 0) |
* CWorkspace: prefer M*M [VT] + M*M [L] + M [tau] + M*NB [work] |
|
* RWorkspace: need 0 |
* |
* |
CALL ZUNGLQ( M, N, M, A, LDA, WORK( ITAU ), |
CALL ZUNGLQ( M, N, M, A, LDA, WORK( ITAU ), |
$ WORK( NWORK ), LWORK-NWORK+1, IERR ) |
$ WORK( NWORK ), LWORK-NWORK+1, IERR ) |
Line 1407
|
Line 1526
|
NWORK = ITAUP + M |
NWORK = ITAUP + M |
* |
* |
* Bidiagonalize L in WORK(IL) |
* Bidiagonalize L in WORK(IL) |
* (CWorkspace: need M*M+3*M, prefer M*M+2*M+2*M*NB) |
* CWorkspace: need M*M [VT] + M*M [L] + 2*M [tauq, taup] + M [work] |
* (RWorkspace: need M) |
* CWorkspace: prefer M*M [VT] + M*M [L] + 2*M [tauq, taup] + 2*M*NB [work] |
|
* RWorkspace: need M [e] |
* |
* |
CALL ZGEBRD( M, M, WORK( IL ), LDWRKL, S, RWORK( IE ), |
CALL ZGEBRD( M, M, WORK( IL ), LDWRKL, S, RWORK( IE ), |
$ WORK( ITAUQ ), WORK( ITAUP ), WORK( NWORK ), |
$ WORK( ITAUQ ), WORK( ITAUP ), WORK( NWORK ), |
Line 1417
|
Line 1537
|
* Perform bidiagonal SVD, computing left singular vectors |
* Perform bidiagonal SVD, computing left singular vectors |
* of bidiagonal matrix in RWORK(IRU) and computing right |
* of bidiagonal matrix in RWORK(IRU) and computing right |
* singular vectors of bidiagonal matrix in RWORK(IRVT) |
* singular vectors of bidiagonal matrix in RWORK(IRVT) |
* (CWorkspace: need 0) |
* CWorkspace: need 0 |
* (RWorkspace: need BDSPAC) |
* RWorkspace: need M [e] + M*M [RU] + M*M [RVT] + BDSPAC |
* |
* |
IRU = IE + M |
IRU = IE + M |
IRVT = IRU + M*M |
IRVT = IRU + M*M |
Line 1429
|
Line 1549
|
* |
* |
* Copy real matrix RWORK(IRU) to complex matrix WORK(IU) |
* Copy real matrix RWORK(IRU) to complex matrix WORK(IU) |
* Overwrite WORK(IU) by the left singular vectors of L |
* Overwrite WORK(IU) by the left singular vectors of L |
* (CWorkspace: need N*N+3*N, prefer M*N+2*N+N*NB) |
* CWorkspace: need M*M [VT] + M*M [L] + 2*M [tauq, taup] + M [work] |
* (RWorkspace: 0) |
* CWorkspace: prefer M*M [VT] + M*M [L] + 2*M [tauq, taup] + M*NB [work] |
|
* RWorkspace: need 0 |
* |
* |
CALL ZLACP2( 'F', M, M, RWORK( IRU ), M, U, LDU ) |
CALL ZLACP2( 'F', M, M, RWORK( IRU ), M, U, LDU ) |
CALL ZUNMBR( 'Q', 'L', 'N', M, M, M, WORK( IL ), LDWRKL, |
CALL ZUNMBR( 'Q', 'L', 'N', M, M, M, WORK( IL ), LDWRKL, |
Line 1439
|
Line 1560
|
* |
* |
* Copy real matrix RWORK(IRVT) to complex matrix WORK(IVT) |
* Copy real matrix RWORK(IRVT) to complex matrix WORK(IVT) |
* Overwrite WORK(IVT) by the right singular vectors of L |
* Overwrite WORK(IVT) by the right singular vectors of L |
* (CWorkspace: need N*N+3*N, prefer M*N+2*N+N*NB) |
* CWorkspace: need M*M [VT] + M*M [L] + 2*M [tauq, taup] + M [work] |
* (RWorkspace: 0) |
* CWorkspace: prefer M*M [VT] + M*M [L] + 2*M [tauq, taup] + M*NB [work] |
|
* RWorkspace: need 0 |
* |
* |
CALL ZLACP2( 'F', M, M, RWORK( IRVT ), M, WORK( IVT ), |
CALL ZLACP2( 'F', M, M, RWORK( IRVT ), M, WORK( IVT ), |
$ LDWKVT ) |
$ LDWKVT ) |
Line 1450
|
Line 1572
|
* |
* |
* Multiply right singular vectors of L in WORK(IL) by Q |
* Multiply right singular vectors of L in WORK(IL) by Q |
* in A, storing result in WORK(IL) and copying to A |
* in A, storing result in WORK(IL) and copying to A |
* (CWorkspace: need 2*M*M, prefer M*M+M*N)) |
* CWorkspace: need M*M [VT] + M*M [L] |
* (RWorkspace: 0) |
* CWorkspace: prefer M*M [VT] + M*N [L] |
|
* RWorkspace: need 0 |
* |
* |
DO 40 I = 1, N, CHUNK |
DO 40 I = 1, N, CHUNK |
BLK = MIN( N-I+1, CHUNK ) |
BLK = MIN( N-I+1, CHUNK ) |
Line 1464
|
Line 1587
|
* |
* |
ELSE IF( WNTQS ) THEN |
ELSE IF( WNTQS ) THEN |
* |
* |
* Path 3t (N much larger than M, JOBZ='S') |
* Path 3t (N >> M, JOBZ='S') |
* M right singular vectors to be computed in VT and |
* M right singular vectors to be computed in VT and |
* M left singular vectors to be computed in U |
* M left singular vectors to be computed in U |
* |
* |
IL = 1 |
IL = 1 |
* |
* |
Line 1477
|
Line 1600
|
NWORK = ITAU + M |
NWORK = ITAU + M |
* |
* |
* Compute A=L*Q |
* Compute A=L*Q |
* (CWorkspace: need 2*M, prefer M+M*NB) |
* CWorkspace: need M*M [L] + M [tau] + M [work] |
* (RWorkspace: 0) |
* CWorkspace: prefer M*M [L] + M [tau] + M*NB [work] |
|
* RWorkspace: need 0 |
* |
* |
CALL ZGELQF( M, N, A, LDA, WORK( ITAU ), WORK( NWORK ), |
CALL ZGELQF( M, N, A, LDA, WORK( ITAU ), WORK( NWORK ), |
$ LWORK-NWORK+1, IERR ) |
$ LWORK-NWORK+1, IERR ) |
Line 1490
|
Line 1614
|
$ WORK( IL+LDWRKL ), LDWRKL ) |
$ WORK( IL+LDWRKL ), LDWRKL ) |
* |
* |
* Generate Q in A |
* Generate Q in A |
* (CWorkspace: need M*M+2*M, prefer M*M+M+M*NB) |
* CWorkspace: need M*M [L] + M [tau] + M [work] |
* (RWorkspace: 0) |
* CWorkspace: prefer M*M [L] + M [tau] + M*NB [work] |
|
* RWorkspace: need 0 |
* |
* |
CALL ZUNGLQ( M, N, M, A, LDA, WORK( ITAU ), |
CALL ZUNGLQ( M, N, M, A, LDA, WORK( ITAU ), |
$ WORK( NWORK ), LWORK-NWORK+1, IERR ) |
$ WORK( NWORK ), LWORK-NWORK+1, IERR ) |
Line 1501
|
Line 1626
|
NWORK = ITAUP + M |
NWORK = ITAUP + M |
* |
* |
* Bidiagonalize L in WORK(IL) |
* Bidiagonalize L in WORK(IL) |
* (CWorkspace: need M*M+3*M, prefer M*M+2*M+2*M*NB) |
* CWorkspace: need M*M [L] + 2*M [tauq, taup] + M [work] |
* (RWorkspace: need M) |
* CWorkspace: prefer M*M [L] + 2*M [tauq, taup] + 2*M*NB [work] |
|
* RWorkspace: need M [e] |
* |
* |
CALL ZGEBRD( M, M, WORK( IL ), LDWRKL, S, RWORK( IE ), |
CALL ZGEBRD( M, M, WORK( IL ), LDWRKL, S, RWORK( IE ), |
$ WORK( ITAUQ ), WORK( ITAUP ), WORK( NWORK ), |
$ WORK( ITAUQ ), WORK( ITAUP ), WORK( NWORK ), |
Line 1511
|
Line 1637
|
* Perform bidiagonal SVD, computing left singular vectors |
* Perform bidiagonal SVD, computing left singular vectors |
* of bidiagonal matrix in RWORK(IRU) and computing right |
* of bidiagonal matrix in RWORK(IRU) and computing right |
* singular vectors of bidiagonal matrix in RWORK(IRVT) |
* singular vectors of bidiagonal matrix in RWORK(IRVT) |
* (CWorkspace: need 0) |
* CWorkspace: need 0 |
* (RWorkspace: need BDSPAC) |
* RWorkspace: need M [e] + M*M [RU] + M*M [RVT] + BDSPAC |
* |
* |
IRU = IE + M |
IRU = IE + M |
IRVT = IRU + M*M |
IRVT = IRU + M*M |
Line 1523
|
Line 1649
|
* |
* |
* Copy real matrix RWORK(IRU) to complex matrix U |
* Copy real matrix RWORK(IRU) to complex matrix U |
* Overwrite U by left singular vectors of L |
* Overwrite U by left singular vectors of L |
* (CWorkspace: need M*M+3*M, prefer M*M+2*M+M*NB) |
* CWorkspace: need M*M [L] + 2*M [tauq, taup] + M [work] |
* (RWorkspace: 0) |
* CWorkspace: prefer M*M [L] + 2*M [tauq, taup] + M*NB [work] |
|
* RWorkspace: need 0 |
* |
* |
CALL ZLACP2( 'F', M, M, RWORK( IRU ), M, U, LDU ) |
CALL ZLACP2( 'F', M, M, RWORK( IRU ), M, U, LDU ) |
CALL ZUNMBR( 'Q', 'L', 'N', M, M, M, WORK( IL ), LDWRKL, |
CALL ZUNMBR( 'Q', 'L', 'N', M, M, M, WORK( IL ), LDWRKL, |
Line 1533
|
Line 1660
|
* |
* |
* Copy real matrix RWORK(IRVT) to complex matrix VT |
* Copy real matrix RWORK(IRVT) to complex matrix VT |
* Overwrite VT by left singular vectors of L |
* Overwrite VT by left singular vectors of L |
* (CWorkspace: need M*M+3*M, prefer M*M+2*M+M*NB) |
* CWorkspace: need M*M [L] + 2*M [tauq, taup] + M [work] |
* (RWorkspace: 0) |
* CWorkspace: prefer M*M [L] + 2*M [tauq, taup] + M*NB [work] |
|
* RWorkspace: need 0 |
* |
* |
CALL ZLACP2( 'F', M, M, RWORK( IRVT ), M, VT, LDVT ) |
CALL ZLACP2( 'F', M, M, RWORK( IRVT ), M, VT, LDVT ) |
CALL ZUNMBR( 'P', 'R', 'C', M, M, M, WORK( IL ), LDWRKL, |
CALL ZUNMBR( 'P', 'R', 'C', M, M, M, WORK( IL ), LDWRKL, |
Line 1543
|
Line 1671
|
* |
* |
* Copy VT to WORK(IL), multiply right singular vectors of L |
* Copy VT to WORK(IL), multiply right singular vectors of L |
* in WORK(IL) by Q in A, storing result in VT |
* in WORK(IL) by Q in A, storing result in VT |
* (CWorkspace: need M*M) |
* CWorkspace: need M*M [L] |
* (RWorkspace: 0) |
* RWorkspace: need 0 |
* |
* |
CALL ZLACPY( 'F', M, M, VT, LDVT, WORK( IL ), LDWRKL ) |
CALL ZLACPY( 'F', M, M, VT, LDVT, WORK( IL ), LDWRKL ) |
CALL ZGEMM( 'N', 'N', M, N, M, CONE, WORK( IL ), LDWRKL, |
CALL ZGEMM( 'N', 'N', M, N, M, CONE, WORK( IL ), LDWRKL, |
Line 1552
|
Line 1680
|
* |
* |
ELSE IF( WNTQA ) THEN |
ELSE IF( WNTQA ) THEN |
* |
* |
* Path 9t (N much larger than M, JOBZ='A') |
* Path 4t (N >> M, JOBZ='A') |
* N right singular vectors to be computed in VT and |
* N right singular vectors to be computed in VT and |
* M left singular vectors to be computed in U |
* M left singular vectors to be computed in U |
* |
* |
Line 1565
|
Line 1693
|
NWORK = ITAU + M |
NWORK = ITAU + M |
* |
* |
* Compute A=L*Q, copying result to VT |
* Compute A=L*Q, copying result to VT |
* (CWorkspace: need 2*M, prefer M+M*NB) |
* CWorkspace: need M*M [VT] + M [tau] + M [work] |
* (RWorkspace: 0) |
* CWorkspace: prefer M*M [VT] + M [tau] + M*NB [work] |
|
* RWorkspace: need 0 |
* |
* |
CALL ZGELQF( M, N, A, LDA, WORK( ITAU ), WORK( NWORK ), |
CALL ZGELQF( M, N, A, LDA, WORK( ITAU ), WORK( NWORK ), |
$ LWORK-NWORK+1, IERR ) |
$ LWORK-NWORK+1, IERR ) |
CALL ZLACPY( 'U', M, N, A, LDA, VT, LDVT ) |
CALL ZLACPY( 'U', M, N, A, LDA, VT, LDVT ) |
* |
* |
* Generate Q in VT |
* Generate Q in VT |
* (CWorkspace: need M+N, prefer M+N*NB) |
* CWorkspace: need M*M [VT] + M [tau] + N [work] |
* (RWorkspace: 0) |
* CWorkspace: prefer M*M [VT] + M [tau] + N*NB [work] |
|
* RWorkspace: need 0 |
* |
* |
CALL ZUNGLQ( N, N, M, VT, LDVT, WORK( ITAU ), |
CALL ZUNGLQ( N, N, M, VT, LDVT, WORK( ITAU ), |
$ WORK( NWORK ), LWORK-NWORK+1, IERR ) |
$ WORK( NWORK ), LWORK-NWORK+1, IERR ) |
Line 1589
|
Line 1719
|
NWORK = ITAUP + M |
NWORK = ITAUP + M |
* |
* |
* Bidiagonalize L in A |
* Bidiagonalize L in A |
* (CWorkspace: need M*M+3*M, prefer M*M+2*M+2*M*NB) |
* CWorkspace: need M*M [VT] + 2*M [tauq, taup] + M [work] |
* (RWorkspace: need M) |
* CWorkspace: prefer M*M [VT] + 2*M [tauq, taup] + 2*M*NB [work] |
|
* RWorkspace: need M [e] |
* |
* |
CALL ZGEBRD( M, M, A, LDA, S, RWORK( IE ), WORK( ITAUQ ), |
CALL ZGEBRD( M, M, A, LDA, S, RWORK( IE ), WORK( ITAUQ ), |
$ WORK( ITAUP ), WORK( NWORK ), LWORK-NWORK+1, |
$ WORK( ITAUP ), WORK( NWORK ), LWORK-NWORK+1, |
Line 1599
|
Line 1730
|
* Perform bidiagonal SVD, computing left singular vectors |
* Perform bidiagonal SVD, computing left singular vectors |
* of bidiagonal matrix in RWORK(IRU) and computing right |
* of bidiagonal matrix in RWORK(IRU) and computing right |
* singular vectors of bidiagonal matrix in RWORK(IRVT) |
* singular vectors of bidiagonal matrix in RWORK(IRVT) |
* (CWorkspace: need 0) |
* CWorkspace: need 0 |
* (RWorkspace: need BDSPAC) |
* RWorkspace: need M [e] + M*M [RU] + M*M [RVT] + BDSPAC |
* |
* |
IRU = IE + M |
IRU = IE + M |
IRVT = IRU + M*M |
IRVT = IRU + M*M |
Line 1611
|
Line 1742
|
* |
* |
* Copy real matrix RWORK(IRU) to complex matrix U |
* Copy real matrix RWORK(IRU) to complex matrix U |
* Overwrite U by left singular vectors of L |
* Overwrite U by left singular vectors of L |
* (CWorkspace: need 3*M, prefer 2*M+M*NB) |
* CWorkspace: need M*M [VT] + 2*M [tauq, taup] + M [work] |
* (RWorkspace: 0) |
* CWorkspace: prefer M*M [VT] + 2*M [tauq, taup] + M*NB [work] |
|
* RWorkspace: need 0 |
* |
* |
CALL ZLACP2( 'F', M, M, RWORK( IRU ), M, U, LDU ) |
CALL ZLACP2( 'F', M, M, RWORK( IRU ), M, U, LDU ) |
CALL ZUNMBR( 'Q', 'L', 'N', M, M, M, A, LDA, |
CALL ZUNMBR( 'Q', 'L', 'N', M, M, M, A, LDA, |
Line 1621
|
Line 1753
|
* |
* |
* Copy real matrix RWORK(IRVT) to complex matrix WORK(IVT) |
* Copy real matrix RWORK(IRVT) to complex matrix WORK(IVT) |
* Overwrite WORK(IVT) by right singular vectors of L |
* Overwrite WORK(IVT) by right singular vectors of L |
* (CWorkspace: need M*M+3*M, prefer M*M+2*M+M*NB) |
* CWorkspace: need M*M [VT] + 2*M [tauq, taup] + M [work] |
* (RWorkspace: 0) |
* CWorkspace: prefer M*M [VT] + 2*M [tauq, taup] + M*NB [work] |
|
* RWorkspace: need 0 |
* |
* |
CALL ZLACP2( 'F', M, M, RWORK( IRVT ), M, WORK( IVT ), |
CALL ZLACP2( 'F', M, M, RWORK( IRVT ), M, WORK( IVT ), |
$ LDWKVT ) |
$ LDWKVT ) |
Line 1632
|
Line 1765
|
* |
* |
* Multiply right singular vectors of L in WORK(IVT) by |
* Multiply right singular vectors of L in WORK(IVT) by |
* Q in VT, storing result in A |
* Q in VT, storing result in A |
* (CWorkspace: need M*M) |
* CWorkspace: need M*M [VT] |
* (RWorkspace: 0) |
* RWorkspace: need 0 |
* |
* |
CALL ZGEMM( 'N', 'N', M, N, M, CONE, WORK( IVT ), LDWKVT, |
CALL ZGEMM( 'N', 'N', M, N, M, CONE, WORK( IVT ), LDWKVT, |
$ VT, LDVT, CZERO, A, LDA ) |
$ VT, LDVT, CZERO, A, LDA ) |
Line 1648
|
Line 1781
|
* |
* |
* MNTHR2 <= N < MNTHR1 |
* MNTHR2 <= N < MNTHR1 |
* |
* |
* Path 5t (N much larger than M, but not as much as MNTHR1) |
* Path 5t (N >> M, but not as much as MNTHR1) |
* Reduce to bidiagonal form without QR decomposition, use |
* Reduce to bidiagonal form without QR decomposition, use |
* ZUNGBR and matrix multiplication to compute singular vectors |
* ZUNGBR and matrix multiplication to compute singular vectors |
* |
* |
* |
|
IE = 1 |
IE = 1 |
NRWORK = IE + M |
NRWORK = IE + M |
ITAUQ = 1 |
ITAUQ = 1 |
Line 1660
|
Line 1792
|
NWORK = ITAUP + M |
NWORK = ITAUP + M |
* |
* |
* Bidiagonalize A |
* Bidiagonalize A |
* (CWorkspace: need 2*M+N, prefer 2*M+(M+N)*NB) |
* CWorkspace: need 2*M [tauq, taup] + N [work] |
* (RWorkspace: M) |
* CWorkspace: prefer 2*M [tauq, taup] + (M+N)*NB [work] |
|
* RWorkspace: need M [e] |
* |
* |
CALL ZGEBRD( M, N, A, LDA, S, RWORK( IE ), WORK( ITAUQ ), |
CALL ZGEBRD( M, N, A, LDA, S, RWORK( IE ), WORK( ITAUQ ), |
$ WORK( ITAUP ), WORK( NWORK ), LWORK-NWORK+1, |
$ WORK( ITAUP ), WORK( NWORK ), LWORK-NWORK+1, |
Line 1669
|
Line 1802
|
* |
* |
IF( WNTQN ) THEN |
IF( WNTQN ) THEN |
* |
* |
|
* Path 5tn (N >> M, JOBZ='N') |
* Compute singular values only |
* Compute singular values only |
* (Cworkspace: 0) |
* CWorkspace: need 0 |
* (Rworkspace: need BDSPAN) |
* RWorkspace: need M [e] + BDSPAC |
* |
* |
CALL DBDSDC( 'L', 'N', M, S, RWORK( IE ), DUM, 1, DUM, 1, |
CALL DBDSDC( 'L', 'N', M, S, RWORK( IE ), DUM,1,DUM,1, |
$ DUM, IDUM, RWORK( NRWORK ), IWORK, INFO ) |
$ DUM, IDUM, RWORK( NRWORK ), IWORK, INFO ) |
ELSE IF( WNTQO ) THEN |
ELSE IF( WNTQO ) THEN |
IRVT = NRWORK |
IRVT = NRWORK |
Line 1681
|
Line 1815
|
NRWORK = IRU + M*M |
NRWORK = IRU + M*M |
IVT = NWORK |
IVT = NWORK |
* |
* |
|
* Path 5to (N >> M, JOBZ='O') |
* Copy A to U, generate Q |
* Copy A to U, generate Q |
* (Cworkspace: need 2*M, prefer M+M*NB) |
* CWorkspace: need 2*M [tauq, taup] + M [work] |
* (Rworkspace: 0) |
* CWorkspace: prefer 2*M [tauq, taup] + M*NB [work] |
|
* RWorkspace: need 0 |
* |
* |
CALL ZLACPY( 'L', M, M, A, LDA, U, LDU ) |
CALL ZLACPY( 'L', M, M, A, LDA, U, LDU ) |
CALL ZUNGBR( 'Q', M, M, N, U, LDU, WORK( ITAUQ ), |
CALL ZUNGBR( 'Q', M, M, N, U, LDU, WORK( ITAUQ ), |
$ WORK( NWORK ), LWORK-NWORK+1, IERR ) |
$ WORK( NWORK ), LWORK-NWORK+1, IERR ) |
* |
* |
* Generate P**H in A |
* Generate P**H in A |
* (Cworkspace: need 2*M, prefer M+M*NB) |
* CWorkspace: need 2*M [tauq, taup] + M [work] |
* (Rworkspace: 0) |
* CWorkspace: prefer 2*M [tauq, taup] + M*NB [work] |
|
* RWorkspace: need 0 |
* |
* |
CALL ZUNGBR( 'P', M, N, M, A, LDA, WORK( ITAUP ), |
CALL ZUNGBR( 'P', M, N, M, A, LDA, WORK( ITAUP ), |
$ WORK( NWORK ), LWORK-NWORK+1, IERR ) |
$ WORK( NWORK ), LWORK-NWORK+1, IERR ) |
* |
* |
LDWKVT = M |
LDWKVT = M |
IF( LWORK.GE.M*N+3*M ) THEN |
IF( LWORK .GE. M*N + 3*M ) THEN |
* |
* |
* WORK( IVT ) is M by N |
* WORK( IVT ) is M by N |
* |
* |
Line 1707
|
Line 1844
|
* |
* |
* WORK( IVT ) is M by CHUNK |
* WORK( IVT ) is M by CHUNK |
* |
* |
CHUNK = ( LWORK-3*M ) / M |
CHUNK = ( LWORK - 3*M ) / M |
NWORK = IVT + LDWKVT*CHUNK |
NWORK = IVT + LDWKVT*CHUNK |
END IF |
END IF |
* |
* |
* Perform bidiagonal SVD, computing left singular vectors |
* Perform bidiagonal SVD, computing left singular vectors |
* of bidiagonal matrix in RWORK(IRU) and computing right |
* of bidiagonal matrix in RWORK(IRU) and computing right |
* singular vectors of bidiagonal matrix in RWORK(IRVT) |
* singular vectors of bidiagonal matrix in RWORK(IRVT) |
* (CWorkspace: need 0) |
* CWorkspace: need 0 |
* (RWorkspace: need BDSPAC) |
* RWorkspace: need M [e] + M*M [RVT] + M*M [RU] + BDSPAC |
* |
* |
CALL DBDSDC( 'L', 'I', M, S, RWORK( IE ), RWORK( IRU ), |
CALL DBDSDC( 'L', 'I', M, S, RWORK( IE ), RWORK( IRU ), |
$ M, RWORK( IRVT ), M, DUM, IDUM, |
$ M, RWORK( IRVT ), M, DUM, IDUM, |
Line 1723
|
Line 1860
|
* |
* |
* Multiply Q in U by real matrix RWORK(IRVT) |
* Multiply Q in U by real matrix RWORK(IRVT) |
* storing the result in WORK(IVT), copying to U |
* storing the result in WORK(IVT), copying to U |
* (Cworkspace: need 0) |
* CWorkspace: need 2*M [tauq, taup] + M*M [VT] |
* (Rworkspace: need 2*M*M) |
* RWorkspace: need M [e] + M*M [RVT] + M*M [RU] + 2*M*M [rwork] |
* |
* |
CALL ZLACRM( M, M, U, LDU, RWORK( IRU ), M, WORK( IVT ), |
CALL ZLACRM( M, M, U, LDU, RWORK( IRU ), M, WORK( IVT ), |
$ LDWKVT, RWORK( NRWORK ) ) |
$ LDWKVT, RWORK( NRWORK ) ) |
Line 1732
|
Line 1869
|
* |
* |
* Multiply RWORK(IRVT) by P**H in A, storing the |
* Multiply RWORK(IRVT) by P**H in A, storing the |
* result in WORK(IVT), copying to A |
* result in WORK(IVT), copying to A |
* (CWorkspace: need M*M, prefer M*N) |
* CWorkspace: need 2*M [tauq, taup] + M*M [VT] |
* (Rworkspace: need 2*M*M, prefer 2*M*N) |
* CWorkspace: prefer 2*M [tauq, taup] + M*N [VT] |
|
* RWorkspace: need M [e] + M*M [RVT] + 2*M*M [rwork] |
|
* RWorkspace: prefer M [e] + M*M [RVT] + 2*M*N [rwork] < M + 5*M*M since N < 2*M here |
* |
* |
NRWORK = IRU |
NRWORK = IRU |
DO 50 I = 1, N, CHUNK |
DO 50 I = 1, N, CHUNK |
Line 1745
|
Line 1884
|
50 CONTINUE |
50 CONTINUE |
ELSE IF( WNTQS ) THEN |
ELSE IF( WNTQS ) THEN |
* |
* |
|
* Path 5ts (N >> M, JOBZ='S') |
* Copy A to U, generate Q |
* Copy A to U, generate Q |
* (Cworkspace: need 2*M, prefer M+M*NB) |
* CWorkspace: need 2*M [tauq, taup] + M [work] |
* (Rworkspace: 0) |
* CWorkspace: prefer 2*M [tauq, taup] + M*NB [work] |
|
* RWorkspace: need 0 |
* |
* |
CALL ZLACPY( 'L', M, M, A, LDA, U, LDU ) |
CALL ZLACPY( 'L', M, M, A, LDA, U, LDU ) |
CALL ZUNGBR( 'Q', M, M, N, U, LDU, WORK( ITAUQ ), |
CALL ZUNGBR( 'Q', M, M, N, U, LDU, WORK( ITAUQ ), |
$ WORK( NWORK ), LWORK-NWORK+1, IERR ) |
$ WORK( NWORK ), LWORK-NWORK+1, IERR ) |
* |
* |
* Copy A to VT, generate P**H |
* Copy A to VT, generate P**H |
* (Cworkspace: need 2*M, prefer M+M*NB) |
* CWorkspace: need 2*M [tauq, taup] + M [work] |
* (Rworkspace: 0) |
* CWorkspace: prefer 2*M [tauq, taup] + M*NB [work] |
|
* RWorkspace: need 0 |
* |
* |
CALL ZLACPY( 'U', M, N, A, LDA, VT, LDVT ) |
CALL ZLACPY( 'U', M, N, A, LDA, VT, LDVT ) |
CALL ZUNGBR( 'P', M, N, M, VT, LDVT, WORK( ITAUP ), |
CALL ZUNGBR( 'P', M, N, M, VT, LDVT, WORK( ITAUP ), |
Line 1764
|
Line 1906
|
* Perform bidiagonal SVD, computing left singular vectors |
* Perform bidiagonal SVD, computing left singular vectors |
* of bidiagonal matrix in RWORK(IRU) and computing right |
* of bidiagonal matrix in RWORK(IRU) and computing right |
* singular vectors of bidiagonal matrix in RWORK(IRVT) |
* singular vectors of bidiagonal matrix in RWORK(IRVT) |
* (CWorkspace: need 0) |
* CWorkspace: need 0 |
* (RWorkspace: need BDSPAC) |
* RWorkspace: need M [e] + M*M [RVT] + M*M [RU] + BDSPAC |
* |
* |
IRVT = NRWORK |
IRVT = NRWORK |
IRU = IRVT + M*M |
IRU = IRVT + M*M |
Line 1776
|
Line 1918
|
* |
* |
* Multiply Q in U by real matrix RWORK(IRU), storing the |
* Multiply Q in U by real matrix RWORK(IRU), storing the |
* result in A, copying to U |
* result in A, copying to U |
* (CWorkspace: need 0) |
* CWorkspace: need 0 |
* (Rworkspace: need 3*M*M) |
* RWorkspace: need M [e] + M*M [RVT] + M*M [RU] + 2*M*M [rwork] |
* |
* |
CALL ZLACRM( M, M, U, LDU, RWORK( IRU ), M, A, LDA, |
CALL ZLACRM( M, M, U, LDU, RWORK( IRU ), M, A, LDA, |
$ RWORK( NRWORK ) ) |
$ RWORK( NRWORK ) ) |
Line 1785
|
Line 1927
|
* |
* |
* Multiply real matrix RWORK(IRVT) by P**H in VT, |
* Multiply real matrix RWORK(IRVT) by P**H in VT, |
* storing the result in A, copying to VT |
* storing the result in A, copying to VT |
* (Cworkspace: need 0) |
* CWorkspace: need 0 |
* (Rworkspace: need M*M+2*M*N) |
* RWorkspace: need M [e] + M*M [RVT] + 2*M*N [rwork] < M + 5*M*M since N < 2*M here |
* |
* |
NRWORK = IRU |
NRWORK = IRU |
CALL ZLARCM( M, N, RWORK( IRVT ), M, VT, LDVT, A, LDA, |
CALL ZLARCM( M, N, RWORK( IRVT ), M, VT, LDVT, A, LDA, |
Line 1794
|
Line 1936
|
CALL ZLACPY( 'F', M, N, A, LDA, VT, LDVT ) |
CALL ZLACPY( 'F', M, N, A, LDA, VT, LDVT ) |
ELSE |
ELSE |
* |
* |
|
* Path 5ta (N >> M, JOBZ='A') |
* Copy A to U, generate Q |
* Copy A to U, generate Q |
* (Cworkspace: need 2*M, prefer M+M*NB) |
* CWorkspace: need 2*M [tauq, taup] + M [work] |
* (Rworkspace: 0) |
* CWorkspace: prefer 2*M [tauq, taup] + M*NB [work] |
|
* RWorkspace: need 0 |
* |
* |
CALL ZLACPY( 'L', M, M, A, LDA, U, LDU ) |
CALL ZLACPY( 'L', M, M, A, LDA, U, LDU ) |
CALL ZUNGBR( 'Q', M, M, N, U, LDU, WORK( ITAUQ ), |
CALL ZUNGBR( 'Q', M, M, N, U, LDU, WORK( ITAUQ ), |
$ WORK( NWORK ), LWORK-NWORK+1, IERR ) |
$ WORK( NWORK ), LWORK-NWORK+1, IERR ) |
* |
* |
* Copy A to VT, generate P**H |
* Copy A to VT, generate P**H |
* (Cworkspace: need 2*M, prefer M+M*NB) |
* CWorkspace: need 2*M [tauq, taup] + N [work] |
* (Rworkspace: 0) |
* CWorkspace: prefer 2*M [tauq, taup] + N*NB [work] |
|
* RWorkspace: need 0 |
* |
* |
CALL ZLACPY( 'U', M, N, A, LDA, VT, LDVT ) |
CALL ZLACPY( 'U', M, N, A, LDA, VT, LDVT ) |
CALL ZUNGBR( 'P', N, N, M, VT, LDVT, WORK( ITAUP ), |
CALL ZUNGBR( 'P', N, N, M, VT, LDVT, WORK( ITAUP ), |
Line 1813
|
Line 1958
|
* Perform bidiagonal SVD, computing left singular vectors |
* Perform bidiagonal SVD, computing left singular vectors |
* of bidiagonal matrix in RWORK(IRU) and computing right |
* of bidiagonal matrix in RWORK(IRU) and computing right |
* singular vectors of bidiagonal matrix in RWORK(IRVT) |
* singular vectors of bidiagonal matrix in RWORK(IRVT) |
* (CWorkspace: need 0) |
* CWorkspace: need 0 |
* (RWorkspace: need BDSPAC) |
* RWorkspace: need M [e] + M*M [RVT] + M*M [RU] + BDSPAC |
* |
* |
IRVT = NRWORK |
IRVT = NRWORK |
IRU = IRVT + M*M |
IRU = IRVT + M*M |
Line 1825
|
Line 1970
|
* |
* |
* Multiply Q in U by real matrix RWORK(IRU), storing the |
* Multiply Q in U by real matrix RWORK(IRU), storing the |
* result in A, copying to U |
* result in A, copying to U |
* (CWorkspace: need 0) |
* CWorkspace: need 0 |
* (Rworkspace: need 3*M*M) |
* RWorkspace: need M [e] + M*M [RVT] + M*M [RU] + 2*M*M [rwork] |
* |
* |
CALL ZLACRM( M, M, U, LDU, RWORK( IRU ), M, A, LDA, |
CALL ZLACRM( M, M, U, LDU, RWORK( IRU ), M, A, LDA, |
$ RWORK( NRWORK ) ) |
$ RWORK( NRWORK ) ) |
Line 1834
|
Line 1979
|
* |
* |
* Multiply real matrix RWORK(IRVT) by P**H in VT, |
* Multiply real matrix RWORK(IRVT) by P**H in VT, |
* storing the result in A, copying to VT |
* storing the result in A, copying to VT |
* (Cworkspace: need 0) |
* CWorkspace: need 0 |
* (Rworkspace: need M*M+2*M*N) |
* RWorkspace: need M [e] + M*M [RVT] + 2*M*N [rwork] < M + 5*M*M since N < 2*M here |
* |
* |
|
NRWORK = IRU |
CALL ZLARCM( M, N, RWORK( IRVT ), M, VT, LDVT, A, LDA, |
CALL ZLARCM( M, N, RWORK( IRVT ), M, VT, LDVT, A, LDA, |
$ RWORK( NRWORK ) ) |
$ RWORK( NRWORK ) ) |
CALL ZLACPY( 'F', M, N, A, LDA, VT, LDVT ) |
CALL ZLACPY( 'F', M, N, A, LDA, VT, LDVT ) |
Line 1846
|
Line 1992
|
* |
* |
* N .LT. MNTHR2 |
* N .LT. MNTHR2 |
* |
* |
* Path 6t (N greater than M, but not much larger) |
* Path 6t (N > M, but not much larger) |
* Reduce to bidiagonal form without LQ decomposition |
* Reduce to bidiagonal form without LQ decomposition |
* Use ZUNMBR to compute singular vectors |
* Use ZUNMBR to compute singular vectors |
* |
* |
Line 1857
|
Line 2003
|
NWORK = ITAUP + M |
NWORK = ITAUP + M |
* |
* |
* Bidiagonalize A |
* Bidiagonalize A |
* (CWorkspace: need 2*M+N, prefer 2*M+(M+N)*NB) |
* CWorkspace: need 2*M [tauq, taup] + N [work] |
* (RWorkspace: M) |
* CWorkspace: prefer 2*M [tauq, taup] + (M+N)*NB [work] |
|
* RWorkspace: need M [e] |
* |
* |
CALL ZGEBRD( M, N, A, LDA, S, RWORK( IE ), WORK( ITAUQ ), |
CALL ZGEBRD( M, N, A, LDA, S, RWORK( IE ), WORK( ITAUQ ), |
$ WORK( ITAUP ), WORK( NWORK ), LWORK-NWORK+1, |
$ WORK( ITAUP ), WORK( NWORK ), LWORK-NWORK+1, |
$ IERR ) |
$ IERR ) |
IF( WNTQN ) THEN |
IF( WNTQN ) THEN |
* |
* |
|
* Path 6tn (N > M, JOBZ='N') |
* Compute singular values only |
* Compute singular values only |
* (Cworkspace: 0) |
* CWorkspace: need 0 |
* (Rworkspace: need BDSPAN) |
* RWorkspace: need M [e] + BDSPAC |
* |
* |
CALL DBDSDC( 'L', 'N', M, S, RWORK( IE ), DUM, 1, DUM, 1, |
CALL DBDSDC( 'L', 'N', M, S, RWORK( IE ), DUM,1,DUM,1, |
$ DUM, IDUM, RWORK( NRWORK ), IWORK, INFO ) |
$ DUM, IDUM, RWORK( NRWORK ), IWORK, INFO ) |
ELSE IF( WNTQO ) THEN |
ELSE IF( WNTQO ) THEN |
|
* Path 6to (N > M, JOBZ='O') |
LDWKVT = M |
LDWKVT = M |
IVT = NWORK |
IVT = NWORK |
IF( LWORK.GE.M*N+3*M ) THEN |
IF( LWORK .GE. M*N + 3*M ) THEN |
* |
* |
* WORK( IVT ) is M by N |
* WORK( IVT ) is M by N |
* |
* |
Line 1885
|
Line 2034
|
* |
* |
* WORK( IVT ) is M by CHUNK |
* WORK( IVT ) is M by CHUNK |
* |
* |
CHUNK = ( LWORK-3*M ) / M |
CHUNK = ( LWORK - 3*M ) / M |
NWORK = IVT + LDWKVT*CHUNK |
NWORK = IVT + LDWKVT*CHUNK |
END IF |
END IF |
* |
* |
* Perform bidiagonal SVD, computing left singular vectors |
* Perform bidiagonal SVD, computing left singular vectors |
* of bidiagonal matrix in RWORK(IRU) and computing right |
* of bidiagonal matrix in RWORK(IRU) and computing right |
* singular vectors of bidiagonal matrix in RWORK(IRVT) |
* singular vectors of bidiagonal matrix in RWORK(IRVT) |
* (CWorkspace: need 0) |
* CWorkspace: need 0 |
* (RWorkspace: need BDSPAC) |
* RWorkspace: need M [e] + M*M [RVT] + M*M [RU] + BDSPAC |
* |
* |
IRVT = NRWORK |
IRVT = NRWORK |
IRU = IRVT + M*M |
IRU = IRVT + M*M |
Line 1904
|
Line 2053
|
* |
* |
* Copy real matrix RWORK(IRU) to complex matrix U |
* Copy real matrix RWORK(IRU) to complex matrix U |
* Overwrite U by left singular vectors of A |
* Overwrite U by left singular vectors of A |
* (Cworkspace: need 2*M, prefer M+M*NB) |
* CWorkspace: need 2*M [tauq, taup] + M*M [VT] + M [work] |
* (Rworkspace: need 0) |
* CWorkspace: prefer 2*M [tauq, taup] + M*M [VT] + M*NB [work] |
|
* RWorkspace: need M [e] + M*M [RVT] + M*M [RU] |
* |
* |
CALL ZLACP2( 'F', M, M, RWORK( IRU ), M, U, LDU ) |
CALL ZLACP2( 'F', M, M, RWORK( IRU ), M, U, LDU ) |
CALL ZUNMBR( 'Q', 'L', 'N', M, M, N, A, LDA, |
CALL ZUNMBR( 'Q', 'L', 'N', M, M, N, A, LDA, |
$ WORK( ITAUQ ), U, LDU, WORK( NWORK ), |
$ WORK( ITAUQ ), U, LDU, WORK( NWORK ), |
$ LWORK-NWORK+1, IERR ) |
$ LWORK-NWORK+1, IERR ) |
* |
* |
IF( LWORK.GE.M*N+3*M ) THEN |
IF( LWORK .GE. M*N + 3*M ) THEN |
* |
* |
* Copy real matrix RWORK(IRVT) to complex matrix WORK(IVT) |
* Path 6to-fast |
* Overwrite WORK(IVT) by right singular vectors of A, |
* Copy real matrix RWORK(IRVT) to complex matrix WORK(IVT) |
* copying to A |
* Overwrite WORK(IVT) by right singular vectors of A, |
* (Cworkspace: need M*N+2*M, prefer M*N+M+M*NB) |
* copying to A |
* (Rworkspace: need 0) |
* CWorkspace: need 2*M [tauq, taup] + M*N [VT] + M [work] |
|
* CWorkspace: prefer 2*M [tauq, taup] + M*N [VT] + M*NB [work] |
|
* RWorkspace: need M [e] + M*M [RVT] |
* |
* |
CALL ZLACP2( 'F', M, M, RWORK( IRVT ), M, WORK( IVT ), |
CALL ZLACP2( 'F', M, M, RWORK( IRVT ), M, WORK( IVT ), |
$ LDWKVT ) |
$ LDWKVT ) |
Line 1928
|
Line 2080
|
CALL ZLACPY( 'F', M, N, WORK( IVT ), LDWKVT, A, LDA ) |
CALL ZLACPY( 'F', M, N, WORK( IVT ), LDWKVT, A, LDA ) |
ELSE |
ELSE |
* |
* |
|
* Path 6to-slow |
* Generate P**H in A |
* Generate P**H in A |
* (Cworkspace: need 2*M, prefer M+M*NB) |
* CWorkspace: need 2*M [tauq, taup] + M*M [VT] + M [work] |
* (Rworkspace: need 0) |
* CWorkspace: prefer 2*M [tauq, taup] + M*M [VT] + M*NB [work] |
|
* RWorkspace: need 0 |
* |
* |
CALL ZUNGBR( 'P', M, N, M, A, LDA, WORK( ITAUP ), |
CALL ZUNGBR( 'P', M, N, M, A, LDA, WORK( ITAUP ), |
$ WORK( NWORK ), LWORK-NWORK+1, IERR ) |
$ WORK( NWORK ), LWORK-NWORK+1, IERR ) |
* |
* |
* Multiply Q in A by real matrix RWORK(IRU), storing the |
* Multiply Q in A by real matrix RWORK(IRU), storing the |
* result in WORK(IU), copying to A |
* result in WORK(IU), copying to A |
* (CWorkspace: need M*M, prefer M*N) |
* CWorkspace: need 2*M [tauq, taup] + M*M [VT] |
* (Rworkspace: need 3*M*M, prefer M*M+2*M*N) |
* CWorkspace: prefer 2*M [tauq, taup] + M*N [VT] |
|
* RWorkspace: need M [e] + M*M [RVT] + 2*M*M [rwork] |
|
* RWorkspace: prefer M [e] + M*M [RVT] + 2*M*N [rwork] < M + 5*M*M since N < 2*M here |
* |
* |
NRWORK = IRU |
NRWORK = IRU |
DO 60 I = 1, N, CHUNK |
DO 60 I = 1, N, CHUNK |
Line 1952
|
Line 2108
|
END IF |
END IF |
ELSE IF( WNTQS ) THEN |
ELSE IF( WNTQS ) THEN |
* |
* |
|
* Path 6ts (N > M, JOBZ='S') |
* Perform bidiagonal SVD, computing left singular vectors |
* Perform bidiagonal SVD, computing left singular vectors |
* of bidiagonal matrix in RWORK(IRU) and computing right |
* of bidiagonal matrix in RWORK(IRU) and computing right |
* singular vectors of bidiagonal matrix in RWORK(IRVT) |
* singular vectors of bidiagonal matrix in RWORK(IRVT) |
* (CWorkspace: need 0) |
* CWorkspace: need 0 |
* (RWorkspace: need BDSPAC) |
* RWorkspace: need M [e] + M*M [RVT] + M*M [RU] + BDSPAC |
* |
* |
IRVT = NRWORK |
IRVT = NRWORK |
IRU = IRVT + M*M |
IRU = IRVT + M*M |
Line 1967
|
Line 2124
|
* |
* |
* Copy real matrix RWORK(IRU) to complex matrix U |
* Copy real matrix RWORK(IRU) to complex matrix U |
* Overwrite U by left singular vectors of A |
* Overwrite U by left singular vectors of A |
* (CWorkspace: need 3*M, prefer 2*M+M*NB) |
* CWorkspace: need 2*M [tauq, taup] + M [work] |
* (RWorkspace: M*M) |
* CWorkspace: prefer 2*M [tauq, taup] + M*NB [work] |
|
* RWorkspace: need M [e] + M*M [RVT] + M*M [RU] |
* |
* |
CALL ZLACP2( 'F', M, M, RWORK( IRU ), M, U, LDU ) |
CALL ZLACP2( 'F', M, M, RWORK( IRU ), M, U, LDU ) |
CALL ZUNMBR( 'Q', 'L', 'N', M, M, N, A, LDA, |
CALL ZUNMBR( 'Q', 'L', 'N', M, M, N, A, LDA, |
Line 1977
|
Line 2135
|
* |
* |
* Copy real matrix RWORK(IRVT) to complex matrix VT |
* Copy real matrix RWORK(IRVT) to complex matrix VT |
* Overwrite VT by right singular vectors of A |
* Overwrite VT by right singular vectors of A |
* (CWorkspace: need 3*M, prefer 2*M+M*NB) |
* CWorkspace: need 2*M [tauq, taup] + M [work] |
* (RWorkspace: M*M) |
* CWorkspace: prefer 2*M [tauq, taup] + M*NB [work] |
|
* RWorkspace: need M [e] + M*M [RVT] |
* |
* |
CALL ZLASET( 'F', M, N, CZERO, CZERO, VT, LDVT ) |
CALL ZLASET( 'F', M, N, CZERO, CZERO, VT, LDVT ) |
CALL ZLACP2( 'F', M, M, RWORK( IRVT ), M, VT, LDVT ) |
CALL ZLACP2( 'F', M, M, RWORK( IRVT ), M, VT, LDVT ) |
Line 1987
|
Line 2146
|
$ LWORK-NWORK+1, IERR ) |
$ LWORK-NWORK+1, IERR ) |
ELSE |
ELSE |
* |
* |
|
* Path 6ta (N > M, JOBZ='A') |
* Perform bidiagonal SVD, computing left singular vectors |
* Perform bidiagonal SVD, computing left singular vectors |
* of bidiagonal matrix in RWORK(IRU) and computing right |
* of bidiagonal matrix in RWORK(IRU) and computing right |
* singular vectors of bidiagonal matrix in RWORK(IRVT) |
* singular vectors of bidiagonal matrix in RWORK(IRVT) |
* (CWorkspace: need 0) |
* CWorkspace: need 0 |
* (RWorkspace: need BDSPAC) |
* RWorkspace: need M [e] + M*M [RVT] + M*M [RU] + BDSPAC |
* |
* |
IRVT = NRWORK |
IRVT = NRWORK |
IRU = IRVT + M*M |
IRU = IRVT + M*M |
Line 2003
|
Line 2163
|
* |
* |
* Copy real matrix RWORK(IRU) to complex matrix U |
* Copy real matrix RWORK(IRU) to complex matrix U |
* Overwrite U by left singular vectors of A |
* Overwrite U by left singular vectors of A |
* (CWorkspace: need 3*M, prefer 2*M+M*NB) |
* CWorkspace: need 2*M [tauq, taup] + M [work] |
* (RWorkspace: M*M) |
* CWorkspace: prefer 2*M [tauq, taup] + M*NB [work] |
|
* RWorkspace: need M [e] + M*M [RVT] + M*M [RU] |
* |
* |
CALL ZLACP2( 'F', M, M, RWORK( IRU ), M, U, LDU ) |
CALL ZLACP2( 'F', M, M, RWORK( IRU ), M, U, LDU ) |
CALL ZUNMBR( 'Q', 'L', 'N', M, M, N, A, LDA, |
CALL ZUNMBR( 'Q', 'L', 'N', M, M, N, A, LDA, |
Line 2017
|
Line 2178
|
* |
* |
* Copy real matrix RWORK(IRVT) to complex matrix VT |
* Copy real matrix RWORK(IRVT) to complex matrix VT |
* Overwrite VT by right singular vectors of A |
* Overwrite VT by right singular vectors of A |
* (CWorkspace: need 2*M+N, prefer 2*M+N*NB) |
* CWorkspace: need 2*M [tauq, taup] + N [work] |
* (RWorkspace: M*M) |
* CWorkspace: prefer 2*M [tauq, taup] + N*NB [work] |
|
* RWorkspace: need M [e] + M*M [RVT] |
* |
* |
CALL ZLACP2( 'F', M, M, RWORK( IRVT ), M, VT, LDVT ) |
CALL ZLACP2( 'F', M, M, RWORK( IRVT ), M, VT, LDVT ) |
CALL ZUNMBR( 'P', 'R', 'C', N, N, M, A, LDA, |
CALL ZUNMBR( 'P', 'R', 'C', N, N, M, A, LDA, |