188 SUBROUTINE zlalsd( UPLO, SMLSIZ, N, NRHS, D, E, B, LDB, RCOND,
189 $ rank, work, rwork, iwork, info )
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 dlasrt(ID, N, D, INFO)
DLASRT sorts numbers in increasing or decreasing order.
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 zlacpy(UPLO, M, N, A, LDA, B, LDB)
ZLACPY copies all or part of one two-dimensional array to another.
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 zdrot(N, CX, INCX, CY, INCY, C, S)
ZDROT
subroutine dlartg(F, G, CS, SN, R)
DLARTG generates a plane rotation with real cosine and real sine.
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 dgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
DGEMM
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 xerbla(SRNAME, INFO)
XERBLA
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.
set ue cd $ADTTMP cat<< EOF > tmp f Program LinearEquations Implicit none Real b(3) integer i
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 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 zcopy(N, ZX, INCX, ZY, INCY)
ZCOPY
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...
set ue cd $ADTTMP cat<< EOF > tmp f Program LinearEquations Implicit none Real j
DOUBLE PRECISION function dlamch(CMACH)
DLAMCH
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...
INTEGER function idamax(N, DX, INCX)
IDAMAX