185 SUBROUTINE zlalsd( UPLO, SMLSIZ, N, NRHS, D, E, B, LDB, RCOND,
186 $ RANK, WORK, RWORK, IWORK, INFO )
194 INTEGER INFO, LDB, N, NRHS, RANK, SMLSIZ
195 DOUBLE PRECISION RCOND
199 DOUBLE PRECISION D( * ), E( * ), RWORK( * )
200 COMPLEX*16 B( LDB, * ), WORK( * )
206 DOUBLE PRECISION ZERO, ONE, TWO
207 parameter( zero = 0.0d0, one = 1.0d0, two = 2.0d0 )
209 parameter( czero = ( 0.0d0, 0.0d0 ) )
212 INTEGER BX, BXST, C, DIFL, DIFR, GIVCOL, GIVNUM,
213 $ givptr, i, icmpq1, icmpq2, irwb, irwib, irwrb,
214 $ irwu, irwvt, irwwrk, iwk, j, jcol, jimag,
215 $ jreal, jrow, k, nlvl, nm1, nrwork, nsize, nsub,
216 $ perm, poles, s, sizei, smlszp, sqre, st, st1,
218 DOUBLE PRECISION CS, EPS, ORGNRM, RCND, R, SN, TOL
222 DOUBLE PRECISION DLAMCH, DLANST
223 EXTERNAL idamax, dlamch, dlanst
231 INTRINSIC abs, dble, dcmplx, dimag, int, log, sign
241 ELSE IF( nrhs.LT.1 )
THEN
243 ELSE IF( ( ldb.LT.1 ) .OR. ( ldb.LT.n ) )
THEN
247 CALL xerbla(
'ZLALSD', -info )
251 eps = dlamch(
'Epsilon' )
255 IF( ( rcond.LE.zero ) .OR. ( rcond.GE.one ) )
THEN
267 ELSE IF( n.EQ.1 )
THEN
268 IF( d( 1 ).EQ.zero )
THEN
269 CALL zlaset(
'A', 1, nrhs, czero, czero, b, ldb )
272 CALL zlascl(
'G', 0, 0, d( 1 ), one, 1, nrhs, b, ldb, info )
273 d( 1 ) = abs( d( 1 ) )
280 IF( uplo.EQ.
'L' )
THEN
282 CALL dlartg( d( i ), e( i ), cs, sn, r )
285 d( i+1 ) = cs*d( i+1 )
287 CALL zdrot( 1, b( i, 1 ), 1, b( i+1, 1 ), 1, cs, sn )
298 CALL zdrot( 1, b( j, i ), 1, b( j+1, i ), 1, cs, sn )
307 orgnrm = dlanst(
'M', n, d, e )
308 IF( orgnrm.EQ.zero )
THEN
309 CALL zlaset(
'A', n, nrhs, czero, czero, b, ldb )
313 CALL dlascl(
'G', 0, 0, orgnrm, one, n, 1, d, n, info )
314 CALL dlascl(
'G', 0, 0, orgnrm, one, nm1, 1, e, nm1, info )
319 IF( n.LE.smlsiz )
THEN
324 irwib = irwrb + n*nrhs
325 irwb = irwib + n*nrhs
326 CALL dlaset(
'A', n, n, zero, one, rwork( irwu ), n )
327 CALL dlaset(
'A', n, n, zero, one, rwork( irwvt ), n )
328 CALL dlasdq(
'U', 0, n, n, n, 0, d, e, rwork( irwvt ), n,
329 $ rwork( irwu ), n, rwork( irwwrk ), 1,
330 $ rwork( irwwrk ), info )
343 rwork( j ) = dble( b( jrow, jcol ) )
346 CALL dgemm(
'T',
'N', n, nrhs, n, one, rwork( irwu ), n,
347 $ rwork( irwb ), n, zero, rwork( irwrb ), n )
352 rwork( j ) = dimag( b( jrow, jcol ) )
355 CALL dgemm(
'T',
'N', n, nrhs, n, one, rwork( irwu ), n,
356 $ rwork( irwb ), n, zero, rwork( irwib ), n )
363 b( jrow, jcol ) = dcmplx( rwork( jreal ),
368 tol = rcnd*abs( d( idamax( n, d, 1 ) ) )
370 IF( d( i ).LE.tol )
THEN
371 CALL zlaset(
'A', 1, nrhs, czero, czero, b( i, 1 ), ldb )
373 CALL zlascl(
'G', 0, 0, d( i ), one, 1, nrhs, b( i, 1 ),
387 DO 120 jcol = 1, nrhs
390 rwork( j ) = dble( b( jrow, jcol ) )
393 CALL dgemm(
'T',
'N', n, nrhs, n, one, rwork( irwvt ), n,
394 $ rwork( irwb ), n, zero, rwork( irwrb ), n )
396 DO 140 jcol = 1, nrhs
399 rwork( j ) = dimag( b( jrow, jcol ) )
402 CALL dgemm(
'T',
'N', n, nrhs, n, one, rwork( irwvt ), n,
403 $ rwork( irwb ), n, zero, rwork( irwib ), n )
406 DO 160 jcol = 1, nrhs
410 b( jrow, jcol ) = dcmplx( rwork( jreal ),
417 CALL dlascl(
'G', 0, 0, one, orgnrm, n, 1, d, n, info )
418 CALL dlasrt(
'D', n, d, info )
419 CALL zlascl(
'G', 0, 0, orgnrm, one, n, nrhs, b, ldb, info )
426 nlvl = int( log( dble( n ) / dble( 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 zcopy( nrhs, b( n, 1 ), ldb, work( bx+nm1 ), n )
499 IF( nsize.EQ.1 )
THEN
504 CALL zcopy( nrhs, b( st, 1 ), ldb, work( bx+st1 ), n )
505 ELSE IF( nsize.LE.smlsiz )
THEN
509 CALL dlaset(
'A', nsize, nsize, zero, one,
510 $ rwork( vt+st1 ), n )
511 CALL dlaset(
'A', nsize, nsize, zero, one,
512 $ rwork( u+st1 ), n )
513 CALL dlasdq(
'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 ) = dble( b( jrow, jcol ) )
532 CALL dgemm(
'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 ) = dimag( b( jrow, jcol ) )
542 CALL dgemm(
'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 ) = dcmplx( rwork( jreal ),
556 CALL zlacpy(
'A', nsize, nrhs, b( st, 1 ), ldb,
557 $ work( bx+st1 ), n )
562 CALL dlasda( 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 zlalsa( 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( idamax( n, d, 1 ) ) )
601 IF( abs( d( i ) ).LE.tol )
THEN
602 CALL zlaset(
'A', 1, nrhs, czero, czero, work( bx+i-1 ), n )
605 CALL zlascl(
'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 zcopy( 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 ) = dble( work( j+jrow ) )
639 CALL dgemm(
'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 ) = dimag( work( j+jrow ) )
651 CALL dgemm(
'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 ) = dcmplx( rwork( jreal ),
665 CALL zlalsa( 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 dlascl(
'G', 0, 0, one, orgnrm, n, 1, d, n, info )
683 CALL dlasrt(
'D', n, d, info )
684 CALL zlascl(
'G', 0, 0, orgnrm, one, n, nrhs, b, ldb, info )
subroutine dlascl(TYPE, KL, KU, CFROM, CTO, M, N, A, LDA, INFO)
DLASCL multiplies a general rectangular matrix by a real scalar defined as cto/cfrom.
subroutine dlartg(f, g, c, s, r)
DLARTG generates a plane rotation with real cosine and real sine.
subroutine dlaset(UPLO, M, N, ALPHA, BETA, A, LDA)
DLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values.
subroutine dlasda(ICOMPQ, SMLSIZ, N, SQRE, D, E, U, LDU, VT, K, DIFL, DIFR, Z, POLES, GIVPTR, GIVCOL, LDGCOL, PERM, GIVNUM, C, S, WORK, IWORK, INFO)
DLASDA computes the singular value decomposition (SVD) of a real upper bidiagonal matrix with diagona...
subroutine dlasdq(UPLO, SQRE, N, NCVT, NRU, NCC, D, E, VT, LDVT, U, LDU, C, LDC, WORK, INFO)
DLASDQ computes the SVD of a real bidiagonal matrix with diagonal d and off-diagonal e....
subroutine xerbla(SRNAME, INFO)
XERBLA
subroutine dlasrt(ID, N, D, INFO)
DLASRT sorts numbers in increasing or decreasing order.
subroutine zdrot(N, ZX, INCX, ZY, INCY, C, S)
ZDROT
subroutine zcopy(N, ZX, INCX, ZY, INCY)
ZCOPY
subroutine zlascl(TYPE, KL, KU, CFROM, CTO, M, N, A, LDA, INFO)
ZLASCL multiplies a general rectangular matrix by a real scalar defined as cto/cfrom.
subroutine zlacpy(UPLO, M, N, A, LDA, B, LDB)
ZLACPY copies all or part of one two-dimensional array to another.
subroutine zlaset(UPLO, M, N, ALPHA, BETA, A, LDA)
ZLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values.
subroutine zlalsd(UPLO, SMLSIZ, N, NRHS, D, E, B, LDB, RCOND, RANK, WORK, RWORK, IWORK, INFO)
ZLALSD uses the singular value decomposition of A to solve the least squares problem.
subroutine zlalsa(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)
ZLALSA computes the SVD of the coefficient matrix in compact form. Used by sgelsd.
subroutine dgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
DGEMM