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 dlacpy(UPLO, M, N, A, LDA, B, LDB)
DLACPY copies all or part of one two-dimensional array to another.
subroutine dgebrd(M, N, A, LDA, D, E, TAUQ, TAUP, WORK, LWORK, INFO)
DGEBRD
subroutine dcopy(N, DX, INCX, DY, INCY)
DCOPY
subroutine dgemv(TRANS, M, N, ALPHA, A, LDA, X, INCX, BETA, Y, INCY)
DGEMV
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 dormqr(SIDE, TRANS, M, N, K, A, LDA, TAU, C, LDC, WORK, LWORK, INFO)
DORMQR
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 dgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
DGEMM
subroutine dgelqf(M, N, A, LDA, TAU, WORK, LWORK, INFO)
DGELQF
subroutine dormlq(SIDE, TRANS, M, N, K, A, LDA, TAU, C, LDC, WORK, LWORK, INFO)
DORMLQ
subroutine dormbr(VECT, SIDE, TRANS, M, N, K, A, LDA, TAU, C, LDC, WORK, LWORK, INFO)
DORMBR
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 xerbla(SRNAME, INFO)
XERBLA
subroutine dlabad(SMALL, LARGE)
DLABAD
subroutine dgeqrf(M, N, A, LDA, TAU, WORK, LWORK, INFO)
DGEQRF
subroutine dbdsqr(UPLO, N, NCVT, NRU, NCC, D, E, VT, LDVT, U, LDU, C, LDC, WORK, INFO)
DBDSQR
subroutine drscl(N, SA, SX, INCX)
DRSCL multiplies a vector by the reciprocal of a real scalar.