--- rpl/lapack/lapack/zggglm.f 2010/12/21 13:53:45 1.7
+++ rpl/lapack/lapack/zggglm.f 2023/08/07 08:39:21 1.19
@@ -1,10 +1,191 @@
+*> \brief \b ZGGGLM
+*
+* =========== DOCUMENTATION ===========
+*
+* Online html documentation available at
+* http://www.netlib.org/lapack/explore-html/
+*
+*> \htmlonly
+*> Download ZGGGLM + dependencies
+*>
+*> [TGZ]
+*>
+*> [ZIP]
+*>
+*> [TXT]
+*> \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.
+*
+*> \ingroup complex16OTHEReigen
+*
+* =====================================================================
SUBROUTINE ZGGGLM( N, M, P, A, LDA, B, LDB, D, X, Y, WORK, LWORK,
$ INFO )
*
-* -- LAPACK driver routine (version 3.2) --
+* -- LAPACK driver routine --
* -- LAPACK is a software package provided by Univ. of Tennessee, --
* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
-* November 2006
*
* .. Scalar Arguments ..
INTEGER INFO, LDA, LDB, LWORK, M, N, P
@@ -14,101 +195,6 @@
$ 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 ..
@@ -127,7 +213,7 @@
* ..
* .. External Functions ..
INTEGER ILAENV
- EXTERNAL ILAENV
+ EXTERNAL ILAENV
* ..
* .. Intrinsic Functions ..
INTRINSIC INT, MAX, MIN
@@ -182,24 +268,31 @@
*
* Quick return if possible
*
- IF( N.EQ.0 )
- $ RETURN
+ IF( N.EQ.0 ) THEN
+ DO I = 1, M
+ X(I) = CZERO
+ END DO
+ DO I = 1, P
+ Y(I) = CZERO
+ END DO
+ RETURN
+ END IF
*
* Compute the GQR factorization of matrices A and B:
*
-* Q'*A = ( R11 ) M, Q'*B*Z' = ( T11 T12 ) M
-* ( 0 ) N-M ( 0 T22 ) N-M
-* M M+P-N N-M
+* Q**H*A = ( R11 ) M, Q**H*B*Z**H = ( T11 T12 ) M
+* ( 0 ) N-M ( 0 T22 ) N-M
+* M M+P-N N-M
*
* where R11 and T22 are upper triangular, and Q and Z are
* unitary.
*
CALL ZGGQRF( N, M, P, A, LDA, WORK, B, LDB, WORK( M+1 ),
$ WORK( M+NP+1 ), LWORK-M-NP, INFO )
- LOPT = WORK( M+NP+1 )
+ LOPT = INT( WORK( M+NP+1 ) )
*
-* Update left-hand-side vector d = Q'*d = ( d1 ) M
-* ( d2 ) N-M
+* Update left-hand-side vector d = Q**H*d = ( d1 ) M
+* ( d2 ) N-M
*
CALL ZUNMQR( 'Left', 'Conjugate transpose', N, 1, M, A, LDA, WORK,
$ D, MAX( 1, N ), WORK( M+NP+1 ), LWORK-M-NP, INFO )
@@ -246,7 +339,7 @@
CALL ZCOPY( M, D, 1, X, 1 )
END IF
*
-* Backward transformation y = Z'*y
+* Backward transformation y = Z**H *y
*
CALL ZUNMRQ( 'Left', 'Conjugate transpose', P, 1, NP,
$ B( MAX( 1, N-P+1 ), 1 ), LDB, WORK( M+1 ), Y,