180 SUBROUTINE zgelsx( M, N, NRHS, A, LDA, B, LDB, JPVT, RCOND,
181 $ RANK, WORK, RWORK, INFO )
188 INTEGER INFO, LDA, LDB, M, N, NRHS, RANK
189 DOUBLE PRECISION RCOND
193 DOUBLE PRECISION RWORK( * )
194 COMPLEX*16 A( LDA, * ), B( LDB, * ), WORK( * )
201 parameter( imax = 1, imin = 2 )
202 DOUBLE PRECISION ZERO, ONE, DONE, NTDONE
203 parameter( zero = 0.0d+0, one = 1.0d+0, done = zero,
205 COMPLEX*16 CZERO, CONE
206 parameter( czero = ( 0.0d+0, 0.0d+0 ),
207 $ cone = ( 1.0d+0, 0.0d+0 ) )
210 INTEGER I, IASCL, IBSCL, ISMAX, ISMIN, J, K, MN
211 DOUBLE PRECISION ANRM, BIGNUM, BNRM, SMAX, SMAXPR, SMIN, SMINPR,
213 COMPLEX*16 C1, C2, S1, S2, T1, T2
220 DOUBLE PRECISION DLAMCH, ZLANGE
221 EXTERNAL dlamch, zlange
224 INTRINSIC abs, dconjg, max, min
237 ELSE IF( n.LT.0 )
THEN
239 ELSE IF( nrhs.LT.0 )
THEN
241 ELSE IF( lda.LT.max( 1, m ) )
THEN
243 ELSE IF( ldb.LT.max( 1, m, n ) )
THEN
248 CALL xerbla(
'ZGELSX', -info )
254 IF( min( m, n, nrhs ).EQ.0 )
THEN
261 smlnum = dlamch(
'S' ) / dlamch(
'P' )
262 bignum = one / smlnum
266 anrm = zlange(
'M', m, n, a, lda, rwork )
268 IF( anrm.GT.zero .AND. anrm.LT.smlnum )
THEN
272 CALL zlascl(
'G', 0, 0, anrm, smlnum, m, n, a, lda, info )
274 ELSE IF( anrm.GT.bignum )
THEN
278 CALL zlascl(
'G', 0, 0, anrm, bignum, m, n, a, lda, info )
280 ELSE IF( anrm.EQ.zero )
THEN
284 CALL zlaset(
'F', max( m, n ), nrhs, czero, czero, b, ldb )
289 bnrm = zlange(
'M', m, nrhs, b, ldb, rwork )
291 IF( bnrm.GT.zero .AND. bnrm.LT.smlnum )
THEN
295 CALL zlascl(
'G', 0, 0, bnrm, smlnum, m, nrhs, b, ldb,
298 ELSE IF( bnrm.GT.bignum )
THEN
302 CALL zlascl(
'G', 0, 0, bnrm, bignum, m, nrhs, b, ldb,
310 CALL zgeqpf( m, n, a, lda, jpvt, work( 1 ), work( mn+1 ),
320 smax = abs( a( 1, 1 ) )
322 IF( abs( a( 1, 1 ) ).EQ.zero )
THEN
324 CALL zlaset(
'F', max( m, n ), nrhs, czero, czero, b, ldb )
331 IF( rank.LT.mn )
THEN
333 CALL zlaic1( imin, rank, work( ismin ), smin, a( 1, i ),
334 $ a( i, i ), sminpr, s1, c1 )
335 CALL zlaic1( imax, rank, work( ismax ), smax, a( 1, i ),
336 $ a( i, i ), smaxpr, s2, c2 )
338 IF( smaxpr*rcond.LE.sminpr )
THEN
340 work( ismin+i-1 ) = s1*work( ismin+i-1 )
341 work( ismax+i-1 ) = s2*work( ismax+i-1 )
343 work( ismin+rank ) = c1
344 work( ismax+rank ) = c2
359 $
CALL ztzrqf( rank, n, a, lda, work( mn+1 ), info )
365 CALL zunm2r(
'Left',
'Conjugate transpose', m, nrhs, mn, a,
366 $ lda, work( 1 ), b, ldb, work( 2*mn+1 ), info )
372 CALL ztrsm(
'Left',
'Upper',
'No transpose',
'Non-unit', rank,
373 $ nrhs, cone, a, lda, b, ldb )
375 DO 40 i = rank + 1, n
385 CALL zlatzm(
'Left', n-rank+1, nrhs, a( i, rank+1 ),
386 $ lda, dconjg( work( mn+i ) ), b( i, 1 ),
387 $ b( rank+1, 1 ), ldb, work( 2*mn+1 ) )
397 work( 2*mn+i ) = ntdone
400 IF( work( 2*mn+i ).EQ.ntdone )
THEN
401 IF( jpvt( i ).NE.i )
THEN
404 t2 = b( jpvt( k ), j )
406 b( jpvt( k ), j ) = t1
407 work( 2*mn+k ) = done
410 t2 = b( jpvt( k ), j )
414 work( 2*mn+k ) = done
422 IF( iascl.EQ.1 )
THEN
423 CALL zlascl(
'G', 0, 0, anrm, smlnum, n, nrhs, b, ldb,
425 CALL zlascl(
'U', 0, 0, smlnum, anrm, rank, rank, a, lda,
427 ELSE IF( iascl.EQ.2 )
THEN
428 CALL zlascl(
'G', 0, 0, anrm, bignum, n, nrhs, b, ldb,
430 CALL zlascl(
'U', 0, 0, bignum, anrm, rank, rank, a, lda,
433 IF( ibscl.EQ.1 )
THEN
434 CALL zlascl(
'G', 0, 0, smlnum, bnrm, n, nrhs, b, ldb,
436 ELSE IF( ibscl.EQ.2 )
THEN
437 CALL zlascl(
'G', 0, 0, bignum, bnrm, n, nrhs, b, ldb,