LAPACK 3.11.0 LAPACK: Linear Algebra PACKage
Searching...
No Matches
sgelst.f
Go to the documentation of this file.
1*> \brief <b> SGELST solves overdetermined or underdetermined systems for GE matrices using QR or LQ factorization with compact WY representation of Q.</b>
2*
3* =========== DOCUMENTATION ===========
4*
5* Online html documentation available at
6* http://www.netlib.org/lapack/explore-html/
7*
8*> \htmlonly
10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/sgelst.f">
11*> [TGZ]</a>
12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/sgelst.f">
13*> [ZIP]</a>
14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/sgelst.f">
15*> [TXT]</a>
16*> \endhtmlonly
17*
18* Definition:
19* ===========
20*
21* SUBROUTINE SGELST( TRANS, M, N, NRHS, A, LDA, B, LDB, WORK, LWORK,
22* INFO )
23*
24* .. Scalar Arguments ..
25* CHARACTER TRANS
26* INTEGER INFO, LDA, LDB, LWORK, M, N, NRHS
27* ..
28* .. Array Arguments ..
29* REAL A( LDA, * ), B( LDB, * ), WORK( * )
30* ..
31*
32*
33*> \par Purpose:
34* =============
35*>
36*> \verbatim
37*>
38*> SGELST solves overdetermined or underdetermined real linear systems
39*> involving an M-by-N matrix A, or its transpose, using a QR or LQ
40*> factorization of A with compact WY representation of Q.
41*> It is assumed that A has full rank.
42*>
43*> The following options are provided:
44*>
45*> 1. If TRANS = 'N' and m >= n: find the least squares solution of
46*> an overdetermined system, i.e., solve the least squares problem
47*> minimize || B - A*X ||.
48*>
49*> 2. If TRANS = 'N' and m < n: find the minimum norm solution of
50*> an underdetermined system A * X = B.
51*>
52*> 3. If TRANS = 'T' and m >= n: find the minimum norm solution of
53*> an underdetermined system A**T * X = B.
54*>
55*> 4. If TRANS = 'T' and m < n: find the least squares solution of
56*> an overdetermined system, i.e., solve the least squares problem
57*> minimize || B - A**T * X ||.
58*>
59*> Several right hand side vectors b and solution vectors x can be
60*> handled in a single call; they are stored as the columns of the
61*> M-by-NRHS right hand side matrix B and the N-by-NRHS solution
62*> matrix X.
63*> \endverbatim
64*
65* Arguments:
66* ==========
67*
68*> \param[in] TRANS
69*> \verbatim
70*> TRANS is CHARACTER*1
71*> = 'N': the linear system involves A;
72*> = 'T': the linear system involves A**T.
73*> \endverbatim
74*>
75*> \param[in] M
76*> \verbatim
77*> M is INTEGER
78*> The number of rows of the matrix A. M >= 0.
79*> \endverbatim
80*>
81*> \param[in] N
82*> \verbatim
83*> N is INTEGER
84*> The number of columns of the matrix A. N >= 0.
85*> \endverbatim
86*>
87*> \param[in] NRHS
88*> \verbatim
89*> NRHS is INTEGER
90*> The number of right hand sides, i.e., the number of
91*> columns of the matrices B and X. NRHS >=0.
92*> \endverbatim
93*>
94*> \param[in,out] A
95*> \verbatim
96*> A is REAL array, dimension (LDA,N)
97*> On entry, the M-by-N matrix A.
98*> On exit,
99*> if M >= N, A is overwritten by details of its QR
100*> factorization as returned by SGEQRT;
101*> if M < N, A is overwritten by details of its LQ
102*> factorization as returned by SGELQT.
103*> \endverbatim
104*>
105*> \param[in] LDA
106*> \verbatim
107*> LDA is INTEGER
108*> The leading dimension of the array A. LDA >= max(1,M).
109*> \endverbatim
110*>
111*> \param[in,out] B
112*> \verbatim
113*> B is REAL array, dimension (LDB,NRHS)
114*> On entry, the matrix B of right hand side vectors, stored
115*> columnwise; B is M-by-NRHS if TRANS = 'N', or N-by-NRHS
116*> if TRANS = 'T'.
117*> On exit, if INFO = 0, B is overwritten by the solution
118*> vectors, stored columnwise:
119*> if TRANS = 'N' and m >= n, rows 1 to n of B contain the least
120*> squares solution vectors; the residual sum of squares for the
121*> solution in each column is given by the sum of squares of
122*> elements N+1 to M in that column;
123*> if TRANS = 'N' and m < n, rows 1 to N of B contain the
124*> minimum norm solution vectors;
125*> if TRANS = 'T' and m >= n, rows 1 to M of B contain the
126*> minimum norm solution vectors;
127*> if TRANS = 'T' and m < n, rows 1 to M of B contain the
128*> least squares solution vectors; the residual sum of squares
129*> for the solution in each column is given by the sum of
130*> squares of elements M+1 to N in that column.
131*> \endverbatim
132*>
133*> \param[in] LDB
134*> \verbatim
135*> LDB is INTEGER
136*> The leading dimension of the array B. LDB >= MAX(1,M,N).
137*> \endverbatim
138*>
139*> \param[out] WORK
140*> \verbatim
141*> WORK is REAL array, dimension (MAX(1,LWORK))
142*> On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
143*> \endverbatim
144*>
145*> \param[in] LWORK
146*> \verbatim
147*> LWORK is INTEGER
148*> The dimension of the array WORK.
149*> LWORK >= max( 1, MN + max( MN, NRHS ) ).
150*> For optimal performance,
151*> LWORK >= max( 1, (MN + max( MN, NRHS ))*NB ).
152*> where MN = min(M,N) and NB is the optimum block size.
153*>
154*> If LWORK = -1, then a workspace query is assumed; the routine
155*> only calculates the optimal size of the WORK array, returns
156*> this value as the first entry of the WORK array, and no error
157*> message related to LWORK is issued by XERBLA.
158*> \endverbatim
159*>
160*> \param[out] INFO
161*> \verbatim
162*> INFO is INTEGER
163*> = 0: successful exit
164*> < 0: if INFO = -i, the i-th argument had an illegal value
165*> > 0: if INFO = i, the i-th diagonal element of the
166*> triangular factor of A is zero, so that A does not have
167*> full rank; the least squares solution could not be
168*> computed.
169*> \endverbatim
170*
171* Authors:
172* ========
173*
174*> \author Univ. of Tennessee
175*> \author Univ. of California Berkeley
176*> \author Univ. of Colorado Denver
177*> \author NAG Ltd.
178*
179*> \ingroup realGEsolve
180*
181*> \par Contributors:
182* ==================
183*>
184*> \verbatim
185*>
186*> November 2022, Igor Kozachenko,
187*> Computer Science Division,
188*> University of California, Berkeley
189*> \endverbatim
190*
191* =====================================================================
192 SUBROUTINE sgelst( TRANS, M, N, NRHS, A, LDA, B, LDB, WORK, LWORK,
193 \$ INFO )
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 CHARACTER TRANS
201 INTEGER INFO, LDA, LDB, LWORK, M, N, NRHS
202* ..
203* .. Array Arguments ..
204 REAL A( LDA, * ), B( LDB, * ), WORK( * )
205* ..
206*
207* =====================================================================
208*
209* .. Parameters ..
210 REAL ZERO, ONE
211 parameter( zero = 0.0e+0, one = 1.0e+0 )
212* ..
213* .. Local Scalars ..
214 LOGICAL LQUERY, TPSD
215 INTEGER BROW, I, IASCL, IBSCL, J, LWOPT, MN, MNNRHS,
216 \$ nb, nbmin, scllen
217 REAL ANRM, BIGNUM, BNRM, SMLNUM
218* ..
219* .. Local Arrays ..
220 REAL RWORK( 1 )
221* ..
222* .. External Functions ..
223 LOGICAL LSAME
224 INTEGER ILAENV
225 REAL SLAMCH, SLANGE
226 EXTERNAL lsame, ilaenv, slamch, slange
227* ..
228* .. External Subroutines ..
229 EXTERNAL sgelqt, sgeqrt, sgemlqt, sgemqrt, slabad,
231* ..
232* .. Intrinsic Functions ..
233 INTRINSIC real, max, min
234* ..
235* .. Executable Statements ..
236*
237* Test the input arguments.
238*
239 info = 0
240 mn = min( m, n )
241 lquery = ( lwork.EQ.-1 )
242 IF( .NOT.( lsame( trans, 'N' ) .OR. lsame( trans, 'T' ) ) ) THEN
243 info = -1
244 ELSE IF( m.LT.0 ) THEN
245 info = -2
246 ELSE IF( n.LT.0 ) THEN
247 info = -3
248 ELSE IF( nrhs.LT.0 ) THEN
249 info = -4
250 ELSE IF( lda.LT.max( 1, m ) ) THEN
251 info = -6
252 ELSE IF( ldb.LT.max( 1, m, n ) ) THEN
253 info = -8
254 ELSE IF( lwork.LT.max( 1, mn+max( mn, nrhs ) ) .AND. .NOT.lquery )
255 \$ THEN
256 info = -10
257 END IF
258*
259* Figure out optimal block size and optimal workspace size
260*
261 IF( info.EQ.0 .OR. info.EQ.-10 ) THEN
262*
263 tpsd = .true.
264 IF( lsame( trans, 'N' ) )
265 \$ tpsd = .false.
266*
267 nb = ilaenv( 1, 'SGELST', ' ', m, n, -1, -1 )
268*
269 mnnrhs = max( mn, nrhs )
270 lwopt = max( 1, (mn+mnnrhs)*nb )
271 work( 1 ) = real( lwopt )
272*
273 END IF
274*
275 IF( info.NE.0 ) THEN
276 CALL xerbla( 'SGELST ', -info )
277 RETURN
278 ELSE IF( lquery ) THEN
279 RETURN
280 END IF
281*
282* Quick return if possible
283*
284 IF( min( m, n, nrhs ).EQ.0 ) THEN
285 CALL slaset( 'Full', max( m, n ), nrhs, zero, zero, b, ldb )
286 work( 1 ) = real( lwopt )
287 RETURN
288 END IF
289*
290* *GEQRT and *GELQT routines cannot accept NB larger than min(M,N)
291*
292 IF( nb.GT.mn ) nb = mn
293*
294* Determine the block size from the supplied LWORK
295* ( at this stage we know that LWORK >= (minimum required workspace,
296* but it may be less than optimal)
297*
298 nb = min( nb, lwork/( mn + mnnrhs ) )
299*
300* The minimum value of NB, when blocked code is used
301*
302 nbmin = max( 2, ilaenv( 2, 'SGELST', ' ', m, n, -1, -1 ) )
303*
304 IF( nb.LT.nbmin ) THEN
305 nb = 1
306 END IF
307*
308* Get machine parameters
309*
310 smlnum = slamch( 'S' ) / slamch( 'P' )
311 bignum = one / smlnum
312 CALL slabad( smlnum, bignum )
313*
314* Scale A, B if max element outside range [SMLNUM,BIGNUM]
315*
316 anrm = slange( 'M', m, n, a, lda, rwork )
317 iascl = 0
318 IF( anrm.GT.zero .AND. anrm.LT.smlnum ) THEN
319*
320* Scale matrix norm up to SMLNUM
321*
322 CALL slascl( 'G', 0, 0, anrm, smlnum, m, n, a, lda, info )
323 iascl = 1
324 ELSE IF( anrm.GT.bignum ) THEN
325*
326* Scale matrix norm down to BIGNUM
327*
328 CALL slascl( 'G', 0, 0, anrm, bignum, m, n, a, lda, info )
329 iascl = 2
330 ELSE IF( anrm.EQ.zero ) THEN
331*
332* Matrix all zero. Return zero solution.
333*
334 CALL slaset( 'Full', max( m, n ), nrhs, zero, zero, b, ldb )
335 work( 1 ) = real( lwopt )
336 RETURN
337 END IF
338*
339 brow = m
340 IF( tpsd )
341 \$ brow = n
342 bnrm = slange( 'M', brow, nrhs, b, ldb, rwork )
343 ibscl = 0
344 IF( bnrm.GT.zero .AND. bnrm.LT.smlnum ) THEN
345*
346* Scale matrix norm up to SMLNUM
347*
348 CALL slascl( 'G', 0, 0, bnrm, smlnum, brow, nrhs, b, ldb,
349 \$ info )
350 ibscl = 1
351 ELSE IF( bnrm.GT.bignum ) THEN
352*
353* Scale matrix norm down to BIGNUM
354*
355 CALL slascl( 'G', 0, 0, bnrm, bignum, brow, nrhs, b, ldb,
356 \$ info )
357 ibscl = 2
358 END IF
359*
360 IF( m.GE.n ) THEN
361*
362* M > N:
363* Compute the blocked QR factorization of A,
364* using the compact WY representation of Q,
365* workspace at least N, optimally N*NB.
366*
367 CALL sgeqrt( m, n, nb, a, lda, work( 1 ), nb,
368 \$ work( mn*nb+1 ), info )
369*
370 IF( .NOT.tpsd ) THEN
371*
372* M > N, A is not transposed:
373* Overdetermined system of equations,
374* least-squares problem, min || A * X - B ||.
375*
376* Compute B(1:M,1:NRHS) := Q**T * B(1:M,1:NRHS),
377* using the compact WY representation of Q,
378* workspace at least NRHS, optimally NRHS*NB.
379*
380 CALL sgemqrt( 'Left', 'Transpose', m, nrhs, n, nb, a, lda,
381 \$ work( 1 ), nb, b, ldb, work( mn*nb+1 ),
382 \$ info )
383*
384* Compute B(1:N,1:NRHS) := inv(R) * B(1:N,1:NRHS)
385*
386 CALL strtrs( 'Upper', 'No transpose', 'Non-unit', n, nrhs,
387 \$ a, lda, b, ldb, info )
388*
389 IF( info.GT.0 ) THEN
390 RETURN
391 END IF
392*
393 scllen = n
394*
395 ELSE
396*
397* M > N, A is transposed:
398* Underdetermined system of equations,
399* minimum norm solution of A**T * X = B.
400*
401* Compute B := inv(R**T) * B in two row blocks of B.
402*
403* Block 1: B(1:N,1:NRHS) := inv(R**T) * B(1:N,1:NRHS)
404*
405 CALL strtrs( 'Upper', 'Transpose', 'Non-unit', n, nrhs,
406 \$ a, lda, b, ldb, info )
407*
408 IF( info.GT.0 ) THEN
409 RETURN
410 END IF
411*
412* Block 2: Zero out all rows below the N-th row in B:
413* B(N+1:M,1:NRHS) = ZERO
414*
415 DO j = 1, nrhs
416 DO i = n + 1, m
417 b( i, j ) = zero
418 END DO
419 END DO
420*
421* Compute B(1:M,1:NRHS) := Q(1:N,:) * B(1:N,1:NRHS),
422* using the compact WY representation of Q,
423* workspace at least NRHS, optimally NRHS*NB.
424*
425 CALL sgemqrt( 'Left', 'No transpose', m, nrhs, n, nb,
426 \$ a, lda, work( 1 ), nb, b, ldb,
427 \$ work( mn*nb+1 ), info )
428*
429 scllen = m
430*
431 END IF
432*
433 ELSE
434*
435* M < N:
436* Compute the blocked LQ factorization of A,
437* using the compact WY representation of Q,
438* workspace at least M, optimally M*NB.
439*
440 CALL sgelqt( m, n, nb, a, lda, work( 1 ), nb,
441 \$ work( mn*nb+1 ), info )
442*
443 IF( .NOT.tpsd ) THEN
444*
445* M < N, A is not transposed:
446* Underdetermined system of equations,
447* minimum norm solution of A * X = B.
448*
449* Compute B := inv(L) * B in two row blocks of B.
450*
451* Block 1: B(1:M,1:NRHS) := inv(L) * B(1:M,1:NRHS)
452*
453 CALL strtrs( 'Lower', 'No transpose', 'Non-unit', m, nrhs,
454 \$ a, lda, b, ldb, info )
455*
456 IF( info.GT.0 ) THEN
457 RETURN
458 END IF
459*
460* Block 2: Zero out all rows below the M-th row in B:
461* B(M+1:N,1:NRHS) = ZERO
462*
463 DO j = 1, nrhs
464 DO i = m + 1, n
465 b( i, j ) = zero
466 END DO
467 END DO
468*
469* Compute B(1:N,1:NRHS) := Q(1:N,:)**T * B(1:M,1:NRHS),
470* using the compact WY representation of Q,
471* workspace at least NRHS, optimally NRHS*NB.
472*
473 CALL sgemlqt( 'Left', 'Transpose', n, nrhs, m, nb, a, lda,
474 \$ work( 1 ), nb, b, ldb,
475 \$ work( mn*nb+1 ), info )
476*
477 scllen = n
478*
479 ELSE
480*
481* M < N, A is transposed:
482* Overdetermined system of equations,
483* least-squares problem, min || A**T * X - B ||.
484*
485* Compute B(1:N,1:NRHS) := Q * B(1:N,1:NRHS),
486* using the compact WY representation of Q,
487* workspace at least NRHS, optimally NRHS*NB.
488*
489 CALL sgemlqt( 'Left', 'No transpose', n, nrhs, m, nb,
490 \$ a, lda, work( 1 ), nb, b, ldb,
491 \$ work( mn*nb+1), info )
492*
493* Compute B(1:M,1:NRHS) := inv(L**T) * B(1:M,1:NRHS)
494*
495 CALL strtrs( 'Lower', 'Transpose', 'Non-unit', m, nrhs,
496 \$ a, lda, b, ldb, info )
497*
498 IF( info.GT.0 ) THEN
499 RETURN
500 END IF
501*
502 scllen = m
503*
504 END IF
505*
506 END IF
507*
508* Undo scaling
509*
510 IF( iascl.EQ.1 ) THEN
511 CALL slascl( 'G', 0, 0, anrm, smlnum, scllen, nrhs, b, ldb,
512 \$ info )
513 ELSE IF( iascl.EQ.2 ) THEN
514 CALL slascl( 'G', 0, 0, anrm, bignum, scllen, nrhs, b, ldb,
515 \$ info )
516 END IF
517 IF( ibscl.EQ.1 ) THEN
518 CALL slascl( 'G', 0, 0, smlnum, bnrm, scllen, nrhs, b, ldb,
519 \$ info )
520 ELSE IF( ibscl.EQ.2 ) THEN
521 CALL slascl( 'G', 0, 0, bignum, bnrm, scllen, nrhs, b, ldb,
522 \$ info )
523 END IF
524*
525 work( 1 ) = real( lwopt )
526*
527 RETURN
528*
529* End of SGELST
530*
531 END
subroutine slascl(TYPE, KL, KU, CFROM, CTO, M, N, A, LDA, INFO)
SLASCL multiplies a general rectangular matrix by a real scalar defined as cto/cfrom.
Definition: slascl.f:143
subroutine slaset(UPLO, M, N, ALPHA, BETA, A, LDA)
SLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values.
Definition: slaset.f:110
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:60
subroutine sgelqt(M, N, MB, A, LDA, T, LDT, WORK, INFO)
SGELQT
Definition: sgelqt.f:124
subroutine sgemlqt(SIDE, TRANS, M, N, K, MB, V, LDV, T, LDT, C, LDC, WORK, INFO)
SGEMLQT
Definition: sgemlqt.f:153
subroutine sgeqrt(M, N, NB, A, LDA, T, LDT, WORK, INFO)
SGEQRT
Definition: sgeqrt.f:141
subroutine sgemqrt(SIDE, TRANS, M, N, K, NB, V, LDV, T, LDT, C, LDC, WORK, INFO)
SGEMQRT
Definition: sgemqrt.f:168
subroutine sgelst(TRANS, M, N, NRHS, A, LDA, B, LDB, WORK, LWORK, INFO)
SGELST solves overdetermined or underdetermined systems for GE matrices using QR or LQ factorization ...
Definition: sgelst.f:194
subroutine strtrs(UPLO, TRANS, DIAG, N, NRHS, A, LDA, B, LDB, INFO)
STRTRS
Definition: strtrs.f:140