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

◆ sgels()

subroutine sgels ( 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 
)

SGELS solves overdetermined or underdetermined systems for GE matrices

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

Purpose:
 SGELS 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.  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 SGEQRF;
            if M <  N, A is overwritten by details of its LQ
                       factorization as returned by SGELQF.
[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.

Definition at line 181 of file sgels.f.

183*
184* -- LAPACK driver routine --
185* -- LAPACK is a software package provided by Univ. of Tennessee, --
186* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
187*
188* .. Scalar Arguments ..
189 CHARACTER TRANS
190 INTEGER INFO, LDA, LDB, LWORK, M, N, NRHS
191* ..
192* .. Array Arguments ..
193 REAL A( LDA, * ), B( LDB, * ), WORK( * )
194* ..
195*
196* =====================================================================
197*
198* .. Parameters ..
199 REAL ZERO, ONE
200 parameter( zero = 0.0e0, one = 1.0e0 )
201* ..
202* .. Local Scalars ..
203 LOGICAL LQUERY, TPSD
204 INTEGER BROW, I, IASCL, IBSCL, J, MN, NB, SCLLEN, WSIZE
205 REAL ANRM, BIGNUM, BNRM, SMLNUM
206* ..
207* .. Local Arrays ..
208 REAL RWORK( 1 )
209* ..
210* .. External Functions ..
211 LOGICAL LSAME
212 INTEGER ILAENV
213 REAL SLAMCH, SLANGE
214 EXTERNAL lsame, ilaenv, slamch, slange
215* ..
216* .. External Subroutines ..
217 EXTERNAL sgelqf, sgeqrf, slabad, slascl, slaset, sormlq,
219* ..
220* .. Intrinsic Functions ..
221 INTRINSIC max, min, real
222* ..
223* .. Executable Statements ..
224*
225* Test the input arguments.
226*
227 info = 0
228 mn = min( m, n )
229 lquery = ( lwork.EQ.-1 )
230 IF( .NOT.( lsame( trans, 'N' ) .OR. lsame( trans, 'T' ) ) ) THEN
231 info = -1
232 ELSE IF( m.LT.0 ) THEN
233 info = -2
234 ELSE IF( n.LT.0 ) THEN
235 info = -3
236 ELSE IF( nrhs.LT.0 ) THEN
237 info = -4
238 ELSE IF( lda.LT.max( 1, m ) ) THEN
239 info = -6
240 ELSE IF( ldb.LT.max( 1, m, n ) ) THEN
241 info = -8
242 ELSE IF( lwork.LT.max( 1, mn + max( mn, nrhs ) ) .AND.
243 $ .NOT.lquery ) THEN
244 info = -10
245 END IF
246*
247* Figure out optimal block size
248*
249 IF( info.EQ.0 .OR. info.EQ.-10 ) THEN
250*
251 tpsd = .true.
252 IF( lsame( trans, 'N' ) )
253 $ tpsd = .false.
254*
255 IF( m.GE.n ) THEN
256 nb = ilaenv( 1, 'SGEQRF', ' ', m, n, -1, -1 )
257 IF( tpsd ) THEN
258 nb = max( nb, ilaenv( 1, 'SORMQR', 'LN', m, nrhs, n,
259 $ -1 ) )
260 ELSE
261 nb = max( nb, ilaenv( 1, 'SORMQR', 'LT', m, nrhs, n,
262 $ -1 ) )
263 END IF
264 ELSE
265 nb = ilaenv( 1, 'SGELQF', ' ', m, n, -1, -1 )
266 IF( tpsd ) THEN
267 nb = max( nb, ilaenv( 1, 'SORMLQ', 'LT', n, nrhs, m,
268 $ -1 ) )
269 ELSE
270 nb = max( nb, ilaenv( 1, 'SORMLQ', 'LN', n, nrhs, m,
271 $ -1 ) )
272 END IF
273 END IF
274*
275 wsize = max( 1, mn + max( mn, nrhs )*nb )
276 work( 1 ) = real( wsize )
277*
278 END IF
279*
280 IF( info.NE.0 ) THEN
281 CALL xerbla( 'SGELS ', -info )
282 RETURN
283 ELSE IF( lquery ) THEN
284 RETURN
285 END IF
286*
287* Quick return if possible
288*
289 IF( min( m, n, nrhs ).EQ.0 ) THEN
290 CALL slaset( 'Full', max( m, n ), nrhs, zero, zero, b, ldb )
291 RETURN
292 END IF
293*
294* Get machine parameters
295*
296 smlnum = slamch( 'S' ) / slamch( 'P' )
297 bignum = one / smlnum
298 CALL slabad( smlnum, bignum )
299*
300* Scale A, B if max element outside range [SMLNUM,BIGNUM]
301*
302 anrm = slange( 'M', m, n, a, lda, rwork )
303 iascl = 0
304 IF( anrm.GT.zero .AND. anrm.LT.smlnum ) THEN
305*
306* Scale matrix norm up to SMLNUM
307*
308 CALL slascl( 'G', 0, 0, anrm, smlnum, m, n, a, lda, info )
309 iascl = 1
310 ELSE IF( anrm.GT.bignum ) THEN
311*
312* Scale matrix norm down to BIGNUM
313*
314 CALL slascl( 'G', 0, 0, anrm, bignum, m, n, a, lda, info )
315 iascl = 2
316 ELSE IF( anrm.EQ.zero ) THEN
317*
318* Matrix all zero. Return zero solution.
319*
320 CALL slaset( 'F', max( m, n ), nrhs, zero, zero, b, ldb )
321 GO TO 50
322 END IF
323*
324 brow = m
325 IF( tpsd )
326 $ brow = n
327 bnrm = slange( 'M', brow, nrhs, b, ldb, rwork )
328 ibscl = 0
329 IF( bnrm.GT.zero .AND. bnrm.LT.smlnum ) THEN
330*
331* Scale matrix norm up to SMLNUM
332*
333 CALL slascl( 'G', 0, 0, bnrm, smlnum, brow, nrhs, b, ldb,
334 $ info )
335 ibscl = 1
336 ELSE IF( bnrm.GT.bignum ) THEN
337*
338* Scale matrix norm down to BIGNUM
339*
340 CALL slascl( 'G', 0, 0, bnrm, bignum, brow, nrhs, b, ldb,
341 $ info )
342 ibscl = 2
343 END IF
344*
345 IF( m.GE.n ) THEN
346*
347* compute QR factorization of A
348*
349 CALL sgeqrf( m, n, a, lda, work( 1 ), work( mn+1 ), lwork-mn,
350 $ info )
351*
352* workspace at least N, optimally N*NB
353*
354 IF( .NOT.tpsd ) THEN
355*
356* Least-Squares Problem min || A * X - B ||
357*
358* B(1:M,1:NRHS) := Q**T * B(1:M,1:NRHS)
359*
360 CALL sormqr( 'Left', 'Transpose', m, nrhs, n, a, lda,
361 $ work( 1 ), b, ldb, work( mn+1 ), lwork-mn,
362 $ info )
363*
364* workspace at least NRHS, optimally NRHS*NB
365*
366* B(1:N,1:NRHS) := inv(R) * B(1:N,1:NRHS)
367*
368 CALL strtrs( 'Upper', 'No transpose', 'Non-unit', n, nrhs,
369 $ a, lda, b, ldb, info )
370*
371 IF( info.GT.0 ) THEN
372 RETURN
373 END IF
374*
375 scllen = n
376*
377 ELSE
378*
379* Underdetermined system of equations A**T * X = B
380*
381* B(1:N,1:NRHS) := inv(R**T) * B(1:N,1:NRHS)
382*
383 CALL strtrs( 'Upper', 'Transpose', 'Non-unit', n, nrhs,
384 $ a, lda, b, ldb, info )
385*
386 IF( info.GT.0 ) THEN
387 RETURN
388 END IF
389*
390* B(N+1:M,1:NRHS) = ZERO
391*
392 DO 20 j = 1, nrhs
393 DO 10 i = n + 1, m
394 b( i, j ) = zero
395 10 CONTINUE
396 20 CONTINUE
397*
398* B(1:M,1:NRHS) := Q(1:N,:) * B(1:N,1:NRHS)
399*
400 CALL sormqr( 'Left', 'No transpose', m, nrhs, n, a, lda,
401 $ work( 1 ), b, ldb, work( mn+1 ), lwork-mn,
402 $ info )
403*
404* workspace at least NRHS, optimally NRHS*NB
405*
406 scllen = m
407*
408 END IF
409*
410 ELSE
411*
412* Compute LQ factorization of A
413*
414 CALL sgelqf( m, n, a, lda, work( 1 ), work( mn+1 ), lwork-mn,
415 $ info )
416*
417* workspace at least M, optimally M*NB.
418*
419 IF( .NOT.tpsd ) THEN
420*
421* underdetermined system of equations A * X = B
422*
423* B(1:M,1:NRHS) := inv(L) * B(1:M,1:NRHS)
424*
425 CALL strtrs( 'Lower', 'No transpose', 'Non-unit', m, nrhs,
426 $ a, lda, b, ldb, info )
427*
428 IF( info.GT.0 ) THEN
429 RETURN
430 END IF
431*
432* B(M+1:N,1:NRHS) = 0
433*
434 DO 40 j = 1, nrhs
435 DO 30 i = m + 1, n
436 b( i, j ) = zero
437 30 CONTINUE
438 40 CONTINUE
439*
440* B(1:N,1:NRHS) := Q(1:N,:)**T * B(1:M,1:NRHS)
441*
442 CALL sormlq( 'Left', 'Transpose', n, nrhs, m, a, lda,
443 $ work( 1 ), b, ldb, work( mn+1 ), lwork-mn,
444 $ info )
445*
446* workspace at least NRHS, optimally NRHS*NB
447*
448 scllen = n
449*
450 ELSE
451*
452* overdetermined system min || A**T * X - B ||
453*
454* B(1:N,1:NRHS) := Q * B(1:N,1:NRHS)
455*
456 CALL sormlq( 'Left', 'No transpose', n, nrhs, m, a, lda,
457 $ work( 1 ), b, ldb, work( mn+1 ), lwork-mn,
458 $ info )
459*
460* workspace at least NRHS, optimally NRHS*NB
461*
462* B(1:M,1:NRHS) := inv(L**T) * B(1:M,1:NRHS)
463*
464 CALL strtrs( 'Lower', 'Transpose', 'Non-unit', m, nrhs,
465 $ a, lda, b, ldb, info )
466*
467 IF( info.GT.0 ) THEN
468 RETURN
469 END IF
470*
471 scllen = m
472*
473 END IF
474*
475 END IF
476*
477* Undo scaling
478*
479 IF( iascl.EQ.1 ) THEN
480 CALL slascl( 'G', 0, 0, anrm, smlnum, scllen, nrhs, b, ldb,
481 $ info )
482 ELSE IF( iascl.EQ.2 ) THEN
483 CALL slascl( 'G', 0, 0, anrm, bignum, scllen, nrhs, b, ldb,
484 $ info )
485 END IF
486 IF( ibscl.EQ.1 ) THEN
487 CALL slascl( 'G', 0, 0, smlnum, bnrm, scllen, nrhs, b, ldb,
488 $ info )
489 ELSE IF( ibscl.EQ.2 ) THEN
490 CALL slascl( 'G', 0, 0, bignum, bnrm, scllen, nrhs, b, ldb,
491 $ info )
492 END IF
493*
494 50 CONTINUE
495 work( 1 ) = real( wsize )
496*
497 RETURN
498*
499* End of SGELS
500*
subroutine slabad(SMALL, LARGE)
SLABAD
Definition: slabad.f:74
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
integer function ilaenv(ISPEC, NAME, OPTS, N1, N2, N3, N4)
ILAENV
Definition: ilaenv.f:162
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:60
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:53
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 sgeqrf(M, N, A, LDA, TAU, WORK, LWORK, INFO)
SGEQRF
Definition: sgeqrf.f:146
subroutine sgelqf(M, N, A, LDA, TAU, WORK, LWORK, INFO)
SGELQF
Definition: sgelqf.f:143
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 sormlq(SIDE, TRANS, M, N, K, A, LDA, TAU, C, LDC, WORK, LWORK, INFO)
SORMLQ
Definition: sormlq.f:168
real function slamch(CMACH)
SLAMCH
Definition: slamch.f:68
Here is the call graph for this function:
Here is the caller graph for this function: