version 1.6, 2010/08/13 21:04:04
|
version 1.14, 2015/11/26 11:44:22
|
Line 1
|
Line 1
|
|
*> \brief \b ZGGGLM |
|
* |
|
* =========== DOCUMENTATION =========== |
|
* |
|
* Online html documentation available at |
|
* http://www.netlib.org/lapack/explore-html/ |
|
* |
|
*> \htmlonly |
|
*> Download ZGGGLM + dependencies |
|
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zggglm.f"> |
|
*> [TGZ]</a> |
|
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zggglm.f"> |
|
*> [ZIP]</a> |
|
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zggglm.f"> |
|
*> [TXT]</a> |
|
*> \endhtmlonly |
|
* |
|
* Definition: |
|
* =========== |
|
* |
|
* SUBROUTINE ZGGGLM( N, M, P, A, LDA, B, LDB, D, X, Y, WORK, LWORK, |
|
* INFO ) |
|
* |
|
* .. Scalar Arguments .. |
|
* INTEGER INFO, LDA, LDB, LWORK, M, N, P |
|
* .. |
|
* .. Array Arguments .. |
|
* COMPLEX*16 A( LDA, * ), B( LDB, * ), D( * ), WORK( * ), |
|
* $ X( * ), Y( * ) |
|
* .. |
|
* |
|
* |
|
*> \par Purpose: |
|
* ============= |
|
*> |
|
*> \verbatim |
|
*> |
|
*> ZGGGLM solves a general Gauss-Markov linear model (GLM) problem: |
|
*> |
|
*> minimize || y ||_2 subject to d = A*x + B*y |
|
*> x |
|
*> |
|
*> where A is an N-by-M matrix, B is an N-by-P matrix, and d is a |
|
*> given N-vector. It is assumed that M <= N <= M+P, and |
|
*> |
|
*> rank(A) = M and rank( A B ) = N. |
|
*> |
|
*> Under these assumptions, the constrained equation is always |
|
*> consistent, and there is a unique solution x and a minimal 2-norm |
|
*> solution y, which is obtained using a generalized QR factorization |
|
*> of the matrices (A, B) given by |
|
*> |
|
*> A = Q*(R), B = Q*T*Z. |
|
*> (0) |
|
*> |
|
*> In particular, if matrix B is square nonsingular, then the problem |
|
*> GLM is equivalent to the following weighted linear least squares |
|
*> problem |
|
*> |
|
*> minimize || inv(B)*(d-A*x) ||_2 |
|
*> x |
|
*> |
|
*> where inv(B) denotes the inverse of B. |
|
*> \endverbatim |
|
* |
|
* Arguments: |
|
* ========== |
|
* |
|
*> \param[in] N |
|
*> \verbatim |
|
*> N is INTEGER |
|
*> The number of rows of the matrices A and B. N >= 0. |
|
*> \endverbatim |
|
*> |
|
*> \param[in] M |
|
*> \verbatim |
|
*> M is INTEGER |
|
*> The number of columns of the matrix A. 0 <= M <= N. |
|
*> \endverbatim |
|
*> |
|
*> \param[in] P |
|
*> \verbatim |
|
*> P is INTEGER |
|
*> The number of columns of the matrix B. P >= N-M. |
|
*> \endverbatim |
|
*> |
|
*> \param[in,out] A |
|
*> \verbatim |
|
*> A is COMPLEX*16 array, dimension (LDA,M) |
|
*> On entry, the N-by-M matrix A. |
|
*> On exit, the upper triangular part of the array A contains |
|
*> the M-by-M upper triangular matrix R. |
|
*> \endverbatim |
|
*> |
|
*> \param[in] LDA |
|
*> \verbatim |
|
*> LDA is INTEGER |
|
*> The leading dimension of the array A. LDA >= max(1,N). |
|
*> \endverbatim |
|
*> |
|
*> \param[in,out] B |
|
*> \verbatim |
|
*> B is COMPLEX*16 array, dimension (LDB,P) |
|
*> On entry, the N-by-P matrix B. |
|
*> On exit, if N <= P, the upper triangle of the subarray |
|
*> B(1:N,P-N+1:P) contains the N-by-N upper triangular matrix T; |
|
*> if N > P, the elements on and above the (N-P)th subdiagonal |
|
*> contain the N-by-P upper trapezoidal matrix T. |
|
*> \endverbatim |
|
*> |
|
*> \param[in] LDB |
|
*> \verbatim |
|
*> LDB is INTEGER |
|
*> The leading dimension of the array B. LDB >= max(1,N). |
|
*> \endverbatim |
|
*> |
|
*> \param[in,out] D |
|
*> \verbatim |
|
*> D is COMPLEX*16 array, dimension (N) |
|
*> On entry, D is the left hand side of the GLM equation. |
|
*> On exit, D is destroyed. |
|
*> \endverbatim |
|
*> |
|
*> \param[out] X |
|
*> \verbatim |
|
*> X is COMPLEX*16 array, dimension (M) |
|
*> \endverbatim |
|
*> |
|
*> \param[out] Y |
|
*> \verbatim |
|
*> Y is COMPLEX*16 array, dimension (P) |
|
*> |
|
*> On exit, X and Y are the solutions of the GLM problem. |
|
*> \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,N+M+P). |
|
*> For optimum performance, LWORK >= M+min(N,P)+max(N,P)*NB, |
|
*> where NB is an upper bound for the optimal blocksizes for |
|
*> ZGEQRF, ZGERQF, ZUNMQR and ZUNMRQ. |
|
*> |
|
*> 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. |
|
*> = 1: the upper triangular factor R associated with A in the |
|
*> generalized QR factorization of the pair (A, B) is |
|
*> singular, so that rank(A) < M; the least squares |
|
*> solution could not be computed. |
|
*> = 2: the bottom (N-M) by (N-M) part of the upper trapezoidal |
|
*> factor T associated with B in the generalized QR |
|
*> factorization of the pair (A, B) is singular, so that |
|
*> rank( A B ) < N; 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 2015 |
|
* |
|
*> \ingroup complex16OTHEReigen |
|
* |
|
* ===================================================================== |
SUBROUTINE ZGGGLM( N, M, P, A, LDA, B, LDB, D, X, Y, WORK, LWORK, |
SUBROUTINE ZGGGLM( N, M, P, A, LDA, B, LDB, D, X, Y, WORK, LWORK, |
$ INFO ) |
$ INFO ) |
* |
* |
* -- LAPACK driver routine (version 3.2) -- |
* -- LAPACK driver routine (version 3.6.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 2015 |
* |
* |
* .. Scalar Arguments .. |
* .. Scalar Arguments .. |
INTEGER INFO, LDA, LDB, LWORK, M, N, P |
INTEGER INFO, LDA, LDB, LWORK, M, N, P |
Line 14
|
Line 198
|
$ X( * ), Y( * ) |
$ X( * ), Y( * ) |
* .. |
* .. |
* |
* |
* Purpose |
|
* ======= |
|
* |
|
* ZGGGLM solves a general Gauss-Markov linear model (GLM) problem: |
|
* |
|
* minimize || y ||_2 subject to d = A*x + B*y |
|
* x |
|
* |
|
* where A is an N-by-M matrix, B is an N-by-P matrix, and d is a |
|
* given N-vector. It is assumed that M <= N <= M+P, and |
|
* |
|
* rank(A) = M and rank( A B ) = N. |
|
* |
|
* Under these assumptions, the constrained equation is always |
|
* consistent, and there is a unique solution x and a minimal 2-norm |
|
* solution y, which is obtained using a generalized QR factorization |
|
* of the matrices (A, B) given by |
|
* |
|
* A = Q*(R), B = Q*T*Z. |
|
* (0) |
|
* |
|
* In particular, if matrix B is square nonsingular, then the problem |
|
* GLM is equivalent to the following weighted linear least squares |
|
* problem |
|
* |
|
* minimize || inv(B)*(d-A*x) ||_2 |
|
* x |
|
* |
|
* where inv(B) denotes the inverse of B. |
|
* |
|
* Arguments |
|
* ========= |
|
* |
|
* N (input) INTEGER |
|
* The number of rows of the matrices A and B. N >= 0. |
|
* |
|
* M (input) INTEGER |
|
* The number of columns of the matrix A. 0 <= M <= N. |
|
* |
|
* P (input) INTEGER |
|
* The number of columns of the matrix B. P >= N-M. |
|
* |
|
* A (input/output) COMPLEX*16 array, dimension (LDA,M) |
|
* On entry, the N-by-M matrix A. |
|
* On exit, the upper triangular part of the array A contains |
|
* the M-by-M upper triangular matrix R. |
|
* |
|
* LDA (input) INTEGER |
|
* The leading dimension of the array A. LDA >= max(1,N). |
|
* |
|
* B (input/output) COMPLEX*16 array, dimension (LDB,P) |
|
* On entry, the N-by-P matrix B. |
|
* On exit, if N <= P, the upper triangle of the subarray |
|
* B(1:N,P-N+1:P) contains the N-by-N upper triangular matrix T; |
|
* if N > P, the elements on and above the (N-P)th subdiagonal |
|
* contain the N-by-P upper trapezoidal matrix T. |
|
* |
|
* LDB (input) INTEGER |
|
* The leading dimension of the array B. LDB >= max(1,N). |
|
* |
|
* D (input/output) COMPLEX*16 array, dimension (N) |
|
* On entry, D is the left hand side of the GLM equation. |
|
* On exit, D is destroyed. |
|
* |
|
* X (output) COMPLEX*16 array, dimension (M) |
|
* Y (output) COMPLEX*16 array, dimension (P) |
|
* On exit, X and Y are the solutions of the GLM problem. |
|
* |
|
* 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,N+M+P). |
|
* For optimum performance, LWORK >= M+min(N,P)+max(N,P)*NB, |
|
* where NB is an upper bound for the optimal blocksizes for |
|
* ZGEQRF, ZGERQF, ZUNMQR and ZUNMRQ. |
|
* |
|
* 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. |
|
* = 1: the upper triangular factor R associated with A in the |
|
* generalized QR factorization of the pair (A, B) is |
|
* singular, so that rank(A) < M; the least squares |
|
* solution could not be computed. |
|
* = 2: the bottom (N-M) by (N-M) part of the upper trapezoidal |
|
* factor T associated with B in the generalized QR |
|
* factorization of the pair (A, B) is singular, so that |
|
* rank( A B ) < N; the least squares solution could not |
|
* be computed. |
|
* |
|
* =================================================================== |
* =================================================================== |
* |
* |
* .. Parameters .. |
* .. Parameters .. |
Line 187
|
Line 276
|
* |
* |
* Compute the GQR factorization of matrices A and B: |
* Compute the GQR factorization of matrices A and B: |
* |
* |
* Q'*A = ( R11 ) M, Q'*B*Z' = ( T11 T12 ) M |
* Q**H*A = ( R11 ) M, Q**H*B*Z**H = ( T11 T12 ) M |
* ( 0 ) N-M ( 0 T22 ) N-M |
* ( 0 ) N-M ( 0 T22 ) N-M |
* M M+P-N N-M |
* M M+P-N N-M |
* |
* |
* where R11 and T22 are upper triangular, and Q and Z are |
* where R11 and T22 are upper triangular, and Q and Z are |
* unitary. |
* unitary. |
Line 198
|
Line 287
|
$ WORK( M+NP+1 ), LWORK-M-NP, INFO ) |
$ WORK( M+NP+1 ), LWORK-M-NP, INFO ) |
LOPT = WORK( M+NP+1 ) |
LOPT = WORK( M+NP+1 ) |
* |
* |
* Update left-hand-side vector d = Q'*d = ( d1 ) M |
* Update left-hand-side vector d = Q**H*d = ( d1 ) M |
* ( d2 ) N-M |
* ( d2 ) N-M |
* |
* |
CALL ZUNMQR( 'Left', 'Conjugate transpose', N, 1, M, A, LDA, WORK, |
CALL ZUNMQR( 'Left', 'Conjugate transpose', N, 1, M, A, LDA, WORK, |
$ D, MAX( 1, N ), WORK( M+NP+1 ), LWORK-M-NP, INFO ) |
$ D, MAX( 1, N ), WORK( M+NP+1 ), LWORK-M-NP, INFO ) |
Line 246
|
Line 335
|
CALL ZCOPY( M, D, 1, X, 1 ) |
CALL ZCOPY( M, D, 1, X, 1 ) |
END IF |
END IF |
* |
* |
* Backward transformation y = Z'*y |
* Backward transformation y = Z**H *y |
* |
* |
CALL ZUNMRQ( 'Left', 'Conjugate transpose', P, 1, NP, |
CALL ZUNMRQ( 'Left', 'Conjugate transpose', P, 1, NP, |
$ B( MAX( 1, N-P+1 ), 1 ), LDB, WORK( M+1 ), Y, |
$ B( MAX( 1, N-P+1 ), 1 ), LDB, WORK( M+1 ), Y, |