LAPACK  3.10.1
LAPACK: Linear Algebra PACKage

◆ 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.
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
                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 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 183 of file dggglm.f.

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 = 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 *
integer function ilaenv(ISPEC, NAME, OPTS, N1, N2, N3, N4)
ILAENV
Definition: ilaenv.f:162
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:60
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:156
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
subroutine dtrtrs(UPLO, TRANS, DIAG, N, NRHS, A, LDA, B, LDB, INFO)
DTRTRS
Definition: dtrtrs.f:140
subroutine dggqrf(N, M, P, A, LDA, TAUA, B, LDB, TAUB, WORK, LWORK, INFO)
DGGQRF
Definition: dggqrf.f:215
Here is the call graph for this function:
Here is the caller graph for this function: