LAPACK 3.12.1
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches

◆ dggglm()

subroutine dggglm ( integer n,
integer m,
integer p,
double precision, dimension( lda, * ) a,
integer lda,
double precision, dimension( ldb, * ) b,
integer ldb,
double precision, dimension( * ) d,
double precision, dimension( * ) x,
double precision, dimension( * ) y,
double precision, dimension( * ) work,
integer lwork,
integer info )

DGGGLM

Download DGGGLM + dependencies [TGZ] [ZIP] [TXT]

Purpose:
!>
!> DGGGLM 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.
!>
!> Callers of this subroutine should note that the singularity/rank-deficiency checks
!> implemented in this subroutine are rudimentary. The DTRTRS subroutine called by this
!> subroutine only signals a failure due to singularity if the problem is exactly singular.
!>
!> It is conceivable for one (or more) of the factors involved in the generalized QR
!> factorization of the pair (A, B) to be subnormally close to singularity without this
!> subroutine signalling an error. The solutions computed for such almost-rank-deficient
!> problems may be less accurate due to a loss of numerical precision.
!> 
!>
Parameters
[in]N
!>          N is INTEGER
!>          The number of rows of the matrices A and B.  N >= 0.
!>
[in]M
!>          M is INTEGER
!>          The number of columns of the matrix A.  0 <= M <= N.
!>
[in]P
!>          P is INTEGER
!>          The number of columns of the matrix B.  P >= N-M.
!>
[in,out]A
!>          A is DOUBLE PRECISION 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.
!>
[in]LDA
!>          LDA is INTEGER
!>          The leading dimension of the array A. LDA >= max(1,N).
!>
[in,out]B
!>          B is DOUBLE PRECISION 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.
!>
[in]LDB
!>          LDB is INTEGER
!>          The leading dimension of the array B. LDB >= max(1,N).
!>
[in,out]D
!>          D is DOUBLE PRECISION array, dimension (N)
!>          On entry, D is the left hand side of the GLM equation.
!>          On exit, D is destroyed.
!>
[out]X
!>          X is DOUBLE PRECISION array, dimension (M)
!>
[out]Y
!>          Y is DOUBLE PRECISION array, dimension (P)
!>
!>          On exit, X and Y are the solutions of the GLM problem.
!>
[out]WORK
!>          WORK is DOUBLE PRECISION array, dimension (MAX(1,LWORK))
!>          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
!>
[in]LWORK
!>          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
!>          DGEQRF, SGERQF, DORMQR and SORMRQ.
!>
!>          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.
!>
[out]INFO
!>          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 exactly
!>                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 exactly singular, so that
!>                rank( A B ) < N; the least squares solution could not
!>                be computed.
!>
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.

