198 INTEGER info, ldb, n, nrhs, rank, smlsiz
199 DOUBLE PRECISION rcond
203 DOUBLE PRECISION d( * ), e( * ), rwork( * )
204 COMPLEX*16 b( ldb, * ), work( * )
210 DOUBLE PRECISION zero, one, two
211 parameter ( zero = 0.0d0, one = 1.0d0, two = 2.0d0 )
213 parameter ( czero = ( 0.0d0, 0.0d0 ) )
216 INTEGER bx, bxst, c, difl, difr, givcol, givnum,
217 $ givptr, i, icmpq1, icmpq2, irwb, irwib, irwrb,
218 $ irwu, irwvt, irwwrk, iwk, j, jcol, jimag,
219 $ jreal, jrow, k, nlvl, nm1, nrwork, nsize, nsub,
220 $ perm, poles, s, sizei, smlszp, sqre, st, st1,
222 DOUBLE PRECISION cs, eps, orgnrm, rcnd, r, sn, tol
235 INTRINSIC abs, dble, dcmplx, dimag, int, log, sign
245 ELSE IF( nrhs.LT.1 )
THEN
247 ELSE IF( ( ldb.LT.1 ) .OR. ( ldb.LT.n ) )
THEN
251 CALL xerbla(
'ZLALSD', -info )
259 IF( ( rcond.LE.zero ) .OR. ( rcond.GE.one ) )
THEN
271 ELSE IF( n.EQ.1 )
THEN
272 IF( d( 1 ).EQ.zero )
THEN
273 CALL zlaset(
'A', 1, nrhs, czero, czero, b, ldb )
276 CALL zlascl(
'G', 0, 0, d( 1 ), one, 1, nrhs, b, ldb, info )
277 d( 1 ) = abs( d( 1 ) )
284 IF( uplo.EQ.
'L' )
THEN
286 CALL dlartg( d( i ), e( i ), cs, sn, r )
289 d( i+1 ) = cs*d( i+1 )
291 CALL zdrot( 1, b( i, 1 ), 1, b( i+1, 1 ), 1, cs, sn )
302 CALL zdrot( 1, b( j, i ), 1, b( j+1, i ), 1, cs, sn )
311 orgnrm =
dlanst(
'M', n, d, e )
312 IF( orgnrm.EQ.zero )
THEN
313 CALL zlaset(
'A', n, nrhs, czero, czero, b, ldb )
317 CALL dlascl(
'G', 0, 0, orgnrm, one, n, 1, d, n, info )
318 CALL dlascl(
'G', 0, 0, orgnrm, one, nm1, 1, e, nm1, info )
323 IF( n.LE.smlsiz )
THEN
328 irwib = irwrb + n*nrhs
329 irwb = irwib + n*nrhs
330 CALL dlaset(
'A', n, n, zero, one, rwork( irwu ), n )
331 CALL dlaset(
'A', n, n, zero, one, rwork( irwvt ), n )
332 CALL dlasdq(
'U', 0, n, n, n, 0, d, e, rwork( irwvt ), n,
333 $ rwork( irwu ), n, rwork( irwwrk ), 1,
334 $ rwork( irwwrk ), info )
347 rwork( j ) = dble( b( jrow, jcol ) )
350 CALL dgemm(
'T',
'N', n, nrhs, n, one, rwork( irwu ), n,
351 $ rwork( irwb ), n, zero, rwork( irwrb ), n )
356 rwork( j ) = dimag( b( jrow, jcol ) )
359 CALL dgemm(
'T',
'N', n, nrhs, n, one, rwork( irwu ), n,
360 $ rwork( irwb ), n, zero, rwork( irwib ), n )
367 b( jrow, jcol ) = dcmplx( rwork( jreal ),
372 tol = rcnd*abs( d(
idamax( n, d, 1 ) ) )
374 IF( d( i ).LE.tol )
THEN
375 CALL zlaset(
'A', 1, nrhs, czero, czero, b( i, 1 ), ldb )
377 CALL zlascl(
'G', 0, 0, d( i ), one, 1, nrhs, b( i, 1 ),
391 DO 120 jcol = 1, nrhs
394 rwork( j ) = dble( b( jrow, jcol ) )
397 CALL dgemm(
'T',
'N', n, nrhs, n, one, rwork( irwvt ), n,
398 $ rwork( irwb ), n, zero, rwork( irwrb ), n )
400 DO 140 jcol = 1, nrhs
403 rwork( j ) = dimag( b( jrow, jcol ) )
406 CALL dgemm(
'T',
'N', n, nrhs, n, one, rwork( irwvt ), n,
407 $ rwork( irwb ), n, zero, rwork( irwib ), n )
410 DO 160 jcol = 1, nrhs
414 b( jrow, jcol ) = dcmplx( rwork( jreal ),
421 CALL dlascl(
'G', 0, 0, one, orgnrm, n, 1, d, n, info )
422 CALL dlasrt(
'D', n, d, info )
423 CALL zlascl(
'G', 0, 0, orgnrm, one, n, nrhs, b, ldb, info )
430 nlvl = int( log( dble( n ) / dble( smlsiz+1 ) ) / log( two ) ) + 1
442 givnum = poles + 2*nlvl*n
443 nrwork = givnum + 2*nlvl*n
447 irwib = irwrb + smlsiz*nrhs
448 irwb = irwib + smlsiz*nrhs
454 givcol = perm + nlvl*n
455 iwk = givcol + nlvl*n*2
464 IF( abs( d( i ) ).LT.eps )
THEN
465 d( i ) = sign( eps, d( i ) )
470 IF( ( abs( e( i ) ).LT.eps ) .OR. ( i.EQ.nm1 ) )
THEN
482 iwork( sizei+nsub-1 ) = nsize
483 ELSE IF( abs( e( i ) ).GE.eps )
THEN
488 iwork( sizei+nsub-1 ) = nsize
496 iwork( sizei+nsub-1 ) = nsize
499 iwork( sizei+nsub-1 ) = 1
500 CALL zcopy( nrhs, b( n, 1 ), ldb, work( bx+nm1 ), n )
503 IF( nsize.EQ.1 )
THEN
508 CALL zcopy( nrhs, b( st, 1 ), ldb, work( bx+st1 ), n )
509 ELSE IF( nsize.LE.smlsiz )
THEN
513 CALL dlaset(
'A', nsize, nsize, zero, one,
514 $ rwork( vt+st1 ), n )
515 CALL dlaset(
'A', nsize, nsize, zero, one,
516 $ rwork( u+st1 ), n )
517 CALL dlasdq(
'U', 0, nsize, nsize, nsize, 0, d( st ),
518 $ e( st ), rwork( vt+st1 ), n, rwork( u+st1 ),
519 $ n, rwork( nrwork ), 1, rwork( nrwork ),
530 DO 190 jcol = 1, nrhs
531 DO 180 jrow = st, st + nsize - 1
533 rwork( j ) = dble( b( jrow, jcol ) )
536 CALL dgemm(
'T',
'N', nsize, nrhs, nsize, one,
537 $ rwork( u+st1 ), n, rwork( irwb ), nsize,
538 $ zero, rwork( irwrb ), nsize )
540 DO 210 jcol = 1, nrhs
541 DO 200 jrow = st, st + nsize - 1
543 rwork( j ) = dimag( b( jrow, jcol ) )
546 CALL dgemm(
'T',
'N', nsize, nrhs, nsize, one,
547 $ rwork( u+st1 ), n, rwork( irwb ), nsize,
548 $ zero, rwork( irwib ), nsize )
551 DO 230 jcol = 1, nrhs
552 DO 220 jrow = st, st + nsize - 1
555 b( jrow, jcol ) = dcmplx( rwork( jreal ),
560 CALL zlacpy(
'A', nsize, nrhs, b( st, 1 ), ldb,
561 $ work( bx+st1 ), n )
566 CALL dlasda( icmpq1, smlsiz, nsize, sqre, d( st ),
567 $ e( st ), rwork( u+st1 ), n, rwork( vt+st1 ),
568 $ iwork( k+st1 ), rwork( difl+st1 ),
569 $ rwork( difr+st1 ), rwork( z+st1 ),
570 $ rwork( poles+st1 ), iwork( givptr+st1 ),
571 $ iwork( givcol+st1 ), n, iwork( perm+st1 ),
572 $ rwork( givnum+st1 ), rwork( c+st1 ),
573 $ rwork( s+st1 ), rwork( nrwork ),
574 $ iwork( iwk ), info )
579 CALL zlalsa( icmpq2, smlsiz, nsize, nrhs, b( st, 1 ),
580 $ ldb, work( bxst ), n, rwork( u+st1 ), n,
581 $ rwork( vt+st1 ), iwork( k+st1 ),
582 $ rwork( difl+st1 ), rwork( difr+st1 ),
583 $ rwork( z+st1 ), rwork( poles+st1 ),
584 $ iwork( givptr+st1 ), iwork( givcol+st1 ), n,
585 $ iwork( perm+st1 ), rwork( givnum+st1 ),
586 $ rwork( c+st1 ), rwork( s+st1 ),
587 $ rwork( nrwork ), iwork( iwk ), info )
598 tol = rcnd*abs( d(
idamax( n, d, 1 ) ) )
605 IF( abs( d( i ) ).LE.tol )
THEN
606 CALL zlaset(
'A', 1, nrhs, czero, czero, work( bx+i-1 ), n )
609 CALL zlascl(
'G', 0, 0, d( i ), one, 1, nrhs,
610 $ work( bx+i-1 ), n, info )
612 d( i ) = abs( d( i ) )
621 nsize = iwork( sizei+i-1 )
623 IF( nsize.EQ.1 )
THEN
624 CALL zcopy( nrhs, work( bxst ), n, b( st, 1 ), ldb )
625 ELSE IF( nsize.LE.smlsiz )
THEN
636 DO 270 jcol = 1, nrhs
638 DO 260 jrow = 1, nsize
640 rwork( jreal ) = dble( work( j+jrow ) )
643 CALL dgemm(
'T',
'N', nsize, nrhs, nsize, one,
644 $ rwork( vt+st1 ), n, rwork( irwb ), nsize, zero,
645 $ rwork( irwrb ), nsize )
648 DO 290 jcol = 1, nrhs
650 DO 280 jrow = 1, nsize
652 rwork( jimag ) = dimag( work( j+jrow ) )
655 CALL dgemm(
'T',
'N', nsize, nrhs, nsize, one,
656 $ rwork( vt+st1 ), n, rwork( irwb ), nsize, zero,
657 $ rwork( irwib ), nsize )
660 DO 310 jcol = 1, nrhs
661 DO 300 jrow = st, st + nsize - 1
664 b( jrow, jcol ) = dcmplx( rwork( jreal ),
669 CALL zlalsa( icmpq2, smlsiz, nsize, nrhs, work( bxst ), n,
670 $ b( st, 1 ), ldb, rwork( u+st1 ), n,
671 $ rwork( vt+st1 ), iwork( k+st1 ),
672 $ rwork( difl+st1 ), rwork( difr+st1 ),
673 $ rwork( z+st1 ), rwork( poles+st1 ),
674 $ iwork( givptr+st1 ), iwork( givcol+st1 ), n,
675 $ iwork( perm+st1 ), rwork( givnum+st1 ),
676 $ rwork( c+st1 ), rwork( s+st1 ),
677 $ rwork( nrwork ), iwork( iwk ), info )
686 CALL dlascl(
'G', 0, 0, one, orgnrm, n, 1, d, n, info )
687 CALL dlasrt(
'D', n, d, info )
688 CALL zlascl(
'G', 0, 0, orgnrm, one, n, nrhs, b, ldb, info )
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 dlasrt(ID, N, D, INFO)
DLASRT sorts numbers in increasing or decreasing order.
subroutine zlacpy(UPLO, M, N, A, LDA, B, LDB)
ZLACPY copies all or part of one two-dimensional array to another.
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.
integer function idamax(N, DX, INCX)
IDAMAX
subroutine zcopy(N, ZX, INCX, ZY, INCY)
ZCOPY
double precision function dlamch(CMACH)
DLAMCH
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 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 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 dgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
DGEMM
subroutine xerbla(SRNAME, INFO)
XERBLA
subroutine zdrot(N, CX, INCX, CY, INCY, C, S)
ZDROT
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...
double precision function dlanst(NORM, N, D, E)
DLANST 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.
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 dlartg(F, G, CS, SN, R)
DLARTG generates a plane rotation with real cosine and real sine.