LAPACK  3.6.1
LAPACK: Linear Algebra PACKage
subroutine cggglm ( integer  N,
integer  M,
integer  P,
complex, dimension( lda, * )  A,
integer  LDA,
complex, dimension( ldb, * )  B,
integer  LDB,
complex, dimension( * )  D,
complex, dimension( * )  X,
complex, dimension( * )  Y,
complex, dimension( * )  WORK,
integer  LWORK,
integer  INFO 
)

CGGGLM

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

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

          On exit, X and Y are the solutions of the GLM problem.
[out]WORK
          WORK is COMPLEX 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
          CGEQRF, CGERQF, CUNMQR and CUNMRQ.

          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.
Date
November 2015

Definition at line 187 of file cggglm.f.

187 *
188 * -- LAPACK driver routine (version 3.6.0) --
189 * -- LAPACK is a software package provided by Univ. of Tennessee, --
190 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
191 * November 2015
192 *
193 * .. Scalar Arguments ..
194  INTEGER info, lda, ldb, lwork, m, n, p
195 * ..
196 * .. Array Arguments ..
197  COMPLEX a( lda, * ), b( ldb, * ), d( * ), work( * ),
198  $ x( * ), y( * )
199 * ..
200 *
201 * ===================================================================
202 *
203 * .. Parameters ..
204  COMPLEX czero, cone
205  parameter ( czero = ( 0.0e+0, 0.0e+0 ),
206  $ cone = ( 1.0e+0, 0.0e+0 ) )
207 * ..
208 * .. Local Scalars ..
209  LOGICAL lquery
210  INTEGER i, lopt, lwkmin, lwkopt, nb, nb1, nb2, nb3,
211  $ nb4, np
212 * ..
213 * .. External Subroutines ..
214  EXTERNAL ccopy, cgemv, cggqrf, ctrtrs, cunmqr, cunmrq,
215  $ xerbla
216 * ..
217 * .. External Functions ..
218  INTEGER ilaenv
219  EXTERNAL ilaenv
220 * ..
221 * .. Intrinsic Functions ..
222  INTRINSIC int, max, min
223 * ..
224 * .. Executable Statements ..
225 *
226 * Test the input parameters
227 *
228  info = 0
229  np = min( n, p )
230  lquery = ( lwork.EQ.-1 )
231  IF( n.LT.0 ) THEN
232  info = -1
233  ELSE IF( m.LT.0 .OR. m.GT.n ) THEN
234  info = -2
235  ELSE IF( p.LT.0 .OR. p.LT.n-m ) THEN
236  info = -3
237  ELSE IF( lda.LT.max( 1, n ) ) THEN
238  info = -5
239  ELSE IF( ldb.LT.max( 1, n ) ) THEN
240  info = -7
241  END IF
242 *
243 * Calculate workspace
244 *
245  IF( info.EQ.0) THEN
246  IF( n.EQ.0 ) THEN
247  lwkmin = 1
248  lwkopt = 1
249  ELSE
250  nb1 = ilaenv( 1, 'CGEQRF', ' ', n, m, -1, -1 )
251  nb2 = ilaenv( 1, 'CGERQF', ' ', n, m, -1, -1 )
252  nb3 = ilaenv( 1, 'CUNMQR', ' ', n, m, p, -1 )
253  nb4 = ilaenv( 1, 'CUNMRQ', ' ', n, m, p, -1 )
254  nb = max( nb1, nb2, nb3, nb4 )
255  lwkmin = m + n + p
256  lwkopt = m + np + max( n, p )*nb
257  END IF
258  work( 1 ) = lwkopt
259 *
260  IF( lwork.LT.lwkmin .AND. .NOT.lquery ) THEN
261  info = -12
262  END IF
263  END IF
264 *
265  IF( info.NE.0 ) THEN
266  CALL xerbla( 'CGGGLM', -info )
267  RETURN
268  ELSE IF( lquery ) THEN
269  RETURN
270  END IF
271 *
272 * Quick return if possible
273 *
274  IF( n.EQ.0 )
275  $ RETURN
276 *
277 * Compute the GQR factorization of matrices A and B:
278 *
279 * Q**H*A = ( R11 ) M, Q**H*B*Z**H = ( T11 T12 ) M
280 * ( 0 ) N-M ( 0 T22 ) N-M
281 * M M+P-N N-M
282 *
283 * where R11 and T22 are upper triangular, and Q and Z are
284 * unitary.
285 *
286  CALL cggqrf( n, m, p, a, lda, work, b, ldb, work( m+1 ),
287  $ work( m+np+1 ), lwork-m-np, info )
288  lopt = work( m+np+1 )
289 *
290 * Update left-hand-side vector d = Q**H*d = ( d1 ) M
291 * ( d2 ) N-M
292 *
293  CALL cunmqr( 'Left', 'Conjugate transpose', n, 1, m, a, lda, work,
294  $ d, max( 1, n ), work( m+np+1 ), lwork-m-np, info )
295  lopt = max( lopt, int( work( m+np+1 ) ) )
296 *
297 * Solve T22*y2 = d2 for y2
298 *
299  IF( n.GT.m ) THEN
300  CALL ctrtrs( 'Upper', 'No transpose', 'Non unit', n-m, 1,
301  $ b( m+1, m+p-n+1 ), ldb, d( m+1 ), n-m, info )
302 *
303  IF( info.GT.0 ) THEN
304  info = 1
305  RETURN
306  END IF
307 *
308  CALL ccopy( n-m, d( m+1 ), 1, y( m+p-n+1 ), 1 )
309  END IF
310 *
311 * Set y1 = 0
312 *
313  DO 10 i = 1, m + p - n
314  y( i ) = czero
315  10 CONTINUE
316 *
317 * Update d1 = d1 - T12*y2
318 *
319  CALL cgemv( 'No transpose', m, n-m, -cone, b( 1, m+p-n+1 ), ldb,
320  $ y( m+p-n+1 ), 1, cone, d, 1 )
321 *
322 * Solve triangular system: R11*x = d1
323 *
324  IF( m.GT.0 ) THEN
325  CALL ctrtrs( 'Upper', 'No Transpose', 'Non unit', m, 1, a, lda,
326  $ d, m, info )
327 *
328  IF( info.GT.0 ) THEN
329  info = 2
330  RETURN
331  END IF
332 *
333 * Copy D to X
334 *
335  CALL ccopy( m, d, 1, x, 1 )
336  END IF
337 *
338 * Backward transformation y = Z**H *y
339 *
340  CALL cunmrq( 'Left', 'Conjugate transpose', p, 1, np,
341  $ b( max( 1, n-p+1 ), 1 ), ldb, work( m+1 ), y,
342  $ max( 1, p ), work( m+np+1 ), lwork-m-np, info )
343  work( 1 ) = m + np + max( lopt, int( work( m+np+1 ) ) )
344 *
345  RETURN
346 *
347 * End of CGGGLM
348 *
subroutine cggqrf(N, M, P, A, LDA, TAUA, B, LDB, TAUB, WORK, LWORK, INFO)
CGGQRF
Definition: cggqrf.f:217
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
subroutine cunmrq(SIDE, TRANS, M, N, K, A, LDA, TAU, C, LDC, WORK, LWORK, INFO)
CUNMRQ
Definition: cunmrq.f:170
subroutine cunmqr(SIDE, TRANS, M, N, K, A, LDA, TAU, C, LDC, WORK, LWORK, INFO)
CUNMQR
Definition: cunmqr.f:170
subroutine cgemv(TRANS, M, N, ALPHA, A, LDA, X, INCX, BETA, Y, INCY)
CGEMV
Definition: cgemv.f:160
subroutine ctrtrs(UPLO, TRANS, DIAG, N, NRHS, A, LDA, B, LDB, INFO)
CTRTRS
Definition: ctrtrs.f:142
integer function ilaenv(ISPEC, NAME, OPTS, N1, N2, N3, N4)
Definition: tstiee.f:83
subroutine ccopy(N, CX, INCX, CY, INCY)
CCOPY
Definition: ccopy.f:52

Here is the call graph for this function:

Here is the caller graph for this function: