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

◆ sgglse()

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

SGGLSE solves overdetermined or underdetermined systems for OTHER matrices

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

Purpose:
 SGGLSE 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 REAL 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 REAL 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 REAL 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 REAL 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 REAL array, dimension (N)
          On exit, X is the solution of the LSE problem.
[out]WORK
          WORK is REAL 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
          SGEQRF, SGERQF, SORMQR 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 sgglse.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 REAL A( LDA, * ), B( LDB, * ), C( * ), D( * ),
190 $ WORK( * ), X( * )
191* ..
192*
193* =====================================================================
194*
195* .. Parameters ..
196 REAL ONE
197 parameter( one = 1.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 saxpy, scopy, sgemv, sggrqf, sormqr, sormrq,
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, 'SGEQRF', ' ', m, n, -1, -1 )
243 nb2 = ilaenv( 1, 'SGERQF', ' ', m, n, -1, -1 )
244 nb3 = ilaenv( 1, 'SORMQR', ' ', m, n, p, -1 )
245 nb4 = ilaenv( 1, 'SORMRQ', ' ', 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( 'SGGLSE', -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**T = ( 0 T12 ) P Z**T*A*Q**T = ( 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* orthogonal.
277*
278 CALL sggrqf( 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**T *c = ( c1 ) N-P
283* ( c2 ) M+P-N
284*
285 CALL sormqr( 'Left', 'Transpose', m, 1, mn, a, lda, work( p+1 ),
286 $ c, max( 1, m ), work( p+mn+1 ), 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 strtrs( '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 scopy( p, d, 1, x( n-p+1 ), 1 )
303*
304* Update c1
305*
306 CALL sgemv( 'No transpose', n-p, p, -one, a( 1, n-p+1 ), lda,
307 $ d, 1, one, c, 1 )
308 END IF
309*
310* Solve R11*x1 = c1 for x1
311*
312 IF( n.GT.p ) THEN
313 CALL strtrs( '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 scopy( 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 sgemv( 'No transpose', nr, n-m, -one, a( n-p+1, m+1 ),
332 $ lda, d( nr+1 ), 1, one, c( n-p+1 ), 1 )
333 ELSE
334 nr = p
335 END IF
336 IF( nr.GT.0 ) THEN
337 CALL strmv( 'Upper', 'No transpose', 'Non unit', nr,
338 $ a( n-p+1, n-p+1 ), lda, d, 1 )
339 CALL saxpy( nr, -one, d, 1, c( n-p+1 ), 1 )
340 END IF
341*
342* Backward transformation x = Q**T*x
343*
344 CALL sormrq( 'Left', 'Transpose', n, 1, p, b, ldb, work( 1 ), x,
345 $ 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 SGGLSE
351*
subroutine xerbla(srname, info)
Definition cblat2.f:3285
subroutine saxpy(n, sa, sx, incx, sy, incy)
SAXPY
Definition saxpy.f:89
subroutine scopy(n, sx, incx, sy, incy)
SCOPY
Definition scopy.f:82
subroutine sgemv(trans, m, n, alpha, a, lda, x, incx, beta, y, incy)
SGEMV
Definition sgemv.f:158
subroutine sggrqf(m, p, n, a, lda, taua, b, ldb, taub, work, lwork, info)
SGGRQF
Definition sggrqf.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 strmv(uplo, trans, diag, n, a, lda, x, incx)
STRMV
Definition strmv.f:147
subroutine strtrs(uplo, trans, diag, n, nrhs, a, lda, b, ldb, info)
STRTRS
Definition strtrs.f:140
subroutine sormqr(side, trans, m, n, k, a, lda, tau, c, ldc, work, lwork, info)
SORMQR
Definition sormqr.f:168
subroutine sormrq(side, trans, m, n, k, a, lda, tau, c, ldc, work, lwork, info)
SORMRQ
Definition sormrq.f:168
Here is the call graph for this function:
Here is the caller graph for this function: