172 SUBROUTINE sgelss( M, N, NRHS, A, LDA, B, LDB, S, RCOND, RANK,
173 $ work, lwork, info )
181 INTEGER info, lda, ldb, lwork, m, n, nrhs, rank
185 REAL a( lda, * ),
b( ldb, * ), s( * ), work( * )
192 parameter( zero = 0.0e+0, one = 1.0e+0 )
196 INTEGER bdspac, bl, chunk, i, iascl, ibscl, ie, il,
197 $ itau, itaup, itauq, iwork, ldwork, maxmn,
198 $ maxwrk, minmn, minwrk, mm, mnthr
199 INTEGER lwork_sgeqrf, lwork_sormqr, lwork_sgebrd,
200 $ lwork_sormbr, lwork_sorgbr, lwork_sormlq
201 REAL anrm, bignum, bnrm, eps, sfmin, smlnum, thr
226 lquery = ( lwork.EQ.-1 )
229 ELSE IF( n.LT.0 )
THEN
231 ELSE IF( nrhs.LT.0 )
THEN
233 ELSE IF( lda.LT.max( 1, m ) )
THEN
235 ELSE IF( ldb.LT.max( 1, maxmn ) )
THEN
249 IF( minmn.GT.0 )
THEN
251 mnthr =
ilaenv( 6,
'SGELSS',
' ', m, n, nrhs, -1 )
252 IF( m.GE.n .AND. m.GE.mnthr )
THEN
258 CALL
sgeqrf( m, n, a, lda, dum(1), dum(1), -1, info )
261 CALL
sormqr(
'L',
'T', m, nrhs, n, a, lda, dum(1),
b,
262 $ ldb, dum(1), -1, info )
265 maxwrk = max( maxwrk, n + lwork_sgeqrf )
266 maxwrk = max( maxwrk, n + lwork_sormqr )
274 bdspac = max( 1, 5*n )
276 CALL
sgebrd( mm, n, a, lda, s, dum(1), dum(1),
277 $ dum(1), dum(1), -1, info )
280 CALL
sormbr(
'Q',
'L',
'T', mm, nrhs, n, a, lda, dum(1),
281 $
b, ldb, dum(1), -1, info )
284 CALL
sorgbr(
'P', n, n, n, a, lda, dum(1),
288 maxwrk = max( maxwrk, 3*n + lwork_sgebrd )
289 maxwrk = max( maxwrk, 3*n + lwork_sormbr )
290 maxwrk = max( maxwrk, 3*n + lwork_sorgbr )
291 maxwrk = max( maxwrk, bdspac )
292 maxwrk = max( maxwrk, n*nrhs )
293 minwrk = max( 3*n + mm, 3*n + nrhs, bdspac )
294 maxwrk = max( minwrk, maxwrk )
300 bdspac = max( 1, 5*m )
301 minwrk = max( 3*m+nrhs, 3*m+n, bdspac )
302 IF( n.GE.mnthr )
THEN
308 CALL
sgebrd( m, m, a, lda, s, dum(1), dum(1),
309 $ dum(1), dum(1), -1, info )
312 CALL
sormbr(
'Q',
'L',
'T', m, nrhs, n, a, lda,
313 $ dum(1),
b, ldb, dum(1), -1, info )
316 CALL
sorgbr(
'P', m, m, m, a, lda, dum(1),
320 CALL
sormlq(
'L',
'T', n, nrhs, m, a, lda, dum(1),
321 $
b, ldb, dum(1), -1, info )
324 maxwrk = m + m*
ilaenv( 1,
'SGELQF',
' ', m, n, -1,
326 maxwrk = max( maxwrk, m*m + 4*m + lwork_sgebrd )
327 maxwrk = max( maxwrk, m*m + 4*m + lwork_sormbr )
328 maxwrk = max( maxwrk, m*m + 4*m + lwork_sorgbr )
329 maxwrk = max( maxwrk, m*m + m + bdspac )
331 maxwrk = max( maxwrk, m*m + m + m*nrhs )
333 maxwrk = max( maxwrk, m*m + 2*m )
335 maxwrk = max( maxwrk, m + lwork_sormlq )
341 CALL
sgebrd( m, n, a, lda, s, dum(1), dum(1),
342 $ dum(1), dum(1), -1, info )
345 CALL
sormbr(
'Q',
'L',
'T', m, nrhs, m, a, lda,
346 $ dum(1),
b, ldb, dum(1), -1, info )
349 CALL
sorgbr(
'P', m, n, m, a, lda, dum(1),
352 maxwrk = 3*m + lwork_sgebrd
353 maxwrk = max( maxwrk, 3*m + lwork_sormbr )
354 maxwrk = max( maxwrk, 3*m + lwork_sorgbr )
355 maxwrk = max( maxwrk, bdspac )
356 maxwrk = max( maxwrk, n*nrhs )
359 maxwrk = max( minwrk, maxwrk )
363 IF( lwork.LT.minwrk .AND. .NOT.lquery )
368 CALL
xerbla(
'SGELSS', -info )
370 ELSE IF( lquery )
THEN
376 IF( m.EQ.0 .OR. n.EQ.0 )
THEN
386 bignum = one / smlnum
387 CALL
slabad( smlnum, bignum )
391 anrm =
slange(
'M', m, n, a, lda, work )
393 IF( anrm.GT.zero .AND. anrm.LT.smlnum )
THEN
397 CALL
slascl(
'G', 0, 0, anrm, smlnum, m, n, a, lda, info )
399 ELSE IF( anrm.GT.bignum )
THEN
403 CALL
slascl(
'G', 0, 0, anrm, bignum, m, n, a, lda, info )
405 ELSE IF( anrm.EQ.zero )
THEN
409 CALL
slaset(
'F', max( m, n ), nrhs, zero, zero,
b, ldb )
410 CALL
slaset(
'F', minmn, 1, zero, zero, s, 1 )
417 bnrm =
slange(
'M', m, nrhs,
b, ldb, work )
419 IF( bnrm.GT.zero .AND. bnrm.LT.smlnum )
THEN
423 CALL
slascl(
'G', 0, 0, bnrm, smlnum, m, nrhs,
b, ldb, info )
425 ELSE IF( bnrm.GT.bignum )
THEN
429 CALL
slascl(
'G', 0, 0, bnrm, bignum, m, nrhs,
b, ldb, info )
440 IF( m.GE.mnthr )
THEN
451 CALL
sgeqrf( m, n, a, lda, work( itau ), work( iwork ),
452 $ lwork-iwork+1, info )
457 CALL
sormqr(
'L',
'T', m, nrhs, n, a, lda, work( itau ),
b,
458 $ ldb, work( iwork ), lwork-iwork+1, info )
463 $ CALL
slaset(
'L', n-1, n-1, zero, zero, a( 2, 1 ), lda )
474 CALL
sgebrd( mm, n, a, lda, s, work( ie ), work( itauq ),
475 $ work( itaup ), work( iwork ), lwork-iwork+1,
481 CALL
sormbr(
'Q',
'L',
'T', mm, nrhs, n, a, lda, work( itauq ),
482 $
b, ldb, work( iwork ), lwork-iwork+1, info )
487 CALL
sorgbr(
'P', n, n, n, a, lda, work( itaup ),
488 $ work( iwork ), lwork-iwork+1, info )
496 CALL
sbdsqr(
'U', n, n, 0, nrhs, s, work( ie ), a, lda, dum,
497 $ 1,
b, ldb, work( iwork ), info )
503 thr = max( rcond*s( 1 ), sfmin )
505 $ thr = max( eps*s( 1 ), sfmin )
508 IF( s( i ).GT.thr )
THEN
509 CALL
srscl( nrhs, s( i ),
b( i, 1 ), ldb )
512 CALL
slaset(
'F', 1, nrhs, zero, zero,
b( i, 1 ), ldb )
519 IF( lwork.GE.ldb*nrhs .AND. nrhs.GT.1 )
THEN
520 CALL
sgemm(
'T',
'N', n, nrhs, n, one, a, lda,
b, ldb, zero,
522 CALL
slacpy(
'G', n, nrhs, work, ldb,
b, ldb )
523 ELSE IF( nrhs.GT.1 )
THEN
525 DO 20 i = 1, nrhs, chunk
526 bl = min( nrhs-i+1, chunk )
527 CALL
sgemm(
'T',
'N', n, bl, n, one, a, lda,
b( 1, i ),
528 $ ldb, zero, work, n )
529 CALL
slacpy(
'G', n, bl, work, n,
b( 1, i ), ldb )
532 CALL
sgemv(
'T', n, n, one, a, lda,
b, 1, zero, work, 1 )
533 CALL
scopy( n, work, 1,
b, 1 )
536 ELSE IF( n.GE.mnthr .AND. lwork.GE.4*m+m*m+
537 $ max( m, 2*m-4, nrhs, n-3*m ) )
THEN
543 IF( lwork.GE.max( 4*m+m*lda+max( m, 2*m-4, nrhs, n-3*m ),
544 $ m*lda+m+m*nrhs ) )ldwork = lda
551 CALL
sgelqf( m, n, a, lda, work( itau ), work( iwork ),
552 $ lwork-iwork+1, info )
557 CALL
slacpy(
'L', m, m, a, lda, work( il ), ldwork )
558 CALL
slaset(
'U', m-1, m-1, zero, zero, work( il+ldwork ),
568 CALL
sgebrd( m, m, work( il ), ldwork, s, work( ie ),
569 $ work( itauq ), work( itaup ), work( iwork ),
570 $ lwork-iwork+1, info )
575 CALL
sormbr(
'Q',
'L',
'T', m, nrhs, m, work( il ), ldwork,
576 $ work( itauq ),
b, ldb, work( iwork ),
577 $ lwork-iwork+1, info )
582 CALL
sorgbr(
'P', m, m, m, work( il ), ldwork, work( itaup ),
583 $ work( iwork ), lwork-iwork+1, info )
591 CALL
sbdsqr(
'U', m, m, 0, nrhs, s, work( ie ), work( il ),
592 $ ldwork, a, lda,
b, ldb, work( iwork ), info )
598 thr = max( rcond*s( 1 ), sfmin )
600 $ thr = max( eps*s( 1 ), sfmin )
603 IF( s( i ).GT.thr )
THEN
604 CALL
srscl( nrhs, s( i ),
b( i, 1 ), ldb )
607 CALL
slaset(
'F', 1, nrhs, zero, zero,
b( i, 1 ), ldb )
615 IF( lwork.GE.ldb*nrhs+iwork-1 .AND. nrhs.GT.1 )
THEN
616 CALL
sgemm(
'T',
'N', m, nrhs, m, one, work( il ), ldwork,
617 $
b, ldb, zero, work( iwork ), ldb )
618 CALL
slacpy(
'G', m, nrhs, work( iwork ), ldb,
b, ldb )
619 ELSE IF( nrhs.GT.1 )
THEN
620 chunk = ( lwork-iwork+1 ) / m
621 DO 40 i = 1, nrhs, chunk
622 bl = min( nrhs-i+1, chunk )
623 CALL
sgemm(
'T',
'N', m, bl, m, one, work( il ), ldwork,
624 $
b( 1, i ), ldb, zero, work( iwork ), m )
625 CALL
slacpy(
'G', m, bl, work( iwork ), m,
b( 1, i ),
629 CALL
sgemv(
'T', m, m, one, work( il ), ldwork,
b( 1, 1 ),
630 $ 1, zero, work( iwork ), 1 )
631 CALL
scopy( m, work( iwork ), 1,
b( 1, 1 ), 1 )
636 CALL
slaset(
'F', n-m, nrhs, zero, zero,
b( m+1, 1 ), ldb )
642 CALL
sormlq(
'L',
'T', n, nrhs, m, a, lda, work( itau ),
b,
643 $ ldb, work( iwork ), lwork-iwork+1, info )
657 CALL
sgebrd( m, n, a, lda, s, work( ie ), work( itauq ),
658 $ work( itaup ), work( iwork ), lwork-iwork+1,
664 CALL
sormbr(
'Q',
'L',
'T', m, nrhs, n, a, lda, work( itauq ),
665 $
b, ldb, work( iwork ), lwork-iwork+1, info )
670 CALL
sorgbr(
'P', m, n, m, a, lda, work( itaup ),
671 $ work( iwork ), lwork-iwork+1, info )
679 CALL
sbdsqr(
'L', m, n, 0, nrhs, s, work( ie ), a, lda, dum,
680 $ 1,
b, ldb, work( iwork ), info )
686 thr = max( rcond*s( 1 ), sfmin )
688 $ thr = max( eps*s( 1 ), sfmin )
691 IF( s( i ).GT.thr )
THEN
692 CALL
srscl( nrhs, s( i ),
b( i, 1 ), ldb )
695 CALL
slaset(
'F', 1, nrhs, zero, zero,
b( i, 1 ), ldb )
702 IF( lwork.GE.ldb*nrhs .AND. nrhs.GT.1 )
THEN
703 CALL
sgemm(
'T',
'N', n, nrhs, m, one, a, lda,
b, ldb, zero,
705 CALL
slacpy(
'F', n, nrhs, work, ldb,
b, ldb )
706 ELSE IF( nrhs.GT.1 )
THEN
708 DO 60 i = 1, nrhs, chunk
709 bl = min( nrhs-i+1, chunk )
710 CALL
sgemm(
'T',
'N', n, bl, m, one, a, lda,
b( 1, i ),
711 $ ldb, zero, work, n )
712 CALL
slacpy(
'F', n, bl, work, n,
b( 1, i ), ldb )
715 CALL
sgemv(
'T', m, n, one, a, lda,
b, 1, zero, work, 1 )
716 CALL
scopy( n, work, 1,
b, 1 )
722 IF( iascl.EQ.1 )
THEN
723 CALL
slascl(
'G', 0, 0, anrm, smlnum, n, nrhs,
b, ldb, info )
724 CALL
slascl(
'G', 0, 0, smlnum, anrm, minmn, 1, s, minmn,
726 ELSE IF( iascl.EQ.2 )
THEN
727 CALL
slascl(
'G', 0, 0, anrm, bignum, n, nrhs,
b, ldb, info )
728 CALL
slascl(
'G', 0, 0, bignum, anrm, minmn, 1, s, minmn,
731 IF( ibscl.EQ.1 )
THEN
732 CALL
slascl(
'G', 0, 0, smlnum, bnrm, n, nrhs,
b, ldb, info )
733 ELSE IF( ibscl.EQ.2 )
THEN
734 CALL
slascl(
'G', 0, 0, bignum, bnrm, n, nrhs,
b, ldb, info )
subroutine sgebrd(M, N, A, LDA, D, E, TAUQ, TAUP, WORK, LWORK, INFO)
SGEBRD
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...
subroutine sormlq(SIDE, TRANS, M, N, K, A, LDA, TAU, C, LDC, WORK, LWORK, INFO)
SORMLQ
subroutine sormqr(SIDE, TRANS, M, N, K, A, LDA, TAU, C, LDC, WORK, LWORK, INFO)
SORMQR
REAL function slamch(CMACH)
SLAMCH
subroutine sorgbr(VECT, M, N, K, A, LDA, TAU, WORK, LWORK, INFO)
SORGBR
subroutine sbdsqr(UPLO, N, NCVT, NRU, NCC, D, E, VT, LDVT, U, LDU, C, LDC, WORK, INFO)
SBDSQR
subroutine sgelss(M, N, NRHS, A, LDA, B, LDB, S, RCOND, RANK, WORK, LWORK, INFO)
SGELSS solves overdetermined or underdetermined systems for GE matrices
subroutine xerbla(SRNAME, INFO)
XERBLA
set ue cd $ADTTMP cat<< EOF > tmp f Program LinearEquations Implicit none Real b(3) integer i
subroutine scopy(N, SX, INCX, SY, INCY)
SCOPY
subroutine sgelqf(M, N, A, LDA, TAU, WORK, LWORK, INFO)
SGELQF
subroutine sgemv(TRANS, M, N, ALPHA, A, LDA, X, INCX, BETA, Y, INCY)
SGEMV
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.
subroutine slacpy(UPLO, M, N, A, LDA, B, LDB)
SLACPY copies all or part of one two-dimensional array to another.
INTEGER function ilaenv(ISPEC, NAME, OPTS, N1, N2, N3, N4)
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 ...
subroutine sgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
SGEMM
subroutine srscl(N, SA, SX, INCX)
SRSCL multiplies a vector by the reciprocal of a real scalar.
subroutine sgeqrf(M, N, A, LDA, TAU, WORK, LWORK, INFO)
SGEQRF
subroutine slabad(SMALL, LARGE)
SLABAD
subroutine sormbr(VECT, SIDE, TRANS, M, N, K, A, LDA, TAU, C, LDC, WORK, LWORK, INFO)
SORMBR