Definition at line 191 of file dggglm.f.

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 DOUBLE PRECISION A( LDA, * ), B( LDB, * ), D( * ), WORK( * ),
204 $ X( * ), Y( * )
205* ..
206*
207* ===================================================================
208*
209* .. Parameters ..
210 DOUBLE PRECISION ZERO, ONE
211 parameter( zero = 0.0d+0, one = 1.0d+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 dcopy, dgemv, dggqrf, dormqr, dormrq,
220 $ dtrtrs,
221 $ xerbla
222* ..
223* .. External Functions ..
224 INTEGER ILAENV
225 EXTERNAL ilaenv
226* ..
227* .. Intrinsic Functions ..
228 INTRINSIC int, max, min
229* ..
230* .. Executable Statements ..
231*
232* Test the input parameters
233*
234 info = 0
235 np = min( n, p )
236 lquery = ( lwork.EQ.-1 )
237 IF( n.LT.0 ) THEN
238 info = -1
239 ELSE IF( m.LT.0 .OR. m.GT.n ) THEN
240 info = -2
241 ELSE IF( p.LT.0 .OR. p.LT.n-m ) THEN
242 info = -3
243 ELSE IF( lda.LT.max( 1, n ) ) THEN
244 info = -5
245 ELSE IF( ldb.LT.max( 1, n ) ) THEN
246 info = -7
247 END IF
248*
249* Calculate workspace
250*
251 IF( info.EQ.0) THEN
252 IF( n.EQ.0 ) THEN
253 lwkmin = 1
254 lwkopt = 1
255 ELSE
256 nb1 = ilaenv( 1, 'DGEQRF', ' ', n, m, -1, -1 )
257 nb2 = ilaenv( 1, 'DGERQF', ' ', n, m, -1, -1 )
258 nb3 = ilaenv( 1, 'DORMQR', ' ', n, m, p, -1 )
259 nb4 = ilaenv( 1, 'DORMRQ', ' ', n, m, p, -1 )
260 nb = max( nb1, nb2, nb3, nb4 )
261 lwkmin = m + n + p
262 lwkopt = m + np + max( n, p )*nb
263 END IF
264 work( 1 ) = lwkopt
265*
266 IF( lwork.LT.lwkmin .AND. .NOT.lquery ) THEN
267 info = -12
268 END IF
269 END IF
270*
271 IF( info.NE.0 ) THEN
272 CALL xerbla( 'DGGGLM', -info )
273 RETURN
274 ELSE IF( lquery ) THEN
275 RETURN
276 END IF
277*
278* Quick return if possible
279*
280 IF( n.EQ.0 ) THEN
281 DO i = 1, m
282 x(i) = zero
283 END DO
284 DO i = 1, p
285 y(i) = zero
286 END DO
287 RETURN
288 END IF
289*
290* Compute the GQR factorization of matrices A and B:
291*
292* Q**T*A = ( R11 ) M, Q**T*B*Z**T = ( T11 T12 ) M
293* ( 0 ) N-M ( 0 T22 ) N-M
294* M M+P-N N-M
295*
296* where R11 and T22 are upper triangular, and Q and Z are
297* orthogonal.
298*
299 CALL dggqrf( n, m, p, a, lda, work, b, ldb, work( m+1 ),
300 $ work( m+np+1 ), lwork-m-np, info )
301 lopt = int( work( m+np+1 ) )
302*
303* Update left-hand-side vector d = Q**T*d = ( d1 ) M
304* ( d2 ) N-M
305*
306 CALL dormqr( 'Left', 'Transpose', n, 1, m, a, lda, work, d,
307 $ max( 1, n ), work( m+np+1 ), lwork-m-np, info )
308 lopt = max( lopt, int( work( m+np+1 ) ) )
309*
310* Solve T22*y2 = d2 for y2
311*
312 IF( n.GT.m ) THEN
313 CALL dtrtrs( 'Upper', 'No transpose', 'Non unit', n-m, 1,
314 $ b( m+1, m+p-n+1 ), ldb, d( m+1 ), n-m, info )
315*
316 IF( info.GT.0 ) THEN
317 info = 1
318 RETURN
319 END IF
320*
321 CALL dcopy( n-m, d( m+1 ), 1, y( m+p-n+1 ), 1 )
322 END IF
323*
324* Set y1 = 0
325*
326 DO 10 i = 1, m + p - n
327 y( i ) = zero
328 10 CONTINUE
329*
330* Update d1 = d1 - T12*y2
331*
332 CALL dgemv( 'No transpose', m, n-m, -one, b( 1, m+p-n+1 ), ldb,
333 $ y( m+p-n+1 ), 1, one, d, 1 )
334*
335* Solve triangular system: R11*x = d1
336*
337 IF( m.GT.0 ) THEN
338 CALL dtrtrs( 'Upper', 'No Transpose', 'Non unit', m, 1, a,
339 $ lda,
340 $ d, m, info )
341*
342 IF( info.GT.0 ) THEN
343 info = 2
344 RETURN
345 END IF
346*
347* Copy D to X
348*
349 CALL dcopy( m, d, 1, x, 1 )
350 END IF
351*
352* Backward transformation y = Z**T *y
353*
354 CALL dormrq( 'Left', 'Transpose', p, 1, np,
355 $ b( max( 1, n-p+1 ), 1 ), ldb, work( m+1 ), y,
356 $ max( 1, p ), work( m+np+1 ), lwork-m-np, info )
357 work( 1 ) = m + np + max( lopt, int( work( m+np+1 ) ) )
358*
359 RETURN
360*
361* End of DGGGLM
362*
xerbla
subroutine xerbla(srname, info)
Definition cblat2.f:3285
dcopy
subroutine dcopy(n, dx, incx, dy, incy)
DCOPY
Definition dcopy.f:82
dgemv
subroutine dgemv(trans, m, n, alpha, a, lda, x, incx, beta, y, incy)
DGEMV
Definition dgemv.f:158
dggqrf
subroutine dggqrf(n, m, p, a, lda, taua, b, ldb, taub, work, lwork, info)
DGGQRF
Definition dggqrf.f:213
ilaenv
integer function ilaenv(ispec, name, opts, n1, n2, n3, n4)
ILAENV
Definition ilaenv.f:160
dtrtrs
subroutine dtrtrs(uplo, trans, diag, n, nrhs, a, lda, b, ldb, info)
DTRTRS
Definition dtrtrs.f:144
dormqr
subroutine dormqr(side, trans, m, n, k, a, lda, tau, c, ldc, work, lwork, info)
DORMQR
Definition dormqr.f:165
dormrq
subroutine dormrq(side, trans, m, n, k, a, lda, tau, c, ldc, work, lwork, info)
DORMRQ
Definition dormrq.f:165
Here is the call graph for this function:
Here is the caller graph for this function: