LAPACK 3.12.0
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches

◆ 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.

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, SROUNDUP_LWORK
215* ..
216* .. External Subroutines ..
217 EXTERNAL sgelqf, sgeqrf, slascl, slaset, sormlq,
219* ..
220* .. Intrinsic Functions ..
221 INTRINSIC max, min
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 ) = sroundup_lwork( 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*
299* Scale A, B if max element outside range [SMLNUM,BIGNUM]
300*
301 anrm = slange( 'M', m, n, a, lda, rwork )
302 iascl = 0
303 IF( anrm.GT.zero .AND. anrm.LT.smlnum ) THEN
304*
305* Scale matrix norm up to SMLNUM
306*
307 CALL slascl( 'G', 0, 0, anrm, smlnum, m, n, a, lda, info )
308 iascl = 1
309 ELSE IF( anrm.GT.bignum ) THEN
310*
311* Scale matrix norm down to BIGNUM
312*
313 CALL slascl( 'G', 0, 0, anrm, bignum, m, n, a, lda, info )
314 iascl = 2
315 ELSE IF( anrm.EQ.zero ) THEN
316*
317* Matrix all zero. Return zero solution.
318*
319 CALL slaset( 'F', max( m, n ), nrhs, zero, zero, b, ldb )
320 GO TO 50
321 END IF
322*
323 brow = m
324 IF( tpsd )
325 $ brow = n
326 bnrm = slange( 'M', brow, nrhs, b, ldb, rwork )
327 ibscl = 0
328 IF( bnrm.GT.zero .AND. bnrm.LT.smlnum ) THEN
329*
330* Scale matrix norm up to SMLNUM
331*
332 CALL slascl( 'G', 0, 0, bnrm, smlnum, brow, nrhs, b, ldb,
333 $ info )
334 ibscl = 1
335 ELSE IF( bnrm.GT.bignum ) THEN
336*
337* Scale matrix norm down to BIGNUM
338*
339 CALL slascl( 'G', 0, 0, bnrm, bignum, brow, nrhs, b, ldb,
340 $ info )
341 ibscl = 2
342 END IF
343*
344 IF( m.GE.n ) THEN
345*
346* compute QR factorization of A
347*
348 CALL sgeqrf( m, n, a, lda, work( 1 ), work( mn+1 ), lwork-mn,
349 $ info )
350*
351* workspace at least N, optimally N*NB
352*
353 IF( .NOT.tpsd ) THEN
354*
355* Least-Squares Problem min || A * X - B ||
356*
357* B(1:M,1:NRHS) := Q**T * B(1:M,1:NRHS)
358*
359 CALL sormqr( 'Left', 'Transpose', m, nrhs, n, a, lda,
360 $ work( 1 ), b, ldb, work( mn+1 ), lwork-mn,
361 $ info )
362*
363* workspace at least NRHS, optimally NRHS*NB
364*
365* B(1:N,1:NRHS) := inv(R) * B(1:N,1:NRHS)
366*
367 CALL strtrs( 'Upper', 'No transpose', 'Non-unit', n, nrhs,
368 $ a, lda, b, ldb, info )
369*
370 IF( info.GT.0 ) THEN
371 RETURN
372 END IF
373*
374 scllen = n
375*
376 ELSE
377*
378* Underdetermined system of equations A**T * X = B
379*
380* B(1:N,1:NRHS) := inv(R**T) * B(1:N,1:NRHS)
381*
382 CALL strtrs( 'Upper', 'Transpose', 'Non-unit', n, nrhs,
383 $ a, lda, b, ldb, info )
384*
385 IF( info.GT.0 ) THEN
386 RETURN
387 END IF
388*
389* B(N+1:M,1:NRHS) = ZERO
390*
391 DO 20 j = 1, nrhs
392 DO 10 i = n + 1, m
393 b( i, j ) = zero
394 10 CONTINUE
395 20 CONTINUE
396*
397* B(1:M,1:NRHS) := Q(1:N,:) * B(1:N,1:NRHS)
398*
399 CALL sormqr( 'Left', 'No transpose', m, nrhs, n, a, lda,
400 $ work( 1 ), b, ldb, work( mn+1 ), lwork-mn,
401 $ info )
402*
403* workspace at least NRHS, optimally NRHS*NB
404*
405 scllen = m
406*
407 END IF
408*
409 ELSE
410*
411* Compute LQ factorization of A
412*
413 CALL sgelqf( m, n, a, lda, work( 1 ), work( mn+1 ), lwork-mn,
414 $ info )
415*
416* workspace at least M, optimally M*NB.
417*
418 IF( .NOT.tpsd ) THEN
419*
420* underdetermined system of equations A * X = B
421*
422* B(1:M,1:NRHS) := inv(L) * B(1:M,1:NRHS)
423*
424 CALL strtrs( 'Lower', 'No transpose', 'Non-unit', m, nrhs,
425 $ a, lda, b, ldb, info )
426*
427 IF( info.GT.0 ) THEN
428 RETURN
429 END IF
430*
431* B(M+1:N,1:NRHS) = 0
432*
433 DO 40 j = 1, nrhs
434 DO 30 i = m + 1, n
435 b( i, j ) = zero
436 30 CONTINUE
437 40 CONTINUE
438*
439* B(1:N,1:NRHS) := Q(1:N,:)**T * B(1:M,1:NRHS)
440*
441 CALL sormlq( 'Left', 'Transpose', n, nrhs, m, a, lda,
442 $ work( 1 ), b, ldb, work( mn+1 ), lwork-mn,
443 $ info )
444*
445* workspace at least NRHS, optimally NRHS*NB
446*
447 scllen = n
448*
449 ELSE
450*
451* overdetermined system min || A**T * X - B ||
452*
453* B(1:N,1:NRHS) := Q * B(1:N,1:NRHS)
454*
455 CALL sormlq( 'Left', 'No transpose', n, nrhs, m, a, lda,
456 $ work( 1 ), b, ldb, work( mn+1 ), lwork-mn,
457 $ info )
458*
459* workspace at least NRHS, optimally NRHS*NB
460*
461* B(1:M,1:NRHS) := inv(L**T) * B(1:M,1:NRHS)
462*
463 CALL strtrs( 'Lower', 'Transpose', 'Non-unit', m, nrhs,
464 $ a, lda, b, ldb, info )
465*
466 IF( info.GT.0 ) THEN
467 RETURN
468 END IF
469*
470 scllen = m
471*
472 END IF
473*
474 END IF
475*
476* Undo scaling
477*
478 IF( iascl.EQ.1 ) THEN
479 CALL slascl( 'G', 0, 0, anrm, smlnum, scllen, nrhs, b, ldb,
480 $ info )
481 ELSE IF( iascl.EQ.2 ) THEN
482 CALL slascl( 'G', 0, 0, anrm, bignum, scllen, nrhs, b, ldb,
483 $ info )
484 END IF
485 IF( ibscl.EQ.1 ) THEN
486 CALL slascl( 'G', 0, 0, smlnum, bnrm, scllen, nrhs, b, ldb,
487 $ info )
488 ELSE IF( ibscl.EQ.2 ) THEN
489 CALL slascl( 'G', 0, 0, bignum, bnrm, scllen, nrhs, b, ldb,
490 $ info )
491 END IF
492*
493 50 CONTINUE
494 work( 1 ) = sroundup_lwork( wsize )
495*
496 RETURN
497*
498* End of SGELS
499*
subroutine xerbla(srname, info)
Definition cblat2.f:3285
subroutine sgelqf(m, n, a, lda, tau, work, lwork, info)
SGELQF
Definition sgelqf.f:143
subroutine sgeqrf(m, n, a, lda, tau, work, lwork, info)
SGEQRF
Definition sgeqrf.f:146
integer function ilaenv(ispec, name, opts, n1, n2, n3, n4)
ILAENV
Definition ilaenv.f:162
real function slamch(cmach)
SLAMCH
Definition slamch.f:68
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 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
logical function lsame(ca, cb)
LSAME
Definition lsame.f:48
real function sroundup_lwork(lwork)
SROUNDUP_LWORK
subroutine strtrs(uplo, trans, diag, n, nrhs, a, lda, b, ldb, info)
STRTRS
Definition strtrs.f:140
subroutine sormlq(side, trans, m, n, k, a, lda, tau, c, ldc, work, lwork, info)
SORMLQ
Definition sormlq.f:168
subroutine sormqr(side, trans, m, n, k, a, lda, tau, c, ldc, work, lwork, info)
SORMQR
Definition sormqr.f:168
Here is the call graph for this function:
Here is the caller graph for this function: