176 SUBROUTINE dgelsx( M, N, NRHS, A, LDA, B, LDB, JPVT, RCOND, RANK,
184 INTEGER INFO, LDA, LDB, M, N, NRHS, RANK
185 DOUBLE PRECISION RCOND
189 DOUBLE PRECISION A( LDA, * ), B( LDB, * ), WORK( * )
196 parameter( imax = 1, imin = 2 )
197 DOUBLE PRECISION ZERO, ONE, DONE, NTDONE
198 parameter( zero = 0.0d0, one = 1.0d0, done = zero,
202 INTEGER I, IASCL, IBSCL, ISMAX, ISMIN, J, K, MN
203 DOUBLE PRECISION ANRM, BIGNUM, BNRM, C1, C2, S1, S2, SMAX,
204 $ smaxpr, smin, sminpr, smlnum, t1, t2
207 DOUBLE PRECISION DLAMCH, DLANGE
208 EXTERNAL dlamch, dlange
215 INTRINSIC abs, max, min
228 ELSE IF( n.LT.0 )
THEN
230 ELSE IF( nrhs.LT.0 )
THEN
232 ELSE IF( lda.LT.max( 1, m ) )
THEN
234 ELSE IF( ldb.LT.max( 1, m, n ) )
THEN
239 CALL xerbla(
'DGELSX', -info )
245 IF( min( m, n, nrhs ).EQ.0 )
THEN
252 smlnum = dlamch(
'S' ) / dlamch(
'P' )
253 bignum = one / smlnum
254 CALL dlabad( smlnum, bignum )
258 anrm = dlange(
'M', m, n, a, lda, work )
260 IF( anrm.GT.zero .AND. anrm.LT.smlnum )
THEN
264 CALL dlascl(
'G', 0, 0, anrm, smlnum, m, n, a, lda, info )
266 ELSE IF( anrm.GT.bignum )
THEN
270 CALL dlascl(
'G', 0, 0, anrm, bignum, m, n, a, lda, info )
272 ELSE IF( anrm.EQ.zero )
THEN
276 CALL dlaset(
'F', max( m, n ), nrhs, zero, zero, b, ldb )
281 bnrm = dlange(
'M', m, nrhs, b, ldb, work )
283 IF( bnrm.GT.zero .AND. bnrm.LT.smlnum )
THEN
287 CALL dlascl(
'G', 0, 0, bnrm, smlnum, m, nrhs, b, ldb, info )
289 ELSE IF( bnrm.GT.bignum )
THEN
293 CALL dlascl(
'G', 0, 0, bnrm, bignum, m, nrhs, b, ldb, info )
300 CALL dgeqpf( m, n, a, lda, jpvt, work( 1 ), work( mn+1 ), info )
309 smax = abs( a( 1, 1 ) )
311 IF( abs( a( 1, 1 ) ).EQ.zero )
THEN
313 CALL dlaset(
'F', max( m, n ), nrhs, zero, zero, b, ldb )
320 IF( rank.LT.mn )
THEN
322 CALL dlaic1( imin, rank, work( ismin ), smin, a( 1, i ),
323 $ a( i, i ), sminpr, s1, c1 )
324 CALL dlaic1( imax, rank, work( ismax ), smax, a( 1, i ),
325 $ a( i, i ), smaxpr, s2, c2 )
327 IF( smaxpr*rcond.LE.sminpr )
THEN
329 work( ismin+i-1 ) = s1*work( ismin+i-1 )
330 work( ismax+i-1 ) = s2*work( ismax+i-1 )
332 work( ismin+rank ) = c1
333 work( ismax+rank ) = c2
348 $
CALL dtzrqf( rank, n, a, lda, work( mn+1 ), info )
354 CALL dorm2r(
'Left',
'Transpose', m, nrhs, mn, a, lda, work( 1 ),
355 $ b, ldb, work( 2*mn+1 ), info )
361 CALL dtrsm(
'Left',
'Upper',
'No transpose',
'Non-unit', rank,
362 $ nrhs, one, a, lda, b, ldb )
364 DO 40 i = rank + 1, n
374 CALL dlatzm(
'Left', n-rank+1, nrhs, a( i, rank+1 ), lda,
375 $ work( mn+i ), b( i, 1 ), b( rank+1, 1 ), ldb,
386 work( 2*mn+i ) = ntdone
389 IF( work( 2*mn+i ).EQ.ntdone )
THEN
390 IF( jpvt( i ).NE.i )
THEN
393 t2 = b( jpvt( k ), j )
395 b( jpvt( k ), j ) = t1
396 work( 2*mn+k ) = done
399 t2 = b( jpvt( k ), j )
403 work( 2*mn+k ) = done
411 IF( iascl.EQ.1 )
THEN
412 CALL dlascl(
'G', 0, 0, anrm, smlnum, n, nrhs, b, ldb, info )
413 CALL dlascl(
'U', 0, 0, smlnum, anrm, rank, rank, a, lda,
415 ELSE IF( iascl.EQ.2 )
THEN
416 CALL dlascl(
'G', 0, 0, anrm, bignum, n, nrhs, b, ldb, info )
417 CALL dlascl(
'U', 0, 0, bignum, anrm, rank, rank, a, lda,
420 IF( ibscl.EQ.1 )
THEN
421 CALL dlascl(
'G', 0, 0, smlnum, bnrm, n, nrhs, b, ldb, info )
422 ELSE IF( ibscl.EQ.2 )
THEN
423 CALL dlascl(
'G', 0, 0, bignum, bnrm, n, nrhs, b, ldb, info )
subroutine dlabad(SMALL, LARGE)
DLABAD
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 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 xerbla(SRNAME, INFO)
XERBLA
subroutine dtrsm(SIDE, UPLO, TRANSA, DIAG, M, N, ALPHA, A, LDA, B, LDB)
DTRSM
subroutine dgeqpf(M, N, A, LDA, JPVT, TAU, WORK, INFO)
DGEQPF
subroutine dgelsx(M, N, NRHS, A, LDA, B, LDB, JPVT, RCOND, RANK, WORK, INFO)
DGELSX solves overdetermined or underdetermined systems for GE matrices
subroutine dlaic1(JOB, J, X, SEST, W, GAMMA, SESTPR, S, C)
DLAIC1 applies one step of incremental condition estimation.
subroutine dorm2r(SIDE, TRANS, M, N, K, A, LDA, TAU, C, LDC, WORK, INFO)
DORM2R multiplies a general matrix by the orthogonal matrix from a QR factorization determined by sge...
subroutine dtzrqf(M, N, A, LDA, TAU, INFO)
DTZRQF
subroutine dlatzm(SIDE, M, N, V, INCV, TAU, C1, C2, LDC, WORK)
DLATZM