LAPACK  3.8.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 )
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.```
Date
June 2017

Definition at line 162 of file sgetsls.f.

162 *
163 * -- LAPACK driver routine (version 3.7.1) --
164 * -- LAPACK is a software package provided by Univ. of Tennessee, --
165 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
166 * June 2017
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( 1 )
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,
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( 1 ) )
237  CALL sgemqr( 'L', trans, m, nrhs, n, a, lda, tq,
238  \$ tszo, b, ldb, workq, -1, info2 )
239  lwo = max( lwo, int( workq( 1 ) ) )
240  CALL sgeqr( m, n, a, lda, tq, -2, workq, -2, info2 )
241  tszm = int( tq( 1 ) )
242  lwm = int( workq( 1 ) )
243  CALL sgemqr( 'L', trans, m, nrhs, n, a, lda, tq,
244  \$ tszm, b, ldb, workq, -1, info2 )
245  lwm = max( lwm, int( workq( 1 ) ) )
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( 1 ) )
252  CALL sgemlq( 'L', trans, n, nrhs, m, a, lda, tq,
253  \$ tszo, b, ldb, workq, -1, info2 )
254  lwo = max( lwo, int( workq( 1 ) ) )
255  CALL sgelq( m, n, a, lda, tq, -2, workq, -2, info2 )
256  tszm = int( tq( 1 ) )
257  lwm = int( workq( 1 ) )
258  CALL sgemlq( 'L', trans, n, nrhs, m, a, lda, tq,
259  \$ tszo, b, ldb, workq, -1, info2 )
260  lwm = max( lwm, int( workq( 1 ) ) )
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 *
integer function ilaenv(ISPEC, NAME, OPTS, N1, N2, N3, N4)
ILAENV
Definition: tstiee.f:83
subroutine sgemlq(SIDE, TRANS, M, N, K, A, LDA, T, TSIZE, C, LDC, WORK, LWORK, INFO)
Definition: sgemlq.f:169
subroutine sgelq(M, N, A, LDA, T, TSIZE, WORK, LWORK, INFO)
Definition: sgelq.f:161
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 sgeqr(M, N, A, LDA, T, TSIZE, WORK, LWORK, INFO)
Definition: sgeqr.f:162
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