LAPACK  3.10.0
LAPACK: Linear Algebra PACKage
zggglm.f
Go to the documentation of this file.
1 *> \brief \b ZGGGLM
2 *
3 * =========== DOCUMENTATION ===========
4 *
5 * Online html documentation available at
6 * http://www.netlib.org/lapack/explore-html/
7 *
8 *> \htmlonly
9 *> Download ZGGGLM + dependencies
10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zggglm.f">
11 *> [TGZ]</a>
12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zggglm.f">
13 *> [ZIP]</a>
14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zggglm.f">
15 *> [TXT]</a>
16 *> \endhtmlonly
17 *
18 * Definition:
19 * ===========
20 *
21 * SUBROUTINE ZGGGLM( N, M, P, A, LDA, B, LDB, D, X, Y, WORK, LWORK,
22 * INFO )
23 *
24 * .. Scalar Arguments ..
25 * INTEGER INFO, LDA, LDB, LWORK, M, N, P
26 * ..
27 * .. Array Arguments ..
28 * COMPLEX*16 A( LDA, * ), B( LDB, * ), D( * ), WORK( * ),
29 * $ X( * ), Y( * )
30 * ..
31 *
32 *
33 *> \par Purpose:
34 * =============
35 *>
36 *> \verbatim
37 *>
38 *> ZGGGLM solves a general Gauss-Markov linear model (GLM) problem:
39 *>
40 *> minimize || y ||_2 subject to d = A*x + B*y
41 *> x
42 *>
43 *> where A is an N-by-M matrix, B is an N-by-P matrix, and d is a
44 *> given N-vector. It is assumed that M <= N <= M+P, and
45 *>
46 *> rank(A) = M and rank( A B ) = N.
47 *>
48 *> Under these assumptions, the constrained equation is always
49 *> consistent, and there is a unique solution x and a minimal 2-norm
50 *> solution y, which is obtained using a generalized QR factorization
51 *> of the matrices (A, B) given by
52 *>
53 *> A = Q*(R), B = Q*T*Z.
54 *> (0)
55 *>
56 *> In particular, if matrix B is square nonsingular, then the problem
57 *> GLM is equivalent to the following weighted linear least squares
58 *> problem
59 *>
60 *> minimize || inv(B)*(d-A*x) ||_2
61 *> x
62 *>
63 *> where inv(B) denotes the inverse of B.
64 *> \endverbatim
65 *
66 * Arguments:
67 * ==========
68 *
69 *> \param[in] N
70 *> \verbatim
71 *> N is INTEGER
72 *> The number of rows of the matrices A and B. N >= 0.
73 *> \endverbatim
74 *>
75 *> \param[in] M
76 *> \verbatim
77 *> M is INTEGER
78 *> The number of columns of the matrix A. 0 <= M <= N.
79 *> \endverbatim
80 *>
81 *> \param[in] P
82 *> \verbatim
83 *> P is INTEGER
84 *> The number of columns of the matrix B. P >= N-M.
85 *> \endverbatim
86 *>
87 *> \param[in,out] A
88 *> \verbatim
89 *> A is COMPLEX*16 array, dimension (LDA,M)
90 *> On entry, the N-by-M matrix A.
91 *> On exit, the upper triangular part of the array A contains
92 *> the M-by-M upper triangular matrix R.
93 *> \endverbatim
94 *>
95 *> \param[in] LDA
96 *> \verbatim
97 *> LDA is INTEGER
98 *> The leading dimension of the array A. LDA >= max(1,N).
99 *> \endverbatim
100 *>
101 *> \param[in,out] B
102 *> \verbatim
103 *> B is COMPLEX*16 array, dimension (LDB,P)
104 *> On entry, the N-by-P matrix B.
105 *> On exit, if N <= P, the upper triangle of the subarray
106 *> B(1:N,P-N+1:P) contains the N-by-N upper triangular matrix T;
107 *> if N > P, the elements on and above the (N-P)th subdiagonal
108 *> contain the N-by-P upper trapezoidal matrix T.
109 *> \endverbatim
110 *>
111 *> \param[in] LDB
112 *> \verbatim
113 *> LDB is INTEGER
114 *> The leading dimension of the array B. LDB >= max(1,N).
115 *> \endverbatim
116 *>
117 *> \param[in,out] D
118 *> \verbatim
119 *> D is COMPLEX*16 array, dimension (N)
120 *> On entry, D is the left hand side of the GLM equation.
121 *> On exit, D is destroyed.
122 *> \endverbatim
123 *>
124 *> \param[out] X
125 *> \verbatim
126 *> X is COMPLEX*16 array, dimension (M)
127 *> \endverbatim
128 *>
129 *> \param[out] Y
130 *> \verbatim
131 *> Y is COMPLEX*16 array, dimension (P)
132 *>
133 *> On exit, X and Y are the solutions of the GLM problem.
134 *> \endverbatim
135 *>
136 *> \param[out] WORK
137 *> \verbatim
138 *> WORK is COMPLEX*16 array, dimension (MAX(1,LWORK))
139 *> On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
140 *> \endverbatim
141 *>
142 *> \param[in] LWORK
143 *> \verbatim
144 *> LWORK is INTEGER
145 *> The dimension of the array WORK. LWORK >= max(1,N+M+P).
146 *> For optimum performance, LWORK >= M+min(N,P)+max(N,P)*NB,
147 *> where NB is an upper bound for the optimal blocksizes for
148 *> ZGEQRF, ZGERQF, ZUNMQR and ZUNMRQ.
149 *>
150 *> If LWORK = -1, then a workspace query is assumed; the routine
151 *> only calculates the optimal size of the WORK array, returns
152 *> this value as the first entry of the WORK array, and no error
153 *> message related to LWORK is issued by XERBLA.
154 *> \endverbatim
155 *>
156 *> \param[out] INFO
157 *> \verbatim
158 *> INFO is INTEGER
159 *> = 0: successful exit.
160 *> < 0: if INFO = -i, the i-th argument had an illegal value.
161 *> = 1: the upper triangular factor R associated with A in the
162 *> generalized QR factorization of the pair (A, B) is
163 *> singular, so that rank(A) < M; the least squares
164 *> solution could not be computed.
165 *> = 2: the bottom (N-M) by (N-M) part of the upper trapezoidal
166 *> factor T associated with B in the generalized QR
167 *> factorization of the pair (A, B) is singular, so that
168 *> rank( A B ) < N; the least squares solution could not
169 *> be computed.
170 *> \endverbatim
171 *
172 * Authors:
173 * ========
174 *
175 *> \author Univ. of Tennessee
176 *> \author Univ. of California Berkeley
177 *> \author Univ. of Colorado Denver
178 *> \author NAG Ltd.
179 *
180 *> \ingroup complex16OTHEReigen
181 *
182 * =====================================================================
183  SUBROUTINE zggglm( N, M, P, A, LDA, B, LDB, D, X, Y, WORK, LWORK,
184  $ INFO )
185 *
186 * -- LAPACK driver routine --
187 * -- LAPACK is a software package provided by Univ. of Tennessee, --
188 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
189 *
190 * .. Scalar Arguments ..
191  INTEGER INFO, LDA, LDB, LWORK, M, N, P
192 * ..
193 * .. Array Arguments ..
194  COMPLEX*16 A( LDA, * ), B( LDB, * ), D( * ), WORK( * ),
195  $ x( * ), y( * )
196 * ..
197 *
198 * ===================================================================
199 *
200 * .. Parameters ..
201  COMPLEX*16 CZERO, CONE
202  parameter( czero = ( 0.0d+0, 0.0d+0 ),
203  $ cone = ( 1.0d+0, 0.0d+0 ) )
204 * ..
205 * .. Local Scalars ..
206  LOGICAL LQUERY
207  INTEGER I, LOPT, LWKMIN, LWKOPT, NB, NB1, NB2, NB3,
208  $ nb4, np
209 * ..
210 * .. External Subroutines ..
211  EXTERNAL xerbla, zcopy, zgemv, zggqrf, ztrtrs, zunmqr,
212  $ zunmrq
213 * ..
214 * .. External Functions ..
215  INTEGER ILAENV
216  EXTERNAL ilaenv
217 * ..
218 * .. Intrinsic Functions ..
219  INTRINSIC int, max, min
220 * ..
221 * .. Executable Statements ..
222 *
223 * Test the input parameters
224 *
225  info = 0
226  np = min( n, p )
227  lquery = ( lwork.EQ.-1 )
228  IF( n.LT.0 ) THEN
229  info = -1
230  ELSE IF( m.LT.0 .OR. m.GT.n ) THEN
231  info = -2
232  ELSE IF( p.LT.0 .OR. p.LT.n-m ) THEN
233  info = -3
234  ELSE IF( lda.LT.max( 1, n ) ) THEN
235  info = -5
236  ELSE IF( ldb.LT.max( 1, n ) ) THEN
237  info = -7
238  END IF
239 *
240 * Calculate workspace
241 *
242  IF( info.EQ.0) THEN
243  IF( n.EQ.0 ) THEN
244  lwkmin = 1
245  lwkopt = 1
246  ELSE
247  nb1 = ilaenv( 1, 'ZGEQRF', ' ', n, m, -1, -1 )
248  nb2 = ilaenv( 1, 'ZGERQF', ' ', n, m, -1, -1 )
249  nb3 = ilaenv( 1, 'ZUNMQR', ' ', n, m, p, -1 )
250  nb4 = ilaenv( 1, 'ZUNMRQ', ' ', n, m, p, -1 )
251  nb = max( nb1, nb2, nb3, nb4 )
252  lwkmin = m + n + p
253  lwkopt = m + np + max( n, p )*nb
254  END IF
255  work( 1 ) = lwkopt
256 *
257  IF( lwork.LT.lwkmin .AND. .NOT.lquery ) THEN
258  info = -12
259  END IF
260  END IF
261 *
262  IF( info.NE.0 ) THEN
263  CALL xerbla( 'ZGGGLM', -info )
264  RETURN
265  ELSE IF( lquery ) THEN
266  RETURN
267  END IF
268 *
269 * Quick return if possible
270 *
271  IF( n.EQ.0 ) THEN
272  DO i = 1, m
273  x(i) = czero
274  END DO
275  DO i = 1, p
276  y(i) = czero
277  END DO
278  RETURN
279  END IF
280 *
281 * Compute the GQR factorization of matrices A and B:
282 *
283 * Q**H*A = ( R11 ) M, Q**H*B*Z**H = ( T11 T12 ) M
284 * ( 0 ) N-M ( 0 T22 ) N-M
285 * M M+P-N N-M
286 *
287 * where R11 and T22 are upper triangular, and Q and Z are
288 * unitary.
289 *
290  CALL zggqrf( n, m, p, a, lda, work, b, ldb, work( m+1 ),
291  $ work( m+np+1 ), lwork-m-np, info )
292  lopt = dble( work( m+np+1 ) )
293 *
294 * Update left-hand-side vector d = Q**H*d = ( d1 ) M
295 * ( d2 ) N-M
296 *
297  CALL zunmqr( 'Left', 'Conjugate transpose', n, 1, m, a, lda, work,
298  $ d, max( 1, n ), work( m+np+1 ), lwork-m-np, info )
299  lopt = max( lopt, int( work( m+np+1 ) ) )
300 *
301 * Solve T22*y2 = d2 for y2
302 *
303  IF( n.GT.m ) THEN
304  CALL ztrtrs( 'Upper', 'No transpose', 'Non unit', n-m, 1,
305  $ b( m+1, m+p-n+1 ), ldb, d( m+1 ), n-m, info )
306 *
307  IF( info.GT.0 ) THEN
308  info = 1
309  RETURN
310  END IF
311 *
312  CALL zcopy( n-m, d( m+1 ), 1, y( m+p-n+1 ), 1 )
313  END IF
314 *
315 * Set y1 = 0
316 *
317  DO 10 i = 1, m + p - n
318  y( i ) = czero
319  10 CONTINUE
320 *
321 * Update d1 = d1 - T12*y2
322 *
323  CALL zgemv( 'No transpose', m, n-m, -cone, b( 1, m+p-n+1 ), ldb,
324  $ y( m+p-n+1 ), 1, cone, d, 1 )
325 *
326 * Solve triangular system: R11*x = d1
327 *
328  IF( m.GT.0 ) THEN
329  CALL ztrtrs( 'Upper', 'No Transpose', 'Non unit', m, 1, a, lda,
330  $ d, m, info )
331 *
332  IF( info.GT.0 ) THEN
333  info = 2
334  RETURN
335  END IF
336 *
337 * Copy D to X
338 *
339  CALL zcopy( m, d, 1, x, 1 )
340  END IF
341 *
342 * Backward transformation y = Z**H *y
343 *
344  CALL zunmrq( 'Left', 'Conjugate transpose', p, 1, np,
345  $ b( max( 1, n-p+1 ), 1 ), ldb, work( m+1 ), y,
346  $ max( 1, p ), work( m+np+1 ), lwork-m-np, info )
347  work( 1 ) = m + np + max( lopt, int( work( m+np+1 ) ) )
348 *
349  RETURN
350 *
351 * End of ZGGGLM
352 *
353  END
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:60
subroutine zcopy(N, ZX, INCX, ZY, INCY)
ZCOPY
Definition: zcopy.f:81
subroutine zgemv(TRANS, M, N, ALPHA, A, LDA, X, INCX, BETA, Y, INCY)
ZGEMV
Definition: zgemv.f:158
subroutine ztrtrs(UPLO, TRANS, DIAG, N, NRHS, A, LDA, B, LDB, INFO)
ZTRTRS
Definition: ztrtrs.f:140
subroutine zggqrf(N, M, P, A, LDA, TAUA, B, LDB, TAUB, WORK, LWORK, INFO)
ZGGQRF
Definition: zggqrf.f:215
subroutine zunmrq(SIDE, TRANS, M, N, K, A, LDA, TAU, C, LDC, WORK, LWORK, INFO)
ZUNMRQ
Definition: zunmrq.f:167
subroutine zunmqr(SIDE, TRANS, M, N, K, A, LDA, TAU, C, LDC, WORK, LWORK, INFO)
ZUNMQR
Definition: zunmqr.f:167
subroutine zggglm(N, M, P, A, LDA, B, LDB, D, X, Y, WORK, LWORK, INFO)
ZGGGLM
Definition: zggglm.f:185