225 SUBROUTINE zgelsd( M, N, NRHS, A, LDA, B, LDB, S, RCOND, RANK,
226 $ work, lwork, rwork, iwork, info )
234 INTEGER INFO, LDA, LDB, LWORK, M, N, NRHS, RANK
235 DOUBLE PRECISION RCOND
239 DOUBLE PRECISION RWORK( * ), S( * )
240 COMPLEX*16 A( lda, * ), B( ldb, * ), WORK( * )
246 DOUBLE PRECISION ZERO, ONE, TWO
247 parameter( zero = 0.0d+0, one = 1.0d+0, two = 2.0d+0 )
249 parameter( czero = ( 0.0d+0, 0.0d+0 ) )
253 INTEGER IASCL, IBSCL, IE, IL, ITAU, ITAUP, ITAUQ,
254 $ ldwork, liwork, lrwork, maxmn, maxwrk, minmn,
255 $ minwrk, mm, mnthr, nlvl, nrwork, nwork, smlsiz
256 DOUBLE PRECISION ANRM, BIGNUM, BNRM, EPS, SFMIN, SMLNUM
265 DOUBLE PRECISION DLAMCH, ZLANGE
266 EXTERNAL ilaenv, dlamch, zlange
269 INTRINSIC int, log, max, min, dble
278 lquery = ( lwork.EQ.-1 )
281 ELSE IF( n.LT.0 )
THEN
283 ELSE IF( nrhs.LT.0 )
THEN
285 ELSE IF( lda.LT.max( 1, m ) )
THEN
287 ELSE IF( ldb.LT.max( 1, maxmn ) )
THEN
303 IF( minmn.GT.0 )
THEN
304 smlsiz = ilaenv( 9,
'ZGELSD',
' ', 0, 0, 0, 0 )
305 mnthr = ilaenv( 6,
'ZGELSD',
' ', m, n, nrhs, -1 )
306 nlvl = max( int( log( dble( minmn ) / dble( smlsiz + 1 ) ) /
307 $ log( two ) ) + 1, 0 )
308 liwork = 3*minmn*nlvl + 11*minmn
310 IF( m.GE.n .AND. m.GE.mnthr )
THEN
316 maxwrk = max( maxwrk, n*ilaenv( 1,
'ZGEQRF',
' ', m, n,
318 maxwrk = max( maxwrk, nrhs*ilaenv( 1,
'ZUNMQR',
'LC', m,
325 lrwork = 10*n + 2*n*smlsiz + 8*n*nlvl + 3*smlsiz*nrhs +
326 $ max( (smlsiz+1)**2, n*(1+nrhs) + 2*nrhs )
327 maxwrk = max( maxwrk, 2*n + ( mm + n )*ilaenv( 1,
328 $
'ZGEBRD',
' ', mm, n, -1, -1 ) )
329 maxwrk = max( maxwrk, 2*n + nrhs*ilaenv( 1,
'ZUNMBR',
330 $
'QLC', mm, nrhs, n, -1 ) )
331 maxwrk = max( maxwrk, 2*n + ( n - 1 )*ilaenv( 1,
332 $
'ZUNMBR',
'PLN', n, nrhs, n, -1 ) )
333 maxwrk = max( maxwrk, 2*n + n*nrhs )
334 minwrk = max( 2*n + mm, 2*n + n*nrhs )
337 lrwork = 10*m + 2*m*smlsiz + 8*m*nlvl + 3*smlsiz*nrhs +
338 $ max( (smlsiz+1)**2, n*(1+nrhs) + 2*nrhs )
339 IF( n.GE.mnthr )
THEN
344 maxwrk = m + m*ilaenv( 1,
'ZGELQF',
' ', m, n, -1,
346 maxwrk = max( maxwrk, m*m + 4*m + 2*m*ilaenv( 1,
347 $
'ZGEBRD',
' ', m, m, -1, -1 ) )
348 maxwrk = max( maxwrk, m*m + 4*m + nrhs*ilaenv( 1,
349 $
'ZUNMBR',
'QLC', m, nrhs, m, -1 ) )
350 maxwrk = max( maxwrk, m*m + 4*m + ( m - 1 )*ilaenv( 1,
351 $
'ZUNMLQ',
'LC', n, nrhs, m, -1 ) )
353 maxwrk = max( maxwrk, m*m + m + m*nrhs )
355 maxwrk = max( maxwrk, m*m + 2*m )
357 maxwrk = max( maxwrk, m*m + 4*m + m*nrhs )
360 maxwrk = max( maxwrk,
361 $ 4*m+m*m+max( m, 2*m-4, nrhs, n-3*m ) )
366 maxwrk = 2*m + ( n + m )*ilaenv( 1,
'ZGEBRD',
' ', m,
368 maxwrk = max( maxwrk, 2*m + nrhs*ilaenv( 1,
'ZUNMBR',
369 $
'QLC', m, nrhs, m, -1 ) )
370 maxwrk = max( maxwrk, 2*m + m*ilaenv( 1,
'ZUNMBR',
371 $
'PLN', n, nrhs, m, -1 ) )
372 maxwrk = max( maxwrk, 2*m + m*nrhs )
374 minwrk = max( 2*m + n, 2*m + m*nrhs )
377 minwrk = min( minwrk, maxwrk )
382 IF( lwork.LT.minwrk .AND. .NOT.lquery )
THEN
388 CALL
xerbla(
'ZGELSD', -info )
390 ELSE IF( lquery )
THEN
396 IF( m.EQ.0 .OR. n.EQ.0 )
THEN
404 sfmin = dlamch(
'S' )
406 bignum = one / smlnum
407 CALL
dlabad( smlnum, bignum )
411 anrm = zlange(
'M', m, n, a, lda, rwork )
413 IF( anrm.GT.zero .AND. anrm.LT.smlnum )
THEN
417 CALL
zlascl(
'G', 0, 0, anrm, smlnum, m, n, a, lda, info )
419 ELSE IF( anrm.GT.bignum )
THEN
423 CALL
zlascl(
'G', 0, 0, anrm, bignum, m, n, a, lda, info )
425 ELSE IF( anrm.EQ.zero )
THEN
429 CALL
zlaset(
'F', max( m, n ), nrhs, czero, czero, b, ldb )
430 CALL
dlaset(
'F', minmn, 1, zero, zero, s, 1 )
437 bnrm = zlange(
'M', m, nrhs, b, ldb, rwork )
439 IF( bnrm.GT.zero .AND. bnrm.LT.smlnum )
THEN
443 CALL
zlascl(
'G', 0, 0, bnrm, smlnum, m, nrhs, b, ldb, info )
445 ELSE IF( bnrm.GT.bignum )
THEN
449 CALL
zlascl(
'G', 0, 0, bnrm, bignum, m, nrhs, b, ldb, info )
456 $ CALL
zlaset(
'F', n-m, nrhs, czero, czero, b( m+1, 1 ), ldb )
465 IF( m.GE.mnthr )
THEN
477 CALL
zgeqrf( m, n, a, lda, work( itau ), work( nwork ),
478 $ lwork-nwork+1, info )
484 CALL
zunmqr(
'L',
'C', m, nrhs, n, a, lda, work( itau ), b,
485 $ ldb, work( nwork ), lwork-nwork+1, info )
490 CALL
zlaset(
'L', n-1, n-1, czero, czero, a( 2, 1 ),
505 CALL
zgebrd( mm, n, a, lda, s, rwork( ie ), work( itauq ),
506 $ work( itaup ), work( nwork ), lwork-nwork+1,
512 CALL
zunmbr(
'Q',
'L',
'C', mm, nrhs, n, a, lda, work( itauq ),
513 $ b, ldb, work( nwork ), lwork-nwork+1, info )
517 CALL
zlalsd(
'U', smlsiz, n, nrhs, s, rwork( ie ), b, ldb,
518 $ rcond, rank, work( nwork ), rwork( nrwork ),
526 CALL
zunmbr(
'P',
'L',
'N', n, nrhs, n, a, lda, work( itaup ),
527 $ b, ldb, work( nwork ), lwork-nwork+1, info )
529 ELSE IF( n.GE.mnthr .AND. lwork.GE.4*m+m*m+
530 $ max( m, 2*m-4, nrhs, n-3*m ) )
THEN
536 IF( lwork.GE.max( 4*m+m*lda+max( m, 2*m-4, nrhs, n-3*m ),
537 $ m*lda+m+m*nrhs ) )ldwork = lda
544 CALL
zgelqf( m, n, a, lda, work( itau ), work( nwork ),
545 $ lwork-nwork+1, info )
550 CALL
zlacpy(
'L', m, m, a, lda, work( il ), ldwork )
551 CALL
zlaset(
'U', m-1, m-1, czero, czero, work( il+ldwork ),
553 itauq = il + ldwork*m
563 CALL
zgebrd( m, m, work( il ), ldwork, s, rwork( ie ),
564 $ work( itauq ), work( itaup ), work( nwork ),
565 $ lwork-nwork+1, info )
570 CALL
zunmbr(
'Q',
'L',
'C', m, nrhs, m, work( il ), ldwork,
571 $ work( itauq ), b, ldb, work( nwork ),
572 $ lwork-nwork+1, info )
576 CALL
zlalsd(
'U', smlsiz, m, nrhs, s, rwork( ie ), b, ldb,
577 $ rcond, rank, work( nwork ), rwork( nrwork ),
585 CALL
zunmbr(
'P',
'L',
'N', m, nrhs, m, work( il ), ldwork,
586 $ work( itaup ), b, ldb, work( nwork ),
587 $ lwork-nwork+1, info )
591 CALL
zlaset(
'F', n-m, nrhs, czero, czero, b( m+1, 1 ), ldb )
597 CALL
zunmlq(
'L',
'C', n, nrhs, m, a, lda, work( itau ), b,
598 $ ldb, work( nwork ), lwork-nwork+1, info )
614 CALL
zgebrd( m, n, a, lda, s, rwork( ie ), work( itauq ),
615 $ work( itaup ), work( nwork ), lwork-nwork+1,
621 CALL
zunmbr(
'Q',
'L',
'C', m, nrhs, n, a, lda, work( itauq ),
622 $ b, ldb, work( nwork ), lwork-nwork+1, info )
626 CALL
zlalsd(
'L', smlsiz, m, nrhs, s, rwork( ie ), b, ldb,
627 $ rcond, rank, work( nwork ), rwork( nrwork ),
635 CALL
zunmbr(
'P',
'L',
'N', n, nrhs, m, a, lda, work( itaup ),
636 $ b, ldb, work( nwork ), lwork-nwork+1, info )
642 IF( iascl.EQ.1 )
THEN
643 CALL
zlascl(
'G', 0, 0, anrm, smlnum, n, nrhs, b, ldb, info )
644 CALL
dlascl(
'G', 0, 0, smlnum, anrm, minmn, 1, s, minmn,
646 ELSE IF( iascl.EQ.2 )
THEN
647 CALL
zlascl(
'G', 0, 0, anrm, bignum, n, nrhs, b, ldb, info )
648 CALL
dlascl(
'G', 0, 0, bignum, anrm, minmn, 1, s, minmn,
651 IF( ibscl.EQ.1 )
THEN
652 CALL
zlascl(
'G', 0, 0, smlnum, bnrm, n, nrhs, b, ldb, info )
653 ELSE IF( ibscl.EQ.2 )
THEN
654 CALL
zlascl(
'G', 0, 0, bignum, bnrm, n, nrhs, b, ldb, info )
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 zlalsd(UPLO, SMLSIZ, N, NRHS, D, E, B, LDB, RCOND, RANK, WORK, RWORK, IWORK, INFO)
ZLALSD uses the singular value decomposition of A to solve the least squares problem.
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
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 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 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 zgelsd(M, N, NRHS, A, LDA, B, LDB, S, RCOND, RANK, WORK, LWORK, RWORK, IWORK, INFO)
ZGELSD computes the minimum-norm solution to a linear least squares problem for GE matrices ...
subroutine zunmqr(SIDE, TRANS, M, N, K, A, LDA, TAU, C, LDC, WORK, LWORK, INFO)
ZUNMQR
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