204 SUBROUTINE dgelsy( M, N, NRHS, A, LDA, B, LDB, JPVT, RCOND, RANK,
205 $ work, lwork, info )
213 INTEGER info, lda, ldb, lwork, m, n, nrhs, rank
214 DOUBLE PRECISION rcond
218 DOUBLE PRECISION a( lda, * ),
b( ldb, * ), work( * )
225 parameter( imax = 1, imin = 2 )
226 DOUBLE PRECISION zero, one
227 parameter( zero = 0.0d+0, one = 1.0d+0 )
231 INTEGER i, iascl, ibscl, ismax, ismin,
j, lwkmin,
232 $ lwkopt, mn, nb, nb1, nb2, nb3, nb4
233 DOUBLE PRECISION anrm, bignum, bnrm, c1, c2, s1, s2, smax,
234 $ smaxpr, smin, sminpr, smlnum, wsize
246 INTRINSIC abs, max, min
257 lquery = ( lwork.EQ.-1 )
260 ELSE IF( n.LT.0 )
THEN
262 ELSE IF( nrhs.LT.0 )
THEN
264 ELSE IF( lda.LT.max( 1, m ) )
THEN
266 ELSE IF( ldb.LT.max( 1, m, n ) )
THEN
273 IF( mn.EQ.0 .OR. nrhs.EQ.0 )
THEN
277 nb1 =
ilaenv( 1,
'DGEQRF',
' ', m, n, -1, -1 )
278 nb2 =
ilaenv( 1,
'DGERQF',
' ', m, n, -1, -1 )
279 nb3 =
ilaenv( 1,
'DORMQR',
' ', m, n, nrhs, -1 )
280 nb4 =
ilaenv( 1,
'DORMRQ',
' ', m, n, nrhs, -1 )
281 nb = max( nb1, nb2, nb3, nb4 )
282 lwkmin = mn + max( 2*mn, n + 1, mn + nrhs )
283 lwkopt = max( lwkmin,
284 $ mn + 2*n + nb*( n + 1 ), 2*mn + nb*nrhs )
288 IF( lwork.LT.lwkmin .AND. .NOT.lquery )
THEN
294 CALL
xerbla(
'DGELSY', -info )
296 ELSE IF( lquery )
THEN
302 IF( mn.EQ.0 .OR. nrhs.EQ.0 )
THEN
310 bignum = one / smlnum
311 CALL
dlabad( smlnum, bignum )
315 anrm =
dlange(
'M', m, n, a, lda, work )
317 IF( anrm.GT.zero .AND. anrm.LT.smlnum )
THEN
321 CALL
dlascl(
'G', 0, 0, anrm, smlnum, m, n, a, lda, info )
323 ELSE IF( anrm.GT.bignum )
THEN
327 CALL
dlascl(
'G', 0, 0, anrm, bignum, m, n, a, lda, info )
329 ELSE IF( anrm.EQ.zero )
THEN
333 CALL
dlaset(
'F', max( m, n ), nrhs, zero, zero,
b, ldb )
338 bnrm =
dlange(
'M', m, nrhs,
b, ldb, work )
340 IF( bnrm.GT.zero .AND. bnrm.LT.smlnum )
THEN
344 CALL
dlascl(
'G', 0, 0, bnrm, smlnum, m, nrhs,
b, ldb, info )
346 ELSE IF( bnrm.GT.bignum )
THEN
350 CALL
dlascl(
'G', 0, 0, bnrm, bignum, m, nrhs,
b, ldb, info )
357 CALL
dgeqp3( m, n, a, lda, jpvt, work( 1 ), work( mn+1 ),
359 wsize = mn + work( mn+1 )
368 smax = abs( a( 1, 1 ) )
370 IF( abs( a( 1, 1 ) ).EQ.zero )
THEN
372 CALL
dlaset(
'F', max( m, n ), nrhs, zero, zero,
b, ldb )
379 IF( rank.LT.mn )
THEN
381 CALL
dlaic1( imin, rank, work( ismin ), smin, a( 1, i ),
382 $ a( i, i ), sminpr, s1, c1 )
383 CALL
dlaic1( imax, rank, work( ismax ), smax, a( 1, i ),
384 $ a( i, i ), smaxpr, s2, c2 )
386 IF( smaxpr*rcond.LE.sminpr )
THEN
388 work( ismin+i-1 ) = s1*work( ismin+i-1 )
389 work( ismax+i-1 ) = s2*work( ismax+i-1 )
391 work( ismin+rank ) = c1
392 work( ismax+rank ) = c2
409 $ CALL
dtzrzf( rank, n, a, lda, work( mn+1 ), work( 2*mn+1 ),
417 CALL
dormqr(
'Left',
'Transpose', m, nrhs, mn, a, lda, work( 1 ),
418 $
b, ldb, work( 2*mn+1 ), lwork-2*mn, info )
419 wsize = max( wsize, 2*mn+work( 2*mn+1 ) )
425 CALL
dtrsm(
'Left',
'Upper',
'No transpose',
'Non-unit', rank,
426 $ nrhs, one, a, lda,
b, ldb )
429 DO 30 i = rank + 1, n
437 CALL
dormrz(
'Left',
'Transpose', n, nrhs, rank, n-rank, a,
438 $ lda, work( mn+1 ),
b, ldb, work( 2*mn+1 ),
448 work( jpvt( i ) ) =
b( i,
j )
450 CALL
dcopy( n, work( 1 ), 1,
b( 1,
j ), 1 )
457 IF( iascl.EQ.1 )
THEN
458 CALL
dlascl(
'G', 0, 0, anrm, smlnum, n, nrhs,
b, ldb, info )
459 CALL
dlascl(
'U', 0, 0, smlnum, anrm, rank, rank, a, lda,
461 ELSE IF( iascl.EQ.2 )
THEN
462 CALL
dlascl(
'G', 0, 0, anrm, bignum, n, nrhs,
b, ldb, info )
463 CALL
dlascl(
'U', 0, 0, bignum, anrm, rank, rank, a, lda,
466 IF( ibscl.EQ.1 )
THEN
467 CALL
dlascl(
'G', 0, 0, smlnum, bnrm, n, nrhs,
b, ldb, info )
468 ELSE IF( ibscl.EQ.2 )
THEN
469 CALL
dlascl(
'G', 0, 0, bignum, bnrm, n, nrhs,
b, ldb, info )
subroutine dtzrzf(M, N, A, LDA, TAU, WORK, LWORK, INFO)
DTZRZF
subroutine dlaic1(JOB, J, X, SEST, W, GAMMA, SESTPR, S, C)
DLAIC1 applies one step of incremental condition estimation.
subroutine dcopy(N, DX, INCX, DY, INCY)
DCOPY
subroutine dormqr(SIDE, TRANS, M, N, K, A, LDA, TAU, C, LDC, WORK, LWORK, INFO)
DORMQR
subroutine xerbla(SRNAME, INFO)
XERBLA
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 dgelsy(M, N, NRHS, A, LDA, B, LDB, JPVT, RCOND, RANK, WORK, LWORK, INFO)
DGELSY solves overdetermined or underdetermined systems for GE matrices
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 dormrz(SIDE, TRANS, M, N, K, L, A, LDA, TAU, C, LDC, WORK, LWORK, INFO)
DORMRZ
subroutine dtrsm(SIDE, UPLO, TRANSA, DIAG, M, N, ALPHA, A, LDA, B, LDB)
DTRSM
INTEGER function ilaenv(ISPEC, NAME, OPTS, N1, N2, N3, N4)
set ue cd $ADTTMP cat<< EOF > tmp f Program LinearEquations Implicit none Real j
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 dgeqp3(M, N, A, LDA, JPVT, TAU, WORK, LWORK, INFO)
DGEQP3