210 SUBROUTINE zgelsy( M, N, NRHS, A, LDA, B, LDB, JPVT, RCOND, RANK,
211 $ work, lwork, rwork, info )
219 INTEGER info, lda, ldb, lwork, m, n, nrhs, rank
220 DOUBLE PRECISION rcond
224 DOUBLE PRECISION rwork( * )
225 COMPLEX*16 a( lda, * ),
b( ldb, * ), work( * )
232 parameter( imax = 1, imin = 2 )
233 DOUBLE PRECISION zero, one
234 parameter( zero = 0.0d+0, one = 1.0d+0 )
235 COMPLEX*16 czero, cone
236 parameter( czero = ( 0.0d+0, 0.0d+0 ),
237 $ cone = ( 1.0d+0, 0.0d+0 ) )
241 INTEGER i, iascl, ibscl, ismax, ismin,
j, lwkopt, mn,
242 $ nb, nb1, nb2, nb3, nb4
243 DOUBLE PRECISION anrm, bignum, bnrm, smax, smaxpr, smin, sminpr,
245 COMPLEX*16 c1, c2, s1, s2
257 INTRINSIC abs, dble, dcmplx, max, min
268 nb1 =
ilaenv( 1,
'ZGEQRF',
' ', m, n, -1, -1 )
269 nb2 =
ilaenv( 1,
'ZGERQF',
' ', m, n, -1, -1 )
270 nb3 =
ilaenv( 1,
'ZUNMQR',
' ', m, n, nrhs, -1 )
271 nb4 =
ilaenv( 1,
'ZUNMRQ',
' ', m, n, nrhs, -1 )
272 nb = max( nb1, nb2, nb3, nb4 )
273 lwkopt = max( 1, mn+2*n+nb*( n+1 ), 2*mn+nb*nrhs )
274 work( 1 ) = dcmplx( lwkopt )
275 lquery = ( lwork.EQ.-1 )
278 ELSE IF( n.LT.0 )
THEN
280 ELSE IF( nrhs.LT.0 )
THEN
282 ELSE IF( lda.LT.max( 1, m ) )
THEN
284 ELSE IF( ldb.LT.max( 1, m, n ) )
THEN
286 ELSE IF( lwork.LT.( mn+max( 2*mn, n+1, mn+nrhs ) ) .AND. .NOT.
292 CALL
xerbla(
'ZGELSY', -info )
294 ELSE IF( lquery )
THEN
300 IF( min( m, n, nrhs ).EQ.0 )
THEN
308 bignum = one / smlnum
309 CALL
dlabad( smlnum, bignum )
313 anrm =
zlange(
'M', m, n, a, lda, rwork )
315 IF( anrm.GT.zero .AND. anrm.LT.smlnum )
THEN
319 CALL
zlascl(
'G', 0, 0, anrm, smlnum, m, n, a, lda, info )
321 ELSE IF( anrm.GT.bignum )
THEN
325 CALL
zlascl(
'G', 0, 0, anrm, bignum, m, n, a, lda, info )
327 ELSE IF( anrm.EQ.zero )
THEN
331 CALL
zlaset(
'F', max( m, n ), nrhs, czero, czero,
b, ldb )
336 bnrm =
zlange(
'M', m, nrhs,
b, ldb, rwork )
338 IF( bnrm.GT.zero .AND. bnrm.LT.smlnum )
THEN
342 CALL
zlascl(
'G', 0, 0, bnrm, smlnum, m, nrhs,
b, ldb, info )
344 ELSE IF( bnrm.GT.bignum )
THEN
348 CALL
zlascl(
'G', 0, 0, bnrm, bignum, m, nrhs,
b, ldb, info )
355 CALL
zgeqp3( m, n, a, lda, jpvt, work( 1 ), work( mn+1 ),
356 $ lwork-mn, rwork, info )
357 wsize = mn + dble( work( mn+1 ) )
366 smax = abs( a( 1, 1 ) )
368 IF( abs( a( 1, 1 ) ).EQ.zero )
THEN
370 CALL
zlaset(
'F', max( m, n ), nrhs, czero, czero,
b, ldb )
377 IF( rank.LT.mn )
THEN
379 CALL
zlaic1( imin, rank, work( ismin ), smin, a( 1, i ),
380 $ a( i, i ), sminpr, s1, c1 )
381 CALL
zlaic1( 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
ztzrzf( rank, n, a, lda, work( mn+1 ), work( 2*mn+1 ),
415 CALL
zunmqr(
'Left',
'Conjugate transpose', m, nrhs, mn, a, lda,
416 $ work( 1 ),
b, ldb, work( 2*mn+1 ), lwork-2*mn, info )
417 wsize = max( wsize, 2*mn+dble( work( 2*mn+1 ) ) )
423 CALL
ztrsm(
'Left',
'Upper',
'No transpose',
'Non-unit', rank,
424 $ nrhs, cone, a, lda,
b, ldb )
427 DO 30 i = rank + 1, n
435 CALL
zunmrz(
'Left',
'Conjugate transpose', n, nrhs, rank,
436 $ n-rank, a, lda, work( mn+1 ),
b, ldb,
437 $ work( 2*mn+1 ), lwork-2*mn, info )
446 work( jpvt( i ) ) =
b( i,
j )
448 CALL
zcopy( n, work( 1 ), 1,
b( 1,
j ), 1 )
455 IF( iascl.EQ.1 )
THEN
456 CALL
zlascl(
'G', 0, 0, anrm, smlnum, n, nrhs,
b, ldb, info )
457 CALL
zlascl(
'U', 0, 0, smlnum, anrm, rank, rank, a, lda,
459 ELSE IF( iascl.EQ.2 )
THEN
460 CALL
zlascl(
'G', 0, 0, anrm, bignum, n, nrhs,
b, ldb, info )
461 CALL
zlascl(
'U', 0, 0, bignum, anrm, rank, rank, a, lda,
464 IF( ibscl.EQ.1 )
THEN
465 CALL
zlascl(
'G', 0, 0, smlnum, bnrm, n, nrhs,
b, ldb, info )
466 ELSE IF( ibscl.EQ.2 )
THEN
467 CALL
zlascl(
'G', 0, 0, bignum, bnrm, n, nrhs,
b, ldb, info )
471 work( 1 ) = dcmplx( lwkopt )
subroutine ztrsm(SIDE, UPLO, TRANSA, DIAG, M, N, ALPHA, A, LDA, B, LDB)
ZTRSM
subroutine zlaic1(JOB, J, X, SEST, W, GAMMA, SESTPR, S, C)
ZLAIC1 applies one step of incremental condition estimation.
subroutine zgeqp3(M, N, A, LDA, JPVT, TAU, WORK, LWORK, RWORK, INFO)
ZGEQP3
subroutine zunmrz(SIDE, TRANS, M, N, K, L, A, LDA, TAU, C, LDC, WORK, LWORK, INFO)
ZUNMRZ
DOUBLE PRECISION function zlange(NORM, M, N, A, LDA, WORK)
ZLANGE returns the value of the 1-norm, Frobenius norm, infinity-norm, or the largest absolute value ...
subroutine zgelsy(M, N, NRHS, A, LDA, B, LDB, JPVT, RCOND, RANK, WORK, LWORK, RWORK, INFO)
ZGELSY solves overdetermined or underdetermined systems for GE matrices
subroutine xerbla(SRNAME, INFO)
XERBLA
subroutine dlabad(SMALL, LARGE)
DLABAD
set ue cd $ADTTMP cat<< EOF > tmp f Program LinearEquations Implicit none Real b(3) integer i
subroutine zlascl(TYPE, KL, KU, CFROM, CTO, M, N, A, LDA, INFO)
ZLASCL multiplies a general rectangular matrix by a real scalar defined as cto/cfrom.
subroutine zcopy(N, ZX, INCX, ZY, INCY)
ZCOPY
subroutine zlaset(UPLO, M, N, ALPHA, BETA, A, LDA)
ZLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values...
INTEGER function ilaenv(ISPEC, NAME, OPTS, N1, N2, N3, N4)
set ue cd $ADTTMP cat<< EOF > tmp f Program LinearEquations Implicit none Real j
subroutine ztzrzf(M, N, A, LDA, TAU, WORK, LWORK, INFO)
ZTZRZF
DOUBLE PRECISION function dlamch(CMACH)
DLAMCH
subroutine zunmqr(SIDE, TRANS, M, N, K, A, LDA, TAU, C, LDC, WORK, LWORK, INFO)
ZUNMQR