version 1.3, 2010/08/06 15:28:51
|
version 1.12, 2012/12/14 14:22:45
|
Line 1
|
Line 1
|
|
*> \brief <b> ZGELS solves overdetermined or underdetermined systems for GE matrices</b> |
|
* |
|
* =========== DOCUMENTATION =========== |
|
* |
|
* Online html documentation available at |
|
* http://www.netlib.org/lapack/explore-html/ |
|
* |
|
*> \htmlonly |
|
*> Download ZGELS + dependencies |
|
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zgels.f"> |
|
*> [TGZ]</a> |
|
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zgels.f"> |
|
*> [ZIP]</a> |
|
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zgels.f"> |
|
*> [TXT]</a> |
|
*> \endhtmlonly |
|
* |
|
* Definition: |
|
* =========== |
|
* |
|
* SUBROUTINE ZGELS( TRANS, M, N, NRHS, A, LDA, B, LDB, WORK, LWORK, |
|
* INFO ) |
|
* |
|
* .. Scalar Arguments .. |
|
* CHARACTER TRANS |
|
* INTEGER INFO, LDA, LDB, LWORK, M, N, NRHS |
|
* .. |
|
* .. Array Arguments .. |
|
* COMPLEX*16 A( LDA, * ), B( LDB, * ), WORK( * ) |
|
* .. |
|
* |
|
* |
|
*> \par Purpose: |
|
* ============= |
|
*> |
|
*> \verbatim |
|
*> |
|
*> ZGELS solves overdetermined or underdetermined complex linear systems |
|
*> involving an M-by-N matrix A, or its conjugate-transpose, using a QR |
|
*> or LQ factorization of A. It is assumed that A has full rank. |
|
*> |
|
*> The following options are provided: |
|
*> |
|
*> 1. If TRANS = 'N' and m >= n: find the least squares solution of |
|
*> an overdetermined system, i.e., solve the least squares problem |
|
*> minimize || B - A*X ||. |
|
*> |
|
*> 2. If TRANS = 'N' and m < n: find the minimum norm solution of |
|
*> an underdetermined system A * X = B. |
|
*> |
|
*> 3. If TRANS = 'C' and m >= n: find the minimum norm solution of |
|
*> an undetermined system A**H * X = B. |
|
*> |
|
*> 4. If TRANS = 'C' and m < n: find the least squares solution of |
|
*> an overdetermined system, i.e., solve the least squares problem |
|
*> minimize || B - A**H * X ||. |
|
*> |
|
*> 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. |
|
*> \endverbatim |
|
* |
|
* Arguments: |
|
* ========== |
|
* |
|
*> \param[in] TRANS |
|
*> \verbatim |
|
*> TRANS is CHARACTER*1 |
|
*> = 'N': the linear system involves A; |
|
*> = 'C': the linear system involves A**H. |
|
*> \endverbatim |
|
*> |
|
*> \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. |
|
*> if M >= N, A is overwritten by details of its QR |
|
*> factorization as returned by ZGEQRF; |
|
*> if M < N, A is overwritten by details of its LQ |
|
*> factorization as returned by ZGELQF. |
|
*> \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 matrix B of right hand side vectors, stored |
|
*> columnwise; B is M-by-NRHS if TRANS = 'N', or N-by-NRHS |
|
*> if TRANS = 'C'. |
|
*> On exit, if INFO = 0, B is overwritten by the solution |
|
*> vectors, stored columnwise: |
|
*> if TRANS = 'N' and m >= n, rows 1 to n of B contain the least |
|
*> squares solution vectors; the residual sum of squares for the |
|
*> solution in each column is given by the sum of squares of the |
|
*> modulus of elements N+1 to M in that column; |
|
*> if TRANS = 'N' and m < n, rows 1 to N of B contain the |
|
*> minimum norm solution vectors; |
|
*> if TRANS = 'C' and m >= n, rows 1 to M of B contain the |
|
*> minimum norm solution vectors; |
|
*> if TRANS = 'C' and m < n, rows 1 to M of B contain the |
|
*> least squares solution vectors; the residual sum of squares |
|
*> for the solution in each column is given by the sum of |
|
*> squares of the modulus of elements M+1 to N 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] 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 >= max( 1, MN + max( MN, NRHS ) ). |
|
*> For optimal performance, |
|
*> LWORK >= max( 1, MN + max( MN, NRHS )*NB ). |
|
*> where MN = min(M,N) and NB is the optimum block size. |
|
*> |
|
*> 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] INFO |
|
*> \verbatim |
|
*> INFO is INTEGER |
|
*> = 0: successful exit |
|
*> < 0: if INFO = -i, the i-th argument had an illegal value |
|
*> > 0: if INFO = i, the i-th diagonal element of the |
|
*> triangular factor of A is zero, so that A does not have |
|
*> full rank; the least squares solution could not be |
|
*> computed. |
|
*> \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 ZGELS( TRANS, M, N, NRHS, A, LDA, B, LDB, WORK, LWORK, |
SUBROUTINE ZGELS( TRANS, M, N, NRHS, A, LDA, B, LDB, WORK, LWORK, |
$ INFO ) |
$ 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 .. |
CHARACTER TRANS |
CHARACTER TRANS |
Line 14
|
Line 195
|
COMPLEX*16 A( LDA, * ), B( LDB, * ), WORK( * ) |
COMPLEX*16 A( LDA, * ), B( LDB, * ), WORK( * ) |
* .. |
* .. |
* |
* |
* Purpose |
|
* ======= |
|
* |
|
* ZGELS solves overdetermined or underdetermined complex linear systems |
|
* involving an M-by-N matrix A, or its conjugate-transpose, using a QR |
|
* or LQ factorization of A. It is assumed that A has full rank. |
|
* |
|
* The following options are provided: |
|
* |
|
* 1. If TRANS = 'N' and m >= n: find the least squares solution of |
|
* an overdetermined system, i.e., solve the least squares problem |
|
* minimize || B - A*X ||. |
|
* |
|
* 2. If TRANS = 'N' and m < n: find the minimum norm solution of |
|
* an underdetermined system A * X = B. |
|
* |
|
* 3. If TRANS = 'C' and m >= n: find the minimum norm solution of |
|
* an undetermined system A**H * X = B. |
|
* |
|
* 4. If TRANS = 'C' and m < n: find the least squares solution of |
|
* an overdetermined system, i.e., solve the least squares problem |
|
* minimize || B - A**H * X ||. |
|
* |
|
* 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. |
|
* |
|
* Arguments |
|
* ========= |
|
* |
|
* TRANS (input) CHARACTER*1 |
|
* = 'N': the linear system involves A; |
|
* = 'C': the linear system involves A**H. |
|
* |
|
* 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. |
|
* if M >= N, A is overwritten by details of its QR |
|
* factorization as returned by ZGEQRF; |
|
* if M < N, A is overwritten by details of its LQ |
|
* factorization as returned by ZGELQF. |
|
* |
|
* 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 matrix B of right hand side vectors, stored |
|
* columnwise; B is M-by-NRHS if TRANS = 'N', or N-by-NRHS |
|
* if TRANS = 'C'. |
|
* On exit, if INFO = 0, B is overwritten by the solution |
|
* vectors, stored columnwise: |
|
* if TRANS = 'N' and m >= n, rows 1 to n of B contain the least |
|
* squares solution vectors; the residual sum of squares for the |
|
* solution in each column is given by the sum of squares of the |
|
* modulus of elements N+1 to M in that column; |
|
* if TRANS = 'N' and m < n, rows 1 to N of B contain the |
|
* minimum norm solution vectors; |
|
* if TRANS = 'C' and m >= n, rows 1 to M of B contain the |
|
* minimum norm solution vectors; |
|
* if TRANS = 'C' and m < n, rows 1 to M of B contain the |
|
* least squares solution vectors; the residual sum of squares |
|
* for the solution in each column is given by the sum of |
|
* squares of the modulus of elements M+1 to N in that column. |
|
* |
|
* LDB (input) INTEGER |
|
* The leading dimension of the array B. LDB >= MAX(1,M,N). |
|
* |
|
* 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 >= max( 1, MN + max( MN, NRHS ) ). |
|
* For optimal performance, |
|
* LWORK >= max( 1, MN + max( MN, NRHS )*NB ). |
|
* where MN = min(M,N) and NB is the optimum block size. |
|
* |
|
* 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. |
|
* |
|
* INFO (output) INTEGER |
|
* = 0: successful exit |
|
* < 0: if INFO = -i, the i-th argument had an illegal value |
|
* > 0: if INFO = i, the i-th diagonal element of the |
|
* triangular factor of A is zero, so that A does not have |
|
* full rank; the least squares solution could not be |
|
* computed. |
|
* |
|
* ===================================================================== |
* ===================================================================== |
* |
* |
* .. Parameters .. |
* .. Parameters .. |
Line 278
|
Line 359
|
* |
* |
* Least-Squares Problem min || A * X - B || |
* Least-Squares Problem min || A * X - B || |
* |
* |
* 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, N, A, |
CALL ZUNMQR( 'Left', 'Conjugate transpose', M, NRHS, N, A, |
$ LDA, WORK( 1 ), B, LDB, WORK( MN+1 ), LWORK-MN, |
$ LDA, WORK( 1 ), B, LDB, WORK( MN+1 ), LWORK-MN, |
Line 299
|
Line 380
|
* |
* |
ELSE |
ELSE |
* |
* |
* Overdetermined system of equations A' * X = B |
* Overdetermined system of equations A**H * X = B |
* |
* |
* B(1:N,1:NRHS) := inv(R') * B(1:N,1:NRHS) |
* B(1:N,1:NRHS) := inv(R**H) * B(1:N,1:NRHS) |
* |
* |
CALL ZTRTRS( 'Upper', 'Conjugate transpose','Non-unit', |
CALL ZTRTRS( 'Upper', 'Conjugate transpose','Non-unit', |
$ N, NRHS, A, LDA, B, LDB, INFO ) |
$ N, NRHS, A, LDA, B, LDB, INFO ) |
Line 360
|
Line 441
|
30 CONTINUE |
30 CONTINUE |
40 CONTINUE |
40 CONTINUE |
* |
* |
* B(1:N,1:NRHS) := Q(1:N,:)' * B(1:M,1:NRHS) |
* B(1:N,1:NRHS) := Q(1:N,:)**H * B(1:M,1:NRHS) |
* |
* |
CALL ZUNMLQ( 'Left', 'Conjugate transpose', N, NRHS, M, A, |
CALL ZUNMLQ( 'Left', 'Conjugate transpose', N, NRHS, M, A, |
$ LDA, WORK( 1 ), B, LDB, WORK( MN+1 ), LWORK-MN, |
$ LDA, WORK( 1 ), B, LDB, WORK( MN+1 ), LWORK-MN, |
Line 372
|
Line 453
|
* |
* |
ELSE |
ELSE |
* |
* |
* overdetermined system min || A' * X - B || |
* overdetermined system min || A**H * X - B || |
* |
* |
* B(1:N,1:NRHS) := Q * B(1:N,1:NRHS) |
* B(1:N,1:NRHS) := Q * B(1:N,1:NRHS) |
* |
* |
Line 382
|
Line 463
|
* |
* |
* workspace at least NRHS, optimally NRHS*NB |
* workspace at least NRHS, optimally NRHS*NB |
* |
* |
* B(1:M,1:NRHS) := inv(L') * B(1:M,1:NRHS) |
* B(1:M,1:NRHS) := inv(L**H) * B(1:M,1:NRHS) |
* |
* |
CALL ZTRTRS( 'Lower', 'Conjugate transpose', 'Non-unit', |
CALL ZTRTRS( 'Lower', 'Conjugate transpose', 'Non-unit', |
$ M, NRHS, A, LDA, B, LDB, INFO ) |
$ M, NRHS, A, LDA, B, LDB, INFO ) |