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 )
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 slaset(UPLO, M, N, ALPHA, BETA, A, LDA)
SLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values...
real function clange(NORM, M, N, A, LDA, WORK)
CLANGE returns the value of the 1-norm, Frobenius norm, infinity-norm, or the largest absolute value ...
subroutine cunmbr(VECT, SIDE, TRANS, M, N, K, A, LDA, TAU, C, LDC, WORK, LWORK, INFO)
CUNMBR
subroutine cgebrd(M, N, A, LDA, D, E, TAUQ, TAUP, WORK, LWORK, INFO)
CGEBRD
subroutine xerbla(SRNAME, INFO)
XERBLA
subroutine cunmqr(SIDE, TRANS, M, N, K, A, LDA, TAU, C, LDC, WORK, LWORK, INFO)
CUNMQR
set ue cd $ADTTMP cat<< EOF > tmp f Program LinearEquations Implicit none Real b(3) integer i
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 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 clacpy(UPLO, M, N, A, LDA, B, LDB)
CLACPY copies all or part of one two-dimensional array to another.
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
real function slamch(CMACH)
SLAMCH
subroutine cunmlq(SIDE, TRANS, M, N, K, A, LDA, TAU, C, LDC, WORK, LWORK, INFO)
CUNMLQ
subroutine cgeqrf(M, N, A, LDA, TAU, WORK, LWORK, INFO)
CGEQRF
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 slabad(SMALL, LARGE)
SLABAD
integer function ilaenv(ISPEC, NAME, OPTS, N1, N2, N3, N4)