196 INTEGER info, ldb, n, nrhs, rank, smlsiz
201 REAL d( * ), e( * ), rwork( * )
202 COMPLEX b( ldb, * ), work( * )
209 parameter ( zero = 0.0e0, one = 1.0e0, two = 2.0e0 )
211 parameter ( czero = ( 0.0e0, 0.0e0 ) )
214 INTEGER bx, bxst, c, difl, difr, givcol, givnum,
215 $ givptr, i, icmpq1, icmpq2, irwb, irwib, irwrb,
216 $ irwu, irwvt, irwwrk, iwk, j, jcol, jimag,
217 $ jreal, jrow, k, nlvl, nm1, nrwork, nsize, nsub,
218 $ perm, poles, s, sizei, smlszp, sqre, st, st1,
220 REAL cs, eps, orgnrm, r, rcnd, sn, tol
233 INTRINSIC abs, aimag, cmplx, int, log,
REAL, sign
243 ELSE IF( nrhs.LT.1 )
THEN
245 ELSE IF( ( ldb.LT.1 ) .OR. ( ldb.LT.n ) )
THEN
249 CALL xerbla(
'CLALSD', -info )
257 IF( ( rcond.LE.zero ) .OR. ( rcond.GE.one ) )
THEN
269 ELSE IF( n.EQ.1 )
THEN
270 IF( d( 1 ).EQ.zero )
THEN
271 CALL claset(
'A', 1, nrhs, czero, czero, b, ldb )
274 CALL clascl(
'G', 0, 0, d( 1 ), one, 1, nrhs, b, ldb, info )
275 d( 1 ) = abs( d( 1 ) )
282 IF( uplo.EQ.
'L' )
THEN
284 CALL slartg( d( i ), e( i ), cs, sn, r )
287 d( i+1 ) = cs*d( i+1 )
289 CALL csrot( 1, b( i, 1 ), 1, b( i+1, 1 ), 1, cs, sn )
300 CALL csrot( 1, b( j, i ), 1, b( j+1, i ), 1, cs, sn )
309 orgnrm =
slanst(
'M', n, d, e )
310 IF( orgnrm.EQ.zero )
THEN
311 CALL claset(
'A', n, nrhs, czero, czero, b, ldb )
315 CALL slascl(
'G', 0, 0, orgnrm, one, n, 1, d, n, info )
316 CALL slascl(
'G', 0, 0, orgnrm, one, nm1, 1, e, nm1, info )
321 IF( n.LE.smlsiz )
THEN
326 irwib = irwrb + n*nrhs
327 irwb = irwib + n*nrhs
328 CALL slaset(
'A', n, n, zero, one, rwork( irwu ), n )
329 CALL slaset(
'A', n, n, zero, one, rwork( irwvt ), n )
330 CALL slasdq(
'U', 0, n, n, n, 0, d, e, rwork( irwvt ), n,
331 $ rwork( irwu ), n, rwork( irwwrk ), 1,
332 $ rwork( irwwrk ), info )
345 rwork( j ) =
REAL( B( JROW, JCOL ) )
348 CALL sgemm(
'T',
'N', n, nrhs, n, one, rwork( irwu ), n,
349 $ rwork( irwb ), n, zero, rwork( irwrb ), n )
354 rwork( j ) = aimag( b( jrow, jcol ) )
357 CALL sgemm(
'T',
'N', n, nrhs, n, one, rwork( irwu ), n,
358 $ rwork( irwb ), n, zero, rwork( irwib ), n )
365 b( jrow, jcol ) = cmplx( rwork( jreal ), rwork( jimag ) )
369 tol = rcnd*abs( d(
isamax( n, d, 1 ) ) )
371 IF( d( i ).LE.tol )
THEN
372 CALL claset(
'A', 1, nrhs, czero, czero, b( i, 1 ), ldb )
374 CALL clascl(
'G', 0, 0, d( i ), one, 1, nrhs, b( i, 1 ),
388 DO 120 jcol = 1, nrhs
391 rwork( j ) =
REAL( B( JROW, JCOL ) )
394 CALL sgemm(
'T',
'N', n, nrhs, n, one, rwork( irwvt ), n,
395 $ rwork( irwb ), n, zero, rwork( irwrb ), n )
397 DO 140 jcol = 1, nrhs
400 rwork( j ) = aimag( b( jrow, jcol ) )
403 CALL sgemm(
'T',
'N', n, nrhs, n, one, rwork( irwvt ), n,
404 $ rwork( irwb ), n, zero, rwork( irwib ), n )
407 DO 160 jcol = 1, nrhs
411 b( jrow, jcol ) = cmplx( rwork( jreal ), rwork( jimag ) )
417 CALL slascl(
'G', 0, 0, one, orgnrm, n, 1, d, n, info )
418 CALL slasrt(
'D', n, d, info )
419 CALL clascl(
'G', 0, 0, orgnrm, one, n, nrhs, b, ldb, info )
426 nlvl = int( log(
REAL( N ) /
REAL( SMLSIZ+1 ) ) / log( two ) ) + 1
438 givnum = poles + 2*nlvl*n
439 nrwork = givnum + 2*nlvl*n
443 irwib = irwrb + smlsiz*nrhs
444 irwb = irwib + smlsiz*nrhs
450 givcol = perm + nlvl*n
451 iwk = givcol + nlvl*n*2
460 IF( abs( d( i ) ).LT.eps )
THEN
461 d( i ) = sign( eps, d( i ) )
466 IF( ( abs( e( i ) ).LT.eps ) .OR. ( i.EQ.nm1 ) )
THEN
478 iwork( sizei+nsub-1 ) = nsize
479 ELSE IF( abs( e( i ) ).GE.eps )
THEN
484 iwork( sizei+nsub-1 ) = nsize
492 iwork( sizei+nsub-1 ) = nsize
495 iwork( sizei+nsub-1 ) = 1
496 CALL ccopy( nrhs, b( n, 1 ), ldb, work( bx+nm1 ), n )
499 IF( nsize.EQ.1 )
THEN
504 CALL ccopy( nrhs, b( st, 1 ), ldb, work( bx+st1 ), n )
505 ELSE IF( nsize.LE.smlsiz )
THEN
509 CALL slaset(
'A', nsize, nsize, zero, one,
510 $ rwork( vt+st1 ), n )
511 CALL slaset(
'A', nsize, nsize, zero, one,
512 $ rwork( u+st1 ), n )
513 CALL slasdq(
'U', 0, nsize, nsize, nsize, 0, d( st ),
514 $ e( st ), rwork( vt+st1 ), n, rwork( u+st1 ),
515 $ n, rwork( nrwork ), 1, rwork( nrwork ),
526 DO 190 jcol = 1, nrhs
527 DO 180 jrow = st, st + nsize - 1
529 rwork( j ) =
REAL( B( JROW, JCOL ) )
532 CALL sgemm(
'T',
'N', nsize, nrhs, nsize, one,
533 $ rwork( u+st1 ), n, rwork( irwb ), nsize,
534 $ zero, rwork( irwrb ), nsize )
536 DO 210 jcol = 1, nrhs
537 DO 200 jrow = st, st + nsize - 1
539 rwork( j ) = aimag( b( jrow, jcol ) )
542 CALL sgemm(
'T',
'N', nsize, nrhs, nsize, one,
543 $ rwork( u+st1 ), n, rwork( irwb ), nsize,
544 $ zero, rwork( irwib ), nsize )
547 DO 230 jcol = 1, nrhs
548 DO 220 jrow = st, st + nsize - 1
551 b( jrow, jcol ) = cmplx( rwork( jreal ),
556 CALL clacpy(
'A', nsize, nrhs, b( st, 1 ), ldb,
557 $ work( bx+st1 ), n )
562 CALL slasda( icmpq1, smlsiz, nsize, sqre, d( st ),
563 $ e( st ), rwork( u+st1 ), n, rwork( vt+st1 ),
564 $ iwork( k+st1 ), rwork( difl+st1 ),
565 $ rwork( difr+st1 ), rwork( z+st1 ),
566 $ rwork( poles+st1 ), iwork( givptr+st1 ),
567 $ iwork( givcol+st1 ), n, iwork( perm+st1 ),
568 $ rwork( givnum+st1 ), rwork( c+st1 ),
569 $ rwork( s+st1 ), rwork( nrwork ),
570 $ iwork( iwk ), info )
575 CALL clalsa( icmpq2, smlsiz, nsize, nrhs, b( st, 1 ),
576 $ ldb, work( bxst ), n, rwork( u+st1 ), n,
577 $ rwork( vt+st1 ), iwork( k+st1 ),
578 $ rwork( difl+st1 ), rwork( difr+st1 ),
579 $ rwork( z+st1 ), rwork( poles+st1 ),
580 $ iwork( givptr+st1 ), iwork( givcol+st1 ), n,
581 $ iwork( perm+st1 ), rwork( givnum+st1 ),
582 $ rwork( c+st1 ), rwork( s+st1 ),
583 $ rwork( nrwork ), iwork( iwk ), info )
594 tol = rcnd*abs( d(
isamax( n, d, 1 ) ) )
601 IF( abs( d( i ) ).LE.tol )
THEN
602 CALL claset(
'A', 1, nrhs, czero, czero, work( bx+i-1 ), n )
605 CALL clascl(
'G', 0, 0, d( i ), one, 1, nrhs,
606 $ work( bx+i-1 ), n, info )
608 d( i ) = abs( d( i ) )
617 nsize = iwork( sizei+i-1 )
619 IF( nsize.EQ.1 )
THEN
620 CALL ccopy( nrhs, work( bxst ), n, b( st, 1 ), ldb )
621 ELSE IF( nsize.LE.smlsiz )
THEN
632 DO 270 jcol = 1, nrhs
634 DO 260 jrow = 1, nsize
636 rwork( jreal ) =
REAL( WORK( J+JROW ) )
639 CALL sgemm(
'T',
'N', nsize, nrhs, nsize, one,
640 $ rwork( vt+st1 ), n, rwork( irwb ), nsize, zero,
641 $ rwork( irwrb ), nsize )
644 DO 290 jcol = 1, nrhs
646 DO 280 jrow = 1, nsize
648 rwork( jimag ) = aimag( work( j+jrow ) )
651 CALL sgemm(
'T',
'N', nsize, nrhs, nsize, one,
652 $ rwork( vt+st1 ), n, rwork( irwb ), nsize, zero,
653 $ rwork( irwib ), nsize )
656 DO 310 jcol = 1, nrhs
657 DO 300 jrow = st, st + nsize - 1
660 b( jrow, jcol ) = cmplx( rwork( jreal ),
665 CALL clalsa( icmpq2, smlsiz, nsize, nrhs, work( bxst ), n,
666 $ b( st, 1 ), ldb, rwork( u+st1 ), n,
667 $ rwork( vt+st1 ), iwork( k+st1 ),
668 $ rwork( difl+st1 ), rwork( difr+st1 ),
669 $ rwork( z+st1 ), rwork( poles+st1 ),
670 $ iwork( givptr+st1 ), iwork( givcol+st1 ), n,
671 $ iwork( perm+st1 ), rwork( givnum+st1 ),
672 $ rwork( c+st1 ), rwork( s+st1 ),
673 $ rwork( nrwork ), iwork( iwk ), info )
682 CALL slascl(
'G', 0, 0, one, orgnrm, n, 1, d, n, info )
683 CALL slasrt(
'D', n, d, info )
684 CALL clascl(
'G', 0, 0, orgnrm, one, n, nrhs, b, ldb, info )
subroutine slasrt(ID, N, D, INFO)
SLASRT sorts numbers in increasing or decreasing order.
subroutine slartg(F, G, CS, SN, R)
SLARTG generates a plane rotation with real cosine and real sine.
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 sgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
SGEMM
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.
integer function isamax(N, SX, INCX)
ISAMAX
subroutine slasdq(UPLO, SQRE, N, NCVT, NRU, NCC, D, E, VT, LDVT, U, LDU, C, LDC, WORK, INFO)
SLASDQ computes the SVD of a real bidiagonal matrix with diagonal d and off-diagonal e...
subroutine slasda(ICOMPQ, SMLSIZ, N, SQRE, D, E, U, LDU, VT, K, DIFL, DIFR, Z, POLES, GIVPTR, GIVCOL, LDGCOL, PERM, GIVNUM, C, S, WORK, IWORK, INFO)
SLASDA computes the singular value decomposition (SVD) of a real upper bidiagonal matrix with diagona...
subroutine xerbla(SRNAME, INFO)
XERBLA
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 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.
real function slamch(CMACH)
SLAMCH
subroutine clacpy(UPLO, M, N, A, LDA, B, LDB)
CLACPY copies all or part of one two-dimensional array to another.
subroutine csrot(N, CX, INCX, CY, INCY, C, S)
CSROT
subroutine ccopy(N, CX, INCX, CY, INCY)
CCOPY
subroutine clalsa(ICOMPQ, SMLSIZ, N, NRHS, B, LDB, BX, LDBX, U, LDU, VT, K, DIFL, DIFR, Z, POLES, GIVPTR, GIVCOL, LDGCOL, PERM, GIVNUM, C, S, RWORK, IWORK, INFO)
CLALSA computes the SVD of the coefficient matrix in compact form. Used by sgelsd.
real function slanst(NORM, N, D, E)
SLANST returns the value of the 1-norm, or the Frobenius norm, or the infinity norm, or the element of largest absolute value of a real symmetric tridiagonal matrix.