178 SUBROUTINE cgelss( M, N, NRHS, A, LDA, B, LDB, S, RCOND, RANK,
179 $ work, lwork, rwork, info )
187 INTEGER info, lda, ldb, lwork, m, n, nrhs, rank
191 REAL rwork( * ), s( * )
192 COMPLEX a( lda, * ),
b( ldb, * ), work( * )
199 parameter( zero = 0.0e+0, one = 1.0e+0 )
201 parameter( czero = ( 0.0e+0, 0.0e+0 ),
202 $ cone = ( 1.0e+0, 0.0e+0 ) )
206 INTEGER bl, chunk, i, iascl, ibscl, ie, il, irwork,
207 $ itau, itaup, itauq, iwork, ldwork, maxmn,
208 $ maxwrk, minmn, minwrk, mm, mnthr
209 INTEGER lwork_cgeqrf, lwork_cunmqr, lwork_cgebrd,
210 $ lwork_cunmbr, lwork_cungbr, lwork_cunmlq,
212 REAL anrm, bignum, bnrm, eps, sfmin, smlnum, thr
238 lquery = ( lwork.EQ.-1 )
241 ELSE IF( n.LT.0 )
THEN
243 ELSE IF( nrhs.LT.0 )
THEN
245 ELSE IF( lda.LT.max( 1, m ) )
THEN
247 ELSE IF( ldb.LT.max( 1, maxmn ) )
THEN
262 IF( minmn.GT.0 )
THEN
264 mnthr =
ilaenv( 6,
'CGELSS',
' ', m, n, nrhs, -1 )
265 IF( m.GE.n .AND. m.GE.mnthr )
THEN
271 CALL
cgeqrf( m, n, a, lda, dum(1), dum(1), -1, info )
274 CALL
cunmqr(
'L',
'C', m, nrhs, n, a, lda, dum(1),
b,
275 $ ldb, dum(1), -1, info )
278 maxwrk = max( maxwrk, n + n*
ilaenv( 1,
'CGEQRF',
' ', m,
280 maxwrk = max( maxwrk, n + nrhs*
ilaenv( 1,
'CUNMQR',
'LC',
288 CALL
cgebrd( mm, n, a, lda, s, dum(1), dum(1),
289 $ dum(1), dum(1), -1, info )
292 CALL
cunmbr(
'Q',
'L',
'C', mm, nrhs, n, a, lda, dum(1),
293 $
b, ldb, dum(1), -1, info )
296 CALL
cungbr(
'P', n, n, n, a, lda, dum(1),
300 maxwrk = max( maxwrk, 2*n + lwork_cgebrd )
301 maxwrk = max( maxwrk, 2*n + lwork_cunmbr )
302 maxwrk = max( maxwrk, 2*n + lwork_cungbr )
303 maxwrk = max( maxwrk, n*nrhs )
304 minwrk = 2*n + max( nrhs, m )
307 minwrk = 2*m + max( nrhs, n )
308 IF( n.GE.mnthr )
THEN
314 CALL
cgelqf( m, n, a, lda, dum(1), dum(1),
318 CALL
cgebrd( m, m, a, lda, s, dum(1), dum(1),
319 $ dum(1), dum(1), -1, info )
322 CALL
cunmbr(
'Q',
'L',
'C', m, nrhs, n, a, lda,
323 $ dum(1),
b, ldb, dum(1), -1, info )
326 CALL
cungbr(
'P', m, m, m, a, lda, dum(1),
330 CALL
cunmlq(
'L',
'C', n, nrhs, m, a, lda, dum(1),
331 $
b, ldb, dum(1), -1, info )
334 maxwrk = m + lwork_cgelqf
335 maxwrk = max( maxwrk, 3*m + m*m + lwork_cgebrd )
336 maxwrk = max( maxwrk, 3*m + m*m + lwork_cunmbr )
337 maxwrk = max( maxwrk, 3*m + m*m + lwork_cungbr )
339 maxwrk = max( maxwrk, m*m + m + m*nrhs )
341 maxwrk = max( maxwrk, m*m + 2*m )
343 maxwrk = max( maxwrk, m + lwork_cunmlq )
349 CALL
cgebrd( m, n, a, lda, s, dum(1), dum(1),
350 $ dum(1), dum(1), -1, info )
353 CALL
cunmbr(
'Q',
'L',
'C', m, nrhs, m, a, lda,
354 $ dum(1),
b, ldb, dum(1), -1, info )
357 CALL
cungbr(
'P', m, n, m, a, lda, dum(1),
360 maxwrk = 2*m + lwork_cgebrd
361 maxwrk = max( maxwrk, 2*m + lwork_cunmbr )
362 maxwrk = max( maxwrk, 2*m + lwork_cungbr )
363 maxwrk = max( maxwrk, n*nrhs )
366 maxwrk = max( minwrk, maxwrk )
370 IF( lwork.LT.minwrk .AND. .NOT.lquery )
375 CALL
xerbla(
'CGELSS', -info )
377 ELSE IF( lquery )
THEN
383 IF( m.EQ.0 .OR. n.EQ.0 )
THEN
393 bignum = one / smlnum
394 CALL
slabad( smlnum, bignum )
398 anrm =
clange(
'M', m, n, a, lda, rwork )
400 IF( anrm.GT.zero .AND. anrm.LT.smlnum )
THEN
404 CALL
clascl(
'G', 0, 0, anrm, smlnum, m, n, a, lda, info )
406 ELSE IF( anrm.GT.bignum )
THEN
410 CALL
clascl(
'G', 0, 0, anrm, bignum, m, n, a, lda, info )
412 ELSE IF( anrm.EQ.zero )
THEN
416 CALL
claset(
'F', max( m, n ), nrhs, czero, czero,
b, ldb )
417 CALL
slaset(
'F', minmn, 1, zero, zero, s, minmn )
424 bnrm =
clange(
'M', m, nrhs,
b, ldb, rwork )
426 IF( bnrm.GT.zero .AND. bnrm.LT.smlnum )
THEN
430 CALL
clascl(
'G', 0, 0, bnrm, smlnum, m, nrhs,
b, ldb, info )
432 ELSE IF( bnrm.GT.bignum )
THEN
436 CALL
clascl(
'G', 0, 0, bnrm, bignum, m, nrhs,
b, ldb, info )
447 IF( m.GE.mnthr )
THEN
459 CALL
cgeqrf( m, n, a, lda, work( itau ), work( iwork ),
460 $ lwork-iwork+1, info )
466 CALL
cunmqr(
'L',
'C', m, nrhs, n, a, lda, work( itau ),
b,
467 $ ldb, work( iwork ), lwork-iwork+1, info )
472 $ CALL
claset(
'L', n-1, n-1, czero, czero, a( 2, 1 ),
485 CALL
cgebrd( mm, n, a, lda, s, rwork( ie ), work( itauq ),
486 $ work( itaup ), work( iwork ), lwork-iwork+1,
493 CALL
cunmbr(
'Q',
'L',
'C', mm, nrhs, n, a, lda, work( itauq ),
494 $
b, ldb, work( iwork ), lwork-iwork+1, info )
500 CALL
cungbr(
'P', n, n, n, a, lda, work( itaup ),
501 $ work( iwork ), lwork-iwork+1, info )
510 CALL
cbdsqr(
'U', n, n, 0, nrhs, s, rwork( ie ), a, lda, dum,
511 $ 1,
b, ldb, rwork( irwork ), info )
517 thr = max( rcond*s( 1 ), sfmin )
519 $ thr = max( eps*s( 1 ), sfmin )
522 IF( s( i ).GT.thr )
THEN
523 CALL
csrscl( nrhs, s( i ),
b( i, 1 ), ldb )
526 CALL
claset(
'F', 1, nrhs, czero, czero,
b( i, 1 ), ldb )
534 IF( lwork.GE.ldb*nrhs .AND. nrhs.GT.1 )
THEN
535 CALL
cgemm(
'C',
'N', n, nrhs, n, cone, a, lda,
b, ldb,
537 CALL
clacpy(
'G', n, nrhs, work, ldb,
b, ldb )
538 ELSE IF( nrhs.GT.1 )
THEN
540 DO 20 i = 1, nrhs, chunk
541 bl = min( nrhs-i+1, chunk )
542 CALL
cgemm(
'C',
'N', n, bl, n, cone, a, lda,
b( 1, i ),
543 $ ldb, czero, work, n )
544 CALL
clacpy(
'G', n, bl, work, n,
b( 1, i ), ldb )
547 CALL
cgemv(
'C', n, n, cone, a, lda,
b, 1, czero, work, 1 )
548 CALL
ccopy( n, work, 1,
b, 1 )
551 ELSE IF( n.GE.mnthr .AND. lwork.GE.3*m+m*m+max( m, nrhs, n-2*m ) )
560 IF( lwork.GE.3*m+m*lda+max( m, nrhs, n-2*m ) )
569 CALL
cgelqf( m, n, a, lda, work( itau ), work( iwork ),
570 $ lwork-iwork+1, info )
575 CALL
clacpy(
'L', m, m, a, lda, work( il ), ldwork )
576 CALL
claset(
'U', m-1, m-1, czero, czero, work( il+ldwork ),
579 itauq = il + ldwork*m
587 CALL
cgebrd( m, m, work( il ), ldwork, s, rwork( ie ),
588 $ work( itauq ), work( itaup ), work( iwork ),
589 $ lwork-iwork+1, info )
595 CALL
cunmbr(
'Q',
'L',
'C', m, nrhs, m, work( il ), ldwork,
596 $ work( itauq ),
b, ldb, work( iwork ),
597 $ lwork-iwork+1, info )
603 CALL
cungbr(
'P', m, m, m, work( il ), ldwork, work( itaup ),
604 $ work( iwork ), lwork-iwork+1, info )
613 CALL
cbdsqr(
'U', m, m, 0, nrhs, s, rwork( ie ), work( il ),
614 $ ldwork, a, lda,
b, ldb, rwork( irwork ), info )
620 thr = max( rcond*s( 1 ), sfmin )
622 $ thr = max( eps*s( 1 ), sfmin )
625 IF( s( i ).GT.thr )
THEN
626 CALL
csrscl( nrhs, s( i ),
b( i, 1 ), ldb )
629 CALL
claset(
'F', 1, nrhs, czero, czero,
b( i, 1 ), ldb )
632 iwork = il + m*ldwork
638 IF( lwork.GE.ldb*nrhs+iwork-1 .AND. nrhs.GT.1 )
THEN
639 CALL
cgemm(
'C',
'N', m, nrhs, m, cone, work( il ), ldwork,
640 $
b, ldb, czero, work( iwork ), ldb )
641 CALL
clacpy(
'G', m, nrhs, work( iwork ), ldb,
b, ldb )
642 ELSE IF( nrhs.GT.1 )
THEN
643 chunk = ( lwork-iwork+1 ) / m
644 DO 40 i = 1, nrhs, chunk
645 bl = min( nrhs-i+1, chunk )
646 CALL
cgemm(
'C',
'N', m, bl, m, cone, work( il ), ldwork,
647 $
b( 1, i ), ldb, czero, work( iwork ), m )
648 CALL
clacpy(
'G', m, bl, work( iwork ), m,
b( 1, i ),
652 CALL
cgemv(
'C', m, m, cone, work( il ), ldwork,
b( 1, 1 ),
653 $ 1, czero, work( iwork ), 1 )
654 CALL
ccopy( m, work( iwork ), 1,
b( 1, 1 ), 1 )
659 CALL
claset(
'F', n-m, nrhs, czero, czero,
b( m+1, 1 ), ldb )
666 CALL
cunmlq(
'L',
'C', n, nrhs, m, a, lda, work( itau ),
b,
667 $ ldb, work( iwork ), lwork-iwork+1, info )
682 CALL
cgebrd( m, n, a, lda, s, rwork( ie ), work( itauq ),
683 $ work( itaup ), work( iwork ), lwork-iwork+1,
690 CALL
cunmbr(
'Q',
'L',
'C', m, nrhs, n, a, lda, work( itauq ),
691 $
b, ldb, work( iwork ), lwork-iwork+1, info )
697 CALL
cungbr(
'P', m, n, m, a, lda, work( itaup ),
698 $ work( iwork ), lwork-iwork+1, info )
707 CALL
cbdsqr(
'L', m, n, 0, nrhs, s, rwork( ie ), a, lda, dum,
708 $ 1,
b, ldb, rwork( irwork ), info )
714 thr = max( rcond*s( 1 ), sfmin )
716 $ thr = max( eps*s( 1 ), sfmin )
719 IF( s( i ).GT.thr )
THEN
720 CALL
csrscl( nrhs, s( i ),
b( i, 1 ), ldb )
723 CALL
claset(
'F', 1, nrhs, czero, czero,
b( i, 1 ), ldb )
731 IF( lwork.GE.ldb*nrhs .AND. nrhs.GT.1 )
THEN
732 CALL
cgemm(
'C',
'N', n, nrhs, m, cone, a, lda,
b, ldb,
734 CALL
clacpy(
'G', n, nrhs, work, ldb,
b, ldb )
735 ELSE IF( nrhs.GT.1 )
THEN
737 DO 60 i = 1, nrhs, chunk
738 bl = min( nrhs-i+1, chunk )
739 CALL
cgemm(
'C',
'N', n, bl, m, cone, a, lda,
b( 1, i ),
740 $ ldb, czero, work, n )
741 CALL
clacpy(
'F', n, bl, work, n,
b( 1, i ), ldb )
744 CALL
cgemv(
'C', m, n, cone, a, lda,
b, 1, czero, work, 1 )
745 CALL
ccopy( n, work, 1,
b, 1 )
751 IF( iascl.EQ.1 )
THEN
752 CALL
clascl(
'G', 0, 0, anrm, smlnum, n, nrhs,
b, ldb, info )
753 CALL
slascl(
'G', 0, 0, smlnum, anrm, minmn, 1, s, minmn,
755 ELSE IF( iascl.EQ.2 )
THEN
756 CALL
clascl(
'G', 0, 0, anrm, bignum, n, nrhs,
b, ldb, info )
757 CALL
slascl(
'G', 0, 0, bignum, anrm, minmn, 1, s, minmn,
760 IF( ibscl.EQ.1 )
THEN
761 CALL
clascl(
'G', 0, 0, smlnum, bnrm, n, nrhs,
b, ldb, info )
762 ELSE IF( ibscl.EQ.2 )
THEN
763 CALL
clascl(
'G', 0, 0, bignum, bnrm, n, nrhs,
b, ldb, info )
subroutine cungbr(VECT, M, N, K, A, LDA, TAU, WORK, LWORK, INFO)
CUNGBR
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 slaset(UPLO, M, N, ALPHA, BETA, A, LDA)
SLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values...
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 ...
subroutine cunmbr(VECT, SIDE, TRANS, M, N, K, A, LDA, TAU, C, LDC, WORK, LWORK, INFO)
CUNMBR
subroutine cgebrd(M, N, A, LDA, D, E, TAUQ, TAUP, WORK, LWORK, INFO)
CGEBRD
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
subroutine csrscl(N, SA, SX, INCX)
CSRSCL multiplies a vector by the reciprocal of a real scalar.
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 cgemv(TRANS, M, N, ALPHA, A, LDA, X, INCX, BETA, Y, INCY)
CGEMV
subroutine clacpy(UPLO, M, N, A, LDA, B, LDB)
CLACPY copies all or part of one two-dimensional array to another.
subroutine cgelss(M, N, NRHS, A, LDA, B, LDB, S, RCOND, RANK, WORK, LWORK, RWORK, INFO)
CGELSS solves overdetermined or underdetermined systems for GE matrices
subroutine cgelqf(M, N, A, LDA, TAU, WORK, LWORK, INFO)
CGELQF
real function slamch(CMACH)
SLAMCH
subroutine cunmlq(SIDE, TRANS, M, N, K, A, LDA, TAU, C, LDC, WORK, LWORK, INFO)
CUNMLQ
subroutine cgeqrf(M, N, A, LDA, TAU, WORK, LWORK, INFO)
CGEQRF
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.
subroutine slabad(SMALL, LARGE)
SLABAD
subroutine ccopy(N, CX, INCX, CY, INCY)
CCOPY
integer function ilaenv(ISPEC, NAME, OPTS, N1, N2, N3, N4)
subroutine cgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
CGEMM
subroutine cbdsqr(UPLO, N, NCVT, NRU, NCC, D, E, VT, LDVT, U, LDU, C, LDC, RWORK, INFO)
CBDSQR