LAPACK  3.7.0
LAPACK: Linear Algebra PACKage
subroutine sgetsls ( 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 
)
Purpose:

SGETSLS solves overdetermined or underdetermined real linear systems involving an M-by-N matrix A, using a tall skinny QR or short wide 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 undetermined 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,
          A is overwritten by details of its QR or LQ
          factorization as returned by SGEQR or SGELQ.
[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.
          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.
[in]LDB
          LDB is INTEGER
          The leading dimension of the array B. LDB >= MAX(1,M,N).
[out]WORK
          (workspace) REAL array, dimension (MAX(1,LWORK))
          On exit, if INFO = 0, WORK(1) contains optimal (or either minimal
          or optimal, if query was assumed) LWORK.
          See LWORK for details.
[in]LWORK
          LWORK is INTEGER
          The dimension of the array WORK.
          If LWORK = -1 or -2, then a workspace query is assumed.
          If LWORK = -1, the routine calculates optimal size of WORK for the
          optimal performance and returns this value in WORK(1).
          If LWORK = -2, the routine calculates minimal size of WORK and 
          returns this value in WORK(1).
[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 162 of file sgetsls.f.

162 *
163 * -- LAPACK driver routine (version 3.7.0) --
164 * -- LAPACK is a software package provided by Univ. of Tennessee, --
165 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
166 * December 2016
167 *
168 * .. Scalar Arguments ..
169  CHARACTER trans
170  INTEGER info, lda, ldb, lwork, m, n, nrhs
171 * ..
172 * .. Array Arguments ..
173  REAL a( lda, * ), b( ldb, * ), work( * )
174 *
175 * ..
176 *
177 * =====================================================================
178 *
179 * .. Parameters ..
180  REAL zero, one
181  parameter ( zero = 0.0e0, one = 1.0e0 )
182 * ..
183 * .. Local Scalars ..
184  LOGICAL lquery, tran
185  INTEGER i, iascl, ibscl, j, minmn, maxmn, brow,
186  $ scllen, mnk, tszo, tszm, lwo, lwm, lw1, lw2,
187  $ wsizeo, wsizem, info2
188  REAL anrm, bignum, bnrm, smlnum, tq( 5 ), workq
189 * ..
190 * .. External Functions ..
191  LOGICAL lsame
192  INTEGER ilaenv
193  REAL slamch, slange
194  EXTERNAL lsame, ilaenv, slabad, slamch, slange
195 * ..
196 * .. External Subroutines ..
197  EXTERNAL sgeqr, sgemqr, slascl, slaset,
198  $ strtrs, xerbla, sgelq, sgemlq
199 * ..
200 * .. Intrinsic Functions ..
201  INTRINSIC REAL, max, min, int
202 * ..
203 * .. Executable Statements ..
204 *
205 * Test the input arguments.
206 *
207  info = 0
208  minmn = min( m, n )
209  maxmn = max( m, n )
210  mnk = max( minmn, nrhs )
211  tran = lsame( trans, 'T' )
212 *
213  lquery = ( lwork.EQ.-1 .OR. lwork.EQ.-2 )
214  IF( .NOT.( lsame( trans, 'N' ) .OR.
215  $ lsame( trans, 'T' ) ) ) THEN
216  info = -1
217  ELSE IF( m.LT.0 ) THEN
218  info = -2
219  ELSE IF( n.LT.0 ) THEN
220  info = -3
221  ELSE IF( nrhs.LT.0 ) THEN
222  info = -4
223  ELSE IF( lda.LT.max( 1, m ) ) THEN
224  info = -6
225  ELSE IF( ldb.LT.max( 1, m, n ) ) THEN
226  info = -8
227  END IF
228 *
229  IF( info.EQ.0 ) THEN
230 *
231 * Determine the block size and minimum LWORK
232 *
233  IF( m.GE.n ) THEN
234  CALL sgeqr( m, n, a, lda, tq, -1, workq, -1, info2 )
235  tszo = int( tq( 1 ) )
236  lwo = int( workq )
237  CALL sgemqr( 'L', trans, m, nrhs, n, a, lda, tq,
238  $ tszo, b, ldb, workq, -1, info2 )
239  lwo = max( lwo, int( workq ) )
240  CALL sgeqr( m, n, a, lda, tq, -2, workq, -2, info2 )
241  tszm = int( tq( 1 ) )
242  lwm = int( workq )
243  CALL sgemqr( 'L', trans, m, nrhs, n, a, lda, tq,
244  $ tszm, b, ldb, workq, -1, info2 )
245  lwm = max( lwm, int( workq ) )
246  wsizeo = tszo + lwo
247  wsizem = tszm + lwm
248  ELSE
249  CALL sgelq( m, n, a, lda, tq, -1, workq, -1, info2 )
250  tszo = int( tq( 1 ) )
251  lwo = int( workq )
252  CALL sgemlq( 'L', trans, n, nrhs, m, a, lda, tq,
253  $ tszo, b, ldb, workq, -1, info2 )
254  lwo = max( lwo, int( workq ) )
255  CALL sgelq( m, n, a, lda, tq, -2, workq, -2, info2 )
256  tszm = int( tq( 1 ) )
257  lwm = int( workq )
258  CALL sgemlq( 'L', trans, n, nrhs, m, a, lda, tq,
259  $ tszo, b, ldb, workq, -1, info2 )
260  lwm = max( lwm, int( workq ) )
261  wsizeo = tszo + lwo
262  wsizem = tszm + lwm
263  END IF
264 *
265  IF( ( lwork.LT.wsizem ).AND.( .NOT.lquery ) ) THEN
266  info = -10
267  END IF
268 *
269  END IF
270 *
271  IF( info.NE.0 ) THEN
272  CALL xerbla( 'SGETSLS', -info )
273  work( 1 ) = REAL( wsizeo )
274  RETURN
275  END IF
276  IF( lquery ) THEN
277  IF( lwork.EQ.-1 ) work( 1 ) = REAL( wsizeo )
278  IF( lwork.EQ.-2 ) work( 1 ) = REAL( wsizem )
279  RETURN
280  END IF
281  IF( lwork.LT.wsizeo ) THEN
282  lw1 = tszm
283  lw2 = lwm
284  ELSE
285  lw1 = tszo
286  lw2 = lwo
287  END IF
288 *
289 * Quick return if possible
290 *
291  IF( min( m, n, nrhs ).EQ.0 ) THEN
292  CALL slaset( 'FULL', max( m, n ), nrhs, zero, zero,
293  $ 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, work )
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', maxmn, nrhs, zero, zero, b, ldb )
324  GO TO 50
325  END IF
326 *
327  brow = m
328  IF ( tran ) THEN
329  brow = n
330  END IF
331  bnrm = slange( 'M', brow, nrhs, b, ldb, work )
332  ibscl = 0
333  IF( bnrm.GT.zero .AND. bnrm.LT.smlnum ) THEN
334 *
335 * Scale matrix norm up to SMLNUM
336 *
337  CALL slascl( 'G', 0, 0, bnrm, smlnum, brow, nrhs, b, ldb,
338  $ info )
339  ibscl = 1
340  ELSE IF( bnrm.GT.bignum ) THEN
341 *
342 * Scale matrix norm down to BIGNUM
343 *
344  CALL slascl( 'G', 0, 0, bnrm, bignum, brow, nrhs, b, ldb,
345  $ info )
346  ibscl = 2
347  END IF
348 *
349  IF ( m.GE.n ) THEN
350 *
351 * compute QR factorization of A
352 *
353  CALL sgeqr( m, n, a, lda, work( lw2+1 ), lw1,
354  $ work( 1 ), lw2, info )
355  IF ( .NOT.tran ) THEN
356 *
357 * Least-Squares Problem min || A * X - B ||
358 *
359 * B(1:M,1:NRHS) := Q**T * B(1:M,1:NRHS)
360 *
361  CALL sgemqr( 'L' , 'T', m, nrhs, n, a, lda,
362  $ work( lw2+1 ), lw1, b, ldb, work( 1 ), lw2,
363  $ info )
364 *
365 * B(1:N,1:NRHS) := inv(R) * B(1:N,1:NRHS)
366 *
367  CALL strtrs( 'U', 'N', 'N', n, nrhs,
368  $ a, lda, b, ldb, info )
369  IF( info.GT.0 ) THEN
370  RETURN
371  END IF
372  scllen = n
373  ELSE
374 *
375 * Overdetermined system of equations A**T * X = B
376 *
377 * B(1:N,1:NRHS) := inv(R**T) * B(1:N,1:NRHS)
378 *
379  CALL strtrs( 'U', 'T', 'N', n, nrhs,
380  $ a, lda, b, ldb, info )
381 *
382  IF( info.GT.0 ) THEN
383  RETURN
384  END IF
385 *
386 * B(N+1:M,1:NRHS) = ZERO
387 *
388  DO 20 j = 1, nrhs
389  DO 10 i = n + 1, m
390  b( i, j ) = zero
391  10 CONTINUE
392  20 CONTINUE
393 *
394 * B(1:M,1:NRHS) := Q(1:N,:) * B(1:N,1:NRHS)
395 *
396  CALL sgemqr( 'L', 'N', m, nrhs, n, a, lda,
397  $ work( lw2+1 ), lw1, b, ldb, work( 1 ), lw2,
398  $ info )
399 *
400  scllen = m
401 *
402  END IF
403 *
404  ELSE
405 *
406 * Compute LQ factorization of A
407 *
408  CALL sgelq( m, n, a, lda, work( lw2+1 ), lw1,
409  $ work( 1 ), lw2, info )
410 *
411 * workspace at least M, optimally M*NB.
412 *
413  IF( .NOT.tran ) THEN
414 *
415 * underdetermined system of equations A * X = B
416 *
417 * B(1:M,1:NRHS) := inv(L) * B(1:M,1:NRHS)
418 *
419  CALL strtrs( 'L', 'N', 'N', m, nrhs,
420  $ a, lda, b, ldb, info )
421 *
422  IF( info.GT.0 ) THEN
423  RETURN
424  END IF
425 *
426 * B(M+1:N,1:NRHS) = 0
427 *
428  DO 40 j = 1, nrhs
429  DO 30 i = m + 1, n
430  b( i, j ) = zero
431  30 CONTINUE
432  40 CONTINUE
433 *
434 * B(1:N,1:NRHS) := Q(1:N,:)**T * B(1:M,1:NRHS)
435 *
436  CALL sgemlq( 'L', 'T', n, nrhs, m, a, lda,
437  $ work( lw2+1 ), lw1, b, ldb, work( 1 ), lw2,
438  $ info )
439 *
440 * workspace at least NRHS, optimally NRHS*NB
441 *
442  scllen = n
443 *
444  ELSE
445 *
446 * overdetermined system min || A**T * X - B ||
447 *
448 * B(1:N,1:NRHS) := Q * B(1:N,1:NRHS)
449 *
450  CALL sgemlq( 'L', 'N', n, nrhs, m, a, lda,
451  $ work( lw2+1 ), lw1, b, ldb, work( 1 ), lw2,
452  $ info )
453 *
454 * workspace at least NRHS, optimally NRHS*NB
455 *
456 * B(1:M,1:NRHS) := inv(L**T) * B(1:M,1:NRHS)
457 *
458  CALL strtrs( 'Lower', 'Transpose', 'Non-unit', m, nrhs,
459  $ a, lda, b, ldb, info )
460 *
461  IF( info.GT.0 ) THEN
462  RETURN
463  END IF
464 *
465  scllen = m
466 *
467  END IF
468 *
469  END IF
470 *
471 * Undo scaling
472 *
473  IF( iascl.EQ.1 ) THEN
474  CALL slascl( 'G', 0, 0, anrm, smlnum, scllen, nrhs, b, ldb,
475  $ info )
476  ELSE IF( iascl.EQ.2 ) THEN
477  CALL slascl( 'G', 0, 0, anrm, bignum, scllen, nrhs, b, ldb,
478  $ info )
479  END IF
480  IF( ibscl.EQ.1 ) THEN
481  CALL slascl( 'G', 0, 0, smlnum, bnrm, scllen, nrhs, b, ldb,
482  $ info )
483  ELSE IF( ibscl.EQ.2 ) THEN
484  CALL slascl( 'G', 0, 0, bignum, bnrm, scllen, nrhs, b, ldb,
485  $ info )
486  END IF
487 *
488  50 CONTINUE
489  work( 1 ) = REAL( tszo + lwo )
490  RETURN
491 *
492 * End of SGETSLS
493 *
ilaenv
integer function ilaenv(ISPEC, NAME, OPTS, N1, N2, N3, N4)
ILAENV
Definition: tstiee.f:83
strtrs
subroutine strtrs(UPLO, TRANS, DIAG, N, NRHS, A, LDA, B, LDB, INFO)
STRTRS
Definition: strtrs.f:142
slange
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
xerbla
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
sgeqr
subroutine sgeqr(M, N, A, LDA, T, TSIZE, WORK, LWORK, INFO)
Definition: sgeqr.f:162
slaset
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
slascl
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
lsame
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:55
slamch
real function slamch(CMACH)
SLAMCH
Definition: slamch.f:69
slabad
subroutine slabad(SMALL, LARGE)
SLABAD
Definition: slabad.f:76
sgemqr
subroutine sgemqr(SIDE, TRANS, M, N, K, A, LDA, T, TSIZE, C, LDC, WORK, LWORK, INFO)
Definition: sgemqr.f:171
sgemlq
subroutine sgemlq(SIDE, TRANS, M, N, K, A, LDA, T, TSIZE, C, LDC, WORK, LWORK, INFO)
Definition: sgemlq.f:169
sgelq
subroutine sgelq(M, N, A, LDA, T, TSIZE, WORK, LWORK, INFO)
Definition: sgelq.f:161

Here is the call graph for this function:

Here is the caller graph for this function: