version 1.6, 2010/08/13 21:04:03
|
version 1.10, 2011/11/21 22:19:45
|
Line 1
|
Line 1
|
|
*> \brief <b> ZGELSY solves overdetermined or underdetermined systems for GE matrices</b> |
|
* |
|
* =========== DOCUMENTATION =========== |
|
* |
|
* Online html documentation available at |
|
* http://www.netlib.org/lapack/explore-html/ |
|
* |
|
*> \htmlonly |
|
*> Download ZGELSY + dependencies |
|
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zgelsy.f"> |
|
*> [TGZ]</a> |
|
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zgelsy.f"> |
|
*> [ZIP]</a> |
|
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zgelsy.f"> |
|
*> [TXT]</a> |
|
*> \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, |
SUBROUTINE ZGELSY( M, N, NRHS, A, LDA, B, LDB, JPVT, RCOND, RANK, |
$ WORK, LWORK, RWORK, INFO ) |
$ 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, -- |
* -- 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 |
* November 2011 |
* |
* |
* .. 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 225
|
COMPLEX*16 A( LDA, * ), B( LDB, * ), WORK( * ) |
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' [ inv(T11)*Q1'*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 .. |
* .. Parameters .. |
Line 319
|
Line 410
|
* complex workspace: 2*MN. |
* complex workspace: 2*MN. |
* Details of Householder rotations stored in WORK(MN+1:2*MN) |
* Details of Householder rotations stored in WORK(MN+1:2*MN) |
* |
* |
* B(1:M,1:NRHS) := Q' * B(1:M,1:NRHS) |
* B(1:M,1:NRHS) := Q**H * B(1:M,1:NRHS) |
* |
* |
CALL ZUNMQR( 'Left', 'Conjugate transpose', M, NRHS, MN, A, LDA, |
CALL ZUNMQR( 'Left', 'Conjugate transpose', M, NRHS, MN, A, LDA, |
$ WORK( 1 ), B, LDB, WORK( 2*MN+1 ), LWORK-2*MN, INFO ) |
$ WORK( 1 ), B, LDB, WORK( 2*MN+1 ), LWORK-2*MN, INFO ) |
Line 338
|
Line 429
|
30 CONTINUE |
30 CONTINUE |
40 CONTINUE |
40 CONTINUE |
* |
* |
* B(1:N,1:NRHS) := Y' * B(1:N,1:NRHS) |
* B(1:N,1:NRHS) := Y**H * B(1:N,1:NRHS) |
* |
* |
IF( RANK.LT.N ) THEN |
IF( RANK.LT.N ) THEN |
CALL ZUNMRZ( 'Left', 'Conjugate transpose', N, NRHS, RANK, |
CALL ZUNMRZ( 'Left', 'Conjugate transpose', N, NRHS, RANK, |