172 SUBROUTINE dgelss( M, N, NRHS, A, LDA, B, LDB, S, RCOND, RANK,
173 $ work, lwork, info )
181 INTEGER INFO, LDA, LDB, LWORK, M, N, NRHS, RANK
182 DOUBLE PRECISION RCOND
185 DOUBLE PRECISION A( lda, * ), B( ldb, * ), S( * ), WORK( * )
191 DOUBLE PRECISION ZERO, ONE
192 parameter( zero = 0.0d+0, one = 1.0d+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_DGEQRF, LWORK_DORMQR, LWORK_DGEBRD,
200 $ lwork_dormbr, lwork_dorgbr, lwork_dormlq,
202 DOUBLE PRECISION ANRM, BIGNUM, BNRM, EPS, SFMIN, SMLNUM, THR
205 DOUBLE PRECISION DUM( 1 )
214 DOUBLE PRECISION DLAMCH, DLANGE
215 EXTERNAL ilaenv, dlamch, dlange
227 lquery = ( lwork.EQ.-1 )
230 ELSE IF( n.LT.0 )
THEN
232 ELSE IF( nrhs.LT.0 )
THEN
234 ELSE IF( lda.LT.max( 1, m ) )
THEN
236 ELSE IF( ldb.LT.max( 1, maxmn ) )
THEN
250 IF( minmn.GT.0 )
THEN
252 mnthr = ilaenv( 6,
'DGELSS',
' ', m, n, nrhs, -1 )
253 IF( m.GE.n .AND. m.GE.mnthr )
THEN
259 CALL
dgeqrf( m, n, a, lda, dum(1), dum(1), -1, info )
262 CALL
dormqr(
'L',
'T', m, nrhs, n, a, lda, dum(1), b,
263 $ ldb, dum(1), -1, info )
266 maxwrk = max( maxwrk, n + lwork_dgeqrf )
267 maxwrk = max( maxwrk, n + lwork_dormqr )
275 bdspac = max( 1, 5*n )
277 CALL
dgebrd( mm, n, a, lda, s, dum(1), dum(1),
278 $ dum(1), dum(1), -1, info )
281 CALL
dormbr(
'Q',
'L',
'T', mm, nrhs, n, a, lda, dum(1),
282 $ b, ldb, dum(1), -1, info )
285 CALL
dorgbr(
'P', n, n, n, a, lda, dum(1),
289 maxwrk = max( maxwrk, 3*n + lwork_dgebrd )
290 maxwrk = max( maxwrk, 3*n + lwork_dormbr )
291 maxwrk = max( maxwrk, 3*n + lwork_dorgbr )
292 maxwrk = max( maxwrk, bdspac )
293 maxwrk = max( maxwrk, n*nrhs )
294 minwrk = max( 3*n + mm, 3*n + nrhs, bdspac )
295 maxwrk = max( minwrk, maxwrk )
301 bdspac = max( 1, 5*m )
302 minwrk = max( 3*m+nrhs, 3*m+n, bdspac )
303 IF( n.GE.mnthr )
THEN
309 CALL
dgelqf( m, n, a, lda, dum(1), dum(1),
313 CALL
dgebrd( m, m, a, lda, s, dum(1), dum(1),
314 $ dum(1), dum(1), -1, info )
317 CALL
dormbr(
'Q',
'L',
'T', m, nrhs, n, a, lda,
318 $ dum(1), b, ldb, dum(1), -1, info )
321 CALL
dorgbr(
'P', m, m, m, a, lda, dum(1),
325 CALL
dormlq(
'L',
'T', n, nrhs, m, a, lda, dum(1),
326 $ b, ldb, dum(1), -1, info )
329 maxwrk = m + lwork_dgelqf
330 maxwrk = max( maxwrk, m*m + 4*m + lwork_dgebrd )
331 maxwrk = max( maxwrk, m*m + 4*m + lwork_dormbr )
332 maxwrk = max( maxwrk, m*m + 4*m + lwork_dorgbr )
333 maxwrk = max( maxwrk, m*m + m + bdspac )
335 maxwrk = max( maxwrk, m*m + m + m*nrhs )
337 maxwrk = max( maxwrk, m*m + 2*m )
339 maxwrk = max( maxwrk, m + lwork_dormlq )
345 CALL
dgebrd( m, n, a, lda, s, dum(1), dum(1),
346 $ dum(1), dum(1), -1, info )
349 CALL
dormbr(
'Q',
'L',
'T', m, nrhs, m, a, lda,
350 $ dum(1), b, ldb, dum(1), -1, info )
353 CALL
dorgbr(
'P', m, n, m, a, lda, dum(1),
356 maxwrk = 3*m + lwork_dgebrd
357 maxwrk = max( maxwrk, 3*m + lwork_dormbr )
358 maxwrk = max( maxwrk, 3*m + lwork_dorgbr )
359 maxwrk = max( maxwrk, bdspac )
360 maxwrk = max( maxwrk, n*nrhs )
363 maxwrk = max( minwrk, maxwrk )
367 IF( lwork.LT.minwrk .AND. .NOT.lquery )
372 CALL
xerbla(
'DGELSS', -info )
374 ELSE IF( lquery )
THEN
380 IF( m.EQ.0 .OR. n.EQ.0 )
THEN
388 sfmin = dlamch(
'S' )
390 bignum = one / smlnum
391 CALL
dlabad( smlnum, bignum )
395 anrm = dlange(
'M', m, n, a, lda, work )
397 IF( anrm.GT.zero .AND. anrm.LT.smlnum )
THEN
401 CALL
dlascl(
'G', 0, 0, anrm, smlnum, m, n, a, lda, info )
403 ELSE IF( anrm.GT.bignum )
THEN
407 CALL
dlascl(
'G', 0, 0, anrm, bignum, m, n, a, lda, info )
409 ELSE IF( anrm.EQ.zero )
THEN
413 CALL
dlaset(
'F', max( m, n ), nrhs, zero, zero, b, ldb )
414 CALL
dlaset(
'F', minmn, 1, zero, zero, s, minmn )
421 bnrm = dlange(
'M', m, nrhs, b, ldb, work )
423 IF( bnrm.GT.zero .AND. bnrm.LT.smlnum )
THEN
427 CALL
dlascl(
'G', 0, 0, bnrm, smlnum, m, nrhs, b, ldb, info )
429 ELSE IF( bnrm.GT.bignum )
THEN
433 CALL
dlascl(
'G', 0, 0, bnrm, bignum, m, nrhs, b, ldb, info )
444 IF( m.GE.mnthr )
THEN
455 CALL
dgeqrf( m, n, a, lda, work( itau ), work( iwork ),
456 $ lwork-iwork+1, info )
461 CALL
dormqr(
'L',
'T', m, nrhs, n, a, lda, work( itau ), b,
462 $ ldb, work( iwork ), lwork-iwork+1, info )
467 $ CALL
dlaset(
'L', n-1, n-1, zero, zero, a( 2, 1 ), lda )
478 CALL
dgebrd( mm, n, a, lda, s, work( ie ), work( itauq ),
479 $ work( itaup ), work( iwork ), lwork-iwork+1,
485 CALL
dormbr(
'Q',
'L',
'T', mm, nrhs, n, a, lda, work( itauq ),
486 $ b, ldb, work( iwork ), lwork-iwork+1, info )
491 CALL
dorgbr(
'P', n, n, n, a, lda, work( itaup ),
492 $ work( iwork ), lwork-iwork+1, info )
500 CALL
dbdsqr(
'U', n, n, 0, nrhs, s, work( ie ), a, lda, dum,
501 $ 1, b, ldb, work( iwork ), info )
507 thr = max( rcond*s( 1 ), sfmin )
509 $ thr = max( eps*s( 1 ), sfmin )
512 IF( s( i ).GT.thr )
THEN
513 CALL
drscl( nrhs, s( i ), b( i, 1 ), ldb )
516 CALL
dlaset(
'F', 1, nrhs, zero, zero, b( i, 1 ), ldb )
523 IF( lwork.GE.ldb*nrhs .AND. nrhs.GT.1 )
THEN
524 CALL
dgemm(
'T',
'N', n, nrhs, n, one, a, lda, b, ldb, zero,
526 CALL
dlacpy(
'G', n, nrhs, work, ldb, b, ldb )
527 ELSE IF( nrhs.GT.1 )
THEN
529 DO 20 i = 1, nrhs, chunk
530 bl = min( nrhs-i+1, chunk )
531 CALL
dgemm(
'T',
'N', n, bl, n, one, a, lda, b( 1, i ),
532 $ ldb, zero, work, n )
533 CALL
dlacpy(
'G', n, bl, work, n, b( 1, i ), ldb )
536 CALL
dgemv(
'T', n, n, one, a, lda, b, 1, zero, work, 1 )
537 CALL
dcopy( n, work, 1, b, 1 )
540 ELSE IF( n.GE.mnthr .AND. lwork.GE.4*m+m*m+
541 $ max( m, 2*m-4, nrhs, n-3*m ) )
THEN
547 IF( lwork.GE.max( 4*m+m*lda+max( m, 2*m-4, nrhs, n-3*m ),
548 $ m*lda+m+m*nrhs ) )ldwork = lda
555 CALL
dgelqf( m, n, a, lda, work( itau ), work( iwork ),
556 $ lwork-iwork+1, info )
561 CALL
dlacpy(
'L', m, m, a, lda, work( il ), ldwork )
562 CALL
dlaset(
'U', m-1, m-1, zero, zero, work( il+ldwork ),
572 CALL
dgebrd( m, m, work( il ), ldwork, s, work( ie ),
573 $ work( itauq ), work( itaup ), work( iwork ),
574 $ lwork-iwork+1, info )
579 CALL
dormbr(
'Q',
'L',
'T', m, nrhs, m, work( il ), ldwork,
580 $ work( itauq ), b, ldb, work( iwork ),
581 $ lwork-iwork+1, info )
586 CALL
dorgbr(
'P', m, m, m, work( il ), ldwork, work( itaup ),
587 $ work( iwork ), lwork-iwork+1, info )
595 CALL
dbdsqr(
'U', m, m, 0, nrhs, s, work( ie ), work( il ),
596 $ ldwork, a, lda, b, ldb, work( iwork ), info )
602 thr = max( rcond*s( 1 ), sfmin )
604 $ thr = max( eps*s( 1 ), sfmin )
607 IF( s( i ).GT.thr )
THEN
608 CALL
drscl( nrhs, s( i ), b( i, 1 ), ldb )
611 CALL
dlaset(
'F', 1, nrhs, zero, zero, b( i, 1 ), ldb )
619 IF( lwork.GE.ldb*nrhs+iwork-1 .AND. nrhs.GT.1 )
THEN
620 CALL
dgemm(
'T',
'N', m, nrhs, m, one, work( il ), ldwork,
621 $ b, ldb, zero, work( iwork ), ldb )
622 CALL
dlacpy(
'G', m, nrhs, work( iwork ), ldb, b, ldb )
623 ELSE IF( nrhs.GT.1 )
THEN
624 chunk = ( lwork-iwork+1 ) / m
625 DO 40 i = 1, nrhs, chunk
626 bl = min( nrhs-i+1, chunk )
627 CALL
dgemm(
'T',
'N', m, bl, m, one, work( il ), ldwork,
628 $ b( 1, i ), ldb, zero, work( iwork ), m )
629 CALL
dlacpy(
'G', m, bl, work( iwork ), m, b( 1, i ),
633 CALL
dgemv(
'T', m, m, one, work( il ), ldwork, b( 1, 1 ),
634 $ 1, zero, work( iwork ), 1 )
635 CALL
dcopy( m, work( iwork ), 1, b( 1, 1 ), 1 )
640 CALL
dlaset(
'F', n-m, nrhs, zero, zero, b( m+1, 1 ), ldb )
646 CALL
dormlq(
'L',
'T', n, nrhs, m, a, lda, work( itau ), b,
647 $ ldb, work( iwork ), lwork-iwork+1, info )
661 CALL
dgebrd( m, n, a, lda, s, work( ie ), work( itauq ),
662 $ work( itaup ), work( iwork ), lwork-iwork+1,
668 CALL
dormbr(
'Q',
'L',
'T', m, nrhs, n, a, lda, work( itauq ),
669 $ b, ldb, work( iwork ), lwork-iwork+1, info )
674 CALL
dorgbr(
'P', m, n, m, a, lda, work( itaup ),
675 $ work( iwork ), lwork-iwork+1, info )
683 CALL
dbdsqr(
'L', m, n, 0, nrhs, s, work( ie ), a, lda, dum,
684 $ 1, b, ldb, work( iwork ), info )
690 thr = max( rcond*s( 1 ), sfmin )
692 $ thr = max( eps*s( 1 ), sfmin )
695 IF( s( i ).GT.thr )
THEN
696 CALL
drscl( nrhs, s( i ), b( i, 1 ), ldb )
699 CALL
dlaset(
'F', 1, nrhs, zero, zero, b( i, 1 ), ldb )
706 IF( lwork.GE.ldb*nrhs .AND. nrhs.GT.1 )
THEN
707 CALL
dgemm(
'T',
'N', n, nrhs, m, one, a, lda, b, ldb, zero,
709 CALL
dlacpy(
'F', n, nrhs, work, ldb, b, ldb )
710 ELSE IF( nrhs.GT.1 )
THEN
712 DO 60 i = 1, nrhs, chunk
713 bl = min( nrhs-i+1, chunk )
714 CALL
dgemm(
'T',
'N', n, bl, m, one, a, lda, b( 1, i ),
715 $ ldb, zero, work, n )
716 CALL
dlacpy(
'F', n, bl, work, n, b( 1, i ), ldb )
719 CALL
dgemv(
'T', m, n, one, a, lda, b, 1, zero, work, 1 )
720 CALL
dcopy( n, work, 1, b, 1 )
726 IF( iascl.EQ.1 )
THEN
727 CALL
dlascl(
'G', 0, 0, anrm, smlnum, n, nrhs, b, ldb, info )
728 CALL
dlascl(
'G', 0, 0, smlnum, anrm, minmn, 1, s, minmn,
730 ELSE IF( iascl.EQ.2 )
THEN
731 CALL
dlascl(
'G', 0, 0, anrm, bignum, n, nrhs, b, ldb, info )
732 CALL
dlascl(
'G', 0, 0, bignum, anrm, minmn, 1, s, minmn,
735 IF( ibscl.EQ.1 )
THEN
736 CALL
dlascl(
'G', 0, 0, smlnum, bnrm, n, nrhs, b, ldb, info )
737 ELSE IF( ibscl.EQ.2 )
THEN
738 CALL
dlascl(
'G', 0, 0, bignum, bnrm, n, nrhs, b, ldb, info )
subroutine dgelqf(M, N, A, LDA, TAU, WORK, LWORK, INFO)
DGELQF
subroutine dgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
DGEMM
subroutine dormbr(VECT, SIDE, TRANS, M, N, K, A, LDA, TAU, C, LDC, WORK, LWORK, INFO)
DORMBR
subroutine dcopy(N, DX, INCX, DY, INCY)
DCOPY
subroutine dgeqrf(M, N, A, LDA, TAU, WORK, LWORK, INFO)
DGEQRF
subroutine dormqr(SIDE, TRANS, M, N, K, A, LDA, TAU, C, LDC, WORK, LWORK, INFO)
DORMQR
subroutine xerbla(SRNAME, INFO)
XERBLA
subroutine dormlq(SIDE, TRANS, M, N, K, A, LDA, TAU, C, LDC, WORK, LWORK, INFO)
DORMLQ
subroutine dlascl(TYPE, KL, KU, CFROM, CTO, M, N, A, LDA, INFO)
DLASCL multiplies a general rectangular matrix by a real scalar defined as cto/cfrom.
subroutine dlabad(SMALL, LARGE)
DLABAD
subroutine dgebrd(M, N, A, LDA, D, E, TAUQ, TAUP, WORK, LWORK, INFO)
DGEBRD
subroutine dgelss(M, N, NRHS, A, LDA, B, LDB, S, RCOND, RANK, WORK, LWORK, INFO)
DGELSS solves overdetermined or underdetermined systems for GE matrices
subroutine dorgbr(VECT, M, N, K, A, LDA, TAU, WORK, LWORK, INFO)
DORGBR
subroutine dlacpy(UPLO, M, N, A, LDA, B, LDB)
DLACPY copies all or part of one two-dimensional array to another.
subroutine dlaset(UPLO, M, N, ALPHA, BETA, A, LDA)
DLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values...
subroutine drscl(N, SA, SX, INCX)
DRSCL multiplies a vector by the reciprocal of a real scalar.
subroutine dgemv(TRANS, M, N, ALPHA, A, LDA, X, INCX, BETA, Y, INCY)
DGEMV
subroutine dbdsqr(UPLO, N, NCVT, NRU, NCC, D, E, VT, LDVT, U, LDU, C, LDC, WORK, INFO)
DBDSQR