178 SUBROUTINE zgelss( 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
188 DOUBLE PRECISION RCOND
191 DOUBLE PRECISION RWORK( * ), S( * )
192 COMPLEX*16 A( lda, * ), B( ldb, * ), WORK( * )
198 DOUBLE PRECISION ZERO, ONE
199 parameter ( zero = 0.0d+0, one = 1.0d+0 )
200 COMPLEX*16 CZERO, CONE
201 parameter ( czero = ( 0.0d+0, 0.0d+0 ),
202 $ cone = ( 1.0d+0, 0.0d+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_ZGEQRF, LWORK_ZUNMQR, LWORK_ZGEBRD,
210 $ lwork_zunmbr, lwork_zungbr, lwork_zunmlq,
212 DOUBLE PRECISION ANRM, BIGNUM, BNRM, EPS, SFMIN, SMLNUM, THR
225 DOUBLE PRECISION DLAMCH, ZLANGE
226 EXTERNAL ilaenv, dlamch, zlange
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,
'ZGELSS',
' ', m, n, nrhs, -1 )
265 IF( m.GE.n .AND. m.GE.mnthr )
THEN
271 CALL zgeqrf( m, n, a, lda, dum(1), dum(1), -1, info )
274 CALL zunmqr(
'L',
'C', m, nrhs, n, a, lda, dum(1), b,
275 $ ldb, dum(1), -1, info )
278 maxwrk = max( maxwrk, n + n*ilaenv( 1,
'ZGEQRF',
' ', m,
280 maxwrk = max( maxwrk, n + nrhs*ilaenv( 1,
'ZUNMQR',
'LC',
288 CALL zgebrd( mm, n, a, lda, s, s, dum(1), dum(1), dum(1),
292 CALL zunmbr(
'Q',
'L',
'C', mm, nrhs, n, a, lda, dum(1),
293 $ b, ldb, dum(1), -1, info )
296 CALL zungbr(
'P', n, n, n, a, lda, dum(1),
300 maxwrk = max( maxwrk, 2*n + lwork_zgebrd )
301 maxwrk = max( maxwrk, 2*n + lwork_zunmbr )
302 maxwrk = max( maxwrk, 2*n + lwork_zungbr )
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 zgelqf( m, n, a, lda, dum(1), dum(1),
318 CALL zgebrd( m, m, a, lda, s, s, dum(1), dum(1),
322 CALL zunmbr(
'Q',
'L',
'C', m, nrhs, n, a, lda,
323 $ dum(1), b, ldb, dum(1), -1, info )
326 CALL zungbr(
'P', m, m, m, a, lda, dum(1),
330 CALL zunmlq(
'L',
'C', n, nrhs, m, a, lda, dum(1),
331 $ b, ldb, dum(1), -1, info )
334 maxwrk = m + lwork_zgelqf
335 maxwrk = max( maxwrk, 3*m + m*m + lwork_zgebrd )
336 maxwrk = max( maxwrk, 3*m + m*m + lwork_zunmbr )
337 maxwrk = max( maxwrk, 3*m + m*m + lwork_zungbr )
339 maxwrk = max( maxwrk, m*m + m + m*nrhs )
341 maxwrk = max( maxwrk, m*m + 2*m )
343 maxwrk = max( maxwrk, m + lwork_zunmlq )
349 CALL zgebrd( m, n, a, lda, s, s, dum(1), dum(1),
353 CALL zunmbr(
'Q',
'L',
'C', m, nrhs, m, a, lda,
354 $ dum(1), b, ldb, dum(1), -1, info )
357 CALL zungbr(
'P', m, n, m, a, lda, dum(1),
360 maxwrk = 2*m + lwork_zgebrd
361 maxwrk = max( maxwrk, 2*m + lwork_zunmbr )
362 maxwrk = max( maxwrk, 2*m + lwork_zungbr )
363 maxwrk = max( maxwrk, n*nrhs )
366 maxwrk = max( minwrk, maxwrk )
370 IF( lwork.LT.minwrk .AND. .NOT.lquery )
375 CALL xerbla(
'ZGELSS', -info )
377 ELSE IF( lquery )
THEN
383 IF( m.EQ.0 .OR. n.EQ.0 )
THEN
391 sfmin = dlamch(
'S' )
393 bignum = one / smlnum
394 CALL dlabad( smlnum, bignum )
398 anrm = zlange(
'M', m, n, a, lda, rwork )
400 IF( anrm.GT.zero .AND. anrm.LT.smlnum )
THEN
404 CALL zlascl(
'G', 0, 0, anrm, smlnum, m, n, a, lda, info )
406 ELSE IF( anrm.GT.bignum )
THEN
410 CALL zlascl(
'G', 0, 0, anrm, bignum, m, n, a, lda, info )
412 ELSE IF( anrm.EQ.zero )
THEN
416 CALL zlaset(
'F', max( m, n ), nrhs, czero, czero, b, ldb )
417 CALL dlaset(
'F', minmn, 1, zero, zero, s, minmn )
424 bnrm = zlange(
'M', m, nrhs, b, ldb, rwork )
426 IF( bnrm.GT.zero .AND. bnrm.LT.smlnum )
THEN
430 CALL zlascl(
'G', 0, 0, bnrm, smlnum, m, nrhs, b, ldb, info )
432 ELSE IF( bnrm.GT.bignum )
THEN
436 CALL zlascl(
'G', 0, 0, bnrm, bignum, m, nrhs, b, ldb, info )
447 IF( m.GE.mnthr )
THEN
459 CALL zgeqrf( m, n, a, lda, work( itau ), work( iwork ),
460 $ lwork-iwork+1, info )
466 CALL zunmqr(
'L',
'C', m, nrhs, n, a, lda, work( itau ), b,
467 $ ldb, work( iwork ), lwork-iwork+1, info )
472 $
CALL zlaset(
'L', n-1, n-1, czero, czero, a( 2, 1 ),
485 CALL zgebrd( mm, n, a, lda, s, rwork( ie ), work( itauq ),
486 $ work( itaup ), work( iwork ), lwork-iwork+1,
493 CALL zunmbr(
'Q',
'L',
'C', mm, nrhs, n, a, lda, work( itauq ),
494 $ b, ldb, work( iwork ), lwork-iwork+1, info )
500 CALL zungbr(
'P', n, n, n, a, lda, work( itaup ),
501 $ work( iwork ), lwork-iwork+1, info )
510 CALL zbdsqr(
'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 zdrscl( nrhs, s( i ), b( i, 1 ), ldb )
526 CALL zlaset(
'F', 1, nrhs, czero, czero, b( i, 1 ), ldb )
534 IF( lwork.GE.ldb*nrhs .AND. nrhs.GT.1 )
THEN
535 CALL zgemm(
'C',
'N', n, nrhs, n, cone, a, lda, b, ldb,
537 CALL zlacpy(
'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 zgemm(
'C',
'N', n, bl, n, cone, a, lda, b( 1, i ),
543 $ ldb, czero, work, n )
544 CALL zlacpy(
'G', n, bl, work, n, b( 1, i ), ldb )
547 CALL zgemv(
'C', n, n, cone, a, lda, b, 1, czero, work, 1 )
548 CALL zcopy( 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 zgelqf( m, n, a, lda, work( itau ), work( iwork ),
570 $ lwork-iwork+1, info )
575 CALL zlacpy(
'L', m, m, a, lda, work( il ), ldwork )
576 CALL zlaset(
'U', m-1, m-1, czero, czero, work( il+ldwork ),
579 itauq = il + ldwork*m
587 CALL zgebrd( m, m, work( il ), ldwork, s, rwork( ie ),
588 $ work( itauq ), work( itaup ), work( iwork ),
589 $ lwork-iwork+1, info )
595 CALL zunmbr(
'Q',
'L',
'C', m, nrhs, m, work( il ), ldwork,
596 $ work( itauq ), b, ldb, work( iwork ),
597 $ lwork-iwork+1, info )
603 CALL zungbr(
'P', m, m, m, work( il ), ldwork, work( itaup ),
604 $ work( iwork ), lwork-iwork+1, info )
613 CALL zbdsqr(
'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 zdrscl( nrhs, s( i ), b( i, 1 ), ldb )
629 CALL zlaset(
'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 zgemm(
'C',
'N', m, nrhs, m, cone, work( il ), ldwork,
640 $ b, ldb, czero, work( iwork ), ldb )
641 CALL zlacpy(
'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 zgemm(
'C',
'N', m, bl, m, cone, work( il ), ldwork,
647 $ b( 1, i ), ldb, czero, work( iwork ), m )
648 CALL zlacpy(
'G', m, bl, work( iwork ), m, b( 1, i ),
652 CALL zgemv(
'C', m, m, cone, work( il ), ldwork, b( 1, 1 ),
653 $ 1, czero, work( iwork ), 1 )
654 CALL zcopy( m, work( iwork ), 1, b( 1, 1 ), 1 )
659 CALL zlaset(
'F', n-m, nrhs, czero, czero, b( m+1, 1 ), ldb )
666 CALL zunmlq(
'L',
'C', n, nrhs, m, a, lda, work( itau ), b,
667 $ ldb, work( iwork ), lwork-iwork+1, info )
682 CALL zgebrd( m, n, a, lda, s, rwork( ie ), work( itauq ),
683 $ work( itaup ), work( iwork ), lwork-iwork+1,
690 CALL zunmbr(
'Q',
'L',
'C', m, nrhs, n, a, lda, work( itauq ),
691 $ b, ldb, work( iwork ), lwork-iwork+1, info )
697 CALL zungbr(
'P', m, n, m, a, lda, work( itaup ),
698 $ work( iwork ), lwork-iwork+1, info )
707 CALL zbdsqr(
'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 zdrscl( nrhs, s( i ), b( i, 1 ), ldb )
723 CALL zlaset(
'F', 1, nrhs, czero, czero, b( i, 1 ), ldb )
731 IF( lwork.GE.ldb*nrhs .AND. nrhs.GT.1 )
THEN
732 CALL zgemm(
'C',
'N', n, nrhs, m, cone, a, lda, b, ldb,
734 CALL zlacpy(
'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 zgemm(
'C',
'N', n, bl, m, cone, a, lda, b( 1, i ),
740 $ ldb, czero, work, n )
741 CALL zlacpy(
'F', n, bl, work, n, b( 1, i ), ldb )
744 CALL zgemv(
'C', m, n, cone, a, lda, b, 1, czero, work, 1 )
745 CALL zcopy( n, work, 1, b, 1 )
751 IF( iascl.EQ.1 )
THEN
752 CALL zlascl(
'G', 0, 0, anrm, smlnum, n, nrhs, b, ldb, info )
753 CALL dlascl(
'G', 0, 0, smlnum, anrm, minmn, 1, s, minmn,
755 ELSE IF( iascl.EQ.2 )
THEN
756 CALL zlascl(
'G', 0, 0, anrm, bignum, n, nrhs, b, ldb, info )
757 CALL dlascl(
'G', 0, 0, bignum, anrm, minmn, 1, s, minmn,
760 IF( ibscl.EQ.1 )
THEN
761 CALL zlascl(
'G', 0, 0, smlnum, bnrm, n, nrhs, b, ldb, info )
762 ELSE IF( ibscl.EQ.2 )
THEN
763 CALL zlascl(
'G', 0, 0, bignum, bnrm, n, nrhs, b, ldb, info )
subroutine zcopy(N, ZX, INCX, ZY, INCY)
ZCOPY
subroutine zunmbr(VECT, SIDE, TRANS, M, N, K, A, LDA, TAU, C, LDC, WORK, LWORK, INFO)
ZUNMBR
subroutine zungbr(VECT, M, N, K, A, LDA, TAU, WORK, LWORK, INFO)
ZUNGBR
subroutine zgelss(M, N, NRHS, A, LDA, B, LDB, S, RCOND, RANK, WORK, LWORK, RWORK, INFO)
ZGELSS solves overdetermined or underdetermined systems for GE matrices
subroutine zgemv(TRANS, M, N, ALPHA, A, LDA, X, INCX, BETA, Y, INCY)
ZGEMV
subroutine dlascl(TYPE, KL, KU, CFROM, CTO, M, N, A, LDA, INFO)
DLASCL multiplies a general rectangular matrix by a real scalar defined as cto/cfrom.
subroutine zunmlq(SIDE, TRANS, M, N, K, A, LDA, TAU, C, LDC, WORK, LWORK, INFO)
ZUNMLQ
subroutine zlascl(TYPE, KL, KU, CFROM, CTO, M, N, A, LDA, INFO)
ZLASCL multiplies a general rectangular matrix by a real scalar defined as cto/cfrom.
subroutine zlacpy(UPLO, M, N, A, LDA, B, LDB)
ZLACPY copies all or part of one two-dimensional array to another.
subroutine zgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
ZGEMM
subroutine dlaset(UPLO, M, N, ALPHA, BETA, A, LDA)
DLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values...
subroutine zgeqrf(M, N, A, LDA, TAU, WORK, LWORK, INFO)
ZGEQRF VARIANT: left-looking Level 3 BLAS of the algorithm.
subroutine zgebrd(M, N, A, LDA, D, E, TAUQ, TAUP, WORK, LWORK, INFO)
ZGEBRD
subroutine zunmqr(SIDE, TRANS, M, N, K, A, LDA, TAU, C, LDC, WORK, LWORK, INFO)
ZUNMQR
subroutine zdrscl(N, SA, SX, INCX)
ZDRSCL multiplies a vector by the reciprocal of a real scalar.
subroutine xerbla(SRNAME, INFO)
XERBLA
subroutine zlaset(UPLO, M, N, ALPHA, BETA, A, LDA)
ZLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values...
subroutine dlabad(SMALL, LARGE)
DLABAD
subroutine zbdsqr(UPLO, N, NCVT, NRU, NCC, D, E, VT, LDVT, U, LDU, C, LDC, RWORK, INFO)
ZBDSQR
subroutine zgelqf(M, N, A, LDA, TAU, WORK, LWORK, INFO)
ZGELQF