209 SUBROUTINE dgelsd( M, N, NRHS, A, LDA, B, LDB, S, RCOND, RANK,
210 $ work, lwork, iwork, info )
218 INTEGER INFO, LDA, LDB, LWORK, M, N, NRHS, RANK
219 DOUBLE PRECISION RCOND
223 DOUBLE PRECISION A( lda, * ), B( ldb, * ), S( * ), WORK( * )
229 DOUBLE PRECISION ZERO, ONE, TWO
230 parameter ( zero = 0.0d0, one = 1.0d0, two = 2.0d0 )
234 INTEGER IASCL, IBSCL, IE, IL, ITAU, ITAUP, ITAUQ,
235 $ ldwork, liwork, maxmn, maxwrk, minmn, minwrk,
236 $ mm, mnthr, nlvl, nwork, smlsiz, wlalsd
237 DOUBLE PRECISION ANRM, BIGNUM, BNRM, EPS, SFMIN, SMLNUM
245 DOUBLE PRECISION DLAMCH, DLANGE
246 EXTERNAL ilaenv, dlamch, dlange
249 INTRINSIC dble, int, log, max, min
258 mnthr = ilaenv( 6,
'DGELSD',
' ', m, n, nrhs, -1 )
259 lquery = ( lwork.EQ.-1 )
262 ELSE IF( n.LT.0 )
THEN
264 ELSE IF( nrhs.LT.0 )
THEN
266 ELSE IF( lda.LT.max( 1, m ) )
THEN
268 ELSE IF( ldb.LT.max( 1, maxmn ) )
THEN
272 smlsiz = ilaenv( 9,
'DGELSD',
' ', 0, 0, 0, 0 )
283 minmn = max( 1, minmn )
284 nlvl = max( int( log( dble( minmn ) / dble( smlsiz+1 ) ) /
285 $ log( two ) ) + 1, 0 )
289 liwork = 3*minmn*nlvl + 11*minmn
291 IF( m.GE.n .AND. m.GE.mnthr )
THEN
296 maxwrk = max( maxwrk, n+n*ilaenv( 1,
'DGEQRF',
' ', m, n,
298 maxwrk = max( maxwrk, n+nrhs*
299 $ ilaenv( 1,
'DORMQR',
'LT', m, nrhs, n, -1 ) )
305 maxwrk = max( maxwrk, 3*n+( mm+n )*
306 $ ilaenv( 1,
'DGEBRD',
' ', mm, n, -1, -1 ) )
307 maxwrk = max( maxwrk, 3*n+nrhs*
308 $ ilaenv( 1,
'DORMBR',
'QLT', mm, nrhs, n, -1 ) )
309 maxwrk = max( maxwrk, 3*n+( n-1 )*
310 $ ilaenv( 1,
'DORMBR',
'PLN', n, nrhs, n, -1 ) )
311 wlalsd = 9*n+2*n*smlsiz+8*n*nlvl+n*nrhs+(smlsiz+1)**2
312 maxwrk = max( maxwrk, 3*n+wlalsd )
313 minwrk = max( 3*n+mm, 3*n+nrhs, 3*n+wlalsd )
316 wlalsd = 9*m+2*m*smlsiz+8*m*nlvl+m*nrhs+(smlsiz+1)**2
317 IF( n.GE.mnthr )
THEN
322 maxwrk = m + m*ilaenv( 1,
'DGELQF',
' ', m, n, -1, -1 )
323 maxwrk = max( maxwrk, m*m+4*m+2*m*
324 $ ilaenv( 1,
'DGEBRD',
' ', m, m, -1, -1 ) )
325 maxwrk = max( maxwrk, m*m+4*m+nrhs*
326 $ ilaenv( 1,
'DORMBR',
'QLT', m, nrhs, m, -1 ) )
327 maxwrk = max( maxwrk, m*m+4*m+( m-1 )*
328 $ ilaenv( 1,
'DORMBR',
'PLN', m, nrhs, m, -1 ) )
330 maxwrk = max( maxwrk, m*m+m+m*nrhs )
332 maxwrk = max( maxwrk, m*m+2*m )
334 maxwrk = max( maxwrk, m+nrhs*
335 $ ilaenv( 1,
'DORMLQ',
'LT', n, nrhs, m, -1 ) )
336 maxwrk = max( maxwrk, m*m+4*m+wlalsd )
339 maxwrk = max( maxwrk,
340 $ 4*m+m*m+max( m, 2*m-4, nrhs, n-3*m ) )
345 maxwrk = 3*m + ( n+m )*ilaenv( 1,
'DGEBRD',
' ', m, n,
347 maxwrk = max( maxwrk, 3*m+nrhs*
348 $ ilaenv( 1,
'DORMBR',
'QLT', m, nrhs, n, -1 ) )
349 maxwrk = max( maxwrk, 3*m+m*
350 $ ilaenv( 1,
'DORMBR',
'PLN', n, nrhs, m, -1 ) )
351 maxwrk = max( maxwrk, 3*m+wlalsd )
353 minwrk = max( 3*m+nrhs, 3*m+m, 3*m+wlalsd )
355 minwrk = min( minwrk, maxwrk )
359 IF( lwork.LT.minwrk .AND. .NOT.lquery )
THEN
365 CALL xerbla(
'DGELSD', -info )
367 ELSE IF( lquery )
THEN
373 IF( m.EQ.0 .OR. n.EQ.0 )
THEN
381 sfmin = dlamch(
'S' )
383 bignum = one / smlnum
384 CALL dlabad( smlnum, bignum )
388 anrm = dlange(
'M', m, n, a, lda, work )
390 IF( anrm.GT.zero .AND. anrm.LT.smlnum )
THEN
394 CALL dlascl(
'G', 0, 0, anrm, smlnum, m, n, a, lda, info )
396 ELSE IF( anrm.GT.bignum )
THEN
400 CALL dlascl(
'G', 0, 0, anrm, bignum, m, n, a, lda, info )
402 ELSE IF( anrm.EQ.zero )
THEN
406 CALL dlaset(
'F', max( m, n ), nrhs, zero, zero, b, ldb )
407 CALL dlaset(
'F', minmn, 1, zero, zero, s, 1 )
414 bnrm = dlange(
'M', m, nrhs, b, ldb, work )
416 IF( bnrm.GT.zero .AND. bnrm.LT.smlnum )
THEN
420 CALL dlascl(
'G', 0, 0, bnrm, smlnum, m, nrhs, b, ldb, info )
422 ELSE IF( bnrm.GT.bignum )
THEN
426 CALL dlascl(
'G', 0, 0, bnrm, bignum, m, nrhs, b, ldb, info )
433 $
CALL dlaset(
'F', n-m, nrhs, zero, zero, b( m+1, 1 ), ldb )
442 IF( m.GE.mnthr )
THEN
453 CALL dgeqrf( m, n, a, lda, work( itau ), work( nwork ),
454 $ lwork-nwork+1, info )
459 CALL dormqr(
'L',
'T', m, nrhs, n, a, lda, work( itau ), b,
460 $ ldb, work( nwork ), lwork-nwork+1, info )
465 CALL dlaset(
'L', n-1, n-1, zero, zero, a( 2, 1 ), lda )
477 CALL dgebrd( mm, n, a, lda, s, work( ie ), work( itauq ),
478 $ work( itaup ), work( nwork ), lwork-nwork+1,
484 CALL dormbr(
'Q',
'L',
'T', mm, nrhs, n, a, lda, work( itauq ),
485 $ b, ldb, work( nwork ), lwork-nwork+1, info )
489 CALL dlalsd(
'U', smlsiz, n, nrhs, s, work( ie ), b, ldb,
490 $ rcond, rank, work( nwork ), iwork, info )
497 CALL dormbr(
'P',
'L',
'N', n, nrhs, n, a, lda, work( itaup ),
498 $ b, ldb, work( nwork ), lwork-nwork+1, info )
500 ELSE IF( n.GE.mnthr .AND. lwork.GE.4*m+m*m+
501 $ max( m, 2*m-4, nrhs, n-3*m, wlalsd ) )
THEN
507 IF( lwork.GE.max( 4*m+m*lda+max( m, 2*m-4, nrhs, n-3*m ),
508 $ m*lda+m+m*nrhs, 4*m+m*lda+wlalsd ) )ldwork = lda
515 CALL dgelqf( m, n, a, lda, work( itau ), work( nwork ),
516 $ lwork-nwork+1, info )
521 CALL dlacpy(
'L', m, m, a, lda, work( il ), ldwork )
522 CALL dlaset(
'U', m-1, m-1, zero, zero, work( il+ldwork ),
532 CALL dgebrd( m, m, work( il ), ldwork, s, work( ie ),
533 $ work( itauq ), work( itaup ), work( nwork ),
534 $ lwork-nwork+1, info )
539 CALL dormbr(
'Q',
'L',
'T', m, nrhs, m, work( il ), ldwork,
540 $ work( itauq ), b, ldb, work( nwork ),
541 $ lwork-nwork+1, info )
545 CALL dlalsd(
'U', smlsiz, m, nrhs, s, work( ie ), b, ldb,
546 $ rcond, rank, work( nwork ), iwork, info )
553 CALL dormbr(
'P',
'L',
'N', m, nrhs, m, work( il ), ldwork,
554 $ work( itaup ), b, ldb, work( nwork ),
555 $ lwork-nwork+1, info )
559 CALL dlaset(
'F', n-m, nrhs, zero, zero, b( m+1, 1 ), ldb )
565 CALL dormlq(
'L',
'T', n, nrhs, m, a, lda, work( itau ), b,
566 $ ldb, work( nwork ), lwork-nwork+1, info )
580 CALL dgebrd( m, n, a, lda, s, work( ie ), work( itauq ),
581 $ work( itaup ), work( nwork ), lwork-nwork+1,
587 CALL dormbr(
'Q',
'L',
'T', m, nrhs, n, a, lda, work( itauq ),
588 $ b, ldb, work( nwork ), lwork-nwork+1, info )
592 CALL dlalsd(
'L', smlsiz, m, nrhs, s, work( ie ), b, ldb,
593 $ rcond, rank, work( nwork ), iwork, info )
600 CALL dormbr(
'P',
'L',
'N', n, nrhs, m, a, lda, work( itaup ),
601 $ b, ldb, work( nwork ), lwork-nwork+1, info )
607 IF( iascl.EQ.1 )
THEN
608 CALL dlascl(
'G', 0, 0, anrm, smlnum, n, nrhs, b, ldb, info )
609 CALL dlascl(
'G', 0, 0, smlnum, anrm, minmn, 1, s, minmn,
611 ELSE IF( iascl.EQ.2 )
THEN
612 CALL dlascl(
'G', 0, 0, anrm, bignum, n, nrhs, b, ldb, info )
613 CALL dlascl(
'G', 0, 0, bignum, anrm, minmn, 1, s, minmn,
616 IF( ibscl.EQ.1 )
THEN
617 CALL dlascl(
'G', 0, 0, smlnum, bnrm, n, nrhs, b, ldb, info )
618 ELSE IF( ibscl.EQ.2 )
THEN
619 CALL dlascl(
'G', 0, 0, bignum, bnrm, n, nrhs, b, ldb, info )
subroutine dgebrd(M, N, A, LDA, D, E, TAUQ, TAUP, WORK, LWORK, INFO)
DGEBRD
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 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 dlacpy(UPLO, M, N, A, LDA, B, LDB)
DLACPY copies all or part of one two-dimensional array to another.
subroutine dgelsd(M, N, NRHS, A, LDA, B, LDB, S, RCOND, RANK, WORK, LWORK, IWORK, INFO)
DGELSD computes the minimum-norm solution to a linear least squares problem for GE matrices ...
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 dlalsd(UPLO, SMLSIZ, N, NRHS, D, E, B, LDB, RCOND, RANK, WORK, IWORK, INFO)
DLALSD uses the singular value decomposition of A to solve the least squares problem.
subroutine dormqr(SIDE, TRANS, M, N, K, A, LDA, TAU, C, LDC, WORK, LWORK, INFO)
DORMQR
subroutine dgelqf(M, N, A, LDA, TAU, WORK, LWORK, INFO)
DGELQF
subroutine xerbla(SRNAME, INFO)
XERBLA
subroutine dlabad(SMALL, LARGE)
DLABAD
subroutine dgeqrf(M, N, A, LDA, TAU, WORK, LWORK, INFO)
DGEQRF