184 SUBROUTINE zgelsx( M, N, NRHS, A, LDA, B, LDB, JPVT, RCOND, RANK,
185 $ WORK, RWORK, INFO )
193 INTEGER INFO, LDA, LDB, M, N, NRHS, RANK
194 DOUBLE PRECISION RCOND
198 DOUBLE PRECISION RWORK( * )
199 COMPLEX*16 A( LDA, * ), B( LDB, * ), WORK( * )
206 parameter( imax = 1, imin = 2 )
207 DOUBLE PRECISION ZERO, ONE, DONE, NTDONE
208 parameter( zero = 0.0d+0, one = 1.0d+0, done = zero,
210 COMPLEX*16 CZERO, CONE
211 parameter( czero = ( 0.0d+0, 0.0d+0 ),
212 $ cone = ( 1.0d+0, 0.0d+0 ) )
215 INTEGER I, IASCL, IBSCL, ISMAX, ISMIN, J, K, MN
216 DOUBLE PRECISION ANRM, BIGNUM, BNRM, SMAX, SMAXPR, SMIN, SMINPR,
218 COMPLEX*16 C1, C2, S1, S2, T1, T2
225 DOUBLE PRECISION DLAMCH, ZLANGE
226 EXTERNAL dlamch, zlange
229 INTRINSIC abs, dconjg, max, min
242 ELSE IF( n.LT.0 )
THEN
244 ELSE IF( nrhs.LT.0 )
THEN
246 ELSE IF( lda.LT.max( 1, m ) )
THEN
248 ELSE IF( ldb.LT.max( 1, m, n ) )
THEN
253 CALL xerbla(
'ZGELSX', -info )
259 IF( min( m, n, nrhs ).EQ.0 )
THEN
266 smlnum = dlamch(
'S' ) / dlamch(
'P' )
267 bignum = one / smlnum
268 CALL dlabad( smlnum, bignum )
272 anrm = zlange(
'M', m, n, a, lda, rwork )
274 IF( anrm.GT.zero .AND. anrm.LT.smlnum )
THEN
278 CALL zlascl(
'G', 0, 0, anrm, smlnum, m, n, a, lda, info )
280 ELSE IF( anrm.GT.bignum )
THEN
284 CALL zlascl(
'G', 0, 0, anrm, bignum, m, n, a, lda, info )
286 ELSE IF( anrm.EQ.zero )
THEN
290 CALL zlaset(
'F', max( m, n ), nrhs, czero, czero, b, ldb )
295 bnrm = zlange(
'M', m, nrhs, b, ldb, rwork )
297 IF( bnrm.GT.zero .AND. bnrm.LT.smlnum )
THEN
301 CALL zlascl(
'G', 0, 0, bnrm, smlnum, m, nrhs, b, ldb, info )
303 ELSE IF( bnrm.GT.bignum )
THEN
307 CALL zlascl(
'G', 0, 0, bnrm, bignum, m, nrhs, b, ldb, info )
314 CALL zgeqpf( m, n, a, lda, jpvt, work( 1 ), work( mn+1 ), rwork,
324 smax = abs( a( 1, 1 ) )
326 IF( abs( a( 1, 1 ) ).EQ.zero )
THEN
328 CALL zlaset(
'F', max( m, n ), nrhs, czero, czero, b, ldb )
335 IF( rank.LT.mn )
THEN
337 CALL zlaic1( imin, rank, work( ismin ), smin, a( 1, i ),
338 $ a( i, i ), sminpr, s1, c1 )
339 CALL zlaic1( imax, rank, work( ismax ), smax, a( 1, i ),
340 $ a( i, i ), smaxpr, s2, c2 )
342 IF( smaxpr*rcond.LE.sminpr )
THEN
344 work( ismin+i-1 ) = s1*work( ismin+i-1 )
345 work( ismax+i-1 ) = s2*work( ismax+i-1 )
347 work( ismin+rank ) = c1
348 work( ismax+rank ) = c2
363 $
CALL ztzrqf( rank, n, a, lda, work( mn+1 ), info )
369 CALL zunm2r(
'Left',
'Conjugate transpose', m, nrhs, mn, a, lda,
370 $ work( 1 ), b, ldb, work( 2*mn+1 ), info )
376 CALL ztrsm(
'Left',
'Upper',
'No transpose',
'Non-unit', rank,
377 $ nrhs, cone, a, lda, b, ldb )
379 DO 40 i = rank + 1, n
389 CALL zlatzm(
'Left', n-rank+1, nrhs, a( i, rank+1 ), lda,
390 $ dconjg( work( mn+i ) ), b( i, 1 ),
391 $ b( rank+1, 1 ), ldb, work( 2*mn+1 ) )
401 work( 2*mn+i ) = ntdone
404 IF( work( 2*mn+i ).EQ.ntdone )
THEN
405 IF( jpvt( i ).NE.i )
THEN
408 t2 = b( jpvt( k ), j )
410 b( jpvt( k ), j ) = t1
411 work( 2*mn+k ) = done
414 t2 = b( jpvt( k ), j )
418 work( 2*mn+k ) = done
426 IF( iascl.EQ.1 )
THEN
427 CALL zlascl(
'G', 0, 0, anrm, smlnum, n, nrhs, b, ldb, info )
428 CALL zlascl(
'U', 0, 0, smlnum, anrm, rank, rank, a, lda,
430 ELSE IF( iascl.EQ.2 )
THEN
431 CALL zlascl(
'G', 0, 0, anrm, bignum, n, nrhs, b, ldb, info )
432 CALL zlascl(
'U', 0, 0, bignum, anrm, rank, rank, a, lda,
435 IF( ibscl.EQ.1 )
THEN
436 CALL zlascl(
'G', 0, 0, smlnum, bnrm, n, nrhs, b, ldb, info )
437 ELSE IF( ibscl.EQ.2 )
THEN
438 CALL zlascl(
'G', 0, 0, bignum, bnrm, n, nrhs, b, ldb, info )