 LAPACK  3.10.1 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.```

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,
218  \$ sormqr, strtrs, xerbla
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 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: