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

◆ dgglse()

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

DGGLSE solves overdetermined or underdetermined systems for OTHER matrices

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

Purpose:
 DGGLSE 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION array, dimension (N)
          On exit, X is the solution of the LSE problem.
[out]WORK
          WORK is DOUBLE PRECISION 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
          DGEQRF, SGERQF, DORMQR and SORMRQ.

          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 dgglse.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 DOUBLE PRECISION A( LDA, * ), B( LDB, * ), C( * ), D( * ),
190 $ WORK( * ), X( * )
191* ..
192*
193* =====================================================================
194*
195* .. Parameters ..
196 DOUBLE PRECISION ONE
197 parameter( one = 1.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 daxpy, dcopy, dgemv, dggrqf, dormqr, dormrq,
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, 'DGEQRF', ' ', m, n, -1, -1 )
242 nb2 = ilaenv( 1, 'DGERQF', ' ', m, n, -1, -1 )
243 nb3 = ilaenv( 1, 'DORMQR', ' ', m, n, p, -1 )
244 nb4 = ilaenv( 1, 'DORMRQ', ' ', 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( 'DGGLSE', -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**T = ( 0 T12 ) P Z**T*A*Q**T = ( 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* orthogonal.
276*
277 CALL dggrqf( p, m, n, b, ldb, work, a, lda, work( p+1 ),
278 $ work( p+mn+1 ), lwork-p-mn, info )
279 lopt = int( work( p+mn+1 ) )
280*
281* Update c = Z**T *c = ( c1 ) N-P
282* ( c2 ) M+P-N
283*
284 CALL dormqr( 'Left', 'Transpose', m, 1, mn, a, lda, work( p+1 ),
285 $ c, max( 1, m ), work( p+mn+1 ), lwork-p-mn, info )
286 lopt = max( lopt, int( work( p+mn+1 ) ) )
287*
288* Solve T12*x2 = d for x2
289*
290 IF( p.GT.0 ) THEN
291 CALL dtrtrs( 'Upper', 'No transpose', 'Non-unit', p, 1,
292 $ b( 1, n-p+1 ), ldb, d, p, info )
293*
294 IF( info.GT.0 ) THEN
295 info = 1
296 RETURN
297 END IF
298*
299* Put the solution in X
300*
301 CALL dcopy( p, d, 1, x( n-p+1 ), 1 )
302*
303* Update c1
304*
305 CALL dgemv( 'No transpose', n-p, p, -one, a( 1, n-p+1 ), lda,
306 $ d, 1, one, c, 1 )
307 END IF
308*
309* Solve R11*x1 = c1 for x1
310*
311 IF( n.GT.p ) THEN
312 CALL dtrtrs( 'Upper', 'No transpose', 'Non-unit', n-p, 1,
313 $ a, lda, c, n-p, info )
314*
315 IF( info.GT.0 ) THEN
316 info = 2
317 RETURN
318 END IF
319*
320* Put the solutions in X
321*
322 CALL dcopy( n-p, c, 1, x, 1 )
323 END IF
324*
325* Compute the residual vector:
326*
327 IF( m.LT.n ) THEN
328 nr = m + p - n
329 IF( nr.GT.0 )
330 $ CALL dgemv( 'No transpose', nr, n-m, -one, a( n-p+1, m+1 ),
331 $ lda, d( nr+1 ), 1, one, c( n-p+1 ), 1 )
332 ELSE
333 nr = p
334 END IF
335 IF( nr.GT.0 ) THEN
336 CALL dtrmv( 'Upper', 'No transpose', 'Non unit', nr,
337 $ a( n-p+1, n-p+1 ), lda, d, 1 )
338 CALL daxpy( nr, -one, d, 1, c( n-p+1 ), 1 )
339 END IF
340*
341* Backward transformation x = Q**T*x
342*
343 CALL dormrq( 'Left', 'Transpose', n, 1, p, b, ldb, work( 1 ), x,
344 $ n, work( p+mn+1 ), lwork-p-mn, info )
345 work( 1 ) = p + mn + max( lopt, int( work( p+mn+1 ) ) )
346*
347 RETURN
348*
349* End of DGGLSE
350*
subroutine xerbla(srname, info)
Definition cblat2.f:3285
subroutine daxpy(n, da, dx, incx, dy, incy)
DAXPY
Definition daxpy.f:89
subroutine dcopy(n, dx, incx, dy, incy)
DCOPY
Definition dcopy.f:82
subroutine dgemv(trans, m, n, alpha, a, lda, x, incx, beta, y, incy)
DGEMV
Definition dgemv.f:158
subroutine dggrqf(m, p, n, a, lda, taua, b, ldb, taub, work, lwork, info)
DGGRQF
Definition dggrqf.f:214
integer function ilaenv(ispec, name, opts, n1, n2, n3, n4)
ILAENV
Definition ilaenv.f:162
subroutine dtrmv(uplo, trans, diag, n, a, lda, x, incx)
DTRMV
Definition dtrmv.f:147
subroutine dtrtrs(uplo, trans, diag, n, nrhs, a, lda, b, ldb, info)
DTRTRS
Definition dtrtrs.f:140
subroutine dormqr(side, trans, m, n, k, a, lda, tau, c, ldc, work, lwork, info)
DORMQR
Definition dormqr.f:167
subroutine dormrq(side, trans, m, n, k, a, lda, tau, c, ldc, work, lwork, info)
DORMRQ
Definition dormrq.f:167
Here is the call graph for this function:
Here is the caller graph for this function: