version 1.2, 2010/04/21 13:45:29
|
version 1.19, 2023/08/07 08:39:17
|
Line 1
|
Line 1
|
|
*> \brief <b> ZGELSD computes the minimum-norm solution to a linear least squares problem for GE matrices</b> |
|
* |
|
* =========== DOCUMENTATION =========== |
|
* |
|
* Online html documentation available at |
|
* http://www.netlib.org/lapack/explore-html/ |
|
* |
|
*> \htmlonly |
|
*> Download ZGELSD + dependencies |
|
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zgelsd.f"> |
|
*> [TGZ]</a> |
|
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zgelsd.f"> |
|
*> [ZIP]</a> |
|
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zgelsd.f"> |
|
*> [TXT]</a> |
|
*> \endhtmlonly |
|
* |
|
* Definition: |
|
* =========== |
|
* |
|
* SUBROUTINE ZGELSD( M, N, NRHS, A, LDA, B, LDB, S, RCOND, RANK, |
|
* WORK, LWORK, RWORK, IWORK, INFO ) |
|
* |
|
* .. Scalar Arguments .. |
|
* INTEGER INFO, LDA, LDB, LWORK, M, N, NRHS, RANK |
|
* DOUBLE PRECISION RCOND |
|
* .. |
|
* .. Array Arguments .. |
|
* INTEGER IWORK( * ) |
|
* DOUBLE PRECISION RWORK( * ), S( * ) |
|
* COMPLEX*16 A( LDA, * ), B( LDB, * ), WORK( * ) |
|
* .. |
|
* |
|
* |
|
*> \par Purpose: |
|
* ============= |
|
*> |
|
*> \verbatim |
|
*> |
|
*> ZGELSD computes the minimum-norm solution to a real 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 problem is solved in three steps: |
|
*> (1) Reduce the coefficient matrix A to bidiagonal form with |
|
*> Householder transformations, reducing the original problem |
|
*> into a "bidiagonal least squares problem" (BLS) |
|
*> (2) Solve the BLS using a divide and conquer approach. |
|
*> (3) Apply back all the Householder transformations to solve |
|
*> the original least squares problem. |
|
*> |
|
*> The effective rank of A is determined by treating as zero those |
|
*> singular values which are less than RCOND times the largest singular |
|
*> value. |
|
*> |
|
*> The divide and conquer algorithm makes very mild assumptions about |
|
*> floating point arithmetic. It will work on machines with a guard |
|
*> digit in add/subtract, or on those binary machines without guard |
|
*> digits which subtract like the Cray X-MP, Cray Y-MP, Cray C-90, or |
|
*> Cray-2. It could conceivably fail on hexadecimal or decimal machines |
|
*> without guard digits, but we know of none. |
|
*> \endverbatim |
|
* |
|
* Arguments: |
|
* ========== |
|
* |
|
*> \param[in] 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, A has been destroyed. |
|
*> \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 must be at least 1. |
|
*> The exact minimum amount of workspace needed depends on M, |
|
*> N and NRHS. As long as LWORK is at least |
|
*> 2*N + N*NRHS |
|
*> if M is greater than or equal to N or |
|
*> 2*M + M*NRHS |
|
*> if M is less than N, the code will execute correctly. |
|
*> 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 array WORK and the |
|
*> minimum sizes of the arrays RWORK and IWORK, and returns |
|
*> these values as the first entries of the WORK, RWORK and |
|
*> IWORK arrays, and no error message related to LWORK is issued |
|
*> by XERBLA. |
|
*> \endverbatim |
|
*> |
|
*> \param[out] RWORK |
|
*> \verbatim |
|
*> RWORK is DOUBLE PRECISION array, dimension (MAX(1,LRWORK)) |
|
*> LRWORK >= |
|
*> 10*N + 2*N*SMLSIZ + 8*N*NLVL + 3*SMLSIZ*NRHS + |
|
*> MAX( (SMLSIZ+1)**2, N*(1+NRHS) + 2*NRHS ) |
|
*> if M is greater than or equal to N or |
|
*> 10*M + 2*M*SMLSIZ + 8*M*NLVL + 3*SMLSIZ*NRHS + |
|
*> MAX( (SMLSIZ+1)**2, N*(1+NRHS) + 2*NRHS ) |
|
*> if M is less than N, the code will execute correctly. |
|
*> SMLSIZ is returned by ILAENV and is equal to the maximum |
|
*> size of the subproblems at the bottom of the computation |
|
*> tree (usually about 25), and |
|
*> NLVL = MAX( 0, INT( LOG_2( MIN( M,N )/(SMLSIZ+1) ) ) + 1 ) |
|
*> On exit, if INFO = 0, RWORK(1) returns the minimum LRWORK. |
|
*> \endverbatim |
|
*> |
|
*> \param[out] IWORK |
|
*> \verbatim |
|
*> IWORK is INTEGER array, dimension (MAX(1,LIWORK)) |
|
*> LIWORK >= max(1, 3*MINMN*NLVL + 11*MINMN), |
|
*> where MINMN = MIN( M,N ). |
|
*> On exit, if INFO = 0, IWORK(1) returns the minimum LIWORK. |
|
*> \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. |
|
* |
|
*> \ingroup complex16GEsolve |
|
* |
|
*> \par Contributors: |
|
* ================== |
|
*> |
|
*> Ming Gu and Ren-Cang Li, Computer Science Division, University of |
|
*> California at Berkeley, USA \n |
|
*> Osni Marques, LBNL/NERSC, USA \n |
|
* |
|
* ===================================================================== |
SUBROUTINE ZGELSD( M, N, NRHS, A, LDA, B, LDB, S, RCOND, RANK, |
SUBROUTINE ZGELSD( M, N, NRHS, A, LDA, B, LDB, S, RCOND, RANK, |
$ WORK, LWORK, RWORK, IWORK, INFO ) |
$ WORK, LWORK, RWORK, IWORK, INFO ) |
* |
* |
* -- LAPACK driver routine (version 3.2) -- |
* -- LAPACK driver routine -- |
* -- LAPACK is a software package provided by Univ. of Tennessee, -- |
* -- LAPACK is a software package provided by Univ. of Tennessee, -- |
* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- |
* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- |
* November 2006 |
|
* |
* |
* .. Scalar Arguments .. |
* .. Scalar Arguments .. |
INTEGER INFO, LDA, LDB, LWORK, M, N, NRHS, RANK |
INTEGER INFO, LDA, LDB, LWORK, M, N, NRHS, RANK |
Line 16
|
Line 237
|
COMPLEX*16 A( LDA, * ), B( LDB, * ), WORK( * ) |
COMPLEX*16 A( LDA, * ), B( LDB, * ), WORK( * ) |
* .. |
* .. |
* |
* |
* Purpose |
|
* ======= |
|
* |
|
* ZGELSD computes the minimum-norm solution to a real 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 problem is solved in three steps: |
|
* (1) Reduce the coefficient matrix A to bidiagonal form with |
|
* Householder tranformations, reducing the original problem |
|
* into a "bidiagonal least squares problem" (BLS) |
|
* (2) Solve the BLS using a divide and conquer approach. |
|
* (3) Apply back all the Householder tranformations to solve |
|
* the original least squares problem. |
|
* |
|
* The effective rank of A is determined by treating as zero those |
|
* singular values which are less than RCOND times the largest singular |
|
* value. |
|
* |
|
* The divide and conquer algorithm makes very mild assumptions about |
|
* floating point arithmetic. It will work on machines with a guard |
|
* digit in add/subtract, or on those binary machines without guard |
|
* digits which subtract like the Cray X-MP, Cray Y-MP, Cray C-90, or |
|
* Cray-2. It could conceivably fail on hexadecimal or decimal machines |
|
* without guard digits, but we know of none. |
|
* |
|
* Arguments |
|
* ========= |
|
* |
|
* 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) COMPLEX*16 array, dimension (LDA,N) |
|
* On entry, the M-by-N matrix A. |
|
* On exit, A has been destroyed. |
|
* |
|
* 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 must be at least 1. |
|
* The exact minimum amount of workspace needed depends on M, |
|
* N and NRHS. As long as LWORK is at least |
|
* 2*N + N*NRHS |
|
* if M is greater than or equal to N or |
|
* 2*M + M*NRHS |
|
* if M is less than N, the code will execute correctly. |
|
* 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 array WORK and the |
|
* minimum sizes of the arrays RWORK and IWORK, and returns |
|
* these values as the first entries of the WORK, RWORK and |
|
* IWORK arrays, and no error message related to LWORK is issued |
|
* by XERBLA. |
|
* |
|
* RWORK (workspace) DOUBLE PRECISION array, dimension (MAX(1,LRWORK)) |
|
* LRWORK >= |
|
* 10*N + 2*N*SMLSIZ + 8*N*NLVL + 3*SMLSIZ*NRHS + |
|
* (SMLSIZ+1)**2 |
|
* if M is greater than or equal to N or |
|
* 10*M + 2*M*SMLSIZ + 8*M*NLVL + 3*SMLSIZ*NRHS + |
|
* (SMLSIZ+1)**2 |
|
* if M is less than N, the code will execute correctly. |
|
* SMLSIZ is returned by ILAENV and is equal to the maximum |
|
* size of the subproblems at the bottom of the computation |
|
* tree (usually about 25), and |
|
* NLVL = MAX( 0, INT( LOG_2( MIN( M,N )/(SMLSIZ+1) ) ) + 1 ) |
|
* On exit, if INFO = 0, RWORK(1) returns the minimum LRWORK. |
|
* |
|
* IWORK (workspace) INTEGER array, dimension (MAX(1,LIWORK)) |
|
* LIWORK >= max(1, 3*MINMN*NLVL + 11*MINMN), |
|
* where MINMN = MIN( M,N ). |
|
* On exit, if INFO = 0, IWORK(1) returns the minimum LIWORK. |
|
* |
|
* 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. |
|
* |
|
* Further Details |
|
* =============== |
|
* |
|
* Based on contributions by |
|
* Ming Gu and Ren-Cang Li, Computer Science Division, University of |
|
* California at Berkeley, USA |
|
* Osni Marques, LBNL/NERSC, USA |
|
* |
|
* ===================================================================== |
* ===================================================================== |
* |
* |
* .. Parameters .. |
* .. Parameters .. |
Line 229
|
Line 320
|
* Path 1 - overdetermined or exactly determined. |
* Path 1 - overdetermined or exactly determined. |
* |
* |
LRWORK = 10*N + 2*N*SMLSIZ + 8*N*NLVL + 3*SMLSIZ*NRHS + |
LRWORK = 10*N + 2*N*SMLSIZ + 8*N*NLVL + 3*SMLSIZ*NRHS + |
$ ( SMLSIZ + 1 )**2 |
$ MAX( (SMLSIZ+1)**2, N*(1+NRHS) + 2*NRHS ) |
MAXWRK = MAX( MAXWRK, 2*N + ( MM + N )*ILAENV( 1, |
MAXWRK = MAX( MAXWRK, 2*N + ( MM + N )*ILAENV( 1, |
$ 'ZGEBRD', ' ', MM, N, -1, -1 ) ) |
$ 'ZGEBRD', ' ', MM, N, -1, -1 ) ) |
MAXWRK = MAX( MAXWRK, 2*N + NRHS*ILAENV( 1, 'ZUNMBR', |
MAXWRK = MAX( MAXWRK, 2*N + NRHS*ILAENV( 1, 'ZUNMBR', |
Line 241
|
Line 332
|
END IF |
END IF |
IF( N.GT.M ) THEN |
IF( N.GT.M ) THEN |
LRWORK = 10*M + 2*M*SMLSIZ + 8*M*NLVL + 3*SMLSIZ*NRHS + |
LRWORK = 10*M + 2*M*SMLSIZ + 8*M*NLVL + 3*SMLSIZ*NRHS + |
$ ( SMLSIZ + 1 )**2 |
$ MAX( (SMLSIZ+1)**2, N*(1+NRHS) + 2*NRHS ) |
IF( N.GE.MNTHR ) THEN |
IF( N.GE.MNTHR ) THEN |
* |
* |
* Path 2a - underdetermined, with many more columns |
* Path 2a - underdetermined, with many more columns |