version 1.6, 2010/08/13 21:03:45
|
version 1.16, 2016/08/27 15:34:22
|
Line 1
|
Line 1
|
SUBROUTINE DGESVD( JOBU, JOBVT, M, N, A, LDA, S, U, LDU, VT, LDVT, |
*> \brief <b> DGESVD computes the singular value decomposition (SVD) for GE matrices</b> |
$ WORK, LWORK, INFO ) |
|
* |
* |
* -- LAPACK driver routine (version 3.2) -- |
* =========== DOCUMENTATION =========== |
|
* |
|
* Online html documentation available at |
|
* http://www.netlib.org/lapack/explore-html/ |
|
* |
|
*> \htmlonly |
|
*> Download DGESVD + dependencies |
|
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dgesvd.f"> |
|
*> [TGZ]</a> |
|
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dgesvd.f"> |
|
*> [ZIP]</a> |
|
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dgesvd.f"> |
|
*> [TXT]</a> |
|
*> \endhtmlonly |
|
* |
|
* Definition: |
|
* =========== |
|
* |
|
* SUBROUTINE DGESVD( JOBU, JOBVT, M, N, A, LDA, S, U, LDU, VT, LDVT, |
|
* WORK, LWORK, INFO ) |
|
* |
|
* .. Scalar Arguments .. |
|
* CHARACTER JOBU, JOBVT |
|
* INTEGER INFO, LDA, LDU, LDVT, LWORK, M, N |
|
* .. |
|
* .. Array Arguments .. |
|
* DOUBLE PRECISION A( LDA, * ), S( * ), U( LDU, * ), |
|
* $ VT( LDVT, * ), WORK( * ) |
|
* .. |
|
* |
|
* |
|
*> \par Purpose: |
|
* ============= |
|
*> |
|
*> \verbatim |
|
*> |
|
*> DGESVD computes the singular value decomposition (SVD) of a real |
|
*> M-by-N matrix A, optionally computing the left and/or right singular |
|
*> vectors. The SVD is written |
|
*> |
|
*> A = U * SIGMA * transpose(V) |
|
*> |
|
*> where SIGMA is an M-by-N matrix which is zero except for its |
|
*> min(m,n) diagonal elements, U is an M-by-M orthogonal matrix, and |
|
*> V is an N-by-N orthogonal 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 V**T, not V. |
|
*> \endverbatim |
|
* |
|
* Arguments: |
|
* ========== |
|
* |
|
*> \param[in] JOBU |
|
*> \verbatim |
|
*> JOBU is CHARACTER*1 |
|
*> Specifies options for computing all or part of the matrix U: |
|
*> = 'A': all M columns of U are returned in array U: |
|
*> = 'S': the first min(m,n) columns of U (the left singular |
|
*> vectors) are returned in the array U; |
|
*> = 'O': the first min(m,n) columns of U (the left singular |
|
*> vectors) are overwritten on the array A; |
|
*> = 'N': no columns of U (no left singular vectors) are |
|
*> computed. |
|
*> \endverbatim |
|
*> |
|
*> \param[in] JOBVT |
|
*> \verbatim |
|
*> JOBVT is CHARACTER*1 |
|
*> Specifies options for computing all or part of the matrix |
|
*> V**T: |
|
*> = 'A': all N rows of V**T are returned in the array VT; |
|
*> = 'S': the first min(m,n) rows of V**T (the right singular |
|
*> vectors) are returned in the array VT; |
|
*> = 'O': the first min(m,n) rows of V**T (the right singular |
|
*> vectors) are overwritten on the array A; |
|
*> = 'N': no rows of V**T (no right singular vectors) are |
|
*> computed. |
|
*> |
|
*> JOBVT and JOBU cannot both be 'O'. |
|
*> \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 DOUBLE PRECISION array, dimension (LDA,N) |
|
*> On entry, the M-by-N matrix A. |
|
*> On exit, |
|
*> if JOBU = 'O', A is overwritten with the first min(m,n) |
|
*> columns of U (the left singular vectors, |
|
*> stored columnwise); |
|
*> if JOBVT = 'O', A is overwritten with the first min(m,n) |
|
*> rows of V**T (the right singular vectors, |
|
*> stored rowwise); |
|
*> if JOBU .ne. 'O' and JOBVT .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 DOUBLE PRECISION array, dimension (LDU,UCOL) |
|
*> (LDU,M) if JOBU = 'A' or (LDU,min(M,N)) if JOBU = 'S'. |
|
*> If JOBU = 'A', U contains the M-by-M orthogonal matrix U; |
|
*> if JOBU = 'S', U contains the first min(m,n) columns of U |
|
*> (the left singular vectors, stored columnwise); |
|
*> if JOBU = 'N' or 'O', U is not referenced. |
|
*> \endverbatim |
|
*> |
|
*> \param[in] LDU |
|
*> \verbatim |
|
*> LDU is INTEGER |
|
*> The leading dimension of the array U. LDU >= 1; if |
|
*> JOBU = 'S' or 'A', LDU >= M. |
|
*> \endverbatim |
|
*> |
|
*> \param[out] VT |
|
*> \verbatim |
|
*> VT is DOUBLE PRECISION array, dimension (LDVT,N) |
|
*> If JOBVT = 'A', VT contains the N-by-N orthogonal matrix |
|
*> V**T; |
|
*> if JOBVT = 'S', VT contains the first min(m,n) rows of |
|
*> V**T (the right singular vectors, stored rowwise); |
|
*> if JOBVT = 'N' or 'O', VT is not referenced. |
|
*> \endverbatim |
|
*> |
|
*> \param[in] LDVT |
|
*> \verbatim |
|
*> LDVT is INTEGER |
|
*> The leading dimension of the array VT. LDVT >= 1; if |
|
*> JOBVT = 'A', LDVT >= N; if JOBVT = 'S', LDVT >= min(M,N). |
|
*> \endverbatim |
|
*> |
|
*> \param[out] WORK |
|
*> \verbatim |
|
*> WORK is DOUBLE PRECISION array, dimension (MAX(1,LWORK)) |
|
*> On exit, if INFO = 0, WORK(1) returns the optimal LWORK; |
|
*> if INFO > 0, WORK(2:MIN(M,N)) contains the unconverged |
|
*> superdiagonal elements of an upper bidiagonal matrix B |
|
*> whose diagonal is in S (not necessarily sorted). B |
|
*> satisfies A = U * B * VT, so it has the same singular values |
|
*> as A, and singular vectors related by U and VT. |
|
*> \endverbatim |
|
*> |
|
*> \param[in] LWORK |
|
*> \verbatim |
|
*> LWORK is INTEGER |
|
*> The dimension of the array WORK. |
|
*> LWORK >= MAX(1,5*MIN(M,N)) for the paths (see comments inside code): |
|
*> - PATH 1 (M much larger than N, JOBU='N') |
|
*> - PATH 1t (N much larger than M, JOBVT='N') |
|
*> LWORK >= MAX(1,3*MIN(M,N) + MAX(M,N),5*MIN(M,N)) for the other paths |
|
*> For good performance, LWORK should generally be larger. |
|
*> |
|
*> If LWORK = -1, then a workspace query is assumed; the routine |
|
*> only calculates the optimal size of the WORK array, returns |
|
*> this value as the first entry of the WORK array, and no error |
|
*> message related to LWORK is issued by XERBLA. |
|
*> \endverbatim |
|
*> |
|
*> \param[out] INFO |
|
*> \verbatim |
|
*> INFO is INTEGER |
|
*> = 0: successful exit. |
|
*> < 0: if INFO = -i, the i-th argument had an illegal value. |
|
*> > 0: if DBDSQR did not converge, INFO specifies how many |
|
*> superdiagonals of an intermediate bidiagonal form B |
|
*> did not converge to zero. See the description of WORK |
|
*> above for details. |
|
*> \endverbatim |
|
* |
|
* Authors: |
|
* ======== |
|
* |
|
*> \author Univ. of Tennessee |
|
*> \author Univ. of California Berkeley |
|
*> \author Univ. of Colorado Denver |
|
*> \author NAG Ltd. |
|
* |
|
*> \date April 2012 |
|
* |
|
*> \ingroup doubleGEsing |
|
* |
|
* ===================================================================== |
|
SUBROUTINE DGESVD( JOBU, JOBVT, M, N, A, LDA, S, U, LDU, |
|
$ VT, LDVT, WORK, LWORK, INFO ) |
|
* |
|
* -- LAPACK driver routine (version 3.6.1) -- |
* -- LAPACK is a software package provided by Univ. of Tennessee, -- |
* -- LAPACK is a software package provided by Univ. of Tennessee, -- |
* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- |
* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- |
* November 2006 |
* April 2012 |
* |
* |
* .. Scalar Arguments .. |
* .. Scalar Arguments .. |
CHARACTER JOBU, JOBVT |
CHARACTER JOBU, JOBVT |
Line 15
|
Line 225
|
$ VT( LDVT, * ), WORK( * ) |
$ VT( LDVT, * ), WORK( * ) |
* .. |
* .. |
* |
* |
* Purpose |
|
* ======= |
|
* |
|
* DGESVD computes the singular value decomposition (SVD) of a real |
|
* M-by-N matrix A, optionally computing the left and/or right singular |
|
* vectors. The SVD is written |
|
* |
|
* A = U * SIGMA * transpose(V) |
|
* |
|
* where SIGMA is an M-by-N matrix which is zero except for its |
|
* min(m,n) diagonal elements, U is an M-by-M orthogonal matrix, and |
|
* V is an N-by-N orthogonal 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 V**T, not V. |
|
* |
|
* Arguments |
|
* ========= |
|
* |
|
* JOBU (input) CHARACTER*1 |
|
* Specifies options for computing all or part of the matrix U: |
|
* = 'A': all M columns of U are returned in array U: |
|
* = 'S': the first min(m,n) columns of U (the left singular |
|
* vectors) are returned in the array U; |
|
* = 'O': the first min(m,n) columns of U (the left singular |
|
* vectors) are overwritten on the array A; |
|
* = 'N': no columns of U (no left singular vectors) are |
|
* computed. |
|
* |
|
* JOBVT (input) CHARACTER*1 |
|
* Specifies options for computing all or part of the matrix |
|
* V**T: |
|
* = 'A': all N rows of V**T are returned in the array VT; |
|
* = 'S': the first min(m,n) rows of V**T (the right singular |
|
* vectors) are returned in the array VT; |
|
* = 'O': the first min(m,n) rows of V**T (the right singular |
|
* vectors) are overwritten on the array A; |
|
* = 'N': no rows of V**T (no right singular vectors) are |
|
* computed. |
|
* |
|
* JOBVT and JOBU cannot both be 'O'. |
|
* |
|
* 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) DOUBLE PRECISION array, dimension (LDA,N) |
|
* On entry, the M-by-N matrix A. |
|
* On exit, |
|
* if JOBU = 'O', A is overwritten with the first min(m,n) |
|
* columns of U (the left singular vectors, |
|
* stored columnwise); |
|
* if JOBVT = 'O', A is overwritten with the first min(m,n) |
|
* rows of V**T (the right singular vectors, |
|
* stored rowwise); |
|
* if JOBU .ne. 'O' and JOBVT .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) DOUBLE PRECISION array, dimension (LDU,UCOL) |
|
* (LDU,M) if JOBU = 'A' or (LDU,min(M,N)) if JOBU = 'S'. |
|
* If JOBU = 'A', U contains the M-by-M orthogonal matrix U; |
|
* if JOBU = 'S', U contains the first min(m,n) columns of U |
|
* (the left singular vectors, stored columnwise); |
|
* if JOBU = 'N' or 'O', U is not referenced. |
|
* |
|
* LDU (input) INTEGER |
|
* The leading dimension of the array U. LDU >= 1; if |
|
* JOBU = 'S' or 'A', LDU >= M. |
|
* |
|
* VT (output) DOUBLE PRECISION array, dimension (LDVT,N) |
|
* If JOBVT = 'A', VT contains the N-by-N orthogonal matrix |
|
* V**T; |
|
* if JOBVT = 'S', VT contains the first min(m,n) rows of |
|
* V**T (the right singular vectors, stored rowwise); |
|
* if JOBVT = 'N' or 'O', VT is not referenced. |
|
* |
|
* LDVT (input) INTEGER |
|
* The leading dimension of the array VT. LDVT >= 1; if |
|
* JOBVT = 'A', LDVT >= N; if JOBVT = 'S', LDVT >= min(M,N). |
|
* |
|
* WORK (workspace/output) DOUBLE PRECISION array, dimension (MAX(1,LWORK)) |
|
* On exit, if INFO = 0, WORK(1) returns the optimal LWORK; |
|
* if INFO > 0, WORK(2:MIN(M,N)) contains the unconverged |
|
* superdiagonal elements of an upper bidiagonal matrix B |
|
* whose diagonal is in S (not necessarily sorted). B |
|
* satisfies A = U * B * VT, so it has the same singular values |
|
* as A, and singular vectors related by U and VT. |
|
* |
|
* LWORK (input) INTEGER |
|
* The dimension of the array WORK. |
|
* LWORK >= MAX(1,3*MIN(M,N)+MAX(M,N),5*MIN(M,N)). |
|
* For good performance, LWORK should generally be larger. |
|
* |
|
* If LWORK = -1, then a workspace query is assumed; the routine |
|
* only calculates the optimal size of the WORK array, returns |
|
* this value as the first entry of the WORK array, and no error |
|
* message related to LWORK is issued by XERBLA. |
|
* |
|
* INFO (output) INTEGER |
|
* = 0: successful exit. |
|
* < 0: if INFO = -i, the i-th argument had an illegal value. |
|
* > 0: if DBDSQR did not converge, INFO specifies how many |
|
* superdiagonals of an intermediate bidiagonal form B |
|
* did not converge to zero. See the description of WORK |
|
* above for details. |
|
* |
|
* ===================================================================== |
* ===================================================================== |
* |
* |
* .. Parameters .. |
* .. Parameters .. |
Line 144
|
Line 238
|
$ ITAU, ITAUP, ITAUQ, IU, IWORK, LDWRKR, LDWRKU, |
$ ITAU, ITAUP, ITAUQ, IU, IWORK, LDWRKR, LDWRKU, |
$ MAXWRK, MINMN, MINWRK, MNTHR, NCU, NCVT, NRU, |
$ MAXWRK, MINMN, MINWRK, MNTHR, NCU, NCVT, NRU, |
$ NRVT, WRKBL |
$ NRVT, WRKBL |
|
INTEGER LWORK_DGEQRF, LWORK_DORGQR_N, LWORK_DORGQR_M, |
|
$ LWORK_DGEBRD, LWORK_DORGBR_P, LWORK_DORGBR_Q, |
|
$ LWORK_DGELQF, LWORK_DORGLQ_N, LWORK_DORGLQ_M |
DOUBLE PRECISION ANRM, BIGNUM, EPS, SMLNUM |
DOUBLE PRECISION ANRM, BIGNUM, EPS, SMLNUM |
* .. |
* .. |
* .. Local Arrays .. |
* .. Local Arrays .. |
Line 215
|
Line 312
|
* |
* |
MNTHR = ILAENV( 6, 'DGESVD', JOBU // JOBVT, M, N, 0, 0 ) |
MNTHR = ILAENV( 6, 'DGESVD', JOBU // JOBVT, M, N, 0, 0 ) |
BDSPAC = 5*N |
BDSPAC = 5*N |
|
* Compute space needed for DGEQRF |
|
CALL DGEQRF( M, N, A, LDA, DUM(1), DUM(1), -1, IERR ) |
|
LWORK_DGEQRF = INT( DUM(1) ) |
|
* Compute space needed for DORGQR |
|
CALL DORGQR( M, N, N, A, LDA, DUM(1), DUM(1), -1, IERR ) |
|
LWORK_DORGQR_N = INT( DUM(1) ) |
|
CALL DORGQR( M, M, N, A, LDA, DUM(1), DUM(1), -1, IERR ) |
|
LWORK_DORGQR_M = INT( DUM(1) ) |
|
* Compute space needed for DGEBRD |
|
CALL DGEBRD( N, N, A, LDA, S, DUM(1), DUM(1), |
|
$ DUM(1), DUM(1), -1, IERR ) |
|
LWORK_DGEBRD = INT( DUM(1) ) |
|
* Compute space needed for DORGBR P |
|
CALL DORGBR( 'P', N, N, N, A, LDA, DUM(1), |
|
$ DUM(1), -1, IERR ) |
|
LWORK_DORGBR_P = INT( DUM(1) ) |
|
* Compute space needed for DORGBR Q |
|
CALL DORGBR( 'Q', N, N, N, A, LDA, DUM(1), |
|
$ DUM(1), -1, IERR ) |
|
LWORK_DORGBR_Q = INT( DUM(1) ) |
|
* |
IF( M.GE.MNTHR ) THEN |
IF( M.GE.MNTHR ) THEN |
IF( WNTUN ) THEN |
IF( WNTUN ) THEN |
* |
* |
* Path 1 (M much larger than N, JOBU='N') |
* Path 1 (M much larger than N, JOBU='N') |
* |
* |
MAXWRK = N + N*ILAENV( 1, 'DGEQRF', ' ', M, N, -1, |
MAXWRK = N + LWORK_DGEQRF |
$ -1 ) |
MAXWRK = MAX( MAXWRK, 3*N + LWORK_DGEBRD ) |
MAXWRK = MAX( MAXWRK, 3*N+2*N* |
|
$ ILAENV( 1, 'DGEBRD', ' ', N, N, -1, -1 ) ) |
|
IF( WNTVO .OR. WNTVAS ) |
IF( WNTVO .OR. WNTVAS ) |
$ MAXWRK = MAX( MAXWRK, 3*N+( N-1 )* |
$ MAXWRK = MAX( MAXWRK, 3*N + LWORK_DORGBR_P ) |
$ ILAENV( 1, 'DORGBR', 'P', N, N, N, -1 ) ) |
|
MAXWRK = MAX( MAXWRK, BDSPAC ) |
MAXWRK = MAX( MAXWRK, BDSPAC ) |
MINWRK = MAX( 4*N, BDSPAC ) |
MINWRK = MAX( 4*N, BDSPAC ) |
ELSE IF( WNTUO .AND. WNTVN ) THEN |
ELSE IF( WNTUO .AND. WNTVN ) THEN |
* |
* |
* Path 2 (M much larger than N, JOBU='O', JOBVT='N') |
* Path 2 (M much larger than N, JOBU='O', JOBVT='N') |
* |
* |
WRKBL = N + N*ILAENV( 1, 'DGEQRF', ' ', M, N, -1, -1 ) |
WRKBL = N + LWORK_DGEQRF |
WRKBL = MAX( WRKBL, N+N*ILAENV( 1, 'DORGQR', ' ', M, |
WRKBL = MAX( WRKBL, N + LWORK_DORGQR_N ) |
$ N, N, -1 ) ) |
WRKBL = MAX( WRKBL, 3*N + LWORK_DGEBRD ) |
WRKBL = MAX( WRKBL, 3*N+2*N* |
WRKBL = MAX( WRKBL, 3*N + LWORK_DORGBR_Q ) |
$ ILAENV( 1, 'DGEBRD', ' ', N, N, -1, -1 ) ) |
|
WRKBL = MAX( WRKBL, 3*N+N* |
|
$ ILAENV( 1, 'DORGBR', 'Q', N, N, N, -1 ) ) |
|
WRKBL = MAX( WRKBL, BDSPAC ) |
WRKBL = MAX( WRKBL, BDSPAC ) |
MAXWRK = MAX( N*N+WRKBL, N*N+M*N+N ) |
MAXWRK = MAX( N*N + WRKBL, N*N + M*N + N ) |
MINWRK = MAX( 3*N+M, BDSPAC ) |
MINWRK = MAX( 3*N + M, BDSPAC ) |
ELSE IF( WNTUO .AND. WNTVAS ) THEN |
ELSE IF( WNTUO .AND. WNTVAS ) THEN |
* |
* |
* Path 3 (M much larger than N, JOBU='O', JOBVT='S' or |
* Path 3 (M much larger than N, JOBU='O', JOBVT='S' or |
* 'A') |
* 'A') |
* |
* |
WRKBL = N + N*ILAENV( 1, 'DGEQRF', ' ', M, N, -1, -1 ) |
WRKBL = N + LWORK_DGEQRF |
WRKBL = MAX( WRKBL, N+N*ILAENV( 1, 'DORGQR', ' ', M, |
WRKBL = MAX( WRKBL, N + LWORK_DORGQR_N ) |
$ N, N, -1 ) ) |
WRKBL = MAX( WRKBL, 3*N + LWORK_DGEBRD ) |
WRKBL = MAX( WRKBL, 3*N+2*N* |
WRKBL = MAX( WRKBL, 3*N + LWORK_DORGBR_Q ) |
$ ILAENV( 1, 'DGEBRD', ' ', N, N, -1, -1 ) ) |
WRKBL = MAX( WRKBL, 3*N + LWORK_DORGBR_P ) |
WRKBL = MAX( WRKBL, 3*N+N* |
|
$ ILAENV( 1, 'DORGBR', 'Q', N, N, N, -1 ) ) |
|
WRKBL = MAX( WRKBL, 3*N+( N-1 )* |
|
$ ILAENV( 1, 'DORGBR', 'P', N, N, N, -1 ) ) |
|
WRKBL = MAX( WRKBL, BDSPAC ) |
WRKBL = MAX( WRKBL, BDSPAC ) |
MAXWRK = MAX( N*N+WRKBL, N*N+M*N+N ) |
MAXWRK = MAX( N*N + WRKBL, N*N + M*N + N ) |
MINWRK = MAX( 3*N+M, BDSPAC ) |
MINWRK = MAX( 3*N + M, BDSPAC ) |
ELSE IF( WNTUS .AND. WNTVN ) THEN |
ELSE IF( WNTUS .AND. WNTVN ) THEN |
* |
* |
* Path 4 (M much larger than N, JOBU='S', JOBVT='N') |
* Path 4 (M much larger than N, JOBU='S', JOBVT='N') |
* |
* |
WRKBL = N + N*ILAENV( 1, 'DGEQRF', ' ', M, N, -1, -1 ) |
WRKBL = N + LWORK_DGEQRF |
WRKBL = MAX( WRKBL, N+N*ILAENV( 1, 'DORGQR', ' ', M, |
WRKBL = MAX( WRKBL, N + LWORK_DORGQR_N ) |
$ N, N, -1 ) ) |
WRKBL = MAX( WRKBL, 3*N + LWORK_DGEBRD ) |
WRKBL = MAX( WRKBL, 3*N+2*N* |
WRKBL = MAX( WRKBL, 3*N + LWORK_DORGBR_Q ) |
$ ILAENV( 1, 'DGEBRD', ' ', N, N, -1, -1 ) ) |
|
WRKBL = MAX( WRKBL, 3*N+N* |
|
$ ILAENV( 1, 'DORGBR', 'Q', N, N, N, -1 ) ) |
|
WRKBL = MAX( WRKBL, BDSPAC ) |
WRKBL = MAX( WRKBL, BDSPAC ) |
MAXWRK = N*N + WRKBL |
MAXWRK = N*N + WRKBL |
MINWRK = MAX( 3*N+M, BDSPAC ) |
MINWRK = MAX( 3*N + M, BDSPAC ) |
ELSE IF( WNTUS .AND. WNTVO ) THEN |
ELSE IF( WNTUS .AND. WNTVO ) THEN |
* |
* |
* Path 5 (M much larger than N, JOBU='S', JOBVT='O') |
* Path 5 (M much larger than N, JOBU='S', JOBVT='O') |
* |
* |
WRKBL = N + N*ILAENV( 1, 'DGEQRF', ' ', M, N, -1, -1 ) |
WRKBL = N + LWORK_DGEQRF |
WRKBL = MAX( WRKBL, N+N*ILAENV( 1, 'DORGQR', ' ', M, |
WRKBL = MAX( WRKBL, N + LWORK_DORGQR_N ) |
$ N, N, -1 ) ) |
WRKBL = MAX( WRKBL, 3*N + LWORK_DGEBRD ) |
WRKBL = MAX( WRKBL, 3*N+2*N* |
WRKBL = MAX( WRKBL, 3*N + LWORK_DORGBR_Q ) |
$ ILAENV( 1, 'DGEBRD', ' ', N, N, -1, -1 ) ) |
WRKBL = MAX( WRKBL, 3*N + LWORK_DORGBR_P ) |
WRKBL = MAX( WRKBL, 3*N+N* |
|
$ ILAENV( 1, 'DORGBR', 'Q', N, N, N, -1 ) ) |
|
WRKBL = MAX( WRKBL, 3*N+( N-1 )* |
|
$ ILAENV( 1, 'DORGBR', 'P', N, N, N, -1 ) ) |
|
WRKBL = MAX( WRKBL, BDSPAC ) |
WRKBL = MAX( WRKBL, BDSPAC ) |
MAXWRK = 2*N*N + WRKBL |
MAXWRK = 2*N*N + WRKBL |
MINWRK = MAX( 3*N+M, BDSPAC ) |
MINWRK = MAX( 3*N + M, BDSPAC ) |
ELSE IF( WNTUS .AND. WNTVAS ) THEN |
ELSE IF( WNTUS .AND. WNTVAS ) THEN |
* |
* |
* Path 6 (M much larger than N, JOBU='S', JOBVT='S' or |
* Path 6 (M much larger than N, JOBU='S', JOBVT='S' or |
* 'A') |
* 'A') |
* |
* |
WRKBL = N + N*ILAENV( 1, 'DGEQRF', ' ', M, N, -1, -1 ) |
WRKBL = N + LWORK_DGEQRF |
WRKBL = MAX( WRKBL, N+N*ILAENV( 1, 'DORGQR', ' ', M, |
WRKBL = MAX( WRKBL, N + LWORK_DORGQR_N ) |
$ N, N, -1 ) ) |
WRKBL = MAX( WRKBL, 3*N + LWORK_DGEBRD ) |
WRKBL = MAX( WRKBL, 3*N+2*N* |
WRKBL = MAX( WRKBL, 3*N + LWORK_DORGBR_Q ) |
$ ILAENV( 1, 'DGEBRD', ' ', N, N, -1, -1 ) ) |
WRKBL = MAX( WRKBL, 3*N + LWORK_DORGBR_P ) |
WRKBL = MAX( WRKBL, 3*N+N* |
|
$ ILAENV( 1, 'DORGBR', 'Q', N, N, N, -1 ) ) |
|
WRKBL = MAX( WRKBL, 3*N+( N-1 )* |
|
$ ILAENV( 1, 'DORGBR', 'P', N, N, N, -1 ) ) |
|
WRKBL = MAX( WRKBL, BDSPAC ) |
WRKBL = MAX( WRKBL, BDSPAC ) |
MAXWRK = N*N + WRKBL |
MAXWRK = N*N + WRKBL |
MINWRK = MAX( 3*N+M, BDSPAC ) |
MINWRK = MAX( 3*N + M, BDSPAC ) |
ELSE IF( WNTUA .AND. WNTVN ) THEN |
ELSE IF( WNTUA .AND. WNTVN ) THEN |
* |
* |
* Path 7 (M much larger than N, JOBU='A', JOBVT='N') |
* Path 7 (M much larger than N, JOBU='A', JOBVT='N') |
* |
* |
WRKBL = N + N*ILAENV( 1, 'DGEQRF', ' ', M, N, -1, -1 ) |
WRKBL = N + LWORK_DGEQRF |
WRKBL = MAX( WRKBL, N+M*ILAENV( 1, 'DORGQR', ' ', M, |
WRKBL = MAX( WRKBL, N + LWORK_DORGQR_M ) |
$ M, N, -1 ) ) |
WRKBL = MAX( WRKBL, 3*N + LWORK_DGEBRD ) |
WRKBL = MAX( WRKBL, 3*N+2*N* |
WRKBL = MAX( WRKBL, 3*N + LWORK_DORGBR_Q ) |
$ ILAENV( 1, 'DGEBRD', ' ', N, N, -1, -1 ) ) |
|
WRKBL = MAX( WRKBL, 3*N+N* |
|
$ ILAENV( 1, 'DORGBR', 'Q', N, N, N, -1 ) ) |
|
WRKBL = MAX( WRKBL, BDSPAC ) |
WRKBL = MAX( WRKBL, BDSPAC ) |
MAXWRK = N*N + WRKBL |
MAXWRK = N*N + WRKBL |
MINWRK = MAX( 3*N+M, BDSPAC ) |
MINWRK = MAX( 3*N + M, BDSPAC ) |
ELSE IF( WNTUA .AND. WNTVO ) THEN |
ELSE IF( WNTUA .AND. WNTVO ) THEN |
* |
* |
* Path 8 (M much larger than N, JOBU='A', JOBVT='O') |
* Path 8 (M much larger than N, JOBU='A', JOBVT='O') |
* |
* |
WRKBL = N + N*ILAENV( 1, 'DGEQRF', ' ', M, N, -1, -1 ) |
WRKBL = N + LWORK_DGEQRF |
WRKBL = MAX( WRKBL, N+M*ILAENV( 1, 'DORGQR', ' ', M, |
WRKBL = MAX( WRKBL, N + LWORK_DORGQR_M ) |
$ M, N, -1 ) ) |
WRKBL = MAX( WRKBL, 3*N + LWORK_DGEBRD ) |
WRKBL = MAX( WRKBL, 3*N+2*N* |
WRKBL = MAX( WRKBL, 3*N + LWORK_DORGBR_Q ) |
$ ILAENV( 1, 'DGEBRD', ' ', N, N, -1, -1 ) ) |
WRKBL = MAX( WRKBL, 3*N + LWORK_DORGBR_P ) |
WRKBL = MAX( WRKBL, 3*N+N* |
|
$ ILAENV( 1, 'DORGBR', 'Q', N, N, N, -1 ) ) |
|
WRKBL = MAX( WRKBL, 3*N+( N-1 )* |
|
$ ILAENV( 1, 'DORGBR', 'P', N, N, N, -1 ) ) |
|
WRKBL = MAX( WRKBL, BDSPAC ) |
WRKBL = MAX( WRKBL, BDSPAC ) |
MAXWRK = 2*N*N + WRKBL |
MAXWRK = 2*N*N + WRKBL |
MINWRK = MAX( 3*N+M, BDSPAC ) |
MINWRK = MAX( 3*N + M, BDSPAC ) |
ELSE IF( WNTUA .AND. WNTVAS ) THEN |
ELSE IF( WNTUA .AND. WNTVAS ) THEN |
* |
* |
* Path 9 (M much larger than N, JOBU='A', JOBVT='S' or |
* Path 9 (M much larger than N, JOBU='A', JOBVT='S' or |
* 'A') |
* 'A') |
* |
* |
WRKBL = N + N*ILAENV( 1, 'DGEQRF', ' ', M, N, -1, -1 ) |
WRKBL = N + LWORK_DGEQRF |
WRKBL = MAX( WRKBL, N+M*ILAENV( 1, 'DORGQR', ' ', M, |
WRKBL = MAX( WRKBL, N + LWORK_DORGQR_M ) |
$ M, N, -1 ) ) |
WRKBL = MAX( WRKBL, 3*N + LWORK_DGEBRD ) |
WRKBL = MAX( WRKBL, 3*N+2*N* |
WRKBL = MAX( WRKBL, 3*N + LWORK_DORGBR_Q ) |
$ ILAENV( 1, 'DGEBRD', ' ', N, N, -1, -1 ) ) |
WRKBL = MAX( WRKBL, 3*N + LWORK_DORGBR_P ) |
WRKBL = MAX( WRKBL, 3*N+N* |
|
$ ILAENV( 1, 'DORGBR', 'Q', N, N, N, -1 ) ) |
|
WRKBL = MAX( WRKBL, 3*N+( N-1 )* |
|
$ ILAENV( 1, 'DORGBR', 'P', N, N, N, -1 ) ) |
|
WRKBL = MAX( WRKBL, BDSPAC ) |
WRKBL = MAX( WRKBL, BDSPAC ) |
MAXWRK = N*N + WRKBL |
MAXWRK = N*N + WRKBL |
MINWRK = MAX( 3*N+M, BDSPAC ) |
MINWRK = MAX( 3*N + M, BDSPAC ) |
END IF |
END IF |
ELSE |
ELSE |
* |
* |
* Path 10 (M at least N, but not much larger) |
* Path 10 (M at least N, but not much larger) |
* |
* |
MAXWRK = 3*N + ( M+N )*ILAENV( 1, 'DGEBRD', ' ', M, N, |
CALL DGEBRD( M, N, A, LDA, S, DUM(1), DUM(1), |
$ -1, -1 ) |
$ DUM(1), DUM(1), -1, IERR ) |
IF( WNTUS .OR. WNTUO ) |
LWORK_DGEBRD = INT( DUM(1) ) |
$ MAXWRK = MAX( MAXWRK, 3*N+N* |
MAXWRK = 3*N + LWORK_DGEBRD |
$ ILAENV( 1, 'DORGBR', 'Q', M, N, N, -1 ) ) |
IF( WNTUS .OR. WNTUO ) THEN |
IF( WNTUA ) |
CALL DORGBR( 'Q', M, N, N, A, LDA, DUM(1), |
$ MAXWRK = MAX( MAXWRK, 3*N+M* |
$ DUM(1), -1, IERR ) |
$ ILAENV( 1, 'DORGBR', 'Q', M, M, N, -1 ) ) |
LWORK_DORGBR_Q = INT( DUM(1) ) |
IF( .NOT.WNTVN ) |
MAXWRK = MAX( MAXWRK, 3*N + LWORK_DORGBR_Q ) |
$ MAXWRK = MAX( MAXWRK, 3*N+( N-1 )* |
END IF |
$ ILAENV( 1, 'DORGBR', 'P', N, N, N, -1 ) ) |
IF( WNTUA ) THEN |
|
CALL DORGBR( 'Q', M, M, N, A, LDA, DUM(1), |
|
$ DUM(1), -1, IERR ) |
|
LWORK_DORGBR_Q = INT( DUM(1) ) |
|
MAXWRK = MAX( MAXWRK, 3*N + LWORK_DORGBR_Q ) |
|
END IF |
|
IF( .NOT.WNTVN ) THEN |
|
MAXWRK = MAX( MAXWRK, 3*N + LWORK_DORGBR_P ) |
|
END IF |
MAXWRK = MAX( MAXWRK, BDSPAC ) |
MAXWRK = MAX( MAXWRK, BDSPAC ) |
MINWRK = MAX( 3*N+M, BDSPAC ) |
MINWRK = MAX( 3*N + M, BDSPAC ) |
END IF |
END IF |
ELSE IF( MINMN.GT.0 ) THEN |
ELSE IF( MINMN.GT.0 ) THEN |
* |
* |
Line 379
|
Line 473
|
* |
* |
MNTHR = ILAENV( 6, 'DGESVD', JOBU // JOBVT, M, N, 0, 0 ) |
MNTHR = ILAENV( 6, 'DGESVD', JOBU // JOBVT, M, N, 0, 0 ) |
BDSPAC = 5*M |
BDSPAC = 5*M |
|
* Compute space needed for DGELQF |
|
CALL DGELQF( M, N, A, LDA, DUM(1), DUM(1), -1, IERR ) |
|
LWORK_DGELQF = INT( DUM(1) ) |
|
* Compute space needed for DORGLQ |
|
CALL DORGLQ( N, N, M, DUM(1), N, DUM(1), DUM(1), -1, IERR ) |
|
LWORK_DORGLQ_N = INT( DUM(1) ) |
|
CALL DORGLQ( M, N, M, A, LDA, DUM(1), DUM(1), -1, IERR ) |
|
LWORK_DORGLQ_M = INT( DUM(1) ) |
|
* Compute space needed for DGEBRD |
|
CALL DGEBRD( M, M, A, LDA, S, DUM(1), DUM(1), |
|
$ DUM(1), DUM(1), -1, IERR ) |
|
LWORK_DGEBRD = INT( DUM(1) ) |
|
* Compute space needed for DORGBR P |
|
CALL DORGBR( 'P', M, M, M, A, N, DUM(1), |
|
$ DUM(1), -1, IERR ) |
|
LWORK_DORGBR_P = INT( DUM(1) ) |
|
* Compute space needed for DORGBR Q |
|
CALL DORGBR( 'Q', M, M, M, A, N, DUM(1), |
|
$ DUM(1), -1, IERR ) |
|
LWORK_DORGBR_Q = INT( DUM(1) ) |
IF( N.GE.MNTHR ) THEN |
IF( N.GE.MNTHR ) THEN |
IF( WNTVN ) THEN |
IF( WNTVN ) THEN |
* |
* |
* Path 1t(N much larger than M, JOBVT='N') |
* Path 1t(N much larger than M, JOBVT='N') |
* |
* |
MAXWRK = M + M*ILAENV( 1, 'DGELQF', ' ', M, N, -1, |
MAXWRK = M + LWORK_DGELQF |
$ -1 ) |
MAXWRK = MAX( MAXWRK, 3*M + LWORK_DGEBRD ) |
MAXWRK = MAX( MAXWRK, 3*M+2*M* |
|
$ ILAENV( 1, 'DGEBRD', ' ', M, M, -1, -1 ) ) |
|
IF( WNTUO .OR. WNTUAS ) |
IF( WNTUO .OR. WNTUAS ) |
$ MAXWRK = MAX( MAXWRK, 3*M+M* |
$ MAXWRK = MAX( MAXWRK, 3*M + LWORK_DORGBR_Q ) |
$ ILAENV( 1, 'DORGBR', 'Q', M, M, M, -1 ) ) |
|
MAXWRK = MAX( MAXWRK, BDSPAC ) |
MAXWRK = MAX( MAXWRK, BDSPAC ) |
MINWRK = MAX( 4*M, BDSPAC ) |
MINWRK = MAX( 4*M, BDSPAC ) |
ELSE IF( WNTVO .AND. WNTUN ) THEN |
ELSE IF( WNTVO .AND. WNTUN ) THEN |
* |
* |
* Path 2t(N much larger than M, JOBU='N', JOBVT='O') |
* Path 2t(N much larger than M, JOBU='N', JOBVT='O') |
* |
* |
WRKBL = M + M*ILAENV( 1, 'DGELQF', ' ', M, N, -1, -1 ) |
WRKBL = M + LWORK_DGELQF |
WRKBL = MAX( WRKBL, M+M*ILAENV( 1, 'DORGLQ', ' ', M, |
WRKBL = MAX( WRKBL, M + LWORK_DORGLQ_M ) |
$ N, M, -1 ) ) |
WRKBL = MAX( WRKBL, 3*M + LWORK_DGEBRD ) |
WRKBL = MAX( WRKBL, 3*M+2*M* |
WRKBL = MAX( WRKBL, 3*M + LWORK_DORGBR_P ) |
$ ILAENV( 1, 'DGEBRD', ' ', M, M, -1, -1 ) ) |
|
WRKBL = MAX( WRKBL, 3*M+( M-1 )* |
|
$ ILAENV( 1, 'DORGBR', 'P', M, M, M, -1 ) ) |
|
WRKBL = MAX( WRKBL, BDSPAC ) |
WRKBL = MAX( WRKBL, BDSPAC ) |
MAXWRK = MAX( M*M+WRKBL, M*M+M*N+M ) |
MAXWRK = MAX( M*M + WRKBL, M*M + M*N + M ) |
MINWRK = MAX( 3*M+N, BDSPAC ) |
MINWRK = MAX( 3*M + N, BDSPAC ) |
ELSE IF( WNTVO .AND. WNTUAS ) THEN |
ELSE IF( WNTVO .AND. WNTUAS ) THEN |
* |
* |
* Path 3t(N much larger than M, JOBU='S' or 'A', |
* Path 3t(N much larger than M, JOBU='S' or 'A', |
* JOBVT='O') |
* JOBVT='O') |
* |
* |
WRKBL = M + M*ILAENV( 1, 'DGELQF', ' ', M, N, -1, -1 ) |
WRKBL = M + LWORK_DGELQF |
WRKBL = MAX( WRKBL, M+M*ILAENV( 1, 'DORGLQ', ' ', M, |
WRKBL = MAX( WRKBL, M + LWORK_DORGLQ_M ) |
$ N, M, -1 ) ) |
WRKBL = MAX( WRKBL, 3*M + LWORK_DGEBRD ) |
WRKBL = MAX( WRKBL, 3*M+2*M* |
WRKBL = MAX( WRKBL, 3*M + LWORK_DORGBR_P ) |
$ ILAENV( 1, 'DGEBRD', ' ', M, M, -1, -1 ) ) |
WRKBL = MAX( WRKBL, 3*M + LWORK_DORGBR_Q ) |
WRKBL = MAX( WRKBL, 3*M+( M-1 )* |
|
$ ILAENV( 1, 'DORGBR', 'P', M, M, M, -1 ) ) |
|
WRKBL = MAX( WRKBL, 3*M+M* |
|
$ ILAENV( 1, 'DORGBR', 'Q', M, M, M, -1 ) ) |
|
WRKBL = MAX( WRKBL, BDSPAC ) |
WRKBL = MAX( WRKBL, BDSPAC ) |
MAXWRK = MAX( M*M+WRKBL, M*M+M*N+M ) |
MAXWRK = MAX( M*M + WRKBL, M*M + M*N + M ) |
MINWRK = MAX( 3*M+N, BDSPAC ) |
MINWRK = MAX( 3*M + N, BDSPAC ) |
ELSE IF( WNTVS .AND. WNTUN ) THEN |
ELSE IF( WNTVS .AND. WNTUN ) THEN |
* |
* |
* Path 4t(N much larger than M, JOBU='N', JOBVT='S') |
* Path 4t(N much larger than M, JOBU='N', JOBVT='S') |
* |
* |
WRKBL = M + M*ILAENV( 1, 'DGELQF', ' ', M, N, -1, -1 ) |
WRKBL = M + LWORK_DGELQF |
WRKBL = MAX( WRKBL, M+M*ILAENV( 1, 'DORGLQ', ' ', M, |
WRKBL = MAX( WRKBL, M + LWORK_DORGLQ_M ) |
$ N, M, -1 ) ) |
WRKBL = MAX( WRKBL, 3*M + LWORK_DGEBRD ) |
WRKBL = MAX( WRKBL, 3*M+2*M* |
WRKBL = MAX( WRKBL, 3*M + LWORK_DORGBR_P ) |
$ ILAENV( 1, 'DGEBRD', ' ', M, M, -1, -1 ) ) |
|
WRKBL = MAX( WRKBL, 3*M+( M-1 )* |
|
$ ILAENV( 1, 'DORGBR', 'P', M, M, M, -1 ) ) |
|
WRKBL = MAX( WRKBL, BDSPAC ) |
WRKBL = MAX( WRKBL, BDSPAC ) |
MAXWRK = M*M + WRKBL |
MAXWRK = M*M + WRKBL |
MINWRK = MAX( 3*M+N, BDSPAC ) |
MINWRK = MAX( 3*M + N, BDSPAC ) |
ELSE IF( WNTVS .AND. WNTUO ) THEN |
ELSE IF( WNTVS .AND. WNTUO ) THEN |
* |
* |
* Path 5t(N much larger than M, JOBU='O', JOBVT='S') |
* Path 5t(N much larger than M, JOBU='O', JOBVT='S') |
* |
* |
WRKBL = M + M*ILAENV( 1, 'DGELQF', ' ', M, N, -1, -1 ) |
WRKBL = M + LWORK_DGELQF |
WRKBL = MAX( WRKBL, M+M*ILAENV( 1, 'DORGLQ', ' ', M, |
WRKBL = MAX( WRKBL, M + LWORK_DORGLQ_M ) |
$ N, M, -1 ) ) |
WRKBL = MAX( WRKBL, 3*M + LWORK_DGEBRD ) |
WRKBL = MAX( WRKBL, 3*M+2*M* |
WRKBL = MAX( WRKBL, 3*M + LWORK_DORGBR_P ) |
$ ILAENV( 1, 'DGEBRD', ' ', M, M, -1, -1 ) ) |
WRKBL = MAX( WRKBL, 3*M + LWORK_DORGBR_Q ) |
WRKBL = MAX( WRKBL, 3*M+( M-1 )* |
|
$ ILAENV( 1, 'DORGBR', 'P', M, M, M, -1 ) ) |
|
WRKBL = MAX( WRKBL, 3*M+M* |
|
$ ILAENV( 1, 'DORGBR', 'Q', M, M, M, -1 ) ) |
|
WRKBL = MAX( WRKBL, BDSPAC ) |
WRKBL = MAX( WRKBL, BDSPAC ) |
MAXWRK = 2*M*M + WRKBL |
MAXWRK = 2*M*M + WRKBL |
MINWRK = MAX( 3*M+N, BDSPAC ) |
MINWRK = MAX( 3*M + N, BDSPAC ) |
ELSE IF( WNTVS .AND. WNTUAS ) THEN |
ELSE IF( WNTVS .AND. WNTUAS ) THEN |
* |
* |
* Path 6t(N much larger than M, JOBU='S' or 'A', |
* Path 6t(N much larger than M, JOBU='S' or 'A', |
* JOBVT='S') |
* JOBVT='S') |
* |
* |
WRKBL = M + M*ILAENV( 1, 'DGELQF', ' ', M, N, -1, -1 ) |
WRKBL = M + LWORK_DGELQF |
WRKBL = MAX( WRKBL, M+M*ILAENV( 1, 'DORGLQ', ' ', M, |
WRKBL = MAX( WRKBL, M + LWORK_DORGLQ_M ) |
$ N, M, -1 ) ) |
WRKBL = MAX( WRKBL, 3*M + LWORK_DGEBRD ) |
WRKBL = MAX( WRKBL, 3*M+2*M* |
WRKBL = MAX( WRKBL, 3*M + LWORK_DORGBR_P ) |
$ ILAENV( 1, 'DGEBRD', ' ', M, M, -1, -1 ) ) |
WRKBL = MAX( WRKBL, 3*M + LWORK_DORGBR_Q ) |
WRKBL = MAX( WRKBL, 3*M+( M-1 )* |
|
$ ILAENV( 1, 'DORGBR', 'P', M, M, M, -1 ) ) |
|
WRKBL = MAX( WRKBL, 3*M+M* |
|
$ ILAENV( 1, 'DORGBR', 'Q', M, M, M, -1 ) ) |
|
WRKBL = MAX( WRKBL, BDSPAC ) |
WRKBL = MAX( WRKBL, BDSPAC ) |
MAXWRK = M*M + WRKBL |
MAXWRK = M*M + WRKBL |
MINWRK = MAX( 3*M+N, BDSPAC ) |
MINWRK = MAX( 3*M + N, BDSPAC ) |
ELSE IF( WNTVA .AND. WNTUN ) THEN |
ELSE IF( WNTVA .AND. WNTUN ) THEN |
* |
* |
* Path 7t(N much larger than M, JOBU='N', JOBVT='A') |
* Path 7t(N much larger than M, JOBU='N', JOBVT='A') |
* |
* |
WRKBL = M + M*ILAENV( 1, 'DGELQF', ' ', M, N, -1, -1 ) |
WRKBL = M + LWORK_DGELQF |
WRKBL = MAX( WRKBL, M+N*ILAENV( 1, 'DORGLQ', ' ', N, |
WRKBL = MAX( WRKBL, M + LWORK_DORGLQ_N ) |
$ N, M, -1 ) ) |
WRKBL = MAX( WRKBL, 3*M + LWORK_DGEBRD ) |
WRKBL = MAX( WRKBL, 3*M+2*M* |
WRKBL = MAX( WRKBL, 3*M + LWORK_DORGBR_P ) |
$ ILAENV( 1, 'DGEBRD', ' ', M, M, -1, -1 ) ) |
|
WRKBL = MAX( WRKBL, 3*M+( M-1 )* |
|
$ ILAENV( 1, 'DORGBR', 'P', M, M, M, -1 ) ) |
|
WRKBL = MAX( WRKBL, BDSPAC ) |
WRKBL = MAX( WRKBL, BDSPAC ) |
MAXWRK = M*M + WRKBL |
MAXWRK = M*M + WRKBL |
MINWRK = MAX( 3*M+N, BDSPAC ) |
MINWRK = MAX( 3*M + N, BDSPAC ) |
ELSE IF( WNTVA .AND. WNTUO ) THEN |
ELSE IF( WNTVA .AND. WNTUO ) THEN |
* |
* |
* Path 8t(N much larger than M, JOBU='O', JOBVT='A') |
* Path 8t(N much larger than M, JOBU='O', JOBVT='A') |
* |
* |
WRKBL = M + M*ILAENV( 1, 'DGELQF', ' ', M, N, -1, -1 ) |
WRKBL = M + LWORK_DGELQF |
WRKBL = MAX( WRKBL, M+N*ILAENV( 1, 'DORGLQ', ' ', N, |
WRKBL = MAX( WRKBL, M + LWORK_DORGLQ_N ) |
$ N, M, -1 ) ) |
WRKBL = MAX( WRKBL, 3*M + LWORK_DGEBRD ) |
WRKBL = MAX( WRKBL, 3*M+2*M* |
WRKBL = MAX( WRKBL, 3*M + LWORK_DORGBR_P ) |
$ ILAENV( 1, 'DGEBRD', ' ', M, M, -1, -1 ) ) |
WRKBL = MAX( WRKBL, 3*M + LWORK_DORGBR_Q ) |
WRKBL = MAX( WRKBL, 3*M+( M-1 )* |
|
$ ILAENV( 1, 'DORGBR', 'P', M, M, M, -1 ) ) |
|
WRKBL = MAX( WRKBL, 3*M+M* |
|
$ ILAENV( 1, 'DORGBR', 'Q', M, M, M, -1 ) ) |
|
WRKBL = MAX( WRKBL, BDSPAC ) |
WRKBL = MAX( WRKBL, BDSPAC ) |
MAXWRK = 2*M*M + WRKBL |
MAXWRK = 2*M*M + WRKBL |
MINWRK = MAX( 3*M+N, BDSPAC ) |
MINWRK = MAX( 3*M + N, BDSPAC ) |
ELSE IF( WNTVA .AND. WNTUAS ) THEN |
ELSE IF( WNTVA .AND. WNTUAS ) THEN |
* |
* |
* Path 9t(N much larger than M, JOBU='S' or 'A', |
* Path 9t(N much larger than M, JOBU='S' or 'A', |
* JOBVT='A') |
* JOBVT='A') |
* |
* |
WRKBL = M + M*ILAENV( 1, 'DGELQF', ' ', M, N, -1, -1 ) |
WRKBL = M + LWORK_DGELQF |
WRKBL = MAX( WRKBL, M+N*ILAENV( 1, 'DORGLQ', ' ', N, |
WRKBL = MAX( WRKBL, M + LWORK_DORGLQ_N ) |
$ N, M, -1 ) ) |
WRKBL = MAX( WRKBL, 3*M + LWORK_DGEBRD ) |
WRKBL = MAX( WRKBL, 3*M+2*M* |
WRKBL = MAX( WRKBL, 3*M + LWORK_DORGBR_P ) |
$ ILAENV( 1, 'DGEBRD', ' ', M, M, -1, -1 ) ) |
WRKBL = MAX( WRKBL, 3*M + LWORK_DORGBR_Q ) |
WRKBL = MAX( WRKBL, 3*M+( M-1 )* |
|
$ ILAENV( 1, 'DORGBR', 'P', M, M, M, -1 ) ) |
|
WRKBL = MAX( WRKBL, 3*M+M* |
|
$ ILAENV( 1, 'DORGBR', 'Q', M, M, M, -1 ) ) |
|
WRKBL = MAX( WRKBL, BDSPAC ) |
WRKBL = MAX( WRKBL, BDSPAC ) |
MAXWRK = M*M + WRKBL |
MAXWRK = M*M + WRKBL |
MINWRK = MAX( 3*M+N, BDSPAC ) |
MINWRK = MAX( 3*M + N, BDSPAC ) |
END IF |
END IF |
ELSE |
ELSE |
* |
* |
* Path 10t(N greater than M, but not much larger) |
* Path 10t(N greater than M, but not much larger) |
* |
* |
MAXWRK = 3*M + ( M+N )*ILAENV( 1, 'DGEBRD', ' ', M, N, |
CALL DGEBRD( M, N, A, LDA, S, DUM(1), DUM(1), |
$ -1, -1 ) |
$ DUM(1), DUM(1), -1, IERR ) |
IF( WNTVS .OR. WNTVO ) |
LWORK_DGEBRD = INT( DUM(1) ) |
$ MAXWRK = MAX( MAXWRK, 3*M+M* |
MAXWRK = 3*M + LWORK_DGEBRD |
$ ILAENV( 1, 'DORGBR', 'P', M, N, M, -1 ) ) |
IF( WNTVS .OR. WNTVO ) THEN |
IF( WNTVA ) |
* Compute space needed for DORGBR P |
$ MAXWRK = MAX( MAXWRK, 3*M+N* |
CALL DORGBR( 'P', M, N, M, A, N, DUM(1), |
$ ILAENV( 1, 'DORGBR', 'P', N, N, M, -1 ) ) |
$ DUM(1), -1, IERR ) |
IF( .NOT.WNTUN ) |
LWORK_DORGBR_P = INT( DUM(1) ) |
$ MAXWRK = MAX( MAXWRK, 3*M+( M-1 )* |
MAXWRK = MAX( MAXWRK, 3*M + LWORK_DORGBR_P ) |
$ ILAENV( 1, 'DORGBR', 'Q', M, M, M, -1 ) ) |
END IF |
|
IF( WNTVA ) THEN |
|
CALL DORGBR( 'P', N, N, M, A, N, DUM(1), |
|
$ DUM(1), -1, IERR ) |
|
LWORK_DORGBR_P = INT( DUM(1) ) |
|
MAXWRK = MAX( MAXWRK, 3*M + LWORK_DORGBR_P ) |
|
END IF |
|
IF( .NOT.WNTUN ) THEN |
|
MAXWRK = MAX( MAXWRK, 3*M + LWORK_DORGBR_Q ) |
|
END IF |
MAXWRK = MAX( MAXWRK, BDSPAC ) |
MAXWRK = MAX( MAXWRK, BDSPAC ) |
MINWRK = MAX( 3*M+N, BDSPAC ) |
MINWRK = MAX( 3*M + N, BDSPAC ) |
END IF |
END IF |
END IF |
END IF |
MAXWRK = MAX( MAXWRK, MINWRK ) |
MAXWRK = MAX( MAXWRK, MINWRK ) |
Line 594
|
Line 685
|
IWORK = ITAU + N |
IWORK = ITAU + N |
* |
* |
* Compute A=Q*R |
* Compute A=Q*R |
* (Workspace: need 2*N, prefer N+N*NB) |
* (Workspace: need 2*N, prefer N + N*NB) |
* |
* |
CALL DGEQRF( M, N, A, LDA, WORK( ITAU ), WORK( IWORK ), |
CALL DGEQRF( M, N, A, LDA, WORK( ITAU ), WORK( IWORK ), |
$ LWORK-IWORK+1, IERR ) |
$ LWORK-IWORK+1, IERR ) |
* |
* |
* Zero out below R |
* Zero out below R |
* |
* |
CALL DLASET( 'L', N-1, N-1, ZERO, ZERO, A( 2, 1 ), LDA ) |
IF( N .GT. 1 ) THEN |
|
CALL DLASET( 'L', N-1, N-1, ZERO, ZERO, A( 2, 1 ), |
|
$ LDA ) |
|
END IF |
IE = 1 |
IE = 1 |
ITAUQ = IE + N |
ITAUQ = IE + N |
ITAUP = ITAUQ + N |
ITAUP = ITAUQ + N |
IWORK = ITAUP + N |
IWORK = ITAUP + N |
* |
* |
* Bidiagonalize R in A |
* Bidiagonalize R in A |
* (Workspace: need 4*N, prefer 3*N+2*N*NB) |
* (Workspace: need 4*N, prefer 3*N + 2*N*NB) |
* |
* |
CALL DGEBRD( N, N, A, LDA, S, WORK( IE ), WORK( ITAUQ ), |
CALL DGEBRD( N, N, A, LDA, S, WORK( IE ), WORK( ITAUQ ), |
$ WORK( ITAUP ), WORK( IWORK ), LWORK-IWORK+1, |
$ WORK( ITAUP ), WORK( IWORK ), LWORK-IWORK+1, |
Line 617
|
Line 711
|
IF( WNTVO .OR. WNTVAS ) THEN |
IF( WNTVO .OR. WNTVAS ) THEN |
* |
* |
* If right singular vectors desired, generate P'. |
* If right singular vectors desired, generate P'. |
* (Workspace: need 4*N-1, prefer 3*N+(N-1)*NB) |
* (Workspace: need 4*N-1, prefer 3*N + (N-1)*NB) |
* |
* |
CALL DORGBR( 'P', N, N, N, A, LDA, WORK( ITAUP ), |
CALL DORGBR( 'P', N, N, N, A, LDA, WORK( ITAUP ), |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
Line 648
|
Line 742
|
* Sufficient workspace for a fast algorithm |
* Sufficient workspace for a fast algorithm |
* |
* |
IR = 1 |
IR = 1 |
IF( LWORK.GE.MAX( WRKBL, LDA*N+N )+LDA*N ) THEN |
IF( LWORK.GE.MAX( WRKBL, LDA*N + N ) + LDA*N ) THEN |
* |
* |
* WORK(IU) is LDA by N, WORK(IR) is LDA by N |
* WORK(IU) is LDA by N, WORK(IR) is LDA by N |
* |
* |
LDWRKU = LDA |
LDWRKU = LDA |
LDWRKR = LDA |
LDWRKR = LDA |
ELSE IF( LWORK.GE.MAX( WRKBL, LDA*N+N )+N*N ) THEN |
ELSE IF( LWORK.GE.MAX( WRKBL, LDA*N + N ) + N*N ) THEN |
* |
* |
* WORK(IU) is LDA by N, WORK(IR) is N by N |
* WORK(IU) is LDA by N, WORK(IR) is N by N |
* |
* |
Line 671
|
Line 765
|
IWORK = ITAU + N |
IWORK = ITAU + N |
* |
* |
* Compute A=Q*R |
* Compute A=Q*R |
* (Workspace: need N*N+2*N, prefer N*N+N+N*NB) |
* (Workspace: need N*N + 2*N, prefer N*N + N + N*NB) |
* |
* |
CALL DGEQRF( M, N, A, LDA, WORK( ITAU ), |
CALL DGEQRF( M, N, A, LDA, WORK( ITAU ), |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
Line 683
|
Line 777
|
$ LDWRKR ) |
$ LDWRKR ) |
* |
* |
* Generate Q in A |
* Generate Q in A |
* (Workspace: need N*N+2*N, prefer N*N+N+N*NB) |
* (Workspace: need N*N + 2*N, prefer N*N + N + N*NB) |
* |
* |
CALL DORGQR( M, N, N, A, LDA, WORK( ITAU ), |
CALL DORGQR( M, N, N, A, LDA, WORK( ITAU ), |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
Line 693
|
Line 787
|
IWORK = ITAUP + N |
IWORK = ITAUP + N |
* |
* |
* Bidiagonalize R in WORK(IR) |
* Bidiagonalize R in WORK(IR) |
* (Workspace: need N*N+4*N, prefer N*N+3*N+2*N*NB) |
* (Workspace: need N*N + 4*N, prefer N*N + 3*N + 2*N*NB) |
* |
* |
CALL DGEBRD( N, N, WORK( IR ), LDWRKR, S, WORK( IE ), |
CALL DGEBRD( N, N, WORK( IR ), LDWRKR, S, WORK( IE ), |
$ WORK( ITAUQ ), WORK( ITAUP ), |
$ WORK( ITAUQ ), WORK( ITAUP ), |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
* |
* |
* Generate left vectors bidiagonalizing R |
* Generate left vectors bidiagonalizing R |
* (Workspace: need N*N+4*N, prefer N*N+3*N+N*NB) |
* (Workspace: need N*N + 4*N, prefer N*N + 3*N + N*NB) |
* |
* |
CALL DORGBR( 'Q', N, N, N, WORK( IR ), LDWRKR, |
CALL DORGBR( 'Q', N, N, N, WORK( IR ), LDWRKR, |
$ WORK( ITAUQ ), WORK( IWORK ), |
$ WORK( ITAUQ ), WORK( IWORK ), |
Line 709
|
Line 803
|
* |
* |
* Perform bidiagonal QR iteration, computing left |
* Perform bidiagonal QR iteration, computing left |
* singular vectors of R in WORK(IR) |
* singular vectors of R in WORK(IR) |
* (Workspace: need N*N+BDSPAC) |
* (Workspace: need N*N + BDSPAC) |
* |
* |
CALL DBDSQR( 'U', N, 0, N, 0, S, WORK( IE ), DUM, 1, |
CALL DBDSQR( 'U', N, 0, N, 0, S, WORK( IE ), DUM, 1, |
$ WORK( IR ), LDWRKR, DUM, 1, |
$ WORK( IR ), LDWRKR, DUM, 1, |
Line 718
|
Line 812
|
* |
* |
* 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 WORK(IU) and copying to A |
* WORK(IR), storing result in WORK(IU) and copying to A |
* (Workspace: need N*N+2*N, prefer N*N+M*N+N) |
* (Workspace: need N*N + 2*N, prefer N*N + M*N + N) |
* |
* |
DO 10 I = 1, M, LDWRKU |
DO 10 I = 1, M, LDWRKU |
CHUNK = MIN( M-I+1, LDWRKU ) |
CHUNK = MIN( M-I+1, LDWRKU ) |
Line 739
|
Line 833
|
IWORK = ITAUP + N |
IWORK = ITAUP + N |
* |
* |
* Bidiagonalize A |
* Bidiagonalize A |
* (Workspace: need 3*N+M, prefer 3*N+(M+N)*NB) |
* (Workspace: need 3*N + M, prefer 3*N + (M + N)*NB) |
* |
* |
CALL DGEBRD( M, N, A, LDA, S, WORK( IE ), |
CALL DGEBRD( M, N, A, LDA, S, WORK( IE ), |
$ WORK( ITAUQ ), WORK( ITAUP ), |
$ WORK( ITAUQ ), WORK( ITAUP ), |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
* |
* |
* Generate left vectors bidiagonalizing A |
* Generate left vectors bidiagonalizing A |
* (Workspace: need 4*N, prefer 3*N+N*NB) |
* (Workspace: need 4*N, prefer 3*N + N*NB) |
* |
* |
CALL DORGBR( 'Q', M, N, N, A, LDA, WORK( ITAUQ ), |
CALL DORGBR( 'Q', M, N, N, A, LDA, WORK( ITAUQ ), |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
Line 772
|
Line 866
|
* Sufficient workspace for a fast algorithm |
* Sufficient workspace for a fast algorithm |
* |
* |
IR = 1 |
IR = 1 |
IF( LWORK.GE.MAX( WRKBL, LDA*N+N )+LDA*N ) THEN |
IF( LWORK.GE.MAX( WRKBL, LDA*N + N ) + LDA*N ) THEN |
* |
* |
* WORK(IU) is LDA by N and WORK(IR) is LDA by N |
* WORK(IU) is LDA by N and WORK(IR) is LDA by N |
* |
* |
LDWRKU = LDA |
LDWRKU = LDA |
LDWRKR = LDA |
LDWRKR = LDA |
ELSE IF( LWORK.GE.MAX( WRKBL, LDA*N+N )+N*N ) THEN |
ELSE IF( LWORK.GE.MAX( WRKBL, LDA*N + N ) + N*N ) THEN |
* |
* |
* WORK(IU) is LDA by N and WORK(IR) is N by N |
* WORK(IU) is LDA by N and WORK(IR) is N by N |
* |
* |
Line 795
|
Line 889
|
IWORK = ITAU + N |
IWORK = ITAU + N |
* |
* |
* Compute A=Q*R |
* Compute A=Q*R |
* (Workspace: need N*N+2*N, prefer N*N+N+N*NB) |
* (Workspace: need N*N + 2*N, prefer N*N + N + N*NB) |
* |
* |
CALL DGEQRF( M, N, A, LDA, WORK( ITAU ), |
CALL DGEQRF( M, N, A, LDA, WORK( ITAU ), |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
Line 808
|
Line 902
|
$ VT( 2, 1 ), LDVT ) |
$ VT( 2, 1 ), LDVT ) |
* |
* |
* Generate Q in A |
* Generate Q in A |
* (Workspace: need N*N+2*N, prefer N*N+N+N*NB) |
* (Workspace: need N*N + 2*N, prefer N*N + N + N*NB) |
* |
* |
CALL DORGQR( M, N, N, A, LDA, WORK( ITAU ), |
CALL DORGQR( M, N, N, A, LDA, WORK( ITAU ), |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
Line 818
|
Line 912
|
IWORK = ITAUP + N |
IWORK = ITAUP + N |
* |
* |
* Bidiagonalize R in VT, copying result to WORK(IR) |
* Bidiagonalize R in VT, copying result to WORK(IR) |
* (Workspace: need N*N+4*N, prefer N*N+3*N+2*N*NB) |
* (Workspace: need N*N + 4*N, prefer N*N + 3*N + 2*N*NB) |
* |
* |
CALL DGEBRD( N, N, VT, LDVT, S, WORK( IE ), |
CALL DGEBRD( N, N, VT, LDVT, S, WORK( IE ), |
$ WORK( ITAUQ ), WORK( ITAUP ), |
$ WORK( ITAUQ ), WORK( ITAUP ), |
Line 826
|
Line 920
|
CALL DLACPY( 'L', N, N, VT, LDVT, WORK( IR ), LDWRKR ) |
CALL DLACPY( 'L', N, N, VT, LDVT, WORK( IR ), LDWRKR ) |
* |
* |
* Generate left vectors bidiagonalizing R in WORK(IR) |
* Generate left vectors bidiagonalizing R in WORK(IR) |
* (Workspace: need N*N+4*N, prefer N*N+3*N+N*NB) |
* (Workspace: need N*N + 4*N, prefer N*N + 3*N + N*NB) |
* |
* |
CALL DORGBR( 'Q', N, N, N, WORK( IR ), LDWRKR, |
CALL DORGBR( 'Q', N, N, N, WORK( IR ), LDWRKR, |
$ WORK( ITAUQ ), WORK( IWORK ), |
$ WORK( ITAUQ ), WORK( IWORK ), |
$ LWORK-IWORK+1, IERR ) |
$ LWORK-IWORK+1, IERR ) |
* |
* |
* Generate right vectors bidiagonalizing R in VT |
* Generate right vectors bidiagonalizing R in VT |
* (Workspace: need N*N+4*N-1, prefer N*N+3*N+(N-1)*NB) |
* (Workspace: need N*N + 4*N-1, prefer N*N + 3*N + (N-1)*NB) |
* |
* |
CALL DORGBR( 'P', N, N, N, VT, LDVT, WORK( ITAUP ), |
CALL DORGBR( 'P', N, N, N, VT, LDVT, WORK( ITAUP ), |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
Line 842
|
Line 936
|
* Perform bidiagonal QR iteration, computing left |
* Perform bidiagonal QR iteration, computing left |
* singular vectors of R in WORK(IR) and computing right |
* singular vectors of R in WORK(IR) and computing right |
* singular vectors of R in VT |
* singular vectors of R in VT |
* (Workspace: need N*N+BDSPAC) |
* (Workspace: need N*N + BDSPAC) |
* |
* |
CALL DBDSQR( 'U', N, N, N, 0, S, WORK( IE ), VT, LDVT, |
CALL DBDSQR( 'U', N, N, N, 0, S, WORK( IE ), VT, LDVT, |
$ WORK( IR ), LDWRKR, DUM, 1, |
$ WORK( IR ), LDWRKR, DUM, 1, |
Line 851
|
Line 945
|
* |
* |
* 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 WORK(IU) and copying to A |
* WORK(IR), storing result in WORK(IU) and copying to A |
* (Workspace: need N*N+2*N, prefer N*N+M*N+N) |
* (Workspace: need N*N + 2*N, prefer N*N + M*N + N) |
* |
* |
DO 20 I = 1, M, LDWRKU |
DO 20 I = 1, M, LDWRKU |
CHUNK = MIN( M-I+1, LDWRKU ) |
CHUNK = MIN( M-I+1, LDWRKU ) |
Line 870
|
Line 964
|
IWORK = ITAU + N |
IWORK = ITAU + N |
* |
* |
* Compute A=Q*R |
* Compute A=Q*R |
* (Workspace: need 2*N, prefer N+N*NB) |
* (Workspace: need 2*N, prefer N + N*NB) |
* |
* |
CALL DGEQRF( M, N, A, LDA, WORK( ITAU ), |
CALL DGEQRF( M, N, A, LDA, WORK( ITAU ), |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
Line 883
|
Line 977
|
$ VT( 2, 1 ), LDVT ) |
$ VT( 2, 1 ), LDVT ) |
* |
* |
* Generate Q in A |
* Generate Q in A |
* (Workspace: need 2*N, prefer N+N*NB) |
* (Workspace: need 2*N, prefer N + N*NB) |
* |
* |
CALL DORGQR( M, N, N, A, LDA, WORK( ITAU ), |
CALL DORGQR( M, N, N, A, LDA, WORK( ITAU ), |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
Line 893
|
Line 987
|
IWORK = ITAUP + N |
IWORK = ITAUP + N |
* |
* |
* Bidiagonalize R in VT |
* Bidiagonalize R in VT |
* (Workspace: need 4*N, prefer 3*N+2*N*NB) |
* (Workspace: need 4*N, prefer 3*N + 2*N*NB) |
* |
* |
CALL DGEBRD( N, N, VT, LDVT, S, WORK( IE ), |
CALL DGEBRD( N, N, VT, LDVT, S, WORK( IE ), |
$ WORK( ITAUQ ), WORK( ITAUP ), |
$ WORK( ITAUQ ), WORK( ITAUP ), |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
* |
* |
* Multiply Q in A by left vectors bidiagonalizing R |
* Multiply Q in A by left vectors bidiagonalizing R |
* (Workspace: need 3*N+M, prefer 3*N+M*NB) |
* (Workspace: need 3*N + M, prefer 3*N + M*NB) |
* |
* |
CALL DORMBR( 'Q', 'R', 'N', M, N, N, VT, LDVT, |
CALL DORMBR( 'Q', 'R', 'N', M, N, N, VT, LDVT, |
$ WORK( ITAUQ ), A, LDA, WORK( IWORK ), |
$ WORK( ITAUQ ), A, LDA, WORK( IWORK ), |
$ LWORK-IWORK+1, IERR ) |
$ LWORK-IWORK+1, IERR ) |
* |
* |
* Generate right vectors bidiagonalizing R in VT |
* Generate right vectors bidiagonalizing R in VT |
* (Workspace: need 4*N-1, prefer 3*N+(N-1)*NB) |
* (Workspace: need 4*N-1, prefer 3*N + (N-1)*NB) |
* |
* |
CALL DORGBR( 'P', N, N, N, VT, LDVT, WORK( ITAUP ), |
CALL DORGBR( 'P', N, N, N, VT, LDVT, WORK( ITAUP ), |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
Line 951
|
Line 1045
|
IWORK = ITAU + N |
IWORK = ITAU + N |
* |
* |
* Compute A=Q*R |
* Compute A=Q*R |
* (Workspace: need N*N+2*N, prefer N*N+N+N*NB) |
* (Workspace: need N*N + 2*N, prefer N*N + N + N*NB) |
* |
* |
CALL DGEQRF( M, N, A, LDA, WORK( ITAU ), |
CALL DGEQRF( M, N, A, LDA, WORK( ITAU ), |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
Line 964
|
Line 1058
|
$ WORK( IR+1 ), LDWRKR ) |
$ WORK( IR+1 ), LDWRKR ) |
* |
* |
* Generate Q in A |
* Generate Q in A |
* (Workspace: need N*N+2*N, prefer N*N+N+N*NB) |
* (Workspace: need N*N + 2*N, prefer N*N + N + N*NB) |
* |
* |
CALL DORGQR( M, N, N, A, LDA, WORK( ITAU ), |
CALL DORGQR( M, N, N, A, LDA, WORK( ITAU ), |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
Line 974
|
Line 1068
|
IWORK = ITAUP + N |
IWORK = ITAUP + N |
* |
* |
* Bidiagonalize R in WORK(IR) |
* Bidiagonalize R in WORK(IR) |
* (Workspace: need N*N+4*N, prefer N*N+3*N+2*N*NB) |
* (Workspace: need N*N + 4*N, prefer N*N + 3*N + 2*N*NB) |
* |
* |
CALL DGEBRD( N, N, WORK( IR ), LDWRKR, S, |
CALL DGEBRD( N, N, WORK( IR ), LDWRKR, S, |
$ WORK( IE ), WORK( ITAUQ ), |
$ WORK( IE ), WORK( ITAUQ ), |
Line 982
|
Line 1076
|
$ LWORK-IWORK+1, IERR ) |
$ LWORK-IWORK+1, IERR ) |
* |
* |
* Generate left vectors bidiagonalizing R in WORK(IR) |
* Generate left vectors bidiagonalizing R in WORK(IR) |
* (Workspace: need N*N+4*N, prefer N*N+3*N+N*NB) |
* (Workspace: need N*N + 4*N, prefer N*N + 3*N + N*NB) |
* |
* |
CALL DORGBR( 'Q', N, N, N, WORK( IR ), LDWRKR, |
CALL DORGBR( 'Q', N, N, N, WORK( IR ), LDWRKR, |
$ WORK( ITAUQ ), WORK( IWORK ), |
$ WORK( ITAUQ ), WORK( IWORK ), |
Line 991
|
Line 1085
|
* |
* |
* Perform bidiagonal QR iteration, computing left |
* Perform bidiagonal QR iteration, computing left |
* singular vectors of R in WORK(IR) |
* singular vectors of R in WORK(IR) |
* (Workspace: need N*N+BDSPAC) |
* (Workspace: need N*N + BDSPAC) |
* |
* |
CALL DBDSQR( 'U', N, 0, N, 0, S, WORK( IE ), DUM, |
CALL DBDSQR( 'U', N, 0, N, 0, S, WORK( IE ), DUM, |
$ 1, WORK( IR ), LDWRKR, DUM, 1, |
$ 1, WORK( IR ), LDWRKR, DUM, 1, |
Line 1012
|
Line 1106
|
IWORK = ITAU + N |
IWORK = ITAU + N |
* |
* |
* Compute A=Q*R, copying result to U |
* Compute A=Q*R, copying result to U |
* (Workspace: need 2*N, prefer N+N*NB) |
* (Workspace: need 2*N, prefer N + N*NB) |
* |
* |
CALL DGEQRF( M, N, A, LDA, WORK( ITAU ), |
CALL DGEQRF( M, N, A, LDA, WORK( ITAU ), |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
CALL DLACPY( 'L', M, N, A, LDA, U, LDU ) |
CALL DLACPY( 'L', M, N, A, LDA, U, LDU ) |
* |
* |
* Generate Q in U |
* Generate Q in U |
* (Workspace: need 2*N, prefer N+N*NB) |
* (Workspace: need 2*N, prefer N + N*NB) |
* |
* |
CALL DORGQR( M, N, N, U, LDU, WORK( ITAU ), |
CALL DORGQR( M, N, N, U, LDU, WORK( ITAU ), |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
Line 1030
|
Line 1124
|
* |
* |
* Zero out below R in A |
* Zero out below R in A |
* |
* |
CALL DLASET( 'L', N-1, N-1, ZERO, ZERO, A( 2, 1 ), |
IF( N .GT. 1 ) THEN |
$ LDA ) |
CALL DLASET( 'L', N-1, N-1, ZERO, ZERO, |
|
$ A( 2, 1 ), LDA ) |
|
END IF |
* |
* |
* Bidiagonalize R in A |
* Bidiagonalize R in A |
* (Workspace: need 4*N, prefer 3*N+2*N*NB) |
* (Workspace: need 4*N, prefer 3*N + 2*N*NB) |
* |
* |
CALL DGEBRD( N, N, A, LDA, S, WORK( IE ), |
CALL DGEBRD( N, N, A, LDA, S, WORK( IE ), |
$ WORK( ITAUQ ), WORK( ITAUP ), |
$ WORK( ITAUQ ), WORK( ITAUP ), |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
* |
* |
* Multiply Q in U by left vectors bidiagonalizing R |
* Multiply Q in U by left vectors bidiagonalizing R |
* (Workspace: need 3*N+M, prefer 3*N+M*NB) |
* (Workspace: need 3*N + M, prefer 3*N + M*NB) |
* |
* |
CALL DORMBR( 'Q', 'R', 'N', M, N, N, A, LDA, |
CALL DORMBR( 'Q', 'R', 'N', M, N, N, A, LDA, |
$ WORK( ITAUQ ), U, LDU, WORK( IWORK ), |
$ WORK( ITAUQ ), U, LDU, WORK( IWORK ), |
Line 1076
|
Line 1172
|
LDWRKU = LDA |
LDWRKU = LDA |
IR = IU + LDWRKU*N |
IR = IU + LDWRKU*N |
LDWRKR = LDA |
LDWRKR = LDA |
ELSE IF( LWORK.GE.WRKBL+( LDA+N )*N ) THEN |
ELSE IF( LWORK.GE.WRKBL+( LDA + N )*N ) THEN |
* |
* |
* WORK(IU) is LDA by N and WORK(IR) is N by N |
* WORK(IU) is LDA by N and WORK(IR) is N by N |
* |
* |
Line 1095
|
Line 1191
|
IWORK = ITAU + N |
IWORK = ITAU + N |
* |
* |
* Compute A=Q*R |
* Compute A=Q*R |
* (Workspace: need 2*N*N+2*N, prefer 2*N*N+N+N*NB) |
* (Workspace: need 2*N*N + 2*N, prefer 2*N*N + N + N*NB) |
* |
* |
CALL DGEQRF( M, N, A, LDA, WORK( ITAU ), |
CALL DGEQRF( M, N, A, LDA, WORK( ITAU ), |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
Line 1108
|
Line 1204
|
$ WORK( IU+1 ), LDWRKU ) |
$ WORK( IU+1 ), LDWRKU ) |
* |
* |
* Generate Q in A |
* Generate Q in A |
* (Workspace: need 2*N*N+2*N, prefer 2*N*N+N+N*NB) |
* (Workspace: need 2*N*N + 2*N, prefer 2*N*N + N + N*NB) |
* |
* |
CALL DORGQR( M, N, N, A, LDA, WORK( ITAU ), |
CALL DORGQR( M, N, N, A, LDA, WORK( ITAU ), |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
Line 1119
|
Line 1215
|
* |
* |
* Bidiagonalize R in WORK(IU), copying result to |
* Bidiagonalize R in WORK(IU), copying result to |
* WORK(IR) |
* WORK(IR) |
* (Workspace: need 2*N*N+4*N, |
* (Workspace: need 2*N*N + 4*N, |
* prefer 2*N*N+3*N+2*N*NB) |
* prefer 2*N*N+3*N+2*N*NB) |
* |
* |
CALL DGEBRD( N, N, WORK( IU ), LDWRKU, S, |
CALL DGEBRD( N, N, WORK( IU ), LDWRKU, S, |
Line 1130
|
Line 1226
|
$ WORK( IR ), LDWRKR ) |
$ WORK( IR ), LDWRKR ) |
* |
* |
* Generate left bidiagonalizing vectors in WORK(IU) |
* Generate left bidiagonalizing vectors in WORK(IU) |
* (Workspace: need 2*N*N+4*N, prefer 2*N*N+3*N+N*NB) |
* (Workspace: need 2*N*N + 4*N, prefer 2*N*N + 3*N + N*NB) |
* |
* |
CALL DORGBR( 'Q', N, N, N, WORK( IU ), LDWRKU, |
CALL DORGBR( 'Q', N, N, N, WORK( IU ), LDWRKU, |
$ WORK( ITAUQ ), WORK( IWORK ), |
$ WORK( ITAUQ ), WORK( IWORK ), |
$ LWORK-IWORK+1, IERR ) |
$ LWORK-IWORK+1, IERR ) |
* |
* |
* Generate right bidiagonalizing vectors in WORK(IR) |
* Generate right bidiagonalizing vectors in WORK(IR) |
* (Workspace: need 2*N*N+4*N-1, |
* (Workspace: need 2*N*N + 4*N-1, |
* prefer 2*N*N+3*N+(N-1)*NB) |
* prefer 2*N*N+3*N+(N-1)*NB) |
* |
* |
CALL DORGBR( 'P', N, N, N, WORK( IR ), LDWRKR, |
CALL DORGBR( 'P', N, N, N, WORK( IR ), LDWRKR, |
Line 1148
|
Line 1244
|
* Perform bidiagonal QR iteration, computing left |
* Perform bidiagonal QR iteration, computing left |
* singular vectors of R in WORK(IU) and computing |
* singular vectors of R in WORK(IU) and computing |
* right singular vectors of R in WORK(IR) |
* right singular vectors of R in WORK(IR) |
* (Workspace: need 2*N*N+BDSPAC) |
* (Workspace: need 2*N*N + BDSPAC) |
* |
* |
CALL DBDSQR( 'U', N, N, N, 0, S, WORK( IE ), |
CALL DBDSQR( 'U', N, N, N, 0, S, WORK( IE ), |
$ WORK( IR ), LDWRKR, WORK( IU ), |
$ WORK( IR ), LDWRKR, WORK( IU ), |
Line 1175
|
Line 1271
|
IWORK = ITAU + N |
IWORK = ITAU + N |
* |
* |
* Compute A=Q*R, copying result to U |
* Compute A=Q*R, copying result to U |
* (Workspace: need 2*N, prefer N+N*NB) |
* (Workspace: need 2*N, prefer N + N*NB) |
* |
* |
CALL DGEQRF( M, N, A, LDA, WORK( ITAU ), |
CALL DGEQRF( M, N, A, LDA, WORK( ITAU ), |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
CALL DLACPY( 'L', M, N, A, LDA, U, LDU ) |
CALL DLACPY( 'L', M, N, A, LDA, U, LDU ) |
* |
* |
* Generate Q in U |
* Generate Q in U |
* (Workspace: need 2*N, prefer N+N*NB) |
* (Workspace: need 2*N, prefer N + N*NB) |
* |
* |
CALL DORGQR( M, N, N, U, LDU, WORK( ITAU ), |
CALL DORGQR( M, N, N, U, LDU, WORK( ITAU ), |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
Line 1193
|
Line 1289
|
* |
* |
* Zero out below R in A |
* Zero out below R in A |
* |
* |
CALL DLASET( 'L', N-1, N-1, ZERO, ZERO, A( 2, 1 ), |
IF( N .GT. 1 ) THEN |
$ LDA ) |
CALL DLASET( 'L', N-1, N-1, ZERO, ZERO, |
|
$ A( 2, 1 ), LDA ) |
|
END IF |
* |
* |
* Bidiagonalize R in A |
* Bidiagonalize R in A |
* (Workspace: need 4*N, prefer 3*N+2*N*NB) |
* (Workspace: need 4*N, prefer 3*N + 2*N*NB) |
* |
* |
CALL DGEBRD( N, N, A, LDA, S, WORK( IE ), |
CALL DGEBRD( N, N, A, LDA, S, WORK( IE ), |
$ WORK( ITAUQ ), WORK( ITAUP ), |
$ WORK( ITAUQ ), WORK( ITAUP ), |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
* |
* |
* Multiply Q in U by left vectors bidiagonalizing R |
* Multiply Q in U by left vectors bidiagonalizing R |
* (Workspace: need 3*N+M, prefer 3*N+M*NB) |
* (Workspace: need 3*N + M, prefer 3*N + M*NB) |
* |
* |
CALL DORMBR( 'Q', 'R', 'N', M, N, N, A, LDA, |
CALL DORMBR( 'Q', 'R', 'N', M, N, N, A, LDA, |
$ WORK( ITAUQ ), U, LDU, WORK( IWORK ), |
$ WORK( ITAUQ ), U, LDU, WORK( IWORK ), |
$ LWORK-IWORK+1, IERR ) |
$ LWORK-IWORK+1, IERR ) |
* |
* |
* Generate right vectors bidiagonalizing R in A |
* Generate right vectors bidiagonalizing R in A |
* (Workspace: need 4*N-1, prefer 3*N+(N-1)*NB) |
* (Workspace: need 4*N-1, prefer 3*N + (N-1)*NB) |
* |
* |
CALL DORGBR( 'P', N, N, N, A, LDA, WORK( ITAUP ), |
CALL DORGBR( 'P', N, N, N, A, LDA, WORK( ITAUP ), |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
Line 1255
|
Line 1353
|
IWORK = ITAU + N |
IWORK = ITAU + N |
* |
* |
* Compute A=Q*R |
* Compute A=Q*R |
* (Workspace: need N*N+2*N, prefer N*N+N+N*NB) |
* (Workspace: need N*N + 2*N, prefer N*N + N + N*NB) |
* |
* |
CALL DGEQRF( M, N, A, LDA, WORK( ITAU ), |
CALL DGEQRF( M, N, A, LDA, WORK( ITAU ), |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
Line 1268
|
Line 1366
|
$ WORK( IU+1 ), LDWRKU ) |
$ WORK( IU+1 ), LDWRKU ) |
* |
* |
* Generate Q in A |
* Generate Q in A |
* (Workspace: need N*N+2*N, prefer N*N+N+N*NB) |
* (Workspace: need N*N + 2*N, prefer N*N + N + N*NB) |
* |
* |
CALL DORGQR( M, N, N, A, LDA, WORK( ITAU ), |
CALL DORGQR( M, N, N, A, LDA, WORK( ITAU ), |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
Line 1278
|
Line 1376
|
IWORK = ITAUP + N |
IWORK = ITAUP + N |
* |
* |
* Bidiagonalize R in WORK(IU), copying result to VT |
* Bidiagonalize R in WORK(IU), copying result to VT |
* (Workspace: need N*N+4*N, prefer N*N+3*N+2*N*NB) |
* (Workspace: need N*N + 4*N, prefer N*N + 3*N + 2*N*NB) |
* |
* |
CALL DGEBRD( N, N, WORK( IU ), LDWRKU, S, |
CALL DGEBRD( N, N, WORK( IU ), LDWRKU, S, |
$ WORK( IE ), WORK( ITAUQ ), |
$ WORK( IE ), WORK( ITAUQ ), |
Line 1288
|
Line 1386
|
$ LDVT ) |
$ LDVT ) |
* |
* |
* Generate left bidiagonalizing vectors in WORK(IU) |
* Generate left bidiagonalizing vectors in WORK(IU) |
* (Workspace: need N*N+4*N, prefer N*N+3*N+N*NB) |
* (Workspace: need N*N + 4*N, prefer N*N + 3*N + N*NB) |
* |
* |
CALL DORGBR( 'Q', N, N, N, WORK( IU ), LDWRKU, |
CALL DORGBR( 'Q', N, N, N, WORK( IU ), LDWRKU, |
$ WORK( ITAUQ ), WORK( IWORK ), |
$ WORK( ITAUQ ), WORK( IWORK ), |
$ LWORK-IWORK+1, IERR ) |
$ LWORK-IWORK+1, IERR ) |
* |
* |
* Generate right bidiagonalizing vectors in VT |
* Generate right bidiagonalizing vectors in VT |
* (Workspace: need N*N+4*N-1, |
* (Workspace: need N*N + 4*N-1, |
* prefer N*N+3*N+(N-1)*NB) |
* prefer N*N+3*N+(N-1)*NB) |
* |
* |
CALL DORGBR( 'P', N, N, N, VT, LDVT, WORK( ITAUP ), |
CALL DORGBR( 'P', N, N, N, VT, LDVT, WORK( ITAUP ), |
Line 1305
|
Line 1403
|
* Perform bidiagonal QR iteration, computing left |
* Perform bidiagonal QR iteration, computing left |
* singular vectors of R in WORK(IU) and computing |
* singular vectors of R in WORK(IU) and computing |
* right singular vectors of R in VT |
* right singular vectors of R in VT |
* (Workspace: need N*N+BDSPAC) |
* (Workspace: need N*N + BDSPAC) |
* |
* |
CALL DBDSQR( 'U', N, N, N, 0, S, WORK( IE ), VT, |
CALL DBDSQR( 'U', N, N, N, 0, S, WORK( IE ), VT, |
$ LDVT, WORK( IU ), LDWRKU, DUM, 1, |
$ LDVT, WORK( IU ), LDWRKU, DUM, 1, |
Line 1326
|
Line 1424
|
IWORK = ITAU + N |
IWORK = ITAU + N |
* |
* |
* Compute A=Q*R, copying result to U |
* Compute A=Q*R, copying result to U |
* (Workspace: need 2*N, prefer N+N*NB) |
* (Workspace: need 2*N, prefer N + N*NB) |
* |
* |
CALL DGEQRF( M, N, A, LDA, WORK( ITAU ), |
CALL DGEQRF( M, N, A, LDA, WORK( ITAU ), |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
CALL DLACPY( 'L', M, N, A, LDA, U, LDU ) |
CALL DLACPY( 'L', M, N, A, LDA, U, LDU ) |
* |
* |
* Generate Q in U |
* Generate Q in U |
* (Workspace: need 2*N, prefer N+N*NB) |
* (Workspace: need 2*N, prefer N + N*NB) |
* |
* |
CALL DORGQR( M, N, N, U, LDU, WORK( ITAU ), |
CALL DORGQR( M, N, N, U, LDU, WORK( ITAU ), |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
Line 1350
|
Line 1448
|
IWORK = ITAUP + N |
IWORK = ITAUP + N |
* |
* |
* Bidiagonalize R in VT |
* Bidiagonalize R in VT |
* (Workspace: need 4*N, prefer 3*N+2*N*NB) |
* (Workspace: need 4*N, prefer 3*N + 2*N*NB) |
* |
* |
CALL DGEBRD( N, N, VT, LDVT, S, WORK( IE ), |
CALL DGEBRD( N, N, VT, LDVT, S, WORK( IE ), |
$ WORK( ITAUQ ), WORK( ITAUP ), |
$ WORK( ITAUQ ), WORK( ITAUP ), |
Line 1358
|
Line 1456
|
* |
* |
* Multiply Q in U by left bidiagonalizing vectors |
* Multiply Q in U by left bidiagonalizing vectors |
* in VT |
* in VT |
* (Workspace: need 3*N+M, prefer 3*N+M*NB) |
* (Workspace: need 3*N + M, prefer 3*N + M*NB) |
* |
* |
CALL DORMBR( 'Q', 'R', 'N', M, N, N, VT, LDVT, |
CALL DORMBR( 'Q', 'R', 'N', M, N, N, VT, LDVT, |
$ WORK( ITAUQ ), U, LDU, WORK( IWORK ), |
$ WORK( ITAUQ ), U, LDU, WORK( IWORK ), |
$ LWORK-IWORK+1, IERR ) |
$ LWORK-IWORK+1, IERR ) |
* |
* |
* Generate right bidiagonalizing vectors in VT |
* Generate right bidiagonalizing vectors in VT |
* (Workspace: need 4*N-1, prefer 3*N+(N-1)*NB) |
* (Workspace: need 4*N-1, prefer 3*N + (N-1)*NB) |
* |
* |
CALL DORGBR( 'P', N, N, N, VT, LDVT, WORK( ITAUP ), |
CALL DORGBR( 'P', N, N, N, VT, LDVT, WORK( ITAUP ), |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
Line 1412
|
Line 1510
|
IWORK = ITAU + N |
IWORK = ITAU + N |
* |
* |
* Compute A=Q*R, copying result to U |
* Compute A=Q*R, copying result to U |
* (Workspace: need N*N+2*N, prefer N*N+N+N*NB) |
* (Workspace: need N*N + 2*N, prefer N*N + N + N*NB) |
* |
* |
CALL DGEQRF( M, N, A, LDA, WORK( ITAU ), |
CALL DGEQRF( M, N, A, LDA, WORK( ITAU ), |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
Line 1426
|
Line 1524
|
$ WORK( IR+1 ), LDWRKR ) |
$ WORK( IR+1 ), LDWRKR ) |
* |
* |
* Generate Q in U |
* Generate Q in U |
* (Workspace: need N*N+N+M, prefer N*N+N+M*NB) |
* (Workspace: need N*N + N + M, prefer N*N + N + M*NB) |
* |
* |
CALL DORGQR( M, M, N, U, LDU, WORK( ITAU ), |
CALL DORGQR( M, M, N, U, LDU, WORK( ITAU ), |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
Line 1436
|
Line 1534
|
IWORK = ITAUP + N |
IWORK = ITAUP + N |
* |
* |
* Bidiagonalize R in WORK(IR) |
* Bidiagonalize R in WORK(IR) |
* (Workspace: need N*N+4*N, prefer N*N+3*N+2*N*NB) |
* (Workspace: need N*N + 4*N, prefer N*N + 3*N + 2*N*NB) |
* |
* |
CALL DGEBRD( N, N, WORK( IR ), LDWRKR, S, |
CALL DGEBRD( N, N, WORK( IR ), LDWRKR, S, |
$ WORK( IE ), WORK( ITAUQ ), |
$ WORK( IE ), WORK( ITAUQ ), |
Line 1444
|
Line 1542
|
$ LWORK-IWORK+1, IERR ) |
$ LWORK-IWORK+1, IERR ) |
* |
* |
* Generate left bidiagonalizing vectors in WORK(IR) |
* Generate left bidiagonalizing vectors in WORK(IR) |
* (Workspace: need N*N+4*N, prefer N*N+3*N+N*NB) |
* (Workspace: need N*N + 4*N, prefer N*N + 3*N + N*NB) |
* |
* |
CALL DORGBR( 'Q', N, N, N, WORK( IR ), LDWRKR, |
CALL DORGBR( 'Q', N, N, N, WORK( IR ), LDWRKR, |
$ WORK( ITAUQ ), WORK( IWORK ), |
$ WORK( ITAUQ ), WORK( IWORK ), |
Line 1453
|
Line 1551
|
* |
* |
* Perform bidiagonal QR iteration, computing left |
* Perform bidiagonal QR iteration, computing left |
* singular vectors of R in WORK(IR) |
* singular vectors of R in WORK(IR) |
* (Workspace: need N*N+BDSPAC) |
* (Workspace: need N*N + BDSPAC) |
* |
* |
CALL DBDSQR( 'U', N, 0, N, 0, S, WORK( IE ), DUM, |
CALL DBDSQR( 'U', N, 0, N, 0, S, WORK( IE ), DUM, |
$ 1, WORK( IR ), LDWRKR, DUM, 1, |
$ 1, WORK( IR ), LDWRKR, DUM, 1, |
Line 1478
|
Line 1576
|
IWORK = ITAU + N |
IWORK = ITAU + N |
* |
* |
* Compute A=Q*R, copying result to U |
* Compute A=Q*R, copying result to U |
* (Workspace: need 2*N, prefer N+N*NB) |
* (Workspace: need 2*N, prefer N + N*NB) |
* |
* |
CALL DGEQRF( M, N, A, LDA, WORK( ITAU ), |
CALL DGEQRF( M, N, A, LDA, WORK( ITAU ), |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
CALL DLACPY( 'L', M, N, A, LDA, U, LDU ) |
CALL DLACPY( 'L', M, N, A, LDA, U, LDU ) |
* |
* |
* Generate Q in U |
* Generate Q in U |
* (Workspace: need N+M, prefer N+M*NB) |
* (Workspace: need N + M, prefer N + M*NB) |
* |
* |
CALL DORGQR( M, M, N, U, LDU, WORK( ITAU ), |
CALL DORGQR( M, M, N, U, LDU, WORK( ITAU ), |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
Line 1496
|
Line 1594
|
* |
* |
* Zero out below R in A |
* Zero out below R in A |
* |
* |
CALL DLASET( 'L', N-1, N-1, ZERO, ZERO, A( 2, 1 ), |
IF( N .GT. 1 ) THEN |
$ LDA ) |
CALL DLASET( 'L', N-1, N-1, ZERO, ZERO, |
|
$ A( 2, 1 ), LDA ) |
|
END IF |
* |
* |
* Bidiagonalize R in A |
* Bidiagonalize R in A |
* (Workspace: need 4*N, prefer 3*N+2*N*NB) |
* (Workspace: need 4*N, prefer 3*N + 2*N*NB) |
* |
* |
CALL DGEBRD( N, N, A, LDA, S, WORK( IE ), |
CALL DGEBRD( N, N, A, LDA, S, WORK( IE ), |
$ WORK( ITAUQ ), WORK( ITAUP ), |
$ WORK( ITAUQ ), WORK( ITAUP ), |
Line 1508
|
Line 1608
|
* |
* |
* Multiply Q in U by left bidiagonalizing vectors |
* Multiply Q in U by left bidiagonalizing vectors |
* in A |
* in A |
* (Workspace: need 3*N+M, prefer 3*N+M*NB) |
* (Workspace: need 3*N + M, prefer 3*N + M*NB) |
* |
* |
CALL DORMBR( 'Q', 'R', 'N', M, N, N, A, LDA, |
CALL DORMBR( 'Q', 'R', 'N', M, N, N, A, LDA, |
$ WORK( ITAUQ ), U, LDU, WORK( IWORK ), |
$ WORK( ITAUQ ), U, LDU, WORK( IWORK ), |
Line 1543
|
Line 1643
|
LDWRKU = LDA |
LDWRKU = LDA |
IR = IU + LDWRKU*N |
IR = IU + LDWRKU*N |
LDWRKR = LDA |
LDWRKR = LDA |
ELSE IF( LWORK.GE.WRKBL+( LDA+N )*N ) THEN |
ELSE IF( LWORK.GE.WRKBL+( LDA + N )*N ) THEN |
* |
* |
* WORK(IU) is LDA by N and WORK(IR) is N by N |
* WORK(IU) is LDA by N and WORK(IR) is N by N |
* |
* |
Line 1562
|
Line 1662
|
IWORK = ITAU + N |
IWORK = ITAU + N |
* |
* |
* Compute A=Q*R, copying result to U |
* Compute A=Q*R, copying result to U |
* (Workspace: need 2*N*N+2*N, prefer 2*N*N+N+N*NB) |
* (Workspace: need 2*N*N + 2*N, prefer 2*N*N + N + N*NB) |
* |
* |
CALL DGEQRF( M, N, A, LDA, WORK( ITAU ), |
CALL DGEQRF( M, N, A, LDA, WORK( ITAU ), |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
CALL DLACPY( 'L', M, N, A, LDA, U, LDU ) |
CALL DLACPY( 'L', M, N, A, LDA, U, LDU ) |
* |
* |
* Generate Q in U |
* Generate Q in U |
* (Workspace: need 2*N*N+N+M, prefer 2*N*N+N+M*NB) |
* (Workspace: need 2*N*N + N + M, prefer 2*N*N + N + M*NB) |
* |
* |
CALL DORGQR( M, M, N, U, LDU, WORK( ITAU ), |
CALL DORGQR( M, M, N, U, LDU, WORK( ITAU ), |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
Line 1587
|
Line 1687
|
* |
* |
* Bidiagonalize R in WORK(IU), copying result to |
* Bidiagonalize R in WORK(IU), copying result to |
* WORK(IR) |
* WORK(IR) |
* (Workspace: need 2*N*N+4*N, |
* (Workspace: need 2*N*N + 4*N, |
* prefer 2*N*N+3*N+2*N*NB) |
* prefer 2*N*N+3*N+2*N*NB) |
* |
* |
CALL DGEBRD( N, N, WORK( IU ), LDWRKU, S, |
CALL DGEBRD( N, N, WORK( IU ), LDWRKU, S, |
Line 1598
|
Line 1698
|
$ WORK( IR ), LDWRKR ) |
$ WORK( IR ), LDWRKR ) |
* |
* |
* Generate left bidiagonalizing vectors in WORK(IU) |
* Generate left bidiagonalizing vectors in WORK(IU) |
* (Workspace: need 2*N*N+4*N, prefer 2*N*N+3*N+N*NB) |
* (Workspace: need 2*N*N + 4*N, prefer 2*N*N + 3*N + N*NB) |
* |
* |
CALL DORGBR( 'Q', N, N, N, WORK( IU ), LDWRKU, |
CALL DORGBR( 'Q', N, N, N, WORK( IU ), LDWRKU, |
$ WORK( ITAUQ ), WORK( IWORK ), |
$ WORK( ITAUQ ), WORK( IWORK ), |
$ LWORK-IWORK+1, IERR ) |
$ LWORK-IWORK+1, IERR ) |
* |
* |
* Generate right bidiagonalizing vectors in WORK(IR) |
* Generate right bidiagonalizing vectors in WORK(IR) |
* (Workspace: need 2*N*N+4*N-1, |
* (Workspace: need 2*N*N + 4*N-1, |
* prefer 2*N*N+3*N+(N-1)*NB) |
* prefer 2*N*N+3*N+(N-1)*NB) |
* |
* |
CALL DORGBR( 'P', N, N, N, WORK( IR ), LDWRKR, |
CALL DORGBR( 'P', N, N, N, WORK( IR ), LDWRKR, |
Line 1616
|
Line 1716
|
* Perform bidiagonal QR iteration, computing left |
* Perform bidiagonal QR iteration, computing left |
* singular vectors of R in WORK(IU) and computing |
* singular vectors of R in WORK(IU) and computing |
* right singular vectors of R in WORK(IR) |
* right singular vectors of R in WORK(IR) |
* (Workspace: need 2*N*N+BDSPAC) |
* (Workspace: need 2*N*N + BDSPAC) |
* |
* |
CALL DBDSQR( 'U', N, N, N, 0, S, WORK( IE ), |
CALL DBDSQR( 'U', N, N, N, 0, S, WORK( IE ), |
$ WORK( IR ), LDWRKR, WORK( IU ), |
$ WORK( IR ), LDWRKR, WORK( IU ), |
Line 1646
|
Line 1746
|
IWORK = ITAU + N |
IWORK = ITAU + N |
* |
* |
* Compute A=Q*R, copying result to U |
* Compute A=Q*R, copying result to U |
* (Workspace: need 2*N, prefer N+N*NB) |
* (Workspace: need 2*N, prefer N + N*NB) |
* |
* |
CALL DGEQRF( M, N, A, LDA, WORK( ITAU ), |
CALL DGEQRF( M, N, A, LDA, WORK( ITAU ), |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
CALL DLACPY( 'L', M, N, A, LDA, U, LDU ) |
CALL DLACPY( 'L', M, N, A, LDA, U, LDU ) |
* |
* |
* Generate Q in U |
* Generate Q in U |
* (Workspace: need N+M, prefer N+M*NB) |
* (Workspace: need N + M, prefer N + M*NB) |
* |
* |
CALL DORGQR( M, M, N, U, LDU, WORK( ITAU ), |
CALL DORGQR( M, M, N, U, LDU, WORK( ITAU ), |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
Line 1664
|
Line 1764
|
* |
* |
* Zero out below R in A |
* Zero out below R in A |
* |
* |
CALL DLASET( 'L', N-1, N-1, ZERO, ZERO, A( 2, 1 ), |
IF( N .GT. 1 ) THEN |
$ LDA ) |
CALL DLASET( 'L', N-1, N-1, ZERO, ZERO, |
|
$ A( 2, 1 ), LDA ) |
|
END IF |
* |
* |
* Bidiagonalize R in A |
* Bidiagonalize R in A |
* (Workspace: need 4*N, prefer 3*N+2*N*NB) |
* (Workspace: need 4*N, prefer 3*N + 2*N*NB) |
* |
* |
CALL DGEBRD( N, N, A, LDA, S, WORK( IE ), |
CALL DGEBRD( N, N, A, LDA, S, WORK( IE ), |
$ WORK( ITAUQ ), WORK( ITAUP ), |
$ WORK( ITAUQ ), WORK( ITAUP ), |
Line 1676
|
Line 1778
|
* |
* |
* Multiply Q in U by left bidiagonalizing vectors |
* Multiply Q in U by left bidiagonalizing vectors |
* in A |
* in A |
* (Workspace: need 3*N+M, prefer 3*N+M*NB) |
* (Workspace: need 3*N + M, prefer 3*N + M*NB) |
* |
* |
CALL DORMBR( 'Q', 'R', 'N', M, N, N, A, LDA, |
CALL DORMBR( 'Q', 'R', 'N', M, N, N, A, LDA, |
$ WORK( ITAUQ ), U, LDU, WORK( IWORK ), |
$ WORK( ITAUQ ), U, LDU, WORK( IWORK ), |
$ LWORK-IWORK+1, IERR ) |
$ LWORK-IWORK+1, IERR ) |
* |
* |
* Generate right bidiagonalizing vectors in A |
* Generate right bidiagonalizing vectors in A |
* (Workspace: need 4*N-1, prefer 3*N+(N-1)*NB) |
* (Workspace: need 4*N-1, prefer 3*N + (N-1)*NB) |
* |
* |
CALL DORGBR( 'P', N, N, N, A, LDA, WORK( ITAUP ), |
CALL DORGBR( 'P', N, N, N, A, LDA, WORK( ITAUP ), |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
Line 1727
|
Line 1829
|
IWORK = ITAU + N |
IWORK = ITAU + N |
* |
* |
* Compute A=Q*R, copying result to U |
* Compute A=Q*R, copying result to U |
* (Workspace: need N*N+2*N, prefer N*N+N+N*NB) |
* (Workspace: need N*N + 2*N, prefer N*N + N + N*NB) |
* |
* |
CALL DGEQRF( M, N, A, LDA, WORK( ITAU ), |
CALL DGEQRF( M, N, A, LDA, WORK( ITAU ), |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
CALL DLACPY( 'L', M, N, A, LDA, U, LDU ) |
CALL DLACPY( 'L', M, N, A, LDA, U, LDU ) |
* |
* |
* Generate Q in U |
* Generate Q in U |
* (Workspace: need N*N+N+M, prefer N*N+N+M*NB) |
* (Workspace: need N*N + N + M, prefer N*N + N + M*NB) |
* |
* |
CALL DORGQR( M, M, N, U, LDU, WORK( ITAU ), |
CALL DORGQR( M, M, N, U, LDU, WORK( ITAU ), |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
Line 1751
|
Line 1853
|
IWORK = ITAUP + N |
IWORK = ITAUP + N |
* |
* |
* Bidiagonalize R in WORK(IU), copying result to VT |
* Bidiagonalize R in WORK(IU), copying result to VT |
* (Workspace: need N*N+4*N, prefer N*N+3*N+2*N*NB) |
* (Workspace: need N*N + 4*N, prefer N*N + 3*N + 2*N*NB) |
* |
* |
CALL DGEBRD( N, N, WORK( IU ), LDWRKU, S, |
CALL DGEBRD( N, N, WORK( IU ), LDWRKU, S, |
$ WORK( IE ), WORK( ITAUQ ), |
$ WORK( IE ), WORK( ITAUQ ), |
Line 1761
|
Line 1863
|
$ LDVT ) |
$ LDVT ) |
* |
* |
* Generate left bidiagonalizing vectors in WORK(IU) |
* Generate left bidiagonalizing vectors in WORK(IU) |
* (Workspace: need N*N+4*N, prefer N*N+3*N+N*NB) |
* (Workspace: need N*N + 4*N, prefer N*N + 3*N + N*NB) |
* |
* |
CALL DORGBR( 'Q', N, N, N, WORK( IU ), LDWRKU, |
CALL DORGBR( 'Q', N, N, N, WORK( IU ), LDWRKU, |
$ WORK( ITAUQ ), WORK( IWORK ), |
$ WORK( ITAUQ ), WORK( IWORK ), |
$ LWORK-IWORK+1, IERR ) |
$ LWORK-IWORK+1, IERR ) |
* |
* |
* Generate right bidiagonalizing vectors in VT |
* Generate right bidiagonalizing vectors in VT |
* (Workspace: need N*N+4*N-1, |
* (Workspace: need N*N + 4*N-1, |
* prefer N*N+3*N+(N-1)*NB) |
* prefer N*N+3*N+(N-1)*NB) |
* |
* |
CALL DORGBR( 'P', N, N, N, VT, LDVT, WORK( ITAUP ), |
CALL DORGBR( 'P', N, N, N, VT, LDVT, WORK( ITAUP ), |
Line 1778
|
Line 1880
|
* Perform bidiagonal QR iteration, computing left |
* Perform bidiagonal QR iteration, computing left |
* singular vectors of R in WORK(IU) and computing |
* singular vectors of R in WORK(IU) and computing |
* right singular vectors of R in VT |
* right singular vectors of R in VT |
* (Workspace: need N*N+BDSPAC) |
* (Workspace: need N*N + BDSPAC) |
* |
* |
CALL DBDSQR( 'U', N, N, N, 0, S, WORK( IE ), VT, |
CALL DBDSQR( 'U', N, N, N, 0, S, WORK( IE ), VT, |
$ LDVT, WORK( IU ), LDWRKU, DUM, 1, |
$ LDVT, WORK( IU ), LDWRKU, DUM, 1, |
Line 1803
|
Line 1905
|
IWORK = ITAU + N |
IWORK = ITAU + N |
* |
* |
* Compute A=Q*R, copying result to U |
* Compute A=Q*R, copying result to U |
* (Workspace: need 2*N, prefer N+N*NB) |
* (Workspace: need 2*N, prefer N + N*NB) |
* |
* |
CALL DGEQRF( M, N, A, LDA, WORK( ITAU ), |
CALL DGEQRF( M, N, A, LDA, WORK( ITAU ), |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
CALL DLACPY( 'L', M, N, A, LDA, U, LDU ) |
CALL DLACPY( 'L', M, N, A, LDA, U, LDU ) |
* |
* |
* Generate Q in U |
* Generate Q in U |
* (Workspace: need N+M, prefer N+M*NB) |
* (Workspace: need N + M, prefer N + M*NB) |
* |
* |
CALL DORGQR( M, M, N, U, LDU, WORK( ITAU ), |
CALL DORGQR( M, M, N, U, LDU, WORK( ITAU ), |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
Line 1827
|
Line 1929
|
IWORK = ITAUP + N |
IWORK = ITAUP + N |
* |
* |
* Bidiagonalize R in VT |
* Bidiagonalize R in VT |
* (Workspace: need 4*N, prefer 3*N+2*N*NB) |
* (Workspace: need 4*N, prefer 3*N + 2*N*NB) |
* |
* |
CALL DGEBRD( N, N, VT, LDVT, S, WORK( IE ), |
CALL DGEBRD( N, N, VT, LDVT, S, WORK( IE ), |
$ WORK( ITAUQ ), WORK( ITAUP ), |
$ WORK( ITAUQ ), WORK( ITAUP ), |
Line 1835
|
Line 1937
|
* |
* |
* Multiply Q in U by left bidiagonalizing vectors |
* Multiply Q in U by left bidiagonalizing vectors |
* in VT |
* in VT |
* (Workspace: need 3*N+M, prefer 3*N+M*NB) |
* (Workspace: need 3*N + M, prefer 3*N + M*NB) |
* |
* |
CALL DORMBR( 'Q', 'R', 'N', M, N, N, VT, LDVT, |
CALL DORMBR( 'Q', 'R', 'N', M, N, N, VT, LDVT, |
$ WORK( ITAUQ ), U, LDU, WORK( IWORK ), |
$ WORK( ITAUQ ), U, LDU, WORK( IWORK ), |
$ LWORK-IWORK+1, IERR ) |
$ LWORK-IWORK+1, IERR ) |
* |
* |
* Generate right bidiagonalizing vectors in VT |
* Generate right bidiagonalizing vectors in VT |
* (Workspace: need 4*N-1, prefer 3*N+(N-1)*NB) |
* (Workspace: need 4*N-1, prefer 3*N + (N-1)*NB) |
* |
* |
CALL DORGBR( 'P', N, N, N, VT, LDVT, WORK( ITAUP ), |
CALL DORGBR( 'P', N, N, N, VT, LDVT, WORK( ITAUP ), |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
Line 1876
|
Line 1978
|
IWORK = ITAUP + N |
IWORK = ITAUP + N |
* |
* |
* Bidiagonalize A |
* Bidiagonalize A |
* (Workspace: need 3*N+M, prefer 3*N+(M+N)*NB) |
* (Workspace: need 3*N + M, prefer 3*N + (M + N)*NB) |
* |
* |
CALL DGEBRD( M, N, A, LDA, S, WORK( IE ), WORK( ITAUQ ), |
CALL DGEBRD( M, N, A, LDA, S, WORK( IE ), WORK( ITAUQ ), |
$ WORK( ITAUP ), WORK( IWORK ), LWORK-IWORK+1, |
$ WORK( ITAUP ), WORK( IWORK ), LWORK-IWORK+1, |
Line 1885
|
Line 1987
|
* |
* |
* If left singular vectors desired in U, copy result to U |
* If left singular vectors desired in U, copy result to U |
* and generate left bidiagonalizing vectors in U |
* and generate left bidiagonalizing vectors in U |
* (Workspace: need 3*N+NCU, prefer 3*N+NCU*NB) |
* (Workspace: need 3*N + NCU, prefer 3*N + NCU*NB) |
* |
* |
CALL DLACPY( 'L', M, N, A, LDA, U, LDU ) |
CALL DLACPY( 'L', M, N, A, LDA, U, LDU ) |
IF( WNTUS ) |
IF( WNTUS ) |
Line 1899
|
Line 2001
|
* |
* |
* If right singular vectors desired in VT, copy result to |
* If right singular vectors desired in VT, copy result to |
* VT and generate right bidiagonalizing vectors in VT |
* VT and generate right bidiagonalizing vectors in VT |
* (Workspace: need 4*N-1, prefer 3*N+(N-1)*NB) |
* (Workspace: need 4*N-1, prefer 3*N + (N-1)*NB) |
* |
* |
CALL DLACPY( 'U', N, N, A, LDA, VT, LDVT ) |
CALL DLACPY( 'U', N, N, A, LDA, VT, LDVT ) |
CALL DORGBR( 'P', N, N, N, VT, LDVT, WORK( ITAUP ), |
CALL DORGBR( 'P', N, N, N, VT, LDVT, WORK( ITAUP ), |
Line 1909
|
Line 2011
|
* |
* |
* If left singular vectors desired in A, generate left |
* If left singular vectors desired in A, generate left |
* bidiagonalizing vectors in A |
* bidiagonalizing vectors in A |
* (Workspace: need 4*N, prefer 3*N+N*NB) |
* (Workspace: need 4*N, prefer 3*N + N*NB) |
* |
* |
CALL DORGBR( 'Q', M, N, N, A, LDA, WORK( ITAUQ ), |
CALL DORGBR( 'Q', M, N, N, A, LDA, WORK( ITAUQ ), |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
Line 1918
|
Line 2020
|
* |
* |
* If right singular vectors desired in A, generate right |
* If right singular vectors desired in A, generate right |
* bidiagonalizing vectors in A |
* bidiagonalizing vectors in A |
* (Workspace: need 4*N-1, prefer 3*N+(N-1)*NB) |
* (Workspace: need 4*N-1, prefer 3*N + (N-1)*NB) |
* |
* |
CALL DORGBR( 'P', N, N, N, A, LDA, WORK( ITAUP ), |
CALL DORGBR( 'P', N, N, N, A, LDA, WORK( ITAUP ), |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
Line 1980
|
Line 2082
|
IWORK = ITAU + M |
IWORK = ITAU + M |
* |
* |
* Compute A=L*Q |
* Compute A=L*Q |
* (Workspace: need 2*M, prefer M+M*NB) |
* (Workspace: need 2*M, prefer M + M*NB) |
* |
* |
CALL DGELQF( M, N, A, LDA, WORK( ITAU ), WORK( IWORK ), |
CALL DGELQF( M, N, A, LDA, WORK( ITAU ), WORK( IWORK ), |
$ LWORK-IWORK+1, IERR ) |
$ LWORK-IWORK+1, IERR ) |
Line 1994
|
Line 2096
|
IWORK = ITAUP + M |
IWORK = ITAUP + M |
* |
* |
* Bidiagonalize L in A |
* Bidiagonalize L in A |
* (Workspace: need 4*M, prefer 3*M+2*M*NB) |
* (Workspace: need 4*M, prefer 3*M + 2*M*NB) |
* |
* |
CALL DGEBRD( M, M, A, LDA, S, WORK( IE ), WORK( ITAUQ ), |
CALL DGEBRD( M, M, A, LDA, S, WORK( IE ), WORK( ITAUQ ), |
$ WORK( ITAUP ), WORK( IWORK ), LWORK-IWORK+1, |
$ WORK( ITAUP ), WORK( IWORK ), LWORK-IWORK+1, |
Line 2002
|
Line 2104
|
IF( WNTUO .OR. WNTUAS ) THEN |
IF( WNTUO .OR. WNTUAS ) THEN |
* |
* |
* If left singular vectors desired, generate Q |
* If left singular vectors desired, generate Q |
* (Workspace: need 4*M, prefer 3*M+M*NB) |
* (Workspace: need 4*M, prefer 3*M + M*NB) |
* |
* |
CALL DORGBR( 'Q', M, M, M, A, LDA, WORK( ITAUQ ), |
CALL DORGBR( 'Q', M, M, M, A, LDA, WORK( ITAUQ ), |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
Line 2035
|
Line 2137
|
* Sufficient workspace for a fast algorithm |
* Sufficient workspace for a fast algorithm |
* |
* |
IR = 1 |
IR = 1 |
IF( LWORK.GE.MAX( WRKBL, LDA*N+M )+LDA*M ) THEN |
IF( LWORK.GE.MAX( WRKBL, LDA*N + M ) + LDA*M ) THEN |
* |
* |
* WORK(IU) is LDA by N and WORK(IR) is LDA by M |
* WORK(IU) is LDA by N and WORK(IR) is LDA by M |
* |
* |
LDWRKU = LDA |
LDWRKU = LDA |
CHUNK = N |
CHUNK = N |
LDWRKR = LDA |
LDWRKR = LDA |
ELSE IF( LWORK.GE.MAX( WRKBL, LDA*N+M )+M*M ) THEN |
ELSE IF( LWORK.GE.MAX( WRKBL, LDA*N + M ) + M*M ) THEN |
* |
* |
* WORK(IU) is LDA by N and WORK(IR) is M by M |
* WORK(IU) is LDA by N and WORK(IR) is M by M |
* |
* |
Line 2061
|
Line 2163
|
IWORK = ITAU + M |
IWORK = ITAU + M |
* |
* |
* Compute A=L*Q |
* Compute A=L*Q |
* (Workspace: need M*M+2*M, prefer M*M+M+M*NB) |
* (Workspace: need M*M + 2*M, prefer M*M + M + M*NB) |
* |
* |
CALL DGELQF( M, N, A, LDA, WORK( ITAU ), |
CALL DGELQF( M, N, A, LDA, WORK( ITAU ), |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
Line 2073
|
Line 2175
|
$ WORK( IR+LDWRKR ), LDWRKR ) |
$ WORK( IR+LDWRKR ), LDWRKR ) |
* |
* |
* Generate Q in A |
* Generate Q in A |
* (Workspace: need M*M+2*M, prefer M*M+M+M*NB) |
* (Workspace: need M*M + 2*M, prefer M*M + M + M*NB) |
* |
* |
CALL DORGLQ( M, N, M, A, LDA, WORK( ITAU ), |
CALL DORGLQ( M, N, M, A, LDA, WORK( ITAU ), |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
Line 2083
|
Line 2185
|
IWORK = ITAUP + M |
IWORK = ITAUP + M |
* |
* |
* Bidiagonalize L in WORK(IR) |
* Bidiagonalize L in WORK(IR) |
* (Workspace: need M*M+4*M, prefer M*M+3*M+2*M*NB) |
* (Workspace: need M*M + 4*M, prefer M*M + 3*M + 2*M*NB) |
* |
* |
CALL DGEBRD( M, M, WORK( IR ), LDWRKR, S, WORK( IE ), |
CALL DGEBRD( M, M, WORK( IR ), LDWRKR, S, WORK( IE ), |
$ WORK( ITAUQ ), WORK( ITAUP ), |
$ WORK( ITAUQ ), WORK( ITAUP ), |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
* |
* |
* Generate right vectors bidiagonalizing L |
* Generate right vectors bidiagonalizing L |
* (Workspace: need M*M+4*M-1, prefer M*M+3*M+(M-1)*NB) |
* (Workspace: need M*M + 4*M-1, prefer M*M + 3*M + (M-1)*NB) |
* |
* |
CALL DORGBR( 'P', M, M, M, WORK( IR ), LDWRKR, |
CALL DORGBR( 'P', M, M, M, WORK( IR ), LDWRKR, |
$ WORK( ITAUP ), WORK( IWORK ), |
$ WORK( ITAUP ), WORK( IWORK ), |
Line 2099
|
Line 2201
|
* |
* |
* Perform bidiagonal QR iteration, computing right |
* Perform bidiagonal QR iteration, computing right |
* singular vectors of L in WORK(IR) |
* singular vectors of L in WORK(IR) |
* (Workspace: need M*M+BDSPAC) |
* (Workspace: need M*M + BDSPAC) |
* |
* |
CALL DBDSQR( 'U', M, M, 0, 0, S, WORK( IE ), |
CALL DBDSQR( 'U', M, M, 0, 0, S, WORK( IE ), |
$ WORK( IR ), LDWRKR, DUM, 1, DUM, 1, |
$ WORK( IR ), LDWRKR, DUM, 1, DUM, 1, |
Line 2108
|
Line 2210
|
* |
* |
* Multiply right singular vectors of L in WORK(IR) by Q |
* Multiply right singular vectors of L in WORK(IR) by Q |
* in A, storing result in WORK(IU) and copying to A |
* in A, storing result in WORK(IU) and copying to A |
* (Workspace: need M*M+2*M, prefer M*M+M*N+M) |
* (Workspace: need M*M + 2*M, prefer M*M + M*N + M) |
* |
* |
DO 30 I = 1, N, CHUNK |
DO 30 I = 1, N, CHUNK |
BLK = MIN( N-I+1, CHUNK ) |
BLK = MIN( N-I+1, CHUNK ) |
Line 2129
|
Line 2231
|
IWORK = ITAUP + M |
IWORK = ITAUP + M |
* |
* |
* Bidiagonalize A |
* Bidiagonalize A |
* (Workspace: need 3*M+N, prefer 3*M+(M+N)*NB) |
* (Workspace: need 3*M + N, prefer 3*M + (M + N)*NB) |
* |
* |
CALL DGEBRD( M, N, A, LDA, S, WORK( IE ), |
CALL DGEBRD( M, N, A, LDA, S, WORK( IE ), |
$ WORK( ITAUQ ), WORK( ITAUP ), |
$ WORK( ITAUQ ), WORK( ITAUP ), |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
* |
* |
* Generate right vectors bidiagonalizing A |
* Generate right vectors bidiagonalizing A |
* (Workspace: need 4*M, prefer 3*M+M*NB) |
* (Workspace: need 4*M, prefer 3*M + M*NB) |
* |
* |
CALL DORGBR( 'P', M, N, M, A, LDA, WORK( ITAUP ), |
CALL DORGBR( 'P', M, N, M, A, LDA, WORK( ITAUP ), |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
Line 2162
|
Line 2264
|
* Sufficient workspace for a fast algorithm |
* Sufficient workspace for a fast algorithm |
* |
* |
IR = 1 |
IR = 1 |
IF( LWORK.GE.MAX( WRKBL, LDA*N+M )+LDA*M ) THEN |
IF( LWORK.GE.MAX( WRKBL, LDA*N + M ) + LDA*M ) THEN |
* |
* |
* WORK(IU) is LDA by N and WORK(IR) is LDA by M |
* WORK(IU) is LDA by N and WORK(IR) is LDA by M |
* |
* |
LDWRKU = LDA |
LDWRKU = LDA |
CHUNK = N |
CHUNK = N |
LDWRKR = LDA |
LDWRKR = LDA |
ELSE IF( LWORK.GE.MAX( WRKBL, LDA*N+M )+M*M ) THEN |
ELSE IF( LWORK.GE.MAX( WRKBL, LDA*N + M ) + M*M ) THEN |
* |
* |
* WORK(IU) is LDA by N and WORK(IR) is M by M |
* WORK(IU) is LDA by N and WORK(IR) is M by M |
* |
* |
Line 2188
|
Line 2290
|
IWORK = ITAU + M |
IWORK = ITAU + M |
* |
* |
* Compute A=L*Q |
* Compute A=L*Q |
* (Workspace: need M*M+2*M, prefer M*M+M+M*NB) |
* (Workspace: need M*M + 2*M, prefer M*M + M + M*NB) |
* |
* |
CALL DGELQF( M, N, A, LDA, WORK( ITAU ), |
CALL DGELQF( M, N, A, LDA, WORK( ITAU ), |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
Line 2200
|
Line 2302
|
$ LDU ) |
$ LDU ) |
* |
* |
* Generate Q in A |
* Generate Q in A |
* (Workspace: need M*M+2*M, prefer M*M+M+M*NB) |
* (Workspace: need M*M + 2*M, prefer M*M + M + M*NB) |
* |
* |
CALL DORGLQ( M, N, M, A, LDA, WORK( ITAU ), |
CALL DORGLQ( M, N, M, A, LDA, WORK( ITAU ), |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
Line 2210
|
Line 2312
|
IWORK = ITAUP + M |
IWORK = ITAUP + M |
* |
* |
* Bidiagonalize L in U, copying result to WORK(IR) |
* Bidiagonalize L in U, copying result to WORK(IR) |
* (Workspace: need M*M+4*M, prefer M*M+3*M+2*M*NB) |
* (Workspace: need M*M + 4*M, prefer M*M + 3*M + 2*M*NB) |
* |
* |
CALL DGEBRD( M, M, U, LDU, S, WORK( IE ), |
CALL DGEBRD( M, M, U, LDU, S, WORK( IE ), |
$ WORK( ITAUQ ), WORK( ITAUP ), |
$ WORK( ITAUQ ), WORK( ITAUP ), |
Line 2218
|
Line 2320
|
CALL DLACPY( 'U', M, M, U, LDU, WORK( IR ), LDWRKR ) |
CALL DLACPY( 'U', M, M, U, LDU, WORK( IR ), LDWRKR ) |
* |
* |
* Generate right vectors bidiagonalizing L in WORK(IR) |
* Generate right vectors bidiagonalizing L in WORK(IR) |
* (Workspace: need M*M+4*M-1, prefer M*M+3*M+(M-1)*NB) |
* (Workspace: need M*M + 4*M-1, prefer M*M + 3*M + (M-1)*NB) |
* |
* |
CALL DORGBR( 'P', M, M, M, WORK( IR ), LDWRKR, |
CALL DORGBR( 'P', M, M, M, WORK( IR ), LDWRKR, |
$ WORK( ITAUP ), WORK( IWORK ), |
$ WORK( ITAUP ), WORK( IWORK ), |
$ LWORK-IWORK+1, IERR ) |
$ LWORK-IWORK+1, IERR ) |
* |
* |
* Generate left vectors bidiagonalizing L in U |
* Generate left vectors bidiagonalizing L in U |
* (Workspace: need M*M+4*M, prefer M*M+3*M+M*NB) |
* (Workspace: need M*M + 4*M, prefer M*M + 3*M + M*NB) |
* |
* |
CALL DORGBR( 'Q', M, M, M, U, LDU, WORK( ITAUQ ), |
CALL DORGBR( 'Q', M, M, M, U, LDU, WORK( ITAUQ ), |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
Line 2234
|
Line 2336
|
* Perform bidiagonal QR iteration, computing left |
* Perform bidiagonal QR iteration, computing left |
* singular vectors of L in U, and computing right |
* singular vectors of L in U, and computing right |
* singular vectors of L in WORK(IR) |
* singular vectors of L in WORK(IR) |
* (Workspace: need M*M+BDSPAC) |
* (Workspace: need M*M + BDSPAC) |
* |
* |
CALL DBDSQR( 'U', M, M, M, 0, S, WORK( IE ), |
CALL DBDSQR( 'U', M, M, M, 0, S, WORK( IE ), |
$ WORK( IR ), LDWRKR, U, LDU, DUM, 1, |
$ WORK( IR ), LDWRKR, U, LDU, DUM, 1, |
Line 2243
|
Line 2345
|
* |
* |
* Multiply right singular vectors of L in WORK(IR) by Q |
* Multiply right singular vectors of L in WORK(IR) by Q |
* in A, storing result in WORK(IU) and copying to A |
* in A, storing result in WORK(IU) and copying to A |
* (Workspace: need M*M+2*M, prefer M*M+M*N+M)) |
* (Workspace: need M*M + 2*M, prefer M*M + M*N + M)) |
* |
* |
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 2262
|
Line 2364
|
IWORK = ITAU + M |
IWORK = ITAU + M |
* |
* |
* Compute A=L*Q |
* Compute A=L*Q |
* (Workspace: need 2*M, prefer M+M*NB) |
* (Workspace: need 2*M, prefer M + M*NB) |
* |
* |
CALL DGELQF( M, N, A, LDA, WORK( ITAU ), |
CALL DGELQF( M, N, A, LDA, WORK( ITAU ), |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
Line 2274
|
Line 2376
|
$ LDU ) |
$ LDU ) |
* |
* |
* Generate Q in A |
* Generate Q in A |
* (Workspace: need 2*M, prefer M+M*NB) |
* (Workspace: need 2*M, prefer M + M*NB) |
* |
* |
CALL DORGLQ( M, N, M, A, LDA, WORK( ITAU ), |
CALL DORGLQ( M, N, M, A, LDA, WORK( ITAU ), |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
Line 2284
|
Line 2386
|
IWORK = ITAUP + M |
IWORK = ITAUP + M |
* |
* |
* Bidiagonalize L in U |
* Bidiagonalize L in U |
* (Workspace: need 4*M, prefer 3*M+2*M*NB) |
* (Workspace: need 4*M, prefer 3*M + 2*M*NB) |
* |
* |
CALL DGEBRD( M, M, U, LDU, S, WORK( IE ), |
CALL DGEBRD( M, M, U, LDU, S, WORK( IE ), |
$ WORK( ITAUQ ), WORK( ITAUP ), |
$ WORK( ITAUQ ), WORK( ITAUP ), |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
* |
* |
* Multiply right vectors bidiagonalizing L by Q in A |
* Multiply right vectors bidiagonalizing L by Q in A |
* (Workspace: need 3*M+N, prefer 3*M+N*NB) |
* (Workspace: need 3*M + N, prefer 3*M + N*NB) |
* |
* |
CALL DORMBR( 'P', 'L', 'T', M, N, M, U, LDU, |
CALL DORMBR( 'P', 'L', 'T', M, N, M, U, LDU, |
$ WORK( ITAUP ), A, LDA, WORK( IWORK ), |
$ WORK( ITAUP ), A, LDA, WORK( IWORK ), |
$ LWORK-IWORK+1, IERR ) |
$ LWORK-IWORK+1, IERR ) |
* |
* |
* Generate left vectors bidiagonalizing L in U |
* Generate left vectors bidiagonalizing L in U |
* (Workspace: need 4*M, prefer 3*M+M*NB) |
* (Workspace: need 4*M, prefer 3*M + M*NB) |
* |
* |
CALL DORGBR( 'Q', M, M, M, U, LDU, WORK( ITAUQ ), |
CALL DORGBR( 'Q', M, M, M, U, LDU, WORK( ITAUQ ), |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
Line 2342
|
Line 2444
|
IWORK = ITAU + M |
IWORK = ITAU + M |
* |
* |
* Compute A=L*Q |
* Compute A=L*Q |
* (Workspace: need M*M+2*M, prefer M*M+M+M*NB) |
* (Workspace: need M*M + 2*M, prefer M*M + M + M*NB) |
* |
* |
CALL DGELQF( M, N, A, LDA, WORK( ITAU ), |
CALL DGELQF( M, N, A, LDA, WORK( ITAU ), |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
Line 2355
|
Line 2457
|
$ WORK( IR+LDWRKR ), LDWRKR ) |
$ WORK( IR+LDWRKR ), LDWRKR ) |
* |
* |
* Generate Q in A |
* Generate Q in A |
* (Workspace: need M*M+2*M, prefer M*M+M+M*NB) |
* (Workspace: need M*M + 2*M, prefer M*M + M + M*NB) |
* |
* |
CALL DORGLQ( M, N, M, A, LDA, WORK( ITAU ), |
CALL DORGLQ( M, N, M, A, LDA, WORK( ITAU ), |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
Line 2365
|
Line 2467
|
IWORK = ITAUP + M |
IWORK = ITAUP + M |
* |
* |
* Bidiagonalize L in WORK(IR) |
* Bidiagonalize L in WORK(IR) |
* (Workspace: need M*M+4*M, prefer M*M+3*M+2*M*NB) |
* (Workspace: need M*M + 4*M, prefer M*M + 3*M + 2*M*NB) |
* |
* |
CALL DGEBRD( M, M, WORK( IR ), LDWRKR, S, |
CALL DGEBRD( M, M, WORK( IR ), LDWRKR, S, |
$ WORK( IE ), WORK( ITAUQ ), |
$ WORK( IE ), WORK( ITAUQ ), |
Line 2374
|
Line 2476
|
* |
* |
* Generate right vectors bidiagonalizing L in |
* Generate right vectors bidiagonalizing L in |
* WORK(IR) |
* WORK(IR) |
* (Workspace: need M*M+4*M, prefer M*M+3*M+(M-1)*NB) |
* (Workspace: need M*M + 4*M, prefer M*M + 3*M + (M-1)*NB) |
* |
* |
CALL DORGBR( 'P', M, M, M, WORK( IR ), LDWRKR, |
CALL DORGBR( 'P', M, M, M, WORK( IR ), LDWRKR, |
$ WORK( ITAUP ), WORK( IWORK ), |
$ WORK( ITAUP ), WORK( IWORK ), |
Line 2383
|
Line 2485
|
* |
* |
* Perform bidiagonal QR iteration, computing right |
* Perform bidiagonal QR iteration, computing right |
* singular vectors of L in WORK(IR) |
* singular vectors of L in WORK(IR) |
* (Workspace: need M*M+BDSPAC) |
* (Workspace: need M*M + BDSPAC) |
* |
* |
CALL DBDSQR( 'U', M, M, 0, 0, S, WORK( IE ), |
CALL DBDSQR( 'U', M, M, 0, 0, S, WORK( IE ), |
$ WORK( IR ), LDWRKR, DUM, 1, DUM, 1, |
$ WORK( IR ), LDWRKR, DUM, 1, DUM, 1, |
Line 2404
|
Line 2506
|
IWORK = ITAU + M |
IWORK = ITAU + M |
* |
* |
* Compute A=L*Q |
* Compute A=L*Q |
* (Workspace: need 2*M, prefer M+M*NB) |
* (Workspace: need 2*M, prefer M + M*NB) |
* |
* |
CALL DGELQF( M, N, A, LDA, WORK( ITAU ), |
CALL DGELQF( M, N, A, LDA, WORK( ITAU ), |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
Line 2414
|
Line 2516
|
CALL DLACPY( 'U', M, N, A, LDA, VT, LDVT ) |
CALL DLACPY( 'U', M, N, A, LDA, VT, LDVT ) |
* |
* |
* Generate Q in VT |
* Generate Q in VT |
* (Workspace: need 2*M, prefer M+M*NB) |
* (Workspace: need 2*M, prefer M + M*NB) |
* |
* |
CALL DORGLQ( M, N, M, VT, LDVT, WORK( ITAU ), |
CALL DORGLQ( M, N, M, VT, LDVT, WORK( ITAU ), |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
Line 2429
|
Line 2531
|
$ LDA ) |
$ LDA ) |
* |
* |
* Bidiagonalize L in A |
* Bidiagonalize L in A |
* (Workspace: need 4*M, prefer 3*M+2*M*NB) |
* (Workspace: need 4*M, prefer 3*M + 2*M*NB) |
* |
* |
CALL DGEBRD( M, M, A, LDA, S, WORK( IE ), |
CALL DGEBRD( M, M, A, LDA, S, WORK( IE ), |
$ WORK( ITAUQ ), WORK( ITAUP ), |
$ WORK( ITAUQ ), WORK( ITAUP ), |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
* |
* |
* Multiply right vectors bidiagonalizing L by Q in VT |
* Multiply right vectors bidiagonalizing L by Q in VT |
* (Workspace: need 3*M+N, prefer 3*M+N*NB) |
* (Workspace: need 3*M + N, prefer 3*M + N*NB) |
* |
* |
CALL DORMBR( 'P', 'L', 'T', M, N, M, A, LDA, |
CALL DORMBR( 'P', 'L', 'T', M, N, M, A, LDA, |
$ WORK( ITAUP ), VT, LDVT, |
$ WORK( ITAUP ), VT, LDVT, |
Line 2471
|
Line 2573
|
LDWRKU = LDA |
LDWRKU = LDA |
IR = IU + LDWRKU*M |
IR = IU + LDWRKU*M |
LDWRKR = LDA |
LDWRKR = LDA |
ELSE IF( LWORK.GE.WRKBL+( LDA+M )*M ) THEN |
ELSE IF( LWORK.GE.WRKBL+( LDA + M )*M ) THEN |
* |
* |
* WORK(IU) is LDA by M and WORK(IR) is M by M |
* WORK(IU) is LDA by M and WORK(IR) is M by M |
* |
* |
Line 2490
|
Line 2592
|
IWORK = ITAU + M |
IWORK = ITAU + M |
* |
* |
* Compute A=L*Q |
* Compute A=L*Q |
* (Workspace: need 2*M*M+2*M, prefer 2*M*M+M+M*NB) |
* (Workspace: need 2*M*M + 2*M, prefer 2*M*M + M + M*NB) |
* |
* |
CALL DGELQF( M, N, A, LDA, WORK( ITAU ), |
CALL DGELQF( M, N, A, LDA, WORK( ITAU ), |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
Line 2503
|
Line 2605
|
$ WORK( IU+LDWRKU ), LDWRKU ) |
$ WORK( IU+LDWRKU ), LDWRKU ) |
* |
* |
* Generate Q in A |
* Generate Q in A |
* (Workspace: need 2*M*M+2*M, prefer 2*M*M+M+M*NB) |
* (Workspace: need 2*M*M + 2*M, prefer 2*M*M + M + M*NB) |
* |
* |
CALL DORGLQ( M, N, M, A, LDA, WORK( ITAU ), |
CALL DORGLQ( M, N, M, A, LDA, WORK( ITAU ), |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
Line 2514
|
Line 2616
|
* |
* |
* Bidiagonalize L in WORK(IU), copying result to |
* Bidiagonalize L in WORK(IU), copying result to |
* WORK(IR) |
* WORK(IR) |
* (Workspace: need 2*M*M+4*M, |
* (Workspace: need 2*M*M + 4*M, |
* prefer 2*M*M+3*M+2*M*NB) |
* prefer 2*M*M+3*M+2*M*NB) |
* |
* |
CALL DGEBRD( M, M, WORK( IU ), LDWRKU, S, |
CALL DGEBRD( M, M, WORK( IU ), LDWRKU, S, |
Line 2525
|
Line 2627
|
$ WORK( IR ), LDWRKR ) |
$ WORK( IR ), LDWRKR ) |
* |
* |
* Generate right bidiagonalizing vectors in WORK(IU) |
* Generate right bidiagonalizing vectors in WORK(IU) |
* (Workspace: need 2*M*M+4*M-1, |
* (Workspace: need 2*M*M + 4*M-1, |
* prefer 2*M*M+3*M+(M-1)*NB) |
* prefer 2*M*M+3*M+(M-1)*NB) |
* |
* |
CALL DORGBR( 'P', M, M, M, WORK( IU ), LDWRKU, |
CALL DORGBR( 'P', M, M, M, WORK( IU ), LDWRKU, |
Line 2533
|
Line 2635
|
$ LWORK-IWORK+1, IERR ) |
$ LWORK-IWORK+1, IERR ) |
* |
* |
* Generate left bidiagonalizing vectors in WORK(IR) |
* Generate left bidiagonalizing vectors in WORK(IR) |
* (Workspace: need 2*M*M+4*M, prefer 2*M*M+3*M+M*NB) |
* (Workspace: need 2*M*M + 4*M, prefer 2*M*M + 3*M + M*NB) |
* |
* |
CALL DORGBR( 'Q', M, M, M, WORK( IR ), LDWRKR, |
CALL DORGBR( 'Q', M, M, M, WORK( IR ), LDWRKR, |
$ WORK( ITAUQ ), WORK( IWORK ), |
$ WORK( ITAUQ ), WORK( IWORK ), |
Line 2543
|
Line 2645
|
* Perform bidiagonal QR iteration, computing left |
* Perform bidiagonal QR iteration, computing left |
* singular vectors of L in WORK(IR) and computing |
* singular vectors of L in WORK(IR) and computing |
* right singular vectors of L in WORK(IU) |
* right singular vectors of L in WORK(IU) |
* (Workspace: need 2*M*M+BDSPAC) |
* (Workspace: need 2*M*M + BDSPAC) |
* |
* |
CALL DBDSQR( 'U', M, M, M, 0, S, WORK( IE ), |
CALL DBDSQR( 'U', M, M, M, 0, S, WORK( IE ), |
$ WORK( IU ), LDWRKU, WORK( IR ), |
$ WORK( IU ), LDWRKU, WORK( IR ), |
Line 2570
|
Line 2672
|
IWORK = ITAU + M |
IWORK = ITAU + M |
* |
* |
* Compute A=L*Q, copying result to VT |
* Compute A=L*Q, copying result to VT |
* (Workspace: need 2*M, prefer M+M*NB) |
* (Workspace: need 2*M, prefer M + M*NB) |
* |
* |
CALL DGELQF( M, N, A, LDA, WORK( ITAU ), |
CALL DGELQF( M, N, A, LDA, WORK( ITAU ), |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
CALL DLACPY( 'U', M, N, A, LDA, VT, LDVT ) |
CALL DLACPY( 'U', M, N, A, LDA, VT, LDVT ) |
* |
* |
* Generate Q in VT |
* Generate Q in VT |
* (Workspace: need 2*M, prefer M+M*NB) |
* (Workspace: need 2*M, prefer M + M*NB) |
* |
* |
CALL DORGLQ( M, N, M, VT, LDVT, WORK( ITAU ), |
CALL DORGLQ( M, N, M, VT, LDVT, WORK( ITAU ), |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
Line 2592
|
Line 2694
|
$ LDA ) |
$ LDA ) |
* |
* |
* Bidiagonalize L in A |
* Bidiagonalize L in A |
* (Workspace: need 4*M, prefer 3*M+2*M*NB) |
* (Workspace: need 4*M, prefer 3*M + 2*M*NB) |
* |
* |
CALL DGEBRD( M, M, A, LDA, S, WORK( IE ), |
CALL DGEBRD( M, M, A, LDA, S, WORK( IE ), |
$ WORK( ITAUQ ), WORK( ITAUP ), |
$ WORK( ITAUQ ), WORK( ITAUP ), |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
* |
* |
* Multiply right vectors bidiagonalizing L by Q in VT |
* Multiply right vectors bidiagonalizing L by Q in VT |
* (Workspace: need 3*M+N, prefer 3*M+N*NB) |
* (Workspace: need 3*M + N, prefer 3*M + N*NB) |
* |
* |
CALL DORMBR( 'P', 'L', 'T', M, N, M, A, LDA, |
CALL DORMBR( 'P', 'L', 'T', M, N, M, A, LDA, |
$ WORK( ITAUP ), VT, LDVT, |
$ WORK( ITAUP ), VT, LDVT, |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
* |
* |
* Generate left bidiagonalizing vectors of L in A |
* Generate left bidiagonalizing vectors of L in A |
* (Workspace: need 4*M, prefer 3*M+M*NB) |
* (Workspace: need 4*M, prefer 3*M + M*NB) |
* |
* |
CALL DORGBR( 'Q', M, M, M, A, LDA, WORK( ITAUQ ), |
CALL DORGBR( 'Q', M, M, M, A, LDA, WORK( ITAUQ ), |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
Line 2650
|
Line 2752
|
IWORK = ITAU + M |
IWORK = ITAU + M |
* |
* |
* Compute A=L*Q |
* Compute A=L*Q |
* (Workspace: need M*M+2*M, prefer M*M+M+M*NB) |
* (Workspace: need M*M + 2*M, prefer M*M + M + M*NB) |
* |
* |
CALL DGELQF( M, N, A, LDA, WORK( ITAU ), |
CALL DGELQF( M, N, A, LDA, WORK( ITAU ), |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
Line 2663
|
Line 2765
|
$ WORK( IU+LDWRKU ), LDWRKU ) |
$ WORK( IU+LDWRKU ), LDWRKU ) |
* |
* |
* Generate Q in A |
* Generate Q in A |
* (Workspace: need M*M+2*M, prefer M*M+M+M*NB) |
* (Workspace: need M*M + 2*M, prefer M*M + M + M*NB) |
* |
* |
CALL DORGLQ( M, N, M, A, LDA, WORK( ITAU ), |
CALL DORGLQ( M, N, M, A, LDA, WORK( ITAU ), |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
Line 2673
|
Line 2775
|
IWORK = ITAUP + M |
IWORK = ITAUP + M |
* |
* |
* Bidiagonalize L in WORK(IU), copying result to U |
* Bidiagonalize L in WORK(IU), copying result to U |
* (Workspace: need M*M+4*M, prefer M*M+3*M+2*M*NB) |
* (Workspace: need M*M + 4*M, prefer M*M + 3*M + 2*M*NB) |
* |
* |
CALL DGEBRD( M, M, WORK( IU ), LDWRKU, S, |
CALL DGEBRD( M, M, WORK( IU ), LDWRKU, S, |
$ WORK( IE ), WORK( ITAUQ ), |
$ WORK( IE ), WORK( ITAUQ ), |
Line 2683
|
Line 2785
|
$ LDU ) |
$ LDU ) |
* |
* |
* Generate right bidiagonalizing vectors in WORK(IU) |
* Generate right bidiagonalizing vectors in WORK(IU) |
* (Workspace: need M*M+4*M-1, |
* (Workspace: need M*M + 4*M-1, |
* prefer M*M+3*M+(M-1)*NB) |
* prefer M*M+3*M+(M-1)*NB) |
* |
* |
CALL DORGBR( 'P', M, M, M, WORK( IU ), LDWRKU, |
CALL DORGBR( 'P', M, M, M, WORK( IU ), LDWRKU, |
Line 2691
|
Line 2793
|
$ LWORK-IWORK+1, IERR ) |
$ LWORK-IWORK+1, IERR ) |
* |
* |
* Generate left bidiagonalizing vectors in U |
* Generate left bidiagonalizing vectors in U |
* (Workspace: need M*M+4*M, prefer M*M+3*M+M*NB) |
* (Workspace: need M*M + 4*M, prefer M*M + 3*M + M*NB) |
* |
* |
CALL DORGBR( 'Q', M, M, M, U, LDU, WORK( ITAUQ ), |
CALL DORGBR( 'Q', M, M, M, U, LDU, WORK( ITAUQ ), |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
Line 2700
|
Line 2802
|
* Perform bidiagonal QR iteration, computing left |
* Perform bidiagonal QR iteration, computing left |
* singular vectors of L in U and computing right |
* singular vectors of L in U and computing right |
* singular vectors of L in WORK(IU) |
* singular vectors of L in WORK(IU) |
* (Workspace: need M*M+BDSPAC) |
* (Workspace: need M*M + BDSPAC) |
* |
* |
CALL DBDSQR( 'U', M, M, M, 0, S, WORK( IE ), |
CALL DBDSQR( 'U', M, M, M, 0, S, WORK( IE ), |
$ WORK( IU ), LDWRKU, U, LDU, DUM, 1, |
$ WORK( IU ), LDWRKU, U, LDU, DUM, 1, |
Line 2721
|
Line 2823
|
IWORK = ITAU + M |
IWORK = ITAU + M |
* |
* |
* Compute A=L*Q, copying result to VT |
* Compute A=L*Q, copying result to VT |
* (Workspace: need 2*M, prefer M+M*NB) |
* (Workspace: need 2*M, prefer M + M*NB) |
* |
* |
CALL DGELQF( M, N, A, LDA, WORK( ITAU ), |
CALL DGELQF( M, N, A, LDA, WORK( ITAU ), |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
CALL DLACPY( 'U', M, N, A, LDA, VT, LDVT ) |
CALL DLACPY( 'U', M, N, A, LDA, VT, LDVT ) |
* |
* |
* Generate Q in VT |
* Generate Q in VT |
* (Workspace: need 2*M, prefer M+M*NB) |
* (Workspace: need 2*M, prefer M + M*NB) |
* |
* |
CALL DORGLQ( M, N, M, VT, LDVT, WORK( ITAU ), |
CALL DORGLQ( M, N, M, VT, LDVT, WORK( ITAU ), |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
Line 2744
|
Line 2846
|
IWORK = ITAUP + M |
IWORK = ITAUP + M |
* |
* |
* Bidiagonalize L in U |
* Bidiagonalize L in U |
* (Workspace: need 4*M, prefer 3*M+2*M*NB) |
* (Workspace: need 4*M, prefer 3*M + 2*M*NB) |
* |
* |
CALL DGEBRD( M, M, U, LDU, S, WORK( IE ), |
CALL DGEBRD( M, M, U, LDU, S, WORK( IE ), |
$ WORK( ITAUQ ), WORK( ITAUP ), |
$ WORK( ITAUQ ), WORK( ITAUP ), |
Line 2752
|
Line 2854
|
* |
* |
* Multiply right bidiagonalizing vectors in U by Q |
* Multiply right bidiagonalizing vectors in U by Q |
* in VT |
* in VT |
* (Workspace: need 3*M+N, prefer 3*M+N*NB) |
* (Workspace: need 3*M + N, prefer 3*M + N*NB) |
* |
* |
CALL DORMBR( 'P', 'L', 'T', M, N, M, U, LDU, |
CALL DORMBR( 'P', 'L', 'T', M, N, M, U, LDU, |
$ WORK( ITAUP ), VT, LDVT, |
$ WORK( ITAUP ), VT, LDVT, |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
* |
* |
* Generate left bidiagonalizing vectors in U |
* Generate left bidiagonalizing vectors in U |
* (Workspace: need 4*M, prefer 3*M+M*NB) |
* (Workspace: need 4*M, prefer 3*M + M*NB) |
* |
* |
CALL DORGBR( 'Q', M, M, M, U, LDU, WORK( ITAUQ ), |
CALL DORGBR( 'Q', M, M, M, U, LDU, WORK( ITAUQ ), |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
Line 2786
|
Line 2888
|
* N right singular vectors to be computed in VT and |
* N right singular vectors to be computed in VT and |
* no left singular vectors to be computed |
* no left singular vectors to be computed |
* |
* |
IF( LWORK.GE.M*M+MAX( N+M, 4*M, BDSPAC ) ) THEN |
IF( LWORK.GE.M*M+MAX( N + M, 4*M, BDSPAC ) ) THEN |
* |
* |
* Sufficient workspace for a fast algorithm |
* Sufficient workspace for a fast algorithm |
* |
* |
Line 2806
|
Line 2908
|
IWORK = ITAU + M |
IWORK = ITAU + M |
* |
* |
* Compute A=L*Q, copying result to VT |
* Compute A=L*Q, copying result to VT |
* (Workspace: need M*M+2*M, prefer M*M+M+M*NB) |
* (Workspace: need M*M + 2*M, prefer M*M + M + M*NB) |
* |
* |
CALL DGELQF( M, N, A, LDA, WORK( ITAU ), |
CALL DGELQF( M, N, A, LDA, WORK( ITAU ), |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
Line 2820
|
Line 2922
|
$ WORK( IR+LDWRKR ), LDWRKR ) |
$ WORK( IR+LDWRKR ), LDWRKR ) |
* |
* |
* Generate Q in VT |
* Generate Q in VT |
* (Workspace: need M*M+M+N, prefer M*M+M+N*NB) |
* (Workspace: need M*M + M + N, prefer M*M + M + N*NB) |
* |
* |
CALL DORGLQ( N, N, M, VT, LDVT, WORK( ITAU ), |
CALL DORGLQ( N, N, M, VT, LDVT, WORK( ITAU ), |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
Line 2830
|
Line 2932
|
IWORK = ITAUP + M |
IWORK = ITAUP + M |
* |
* |
* Bidiagonalize L in WORK(IR) |
* Bidiagonalize L in WORK(IR) |
* (Workspace: need M*M+4*M, prefer M*M+3*M+2*M*NB) |
* (Workspace: need M*M + 4*M, prefer M*M + 3*M + 2*M*NB) |
* |
* |
CALL DGEBRD( M, M, WORK( IR ), LDWRKR, S, |
CALL DGEBRD( M, M, WORK( IR ), LDWRKR, S, |
$ WORK( IE ), WORK( ITAUQ ), |
$ WORK( IE ), WORK( ITAUQ ), |
Line 2838
|
Line 2940
|
$ LWORK-IWORK+1, IERR ) |
$ LWORK-IWORK+1, IERR ) |
* |
* |
* Generate right bidiagonalizing vectors in WORK(IR) |
* Generate right bidiagonalizing vectors in WORK(IR) |
* (Workspace: need M*M+4*M-1, |
* (Workspace: need M*M + 4*M-1, |
* prefer M*M+3*M+(M-1)*NB) |
* prefer M*M+3*M+(M-1)*NB) |
* |
* |
CALL DORGBR( 'P', M, M, M, WORK( IR ), LDWRKR, |
CALL DORGBR( 'P', M, M, M, WORK( IR ), LDWRKR, |
Line 2848
|
Line 2950
|
* |
* |
* Perform bidiagonal QR iteration, computing right |
* Perform bidiagonal QR iteration, computing right |
* singular vectors of L in WORK(IR) |
* singular vectors of L in WORK(IR) |
* (Workspace: need M*M+BDSPAC) |
* (Workspace: need M*M + BDSPAC) |
* |
* |
CALL DBDSQR( 'U', M, M, 0, 0, S, WORK( IE ), |
CALL DBDSQR( 'U', M, M, 0, 0, S, WORK( IE ), |
$ WORK( IR ), LDWRKR, DUM, 1, DUM, 1, |
$ WORK( IR ), LDWRKR, DUM, 1, DUM, 1, |
Line 2873
|
Line 2975
|
IWORK = ITAU + M |
IWORK = ITAU + M |
* |
* |
* Compute A=L*Q, copying result to VT |
* Compute A=L*Q, copying result to VT |
* (Workspace: need 2*M, prefer M+M*NB) |
* (Workspace: need 2*M, prefer M + M*NB) |
* |
* |
CALL DGELQF( M, N, A, LDA, WORK( ITAU ), |
CALL DGELQF( M, N, A, LDA, WORK( ITAU ), |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
CALL DLACPY( 'U', M, N, A, LDA, VT, LDVT ) |
CALL DLACPY( 'U', M, N, A, LDA, VT, LDVT ) |
* |
* |
* Generate Q in VT |
* Generate Q in VT |
* (Workspace: need M+N, prefer M+N*NB) |
* (Workspace: need M + N, prefer M + N*NB) |
* |
* |
CALL DORGLQ( N, N, M, VT, LDVT, WORK( ITAU ), |
CALL DORGLQ( N, N, M, VT, LDVT, WORK( ITAU ), |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
Line 2895
|
Line 2997
|
$ LDA ) |
$ LDA ) |
* |
* |
* Bidiagonalize L in A |
* Bidiagonalize L in A |
* (Workspace: need 4*M, prefer 3*M+2*M*NB) |
* (Workspace: need 4*M, prefer 3*M + 2*M*NB) |
* |
* |
CALL DGEBRD( M, M, A, LDA, S, WORK( IE ), |
CALL DGEBRD( M, M, A, LDA, S, WORK( IE ), |
$ WORK( ITAUQ ), WORK( ITAUP ), |
$ WORK( ITAUQ ), WORK( ITAUP ), |
Line 2903
|
Line 3005
|
* |
* |
* Multiply right bidiagonalizing vectors in A by Q |
* Multiply right bidiagonalizing vectors in A by Q |
* in VT |
* in VT |
* (Workspace: need 3*M+N, prefer 3*M+N*NB) |
* (Workspace: need 3*M + N, prefer 3*M + N*NB) |
* |
* |
CALL DORMBR( 'P', 'L', 'T', M, N, M, A, LDA, |
CALL DORMBR( 'P', 'L', 'T', M, N, M, A, LDA, |
$ WORK( ITAUP ), VT, LDVT, |
$ WORK( ITAUP ), VT, LDVT, |
Line 2926
|
Line 3028
|
* 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 overwritten on A |
* M left singular vectors to be overwritten on A |
* |
* |
IF( LWORK.GE.2*M*M+MAX( N+M, 4*M, BDSPAC ) ) THEN |
IF( LWORK.GE.2*M*M+MAX( N + M, 4*M, BDSPAC ) ) THEN |
* |
* |
* Sufficient workspace for a fast algorithm |
* Sufficient workspace for a fast algorithm |
* |
* |
Line 2938
|
Line 3040
|
LDWRKU = LDA |
LDWRKU = LDA |
IR = IU + LDWRKU*M |
IR = IU + LDWRKU*M |
LDWRKR = LDA |
LDWRKR = LDA |
ELSE IF( LWORK.GE.WRKBL+( LDA+M )*M ) THEN |
ELSE IF( LWORK.GE.WRKBL+( LDA + M )*M ) THEN |
* |
* |
* WORK(IU) is LDA by M and WORK(IR) is M by M |
* WORK(IU) is LDA by M and WORK(IR) is M by M |
* |
* |
Line 2957
|
Line 3059
|
IWORK = ITAU + M |
IWORK = ITAU + M |
* |
* |
* Compute A=L*Q, copying result to VT |
* Compute A=L*Q, copying result to VT |
* (Workspace: need 2*M*M+2*M, prefer 2*M*M+M+M*NB) |
* (Workspace: need 2*M*M + 2*M, prefer 2*M*M + M + M*NB) |
* |
* |
CALL DGELQF( M, N, A, LDA, WORK( ITAU ), |
CALL DGELQF( M, N, A, LDA, WORK( ITAU ), |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
CALL DLACPY( 'U', M, N, A, LDA, VT, LDVT ) |
CALL DLACPY( 'U', M, N, A, LDA, VT, LDVT ) |
* |
* |
* Generate Q in VT |
* Generate Q in VT |
* (Workspace: need 2*M*M+M+N, prefer 2*M*M+M+N*NB) |
* (Workspace: need 2*M*M + M + N, prefer 2*M*M + M + N*NB) |
* |
* |
CALL DORGLQ( N, N, M, VT, LDVT, WORK( ITAU ), |
CALL DORGLQ( N, N, M, VT, LDVT, WORK( ITAU ), |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
Line 2982
|
Line 3084
|
* |
* |
* Bidiagonalize L in WORK(IU), copying result to |
* Bidiagonalize L in WORK(IU), copying result to |
* WORK(IR) |
* WORK(IR) |
* (Workspace: need 2*M*M+4*M, |
* (Workspace: need 2*M*M + 4*M, |
* prefer 2*M*M+3*M+2*M*NB) |
* prefer 2*M*M+3*M+2*M*NB) |
* |
* |
CALL DGEBRD( M, M, WORK( IU ), LDWRKU, S, |
CALL DGEBRD( M, M, WORK( IU ), LDWRKU, S, |
Line 2993
|
Line 3095
|
$ WORK( IR ), LDWRKR ) |
$ WORK( IR ), LDWRKR ) |
* |
* |
* Generate right bidiagonalizing vectors in WORK(IU) |
* Generate right bidiagonalizing vectors in WORK(IU) |
* (Workspace: need 2*M*M+4*M-1, |
* (Workspace: need 2*M*M + 4*M-1, |
* prefer 2*M*M+3*M+(M-1)*NB) |
* prefer 2*M*M+3*M+(M-1)*NB) |
* |
* |
CALL DORGBR( 'P', M, M, M, WORK( IU ), LDWRKU, |
CALL DORGBR( 'P', M, M, M, WORK( IU ), LDWRKU, |
Line 3001
|
Line 3103
|
$ LWORK-IWORK+1, IERR ) |
$ LWORK-IWORK+1, IERR ) |
* |
* |
* Generate left bidiagonalizing vectors in WORK(IR) |
* Generate left bidiagonalizing vectors in WORK(IR) |
* (Workspace: need 2*M*M+4*M, prefer 2*M*M+3*M+M*NB) |
* (Workspace: need 2*M*M + 4*M, prefer 2*M*M + 3*M + M*NB) |
* |
* |
CALL DORGBR( 'Q', M, M, M, WORK( IR ), LDWRKR, |
CALL DORGBR( 'Q', M, M, M, WORK( IR ), LDWRKR, |
$ WORK( ITAUQ ), WORK( IWORK ), |
$ WORK( ITAUQ ), WORK( IWORK ), |
Line 3011
|
Line 3113
|
* Perform bidiagonal QR iteration, computing left |
* Perform bidiagonal QR iteration, computing left |
* singular vectors of L in WORK(IR) and computing |
* singular vectors of L in WORK(IR) and computing |
* right singular vectors of L in WORK(IU) |
* right singular vectors of L in WORK(IU) |
* (Workspace: need 2*M*M+BDSPAC) |
* (Workspace: need 2*M*M + BDSPAC) |
* |
* |
CALL DBDSQR( 'U', M, M, M, 0, S, WORK( IE ), |
CALL DBDSQR( 'U', M, M, M, 0, S, WORK( IE ), |
$ WORK( IU ), LDWRKU, WORK( IR ), |
$ WORK( IU ), LDWRKU, WORK( IR ), |
Line 3041
|
Line 3143
|
IWORK = ITAU + M |
IWORK = ITAU + M |
* |
* |
* Compute A=L*Q, copying result to VT |
* Compute A=L*Q, copying result to VT |
* (Workspace: need 2*M, prefer M+M*NB) |
* (Workspace: need 2*M, prefer M + M*NB) |
* |
* |
CALL DGELQF( M, N, A, LDA, WORK( ITAU ), |
CALL DGELQF( M, N, A, LDA, WORK( ITAU ), |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
CALL DLACPY( 'U', M, N, A, LDA, VT, LDVT ) |
CALL DLACPY( 'U', M, N, A, LDA, VT, LDVT ) |
* |
* |
* Generate Q in VT |
* Generate Q in VT |
* (Workspace: need M+N, prefer M+N*NB) |
* (Workspace: need M + N, prefer M + N*NB) |
* |
* |
CALL DORGLQ( N, N, M, VT, LDVT, WORK( ITAU ), |
CALL DORGLQ( N, N, M, VT, LDVT, WORK( ITAU ), |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
Line 3063
|
Line 3165
|
$ LDA ) |
$ LDA ) |
* |
* |
* Bidiagonalize L in A |
* Bidiagonalize L in A |
* (Workspace: need 4*M, prefer 3*M+2*M*NB) |
* (Workspace: need 4*M, prefer 3*M + 2*M*NB) |
* |
* |
CALL DGEBRD( M, M, A, LDA, S, WORK( IE ), |
CALL DGEBRD( M, M, A, LDA, S, WORK( IE ), |
$ WORK( ITAUQ ), WORK( ITAUP ), |
$ WORK( ITAUQ ), WORK( ITAUP ), |
Line 3071
|
Line 3173
|
* |
* |
* Multiply right bidiagonalizing vectors in A by Q |
* Multiply right bidiagonalizing vectors in A by Q |
* in VT |
* in VT |
* (Workspace: need 3*M+N, prefer 3*M+N*NB) |
* (Workspace: need 3*M + N, prefer 3*M + N*NB) |
* |
* |
CALL DORMBR( 'P', 'L', 'T', M, N, M, A, LDA, |
CALL DORMBR( 'P', 'L', 'T', M, N, M, A, LDA, |
$ WORK( ITAUP ), VT, LDVT, |
$ WORK( ITAUP ), VT, LDVT, |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
* |
* |
* Generate left bidiagonalizing vectors in A |
* Generate left bidiagonalizing vectors in A |
* (Workspace: need 4*M, prefer 3*M+M*NB) |
* (Workspace: need 4*M, prefer 3*M + M*NB) |
* |
* |
CALL DORGBR( 'Q', M, M, M, A, LDA, WORK( ITAUQ ), |
CALL DORGBR( 'Q', M, M, M, A, LDA, WORK( ITAUQ ), |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
Line 3102
|
Line 3204
|
* 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 |
* |
* |
IF( LWORK.GE.M*M+MAX( N+M, 4*M, BDSPAC ) ) THEN |
IF( LWORK.GE.M*M+MAX( N + M, 4*M, BDSPAC ) ) THEN |
* |
* |
* Sufficient workspace for a fast algorithm |
* Sufficient workspace for a fast algorithm |
* |
* |
Line 3122
|
Line 3224
|
IWORK = ITAU + M |
IWORK = ITAU + M |
* |
* |
* Compute A=L*Q, copying result to VT |
* Compute A=L*Q, copying result to VT |
* (Workspace: need M*M+2*M, prefer M*M+M+M*NB) |
* (Workspace: need M*M + 2*M, prefer M*M + M + M*NB) |
* |
* |
CALL DGELQF( M, N, A, LDA, WORK( ITAU ), |
CALL DGELQF( M, N, A, LDA, WORK( ITAU ), |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
CALL DLACPY( 'U', M, N, A, LDA, VT, LDVT ) |
CALL DLACPY( 'U', M, N, A, LDA, VT, LDVT ) |
* |
* |
* Generate Q in VT |
* Generate Q in VT |
* (Workspace: need M*M+M+N, prefer M*M+M+N*NB) |
* (Workspace: need M*M + M + N, prefer M*M + M + N*NB) |
* |
* |
CALL DORGLQ( N, N, M, VT, LDVT, WORK( ITAU ), |
CALL DORGLQ( N, N, M, VT, LDVT, WORK( ITAU ), |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
Line 3146
|
Line 3248
|
IWORK = ITAUP + M |
IWORK = ITAUP + M |
* |
* |
* Bidiagonalize L in WORK(IU), copying result to U |
* Bidiagonalize L in WORK(IU), copying result to U |
* (Workspace: need M*M+4*M, prefer M*M+3*M+2*M*NB) |
* (Workspace: need M*M + 4*M, prefer M*M + 3*M + 2*M*NB) |
* |
* |
CALL DGEBRD( M, M, WORK( IU ), LDWRKU, S, |
CALL DGEBRD( M, M, WORK( IU ), LDWRKU, S, |
$ WORK( IE ), WORK( ITAUQ ), |
$ WORK( IE ), WORK( ITAUQ ), |
Line 3156
|
Line 3258
|
$ LDU ) |
$ LDU ) |
* |
* |
* Generate right bidiagonalizing vectors in WORK(IU) |
* Generate right bidiagonalizing vectors in WORK(IU) |
* (Workspace: need M*M+4*M, prefer M*M+3*M+(M-1)*NB) |
* (Workspace: need M*M + 4*M, prefer M*M + 3*M + (M-1)*NB) |
* |
* |
CALL DORGBR( 'P', M, M, M, WORK( IU ), LDWRKU, |
CALL DORGBR( 'P', M, M, M, WORK( IU ), LDWRKU, |
$ WORK( ITAUP ), WORK( IWORK ), |
$ WORK( ITAUP ), WORK( IWORK ), |
$ LWORK-IWORK+1, IERR ) |
$ LWORK-IWORK+1, IERR ) |
* |
* |
* Generate left bidiagonalizing vectors in U |
* Generate left bidiagonalizing vectors in U |
* (Workspace: need M*M+4*M, prefer M*M+3*M+M*NB) |
* (Workspace: need M*M + 4*M, prefer M*M + 3*M + M*NB) |
* |
* |
CALL DORGBR( 'Q', M, M, M, U, LDU, WORK( ITAUQ ), |
CALL DORGBR( 'Q', M, M, M, U, LDU, WORK( ITAUQ ), |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
Line 3172
|
Line 3274
|
* Perform bidiagonal QR iteration, computing left |
* Perform bidiagonal QR iteration, computing left |
* singular vectors of L in U and computing right |
* singular vectors of L in U and computing right |
* singular vectors of L in WORK(IU) |
* singular vectors of L in WORK(IU) |
* (Workspace: need M*M+BDSPAC) |
* (Workspace: need M*M + BDSPAC) |
* |
* |
CALL DBDSQR( 'U', M, M, M, 0, S, WORK( IE ), |
CALL DBDSQR( 'U', M, M, M, 0, S, WORK( IE ), |
$ WORK( IU ), LDWRKU, U, LDU, DUM, 1, |
$ WORK( IU ), LDWRKU, U, LDU, DUM, 1, |
Line 3197
|
Line 3299
|
IWORK = ITAU + M |
IWORK = ITAU + M |
* |
* |
* Compute A=L*Q, copying result to VT |
* Compute A=L*Q, copying result to VT |
* (Workspace: need 2*M, prefer M+M*NB) |
* (Workspace: need 2*M, prefer M + M*NB) |
* |
* |
CALL DGELQF( M, N, A, LDA, WORK( ITAU ), |
CALL DGELQF( M, N, A, LDA, WORK( ITAU ), |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
CALL DLACPY( 'U', M, N, A, LDA, VT, LDVT ) |
CALL DLACPY( 'U', M, N, A, LDA, VT, LDVT ) |
* |
* |
* Generate Q in VT |
* Generate Q in VT |
* (Workspace: need M+N, prefer M+N*NB) |
* (Workspace: need M + N, prefer M + N*NB) |
* |
* |
CALL DORGLQ( N, N, M, VT, LDVT, WORK( ITAU ), |
CALL DORGLQ( N, N, M, VT, LDVT, WORK( ITAU ), |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
Line 3220
|
Line 3322
|
IWORK = ITAUP + M |
IWORK = ITAUP + M |
* |
* |
* Bidiagonalize L in U |
* Bidiagonalize L in U |
* (Workspace: need 4*M, prefer 3*M+2*M*NB) |
* (Workspace: need 4*M, prefer 3*M + 2*M*NB) |
* |
* |
CALL DGEBRD( M, M, U, LDU, S, WORK( IE ), |
CALL DGEBRD( M, M, U, LDU, S, WORK( IE ), |
$ WORK( ITAUQ ), WORK( ITAUP ), |
$ WORK( ITAUQ ), WORK( ITAUP ), |
Line 3228
|
Line 3330
|
* |
* |
* Multiply right bidiagonalizing vectors in U by Q |
* Multiply right bidiagonalizing vectors in U by Q |
* in VT |
* in VT |
* (Workspace: need 3*M+N, prefer 3*M+N*NB) |
* (Workspace: need 3*M + N, prefer 3*M + N*NB) |
* |
* |
CALL DORMBR( 'P', 'L', 'T', M, N, M, U, LDU, |
CALL DORMBR( 'P', 'L', 'T', M, N, M, U, LDU, |
$ WORK( ITAUP ), VT, LDVT, |
$ WORK( ITAUP ), VT, LDVT, |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
* |
* |
* Generate left bidiagonalizing vectors in U |
* Generate left bidiagonalizing vectors in U |
* (Workspace: need 4*M, prefer 3*M+M*NB) |
* (Workspace: need 4*M, prefer 3*M + M*NB) |
* |
* |
CALL DORGBR( 'Q', M, M, M, U, LDU, WORK( ITAUQ ), |
CALL DORGBR( 'Q', M, M, M, U, LDU, WORK( ITAUQ ), |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
Line 3269
|
Line 3371
|
IWORK = ITAUP + M |
IWORK = ITAUP + M |
* |
* |
* Bidiagonalize A |
* Bidiagonalize A |
* (Workspace: need 3*M+N, prefer 3*M+(M+N)*NB) |
* (Workspace: need 3*M + N, prefer 3*M + (M + N)*NB) |
* |
* |
CALL DGEBRD( M, N, A, LDA, S, WORK( IE ), WORK( ITAUQ ), |
CALL DGEBRD( M, N, A, LDA, S, WORK( IE ), WORK( ITAUQ ), |
$ WORK( ITAUP ), WORK( IWORK ), LWORK-IWORK+1, |
$ WORK( ITAUP ), WORK( IWORK ), LWORK-IWORK+1, |
Line 3278
|
Line 3380
|
* |
* |
* If left singular vectors desired in U, copy result to U |
* If left singular vectors desired in U, copy result to U |
* and generate left bidiagonalizing vectors in U |
* and generate left bidiagonalizing vectors in U |
* (Workspace: need 4*M-1, prefer 3*M+(M-1)*NB) |
* (Workspace: need 4*M-1, prefer 3*M + (M-1)*NB) |
* |
* |
CALL DLACPY( 'L', M, M, A, LDA, U, LDU ) |
CALL DLACPY( 'L', M, M, A, LDA, U, LDU ) |
CALL DORGBR( 'Q', M, M, N, U, LDU, WORK( ITAUQ ), |
CALL DORGBR( 'Q', M, M, N, U, LDU, WORK( ITAUQ ), |
Line 3288
|
Line 3390
|
* |
* |
* If right singular vectors desired in VT, copy result to |
* If right singular vectors desired in VT, copy result to |
* VT and generate right bidiagonalizing vectors in VT |
* VT and generate right bidiagonalizing vectors in VT |
* (Workspace: need 3*M+NRVT, prefer 3*M+NRVT*NB) |
* (Workspace: need 3*M + NRVT, prefer 3*M + NRVT*NB) |
* |
* |
CALL DLACPY( 'U', M, N, A, LDA, VT, LDVT ) |
CALL DLACPY( 'U', M, N, A, LDA, VT, LDVT ) |
IF( WNTVA ) |
IF( WNTVA ) |
Line 3302
|
Line 3404
|
* |
* |
* If left singular vectors desired in A, generate left |
* If left singular vectors desired in A, generate left |
* bidiagonalizing vectors in A |
* bidiagonalizing vectors in A |
* (Workspace: need 4*M-1, prefer 3*M+(M-1)*NB) |
* (Workspace: need 4*M-1, prefer 3*M + (M-1)*NB) |
* |
* |
CALL DORGBR( 'Q', M, M, N, A, LDA, WORK( ITAUQ ), |
CALL DORGBR( 'Q', M, M, N, A, LDA, WORK( ITAUQ ), |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
Line 3311
|
Line 3413
|
* |
* |
* If right singular vectors desired in A, generate right |
* If right singular vectors desired in A, generate right |
* bidiagonalizing vectors in A |
* bidiagonalizing vectors in A |
* (Workspace: need 4*M, prefer 3*M+M*NB) |
* (Workspace: need 4*M, prefer 3*M + M*NB) |
* |
* |
CALL DORGBR( 'P', M, N, M, A, LDA, WORK( ITAUP ), |
CALL DORGBR( 'P', M, N, M, A, LDA, WORK( ITAUP ), |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |
$ WORK( IWORK ), LWORK-IWORK+1, IERR ) |