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
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
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 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 sormqr(SIDE, TRANS, M, N, K, A, LDA, TAU, C, LDC, WORK, LWORK, INFO)
SORMQR
REAL function slamch(CMACH)
SLAMCH
subroutine sgeqp3(M, N, A, LDA, JPVT, TAU, WORK, LWORK, INFO)
SGEQP3
subroutine xerbla(SRNAME, INFO)
XERBLA
set ue cd $ADTTMP cat<< EOF > tmp f Program LinearEquations Implicit none Real b(3) integer i
subroutine scopy(N, SX, INCX, SY, INCY)
SCOPY
subroutine sormrz(SIDE, TRANS, M, N, K, L, A, LDA, TAU, C, LDC, WORK, LWORK, INFO)
SORMRZ
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 sgelsy(M, N, NRHS, A, LDA, B, LDB, JPVT, RCOND, RANK, WORK, LWORK, INFO)
SGELSY solves overdetermined or underdetermined systems for GE matrices
INTEGER function ilaenv(ISPEC, NAME, OPTS, N1, N2, N3, N4)
REAL function slange(NORM, M, N, A, LDA, WORK)
SLANGE returns the value of the 1-norm, Frobenius norm, infinity-norm, or the largest absolute value ...
set ue cd $ADTTMP cat<< EOF > tmp f Program LinearEquations Implicit none Real j
subroutine slabad(SMALL, LARGE)
SLABAD
subroutine slaic1(JOB, J, X, SEST, W, GAMMA, SESTPR, S, C)
SLAIC1 applies one step of incremental condition estimation.
subroutine strsm(SIDE, UPLO, TRANSA, DIAG, M, N, ALPHA, A, LDA, B, LDB)
STRSM
subroutine stzrzf(M, N, A, LDA, TAU, WORK, LWORK, INFO)
STZRZF