210 SUBROUTINE sgelsd( M, N, NRHS, A, LDA, B, LDB, S, RCOND,
211 $ rank, work, lwork, iwork, info )
219 INTEGER info, lda, ldb, lwork, m, n, nrhs, rank
224 REAL a( lda, * ),
b( ldb, * ), s( * ), work( * )
231 parameter( zero = 0.0e0, one = 1.0e0, two = 2.0e0 )
235 INTEGER iascl, ibscl, ie, il, itau, itaup, itauq,
236 $ ldwork, liwork, maxmn, maxwrk, minmn, minwrk,
237 $ mm, mnthr, nlvl, nwork, smlsiz, wlalsd
238 REAL anrm, bignum, bnrm, eps, sfmin, smlnum
250 INTRINSIC int, log, max, min, real
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
283 IF( minmn.GT.0 )
THEN
284 smlsiz =
ilaenv( 9,
'SGELSD',
' ', 0, 0, 0, 0 )
285 mnthr =
ilaenv( 6,
'SGELSD',
' ', m, n, nrhs, -1 )
286 nlvl = max( int( log(
REAL( MINMN ) /
REAL( SMLSIZ + 1 ) ) /
287 $ log( two ) ) + 1, 0 )
288 liwork = 3*minmn*nlvl + 11*minmn
290 IF( m.GE.n .AND. m.GE.mnthr )
THEN
296 maxwrk = max( maxwrk, n + n*
ilaenv( 1,
'SGEQRF',
' ', m,
298 maxwrk = max( maxwrk, n + nrhs*
ilaenv( 1,
'SORMQR',
'LT',
305 maxwrk = max( maxwrk, 3*n + ( mm + n )*
ilaenv( 1,
306 $
'SGEBRD',
' ', mm, n, -1, -1 ) )
307 maxwrk = max( maxwrk, 3*n + nrhs*
ilaenv( 1,
'SORMBR',
308 $
'QLT', mm, nrhs, n, -1 ) )
309 maxwrk = max( maxwrk, 3*n + ( n - 1 )*
ilaenv( 1,
310 $
'SORMBR',
'PLN', n, nrhs, n, -1 ) )
311 wlalsd = 9*n + 2*n*smlsiz + 8*n*nlvl + n*nrhs +
313 maxwrk = max( maxwrk, 3*n + wlalsd )
314 minwrk = max( 3*n + mm, 3*n + nrhs, 3*n + wlalsd )
317 wlalsd = 9*m + 2*m*smlsiz + 8*m*nlvl + m*nrhs +
319 IF( n.GE.mnthr )
THEN
324 maxwrk = m + m*
ilaenv( 1,
'SGELQF',
' ', m, n, -1,
326 maxwrk = max( maxwrk, m*m + 4*m + 2*m*
ilaenv( 1,
327 $
'SGEBRD',
' ', m, m, -1, -1 ) )
328 maxwrk = max( maxwrk, m*m + 4*m + nrhs*
ilaenv( 1,
329 $
'SORMBR',
'QLT', m, nrhs, m, -1 ) )
330 maxwrk = max( maxwrk, m*m + 4*m + ( m - 1 )*
ilaenv( 1,
331 $
'SORMBR',
'PLN', m, nrhs, m, -1 ) )
333 maxwrk = max( maxwrk, m*m + m + m*nrhs )
335 maxwrk = max( maxwrk, m*m + 2*m )
337 maxwrk = max( maxwrk, m + nrhs*
ilaenv( 1,
'SORMLQ',
338 $
'LT', n, nrhs, m, -1 ) )
339 maxwrk = max( maxwrk, m*m + 4*m + wlalsd )
342 maxwrk = max( maxwrk,
343 $ 4*m+m*m+max( m, 2*m-4, nrhs, n-3*m ) )
348 maxwrk = 3*m + ( n + m )*
ilaenv( 1,
'SGEBRD',
' ', m,
350 maxwrk = max( maxwrk, 3*m + nrhs*
ilaenv( 1,
'SORMBR',
351 $
'QLT', m, nrhs, n, -1 ) )
352 maxwrk = max( maxwrk, 3*m + m*
ilaenv( 1,
'SORMBR',
353 $
'PLN', n, nrhs, m, -1 ) )
354 maxwrk = max( maxwrk, 3*m + wlalsd )
356 minwrk = max( 3*m + nrhs, 3*m + m, 3*m + wlalsd )
359 minwrk = min( minwrk, maxwrk )
363 IF( lwork.LT.minwrk .AND. .NOT.lquery )
THEN
369 CALL
xerbla(
'SGELSD', -info )
371 ELSE IF( lquery )
THEN
377 IF( m.EQ.0 .OR. n.EQ.0 )
THEN
387 bignum = one / smlnum
388 CALL
slabad( smlnum, bignum )
392 anrm =
slange(
'M', m, n, a, lda, work )
394 IF( anrm.GT.zero .AND. anrm.LT.smlnum )
THEN
398 CALL
slascl(
'G', 0, 0, anrm, smlnum, m, n, a, lda, info )
400 ELSE IF( anrm.GT.bignum )
THEN
404 CALL
slascl(
'G', 0, 0, anrm, bignum, m, n, a, lda, info )
406 ELSE IF( anrm.EQ.zero )
THEN
410 CALL
slaset(
'F', max( m, n ), nrhs, zero, zero,
b, ldb )
411 CALL
slaset(
'F', minmn, 1, zero, zero, s, 1 )
418 bnrm =
slange(
'M', m, nrhs,
b, ldb, work )
420 IF( bnrm.GT.zero .AND. bnrm.LT.smlnum )
THEN
424 CALL
slascl(
'G', 0, 0, bnrm, smlnum, m, nrhs,
b, ldb, info )
426 ELSE IF( bnrm.GT.bignum )
THEN
430 CALL
slascl(
'G', 0, 0, bnrm, bignum, m, nrhs,
b, ldb, info )
437 $ CALL
slaset(
'F', n-m, nrhs, zero, zero,
b( m+1, 1 ), ldb )
446 IF( m.GE.mnthr )
THEN
457 CALL
sgeqrf( m, n, a, lda, work( itau ), work( nwork ),
458 $ lwork-nwork+1, info )
463 CALL
sormqr(
'L',
'T', m, nrhs, n, a, lda, work( itau ),
b,
464 $ ldb, work( nwork ), lwork-nwork+1, info )
469 CALL
slaset(
'L', n-1, n-1, zero, zero, a( 2, 1 ), lda )
481 CALL
sgebrd( mm, n, a, lda, s, work( ie ), work( itauq ),
482 $ work( itaup ), work( nwork ), lwork-nwork+1,
488 CALL
sormbr(
'Q',
'L',
'T', mm, nrhs, n, a, lda, work( itauq ),
489 $
b, ldb, work( nwork ), lwork-nwork+1, info )
493 CALL
slalsd(
'U', smlsiz, n, nrhs, s, work( ie ),
b, ldb,
494 $ rcond, rank, work( nwork ), iwork, info )
501 CALL
sormbr(
'P',
'L',
'N', n, nrhs, n, a, lda, work( itaup ),
502 $
b, ldb, work( nwork ), lwork-nwork+1, info )
504 ELSE IF( n.GE.mnthr .AND. lwork.GE.4*m+m*m+
505 $ max( m, 2*m-4, nrhs, n-3*m, wlalsd ) )
THEN
511 IF( lwork.GE.max( 4*m+m*lda+max( m, 2*m-4, nrhs, n-3*m ),
512 $ m*lda+m+m*nrhs, 4*m+m*lda+wlalsd ) )ldwork = lda
519 CALL
sgelqf( m, n, a, lda, work( itau ), work( nwork ),
520 $ lwork-nwork+1, info )
525 CALL
slacpy(
'L', m, m, a, lda, work( il ), ldwork )
526 CALL
slaset(
'U', m-1, m-1, zero, zero, work( il+ldwork ),
536 CALL
sgebrd( m, m, work( il ), ldwork, s, work( ie ),
537 $ work( itauq ), work( itaup ), work( nwork ),
538 $ lwork-nwork+1, info )
543 CALL
sormbr(
'Q',
'L',
'T', m, nrhs, m, work( il ), ldwork,
544 $ work( itauq ),
b, ldb, work( nwork ),
545 $ lwork-nwork+1, info )
549 CALL
slalsd(
'U', smlsiz, m, nrhs, s, work( ie ),
b, ldb,
550 $ rcond, rank, work( nwork ), iwork, info )
557 CALL
sormbr(
'P',
'L',
'N', m, nrhs, m, work( il ), ldwork,
558 $ work( itaup ),
b, ldb, work( nwork ),
559 $ lwork-nwork+1, info )
563 CALL
slaset(
'F', n-m, nrhs, zero, zero,
b( m+1, 1 ), ldb )
569 CALL
sormlq(
'L',
'T', n, nrhs, m, a, lda, work( itau ),
b,
570 $ ldb, work( nwork ), lwork-nwork+1, info )
584 CALL
sgebrd( m, n, a, lda, s, work( ie ), work( itauq ),
585 $ work( itaup ), work( nwork ), lwork-nwork+1,
591 CALL
sormbr(
'Q',
'L',
'T', m, nrhs, n, a, lda, work( itauq ),
592 $
b, ldb, work( nwork ), lwork-nwork+1, info )
596 CALL
slalsd(
'L', smlsiz, m, nrhs, s, work( ie ),
b, ldb,
597 $ rcond, rank, work( nwork ), iwork, info )
604 CALL
sormbr(
'P',
'L',
'N', n, nrhs, m, a, lda, work( itaup ),
605 $
b, ldb, work( nwork ), lwork-nwork+1, info )
611 IF( iascl.EQ.1 )
THEN
612 CALL
slascl(
'G', 0, 0, anrm, smlnum, n, nrhs,
b, ldb, info )
613 CALL
slascl(
'G', 0, 0, smlnum, anrm, minmn, 1, s, minmn,
615 ELSE IF( iascl.EQ.2 )
THEN
616 CALL
slascl(
'G', 0, 0, anrm, bignum, n, nrhs,
b, ldb, info )
617 CALL
slascl(
'G', 0, 0, bignum, anrm, minmn, 1, s, minmn,
620 IF( ibscl.EQ.1 )
THEN
621 CALL
slascl(
'G', 0, 0, smlnum, bnrm, n, nrhs,
b, ldb, info )
622 ELSE IF( ibscl.EQ.2 )
THEN
623 CALL
slascl(
'G', 0, 0, bignum, bnrm, n, nrhs,
b, ldb, info )
subroutine sgebrd(M, N, A, LDA, D, E, TAUQ, TAUP, WORK, LWORK, INFO)
SGEBRD
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 sormlq(SIDE, TRANS, M, N, K, A, LDA, TAU, C, LDC, WORK, LWORK, INFO)
SORMLQ
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 slalsd(UPLO, SMLSIZ, N, NRHS, D, E, B, LDB, RCOND, RANK, WORK, IWORK, INFO)
SLALSD uses the singular value decomposition of A to solve the least squares problem.
subroutine sormqr(SIDE, TRANS, M, N, K, A, LDA, TAU, C, LDC, WORK, LWORK, INFO)
SORMQR
subroutine sgelsd(M, N, NRHS, A, LDA, B, LDB, S, RCOND, RANK, WORK, LWORK, IWORK, INFO)
SGELSD computes the minimum-norm solution to a linear least squares problem for GE matrices ...
subroutine xerbla(SRNAME, INFO)
XERBLA
set ue cd $ADTTMP cat<< EOF > tmp f Program LinearEquations Implicit none Real b(3) integer i
subroutine sgelqf(M, N, A, LDA, TAU, WORK, LWORK, INFO)
SGELQF
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 slacpy(UPLO, M, N, A, LDA, B, LDB)
SLACPY copies all or part of one two-dimensional array to another.
real function slamch(CMACH)
SLAMCH
subroutine sgeqrf(M, N, A, LDA, TAU, WORK, LWORK, INFO)
SGEQRF
subroutine slabad(SMALL, LARGE)
SLABAD
integer function ilaenv(ISPEC, NAME, OPTS, N1, N2, N3, N4)
subroutine sormbr(VECT, SIDE, TRANS, M, N, K, A, LDA, TAU, C, LDC, WORK, LWORK, INFO)
SORMBR