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
267 EXTERNAL clange, slamch, ilaenv
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
405 sfmin = slamch(
'S' )
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 )
subroutine clalsd(UPLO, SMLSIZ, N, NRHS, D, E, B, LDB, RCOND, RANK, WORK, RWORK, IWORK, INFO)
CLALSD uses the singular value decomposition of A to solve the least squares problem.
subroutine claset(UPLO, M, N, ALPHA, BETA, A, LDA)
CLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values...
subroutine clascl(TYPE, KL, KU, CFROM, CTO, M, N, A, LDA, INFO)
CLASCL multiplies a general rectangular matrix by a real scalar defined as cto/cfrom.
subroutine cgelsd(M, N, NRHS, A, LDA, B, LDB, S, RCOND, RANK, WORK, LWORK, RWORK, IWORK, INFO)
CGELSD computes the minimum-norm solution to a linear least squares problem for GE matrices ...
subroutine cgelqf(M, N, A, LDA, TAU, WORK, LWORK, INFO)
CGELQF
subroutine xerbla(SRNAME, INFO)
XERBLA
subroutine cgeqrf(M, N, A, LDA, TAU, WORK, LWORK, INFO)
CGEQRF
subroutine slaset(UPLO, M, N, ALPHA, BETA, A, LDA)
SLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values...
subroutine cgebrd(M, N, A, LDA, D, E, TAUQ, TAUP, WORK, LWORK, INFO)
CGEBRD
subroutine slascl(TYPE, KL, KU, CFROM, CTO, M, N, A, LDA, INFO)
SLASCL multiplies a general rectangular matrix by a real scalar defined as cto/cfrom.
subroutine clacpy(UPLO, M, N, A, LDA, B, LDB)
CLACPY copies all or part of one two-dimensional array to another.
subroutine slabad(SMALL, LARGE)
SLABAD
subroutine cunmqr(SIDE, TRANS, M, N, K, A, LDA, TAU, C, LDC, WORK, LWORK, INFO)
CUNMQR
subroutine cunmbr(VECT, SIDE, TRANS, M, N, K, A, LDA, TAU, C, LDC, WORK, LWORK, INFO)
CUNMBR
subroutine cunmlq(SIDE, TRANS, M, N, K, A, LDA, TAU, C, LDC, WORK, LWORK, INFO)
CUNMLQ