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

◆ zggglm()

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

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