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, s, dum(1), dum(1), dum(1),
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, s, dum(1), dum(1),
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, s, dum(1), dum(1),
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 cbdsqr(UPLO, N, NCVT, NRU, NCC, D, E, VT, LDVT, U, LDU, C, LDC, RWORK, INFO)
CBDSQR
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 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 csrscl(N, SA, SX, INCX)
CSRSCL multiplies a vector by the reciprocal of a real scalar.
integer function ilaenv(ISPEC, NAME, OPTS, N1, N2, N3, N4)
ILAENV
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 cgelqf(M, N, A, LDA, TAU, WORK, LWORK, INFO)
CGELQF
subroutine cgemv(TRANS, M, N, ALPHA, A, LDA, X, INCX, BETA, Y, INCY)
CGEMV
subroutine xerbla(SRNAME, INFO)
XERBLA
subroutine cgeqrf(M, N, A, LDA, TAU, WORK, LWORK, INFO)
CGEQRF
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...
subroutine cgebrd(M, N, A, LDA, D, E, TAUQ, TAUP, WORK, LWORK, INFO)
CGEBRD
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.
real function slamch(CMACH)
SLAMCH
subroutine clacpy(UPLO, M, N, A, LDA, B, LDB)
CLACPY copies all or part of one two-dimensional array to another.
subroutine slabad(SMALL, LARGE)
SLABAD
subroutine cunmqr(SIDE, TRANS, M, N, K, A, LDA, TAU, C, LDC, WORK, LWORK, INFO)
CUNMQR
subroutine ccopy(N, CX, INCX, CY, INCY)
CCOPY
subroutine cunmbr(VECT, SIDE, TRANS, M, N, K, A, LDA, TAU, C, LDC, WORK, LWORK, INFO)
CUNMBR
subroutine cgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
CGEMM
subroutine cunmlq(SIDE, TRANS, M, N, K, A, LDA, TAU, C, LDC, WORK, LWORK, INFO)
CUNMLQ