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
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, dum(1), dum(1),
289 $ dum(1), dum(1), -1, info )
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, dum(1), dum(1),
319 $ dum(1), dum(1), -1, info )
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, dum(1), dum(1),
350 $ dum(1), dum(1), -1, info )
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
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 zungbr(VECT, M, N, K, A, LDA, TAU, WORK, LWORK, INFO)
ZUNGBR
subroutine zbdsqr(UPLO, N, NCVT, NRU, NCC, D, E, VT, LDVT, U, LDU, C, LDC, RWORK, INFO)
ZBDSQR
subroutine zgemv(TRANS, M, N, ALPHA, A, LDA, X, INCX, BETA, Y, INCY)
ZGEMV
subroutine zlacpy(UPLO, M, N, A, LDA, B, LDB)
ZLACPY copies all or part of one two-dimensional array to another.
subroutine zunmlq(SIDE, TRANS, M, N, K, A, LDA, TAU, C, LDC, WORK, LWORK, INFO)
ZUNMLQ
subroutine zgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
ZGEMM
subroutine xerbla(SRNAME, INFO)
XERBLA
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 dlabad(SMALL, LARGE)
DLABAD
set ue cd $ADTTMP cat<< EOF > tmp f Program LinearEquations Implicit none Real b(3) integer i
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 zunmbr(VECT, SIDE, TRANS, M, N, K, A, LDA, TAU, C, LDC, WORK, LWORK, INFO)
ZUNMBR
subroutine zcopy(N, ZX, INCX, ZY, INCY)
ZCOPY
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 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 zgelqf(M, N, A, LDA, TAU, WORK, LWORK, INFO)
ZGELQF
subroutine zgeqrf(M, N, A, LDA, TAU, WORK, LWORK, INFO)
ZGEQRF VARIANT: left-looking Level 3 BLAS of the algorithm.
subroutine zdrscl(N, SA, SX, INCX)
ZDRSCL multiplies a vector by the reciprocal of a real scalar.
double precision function zlange(NORM, M, N, A, LDA, WORK)
ZLANGE returns the value of the 1-norm, Frobenius norm, infinity-norm, or the largest absolute value ...
double precision function dlamch(CMACH)
DLAMCH
subroutine zunmqr(SIDE, TRANS, M, N, K, A, LDA, TAU, C, LDC, WORK, LWORK, INFO)
ZUNMQR
integer function ilaenv(ISPEC, NAME, OPTS, N1, N2, N3, N4)
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 zgebrd(M, N, A, LDA, D, E, TAUQ, TAUP, WORK, LWORK, INFO)
ZGEBRD