LAPACK  3.10.0
LAPACK: Linear Algebra PACKage

◆ zgglse()

subroutine zgglse ( integer  M,
integer  N,
integer  P,
complex*16, dimension( lda, * )  A,
integer  LDA,
complex*16, dimension( ldb, * )  B,
integer  LDB,
complex*16, dimension( * )  C,
complex*16, dimension( * )  D,
complex*16, dimension( * )  X,
complex*16, dimension( * )  WORK,
integer  LWORK,
integer  INFO 
)

ZGGLSE solves overdetermined or underdetermined systems for OTHER matrices

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

Purpose:
 ZGGLSE solves the linear equality-constrained least squares (LSE)
 problem:

         minimize || c - A*x ||_2   subject to   B*x = d

 where A is an M-by-N matrix, B is a P-by-N matrix, c is a given
 M-vector, and d is a given P-vector. It is assumed that
 P <= N <= M+P, and

          rank(B) = P and  rank( (A) ) = N.
                               ( (B) )

 These conditions ensure that the LSE problem has a unique solution,
 which is obtained using a generalized RQ factorization of the
 matrices (B, A) given by

    B = (0 R)*Q,   A = Z*T*Q.
Parameters
[in]M
          M is INTEGER
          The number of rows of the matrix A.  M >= 0.
[in]N
          N is INTEGER
          The number of columns of the matrices A and B. N >= 0.
[in]P
          P is INTEGER
          The number of rows of the matrix B. 0 <= P <= N <= M+P.
[in,out]A
          A is COMPLEX*16 array, dimension (LDA,N)
          On entry, the M-by-N matrix A.
          On exit, the elements on and above the diagonal of the array
          contain the min(M,N)-by-N upper trapezoidal matrix T.
[in]LDA
          LDA is INTEGER
          The leading dimension of the array A. LDA >= max(1,M).
[in,out]B
          B is COMPLEX*16 array, dimension (LDB,N)
          On entry, the P-by-N matrix B.
          On exit, the upper triangle of the subarray B(1:P,N-P+1:N)
          contains the P-by-P upper triangular matrix R.
[in]LDB
          LDB is INTEGER
          The leading dimension of the array B. LDB >= max(1,P).
[in,out]C
          C is COMPLEX*16 array, dimension (M)
          On entry, C contains the right hand side vector for the
          least squares part of the LSE problem.
          On exit, the residual sum of squares for the solution
          is given by the sum of squares of elements N-P+1 to M of
          vector C.
[in,out]D
          D is COMPLEX*16 array, dimension (P)
          On entry, D contains the right hand side vector for the
          constrained equation.
          On exit, D is destroyed.
[out]X
          X is COMPLEX*16 array, dimension (N)
          On exit, X is the solution of the LSE 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,M+N+P).
          For optimum performance LWORK >= P+min(M,N)+max(M,N)*NB,
          where NB is an upper bound for the optimal blocksizes for
          ZGEQRF, CGERQF, ZUNMQR 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 B in the
                generalized RQ factorization of the pair (B, A) is
                singular, so that rank(B) < P; the least squares
                solution could not be computed.
          = 2:  the (N-P) by (N-P) part of the upper trapezoidal factor
                T associated with A in the generalized RQ factorization
                of the pair (B, A) is singular, so that
                rank( (A) ) < N; the least squares solution could not
                    ( (B) )
                be computed.
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.

Definition at line 178 of file zgglse.f.

180 *
181 * -- LAPACK driver routine --
182 * -- LAPACK is a software package provided by Univ. of Tennessee, --
183 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
184 *
185 * .. Scalar Arguments ..
186  INTEGER INFO, LDA, LDB, LWORK, M, N, P
187 * ..
188 * .. Array Arguments ..
189  COMPLEX*16 A( LDA, * ), B( LDB, * ), C( * ), D( * ),
190  $ WORK( * ), X( * )
191 * ..
192 *
193 * =====================================================================
194 *
195 * .. Parameters ..
196  COMPLEX*16 CONE
197  parameter( cone = ( 1.0d+0, 0.0d+0 ) )
198 * ..
199 * .. Local Scalars ..
200  LOGICAL LQUERY
201  INTEGER LOPT, LWKMIN, LWKOPT, MN, NB, NB1, NB2, NB3,
202  $ NB4, NR
203 * ..
204 * .. External Subroutines ..
205  EXTERNAL xerbla, zaxpy, zcopy, zgemv, zggrqf, ztrmv,
206  $ ztrtrs, zunmqr, zunmrq
207 * ..
208 * .. External Functions ..
209  INTEGER ILAENV
210  EXTERNAL ilaenv
211 * ..
212 * .. Intrinsic Functions ..
213  INTRINSIC int, max, min
214 * ..
215 * .. Executable Statements ..
216 *
217 * Test the input parameters
218 *
219  info = 0
220  mn = min( m, n )
221  lquery = ( lwork.EQ.-1 )
222  IF( m.LT.0 ) THEN
223  info = -1
224  ELSE IF( n.LT.0 ) THEN
225  info = -2
226  ELSE IF( p.LT.0 .OR. p.GT.n .OR. p.LT.n-m ) THEN
227  info = -3
228  ELSE IF( lda.LT.max( 1, m ) ) THEN
229  info = -5
230  ELSE IF( ldb.LT.max( 1, p ) ) THEN
231  info = -7
232  END IF
233 *
234 * Calculate workspace
235 *
236  IF( info.EQ.0) THEN
237  IF( n.EQ.0 ) THEN
238  lwkmin = 1
239  lwkopt = 1
240  ELSE
241  nb1 = ilaenv( 1, 'ZGEQRF', ' ', m, n, -1, -1 )
242  nb2 = ilaenv( 1, 'ZGERQF', ' ', m, n, -1, -1 )
243  nb3 = ilaenv( 1, 'ZUNMQR', ' ', m, n, p, -1 )
244  nb4 = ilaenv( 1, 'ZUNMRQ', ' ', m, n, p, -1 )
245  nb = max( nb1, nb2, nb3, nb4 )
246  lwkmin = m + n + p
247  lwkopt = p + mn + max( m, n )*nb
248  END IF
249  work( 1 ) = lwkopt
250 *
251  IF( lwork.LT.lwkmin .AND. .NOT.lquery ) THEN
252  info = -12
253  END IF
254  END IF
255 *
256  IF( info.NE.0 ) THEN
257  CALL xerbla( 'ZGGLSE', -info )
258  RETURN
259  ELSE IF( lquery ) THEN
260  RETURN
261  END IF
262 *
263 * Quick return if possible
264 *
265  IF( n.EQ.0 )
266  $ RETURN
267 *
268 * Compute the GRQ factorization of matrices B and A:
269 *
270 * B*Q**H = ( 0 T12 ) P Z**H*A*Q**H = ( R11 R12 ) N-P
271 * N-P P ( 0 R22 ) M+P-N
272 * N-P P
273 *
274 * where T12 and R11 are upper triangular, and Q and Z are
275 * unitary.
276 *
277  CALL zggrqf( p, m, n, b, ldb, work, a, lda, work( p+1 ),
278  $ work( p+mn+1 ), lwork-p-mn, info )
279  lopt = dble( work( p+mn+1 ) )
280 *
281 * Update c = Z**H *c = ( c1 ) N-P
282 * ( c2 ) M+P-N
283 *
284  CALL zunmqr( 'Left', 'Conjugate Transpose', m, 1, mn, a, lda,
285  $ work( p+1 ), c, max( 1, m ), work( p+mn+1 ),
286  $ lwork-p-mn, info )
287  lopt = max( lopt, int( work( p+mn+1 ) ) )
288 *
289 * Solve T12*x2 = d for x2
290 *
291  IF( p.GT.0 ) THEN
292  CALL ztrtrs( 'Upper', 'No transpose', 'Non-unit', p, 1,
293  $ b( 1, n-p+1 ), ldb, d, p, info )
294 *
295  IF( info.GT.0 ) THEN
296  info = 1
297  RETURN
298  END IF
299 *
300 * Put the solution in X
301 *
302  CALL zcopy( p, d, 1, x( n-p+1 ), 1 )
303 *
304 * Update c1
305 *
306  CALL zgemv( 'No transpose', n-p, p, -cone, a( 1, n-p+1 ), lda,
307  $ d, 1, cone, c, 1 )
308  END IF
309 *
310 * Solve R11*x1 = c1 for x1
311 *
312  IF( n.GT.p ) THEN
313  CALL ztrtrs( 'Upper', 'No transpose', 'Non-unit', n-p, 1,
314  $ a, lda, c, n-p, info )
315 *
316  IF( info.GT.0 ) THEN
317  info = 2
318  RETURN
319  END IF
320 *
321 * Put the solutions in X
322 *
323  CALL zcopy( n-p, c, 1, x, 1 )
324  END IF
325 *
326 * Compute the residual vector:
327 *
328  IF( m.LT.n ) THEN
329  nr = m + p - n
330  IF( nr.GT.0 )
331  $ CALL zgemv( 'No transpose', nr, n-m, -cone, a( n-p+1, m+1 ),
332  $ lda, d( nr+1 ), 1, cone, c( n-p+1 ), 1 )
333  ELSE
334  nr = p
335  END IF
336  IF( nr.GT.0 ) THEN
337  CALL ztrmv( 'Upper', 'No transpose', 'Non unit', nr,
338  $ a( n-p+1, n-p+1 ), lda, d, 1 )
339  CALL zaxpy( nr, -cone, d, 1, c( n-p+1 ), 1 )
340  END IF
341 *
342 * Backward transformation x = Q**H*x
343 *
344  CALL zunmrq( 'Left', 'Conjugate Transpose', n, 1, p, b, ldb,
345  $ work( 1 ), x, n, work( p+mn+1 ), lwork-p-mn, info )
346  work( 1 ) = p + mn + max( lopt, int( work( p+mn+1 ) ) )
347 *
348  RETURN
349 *
350 * End of ZGGLSE
351 *
integer function ilaenv(ISPEC, NAME, OPTS, N1, N2, N3, N4)
ILAENV
Definition: ilaenv.f:162
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:60
subroutine zaxpy(N, ZA, ZX, INCX, ZY, INCY)
ZAXPY
Definition: zaxpy.f:88
subroutine zcopy(N, ZX, INCX, ZY, INCY)
ZCOPY
Definition: zcopy.f:81
subroutine ztrmv(UPLO, TRANS, DIAG, N, A, LDA, X, INCX)
ZTRMV
Definition: ztrmv.f:147
subroutine zgemv(TRANS, M, N, ALPHA, A, LDA, X, INCX, BETA, Y, INCY)
ZGEMV
Definition: zgemv.f:158
subroutine ztrtrs(UPLO, TRANS, DIAG, N, NRHS, A, LDA, B, LDB, INFO)
ZTRTRS
Definition: ztrtrs.f:140
subroutine zggrqf(M, P, N, A, LDA, TAUA, B, LDB, TAUB, WORK, LWORK, INFO)
ZGGRQF
Definition: zggrqf.f:214
subroutine zunmrq(SIDE, TRANS, M, N, K, A, LDA, TAU, C, LDC, WORK, LWORK, INFO)
ZUNMRQ
Definition: zunmrq.f:167
subroutine zunmqr(SIDE, TRANS, M, N, K, A, LDA, TAU, C, LDC, WORK, LWORK, INFO)
ZUNMQR
Definition: zunmqr.f:167
Here is the call graph for this function:
Here is the caller graph for this function: