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

◆ cggglm()

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.

Definition at line 183 of file cggglm.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 COMPLEX A( LDA, * ), B( LDB, * ), D( * ), WORK( * ),
195 $ X( * ), Y( * )
196* ..
197*
198* ===================================================================
199*
200* .. Parameters ..
201 COMPLEX CZERO, CONE
202 parameter( czero = ( 0.0e+0, 0.0e+0 ),
203 $ cone = ( 1.0e+0, 0.0e+0 ) )
204* ..
205* .. Local Scalars ..
206 LOGICAL LQUERY
207 INTEGER I, LOPT, LWKMIN, LWKOPT, NB, NB1, NB2, NB3,
208 $ NB4, NP
209* ..
210* .. External Subroutines ..
211 EXTERNAL ccopy, cgemv, cggqrf, ctrtrs, cunmqr, cunmrq,
212 $ xerbla
213* ..
214* .. External Functions ..
215 INTEGER ILAENV
216 REAL SROUNDUP_LWORK
217 EXTERNAL ilaenv, sroundup_lwork
218* ..
219* .. Intrinsic Functions ..
220 INTRINSIC int, max, min
221* ..
222* .. Executable Statements ..
223*
224* Test the input parameters
225*
226 info = 0
227 np = min( n, p )
228 lquery = ( lwork.EQ.-1 )
229 IF( n.LT.0 ) THEN
230 info = -1
231 ELSE IF( m.LT.0 .OR. m.GT.n ) THEN
232 info = -2
233 ELSE IF( p.LT.0 .OR. p.LT.n-m ) THEN
234 info = -3
235 ELSE IF( lda.LT.max( 1, n ) ) THEN
236 info = -5
237 ELSE IF( ldb.LT.max( 1, n ) ) THEN
238 info = -7
239 END IF
240*
241* Calculate workspace
242*
243 IF( info.EQ.0) THEN
244 IF( n.EQ.0 ) THEN
245 lwkmin = 1
246 lwkopt = 1
247 ELSE
248 nb1 = ilaenv( 1, 'CGEQRF', ' ', n, m, -1, -1 )
249 nb2 = ilaenv( 1, 'CGERQF', ' ', n, m, -1, -1 )
250 nb3 = ilaenv( 1, 'CUNMQR', ' ', n, m, p, -1 )
251 nb4 = ilaenv( 1, 'CUNMRQ', ' ', n, m, p, -1 )
252 nb = max( nb1, nb2, nb3, nb4 )
253 lwkmin = m + n + p
254 lwkopt = m + np + max( n, p )*nb
255 END IF
256 work( 1 ) = sroundup_lwork(lwkopt)
257*
258 IF( lwork.LT.lwkmin .AND. .NOT.lquery ) THEN
259 info = -12
260 END IF
261 END IF
262*
263 IF( info.NE.0 ) THEN
264 CALL xerbla( 'CGGGLM', -info )
265 RETURN
266 ELSE IF( lquery ) THEN
267 RETURN
268 END IF
269*
270* Quick return if possible
271*
272 IF( n.EQ.0 ) THEN
273 DO i = 1, m
274 x(i) = czero
275 END DO
276 DO i = 1, p
277 y(i) = czero
278 END DO
279 RETURN
280 END IF
281*
282* Compute the GQR factorization of matrices A and B:
283*
284* Q**H*A = ( R11 ) M, Q**H*B*Z**H = ( T11 T12 ) M
285* ( 0 ) N-M ( 0 T22 ) N-M
286* M M+P-N N-M
287*
288* where R11 and T22 are upper triangular, and Q and Z are
289* unitary.
290*
291 CALL cggqrf( n, m, p, a, lda, work, b, ldb, work( m+1 ),
292 $ work( m+np+1 ), lwork-m-np, info )
293 lopt = int( work( m+np+1 ) )
294*
295* Update left-hand-side vector d = Q**H*d = ( d1 ) M
296* ( d2 ) N-M
297*
298 CALL cunmqr( 'Left', 'Conjugate transpose', n, 1, m, a, lda, work,
299 $ d, max( 1, n ), work( m+np+1 ), lwork-m-np, info )
300 lopt = max( lopt, int( work( m+np+1 ) ) )
301*
302* Solve T22*y2 = d2 for y2
303*
304 IF( n.GT.m ) THEN
305 CALL ctrtrs( 'Upper', 'No transpose', 'Non unit', n-m, 1,
306 $ b( m+1, m+p-n+1 ), ldb, d( m+1 ), n-m, info )
307*
308 IF( info.GT.0 ) THEN
309 info = 1
310 RETURN
311 END IF
312*
313 CALL ccopy( n-m, d( m+1 ), 1, y( m+p-n+1 ), 1 )
314 END IF
315*
316* Set y1 = 0
317*
318 DO 10 i = 1, m + p - n
319 y( i ) = czero
320 10 CONTINUE
321*
322* Update d1 = d1 - T12*y2
323*
324 CALL cgemv( 'No transpose', m, n-m, -cone, b( 1, m+p-n+1 ), ldb,
325 $ y( m+p-n+1 ), 1, cone, d, 1 )
326*
327* Solve triangular system: R11*x = d1
328*
329 IF( m.GT.0 ) THEN
330 CALL ctrtrs( 'Upper', 'No Transpose', 'Non unit', m, 1, a, lda,
331 $ d, m, info )
332*
333 IF( info.GT.0 ) THEN
334 info = 2
335 RETURN
336 END IF
337*
338* Copy D to X
339*
340 CALL ccopy( m, d, 1, x, 1 )
341 END IF
342*
343* Backward transformation y = Z**H *y
344*
345 CALL cunmrq( 'Left', 'Conjugate transpose', p, 1, np,
346 $ b( max( 1, n-p+1 ), 1 ), ldb, work( m+1 ), y,
347 $ max( 1, p ), work( m+np+1 ), lwork-m-np, info )
348 work( 1 ) = m + np + max( lopt, int( work( m+np+1 ) ) )
349*
350 RETURN
351*
352* End of CGGGLM
353*
subroutine xerbla(srname, info)
Definition cblat2.f:3285
subroutine ccopy(n, cx, incx, cy, incy)
CCOPY
Definition ccopy.f:81
subroutine cgemv(trans, m, n, alpha, a, lda, x, incx, beta, y, incy)
CGEMV
Definition cgemv.f:160
subroutine cggqrf(n, m, p, a, lda, taua, b, ldb, taub, work, lwork, info)
CGGQRF
Definition cggqrf.f:215
integer function ilaenv(ispec, name, opts, n1, n2, n3, n4)
ILAENV
Definition ilaenv.f:162
real function sroundup_lwork(lwork)
SROUNDUP_LWORK
subroutine ctrtrs(uplo, trans, diag, n, nrhs, a, lda, b, ldb, info)
CTRTRS
Definition ctrtrs.f:140
subroutine cunmqr(side, trans, m, n, k, a, lda, tau, c, ldc, work, lwork, info)
CUNMQR
Definition cunmqr.f:168
subroutine cunmrq(side, trans, m, n, k, a, lda, tau, c, ldc, work, lwork, info)
CUNMRQ
Definition cunmrq.f:168
Here is the call graph for this function:
Here is the caller graph for this function: