version 1.5, 2010/08/07 13:22:31
|
version 1.14, 2016/08/27 15:34:45
|
Line 1
|
Line 1
|
|
*> \brief <b> ZGELSS solves overdetermined or underdetermined systems for GE matrices</b> |
|
* |
|
* =========== DOCUMENTATION =========== |
|
* |
|
* Online html documentation available at |
|
* http://www.netlib.org/lapack/explore-html/ |
|
* |
|
*> \htmlonly |
|
*> Download ZGELSS + dependencies |
|
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zgelss.f"> |
|
*> [TGZ]</a> |
|
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zgelss.f"> |
|
*> [ZIP]</a> |
|
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zgelss.f"> |
|
*> [TXT]</a> |
|
*> \endhtmlonly |
|
* |
|
* Definition: |
|
* =========== |
|
* |
|
* SUBROUTINE ZGELSS( M, N, NRHS, A, LDA, B, LDB, S, RCOND, RANK, |
|
* WORK, LWORK, RWORK, INFO ) |
|
* |
|
* .. Scalar Arguments .. |
|
* INTEGER INFO, LDA, LDB, LWORK, M, N, NRHS, RANK |
|
* DOUBLE PRECISION RCOND |
|
* .. |
|
* .. Array Arguments .. |
|
* DOUBLE PRECISION RWORK( * ), S( * ) |
|
* COMPLEX*16 A( LDA, * ), B( LDB, * ), WORK( * ) |
|
* .. |
|
* |
|
* |
|
*> \par Purpose: |
|
* ============= |
|
*> |
|
*> \verbatim |
|
*> |
|
*> ZGELSS computes the minimum norm solution to a complex linear |
|
*> least squares problem: |
|
*> |
|
*> Minimize 2-norm(| b - A*x |). |
|
*> |
|
*> using the singular value decomposition (SVD) of A. A is an M-by-N |
|
*> matrix which may be rank-deficient. |
|
*> |
|
*> Several right hand side vectors b and solution vectors x can be |
|
*> handled in a single call; they are stored as the columns of the |
|
*> M-by-NRHS right hand side matrix B and the N-by-NRHS solution matrix |
|
*> X. |
|
*> |
|
*> The effective rank of A is determined by treating as zero those |
|
*> singular values which are less than RCOND times the largest singular |
|
*> value. |
|
*> \endverbatim |
|
* |
|
* Arguments: |
|
* ========== |
|
* |
|
*> \param[in] M |
|
*> \verbatim |
|
*> M is INTEGER |
|
*> The number of rows of the matrix A. M >= 0. |
|
*> \endverbatim |
|
*> |
|
*> \param[in] N |
|
*> \verbatim |
|
*> N is INTEGER |
|
*> The number of columns of the matrix A. N >= 0. |
|
*> \endverbatim |
|
*> |
|
*> \param[in] NRHS |
|
*> \verbatim |
|
*> NRHS is INTEGER |
|
*> The number of right hand sides, i.e., the number of columns |
|
*> of the matrices B and X. NRHS >= 0. |
|
*> \endverbatim |
|
*> |
|
*> \param[in,out] A |
|
*> \verbatim |
|
*> A is COMPLEX*16 array, dimension (LDA,N) |
|
*> On entry, the M-by-N matrix A. |
|
*> On exit, the first min(m,n) rows of A are overwritten with |
|
*> its right singular vectors, stored rowwise. |
|
*> \endverbatim |
|
*> |
|
*> \param[in] LDA |
|
*> \verbatim |
|
*> LDA is INTEGER |
|
*> The leading dimension of the array A. LDA >= max(1,M). |
|
*> \endverbatim |
|
*> |
|
*> \param[in,out] B |
|
*> \verbatim |
|
*> B is COMPLEX*16 array, dimension (LDB,NRHS) |
|
*> On entry, the M-by-NRHS right hand side matrix B. |
|
*> On exit, B is overwritten by the N-by-NRHS solution matrix X. |
|
*> If m >= n and RANK = n, the residual sum-of-squares for |
|
*> the solution in the i-th column is given by the sum of |
|
*> squares of the modulus of elements n+1:m in that column. |
|
*> \endverbatim |
|
*> |
|
*> \param[in] LDB |
|
*> \verbatim |
|
*> LDB is INTEGER |
|
*> The leading dimension of the array B. LDB >= max(1,M,N). |
|
*> \endverbatim |
|
*> |
|
*> \param[out] S |
|
*> \verbatim |
|
*> S is DOUBLE PRECISION array, dimension (min(M,N)) |
|
*> The singular values of A in decreasing order. |
|
*> The condition number of A in the 2-norm = S(1)/S(min(m,n)). |
|
*> \endverbatim |
|
*> |
|
*> \param[in] RCOND |
|
*> \verbatim |
|
*> RCOND is DOUBLE PRECISION |
|
*> RCOND is used to determine the effective rank of A. |
|
*> Singular values S(i) <= RCOND*S(1) are treated as zero. |
|
*> If RCOND < 0, machine precision is used instead. |
|
*> \endverbatim |
|
*> |
|
*> \param[out] RANK |
|
*> \verbatim |
|
*> RANK is INTEGER |
|
*> The effective rank of A, i.e., the number of singular values |
|
*> which are greater than RCOND*S(1). |
|
*> \endverbatim |
|
*> |
|
*> \param[out] WORK |
|
*> \verbatim |
|
*> WORK is COMPLEX*16 array, dimension (MAX(1,LWORK)) |
|
*> On exit, if INFO = 0, WORK(1) returns the optimal LWORK. |
|
*> \endverbatim |
|
*> |
|
*> \param[in] LWORK |
|
*> \verbatim |
|
*> LWORK is INTEGER |
|
*> The dimension of the array WORK. LWORK >= 1, and also: |
|
*> LWORK >= 2*min(M,N) + max(M,N,NRHS) |
|
*> 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] RWORK |
|
*> \verbatim |
|
*> RWORK is DOUBLE PRECISION array, dimension (5*min(M,N)) |
|
*> \endverbatim |
|
*> |
|
*> \param[out] INFO |
|
*> \verbatim |
|
*> INFO is INTEGER |
|
*> = 0: successful exit |
|
*> < 0: if INFO = -i, the i-th argument had an illegal value. |
|
*> > 0: the algorithm for computing the SVD failed to converge; |
|
*> if INFO = i, i off-diagonal elements of an intermediate |
|
*> bidiagonal form did not converge to zero. |
|
*> \endverbatim |
|
* |
|
* Authors: |
|
* ======== |
|
* |
|
*> \author Univ. of Tennessee |
|
*> \author Univ. of California Berkeley |
|
*> \author Univ. of Colorado Denver |
|
*> \author NAG Ltd. |
|
* |
|
*> \date June 2016 |
|
* |
|
*> \ingroup complex16GEsolve |
|
* |
|
* ===================================================================== |
SUBROUTINE ZGELSS( M, N, NRHS, A, LDA, B, LDB, S, RCOND, RANK, |
SUBROUTINE ZGELSS( M, N, NRHS, A, LDA, B, LDB, S, RCOND, RANK, |
$ WORK, LWORK, RWORK, INFO ) |
$ WORK, LWORK, RWORK, INFO ) |
* |
* |
* -- LAPACK driver routine (version 3.2) -- |
* -- 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 |
* June 2016 |
* |
* |
* .. Scalar Arguments .. |
* .. Scalar Arguments .. |
INTEGER INFO, LDA, LDB, LWORK, M, N, NRHS, RANK |
INTEGER INFO, LDA, LDB, LWORK, M, N, NRHS, RANK |
Line 15
|
Line 192
|
COMPLEX*16 A( LDA, * ), B( LDB, * ), WORK( * ) |
COMPLEX*16 A( LDA, * ), B( LDB, * ), WORK( * ) |
* .. |
* .. |
* |
* |
* Purpose |
|
* ======= |
|
* |
|
* ZGELSS computes the minimum norm solution to a complex linear |
|
* least squares problem: |
|
* |
|
* Minimize 2-norm(| b - A*x |). |
|
* |
|
* using the singular value decomposition (SVD) of A. A is an M-by-N |
|
* matrix which may be rank-deficient. |
|
* |
|
* Several right hand side vectors b and solution vectors x can be |
|
* handled in a single call; they are stored as the columns of the |
|
* M-by-NRHS right hand side matrix B and the N-by-NRHS solution matrix |
|
* X. |
|
* |
|
* The effective rank of A is determined by treating as zero those |
|
* singular values which are less than RCOND times the largest singular |
|
* value. |
|
* |
|
* Arguments |
|
* ========= |
|
* |
|
* M (input) INTEGER |
|
* The number of rows of the matrix A. M >= 0. |
|
* |
|
* N (input) INTEGER |
|
* The number of columns of the matrix A. N >= 0. |
|
* |
|
* NRHS (input) INTEGER |
|
* The number of right hand sides, i.e., the number of columns |
|
* of the matrices B and X. NRHS >= 0. |
|
* |
|
* A (input/output) COMPLEX*16 array, dimension (LDA,N) |
|
* On entry, the M-by-N matrix A. |
|
* On exit, the first min(m,n) rows of A are overwritten with |
|
* its right singular vectors, stored rowwise. |
|
* |
|
* LDA (input) INTEGER |
|
* The leading dimension of the array A. LDA >= max(1,M). |
|
* |
|
* B (input/output) COMPLEX*16 array, dimension (LDB,NRHS) |
|
* On entry, the M-by-NRHS right hand side matrix B. |
|
* On exit, B is overwritten by the N-by-NRHS solution matrix X. |
|
* If m >= n and RANK = n, the residual sum-of-squares for |
|
* the solution in the i-th column is given by the sum of |
|
* squares of the modulus of elements n+1:m in that column. |
|
* |
|
* LDB (input) INTEGER |
|
* The leading dimension of the array B. LDB >= max(1,M,N). |
|
* |
|
* S (output) DOUBLE PRECISION array, dimension (min(M,N)) |
|
* The singular values of A in decreasing order. |
|
* The condition number of A in the 2-norm = S(1)/S(min(m,n)). |
|
* |
|
* RCOND (input) DOUBLE PRECISION |
|
* RCOND is used to determine the effective rank of A. |
|
* Singular values S(i) <= RCOND*S(1) are treated as zero. |
|
* If RCOND < 0, machine precision is used instead. |
|
* |
|
* RANK (output) INTEGER |
|
* The effective rank of A, i.e., the number of singular values |
|
* which are greater than RCOND*S(1). |
|
* |
|
* WORK (workspace/output) COMPLEX*16 array, dimension (MAX(1,LWORK)) |
|
* On exit, if INFO = 0, WORK(1) returns the optimal LWORK. |
|
* |
|
* LWORK (input) INTEGER |
|
* The dimension of the array WORK. LWORK >= 1, and also: |
|
* LWORK >= 2*min(M,N) + max(M,N,NRHS) |
|
* 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. |
|
* |
|
* RWORK (workspace) DOUBLE PRECISION array, dimension (5*min(M,N)) |
|
* |
|
* INFO (output) INTEGER |
|
* = 0: successful exit |
|
* < 0: if INFO = -i, the i-th argument had an illegal value. |
|
* > 0: the algorithm for computing the SVD failed to converge; |
|
* if INFO = i, i off-diagonal elements of an intermediate |
|
* bidiagonal form did not converge to zero. |
|
* |
|
* ===================================================================== |
* ===================================================================== |
* |
* |
* .. Parameters .. |
* .. Parameters .. |
Line 115
|
Line 206
|
INTEGER BL, CHUNK, I, IASCL, IBSCL, IE, IL, IRWORK, |
INTEGER BL, CHUNK, I, IASCL, IBSCL, IE, IL, IRWORK, |
$ ITAU, ITAUP, ITAUQ, IWORK, LDWORK, MAXMN, |
$ ITAU, ITAUP, ITAUQ, IWORK, LDWORK, MAXMN, |
$ MAXWRK, MINMN, MINWRK, MM, MNTHR |
$ MAXWRK, MINMN, MINWRK, MM, MNTHR |
|
INTEGER LWORK_ZGEQRF, LWORK_ZUNMQR, LWORK_ZGEBRD, |
|
$ LWORK_ZUNMBR, LWORK_ZUNGBR, LWORK_ZUNMLQ, |
|
$ LWORK_ZGELQF |
DOUBLE PRECISION ANRM, BIGNUM, BNRM, EPS, SFMIN, SMLNUM, THR |
DOUBLE PRECISION ANRM, BIGNUM, BNRM, EPS, SFMIN, SMLNUM, THR |
* .. |
* .. |
* .. Local Arrays .. |
* .. Local Arrays .. |
COMPLEX*16 VDUM( 1 ) |
COMPLEX*16 DUM( 1 ) |
* .. |
* .. |
* .. External Subroutines .. |
* .. External Subroutines .. |
EXTERNAL DLABAD, DLASCL, DLASET, XERBLA, ZBDSQR, ZCOPY, |
EXTERNAL DLABAD, DLASCL, DLASET, XERBLA, ZBDSQR, ZCOPY, |
Line 173
|
Line 267
|
* Path 1a - overdetermined, with many more rows than |
* Path 1a - overdetermined, with many more rows than |
* columns |
* columns |
* |
* |
|
* Compute space needed for ZGEQRF |
|
CALL ZGEQRF( M, N, A, LDA, DUM(1), DUM(1), -1, INFO ) |
|
LWORK_ZGEQRF=DUM(1) |
|
* Compute space needed for ZUNMQR |
|
CALL ZUNMQR( 'L', 'C', M, NRHS, N, A, LDA, DUM(1), B, |
|
$ LDB, DUM(1), -1, INFO ) |
|
LWORK_ZUNMQR=DUM(1) |
MM = N |
MM = N |
MAXWRK = MAX( MAXWRK, N + N*ILAENV( 1, 'ZGEQRF', ' ', M, |
MAXWRK = MAX( MAXWRK, N + N*ILAENV( 1, 'ZGEQRF', ' ', M, |
$ N, -1, -1 ) ) |
$ N, -1, -1 ) ) |
Line 183
|
Line 284
|
* |
* |
* Path 1 - overdetermined or exactly determined |
* Path 1 - overdetermined or exactly determined |
* |
* |
MAXWRK = MAX( MAXWRK, 2*N + ( MM + N )*ILAENV( 1, |
* Compute space needed for ZGEBRD |
$ 'ZGEBRD', ' ', MM, N, -1, -1 ) ) |
CALL ZGEBRD( MM, N, A, LDA, S, S, DUM(1), DUM(1), DUM(1), |
MAXWRK = MAX( MAXWRK, 2*N + NRHS*ILAENV( 1, 'ZUNMBR', |
$ -1, INFO ) |
$ 'QLC', MM, NRHS, N, -1 ) ) |
LWORK_ZGEBRD=DUM(1) |
MAXWRK = MAX( MAXWRK, 2*N + ( N - 1 )*ILAENV( 1, |
* Compute space needed for ZUNMBR |
$ 'ZUNGBR', 'P', N, N, N, -1 ) ) |
CALL ZUNMBR( 'Q', 'L', 'C', MM, NRHS, N, A, LDA, DUM(1), |
|
$ B, LDB, DUM(1), -1, INFO ) |
|
LWORK_ZUNMBR=DUM(1) |
|
* Compute space needed for ZUNGBR |
|
CALL ZUNGBR( 'P', N, N, N, A, LDA, DUM(1), |
|
$ DUM(1), -1, INFO ) |
|
LWORK_ZUNGBR=DUM(1) |
|
* Compute total workspace needed |
|
MAXWRK = MAX( MAXWRK, 2*N + LWORK_ZGEBRD ) |
|
MAXWRK = MAX( MAXWRK, 2*N + LWORK_ZUNMBR ) |
|
MAXWRK = MAX( MAXWRK, 2*N + LWORK_ZUNGBR ) |
MAXWRK = MAX( MAXWRK, N*NRHS ) |
MAXWRK = MAX( MAXWRK, N*NRHS ) |
MINWRK = 2*N + MAX( NRHS, M ) |
MINWRK = 2*N + MAX( NRHS, M ) |
END IF |
END IF |
Line 199
|
Line 310
|
* Path 2a - underdetermined, with many more columns |
* Path 2a - underdetermined, with many more columns |
* than rows |
* than rows |
* |
* |
MAXWRK = M + M*ILAENV( 1, 'ZGELQF', ' ', M, N, -1, |
* Compute space needed for ZGELQF |
$ -1 ) |
CALL ZGELQF( M, N, A, LDA, DUM(1), DUM(1), |
MAXWRK = MAX( MAXWRK, 3*M + M*M + 2*M*ILAENV( 1, |
$ -1, INFO ) |
$ 'ZGEBRD', ' ', M, M, -1, -1 ) ) |
LWORK_ZGELQF=DUM(1) |
MAXWRK = MAX( MAXWRK, 3*M + M*M + NRHS*ILAENV( 1, |
* Compute space needed for ZGEBRD |
$ 'ZUNMBR', 'QLC', M, NRHS, M, -1 ) ) |
CALL ZGEBRD( M, M, A, LDA, S, S, DUM(1), DUM(1), |
MAXWRK = MAX( MAXWRK, 3*M + M*M + ( M - 1 )*ILAENV( 1, |
$ DUM(1), -1, INFO ) |
$ 'ZUNGBR', 'P', M, M, M, -1 ) ) |
LWORK_ZGEBRD=DUM(1) |
|
* Compute space needed for ZUNMBR |
|
CALL ZUNMBR( 'Q', 'L', 'C', M, NRHS, N, A, LDA, |
|
$ DUM(1), B, LDB, DUM(1), -1, INFO ) |
|
LWORK_ZUNMBR=DUM(1) |
|
* Compute space needed for ZUNGBR |
|
CALL ZUNGBR( 'P', M, M, M, A, LDA, DUM(1), |
|
$ DUM(1), -1, INFO ) |
|
LWORK_ZUNGBR=DUM(1) |
|
* Compute space needed for ZUNMLQ |
|
CALL ZUNMLQ( 'L', 'C', N, NRHS, M, A, LDA, DUM(1), |
|
$ B, LDB, DUM(1), -1, INFO ) |
|
LWORK_ZUNMLQ=DUM(1) |
|
* Compute total workspace needed |
|
MAXWRK = M + LWORK_ZGELQF |
|
MAXWRK = MAX( MAXWRK, 3*M + M*M + LWORK_ZGEBRD ) |
|
MAXWRK = MAX( MAXWRK, 3*M + M*M + LWORK_ZUNMBR ) |
|
MAXWRK = MAX( MAXWRK, 3*M + M*M + LWORK_ZUNGBR ) |
IF( NRHS.GT.1 ) THEN |
IF( NRHS.GT.1 ) THEN |
MAXWRK = MAX( MAXWRK, M*M + M + M*NRHS ) |
MAXWRK = MAX( MAXWRK, M*M + M + M*NRHS ) |
ELSE |
ELSE |
MAXWRK = MAX( MAXWRK, M*M + 2*M ) |
MAXWRK = MAX( MAXWRK, M*M + 2*M ) |
END IF |
END IF |
MAXWRK = MAX( MAXWRK, M + NRHS*ILAENV( 1, 'ZUNMLQ', |
MAXWRK = MAX( MAXWRK, M + LWORK_ZUNMLQ ) |
$ 'LC', N, NRHS, M, -1 ) ) |
|
ELSE |
ELSE |
* |
* |
* Path 2 - underdetermined |
* Path 2 - underdetermined |
* |
* |
MAXWRK = 2*M + ( N + M )*ILAENV( 1, 'ZGEBRD', ' ', M, |
* Compute space needed for ZGEBRD |
$ N, -1, -1 ) |
CALL ZGEBRD( M, N, A, LDA, S, S, DUM(1), DUM(1), |
MAXWRK = MAX( MAXWRK, 2*M + NRHS*ILAENV( 1, 'ZUNMBR', |
$ DUM(1), -1, INFO ) |
$ 'QLC', M, NRHS, M, -1 ) ) |
LWORK_ZGEBRD=DUM(1) |
MAXWRK = MAX( MAXWRK, 2*M + M*ILAENV( 1, 'ZUNGBR', |
* Compute space needed for ZUNMBR |
$ 'P', M, N, M, -1 ) ) |
CALL ZUNMBR( 'Q', 'L', 'C', M, NRHS, M, A, LDA, |
|
$ DUM(1), B, LDB, DUM(1), -1, INFO ) |
|
LWORK_ZUNMBR=DUM(1) |
|
* Compute space needed for ZUNGBR |
|
CALL ZUNGBR( 'P', M, N, M, A, LDA, DUM(1), |
|
$ DUM(1), -1, INFO ) |
|
LWORK_ZUNGBR=DUM(1) |
|
MAXWRK = 2*M + LWORK_ZGEBRD |
|
MAXWRK = MAX( MAXWRK, 2*M + LWORK_ZUNMBR ) |
|
MAXWRK = MAX( MAXWRK, 2*M + LWORK_ZUNGBR ) |
MAXWRK = MAX( MAXWRK, N*NRHS ) |
MAXWRK = MAX( MAXWRK, N*NRHS ) |
END IF |
END IF |
END IF |
END IF |
Line 371
|
Line 507
|
* (CWorkspace: none) |
* (CWorkspace: none) |
* (RWorkspace: need BDSPAC) |
* (RWorkspace: need BDSPAC) |
* |
* |
CALL ZBDSQR( 'U', N, N, 0, NRHS, S, RWORK( IE ), A, LDA, VDUM, |
CALL ZBDSQR( 'U', N, N, 0, NRHS, S, RWORK( IE ), A, LDA, DUM, |
$ 1, B, LDB, RWORK( IRWORK ), INFO ) |
$ 1, B, LDB, RWORK( IRWORK ), INFO ) |
IF( INFO.NE.0 ) |
IF( INFO.NE.0 ) |
$ GO TO 70 |
$ GO TO 70 |
Line 568
|
Line 704
|
* (CWorkspace: none) |
* (CWorkspace: none) |
* (RWorkspace: need BDSPAC) |
* (RWorkspace: need BDSPAC) |
* |
* |
CALL ZBDSQR( 'L', M, N, 0, NRHS, S, RWORK( IE ), A, LDA, VDUM, |
CALL ZBDSQR( 'L', M, N, 0, NRHS, S, RWORK( IE ), A, LDA, DUM, |
$ 1, B, LDB, RWORK( IRWORK ), INFO ) |
$ 1, B, LDB, RWORK( IRWORK ), INFO ) |
IF( INFO.NE.0 ) |
IF( INFO.NE.0 ) |
$ GO TO 70 |
$ GO TO 70 |