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
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
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 dgelqf(M, N, A, LDA, TAU, WORK, LWORK, INFO)
DGELQF
subroutine dormbr(VECT, SIDE, TRANS, M, N, K, A, LDA, TAU, C, LDC, WORK, LWORK, INFO)
DORMBR
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
set ue cd $ADTTMP cat<< EOF > tmp f Program LinearEquations Implicit none Real b(3) integer i
DOUBLE PRECISION function dlange(NORM, M, N, A, LDA, WORK)
DLANGE returns the value of the 1-norm, Frobenius norm, infinity-norm, or the largest absolute value ...
subroutine dgebrd(M, N, A, LDA, D, E, TAUQ, TAUP, WORK, LWORK, INFO)
DGEBRD
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.
INTEGER function ilaenv(ISPEC, NAME, OPTS, N1, N2, N3, N4)
subroutine dlacpy(UPLO, M, N, A, LDA, B, LDB)
DLACPY copies all or part of one two-dimensional array to another.
DOUBLE PRECISION function dlamch(CMACH)
DLAMCH
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 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 ...