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