215 SUBROUTINE zgelsd( M, N, NRHS, A, LDA, B, LDB, S, RCOND, RANK,
216 $ WORK, LWORK, RWORK, IWORK, INFO )
223 INTEGER INFO, LDA, LDB, LWORK, M, N, NRHS, RANK
224 DOUBLE PRECISION RCOND
228 DOUBLE PRECISION RWORK( * ), S( * )
229 COMPLEX*16 A( LDA, * ), B( LDB, * ), WORK( * )
235 DOUBLE PRECISION ZERO, ONE, TWO
236 parameter( zero = 0.0d+0, one = 1.0d+0, two = 2.0d+0 )
238 parameter( czero = ( 0.0d+0, 0.0d+0 ) )
242 INTEGER IASCL, IBSCL, IE, IL, ITAU, ITAUP, ITAUQ,
243 $ ldwork, liwork, lrwork, maxmn, maxwrk, minmn,
244 $ minwrk, mm, mnthr, nlvl, nrwork, nwork, smlsiz
245 DOUBLE PRECISION ANRM, BIGNUM, BNRM, EPS, SFMIN, SMLNUM
255 DOUBLE PRECISION DLAMCH, ZLANGE
256 EXTERNAL ilaenv, dlamch, zlange
259 INTRINSIC int, log, max, min, dble
268 lquery = ( lwork.EQ.-1 )
271 ELSE IF( n.LT.0 )
THEN
273 ELSE IF( nrhs.LT.0 )
THEN
275 ELSE IF( lda.LT.max( 1, m ) )
THEN
277 ELSE IF( ldb.LT.max( 1, maxmn ) )
THEN
293 IF( minmn.GT.0 )
THEN
294 smlsiz = ilaenv( 9,
'ZGELSD',
' ', 0, 0, 0, 0 )
295 mnthr = ilaenv( 6,
'ZGELSD',
' ', m, n, nrhs, -1 )
296 nlvl = max( int( log( dble( minmn ) / dble( smlsiz + 1 ) ) /
297 $ log( two ) ) + 1, 0 )
298 liwork = 3*minmn*nlvl + 11*minmn
300 IF( m.GE.n .AND. m.GE.mnthr )
THEN
306 maxwrk = max( maxwrk, n*ilaenv( 1,
'ZGEQRF',
' ', m,
309 maxwrk = max( maxwrk, nrhs*ilaenv( 1,
'ZUNMQR',
'LC',
317 lrwork = 10*n + 2*n*smlsiz + 8*n*nlvl + 3*smlsiz*nrhs +
318 $ max( (smlsiz+1)**2, n*(1+nrhs) + 2*nrhs )
319 maxwrk = max( maxwrk, 2*n + ( mm + n )*ilaenv( 1,
320 $
'ZGEBRD',
' ', mm, n, -1, -1 ) )
321 maxwrk = max( maxwrk, 2*n + nrhs*ilaenv( 1,
'ZUNMBR',
322 $
'QLC', mm, nrhs, n, -1 ) )
323 maxwrk = max( maxwrk, 2*n + ( n - 1 )*ilaenv( 1,
324 $
'ZUNMBR',
'PLN', n, nrhs, n, -1 ) )
325 maxwrk = max( maxwrk, 2*n + n*nrhs )
326 minwrk = max( 2*n + mm, 2*n + n*nrhs )
329 lrwork = 10*m + 2*m*smlsiz + 8*m*nlvl + 3*smlsiz*nrhs +
330 $ max( (smlsiz+1)**2, n*(1+nrhs) + 2*nrhs )
331 IF( n.GE.mnthr )
THEN
336 maxwrk = m + m*ilaenv( 1,
'ZGELQF',
' ', m, n, -1,
338 maxwrk = max( maxwrk, m*m + 4*m + 2*m*ilaenv( 1,
339 $
'ZGEBRD',
' ', m, m, -1, -1 ) )
340 maxwrk = max( maxwrk, m*m + 4*m + nrhs*ilaenv( 1,
341 $
'ZUNMBR',
'QLC', m, nrhs, m, -1 ) )
342 maxwrk = max( maxwrk,
343 $ m*m + 4*m + ( m - 1 )*ilaenv( 1,
344 $
'ZUNMLQ',
'LC', n, nrhs, m, -1 ) )
346 maxwrk = max( maxwrk, m*m + m + m*nrhs )
348 maxwrk = max( maxwrk, m*m + 2*m )
350 maxwrk = max( maxwrk, m*m + 4*m + m*nrhs )
353 maxwrk = max( maxwrk,
354 $ 4*m+m*m+max( m, 2*m-4, nrhs, n-3*m ) )
359 maxwrk = 2*m + ( n + m )*ilaenv( 1,
'ZGEBRD',
' ',
362 maxwrk = max( maxwrk, 2*m + nrhs*ilaenv( 1,
364 $
'QLC', m, nrhs, m, -1 ) )
365 maxwrk = max( maxwrk, 2*m + m*ilaenv( 1,
'ZUNMBR',
366 $
'PLN', n, nrhs, m, -1 ) )
367 maxwrk = max( maxwrk, 2*m + m*nrhs )
369 minwrk = max( 2*m + n, 2*m + m*nrhs )
372 minwrk = min( minwrk, maxwrk )
377 IF( lwork.LT.minwrk .AND. .NOT.lquery )
THEN
383 CALL xerbla(
'ZGELSD', -info )
385 ELSE IF( lquery )
THEN
391 IF( m.EQ.0 .OR. n.EQ.0 )
THEN
399 sfmin = dlamch(
'S' )
401 bignum = one / smlnum
405 anrm = zlange(
'M', m, n, a, lda, rwork )
407 IF( anrm.GT.zero .AND. anrm.LT.smlnum )
THEN
411 CALL zlascl(
'G', 0, 0, anrm, smlnum, m, n, a, lda, info )
413 ELSE IF( anrm.GT.bignum )
THEN
417 CALL zlascl(
'G', 0, 0, anrm, bignum, m, n, a, lda, info )
419 ELSE IF( anrm.EQ.zero )
THEN
423 CALL zlaset(
'F', max( m, n ), nrhs, czero, czero, b, ldb )
424 CALL dlaset(
'F', minmn, 1, zero, zero, s, 1 )
431 bnrm = zlange(
'M', m, nrhs, b, ldb, rwork )
433 IF( bnrm.GT.zero .AND. bnrm.LT.smlnum )
THEN
437 CALL zlascl(
'G', 0, 0, bnrm, smlnum, m, nrhs, b, ldb,
440 ELSE IF( bnrm.GT.bignum )
THEN
444 CALL zlascl(
'G', 0, 0, bnrm, bignum, m, nrhs, b, ldb,
452 $
CALL zlaset(
'F', n-m, nrhs, czero, czero, b( m+1, 1 ),
462 IF( m.GE.mnthr )
THEN
474 CALL zgeqrf( m, n, a, lda, work( itau ), work( nwork ),
475 $ lwork-nwork+1, info )
481 CALL zunmqr(
'L',
'C', m, nrhs, n, a, lda, work( itau ),
483 $ ldb, work( nwork ), lwork-nwork+1, info )
488 CALL zlaset(
'L', n-1, n-1, czero, czero, a( 2, 1 ),
503 CALL zgebrd( mm, n, a, lda, s, rwork( ie ), work( itauq ),
504 $ work( itaup ), work( nwork ), lwork-nwork+1,
510 CALL zunmbr(
'Q',
'L',
'C', mm, nrhs, n, a, lda,
512 $ b, ldb, work( nwork ), lwork-nwork+1, info )
516 CALL zlalsd(
'U', smlsiz, n, nrhs, s, rwork( ie ), b, ldb,
517 $ rcond, rank, work( nwork ), rwork( nrwork ),
525 CALL zunmbr(
'P',
'L',
'N', n, nrhs, n, a, lda,
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 ),
598 CALL zunmlq(
'L',
'C', n, nrhs, m, a, lda, work( itau ), b,
599 $ ldb, work( nwork ), lwork-nwork+1, info )
615 CALL zgebrd( m, n, a, lda, s, rwork( ie ), work( itauq ),
616 $ work( itaup ), work( nwork ), lwork-nwork+1,
622 CALL zunmbr(
'Q',
'L',
'C', m, nrhs, n, a, lda,
624 $ b, ldb, work( nwork ), lwork-nwork+1, info )
628 CALL zlalsd(
'L', smlsiz, m, nrhs, s, rwork( ie ), b, ldb,
629 $ rcond, rank, work( nwork ), rwork( nrwork ),
637 CALL zunmbr(
'P',
'L',
'N', n, nrhs, m, a, lda,
639 $ b, ldb, work( nwork ), lwork-nwork+1, info )
645 IF( iascl.EQ.1 )
THEN
646 CALL zlascl(
'G', 0, 0, anrm, smlnum, n, nrhs, b, ldb,
648 CALL dlascl(
'G', 0, 0, smlnum, anrm, minmn, 1, s, minmn,
650 ELSE IF( iascl.EQ.2 )
THEN
651 CALL zlascl(
'G', 0, 0, anrm, bignum, n, nrhs, b, ldb,
653 CALL dlascl(
'G', 0, 0, bignum, anrm, minmn, 1, s, minmn,
656 IF( ibscl.EQ.1 )
THEN
657 CALL zlascl(
'G', 0, 0, smlnum, bnrm, n, nrhs, b, ldb,
659 ELSE IF( ibscl.EQ.2 )
THEN
660 CALL zlascl(
'G', 0, 0, bignum, bnrm, n, nrhs, b, ldb,