210 SUBROUTINE cgelsy( 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
225 COMPLEX a( lda, * ),
b( ldb, * ), work( * )
232 parameter( imax = 1, imin = 2 )
234 parameter( zero = 0.0e+0, one = 1.0e+0 )
236 parameter( czero = ( 0.0e+0, 0.0e+0 ),
237 $ cone = ( 1.0e+0, 0.0e+0 ) )
241 INTEGER i, iascl, ibscl, ismax, ismin,
j, lwkopt, mn,
242 $ nb, nb1, nb2, nb3, nb4
243 REAL anrm, bignum, bnrm, smax, smaxpr, smin, sminpr,
245 COMPLEX c1, c2, s1, s2
257 INTRINSIC abs, max, min,
REAL, cmplx
268 nb1 =
ilaenv( 1,
'CGEQRF',
' ', m, n, -1, -1 )
269 nb2 =
ilaenv( 1,
'CGERQF',
' ', m, n, -1, -1 )
270 nb3 =
ilaenv( 1,
'CUNMQR',
' ', m, n, nrhs, -1 )
271 nb4 =
ilaenv( 1,
'CUNMRQ',
' ', 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 ) = cmplx( 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.
292 CALL
xerbla(
'CGELSY', -info )
294 ELSE IF( lquery )
THEN
300 IF( min( m, n, nrhs ).EQ.0 )
THEN
308 bignum = one / smlnum
309 CALL
slabad( smlnum, bignum )
313 anrm =
clange(
'M', m, n, a, lda, rwork )
315 IF( anrm.GT.zero .AND. anrm.LT.smlnum )
THEN
319 CALL
clascl(
'G', 0, 0, anrm, smlnum, m, n, a, lda, info )
321 ELSE IF( anrm.GT.bignum )
THEN
325 CALL
clascl(
'G', 0, 0, anrm, bignum, m, n, a, lda, info )
327 ELSE IF( anrm.EQ.zero )
THEN
331 CALL
claset(
'F', max( m, n ), nrhs, czero, czero,
b, ldb )
336 bnrm =
clange(
'M', m, nrhs,
b, ldb, rwork )
338 IF( bnrm.GT.zero .AND. bnrm.LT.smlnum )
THEN
342 CALL
clascl(
'G', 0, 0, bnrm, smlnum, m, nrhs,
b, ldb, info )
344 ELSE IF( bnrm.GT.bignum )
THEN
348 CALL
clascl(
'G', 0, 0, bnrm, bignum, m, nrhs,
b, ldb, info )
355 CALL
cgeqp3( m, n, a, lda, jpvt, work( 1 ), work( mn+1 ),
356 $ lwork-mn, rwork, info )
357 wsize = mn +
REAL( WORK( MN+1 ) )
366 smax = abs( a( 1, 1 ) )
368 IF( abs( a( 1, 1 ) ).EQ.zero )
THEN
370 CALL
claset(
'F', max( m, n ), nrhs, czero, czero,
b, ldb )
377 IF( rank.LT.mn )
THEN
379 CALL
claic1( imin, rank, work( ismin ), smin, a( 1, i ),
380 $ a( i, i ), sminpr, s1, c1 )
381 CALL
claic1( 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
ctzrzf( rank, n, a, lda, work( mn+1 ), work( 2*mn+1 ),
415 CALL
cunmqr(
'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+
REAL( WORK( 2*MN+1 ) ) )
423 CALL
ctrsm(
'Left',
'Upper',
'No transpose',
'Non-unit', rank,
424 $ nrhs, cone, a, lda,
b, ldb )
427 DO 30 i = rank + 1, n
435 CALL
cunmrz(
'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
ccopy( n, work( 1 ), 1,
b( 1,
j ), 1 )
455 IF( iascl.EQ.1 )
THEN
456 CALL
clascl(
'G', 0, 0, anrm, smlnum, n, nrhs,
b, ldb, info )
457 CALL
clascl(
'U', 0, 0, smlnum, anrm, rank, rank, a, lda,
459 ELSE IF( iascl.EQ.2 )
THEN
460 CALL
clascl(
'G', 0, 0, anrm, bignum, n, nrhs,
b, ldb, info )
461 CALL
clascl(
'U', 0, 0, bignum, anrm, rank, rank, a, lda,
464 IF( ibscl.EQ.1 )
THEN
465 CALL
clascl(
'G', 0, 0, smlnum, bnrm, n, nrhs,
b, ldb, info )
466 ELSE IF( ibscl.EQ.2 )
THEN
467 CALL
clascl(
'G', 0, 0, bignum, bnrm, n, nrhs,
b, ldb, info )
471 work( 1 ) = cmplx( lwkopt )
subroutine claset(UPLO, M, N, ALPHA, BETA, A, LDA)
CLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values...
subroutine ctrsm(SIDE, UPLO, TRANSA, DIAG, M, N, ALPHA, A, LDA, B, LDB)
CTRSM
subroutine cgelsy(M, N, NRHS, A, LDA, B, LDB, JPVT, RCOND, RANK, WORK, LWORK, RWORK, INFO)
CGELSY solves overdetermined or underdetermined systems for GE matrices
REAL function slamch(CMACH)
SLAMCH
subroutine cgeqp3(M, N, A, LDA, JPVT, TAU, WORK, LWORK, RWORK, INFO)
CGEQP3
subroutine claic1(JOB, J, X, SEST, W, GAMMA, SESTPR, S, C)
CLAIC1 applies one step of incremental condition estimation.
subroutine xerbla(SRNAME, INFO)
XERBLA
subroutine cunmqr(SIDE, TRANS, M, N, K, A, LDA, TAU, C, LDC, WORK, LWORK, INFO)
CUNMQR
set ue cd $ADTTMP cat<< EOF > tmp f Program LinearEquations Implicit none Real b(3) integer i
REAL function clange(NORM, M, N, A, LDA, WORK)
CLANGE returns the value of the 1-norm, Frobenius norm, infinity-norm, or the largest absolute value ...
INTEGER function ilaenv(ISPEC, NAME, OPTS, N1, N2, N3, N4)
subroutine clascl(TYPE, KL, KU, CFROM, CTO, M, N, A, LDA, INFO)
CLASCL multiplies a general rectangular matrix by a real scalar defined as cto/cfrom.
set ue cd $ADTTMP cat<< EOF > tmp f Program LinearEquations Implicit none Real j
subroutine slabad(SMALL, LARGE)
SLABAD
subroutine ccopy(N, CX, INCX, CY, INCY)
CCOPY
subroutine cunmrz(SIDE, TRANS, M, N, K, L, A, LDA, TAU, C, LDC, WORK, LWORK, INFO)
CUNMRZ
subroutine ctzrzf(M, N, A, LDA, TAU, WORK, LWORK, INFO)
CTZRZF