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

ZGGGLM

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

Purpose:
 ZGGGLM 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*16 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*16 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*16 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*16 array, dimension (M)
[out]Y
          Y is COMPLEX*16 array, dimension (P)

          On exit, X and Y are the solutions of the GLM problem.
[out]WORK
          WORK is COMPLEX*16 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
          ZGEQRF, ZGERQF, ZUNMQR and ZUNMRQ.

          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 zggglm.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*16 a( lda, * ), b( ldb, * ), d( * ), work( * ),
198  $ x( * ), y( * )
199 * ..
200 *
201 * ===================================================================
202 *
203 * .. Parameters ..
204  COMPLEX*16 czero, cone
205  parameter ( czero = ( 0.0d+0, 0.0d+0 ),
206  $ cone = ( 1.0d+0, 0.0d+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 xerbla, zcopy, zgemv, zggqrf, ztrtrs, zunmqr,
215  $ zunmrq
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, 'ZGEQRF', ' ', n, m, -1, -1 )
251  nb2 = ilaenv( 1, 'ZGERQF', ' ', n, m, -1, -1 )
252  nb3 = ilaenv( 1, 'ZUNMQR', ' ', n, m, p, -1 )
253  nb4 = ilaenv( 1, 'ZUNMRQ', ' ', 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( 'ZGGGLM', -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 zggqrf( 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 zunmqr( '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 ztrtrs( '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 zcopy( 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 zgemv( '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 ztrtrs( '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 zcopy( m, d, 1, x, 1 )
336  END IF
337 *
338 * Backward transformation y = Z**H *y
339 *
340  CALL zunmrq( '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 ZGGGLM
348 *
subroutine ztrtrs(UPLO, TRANS, DIAG, N, NRHS, A, LDA, B, LDB, INFO)
ZTRTRS
Definition: ztrtrs.f:142
subroutine zcopy(N, ZX, INCX, ZY, INCY)
ZCOPY
Definition: zcopy.f:52
subroutine zgemv(TRANS, M, N, ALPHA, A, LDA, X, INCX, BETA, Y, INCY)
ZGEMV
Definition: zgemv.f:160
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
subroutine zunmqr(SIDE, TRANS, M, N, K, A, LDA, TAU, C, LDC, WORK, LWORK, INFO)
ZUNMQR
Definition: zunmqr.f:169
subroutine zggqrf(N, M, P, A, LDA, TAUA, B, LDB, TAUB, WORK, LWORK, INFO)
ZGGQRF
Definition: zggqrf.f:217
integer function ilaenv(ISPEC, NAME, OPTS, N1, N2, N3, N4)
Definition: tstiee.f:83
subroutine zunmrq(SIDE, TRANS, M, N, K, A, LDA, TAU, C, LDC, WORK, LWORK, INFO)
ZUNMRQ
Definition: zunmrq.f:169

Here is the call graph for this function:

Here is the caller graph for this function: