LAPACK  3.10.0
LAPACK: Linear Algebra PACKage

◆ sggglm()

subroutine sggglm ( integer  N,
integer  M,
integer  P,
real, dimension( lda, * )  A,
integer  LDA,
real, dimension( ldb, * )  B,
integer  LDB,
real, dimension( * )  D,
real, dimension( * )  X,
real, dimension( * )  Y,
real, dimension( * )  WORK,
integer  LWORK,
integer  INFO 
)

SGGGLM

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

Purpose:
 SGGGLM 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 REAL 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 REAL 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 REAL array, dimension (N)
          On entry, D is the left hand side of the GLM equation.
          On exit, D is destroyed.
[out]X
          X is REAL array, dimension (M)
[out]Y
          Y is REAL array, dimension (P)

          On exit, X and Y are the solutions of the GLM problem.
[out]WORK
          WORK is REAL 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
          SGEQRF, SGERQF, SORMQR 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 sggglm.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  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  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, 'SGEQRF', ' ', n, m, -1, -1 )
247  nb2 = ilaenv( 1, 'SGERQF', ' ', n, m, -1, -1 )
248  nb3 = ilaenv( 1, 'SORMQR', ' ', n, m, p, -1 )
249  nb4 = ilaenv( 1, 'SORMRQ', ' ', 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( 'SGGGLM', -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 sggqrf( 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 sormqr( '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 strtrs( '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 scopy( 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 sgemv( '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 strtrs( '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 scopy( m, d, 1, x, 1 )
339  END IF
340 *
341 * Backward transformation y = Z**T *y
342 *
343  CALL sormrq( '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 SGGGLM
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 sormrq(SIDE, TRANS, M, N, K, A, LDA, TAU, C, LDC, WORK, LWORK, INFO)
SORMRQ
Definition: sormrq.f:168
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 sggqrf(N, M, P, A, LDA, TAUA, B, LDB, TAUB, WORK, LWORK, INFO)
SGGQRF
Definition: sggqrf.f:215
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:156
Here is the call graph for this function:
Here is the caller graph for this function: