LAPACK  3.8.0
LAPACK: Linear Algebra PACKage

◆ 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.
Date
December 2016

Definition at line 185 of file sgels.f.

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