version 1.7, 2010/08/13 21:04:03
|
version 1.20, 2023/08/07 08:39:19
|
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.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: if INFO = -i, the i-th argument had an illegal value. |
|
*> = -4: if A had a NAN entry. |
|
*> > 0: The updating process of DBDSDC did not converge. |
|
*> = 0: successful exit. |
|
*> \endverbatim |
|
* |
|
* Authors: |
|
* ======== |
|
* |
|
*> \author Univ. of Tennessee |
|
*> \author Univ. of California Berkeley |
|
*> \author Univ. of Colorado Denver |
|
*> \author NAG Ltd. |
|
* |
|
*> \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 -- |
* -- 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..-- |
* June 2010 |
|
* 8-15-00: Improve consistency of WS calculations (eca) |
|
* |
* |
* .. Scalar Arguments .. |
* .. Scalar Arguments .. |
CHARACTER JOBZ |
CHARACTER JOBZ |
Line 18
|
Line 241
|
$ 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 >= min(M,N)*max(5*min(M,N)+7,2*max(M,N)+2*min(M,N)+1) |
|
* |
|
* 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 156
|
Line 251
|
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 279
|
$ ZLASET, ZUNGBR, ZUNGLQ, ZUNGQR, ZUNMBR |
$ ZLASET, ZUNGBR, ZUNGLQ, ZUNGQR, ZUNMBR |
* .. |
* .. |
* .. External Functions .. |
* .. External Functions .. |
LOGICAL LSAME |
LOGICAL LSAME, DISNAN |
INTEGER ILAENV |
DOUBLE PRECISION DLAMCH, ZLANGE, DROUNDUP_LWORK |
DOUBLE PRECISION DLAMCH, ZLANGE |
EXTERNAL LSAME, DLAMCH, ZLANGE, DISNAN, |
EXTERNAL LSAME, ILAENV, DLAMCH, ZLANGE |
$ DROUNDUP_LWORK |
* .. |
* .. |
* .. Intrinsic Functions .. |
* .. Intrinsic Functions .. |
INTRINSIC INT, MAX, MIN, SQRT |
INTRINSIC INT, MAX, MIN, SQRT |
Line 185
|
Line 291
|
* |
* |
* 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 215
|
Line 322
|
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 |
MAXWRK = MAX( MAXWRK, MINWRK ) |
MAXWRK = MAX( MAXWRK, MINWRK ) |
END IF |
END IF |
IF( INFO.EQ.0 ) THEN |
IF( INFO.EQ.0 ) THEN |
WORK( 1 ) = MAXWRK |
WORK( 1 ) = DROUNDUP_LWORK( 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 485
|
Line 646
|
* Scale A if max element outside range [SMLNUM,BIGNUM] |
* Scale A if max element outside range [SMLNUM,BIGNUM] |
* |
* |
ANRM = ZLANGE( 'M', M, N, A, LDA, DUM ) |
ANRM = ZLANGE( 'M', M, N, A, LDA, DUM ) |
|
IF( DISNAN( ANRM ) ) THEN |
|
INFO = -4 |
|
RETURN |
|
END IF |
ISCL = 0 |
ISCL = 0 |
IF( ANRM.GT.ZERO .AND. ANRM.LT.SMLNUM ) THEN |
IF( ANRM.GT.ZERO .AND. ANRM.LT.SMLNUM ) THEN |
ISCL = 1 |
ISCL = 1 |
Line 504
|
Line 669
|
* |
* |
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 527
|
Line 693
|
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 536
|
Line 703
|
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 554
|
Line 721
|
* |
* |
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 579
|
Line 747
|
$ 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 590
|
Line 759
|
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 600
|
Line 770
|
* 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 612
|
Line 782
|
* |
* |
* 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 623
|
Line 794
|
* |
* |
* 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 633
|
Line 805
|
* |
* |
* 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 647
|
Line 820
|
* |
* |
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 660
|
Line 833
|
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 673
|
Line 847
|
$ 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 684
|
Line 859
|
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 694
|
Line 870
|
* 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 706
|
Line 882
|
* |
* |
* 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 716
|
Line 893
|
* |
* |
* 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 726
|
Line 904
|
* |
* |
* 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 735
|
Line 913
|
* |
* |
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 748
|
Line 926
|
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 772
|
Line 952
|
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 785
|
Line 966
|
* 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 794
|
Line 975
|
* |
* |
* 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 805
|
Line 987
|
* |
* |
* 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 815
|
Line 998
|
* |
* |
* 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 831
|
Line 1014
|
* |
* |
* 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 842
|
Line 1025
|
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 862
|
Line 1047
|
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 886
|
Line 1074
|
* |
* |
* 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 902
|
Line 1090
|
* |
* |
* 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 911
|
Line 1099
|
* |
* |
* 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 925
|
Line 1115
|
* |
* |
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 944
|
Line 1137
|
* 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 956
|
Line 1149
|
* |
* |
* 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 965
|
Line 1158
|
* |
* |
* 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 974
|
Line 1167
|
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 993
|
Line 1189
|
* 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 1005
|
Line 1201
|
* |
* |
* 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 1014
|
Line 1210
|
* |
* |
* 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 1027
|
Line 1223
|
* |
* |
* 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 1038
|
Line 1234
|
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 1066
|
Line 1264
|
* |
* |
* 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 1082
|
Line 1281
|
* |
* |
* 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 1108
|
Line 1310
|
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 1133
|
Line 1339
|
* |
* |
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 1148
|
Line 1355
|
* |
* |
* 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 1159
|
Line 1367
|
* |
* |
* 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 1168
|
Line 1377
|
$ 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 1191
|
Line 1401
|
* |
* |
* 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 1201
|
Line 1412
|
* |
* |
* 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 1222
|
Line 1434
|
* |
* |
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 1245
|
Line 1458
|
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 1254
|
Line 1468
|
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 1272
|
Line 1486
|
* 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 1283
|
Line 1497
|
* 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 1302
|
Line 1517
|
$ 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 1313
|
Line 1529
|
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 1323
|
Line 1540
|
* 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 1335
|
Line 1552
|
* |
* |
* 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 1345
|
Line 1563
|
* |
* |
* 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 1356
|
Line 1575
|
* |
* |
* 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 1370
|
Line 1590
|
* |
* |
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 1383
|
Line 1603
|
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 1396
|
Line 1617
|
$ 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 1407
|
Line 1629
|
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 1417
|
Line 1640
|
* Perform bidiagonal SVD, computing left singular vectors |
* Perform bidiagonal SVD, computing left singular vectors |
* of bidiagonal matrix in RWORK(IRU) and computing right |
* of bidiagonal matrix in RWORK(IRU) and computing right |
* singular vectors of bidiagonal matrix in RWORK(IRVT) |
* singular vectors of bidiagonal matrix in RWORK(IRVT) |
* (CWorkspace: need 0) |
* CWorkspace: need 0 |
* (RWorkspace: need BDSPAC) |
* RWorkspace: need M [e] + M*M [RU] + M*M [RVT] + BDSPAC |
* |
* |
IRU = IE + M |
IRU = IE + M |
IRVT = IRU + M*M |
IRVT = IRU + M*M |
Line 1429
|
Line 1652
|
* |
* |
* 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 1439
|
Line 1663
|
* |
* |
* 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 1449
|
Line 1674
|
* |
* |
* 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 1458
|
Line 1683
|
* |
* |
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 1471
|
Line 1696
|
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 1495
|
Line 1722
|
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 1505
|
Line 1733
|
* 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 1517
|
Line 1745
|
* |
* |
* 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 1527
|
Line 1756
|
* |
* |
* 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 1538
|
Line 1768
|
* |
* |
* 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 1554
|
Line 1784
|
* |
* |
* 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 1566
|
Line 1795
|
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 1575
|
Line 1805
|
* |
* |
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 1587
|
Line 1818
|
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 1613
|
Line 1847
|
* |
* |
* 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 1629
|
Line 1863
|
* |
* |
* 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 1638
|
Line 1872
|
* |
* |
* 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 1651
|
Line 1887
|
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 1670
|
Line 1909
|
* 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 1682
|
Line 1921
|
* |
* |
* 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 1691
|
Line 1930
|
* |
* |
* 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 1700
|
Line 1939
|
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 1719
|
Line 1961
|
* 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 1731
|
Line 1973
|
* |
* |
* 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 1740
|
Line 1982
|
* |
* |
* 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 1752
|
Line 1995
|
* |
* |
* 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 1763
|
Line 2006
|
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 1791
|
Line 2037
|
* |
* |
* 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 1810
|
Line 2056
|
* |
* |
* 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 1834
|
Line 2083
|
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 1858
|
Line 2111
|
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 1873
|
Line 2127
|
* |
* |
* 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 1883
|
Line 2138
|
* |
* |
* 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 1893
|
Line 2149
|
$ 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 1909
|
Line 2166
|
* |
* |
* 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 1923
|
Line 2181
|
* |
* |
* 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, |
Line 1955
|
Line 2214
|
* |
* |
* Return optimal workspace in WORK(1) |
* Return optimal workspace in WORK(1) |
* |
* |
WORK( 1 ) = MAXWRK |
WORK( 1 ) = DROUNDUP_LWORK( MAXWRK ) |
* |
* |
RETURN |
RETURN |
* |
* |