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

◆ sgelst()

subroutine sgelst ( character  trans,
integer  m,
integer  n,
integer  nrhs,
real, dimension( lda, * )  a,
integer  lda,
real, dimension( ldb, * )  b,
integer  ldb,
real, dimension( * )  work,
integer  lwork,
integer  info 
)

SGELST solves overdetermined or underdetermined systems for GE matrices using QR or LQ factorization with compact WY representation of Q.

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

Purpose:
 SGELST solves overdetermined or underdetermined real linear systems
 involving an M-by-N matrix A, or its transpose, using a QR or LQ
 factorization of A with compact WY representation of Q.
 It is assumed that A has full rank.

 The following options are provided:

 1. If TRANS = 'N' and m >= n:  find the least squares solution of
    an overdetermined system, i.e., solve the least squares problem
                 minimize || B - A*X ||.

 2. If TRANS = 'N' and m < n:  find the minimum norm solution of
    an underdetermined system A * X = B.

 3. If TRANS = 'T' and m >= n:  find the minimum norm solution of
    an underdetermined system A**T * X = B.

 4. If TRANS = 'T' and m < n:  find the least squares solution of
    an overdetermined system, i.e., solve the least squares problem
                 minimize || B - A**T * X ||.

 Several right hand side vectors b and solution vectors x can be
 handled in a single call; they are stored as the columns of the
 M-by-NRHS right hand side matrix B and the N-by-NRHS solution
 matrix X.
Parameters
[in]TRANS
          TRANS is CHARACTER*1
          = 'N': the linear system involves A;
          = 'T': the linear system involves A**T.
[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 matrix A.  N >= 0.
[in]NRHS
          NRHS is INTEGER
          The number of right hand sides, i.e., the number of
          columns of the matrices B and X. NRHS >=0.
[in,out]A
          A is REAL array, dimension (LDA,N)
          On entry, the M-by-N matrix A.
          On exit,
            if M >= N, A is overwritten by details of its QR
                       factorization as returned by SGEQRT;
            if M <  N, A is overwritten by details of its LQ
                       factorization as returned by SGELQT.
[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,NRHS)
          On entry, the matrix B of right hand side vectors, stored
          columnwise; B is M-by-NRHS if TRANS = 'N', or N-by-NRHS
          if TRANS = 'T'.
          On exit, if INFO = 0, B is overwritten by the solution
          vectors, stored columnwise:
          if TRANS = 'N' and m >= n, rows 1 to n of B contain the least
          squares solution vectors; the residual sum of squares for the
          solution in each column is given by the sum of squares of
          elements N+1 to M in that column;
          if TRANS = 'N' and m < n, rows 1 to N of B contain the
          minimum norm solution vectors;
          if TRANS = 'T' and m >= n, rows 1 to M of B contain the
          minimum norm solution vectors;
          if TRANS = 'T' and m < n, rows 1 to M of B contain the
          least squares solution vectors; the residual sum of squares
          for the solution in each column is given by the sum of
          squares of elements M+1 to N in that column.
[in]LDB
          LDB is INTEGER
          The leading dimension of the array B. LDB >= MAX(1,M,N).
[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, MN + max( MN, NRHS ) ).
          For optimal performance,
          LWORK >= max( 1, (MN + max( MN, NRHS ))*NB ).
          where MN = min(M,N) and NB is the optimum block size.

          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
          > 0:  if INFO =  i, the i-th diagonal element of the
                triangular factor of A is zero, so that A does not have
                full rank; the least squares solution could not be
                computed.
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Contributors:
  November 2022,  Igor Kozachenko,
                  Computer Science Division,
                  University of California, Berkeley

Definition at line 192 of file sgelst.f.

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, SROUNDUP_LWORK
227* ..
228* .. External Subroutines ..
229 EXTERNAL sgelqt, sgeqrt, sgemlqt, sgemqrt,
231* ..
232* .. Intrinsic Functions ..
233 INTRINSIC 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 ) = sroundup_lwork( 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 ) = sroundup_lwork( 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*
313* Scale A, B if max element outside range [SMLNUM,BIGNUM]
314*
315 anrm = slange( 'M', m, n, a, lda, rwork )
316 iascl = 0
317 IF( anrm.GT.zero .AND. anrm.LT.smlnum ) THEN
318*
319* Scale matrix norm up to SMLNUM
320*
321 CALL slascl( 'G', 0, 0, anrm, smlnum, m, n, a, lda, info )
322 iascl = 1
323 ELSE IF( anrm.GT.bignum ) THEN
324*
325* Scale matrix norm down to BIGNUM
326*
327 CALL slascl( 'G', 0, 0, anrm, bignum, m, n, a, lda, info )
328 iascl = 2
329 ELSE IF( anrm.EQ.zero ) THEN
330*
331* Matrix all zero. Return zero solution.
332*
333 CALL slaset( 'Full', max( m, n ), nrhs, zero, zero, b, ldb )
334 work( 1 ) = sroundup_lwork( lwopt )
335 RETURN
336 END IF
337*
338 brow = m
339 IF( tpsd )
340 $ brow = n
341 bnrm = slange( 'M', brow, nrhs, b, ldb, rwork )
342 ibscl = 0
343 IF( bnrm.GT.zero .AND. bnrm.LT.smlnum ) THEN
344*
345* Scale matrix norm up to SMLNUM
346*
347 CALL slascl( 'G', 0, 0, bnrm, smlnum, brow, nrhs, b, ldb,
348 $ info )
349 ibscl = 1
350 ELSE IF( bnrm.GT.bignum ) THEN
351*
352* Scale matrix norm down to BIGNUM
353*
354 CALL slascl( 'G', 0, 0, bnrm, bignum, brow, nrhs, b, ldb,
355 $ info )
356 ibscl = 2
357 END IF
358*
359 IF( m.GE.n ) THEN
360*
361* M > N:
362* Compute the blocked QR factorization of A,
363* using the compact WY representation of Q,
364* workspace at least N, optimally N*NB.
365*
366 CALL sgeqrt( m, n, nb, a, lda, work( 1 ), nb,
367 $ work( mn*nb+1 ), info )
368*
369 IF( .NOT.tpsd ) THEN
370*
371* M > N, A is not transposed:
372* Overdetermined system of equations,
373* least-squares problem, min || A * X - B ||.
374*
375* Compute B(1:M,1:NRHS) := Q**T * B(1:M,1:NRHS),
376* using the compact WY representation of Q,
377* workspace at least NRHS, optimally NRHS*NB.
378*
379 CALL sgemqrt( 'Left', 'Transpose', m, nrhs, n, nb, a, lda,
380 $ work( 1 ), nb, b, ldb, work( mn*nb+1 ),
381 $ info )
382*
383* Compute B(1:N,1:NRHS) := inv(R) * B(1:N,1:NRHS)
384*
385 CALL strtrs( 'Upper', 'No transpose', 'Non-unit', n, nrhs,
386 $ a, lda, b, ldb, info )
387*
388 IF( info.GT.0 ) THEN
389 RETURN
390 END IF
391*
392 scllen = n
393*
394 ELSE
395*
396* M > N, A is transposed:
397* Underdetermined system of equations,
398* minimum norm solution of A**T * X = B.
399*
400* Compute B := inv(R**T) * B in two row blocks of B.
401*
402* Block 1: B(1:N,1:NRHS) := inv(R**T) * B(1:N,1:NRHS)
403*
404 CALL strtrs( 'Upper', 'Transpose', 'Non-unit', n, nrhs,
405 $ a, lda, b, ldb, info )
406*
407 IF( info.GT.0 ) THEN
408 RETURN
409 END IF
410*
411* Block 2: Zero out all rows below the N-th row in B:
412* B(N+1:M,1:NRHS) = ZERO
413*
414 DO j = 1, nrhs
415 DO i = n + 1, m
416 b( i, j ) = zero
417 END DO
418 END DO
419*
420* Compute B(1:M,1:NRHS) := Q(1:N,:) * B(1:N,1:NRHS),
421* using the compact WY representation of Q,
422* workspace at least NRHS, optimally NRHS*NB.
423*
424 CALL sgemqrt( 'Left', 'No transpose', m, nrhs, n, nb,
425 $ a, lda, work( 1 ), nb, b, ldb,
426 $ work( mn*nb+1 ), info )
427*
428 scllen = m
429*
430 END IF
431*
432 ELSE
433*
434* M < N:
435* Compute the blocked LQ factorization of A,
436* using the compact WY representation of Q,
437* workspace at least M, optimally M*NB.
438*
439 CALL sgelqt( m, n, nb, a, lda, work( 1 ), nb,
440 $ work( mn*nb+1 ), info )
441*
442 IF( .NOT.tpsd ) THEN
443*
444* M < N, A is not transposed:
445* Underdetermined system of equations,
446* minimum norm solution of A * X = B.
447*
448* Compute B := inv(L) * B in two row blocks of B.
449*
450* Block 1: B(1:M,1:NRHS) := inv(L) * B(1:M,1:NRHS)
451*
452 CALL strtrs( 'Lower', 'No transpose', 'Non-unit', m, nrhs,
453 $ a, lda, b, ldb, info )
454*
455 IF( info.GT.0 ) THEN
456 RETURN
457 END IF
458*
459* Block 2: Zero out all rows below the M-th row in B:
460* B(M+1:N,1:NRHS) = ZERO
461*
462 DO j = 1, nrhs
463 DO i = m + 1, n
464 b( i, j ) = zero
465 END DO
466 END DO
467*
468* Compute B(1:N,1:NRHS) := Q(1:N,:)**T * B(1:M,1:NRHS),
469* using the compact WY representation of Q,
470* workspace at least NRHS, optimally NRHS*NB.
471*
472 CALL sgemlqt( 'Left', 'Transpose', n, nrhs, m, nb, a, lda,
473 $ work( 1 ), nb, b, ldb,
474 $ work( mn*nb+1 ), info )
475*
476 scllen = n
477*
478 ELSE
479*
480* M < N, A is transposed:
481* Overdetermined system of equations,
482* least-squares problem, min || A**T * X - B ||.
483*
484* Compute B(1:N,1:NRHS) := Q * B(1:N,1:NRHS),
485* using the compact WY representation of Q,
486* workspace at least NRHS, optimally NRHS*NB.
487*
488 CALL sgemlqt( 'Left', 'No transpose', n, nrhs, m, nb,
489 $ a, lda, work( 1 ), nb, b, ldb,
490 $ work( mn*nb+1), info )
491*
492* Compute B(1:M,1:NRHS) := inv(L**T) * B(1:M,1:NRHS)
493*
494 CALL strtrs( 'Lower', 'Transpose', 'Non-unit', m, nrhs,
495 $ a, lda, b, ldb, info )
496*
497 IF( info.GT.0 ) THEN
498 RETURN
499 END IF
500*
501 scllen = m
502*
503 END IF
504*
505 END IF
506*
507* Undo scaling
508*
509 IF( iascl.EQ.1 ) THEN
510 CALL slascl( 'G', 0, 0, anrm, smlnum, scllen, nrhs, b, ldb,
511 $ info )
512 ELSE IF( iascl.EQ.2 ) THEN
513 CALL slascl( 'G', 0, 0, anrm, bignum, scllen, nrhs, b, ldb,
514 $ info )
515 END IF
516 IF( ibscl.EQ.1 ) THEN
517 CALL slascl( 'G', 0, 0, smlnum, bnrm, scllen, nrhs, b, ldb,
518 $ info )
519 ELSE IF( ibscl.EQ.2 ) THEN
520 CALL slascl( 'G', 0, 0, bignum, bnrm, scllen, nrhs, b, ldb,
521 $ info )
522 END IF
523*
524 work( 1 ) = sroundup_lwork( lwopt )
525*
526 RETURN
527*
528* End of SGELST
529*
subroutine xerbla(srname, info)
Definition cblat2.f:3285
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 sgemqrt(side, trans, m, n, k, nb, v, ldv, t, ldt, c, ldc, work, info)
SGEMQRT
Definition sgemqrt.f:168
subroutine sgeqrt(m, n, nb, a, lda, t, ldt, work, info)
SGEQRT
Definition sgeqrt.f:141
integer function ilaenv(ispec, name, opts, n1, n2, n3, n4)
ILAENV
Definition ilaenv.f:162
real function slamch(cmach)
SLAMCH
Definition slamch.f:68
real function slange(norm, m, n, a, lda, work)
SLANGE returns the value of the 1-norm, Frobenius norm, infinity-norm, or the largest absolute value ...
Definition slange.f:114
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
logical function lsame(ca, cb)
LSAME
Definition lsame.f:48
real function sroundup_lwork(lwork)
SROUNDUP_LWORK
subroutine strtrs(uplo, trans, diag, n, nrhs, a, lda, b, ldb, info)
STRTRS
Definition strtrs.f:140
Here is the call graph for this function:
Here is the caller graph for this function: