178 INTEGER INFO, LDA, LDB, LWORK, M, N, NRHS, RANK
182 REAL A( LDA, * ), B( LDB, * ), S( * ), WORK( * )
189 parameter( zero = 0.0e+0, one = 1.0e+0 )
193 INTEGER BDSPAC, BL, CHUNK, I, IASCL, IBSCL, IE, IL,
194 $ ITAU, ITAUP, ITAUQ, IWORK, LDWORK, MAXMN,
195 $ MAXWRK, MINMN, MINWRK, MM, MNTHR
196 INTEGER LWORK_SGEQRF, LWORK_SORMQR, LWORK_SGEBRD,
197 $ LWORK_SORMBR, LWORK_SORGBR, LWORK_SORMLQ
198 REAL ANRM, BIGNUM, BNRM, EPS, SFMIN, SMLNUM, THR
223 lquery = ( lwork.EQ.-1 )
226 ELSE IF( n.LT.0 )
THEN
228 ELSE IF( nrhs.LT.0 )
THEN
230 ELSE IF( lda.LT.max( 1, m ) )
THEN
232 ELSE IF( ldb.LT.max( 1, maxmn ) )
THEN
246 IF( minmn.GT.0 )
THEN
248 mnthr =
ilaenv( 6,
'SGELSS',
' ', m, n, nrhs, -1 )
249 IF( m.GE.n .AND. m.GE.mnthr )
THEN
255 CALL sgeqrf( m, n, a, lda, dum(1), dum(1), -1, info )
258 CALL sormqr(
'L',
'T', m, nrhs, n, a, lda, dum(1), b,
259 $ ldb, dum(1), -1, info )
262 maxwrk = max( maxwrk, n + lwork_sgeqrf )
263 maxwrk = max( maxwrk, n + lwork_sormqr )
271 bdspac = max( 1, 5*n )
273 CALL sgebrd( mm, n, a, lda, s, dum(1), dum(1),
274 $ dum(1), dum(1), -1, info )
277 CALL sormbr(
'Q',
'L',
'T', mm, nrhs, n, a, lda, dum(1),
278 $ b, ldb, dum(1), -1, info )
281 CALL sorgbr(
'P', n, n, n, a, lda, dum(1),
285 maxwrk = max( maxwrk, 3*n + lwork_sgebrd )
286 maxwrk = max( maxwrk, 3*n + lwork_sormbr )
287 maxwrk = max( maxwrk, 3*n + lwork_sorgbr )
288 maxwrk = max( maxwrk, bdspac )
289 maxwrk = max( maxwrk, n*nrhs )
290 minwrk = max( 3*n + mm, 3*n + nrhs, bdspac )
291 maxwrk = max( minwrk, maxwrk )
297 bdspac = max( 1, 5*m )
298 minwrk = max( 3*m+nrhs, 3*m+n, bdspac )
299 IF( n.GE.mnthr )
THEN
305 CALL sgebrd( m, m, a, lda, s, dum(1), dum(1),
306 $ dum(1), dum(1), -1, info )
309 CALL sormbr(
'Q',
'L',
'T', m, nrhs, n, a, lda,
310 $ dum(1), b, ldb, dum(1), -1, info )
313 CALL sorgbr(
'P', m, m, m, a, lda, dum(1),
317 CALL sormlq(
'L',
'T', n, nrhs, m, a, lda, dum(1),
318 $ b, ldb, dum(1), -1, info )
321 maxwrk = m + m*
ilaenv( 1,
'SGELQF',
' ', m, n, -1,
323 maxwrk = max( maxwrk, m*m + 4*m + lwork_sgebrd )
324 maxwrk = max( maxwrk, m*m + 4*m + lwork_sormbr )
325 maxwrk = max( maxwrk, m*m + 4*m + lwork_sorgbr )
326 maxwrk = max( maxwrk, m*m + m + bdspac )
328 maxwrk = max( maxwrk, m*m + m + m*nrhs )
330 maxwrk = max( maxwrk, m*m + 2*m )
332 maxwrk = max( maxwrk, m + lwork_sormlq )
338 CALL sgebrd( m, n, a, lda, s, dum(1), dum(1),
339 $ dum(1), dum(1), -1, info )
342 CALL sormbr(
'Q',
'L',
'T', m, nrhs, m, a, lda,
343 $ dum(1), b, ldb, dum(1), -1, info )
346 CALL sorgbr(
'P', m, n, m, a, lda, dum(1),
349 maxwrk = 3*m + lwork_sgebrd
350 maxwrk = max( maxwrk, 3*m + lwork_sormbr )
351 maxwrk = max( maxwrk, 3*m + lwork_sorgbr )
352 maxwrk = max( maxwrk, bdspac )
353 maxwrk = max( maxwrk, n*nrhs )
356 maxwrk = max( minwrk, maxwrk )
360 IF( lwork.LT.minwrk .AND. .NOT.lquery )
365 CALL xerbla(
'SGELSS', -info )
367 ELSE IF( lquery )
THEN
373 IF( m.EQ.0 .OR. n.EQ.0 )
THEN
383 bignum = one / smlnum
384 CALL slabad( smlnum, bignum )
388 anrm =
slange(
'M', m, n, a, lda, work )
390 IF( anrm.GT.zero .AND. anrm.LT.smlnum )
THEN
394 CALL slascl(
'G', 0, 0, anrm, smlnum, m, n, a, lda, info )
396 ELSE IF( anrm.GT.bignum )
THEN
400 CALL slascl(
'G', 0, 0, anrm, bignum, m, n, a, lda, info )
402 ELSE IF( anrm.EQ.zero )
THEN
406 CALL slaset(
'F', max( m, n ), nrhs, zero, zero, b, ldb )
407 CALL slaset(
'F', minmn, 1, zero, zero, s, minmn )
414 bnrm =
slange(
'M', m, nrhs, b, ldb, work )
416 IF( bnrm.GT.zero .AND. bnrm.LT.smlnum )
THEN
420 CALL slascl(
'G', 0, 0, bnrm, smlnum, m, nrhs, b, ldb, info )
422 ELSE IF( bnrm.GT.bignum )
THEN
426 CALL slascl(
'G', 0, 0, bnrm, bignum, m, nrhs, b, ldb, info )
437 IF( m.GE.mnthr )
THEN
448 CALL sgeqrf( m, n, a, lda, work( itau ), work( iwork ),
449 $ lwork-iwork+1, info )
454 CALL sormqr(
'L',
'T', m, nrhs, n, a, lda, work( itau ), b,
455 $ ldb, work( iwork ), lwork-iwork+1, info )
460 $
CALL slaset(
'L', n-1, n-1, zero, zero, a( 2, 1 ), lda )
471 CALL sgebrd( mm, n, a, lda, s, work( ie ), work( itauq ),
472 $ work( itaup ), work( iwork ), lwork-iwork+1,
478 CALL sormbr(
'Q',
'L',
'T', mm, nrhs, n, a, lda, work( itauq ),
479 $ b, ldb, work( iwork ), lwork-iwork+1, info )
484 CALL sorgbr(
'P', n, n, n, a, lda, work( itaup ),
485 $ work( iwork ), lwork-iwork+1, info )
493 CALL sbdsqr(
'U', n, n, 0, nrhs, s, work( ie ), a, lda, dum,
494 $ 1, b, ldb, work( iwork ), info )
500 thr = max( rcond*s( 1 ), sfmin )
502 $ thr = max( eps*s( 1 ), sfmin )
505 IF( s( i ).GT.thr )
THEN
506 CALL srscl( nrhs, s( i ), b( i, 1 ), ldb )
509 CALL slaset(
'F', 1, nrhs, zero, zero, b( i, 1 ), ldb )
516 IF( lwork.GE.ldb*nrhs .AND. nrhs.GT.1 )
THEN
517 CALL sgemm(
'T',
'N', n, nrhs, n, one, a, lda, b, ldb, zero,
519 CALL slacpy(
'G', n, nrhs, work, ldb, b, ldb )
520 ELSE IF( nrhs.GT.1 )
THEN
522 DO 20 i = 1, nrhs, chunk
523 bl = min( nrhs-i+1, chunk )
524 CALL sgemm(
'T',
'N', n, bl, n, one, a, lda, b( 1, i ),
525 $ ldb, zero, work, n )
526 CALL slacpy(
'G', n, bl, work, n, b( 1, i ), ldb )
529 CALL sgemv(
'T', n, n, one, a, lda, b, 1, zero, work, 1 )
530 CALL scopy( n, work, 1, b, 1 )
533 ELSE IF( n.GE.mnthr .AND. lwork.GE.4*m+m*m+
534 $ max( m, 2*m-4, nrhs, n-3*m ) )
THEN
540 IF( lwork.GE.max( 4*m+m*lda+max( m, 2*m-4, nrhs, n-3*m ),
541 $ m*lda+m+m*nrhs ) )ldwork = lda
548 CALL sgelqf( m, n, a, lda, work( itau ), work( iwork ),
549 $ lwork-iwork+1, info )
554 CALL slacpy(
'L', m, m, a, lda, work( il ), ldwork )
555 CALL slaset(
'U', m-1, m-1, zero, zero, work( il+ldwork ),
565 CALL sgebrd( m, m, work( il ), ldwork, s, work( ie ),
566 $ work( itauq ), work( itaup ), work( iwork ),
567 $ lwork-iwork+1, info )
572 CALL sormbr(
'Q',
'L',
'T', m, nrhs, m, work( il ), ldwork,
573 $ work( itauq ), b, ldb, work( iwork ),
574 $ lwork-iwork+1, info )
579 CALL sorgbr(
'P', m, m, m, work( il ), ldwork, work( itaup ),
580 $ work( iwork ), lwork-iwork+1, info )
588 CALL sbdsqr(
'U', m, m, 0, nrhs, s, work( ie ), work( il ),
589 $ ldwork, a, lda, b, ldb, work( iwork ), info )
595 thr = max( rcond*s( 1 ), sfmin )
597 $ thr = max( eps*s( 1 ), sfmin )
600 IF( s( i ).GT.thr )
THEN
601 CALL srscl( nrhs, s( i ), b( i, 1 ), ldb )
604 CALL slaset(
'F', 1, nrhs, zero, zero, b( i, 1 ), ldb )
612 IF( lwork.GE.ldb*nrhs+iwork-1 .AND. nrhs.GT.1 )
THEN
613 CALL sgemm(
'T',
'N', m, nrhs, m, one, work( il ), ldwork,
614 $ b, ldb, zero, work( iwork ), ldb )
615 CALL slacpy(
'G', m, nrhs, work( iwork ), ldb, b, ldb )
616 ELSE IF( nrhs.GT.1 )
THEN
617 chunk = ( lwork-iwork+1 ) / m
618 DO 40 i = 1, nrhs, chunk
619 bl = min( nrhs-i+1, chunk )
620 CALL sgemm(
'T',
'N', m, bl, m, one, work( il ), ldwork,
621 $ b( 1, i ), ldb, zero, work( iwork ), m )
622 CALL slacpy(
'G', m, bl, work( iwork ), m, b( 1, i ),
626 CALL sgemv(
'T', m, m, one, work( il ), ldwork, b( 1, 1 ),
627 $ 1, zero, work( iwork ), 1 )
628 CALL scopy( m, work( iwork ), 1, b( 1, 1 ), 1 )
633 CALL slaset(
'F', n-m, nrhs, zero, zero, b( m+1, 1 ), ldb )
639 CALL sormlq(
'L',
'T', n, nrhs, m, a, lda, work( itau ), b,
640 $ ldb, work( iwork ), lwork-iwork+1, info )
654 CALL sgebrd( m, n, a, lda, s, work( ie ), work( itauq ),
655 $ work( itaup ), work( iwork ), lwork-iwork+1,
661 CALL sormbr(
'Q',
'L',
'T', m, nrhs, n, a, lda, work( itauq ),
662 $ b, ldb, work( iwork ), lwork-iwork+1, info )
667 CALL sorgbr(
'P', m, n, m, a, lda, work( itaup ),
668 $ work( iwork ), lwork-iwork+1, info )
676 CALL sbdsqr(
'L', m, n, 0, nrhs, s, work( ie ), a, lda, dum,
677 $ 1, b, ldb, work( iwork ), info )
683 thr = max( rcond*s( 1 ), sfmin )
685 $ thr = max( eps*s( 1 ), sfmin )
688 IF( s( i ).GT.thr )
THEN
689 CALL srscl( nrhs, s( i ), b( i, 1 ), ldb )
692 CALL slaset(
'F', 1, nrhs, zero, zero, b( i, 1 ), ldb )
699 IF( lwork.GE.ldb*nrhs .AND. nrhs.GT.1 )
THEN
700 CALL sgemm(
'T',
'N', n, nrhs, m, one, a, lda, b, ldb, zero,
702 CALL slacpy(
'F', n, nrhs, work, ldb, b, ldb )
703 ELSE IF( nrhs.GT.1 )
THEN
705 DO 60 i = 1, nrhs, chunk
706 bl = min( nrhs-i+1, chunk )
707 CALL sgemm(
'T',
'N', n, bl, m, one, a, lda, b( 1, i ),
708 $ ldb, zero, work, n )
709 CALL slacpy(
'F', n, bl, work, n, b( 1, i ), ldb )
712 CALL sgemv(
'T', m, n, one, a, lda, b, 1, zero, work, 1 )
713 CALL scopy( n, work, 1, b, 1 )
719 IF( iascl.EQ.1 )
THEN
720 CALL slascl(
'G', 0, 0, anrm, smlnum, n, nrhs, b, ldb, info )
721 CALL slascl(
'G', 0, 0, smlnum, anrm, minmn, 1, s, minmn,
723 ELSE IF( iascl.EQ.2 )
THEN
724 CALL slascl(
'G', 0, 0, anrm, bignum, n, nrhs, b, ldb, info )
725 CALL slascl(
'G', 0, 0, bignum, anrm, minmn, 1, s, minmn,
728 IF( ibscl.EQ.1 )
THEN
729 CALL slascl(
'G', 0, 0, smlnum, bnrm, n, nrhs, b, ldb, info )
730 ELSE IF( ibscl.EQ.2 )
THEN
731 CALL slascl(
'G', 0, 0, bignum, bnrm, n, nrhs, b, ldb, info )
subroutine slabad(SMALL, LARGE)
SLABAD
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 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 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)
ILAENV
subroutine xerbla(SRNAME, INFO)
XERBLA
subroutine sbdsqr(UPLO, N, NCVT, NRU, NCC, D, E, VT, LDVT, U, LDU, C, LDC, WORK, INFO)
SBDSQR
subroutine sorgbr(VECT, M, N, K, A, LDA, TAU, WORK, LWORK, INFO)
SORGBR
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 sgeqrf(M, N, A, LDA, TAU, WORK, LWORK, INFO)
SGEQRF
subroutine sgebrd(M, N, A, LDA, D, E, TAUQ, TAUP, WORK, LWORK, INFO)
SGEBRD
subroutine sgelqf(M, N, A, LDA, TAU, WORK, LWORK, INFO)
SGELQF
subroutine srscl(N, SA, SX, INCX)
SRSCL multiplies a vector by the reciprocal of a real scalar.
subroutine sormbr(VECT, SIDE, TRANS, M, N, K, A, LDA, TAU, C, LDC, WORK, LWORK, INFO)
SORMBR
subroutine sormqr(SIDE, TRANS, M, N, K, A, LDA, TAU, C, LDC, WORK, LWORK, INFO)
SORMQR
subroutine sormlq(SIDE, TRANS, M, N, K, A, LDA, TAU, C, LDC, WORK, LWORK, INFO)
SORMLQ
subroutine scopy(N, SX, INCX, SY, INCY)
SCOPY
subroutine sgemv(TRANS, M, N, ALPHA, A, LDA, X, INCX, BETA, Y, INCY)
SGEMV
subroutine sgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
SGEMM
real function slamch(CMACH)
SLAMCH