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

◆ cgglse()

subroutine cgglse ( integer  m,
integer  n,
integer  p,
complex, dimension( lda, * )  a,
integer  lda,
complex, dimension( ldb, * )  b,
integer  ldb,
complex, dimension( * )  c,
complex, dimension( * )  d,
complex, dimension( * )  x,
complex, dimension( * )  work,
integer  lwork,
integer  info 
)

CGGLSE solves overdetermined or underdetermined systems for OTHER matrices

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

Purpose:
 CGGLSE 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 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 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 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 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 array, dimension (N)
          On exit, X is the solution of the LSE 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,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
          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 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 cgglse.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 A( LDA, * ), B( LDB, * ), C( * ), D( * ),
190 $ WORK( * ), X( * )
191* ..
192*
193* =====================================================================
194*
195* .. Parameters ..
196 COMPLEX CONE
197 parameter( cone = ( 1.0e+0, 0.0e+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 caxpy, ccopy, cgemv, cggrqf, ctrmv, ctrtrs,
207* ..
208* .. External Functions ..
209 INTEGER ILAENV
210 REAL SROUNDUP_LWORK
211 EXTERNAL ilaenv, sroundup_lwork
212* ..
213* .. Intrinsic Functions ..
214 INTRINSIC int, max, min
215* ..
216* .. Executable Statements ..
217*
218* Test the input parameters
219*
220 info = 0
221 mn = min( m, n )
222 lquery = ( lwork.EQ.-1 )
223 IF( m.LT.0 ) THEN
224 info = -1
225 ELSE IF( n.LT.0 ) THEN
226 info = -2
227 ELSE IF( p.LT.0 .OR. p.GT.n .OR. p.LT.n-m ) THEN
228 info = -3
229 ELSE IF( lda.LT.max( 1, m ) ) THEN
230 info = -5
231 ELSE IF( ldb.LT.max( 1, p ) ) THEN
232 info = -7
233 END IF
234*
235* Calculate workspace
236*
237 IF( info.EQ.0) THEN
238 IF( n.EQ.0 ) THEN
239 lwkmin = 1
240 lwkopt = 1
241 ELSE
242 nb1 = ilaenv( 1, 'CGEQRF', ' ', m, n, -1, -1 )
243 nb2 = ilaenv( 1, 'CGERQF', ' ', m, n, -1, -1 )
244 nb3 = ilaenv( 1, 'CUNMQR', ' ', m, n, p, -1 )
245 nb4 = ilaenv( 1, 'CUNMRQ', ' ', m, n, p, -1 )
246 nb = max( nb1, nb2, nb3, nb4 )
247 lwkmin = m + n + p
248 lwkopt = p + mn + max( m, n )*nb
249 END IF
250 work( 1 ) = sroundup_lwork(lwkopt)
251*
252 IF( lwork.LT.lwkmin .AND. .NOT.lquery ) THEN
253 info = -12
254 END IF
255 END IF
256*
257 IF( info.NE.0 ) THEN
258 CALL xerbla( 'CGGLSE', -info )
259 RETURN
260 ELSE IF( lquery ) THEN
261 RETURN
262 END IF
263*
264* Quick return if possible
265*
266 IF( n.EQ.0 )
267 $ RETURN
268*
269* Compute the GRQ factorization of matrices B and A:
270*
271* B*Q**H = ( 0 T12 ) P Z**H*A*Q**H = ( R11 R12 ) N-P
272* N-P P ( 0 R22 ) M+P-N
273* N-P P
274*
275* where T12 and R11 are upper triangular, and Q and Z are
276* unitary.
277*
278 CALL cggrqf( p, m, n, b, ldb, work, a, lda, work( p+1 ),
279 $ work( p+mn+1 ), lwork-p-mn, info )
280 lopt = int( work( p+mn+1 ) )
281*
282* Update c = Z**H *c = ( c1 ) N-P
283* ( c2 ) M+P-N
284*
285 CALL cunmqr( 'Left', 'Conjugate Transpose', m, 1, mn, a, lda,
286 $ work( p+1 ), c, max( 1, m ), work( p+mn+1 ),
287 $ lwork-p-mn, info )
288 lopt = max( lopt, int( work( p+mn+1 ) ) )
289*
290* Solve T12*x2 = d for x2
291*
292 IF( p.GT.0 ) THEN
293 CALL ctrtrs( 'Upper', 'No transpose', 'Non-unit', p, 1,
294 $ b( 1, n-p+1 ), ldb, d, p, info )
295*
296 IF( info.GT.0 ) THEN
297 info = 1
298 RETURN
299 END IF
300*
301* Put the solution in X
302*
303 CALL ccopy( p, d, 1, x( n-p+1 ), 1 )
304*
305* Update c1
306*
307 CALL cgemv( 'No transpose', n-p, p, -cone, a( 1, n-p+1 ), lda,
308 $ d, 1, cone, c, 1 )
309 END IF
310*
311* Solve R11*x1 = c1 for x1
312*
313 IF( n.GT.p ) THEN
314 CALL ctrtrs( 'Upper', 'No transpose', 'Non-unit', n-p, 1,
315 $ a, lda, c, n-p, info )
316*
317 IF( info.GT.0 ) THEN
318 info = 2
319 RETURN
320 END IF
321*
322* Put the solutions in X
323*
324 CALL ccopy( n-p, c, 1, x, 1 )
325 END IF
326*
327* Compute the residual vector:
328*
329 IF( m.LT.n ) THEN
330 nr = m + p - n
331 IF( nr.GT.0 )
332 $ CALL cgemv( 'No transpose', nr, n-m, -cone, a( n-p+1, m+1 ),
333 $ lda, d( nr+1 ), 1, cone, c( n-p+1 ), 1 )
334 ELSE
335 nr = p
336 END IF
337 IF( nr.GT.0 ) THEN
338 CALL ctrmv( 'Upper', 'No transpose', 'Non unit', nr,
339 $ a( n-p+1, n-p+1 ), lda, d, 1 )
340 CALL caxpy( nr, -cone, d, 1, c( n-p+1 ), 1 )
341 END IF
342*
343* Backward transformation x = Q**H*x
344*
345 CALL cunmrq( 'Left', 'Conjugate Transpose', n, 1, p, b, ldb,
346 $ work( 1 ), x, n, work( p+mn+1 ), lwork-p-mn, info )
347 work( 1 ) = p + mn + max( lopt, int( work( p+mn+1 ) ) )
348*
349 RETURN
350*
351* End of CGGLSE
352*
subroutine xerbla(srname, info)
Definition cblat2.f:3285
subroutine caxpy(n, ca, cx, incx, cy, incy)
CAXPY
Definition caxpy.f:88
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 cggrqf(m, p, n, a, lda, taua, b, ldb, taub, work, lwork, info)
CGGRQF
Definition cggrqf.f:214
integer function ilaenv(ispec, name, opts, n1, n2, n3, n4)
ILAENV
Definition ilaenv.f:162
real function sroundup_lwork(lwork)
SROUNDUP_LWORK
subroutine ctrmv(uplo, trans, diag, n, a, lda, x, incx)
CTRMV
Definition ctrmv.f:147
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: