LAPACK 3.12.0 LAPACK: Linear Algebra PACKage
Searching...
No Matches
sggglm.f
Go to the documentation of this file.
1*> \brief \b SGGGLM
2*
3* =========== DOCUMENTATION ===========
4*
5* Online html documentation available at
6* http://www.netlib.org/lapack/explore-html/
7*
8*> \htmlonly
10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/sggglm.f">
11*> [TGZ]</a>
12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/sggglm.f">
13*> [ZIP]</a>
14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/sggglm.f">
15*> [TXT]</a>
16*> \endhtmlonly
17*
18* Definition:
19* ===========
20*
21* SUBROUTINE SGGGLM( 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* REAL A( LDA, * ), B( LDB, * ), D( * ), WORK( * ),
29* \$ X( * ), Y( * )
30* ..
31*
32*
33*> \par Purpose:
34* =============
35*>
36*> \verbatim
37*>
38*> SGGGLM 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 REAL 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 REAL 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 REAL 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 REAL array, dimension (M)
127*> \endverbatim
128*>
129*> \param[out] Y
130*> \verbatim
131*> Y is REAL 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 REAL 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*> SGEQRF, SGERQF, SORMQR 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 sggglm( 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 REAL A( LDA, * ), B( LDB, * ), D( * ), WORK( * ),
195 \$ x( * ), y( * )
196* ..
197*
198* ===================================================================
199*
200* .. Parameters ..
201 REAL ZERO, ONE
202 parameter( zero = 0.0e+0, one = 1.0e+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 scopy, sgemv, sggqrf, sormqr, sormrq, strtrs,
211 \$ xerbla
212* ..
213* .. External Functions ..
214 INTEGER ILAENV
215 REAL SROUNDUP_LWORK
216 EXTERNAL ilaenv, sroundup_lwork
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, 'SGEQRF', ' ', n, m, -1, -1 )
248 nb2 = ilaenv( 1, 'SGERQF', ' ', n, m, -1, -1 )
249 nb3 = ilaenv( 1, 'SORMQR', ' ', n, m, p, -1 )
250 nb4 = ilaenv( 1, 'SORMRQ', ' ', 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 ) = sroundup_lwork(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( 'SGGGLM', -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) = zero
274 END DO
275 DO i = 1, p
276 y(i) = zero
277 END DO
278 RETURN
279 END IF
280*
281* Compute the GQR factorization of matrices A and B:
282*
283* Q**T*A = ( R11 ) M, Q**T*B*Z**T = ( 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* orthogonal.
289*
290 CALL sggqrf( n, m, p, a, lda, work, b, ldb, work( m+1 ),
291 \$ work( m+np+1 ), lwork-m-np, info )
292 lopt = int( work( m+np+1 ) )
293*
294* Update left-hand-side vector d = Q**T*d = ( d1 ) M
295* ( d2 ) N-M
296*
297 CALL sormqr( 'Left', 'Transpose', n, 1, m, a, lda, work, d,
298 \$ 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 strtrs( '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 scopy( 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 ) = zero
319 10 CONTINUE
320*
321* Update d1 = d1 - T12*y2
322*
323 CALL sgemv( 'No transpose', m, n-m, -one, b( 1, m+p-n+1 ), ldb,
324 \$ y( m+p-n+1 ), 1, one, d, 1 )
325*
326* Solve triangular system: R11*x = d1
327*
328 IF( m.GT.0 ) THEN
329 CALL strtrs( '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 scopy( m, d, 1, x, 1 )
340 END IF
341*
342* Backward transformation y = Z**T *y
343*
344 CALL sormrq( 'Left', '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 SGGGLM
352*
353 END
subroutine xerbla(srname, info)
Definition cblat2.f:3285
subroutine scopy(n, sx, incx, sy, incy)
SCOPY
Definition scopy.f:82
subroutine sgemv(trans, m, n, alpha, a, lda, x, incx, beta, y, incy)
SGEMV
Definition sgemv.f:158
subroutine sggglm(n, m, p, a, lda, b, ldb, d, x, y, work, lwork, info)
SGGGLM
Definition sggglm.f:185
subroutine sggqrf(n, m, p, a, lda, taua, b, ldb, taub, work, lwork, info)
SGGQRF
Definition sggqrf.f:215
subroutine strtrs(uplo, trans, diag, n, nrhs, a, lda, b, ldb, info)
STRTRS
Definition strtrs.f:140
subroutine sormqr(side, trans, m, n, k, a, lda, tau, c, ldc, work, lwork, info)
SORMQR
Definition sormqr.f:168
subroutine sormrq(side, trans, m, n, k, a, lda, tau, c, ldc, work, lwork, info)
SORMRQ
Definition sormrq.f:168