225 SUBROUTINE cgelsd( 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
239 REAL rwork( * ), s( * )
240 COMPLEX a( lda, * ), b( ldb, * ), work( * )
247 parameter( zero = 0.0e+0, one = 1.0e+0, two = 2.0e+0 )
249 parameter( czero = ( 0.0e+0, 0.0e+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 REAL anrm, bignum, bnrm, eps, sfmin, smlnum
270 INTRINSIC int, log, max, min, real
279 lquery = ( lwork.EQ.-1 )
282 ELSE IF( n.LT.0 )
THEN
284 ELSE IF( nrhs.LT.0 )
THEN
286 ELSE IF( lda.LT.max( 1, m ) )
THEN
288 ELSE IF( ldb.LT.max( 1, maxmn ) )
THEN
304 IF( minmn.GT.0 )
THEN
305 smlsiz =
ilaenv( 9,
'CGELSD',
' ', 0, 0, 0, 0 )
306 mnthr =
ilaenv( 6,
'CGELSD',
' ', m, n, nrhs, -1 )
307 nlvl = max( int( log(
REAL( MINMN ) /
REAL( SMLSIZ + 1 ) ) /
308 $ log( two ) ) + 1, 0 )
309 liwork = 3*minmn*nlvl + 11*minmn
311 IF( m.GE.n .AND. m.GE.mnthr )
THEN
317 maxwrk = max( maxwrk, n*
ilaenv( 1,
'CGEQRF',
' ', m, n,
319 maxwrk = max( maxwrk, nrhs*
ilaenv( 1,
'CUNMQR',
'LC', m,
326 lrwork = 10*n + 2*n*smlsiz + 8*n*nlvl + 3*smlsiz*nrhs +
327 $ max( (smlsiz+1)**2, n*(1+nrhs) + 2*nrhs )
328 maxwrk = max( maxwrk, 2*n + ( mm + n )*
ilaenv( 1,
329 $
'CGEBRD',
' ', mm, n, -1, -1 ) )
330 maxwrk = max( maxwrk, 2*n + nrhs*
ilaenv( 1,
'CUNMBR',
331 $
'QLC', mm, nrhs, n, -1 ) )
332 maxwrk = max( maxwrk, 2*n + ( n - 1 )*
ilaenv( 1,
333 $
'CUNMBR',
'PLN', n, nrhs, n, -1 ) )
334 maxwrk = max( maxwrk, 2*n + n*nrhs )
335 minwrk = max( 2*n + mm, 2*n + n*nrhs )
338 lrwork = 10*m + 2*m*smlsiz + 8*m*nlvl + 3*smlsiz*nrhs +
339 $ max( (smlsiz+1)**2, n*(1+nrhs) + 2*nrhs )
340 IF( n.GE.mnthr )
THEN
345 maxwrk = m + m*
ilaenv( 1,
'CGELQF',
' ', m, n, -1,
347 maxwrk = max( maxwrk, m*m + 4*m + 2*m*
ilaenv( 1,
348 $
'CGEBRD',
' ', m, m, -1, -1 ) )
349 maxwrk = max( maxwrk, m*m + 4*m + nrhs*
ilaenv( 1,
350 $
'CUNMBR',
'QLC', m, nrhs, m, -1 ) )
351 maxwrk = max( maxwrk, m*m + 4*m + ( m - 1 )*
ilaenv( 1,
352 $
'CUNMLQ',
'LC', n, nrhs, m, -1 ) )
354 maxwrk = max( maxwrk, m*m + m + m*nrhs )
356 maxwrk = max( maxwrk, m*m + 2*m )
358 maxwrk = max( maxwrk, m*m + 4*m + m*nrhs )
361 maxwrk = max( maxwrk,
362 $ 4*m+m*m+max( m, 2*m-4, nrhs, n-3*m ) )
367 maxwrk = 2*m + ( n + m )*
ilaenv( 1,
'CGEBRD',
' ', m,
369 maxwrk = max( maxwrk, 2*m + nrhs*
ilaenv( 1,
'CUNMBR',
370 $
'QLC', m, nrhs, m, -1 ) )
371 maxwrk = max( maxwrk, 2*m + m*
ilaenv( 1,
'CUNMBR',
372 $
'PLN', n, nrhs, m, -1 ) )
373 maxwrk = max( maxwrk, 2*m + m*nrhs )
375 minwrk = max( 2*m + n, 2*m + m*nrhs )
378 minwrk = min( minwrk, maxwrk )
383 IF( lwork.LT.minwrk .AND. .NOT.lquery )
THEN
389 CALL
xerbla(
'CGELSD', -info )
391 ELSE IF( lquery )
THEN
397 IF( m.EQ.0 .OR. n.EQ.0 )
THEN
407 bignum = one / smlnum
408 CALL
slabad( smlnum, bignum )
412 anrm =
clange(
'M', m, n, a, lda, rwork )
414 IF( anrm.GT.zero .AND. anrm.LT.smlnum )
THEN
418 CALL
clascl(
'G', 0, 0, anrm, smlnum, m, n, a, lda, info )
420 ELSE IF( anrm.GT.bignum )
THEN
424 CALL
clascl(
'G', 0, 0, anrm, bignum, m, n, a, lda, info )
426 ELSE IF( anrm.EQ.zero )
THEN
430 CALL
claset(
'F', max( m, n ), nrhs, czero, czero, b, ldb )
431 CALL
slaset(
'F', minmn, 1, zero, zero, s, 1 )
438 bnrm =
clange(
'M', m, nrhs, b, ldb, rwork )
440 IF( bnrm.GT.zero .AND. bnrm.LT.smlnum )
THEN
444 CALL
clascl(
'G', 0, 0, bnrm, smlnum, m, nrhs, b, ldb, info )
446 ELSE IF( bnrm.GT.bignum )
THEN
450 CALL
clascl(
'G', 0, 0, bnrm, bignum, m, nrhs, b, ldb, info )
457 $ CALL
claset(
'F', n-m, nrhs, czero, czero, b( m+1, 1 ), ldb )
466 IF( m.GE.mnthr )
THEN
478 CALL
cgeqrf( m, n, a, lda, work( itau ), work( nwork ),
479 $ lwork-nwork+1, info )
485 CALL
cunmqr(
'L',
'C', m, nrhs, n, a, lda, work( itau ), b,
486 $ ldb, work( nwork ), lwork-nwork+1, info )
491 CALL
claset(
'L', n-1, n-1, czero, czero, a( 2, 1 ),
506 CALL
cgebrd( mm, n, a, lda, s, rwork( ie ), work( itauq ),
507 $ work( itaup ), work( nwork ), lwork-nwork+1,
513 CALL
cunmbr(
'Q',
'L',
'C', mm, nrhs, n, a, lda, work( itauq ),
514 $ b, ldb, work( nwork ), lwork-nwork+1, info )
518 CALL
clalsd(
'U', smlsiz, n, nrhs, s, rwork( ie ), b, ldb,
519 $ rcond, rank, work( nwork ), rwork( nrwork ),
527 CALL
cunmbr(
'P',
'L',
'N', n, nrhs, n, a, lda, work( itaup ),
528 $ b, ldb, work( nwork ), lwork-nwork+1, info )
530 ELSE IF( n.GE.mnthr .AND. lwork.GE.4*m+m*m+
531 $ max( m, 2*m-4, nrhs, n-3*m ) )
THEN
537 IF( lwork.GE.max( 4*m+m*lda+max( m, 2*m-4, nrhs, n-3*m ),
538 $ m*lda+m+m*nrhs ) )ldwork = lda
545 CALL
cgelqf( m, n, a, lda, work( itau ), work( nwork ),
546 $ lwork-nwork+1, info )
551 CALL
clacpy(
'L', m, m, a, lda, work( il ), ldwork )
552 CALL
claset(
'U', m-1, m-1, czero, czero, work( il+ldwork ),
554 itauq = il + ldwork*m
564 CALL
cgebrd( m, m, work( il ), ldwork, s, rwork( ie ),
565 $ work( itauq ), work( itaup ), work( nwork ),
566 $ lwork-nwork+1, info )
571 CALL
cunmbr(
'Q',
'L',
'C', m, nrhs, m, work( il ), ldwork,
572 $ work( itauq ), b, ldb, work( nwork ),
573 $ lwork-nwork+1, info )
577 CALL
clalsd(
'U', smlsiz, m, nrhs, s, rwork( ie ), b, ldb,
578 $ rcond, rank, work( nwork ), rwork( nrwork ),
586 CALL
cunmbr(
'P',
'L',
'N', m, nrhs, m, work( il ), ldwork,
587 $ work( itaup ), b, ldb, work( nwork ),
588 $ lwork-nwork+1, info )
592 CALL
claset(
'F', n-m, nrhs, czero, czero, b( m+1, 1 ), ldb )
598 CALL
cunmlq(
'L',
'C', n, nrhs, m, a, lda, work( itau ), b,
599 $ ldb, work( nwork ), lwork-nwork+1, info )
615 CALL
cgebrd( m, n, a, lda, s, rwork( ie ), work( itauq ),
616 $ work( itaup ), work( nwork ), lwork-nwork+1,
622 CALL
cunmbr(
'Q',
'L',
'C', m, nrhs, n, a, lda, work( itauq ),
623 $ b, ldb, work( nwork ), lwork-nwork+1, info )
627 CALL
clalsd(
'L', smlsiz, m, nrhs, s, rwork( ie ), b, ldb,
628 $ rcond, rank, work( nwork ), rwork( nrwork ),
636 CALL
cunmbr(
'P',
'L',
'N', n, nrhs, m, a, lda, work( itaup ),
637 $ b, ldb, work( nwork ), lwork-nwork+1, info )
643 IF( iascl.EQ.1 )
THEN
644 CALL
clascl(
'G', 0, 0, anrm, smlnum, n, nrhs, b, ldb, info )
645 CALL
slascl(
'G', 0, 0, smlnum, anrm, minmn, 1, s, minmn,
647 ELSE IF( iascl.EQ.2 )
THEN
648 CALL
clascl(
'G', 0, 0, anrm, bignum, n, nrhs, b, ldb, info )
649 CALL
slascl(
'G', 0, 0, bignum, anrm, minmn, 1, s, minmn,
652 IF( ibscl.EQ.1 )
THEN
653 CALL
clascl(
'G', 0, 0, smlnum, bnrm, n, nrhs, b, ldb, info )
654 ELSE IF( ibscl.EQ.2 )
THEN
655 CALL
clascl(
'G', 0, 0, bignum, bnrm, n, nrhs, b, ldb, info )