204 SUBROUTINE dgelsy( M, N, NRHS, A, LDA, B, LDB, JPVT, RCOND, RANK,
205 $ WORK, LWORK, INFO )
212 INTEGER INFO, LDA, LDB, LWORK, M, N, NRHS, RANK
213 DOUBLE PRECISION RCOND
217 DOUBLE PRECISION A( LDA, * ), B( LDB, * ), WORK( * )
224 parameter( imax = 1, imin = 2 )
225 DOUBLE PRECISION ZERO, ONE
226 parameter( zero = 0.0d+0, one = 1.0d+0 )
230 INTEGER I, IASCL, IBSCL, ISMAX, ISMIN, J, LWKMIN,
231 $ lwkopt, mn, nb, nb1, nb2, nb3, nb4
232 DOUBLE PRECISION ANRM, BIGNUM, BNRM, C1, C2, S1, S2, SMAX,
233 $ smaxpr, smin, sminpr, smlnum, wsize
237 DOUBLE PRECISION DLAMCH, DLANGE
238 EXTERNAL ilaenv, dlamch, dlange
245 INTRINSIC abs, max, min
256 lquery = ( lwork.EQ.-1 )
259 ELSE IF( n.LT.0 )
THEN
261 ELSE IF( nrhs.LT.0 )
THEN
263 ELSE IF( lda.LT.max( 1, m ) )
THEN
265 ELSE IF( ldb.LT.max( 1, m, n ) )
THEN
272 IF( mn.EQ.0 .OR. nrhs.EQ.0 )
THEN
276 nb1 = ilaenv( 1,
'DGEQRF',
' ', m, n, -1, -1 )
277 nb2 = ilaenv( 1,
'DGERQF',
' ', m, n, -1, -1 )
278 nb3 = ilaenv( 1,
'DORMQR',
' ', m, n, nrhs, -1 )
279 nb4 = ilaenv( 1,
'DORMRQ',
' ', m, n, nrhs, -1 )
280 nb = max( nb1, nb2, nb3, nb4 )
281 lwkmin = mn + max( 2*mn, n + 1, mn + nrhs )
282 lwkopt = max( lwkmin,
283 $ mn + 2*n + nb*( n + 1 ), 2*mn + nb*nrhs )
287 IF( lwork.LT.lwkmin .AND. .NOT.lquery )
THEN
293 CALL xerbla(
'DGELSY', -info )
295 ELSE IF( lquery )
THEN
301 IF( mn.EQ.0 .OR. nrhs.EQ.0 )
THEN
308 smlnum = dlamch(
'S' ) / dlamch(
'P' )
309 bignum = one / smlnum
313 anrm = dlange(
'M', m, n, a, lda, work )
315 IF( anrm.GT.zero .AND. anrm.LT.smlnum )
THEN
319 CALL dlascl(
'G', 0, 0, anrm, smlnum, m, n, a, lda, info )
321 ELSE IF( anrm.GT.bignum )
THEN
325 CALL dlascl(
'G', 0, 0, anrm, bignum, m, n, a, lda, info )
327 ELSE IF( anrm.EQ.zero )
THEN
331 CALL dlaset(
'F', max( m, n ), nrhs, zero, zero, b, ldb )
336 bnrm = dlange(
'M', m, nrhs, b, ldb, work )
338 IF( bnrm.GT.zero .AND. bnrm.LT.smlnum )
THEN
342 CALL dlascl(
'G', 0, 0, bnrm, smlnum, m, nrhs, b, ldb, info )
344 ELSE IF( bnrm.GT.bignum )
THEN
348 CALL dlascl(
'G', 0, 0, bnrm, bignum, m, nrhs, b, ldb, info )
355 CALL dgeqp3( m, n, a, lda, jpvt, work( 1 ), work( mn+1 ),
357 wsize = mn + work( mn+1 )
366 smax = abs( a( 1, 1 ) )
368 IF( abs( a( 1, 1 ) ).EQ.zero )
THEN
370 CALL dlaset(
'F', max( m, n ), nrhs, zero, zero, b, ldb )
377 IF( rank.LT.mn )
THEN
379 CALL dlaic1( imin, rank, work( ismin ), smin, a( 1, i ),
380 $ a( i, i ), sminpr, s1, c1 )
381 CALL dlaic1( imax, rank, work( ismax ), smax, a( 1, i ),
382 $ a( i, i ), smaxpr, s2, c2 )
384 IF( smaxpr*rcond.LE.sminpr )
THEN
386 work( ismin+i-1 ) = s1*work( ismin+i-1 )
387 work( ismax+i-1 ) = s2*work( ismax+i-1 )
389 work( ismin+rank ) = c1
390 work( ismax+rank ) = c2
407 $
CALL dtzrzf( rank, n, a, lda, work( mn+1 ), work( 2*mn+1 ),
415 CALL dormqr(
'Left',
'Transpose', m, nrhs, mn, a, lda, work( 1 ),
416 $ b, ldb, work( 2*mn+1 ), lwork-2*mn, info )
417 wsize = max( wsize, 2*mn+work( 2*mn+1 ) )
423 CALL dtrsm(
'Left',
'Upper',
'No transpose',
'Non-unit', rank,
424 $ nrhs, one, a, lda, b, ldb )
427 DO 30 i = rank + 1, n
435 CALL dormrz(
'Left',
'Transpose', n, nrhs, rank, n-rank, a,
436 $ lda, work( mn+1 ), b, ldb, work( 2*mn+1 ),
446 work( jpvt( i ) ) = b( i, j )
448 CALL dcopy( n, work( 1 ), 1, b( 1, j ), 1 )
455 IF( iascl.EQ.1 )
THEN
456 CALL dlascl(
'G', 0, 0, anrm, smlnum, n, nrhs, b, ldb, info )
457 CALL dlascl(
'U', 0, 0, smlnum, anrm, rank, rank, a, lda,
459 ELSE IF( iascl.EQ.2 )
THEN
460 CALL dlascl(
'G', 0, 0, anrm, bignum, n, nrhs, b, ldb, info )
461 CALL dlascl(
'U', 0, 0, bignum, anrm, rank, rank, a, lda,
464 IF( ibscl.EQ.1 )
THEN
465 CALL dlascl(
'G', 0, 0, smlnum, bnrm, n, nrhs, b, ldb, info )
466 ELSE IF( ibscl.EQ.2 )
THEN
467 CALL dlascl(
'G', 0, 0, bignum, bnrm, n, nrhs, b, ldb, info )
subroutine xerbla(srname, info)
subroutine dcopy(n, dx, incx, dy, incy)
DCOPY
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 dgeqp3(m, n, a, lda, jpvt, tau, work, lwork, info)
DGEQP3
subroutine dlaic1(job, j, x, sest, w, gamma, sestpr, s, c)
DLAIC1 applies one step of incremental condition estimation.
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 dtrsm(side, uplo, transa, diag, m, n, alpha, a, lda, b, ldb)
DTRSM
subroutine dtzrzf(m, n, a, lda, tau, work, lwork, info)
DTZRZF
subroutine dormqr(side, trans, m, n, k, a, lda, tau, c, ldc, work, lwork, info)
DORMQR
subroutine dormrz(side, trans, m, n, k, l, a, lda, tau, c, ldc, work, lwork, info)
DORMRZ