204 SUBROUTINE sgelsy( M, N, NRHS, A, LDA, B, LDB, JPVT, RCOND, RANK,
205 $ work, lwork, info )
213 INTEGER INFO, LDA, LDB, LWORK, M, N, NRHS, RANK
218 REAL A( lda, * ), B( ldb, * ), WORK( * )
225 parameter ( imax = 1, imin = 2 )
227 parameter ( zero = 0.0e+0, one = 1.0e+0 )
231 INTEGER I, IASCL, IBSCL, ISMAX, ISMIN, J, LWKMIN,
232 $ lwkopt, mn, nb, nb1, nb2, nb3, nb4
233 REAL ANRM, BIGNUM, BNRM, C1, C2, S1, S2, SMAX,
234 $ smaxpr, smin, sminpr, smlnum, wsize
239 EXTERNAL ilaenv, slamch, slange
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,
'SGEQRF',
' ', m, n, -1, -1 )
278 nb2 = ilaenv( 1,
'SGERQF',
' ', m, n, -1, -1 )
279 nb3 = ilaenv( 1,
'SORMQR',
' ', m, n, nrhs, -1 )
280 nb4 = ilaenv( 1,
'SORMRQ',
' ', 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(
'SGELSY', -info )
296 ELSE IF( lquery )
THEN
302 IF( mn.EQ.0 .OR. nrhs.EQ.0 )
THEN
309 smlnum = slamch(
'S' ) / slamch(
'P' )
310 bignum = one / smlnum
311 CALL slabad( smlnum, bignum )
315 anrm = slange(
'M', m, n, a, lda, work )
317 IF( anrm.GT.zero .AND. anrm.LT.smlnum )
THEN
321 CALL slascl(
'G', 0, 0, anrm, smlnum, m, n, a, lda, info )
323 ELSE IF( anrm.GT.bignum )
THEN
327 CALL slascl(
'G', 0, 0, anrm, bignum, m, n, a, lda, info )
329 ELSE IF( anrm.EQ.zero )
THEN
333 CALL slaset(
'F', max( m, n ), nrhs, zero, zero, b, ldb )
338 bnrm = slange(
'M', m, nrhs, b, ldb, work )
340 IF( bnrm.GT.zero .AND. bnrm.LT.smlnum )
THEN
344 CALL slascl(
'G', 0, 0, bnrm, smlnum, m, nrhs, b, ldb, info )
346 ELSE IF( bnrm.GT.bignum )
THEN
350 CALL slascl(
'G', 0, 0, bnrm, bignum, m, nrhs, b, ldb, info )
357 CALL sgeqp3( 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 slaset(
'F', max( m, n ), nrhs, zero, zero, b, ldb )
379 IF( rank.LT.mn )
THEN
381 CALL slaic1( imin, rank, work( ismin ), smin, a( 1, i ),
382 $ a( i, i ), sminpr, s1, c1 )
383 CALL slaic1( 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 stzrzf( rank, n, a, lda, work( mn+1 ), work( 2*mn+1 ),
417 CALL sormqr(
'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 strsm(
'Left',
'Upper',
'No transpose',
'Non-unit', rank,
426 $ nrhs, one, a, lda, b, ldb )
429 DO 30 i = rank + 1, n
437 CALL sormrz(
'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 scopy( n, work( 1 ), 1, b( 1, j ), 1 )
457 IF( iascl.EQ.1 )
THEN
458 CALL slascl(
'G', 0, 0, anrm, smlnum, n, nrhs, b, ldb, info )
459 CALL slascl(
'U', 0, 0, smlnum, anrm, rank, rank, a, lda,
461 ELSE IF( iascl.EQ.2 )
THEN
462 CALL slascl(
'G', 0, 0, anrm, bignum, n, nrhs, b, ldb, info )
463 CALL slascl(
'U', 0, 0, bignum, anrm, rank, rank, a, lda,
466 IF( ibscl.EQ.1 )
THEN
467 CALL slascl(
'G', 0, 0, smlnum, bnrm, n, nrhs, b, ldb, info )
468 ELSE IF( ibscl.EQ.2 )
THEN
469 CALL slascl(
'G', 0, 0, bignum, bnrm, n, nrhs, b, ldb, info )
subroutine strsm(SIDE, UPLO, TRANSA, DIAG, M, N, ALPHA, A, LDA, B, LDB)
STRSM
subroutine slaic1(JOB, J, X, SEST, W, GAMMA, SESTPR, S, C)
SLAIC1 applies one step of incremental condition estimation.
subroutine sormqr(SIDE, TRANS, M, N, K, A, LDA, TAU, C, LDC, WORK, LWORK, INFO)
SORMQR
subroutine xerbla(SRNAME, INFO)
XERBLA
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 sgelsy(M, N, NRHS, A, LDA, B, LDB, JPVT, RCOND, RANK, WORK, LWORK, INFO)
SGELSY solves overdetermined or underdetermined systems for GE matrices
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 sgeqp3(M, N, A, LDA, JPVT, TAU, WORK, LWORK, INFO)
SGEQP3
subroutine stzrzf(M, N, A, LDA, TAU, WORK, LWORK, INFO)
STZRZF
subroutine slabad(SMALL, LARGE)
SLABAD
subroutine sormrz(SIDE, TRANS, M, N, K, L, A, LDA, TAU, C, LDC, WORK, LWORK, INFO)
SORMRZ
subroutine scopy(N, SX, INCX, SY, INCY)
SCOPY