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
238 DOUBLE PRECISION DLAMCH, DLANGE
239 EXTERNAL ilaenv, dlamch, dlange
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
309 smlnum = dlamch(
'S' ) / dlamch(
'P' )
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 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 dtrsm(SIDE, UPLO, TRANSA, DIAG, M, N, ALPHA, A, LDA, B, LDB)
DTRSM
subroutine dcopy(N, DX, INCX, DY, INCY)
DCOPY
subroutine dtzrzf(M, N, A, LDA, TAU, WORK, LWORK, INFO)
DTZRZF
subroutine dgeqp3(M, N, A, LDA, JPVT, TAU, WORK, LWORK, INFO)
DGEQP3
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 dlaic1(JOB, J, X, SEST, W, GAMMA, SESTPR, S, C)
DLAIC1 applies one step of incremental condition estimation.
subroutine dormqr(SIDE, TRANS, M, N, K, A, LDA, TAU, C, LDC, WORK, LWORK, INFO)
DORMQR
subroutine xerbla(SRNAME, INFO)
XERBLA
subroutine dlabad(SMALL, LARGE)
DLABAD
subroutine dormrz(SIDE, TRANS, M, N, K, L, A, LDA, TAU, C, LDC, WORK, LWORK, INFO)
DORMRZ
subroutine dgelsy(M, N, NRHS, A, LDA, B, LDB, JPVT, RCOND, RANK, WORK, LWORK, INFO)
DGELSY solves overdetermined or underdetermined systems for GE matrices