LAPACK  3.9.0
LAPACK: Linear Algebra PACKage

◆ sgetsls()

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 
)

SGETSLS

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
June 2017

Definition at line 164 of file sgetsls.f.

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