--- rpl/lapack/lapack/zgelsd.f 2010/12/21 13:53:43 1.8 +++ rpl/lapack/lapack/zgelsd.f 2011/11/21 20:43:09 1.9 @@ -1,10 +1,234 @@ +*> \brief ZGELSD computes the minimum-norm solution to a linear least squares problem for GE matrices +* +* =========== DOCUMENTATION =========== +* +* Online html documentation available at +* http://www.netlib.org/lapack/explore-html/ +* +*> \htmlonly +*> Download ZGELSD + dependencies +*> +*> [TGZ] +*> +*> [ZIP] +*> +*> [TXT] +*> \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 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. +*> \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] 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. +* +*> \date November 2011 +* +*> \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, $ WORK, LWORK, RWORK, IWORK, 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 @@ -16,136 +240,6 @@ 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 + -* 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. -* -* 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 ..