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.
!>
!> Callers of this subroutine should note that the singularity/rank-deficiency checks
!> implemented in this subroutine are rudimentary. The STRTRS 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 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 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 sggglm.f.

*
*  -- LAPACK driver routine --
*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
*
*     .. Scalar Arguments ..
      INTEGER            INFO, LDA, LDB, LWORK, M, N, P
*     ..
*     .. Array Arguments ..
      REAL               A( LDA, * ), B( LDB, * ), D( * ), WORK( * ),
     $                   X( * ), Y( * )
*     ..
*
*  ===================================================================
*
*     .. Parameters ..
      REAL               ZERO, ONE
      parameter( zero = 0.0e+0, one = 1.0e+0 )
*     ..
*     .. Local Scalars ..
      LOGICAL            LQUERY
      INTEGER            I, LOPT, LWKMIN, LWKOPT, NB, NB1, NB2, NB3,
     $                   NB4, NP
*     ..
*     .. External Subroutines ..
      EXTERNAL           scopy, sgemv, sggqrf, sormqr, sormrq,
     $                   strtrs,
     $                   xerbla
*     ..
*     .. External Functions ..
      INTEGER            ILAENV
      REAL               SROUNDUP_LWORK
      EXTERNAL           ilaenv, sroundup_lwork
*     ..
*     .. Intrinsic Functions ..
      INTRINSIC          int, max, min
*     ..
*     .. Executable Statements ..
*
*     Test the input parameters
*
      info = 0
      np = min( n, p )
      lquery = ( lwork.EQ.-1 )
      IF( n.LT.0 ) THEN
         info = -1
      ELSE IF( m.LT.0 .OR. m.GT.n ) THEN
         info = -2
      ELSE IF( p.LT.0 .OR. p.LT.n-m ) THEN
         info = -3
      ELSE IF( lda.LT.max( 1, n ) ) THEN
         info = -5
      ELSE IF( ldb.LT.max( 1, n ) ) THEN
         info = -7
      END IF
*
*     Calculate workspace
*
      IF( info.EQ.0) THEN
         IF( n.EQ.0 ) THEN
            lwkmin = 1
            lwkopt = 1
         ELSE
            nb1 = ilaenv( 1, 'SGEQRF', ' ', n, m, -1, -1 )
            nb2 = ilaenv( 1, 'SGERQF', ' ', n, m, -1, -1 )
            nb3 = ilaenv( 1, 'SORMQR', ' ', n, m, p, -1 )
            nb4 = ilaenv( 1, 'SORMRQ', ' ', n, m, p, -1 )
            nb = max( nb1, nb2, nb3, nb4 )
            lwkmin = m + n + p
            lwkopt = m + np + max( n, p )*nb
         END IF
         work( 1 ) = sroundup_lwork(lwkopt)
*
         IF( lwork.LT.lwkmin .AND. .NOT.lquery ) THEN
            info = -12
         END IF
      END IF
*
      IF( info.NE.0 ) THEN
         CALL xerbla( 'SGGGLM', -info )
         RETURN
      ELSE IF( lquery ) THEN
         RETURN
      END IF
*
*     Quick return if possible
*
      IF( n.EQ.0 ) THEN
         DO i = 1, m
            x(i) = zero
         END DO
         DO i = 1, p
            y(i) = zero
         END DO
         RETURN
      END IF
*
*     Compute the GQR factorization of matrices A and B:
*
*          Q**T*A = ( R11 ) M,    Q**T*B*Z**T = ( T11   T12 ) M
*                   (  0  ) N-M                 (  0    T22 ) N-M
*                      M                         M+P-N  N-M
*
*     where R11 and T22 are upper triangular, and Q and Z are
*     orthogonal.
*
      CALL sggqrf( n, m, p, a, lda, work, b, ldb, work( m+1 ),
     $             work( m+np+1 ), lwork-m-np, info )
      lopt = int( work( m+np+1 ) )
*
*     Update left-hand-side vector d = Q**T*d = ( d1 ) M
*                                               ( d2 ) N-M
*
      CALL sormqr( 'Left', 'Transpose', n, 1, m, a, lda, work, d,
     $             max( 1, n ), work( m+np+1 ), lwork-m-np, info )
      lopt = max( lopt, int( work( m+np+1 ) ) )
*
*     Solve T22*y2 = d2 for y2
*
      IF( n.GT.m ) THEN
         CALL strtrs( 'Upper', 'No transpose', 'Non unit', n-m, 1,
     $                b( m+1, m+p-n+1 ), ldb, d( m+1 ), n-m, info )
*
         IF( info.GT.0 ) THEN
            info = 1
            RETURN
         END IF
*
         CALL scopy( n-m, d( m+1 ), 1, y( m+p-n+1 ), 1 )
      END IF
*
*     Set y1 = 0
*
      DO 10 i = 1, m + p - n
         y( i ) = zero
   10 CONTINUE
*
*     Update d1 = d1 - T12*y2
*
      CALL sgemv( 'No transpose', m, n-m, -one, b( 1, m+p-n+1 ), ldb,
     $            y( m+p-n+1 ), 1, one, d, 1 )
*
*     Solve triangular system: R11*x = d1
*
      IF( m.GT.0 ) THEN
         CALL strtrs( 'Upper', 'No Transpose', 'Non unit', m, 1, a,
     $                lda,
     $                d, m, info )
*
         IF( info.GT.0 ) THEN
            info = 2
            RETURN
         END IF
*
*        Copy D to X
*
         CALL scopy( m, d, 1, x, 1 )
      END IF
*
*     Backward transformation y = Z**T *y
*
      CALL sormrq( 'Left', 'Transpose', p, 1, np,
     $             b( max( 1, n-p+1 ), 1 ), ldb, work( m+1 ), y,
     $             max( 1, p ), work( m+np+1 ), lwork-m-np, info )
      work( 1 ) = real( m + np + max( lopt, int( work( m+np+1 ) ) ) )
*
      RETURN
*
*     End of SGGGLM
*

Here is the call graph for this function:

Here is the caller graph for this function: