version 1.1, 2010/01/26 15:22:46
|
version 1.17, 2017/06/17 10:54:11
|
Line 1
|
Line 1
|
SUBROUTINE ZGESDD( JOBZ, M, N, A, LDA, S, U, LDU, VT, LDVT, WORK, |
*> \brief \b ZGESDD |
$ LWORK, RWORK, IWORK, INFO ) |
|
* |
* |
* -- LAPACK driver routine (version 3.2) -- |
* =========== DOCUMENTATION =========== |
|
* |
|
* Online html documentation available at |
|
* http://www.netlib.org/lapack/explore-html/ |
|
* |
|
*> \htmlonly |
|
*> Download ZGESDD + dependencies |
|
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zgesdd.f"> |
|
*> [TGZ]</a> |
|
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zgesdd.f"> |
|
*> [ZIP]</a> |
|
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zgesdd.f"> |
|
*> [TXT]</a> |
|
*> \endhtmlonly |
|
* |
|
* Definition: |
|
* =========== |
|
* |
|
* SUBROUTINE ZGESDD( JOBZ, M, N, A, LDA, S, U, LDU, VT, LDVT, |
|
* WORK, LWORK, RWORK, IWORK, INFO ) |
|
* |
|
* .. Scalar Arguments .. |
|
* CHARACTER JOBZ |
|
* INTEGER INFO, LDA, LDU, LDVT, LWORK, M, N |
|
* .. |
|
* .. Array Arguments .. |
|
* INTEGER IWORK( * ) |
|
* DOUBLE PRECISION RWORK( * ), S( * ) |
|
* COMPLEX*16 A( LDA, * ), U( LDU, * ), VT( LDVT, * ), |
|
* $ WORK( * ) |
|
* .. |
|
* |
|
* |
|
*> \par Purpose: |
|
* ============= |
|
*> |
|
*> \verbatim |
|
*> |
|
*> ZGESDD computes the singular value decomposition (SVD) of a complex |
|
*> M-by-N matrix A, optionally computing the left and/or right singular |
|
*> vectors, by using divide-and-conquer method. The SVD is written |
|
*> |
|
*> A = U * SIGMA * conjugate-transpose(V) |
|
*> |
|
*> where SIGMA is an M-by-N matrix which is zero except for its |
|
*> min(m,n) diagonal elements, U is an M-by-M unitary matrix, and |
|
*> V is an N-by-N unitary matrix. The diagonal elements of SIGMA |
|
*> are the singular values of A; they are real and non-negative, and |
|
*> are returned in descending order. The first min(m,n) columns of |
|
*> U and V are the left and right singular vectors of A. |
|
*> |
|
*> Note that the routine returns VT = V**H, not V. |
|
*> |
|
*> The divide and conquer algorithm makes very mild assumptions about |
|
*> floating point arithmetic. It will work on machines with a guard |
|
*> digit in add/subtract, or on those binary machines without guard |
|
*> digits which subtract like the Cray X-MP, Cray Y-MP, Cray C-90, or |
|
*> Cray-2. It could conceivably fail on hexadecimal or decimal machines |
|
*> without guard digits, but we know of none. |
|
*> \endverbatim |
|
* |
|
* Arguments: |
|
* ========== |
|
* |
|
*> \param[in] JOBZ |
|
*> \verbatim |
|
*> JOBZ is CHARACTER*1 |
|
*> Specifies options for computing all or part of the matrix U: |
|
*> = 'A': all M columns of U and all N rows of V**H are |
|
*> returned in the arrays U and VT; |
|
*> = 'S': the first min(M,N) columns of U and the first |
|
*> min(M,N) rows of V**H are returned in the arrays U |
|
*> and VT; |
|
*> = 'O': If M >= N, the first N columns of U are overwritten |
|
*> in the array A and all rows of V**H are returned in |
|
*> the array VT; |
|
*> otherwise, all columns of U are returned in the |
|
*> array U and the first M rows of V**H are overwritten |
|
*> in the array A; |
|
*> = 'N': no columns of U or rows of V**H are computed. |
|
*> \endverbatim |
|
*> |
|
*> \param[in] M |
|
*> \verbatim |
|
*> M is INTEGER |
|
*> The number of rows of the input matrix A. M >= 0. |
|
*> \endverbatim |
|
*> |
|
*> \param[in] N |
|
*> \verbatim |
|
*> N is INTEGER |
|
*> The number of columns of the input matrix A. N >= 0. |
|
*> \endverbatim |
|
*> |
|
*> \param[in,out] A |
|
*> \verbatim |
|
*> A is COMPLEX*16 array, dimension (LDA,N) |
|
*> On entry, the M-by-N matrix A. |
|
*> On exit, |
|
*> if JOBZ = 'O', A is overwritten with the first N columns |
|
*> of U (the left singular vectors, stored |
|
*> columnwise) if M >= N; |
|
*> A is overwritten with the first M rows |
|
*> of V**H (the right singular vectors, stored |
|
*> rowwise) otherwise. |
|
*> if JOBZ .ne. 'O', the contents of A are destroyed. |
|
*> \endverbatim |
|
*> |
|
*> \param[in] LDA |
|
*> \verbatim |
|
*> LDA is INTEGER |
|
*> The leading dimension of the array A. LDA >= max(1,M). |
|
*> \endverbatim |
|
*> |
|
*> \param[out] S |
|
*> \verbatim |
|
*> S is DOUBLE PRECISION array, dimension (min(M,N)) |
|
*> The singular values of A, sorted so that S(i) >= S(i+1). |
|
*> \endverbatim |
|
*> |
|
*> \param[out] U |
|
*> \verbatim |
|
*> U is COMPLEX*16 array, dimension (LDU,UCOL) |
|
*> UCOL = M if JOBZ = 'A' or JOBZ = 'O' and M < N; |
|
*> UCOL = min(M,N) if JOBZ = 'S'. |
|
*> If JOBZ = 'A' or JOBZ = 'O' and M < N, U contains the M-by-M |
|
*> unitary matrix U; |
|
*> if JOBZ = 'S', U contains the first min(M,N) columns of U |
|
*> (the left singular vectors, stored columnwise); |
|
*> if JOBZ = 'O' and M >= N, or JOBZ = 'N', U is not referenced. |
|
*> \endverbatim |
|
*> |
|
*> \param[in] LDU |
|
*> \verbatim |
|
*> LDU is INTEGER |
|
*> The leading dimension of the array U. LDU >= 1; |
|
*> if JOBZ = 'S' or 'A' or JOBZ = 'O' and M < N, LDU >= M. |
|
*> \endverbatim |
|
*> |
|
*> \param[out] VT |
|
*> \verbatim |
|
*> VT is COMPLEX*16 array, dimension (LDVT,N) |
|
*> If JOBZ = 'A' or JOBZ = 'O' and M >= N, VT contains the |
|
*> N-by-N unitary matrix V**H; |
|
*> if JOBZ = 'S', VT contains the first min(M,N) rows of |
|
*> V**H (the right singular vectors, stored rowwise); |
|
*> if JOBZ = 'O' and M < N, or JOBZ = 'N', VT is not referenced. |
|
*> \endverbatim |
|
*> |
|
*> \param[in] LDVT |
|
*> \verbatim |
|
*> LDVT is INTEGER |
|
*> The leading dimension of the array VT. LDVT >= 1; |
|
*> if JOBZ = 'A' or JOBZ = 'O' and M >= N, LDVT >= N; |
|
*> if JOBZ = 'S', LDVT >= min(M,N). |
|
*> \endverbatim |
|
*> |
|
*> \param[out] WORK |
|
*> \verbatim |
|
*> WORK is COMPLEX*16 array, dimension (MAX(1,LWORK)) |
|
*> On exit, if INFO = 0, WORK(1) returns the optimal LWORK. |
|
*> \endverbatim |
|
*> |
|
*> \param[in] LWORK |
|
*> \verbatim |
|
*> LWORK is INTEGER |
|
*> The dimension of the array WORK. LWORK >= 1. |
|
*> If LWORK = -1, a workspace query is assumed. The optimal |
|
*> size for the WORK array is calculated and stored in WORK(1), |
|
*> 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 |
|
*> |
|
*> \param[out] RWORK |
|
*> \verbatim |
|
*> RWORK is DOUBLE PRECISION array, dimension (MAX(1,LRWORK)) |
|
*> Let mx = max(M,N) and mn = min(M,N). |
|
*> If JOBZ = 'N', LRWORK >= 5*mn (LAPACK <= 3.6 needs 7*mn); |
|
*> 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 |
|
*> |
|
*> \param[out] IWORK |
|
*> \verbatim |
|
*> IWORK is INTEGER array, dimension (8*min(M,N)) |
|
*> \endverbatim |
|
*> |
|
*> \param[out] INFO |
|
*> \verbatim |
|
*> INFO is INTEGER |
|
*> = 0: successful exit. |
|
*> < 0: if INFO = -i, the i-th argument had an illegal value. |
|
*> > 0: The updating process of DBDSDC did not converge. |
|
*> \endverbatim |
|
* |
|
* Authors: |
|
* ======== |
|
* |
|
*> \author Univ. of Tennessee |
|
*> \author Univ. of California Berkeley |
|
*> \author Univ. of Colorado Denver |
|
*> \author NAG Ltd. |
|
* |
|
*> \date June 2016 |
|
* |
|
*> \ingroup complex16GEsing |
|
* |
|
*> \par Contributors: |
|
* ================== |
|
*> |
|
*> Ming Gu and Huan Ren, Computer Science Division, University of |
|
*> California at Berkeley, USA |
|
*> |
|
* ===================================================================== |
|
SUBROUTINE ZGESDD( JOBZ, M, N, A, LDA, S, U, LDU, VT, LDVT, |
|
$ WORK, LWORK, RWORK, IWORK, INFO ) |
|
implicit none |
|
* |
|
* -- 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 2006 |
* June 2016 |
* 8-15-00: Improve consistency of WS calculations (eca) |
|
* |
* |
* .. Scalar Arguments .. |
* .. Scalar Arguments .. |
CHARACTER JOBZ |
CHARACTER JOBZ |
Line 18
|
Line 243
|
$ WORK( * ) |
$ WORK( * ) |
* .. |
* .. |
* |
* |
* Purpose |
|
* ======= |
|
* |
|
* ZGESDD computes the singular value decomposition (SVD) of a complex |
|
* M-by-N matrix A, optionally computing the left and/or right singular |
|
* vectors, by using divide-and-conquer method. The SVD is written |
|
* |
|
* A = U * SIGMA * conjugate-transpose(V) |
|
* |
|
* where SIGMA is an M-by-N matrix which is zero except for its |
|
* min(m,n) diagonal elements, U is an M-by-M unitary matrix, and |
|
* V is an N-by-N unitary matrix. The diagonal elements of SIGMA |
|
* are the singular values of A; they are real and non-negative, and |
|
* are returned in descending order. The first min(m,n) columns of |
|
* U and V are the left and right singular vectors of A. |
|
* |
|
* Note that the routine returns VT = V**H, not V. |
|
* |
|
* The divide and conquer algorithm makes very mild assumptions about |
|
* floating point arithmetic. It will work on machines with a guard |
|
* digit in add/subtract, or on those binary machines without guard |
|
* digits which subtract like the Cray X-MP, Cray Y-MP, Cray C-90, or |
|
* Cray-2. It could conceivably fail on hexadecimal or decimal machines |
|
* without guard digits, but we know of none. |
|
* |
|
* Arguments |
|
* ========= |
|
* |
|
* JOBZ (input) CHARACTER*1 |
|
* Specifies options for computing all or part of the matrix U: |
|
* = 'A': all M columns of U and all N rows of V**H are |
|
* returned in the arrays U and VT; |
|
* = 'S': the first min(M,N) columns of U and the first |
|
* min(M,N) rows of V**H are returned in the arrays U |
|
* and VT; |
|
* = 'O': If M >= N, the first N columns of U are overwritten |
|
* in the array A and all rows of V**H are returned in |
|
* the array VT; |
|
* otherwise, all columns of U are returned in the |
|
* array U and the first M rows of V**H are overwritten |
|
* in the array A; |
|
* = 'N': no columns of U or rows of V**H are computed. |
|
* |
|
* M (input) INTEGER |
|
* The number of rows of the input matrix A. M >= 0. |
|
* |
|
* N (input) INTEGER |
|
* The number of columns of the input matrix A. N >= 0. |
|
* |
|
* A (input/output) COMPLEX*16 array, dimension (LDA,N) |
|
* On entry, the M-by-N matrix A. |
|
* On exit, |
|
* if JOBZ = 'O', A is overwritten with the first N columns |
|
* of U (the left singular vectors, stored |
|
* columnwise) if M >= N; |
|
* A is overwritten with the first M rows |
|
* of V**H (the right singular vectors, stored |
|
* rowwise) otherwise. |
|
* if JOBZ .ne. 'O', the contents of A are destroyed. |
|
* |
|
* LDA (input) INTEGER |
|
* The leading dimension of the array A. LDA >= max(1,M). |
|
* |
|
* S (output) DOUBLE PRECISION array, dimension (min(M,N)) |
|
* The singular values of A, sorted so that S(i) >= S(i+1). |
|
* |
|
* U (output) COMPLEX*16 array, dimension (LDU,UCOL) |
|
* UCOL = M if JOBZ = 'A' or JOBZ = 'O' and M < N; |
|
* UCOL = min(M,N) if JOBZ = 'S'. |
|
* If JOBZ = 'A' or JOBZ = 'O' and M < N, U contains the M-by-M |
|
* unitary matrix U; |
|
* if JOBZ = 'S', U contains the first min(M,N) columns of U |
|
* (the left singular vectors, stored columnwise); |
|
* if JOBZ = 'O' and M >= N, or JOBZ = 'N', U is not referenced. |
|
* |
|
* LDU (input) INTEGER |
|
* The leading dimension of the array U. LDU >= 1; if |
|
* JOBZ = 'S' or 'A' or JOBZ = 'O' and M < N, LDU >= M. |
|
* |
|
* VT (output) COMPLEX*16 array, dimension (LDVT,N) |
|
* If JOBZ = 'A' or JOBZ = 'O' and M >= N, VT contains the |
|
* N-by-N unitary matrix V**H; |
|
* if JOBZ = 'S', VT contains the first min(M,N) rows of |
|
* V**H (the right singular vectors, stored rowwise); |
|
* if JOBZ = 'O' and M < N, or JOBZ = 'N', VT is not referenced. |
|
* |
|
* LDVT (input) INTEGER |
|
* The leading dimension of the array VT. LDVT >= 1; if |
|
* JOBZ = 'A' or JOBZ = 'O' and M >= N, LDVT >= N; |
|
* if JOBZ = 'S', LDVT >= min(M,N). |
|
* |
|
* WORK (workspace/output) COMPLEX*16 array, dimension (MAX(1,LWORK)) |
|
* On exit, if INFO = 0, WORK(1) returns the optimal LWORK. |
|
* |
|
* LWORK (input) INTEGER |
|
* 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 |
|
* size for the WORK array is calculated and stored in WORK(1), |
|
* and no other work except argument checking is performed. |
|
* |
|
* RWORK (workspace) DOUBLE PRECISION array, dimension (MAX(1,LRWORK)) |
|
* If JOBZ = 'N', LRWORK >= 5*min(M,N). |
|
* Otherwise, LRWORK >= 5*min(M,N)*min(M,N) + 7*min(M,N) |
|
* |
|
* IWORK (workspace) INTEGER array, dimension (8*min(M,N)) |
|
* |
|
* INFO (output) INTEGER |
|
* = 0: successful exit. |
|
* < 0: if INFO = -i, the i-th argument had an illegal value. |
|
* > 0: The updating process of DBDSDC did not converge. |
|
* |
|
* Further Details |
|
* =============== |
|
* |
|
* Based on contributions by |
|
* Ming Gu and Huan Ren, Computer Science Division, University of |
|
* California at Berkeley, USA |
|
* |
|
* ===================================================================== |
* ===================================================================== |
* |
* |
* .. 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 155
|
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 173
|
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 184
|
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 214
|
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 459
|
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 503
|
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 526
|
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 535
|
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 553
|
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 578
|
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 589
|
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 599
|
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 611
|
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 622
|
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 632
|
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 646
|
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 659
|
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 672
|
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 683
|
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 693
|
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 705
|
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 715
|
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 725
|
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 734
|
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 747
|
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 771
|
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 784
|
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 793
|
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 804
|
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 814
|
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 830
|
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 841
|
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 861
|
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 885
|
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 901
|
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 910
|
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 924
|
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 943
|
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 955
|
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 964
|
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 973
|
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 992
|
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 1004
|
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 1013
|
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 1026
|
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 1037
|
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 1065
|
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 1081
|
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 1107
|
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 1132
|
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 1147
|
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 1158
|
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 1167
|
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 1190
|
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 1200
|
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 1221
|
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 1244
|
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 1253
|
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 1271
|
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 1282
|
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 1301
|
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 1312
|
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 1322
|
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 1334
|
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 1344
|
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 1355
|
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 1369
|
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 1382
|
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 1395
|
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 1406
|
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 1416
|
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 1428
|
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 1438
|
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 1448
|
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 1457
|
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 1470
|
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 1494
|
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 1504
|
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 1516
|
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 1526
|
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 1537
|
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 1553
|
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 1565
|
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 1574
|
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 1586
|
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 1612
|
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 1628
|
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 1637
|
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 1650
|
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 1669
|
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 1681
|
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 1690
|
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 1699
|
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 1718
|
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 1730
|
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 1739
|
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 1751
|
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 1762
|
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 1790
|
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 1809
|
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 1833
|
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 1857
|
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 1872
|
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 1882
|
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 1892
|
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 1908
|
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 1922
|
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, |