LAPACK 3.12.0
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches
dggglm.f
Go to the documentation of this file.
1*> \brief \b DGGGLM
2*
3* =========== DOCUMENTATION ===========
4*
5* Online html documentation available at
6* http://www.netlib.org/lapack/explore-html/
7*
8*> \htmlonly
9*> Download DGGGLM + dependencies
10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dggglm.f">
11*> [TGZ]</a>
12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dggglm.f">
13*> [ZIP]</a>
14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dggglm.f">
15*> [TXT]</a>
16*> \endhtmlonly
17*
18* Definition:
19* ===========
20*
21* SUBROUTINE DGGGLM( 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* DOUBLE PRECISION A( LDA, * ), B( LDB, * ), D( * ), WORK( * ),
29* $ X( * ), Y( * )
30* ..
31*
32*
33*> \par Purpose:
34* =============
35*>
36*> \verbatim
37*>
38*> DGGGLM 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION array, dimension (M)
127*> \endverbatim
128*>
129*> \param[out] Y
130*> \verbatim
131*> Y is DOUBLE PRECISION 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 DOUBLE PRECISION 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*> DGEQRF, SGERQF, DORMQR and SORMRQ.
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 ggglm
181*
182* =====================================================================
183 SUBROUTINE dggglm( 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 DOUBLE PRECISION A( LDA, * ), B( LDB, * ), D( * ), WORK( * ),
195 $ x( * ), y( * )
196* ..
197*
198* ===================================================================
199*
200* .. Parameters ..
201 DOUBLE PRECISION ZERO, ONE
202 parameter( zero = 0.0d+0, one = 1.0d+0 )
203* ..
204* .. Local Scalars ..
205 LOGICAL LQUERY
206 INTEGER I, LOPT, LWKMIN, LWKOPT, NB, NB1, NB2, NB3,
207 $ nb4, np
208* ..
209* .. External Subroutines ..
210 EXTERNAL dcopy, dgemv, dggqrf, dormqr, dormrq, dtrtrs,
211 $ xerbla
212* ..
213* .. External Functions ..
214 INTEGER ILAENV
215 EXTERNAL ilaenv
216* ..
217* .. Intrinsic Functions ..
218 INTRINSIC int, max, min
219* ..
220* .. Executable Statements ..
221*
222* Test the input parameters
223*
224 info = 0
225 np = min( n, p )
226 lquery = ( lwork.EQ.-1 )
227 IF( n.LT.0 ) THEN
228 info = -1
229 ELSE IF( m.LT.0 .OR. m.GT.n ) THEN
230 info = -2
231 ELSE IF( p.LT.0 .OR. p.LT.n-m ) THEN
232 info = -3
233 ELSE IF( lda.LT.max( 1, n ) ) THEN
234 info = -5
235 ELSE IF( ldb.LT.max( 1, n ) ) THEN
236 info = -7
237 END IF
238*
239* Calculate workspace
240*
241 IF( info.EQ.0) THEN
242 IF( n.EQ.0 ) THEN
243 lwkmin = 1
244 lwkopt = 1
245 ELSE
246 nb1 = ilaenv( 1, 'DGEQRF', ' ', n, m, -1, -1 )
247 nb2 = ilaenv( 1, 'DGERQF', ' ', n, m, -1, -1 )
248 nb3 = ilaenv( 1, 'DORMQR', ' ', n, m, p, -1 )
249 nb4 = ilaenv( 1, 'DORMRQ', ' ', n, m, p, -1 )
250 nb = max( nb1, nb2, nb3, nb4 )
251 lwkmin = m + n + p
252 lwkopt = m + np + max( n, p )*nb
253 END IF
254 work( 1 ) = lwkopt
255*
256 IF( lwork.LT.lwkmin .AND. .NOT.lquery ) THEN
257 info = -12
258 END IF
259 END IF
260*
261 IF( info.NE.0 ) THEN
262 CALL xerbla( 'DGGGLM', -info )
263 RETURN
264 ELSE IF( lquery ) THEN
265 RETURN
266 END IF
267*
268* Quick return if possible
269*
270 IF( n.EQ.0 ) THEN
271 DO i = 1, m
272 x(i) = zero
273 END DO
274 DO i = 1, p
275 y(i) = zero
276 END DO
277 RETURN
278 END IF
279*
280* Compute the GQR factorization of matrices A and B:
281*
282* Q**T*A = ( R11 ) M, Q**T*B*Z**T = ( T11 T12 ) M
283* ( 0 ) N-M ( 0 T22 ) N-M
284* M M+P-N N-M
285*
286* where R11 and T22 are upper triangular, and Q and Z are
287* orthogonal.
288*
289 CALL dggqrf( n, m, p, a, lda, work, b, ldb, work( m+1 ),
290 $ work( m+np+1 ), lwork-m-np, info )
291 lopt = int( work( m+np+1 ) )
292*
293* Update left-hand-side vector d = Q**T*d = ( d1 ) M
294* ( d2 ) N-M
295*
296 CALL dormqr( 'Left', 'Transpose', n, 1, m, a, lda, work, d,
297 $ max( 1, n ), work( m+np+1 ), lwork-m-np, info )
298 lopt = max( lopt, int( work( m+np+1 ) ) )
299*
300* Solve T22*y2 = d2 for y2
301*
302 IF( n.GT.m ) THEN
303 CALL dtrtrs( 'Upper', 'No transpose', 'Non unit', n-m, 1,
304 $ b( m+1, m+p-n+1 ), ldb, d( m+1 ), n-m, info )
305*
306 IF( info.GT.0 ) THEN
307 info = 1
308 RETURN
309 END IF
310*
311 CALL dcopy( n-m, d( m+1 ), 1, y( m+p-n+1 ), 1 )
312 END IF
313*
314* Set y1 = 0
315*
316 DO 10 i = 1, m + p - n
317 y( i ) = zero
318 10 CONTINUE
319*
320* Update d1 = d1 - T12*y2
321*
322 CALL dgemv( 'No transpose', m, n-m, -one, b( 1, m+p-n+1 ), ldb,
323 $ y( m+p-n+1 ), 1, one, d, 1 )
324*
325* Solve triangular system: R11*x = d1
326*
327 IF( m.GT.0 ) THEN
328 CALL dtrtrs( 'Upper', 'No Transpose', 'Non unit', m, 1, a, lda,
329 $ d, m, info )
330*
331 IF( info.GT.0 ) THEN
332 info = 2
333 RETURN
334 END IF
335*
336* Copy D to X
337*
338 CALL dcopy( m, d, 1, x, 1 )
339 END IF
340*
341* Backward transformation y = Z**T *y
342*
343 CALL dormrq( 'Left', 'Transpose', p, 1, np,
344 $ b( max( 1, n-p+1 ), 1 ), ldb, work( m+1 ), y,
345 $ max( 1, p ), work( m+np+1 ), lwork-m-np, info )
346 work( 1 ) = m + np + max( lopt, int( work( m+np+1 ) ) )
347*
348 RETURN
349*
350* End of DGGGLM
351*
352 END
subroutine xerbla(srname, info)
Definition cblat2.f:3285
subroutine dcopy(n, dx, incx, dy, incy)
DCOPY
Definition dcopy.f:82
subroutine dgemv(trans, m, n, alpha, a, lda, x, incx, beta, y, incy)
DGEMV
Definition dgemv.f:158
subroutine dggglm(n, m, p, a, lda, b, ldb, d, x, y, work, lwork, info)
DGGGLM
Definition dggglm.f:185
subroutine dggqrf(n, m, p, a, lda, taua, b, ldb, taub, work, lwork, info)
DGGQRF
Definition dggqrf.f:215
subroutine dtrtrs(uplo, trans, diag, n, nrhs, a, lda, b, ldb, info)
DTRTRS
Definition dtrtrs.f:140
subroutine dormqr(side, trans, m, n, k, a, lda, tau, c, ldc, work, lwork, info)
DORMQR
Definition dormqr.f:167
subroutine dormrq(side, trans, m, n, k, a, lda, tau, c, ldc, work, lwork, info)
DORMRQ
Definition dormrq.f:167