--- rpl/lapack/lapack/zgelsy.f 2011/07/22 07:38:14 1.8 +++ rpl/lapack/lapack/zgelsy.f 2011/11/21 20:43:09 1.9 @@ -1,10 +1,219 @@ +*> \brief ZGELSY solves overdetermined or underdetermined systems for GE matrices +* +* =========== DOCUMENTATION =========== +* +* Online html documentation available at +* http://www.netlib.org/lapack/explore-html/ +* +*> \htmlonly +*> Download ZGELSY + dependencies +*> +*> [TGZ] +*> +*> [ZIP] +*> +*> [TXT] +*> \endhtmlonly +* +* Definition: +* =========== +* +* SUBROUTINE ZGELSY( M, N, NRHS, A, LDA, B, LDB, JPVT, RCOND, RANK, +* WORK, LWORK, RWORK, INFO ) +* +* .. Scalar Arguments .. +* INTEGER INFO, LDA, LDB, LWORK, M, N, NRHS, RANK +* DOUBLE PRECISION RCOND +* .. +* .. Array Arguments .. +* INTEGER JPVT( * ) +* DOUBLE PRECISION RWORK( * ) +* COMPLEX*16 A( LDA, * ), B( LDB, * ), WORK( * ) +* .. +* +* +*> \par Purpose: +* ============= +*> +*> \verbatim +*> +*> ZGELSY computes the minimum-norm solution to a complex linear least +*> squares problem: +*> minimize || A * X - B || +*> using a complete orthogonal factorization 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 routine first computes a QR factorization with column pivoting: +*> A * P = Q * [ R11 R12 ] +*> [ 0 R22 ] +*> with R11 defined as the largest leading submatrix whose estimated +*> condition number is less than 1/RCOND. The order of R11, RANK, +*> is the effective rank of A. +*> +*> Then, R22 is considered to be negligible, and R12 is annihilated +*> by unitary transformations from the right, arriving at the +*> complete orthogonal factorization: +*> A * P = Q * [ T11 0 ] * Z +*> [ 0 0 ] +*> The minimum-norm solution is then +*> X = P * Z**H [ inv(T11)*Q1**H*B ] +*> [ 0 ] +*> where Q1 consists of the first RANK columns of Q. +*> +*> This routine is basically identical to the original xGELSX except +*> three differences: +*> o The permutation of matrix B (the right hand side) is faster and +*> more simple. +*> o The call to the subroutine xGEQPF has been substituted by the +*> the call to the subroutine xGEQP3. This subroutine is a Blas-3 +*> version of the QR factorization with column pivoting. +*> o Matrix B (the right hand side) is updated with Blas-3. +*> \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 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 overwritten by details of its +*> complete orthogonal factorization. +*> \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, the N-by-NRHS solution matrix X. +*> \endverbatim +*> +*> \param[in] LDB +*> \verbatim +*> LDB is INTEGER +*> The leading dimension of the array B. LDB >= max(1,M,N). +*> \endverbatim +*> +*> \param[in,out] JPVT +*> \verbatim +*> JPVT is INTEGER array, dimension (N) +*> On entry, if JPVT(i) .ne. 0, the i-th column of A is permuted +*> to the front of AP, otherwise column i is a free column. +*> On exit, if JPVT(i) = k, then the i-th column of A*P +*> was the k-th column of A. +*> \endverbatim +*> +*> \param[in] RCOND +*> \verbatim +*> RCOND is DOUBLE PRECISION +*> RCOND is used to determine the effective rank of A, which +*> is defined as the order of the largest leading triangular +*> submatrix R11 in the QR factorization with pivoting of A, +*> whose estimated condition number < 1/RCOND. +*> \endverbatim +*> +*> \param[out] RANK +*> \verbatim +*> RANK is INTEGER +*> The effective rank of A, i.e., the order of the submatrix +*> R11. This is the same as the order of the submatrix T11 +*> in the complete orthogonal factorization of A. +*> \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. +*> The unblocked strategy requires that: +*> LWORK >= MN + MAX( 2*MN, N+1, MN+NRHS ) +*> where MN = min(M,N). +*> The block algorithm requires that: +*> LWORK >= MN + MAX( 2*MN, NB*(N+1), MN+MN*NB, MN+NB*NRHS ) +*> where NB is an upper bound on the blocksize returned +*> by ILAENV for the routines ZGEQP3, ZTZRZF, CTZRQF, ZUNMQR, +*> and ZUNMRZ. +*> +*> 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 (2*N) +*> \endverbatim +*> +*> \param[out] INFO +*> \verbatim +*> INFO is INTEGER +*> = 0: successful exit +*> < 0: if INFO = -i, the i-th argument had an illegal value +*> \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: +* ================== +*> +*> A. Petitet, Computer Science Dept., Univ. of Tenn., Knoxville, USA \n +*> E. Quintana-Orti, Depto. de Informatica, Universidad Jaime I, Spain \n +*> G. Quintana-Orti, Depto. de Informatica, Universidad Jaime I, Spain \n +*> +* ===================================================================== SUBROUTINE ZGELSY( M, N, NRHS, A, LDA, B, LDB, JPVT, RCOND, RANK, $ WORK, LWORK, RWORK, INFO ) * -* -- LAPACK driver routine (version 3.3.1) -- +* -- 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..-- -* -- April 2011 -- +* November 2011 * * .. Scalar Arguments .. INTEGER INFO, LDA, LDB, LWORK, M, N, NRHS, RANK @@ -16,124 +225,6 @@ COMPLEX*16 A( LDA, * ), B( LDB, * ), WORK( * ) * .. * -* Purpose -* ======= -* -* ZGELSY computes the minimum-norm solution to a complex linear least -* squares problem: -* minimize || A * X - B || -* using a complete orthogonal factorization 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 routine first computes a QR factorization with column pivoting: -* A * P = Q * [ R11 R12 ] -* [ 0 R22 ] -* with R11 defined as the largest leading submatrix whose estimated -* condition number is less than 1/RCOND. The order of R11, RANK, -* is the effective rank of A. -* -* Then, R22 is considered to be negligible, and R12 is annihilated -* by unitary transformations from the right, arriving at the -* complete orthogonal factorization: -* A * P = Q * [ T11 0 ] * Z -* [ 0 0 ] -* The minimum-norm solution is then -* X = P * Z**H [ inv(T11)*Q1**H*B ] -* [ 0 ] -* where Q1 consists of the first RANK columns of Q. -* -* This routine is basically identical to the original xGELSX except -* three differences: -* o The permutation of matrix B (the right hand side) is faster and -* more simple. -* o The call to the subroutine xGEQPF has been substituted by the -* the call to the subroutine xGEQP3. This subroutine is a Blas-3 -* version of the QR factorization with column pivoting. -* o Matrix B (the right hand side) is updated with Blas-3. -* -* 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 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, A has been overwritten by details of its -* complete orthogonal factorization. -* -* 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, the N-by-NRHS solution matrix X. -* -* LDB (input) INTEGER -* The leading dimension of the array B. LDB >= max(1,M,N). -* -* JPVT (input/output) INTEGER array, dimension (N) -* On entry, if JPVT(i) .ne. 0, the i-th column of A is permuted -* to the front of AP, otherwise column i is a free column. -* On exit, if JPVT(i) = k, then the i-th column of A*P -* was the k-th column of A. -* -* RCOND (input) DOUBLE PRECISION -* RCOND is used to determine the effective rank of A, which -* is defined as the order of the largest leading triangular -* submatrix R11 in the QR factorization with pivoting of A, -* whose estimated condition number < 1/RCOND. -* -* RANK (output) INTEGER -* The effective rank of A, i.e., the order of the submatrix -* R11. This is the same as the order of the submatrix T11 -* in the complete orthogonal factorization of A. -* -* 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. -* The unblocked strategy requires that: -* LWORK >= MN + MAX( 2*MN, N+1, MN+NRHS ) -* where MN = min(M,N). -* The block algorithm requires that: -* LWORK >= MN + MAX( 2*MN, NB*(N+1), MN+MN*NB, MN+NB*NRHS ) -* where NB is an upper bound on the blocksize returned -* by ILAENV for the routines ZGEQP3, ZTZRZF, CTZRQF, ZUNMQR, -* and ZUNMRZ. -* -* 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 (2*N) -* -* INFO (output) INTEGER -* = 0: successful exit -* < 0: if INFO = -i, the i-th argument had an illegal value -* -* Further Details -* =============== -* -* Based on contributions by -* A. Petitet, Computer Science Dept., Univ. of Tenn., Knoxville, USA -* E. Quintana-Orti, Depto. de Informatica, Universidad Jaime I, Spain -* G. Quintana-Orti, Depto. de Informatica, Universidad Jaime I, Spain -* * ===================================================================== * * .. Parameters ..