--- rpl/lapack/lapack/zgelss.f 2010/12/21 13:53:43 1.7 +++ rpl/lapack/lapack/zgelss.f 2011/11/21 20:43:09 1.8 @@ -1,10 +1,187 @@ +*> \brief ZGELSS solves overdetermined or underdetermined systems for GE matrices +* +* =========== DOCUMENTATION =========== +* +* Online html documentation available at +* http://www.netlib.org/lapack/explore-html/ +* +*> \htmlonly +*> Download ZGELSS + dependencies +*> +*> [TGZ] +*> +*> [ZIP] +*> +*> [TXT] +*> \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 November 2011 +* +*> \ingroup complex16GEsolve +* +* ===================================================================== SUBROUTINE ZGELSS( M, N, NRHS, A, LDA, B, LDB, S, RCOND, RANK, $ WORK, LWORK, RWORK, INFO ) * -* -- LAPACK driver routine (version 3.2) -- +* -- LAPACK driver routine (version 3.4.0) -- * -- LAPACK is a software package provided by Univ. of Tennessee, -- * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- -* November 2006 +* November 2011 * * .. Scalar Arguments .. INTEGER INFO, LDA, LDB, LWORK, M, N, NRHS, RANK @@ -15,92 +192,6 @@ 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 .. @@ -115,10 +206,13 @@ INTEGER BL, CHUNK, I, IASCL, IBSCL, IE, IL, IRWORK, $ ITAU, ITAUP, ITAUQ, IWORK, LDWORK, MAXMN, $ 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 * .. * .. Local Arrays .. - COMPLEX*16 VDUM( 1 ) + COMPLEX*16 DUM( 1 ) * .. * .. External Subroutines .. EXTERNAL DLABAD, DLASCL, DLASET, XERBLA, ZBDSQR, ZCOPY, @@ -173,6 +267,13 @@ * Path 1a - overdetermined, with many more rows than * 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 MAXWRK = MAX( MAXWRK, N + N*ILAENV( 1, 'ZGEQRF', ' ', M, $ N, -1, -1 ) ) @@ -183,12 +284,22 @@ * * Path 1 - overdetermined or exactly determined * - MAXWRK = MAX( MAXWRK, 2*N + ( MM + N )*ILAENV( 1, - $ 'ZGEBRD', ' ', MM, N, -1, -1 ) ) - MAXWRK = MAX( MAXWRK, 2*N + NRHS*ILAENV( 1, 'ZUNMBR', - $ 'QLC', MM, NRHS, N, -1 ) ) - MAXWRK = MAX( MAXWRK, 2*N + ( N - 1 )*ILAENV( 1, - $ 'ZUNGBR', 'P', N, N, N, -1 ) ) +* Compute space needed for ZGEBRD + CALL ZGEBRD( MM, N, A, LDA, S, DUM(1), DUM(1), + $ DUM(1), DUM(1), -1, INFO ) + LWORK_ZGEBRD=DUM(1) +* Compute space needed for ZUNMBR + 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 ) MINWRK = 2*N + MAX( NRHS, M ) END IF @@ -199,31 +310,56 @@ * Path 2a - underdetermined, with many more columns * than rows * - MAXWRK = M + M*ILAENV( 1, 'ZGELQF', ' ', M, N, -1, - $ -1 ) - MAXWRK = MAX( MAXWRK, 3*M + M*M + 2*M*ILAENV( 1, - $ 'ZGEBRD', ' ', M, M, -1, -1 ) ) - MAXWRK = MAX( MAXWRK, 3*M + M*M + NRHS*ILAENV( 1, - $ 'ZUNMBR', 'QLC', M, NRHS, M, -1 ) ) - MAXWRK = MAX( MAXWRK, 3*M + M*M + ( M - 1 )*ILAENV( 1, - $ 'ZUNGBR', 'P', M, M, M, -1 ) ) +* Compute space needed for ZGELQF + CALL ZGELQF( M, N, A, LDA, DUM(1), DUM(1), + $ -1, INFO ) + LWORK_ZGELQF=DUM(1) +* Compute space needed for ZGEBRD + CALL ZGEBRD( M, M, A, LDA, S, DUM(1), DUM(1), + $ DUM(1), DUM(1), -1, INFO ) + 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 MAXWRK = MAX( MAXWRK, M*M + M + M*NRHS ) ELSE MAXWRK = MAX( MAXWRK, M*M + 2*M ) END IF - MAXWRK = MAX( MAXWRK, M + NRHS*ILAENV( 1, 'ZUNMLQ', - $ 'LC', N, NRHS, M, -1 ) ) + MAXWRK = MAX( MAXWRK, M + LWORK_ZUNMLQ ) ELSE * * Path 2 - underdetermined * - MAXWRK = 2*M + ( N + M )*ILAENV( 1, 'ZGEBRD', ' ', M, - $ N, -1, -1 ) - MAXWRK = MAX( MAXWRK, 2*M + NRHS*ILAENV( 1, 'ZUNMBR', - $ 'QLC', M, NRHS, M, -1 ) ) - MAXWRK = MAX( MAXWRK, 2*M + M*ILAENV( 1, 'ZUNGBR', - $ 'P', M, N, M, -1 ) ) +* Compute space needed for ZGEBRD + CALL ZGEBRD( M, N, A, LDA, S, DUM(1), DUM(1), + $ DUM(1), DUM(1), -1, INFO ) + LWORK_ZGEBRD=DUM(1) +* Compute space needed for ZUNMBR + 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 ) END IF END IF @@ -371,7 +507,7 @@ * (CWorkspace: none) * (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 ) IF( INFO.NE.0 ) $ GO TO 70 @@ -568,7 +704,7 @@ * (CWorkspace: none) * (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 ) IF( INFO.NE.0 ) $ GO TO 70