179 SUBROUTINE dlalsd( UPLO, SMLSIZ, N, NRHS, D, E, B, LDB, RCOND,
180 $ rank, work, iwork, info )
189 INTEGER INFO, LDB, N, NRHS, RANK, SMLSIZ
190 DOUBLE PRECISION RCOND
194 DOUBLE PRECISION B( ldb, * ), D( * ), E( * ), WORK( * )
200 DOUBLE PRECISION ZERO, ONE, TWO
201 parameter ( zero = 0.0d0, one = 1.0d0, two = 2.0d0 )
204 INTEGER BX, BXST, C, DIFL, DIFR, GIVCOL, GIVNUM,
205 $ givptr, i, icmpq1, icmpq2, iwk, j, k, nlvl,
206 $ nm1, nsize, nsub, nwork, perm, poles, s, sizei,
207 $ smlszp, sqre, st, st1, u, vt, z
208 DOUBLE PRECISION CS, EPS, ORGNRM, R, RCND, SN, TOL
212 DOUBLE PRECISION DLAMCH, DLANST
213 EXTERNAL idamax, dlamch, dlanst
220 INTRINSIC abs, dble, int, log, sign
230 ELSE IF( nrhs.LT.1 )
THEN
232 ELSE IF( ( ldb.LT.1 ) .OR. ( ldb.LT.n ) )
THEN
236 CALL xerbla(
'DLALSD', -info )
240 eps = dlamch(
'Epsilon' )
244 IF( ( rcond.LE.zero ) .OR. ( rcond.GE.one ) )
THEN
256 ELSE IF( n.EQ.1 )
THEN
257 IF( d( 1 ).EQ.zero )
THEN
258 CALL dlaset(
'A', 1, nrhs, zero, zero, b, ldb )
261 CALL dlascl(
'G', 0, 0, d( 1 ), one, 1, nrhs, b, ldb, info )
262 d( 1 ) = abs( d( 1 ) )
269 IF( uplo.EQ.
'L' )
THEN
271 CALL dlartg( d( i ), e( i ), cs, sn, r )
274 d( i+1 ) = cs*d( i+1 )
276 CALL drot( 1, b( i, 1 ), 1, b( i+1, 1 ), 1, cs, sn )
287 CALL drot( 1, b( j, i ), 1, b( j+1, i ), 1, cs, sn )
296 orgnrm = dlanst(
'M', n, d, e )
297 IF( orgnrm.EQ.zero )
THEN
298 CALL dlaset(
'A', n, nrhs, zero, zero, b, ldb )
302 CALL dlascl(
'G', 0, 0, orgnrm, one, n, 1, d, n, info )
303 CALL dlascl(
'G', 0, 0, orgnrm, one, nm1, 1, e, nm1, info )
308 IF( n.LE.smlsiz )
THEN
310 CALL dlaset(
'A', n, n, zero, one, work, n )
311 CALL dlasdq(
'U', 0, n, n, 0, nrhs, d, e, work, n, work, n, b,
312 $ ldb, work( nwork ), info )
316 tol = rcnd*abs( d( idamax( n, d, 1 ) ) )
318 IF( d( i ).LE.tol )
THEN
319 CALL dlaset(
'A', 1, nrhs, zero, zero, b( i, 1 ), ldb )
321 CALL dlascl(
'G', 0, 0, d( i ), one, 1, nrhs, b( i, 1 ),
326 CALL dgemm(
'T',
'N', n, nrhs, n, one, work, n, b, ldb, zero,
328 CALL dlacpy(
'A', n, nrhs, work( nwork ), n, b, ldb )
332 CALL dlascl(
'G', 0, 0, one, orgnrm, n, 1, d, n, info )
333 CALL dlasrt(
'D', n, d, info )
334 CALL dlascl(
'G', 0, 0, orgnrm, one, n, nrhs, b, ldb, info )
341 nlvl = int( log( dble( n ) / dble( smlsiz+1 ) ) / log( two ) ) + 1
353 givnum = poles + 2*nlvl*n
354 bx = givnum + 2*nlvl*n
361 givcol = perm + nlvl*n
362 iwk = givcol + nlvl*n*2
371 IF( abs( d( i ) ).LT.eps )
THEN
372 d( i ) = sign( eps, d( i ) )
377 IF( ( abs( e( i ) ).LT.eps ) .OR. ( i.EQ.nm1 ) )
THEN
389 iwork( sizei+nsub-1 ) = nsize
390 ELSE IF( abs( e( i ) ).GE.eps )
THEN
395 iwork( sizei+nsub-1 ) = nsize
403 iwork( sizei+nsub-1 ) = nsize
406 iwork( sizei+nsub-1 ) = 1
407 CALL dcopy( nrhs, b( n, 1 ), ldb, work( bx+nm1 ), n )
410 IF( nsize.EQ.1 )
THEN
415 CALL dcopy( nrhs, b( st, 1 ), ldb, work( bx+st1 ), n )
416 ELSE IF( nsize.LE.smlsiz )
THEN
420 CALL dlaset(
'A', nsize, nsize, zero, one,
421 $ work( vt+st1 ), n )
422 CALL dlasdq(
'U', 0, nsize, nsize, 0, nrhs, d( st ),
423 $ e( st ), work( vt+st1 ), n, work( nwork ),
424 $ n, b( st, 1 ), ldb, work( nwork ), info )
428 CALL dlacpy(
'A', nsize, nrhs, b( st, 1 ), ldb,
429 $ work( bx+st1 ), n )
434 CALL dlasda( icmpq1, smlsiz, nsize, sqre, d( st ),
435 $ e( st ), work( u+st1 ), n, work( vt+st1 ),
436 $ iwork( k+st1 ), work( difl+st1 ),
437 $ work( difr+st1 ), work( z+st1 ),
438 $ work( poles+st1 ), iwork( givptr+st1 ),
439 $ iwork( givcol+st1 ), n, iwork( perm+st1 ),
440 $ work( givnum+st1 ), work( c+st1 ),
441 $ work( s+st1 ), work( nwork ), iwork( iwk ),
447 CALL dlalsa( icmpq2, smlsiz, nsize, nrhs, b( st, 1 ),
448 $ ldb, work( bxst ), n, work( u+st1 ), n,
449 $ work( vt+st1 ), iwork( k+st1 ),
450 $ work( difl+st1 ), work( difr+st1 ),
451 $ work( z+st1 ), work( poles+st1 ),
452 $ iwork( givptr+st1 ), iwork( givcol+st1 ), n,
453 $ iwork( perm+st1 ), work( givnum+st1 ),
454 $ work( c+st1 ), work( s+st1 ), work( nwork ),
455 $ iwork( iwk ), info )
466 tol = rcnd*abs( d( idamax( n, d, 1 ) ) )
473 IF( abs( d( i ) ).LE.tol )
THEN
474 CALL dlaset(
'A', 1, nrhs, zero, zero, work( bx+i-1 ), n )
477 CALL dlascl(
'G', 0, 0, d( i ), one, 1, nrhs,
478 $ work( bx+i-1 ), n, info )
480 d( i ) = abs( d( i ) )
489 nsize = iwork( sizei+i-1 )
491 IF( nsize.EQ.1 )
THEN
492 CALL dcopy( nrhs, work( bxst ), n, b( st, 1 ), ldb )
493 ELSE IF( nsize.LE.smlsiz )
THEN
494 CALL dgemm(
'T',
'N', nsize, nrhs, nsize, one,
495 $ work( vt+st1 ), n, work( bxst ), n, zero,
498 CALL dlalsa( icmpq2, smlsiz, nsize, nrhs, work( bxst ), n,
499 $ b( st, 1 ), ldb, work( u+st1 ), n,
500 $ work( vt+st1 ), iwork( k+st1 ),
501 $ work( difl+st1 ), work( difr+st1 ),
502 $ work( z+st1 ), work( poles+st1 ),
503 $ iwork( givptr+st1 ), iwork( givcol+st1 ), n,
504 $ iwork( perm+st1 ), work( givnum+st1 ),
505 $ work( c+st1 ), work( s+st1 ), work( nwork ),
506 $ iwork( iwk ), info )
515 CALL dlascl(
'G', 0, 0, one, orgnrm, n, 1, d, n, info )
516 CALL dlasrt(
'D', n, d, info )
517 CALL dlascl(
'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 dcopy(N, DX, INCX, DY, INCY)
DCOPY
subroutine dlacpy(UPLO, M, N, A, LDA, B, LDB)
DLACPY copies all or part of one two-dimensional array to another.
subroutine drot(N, DX, INCX, DY, INCY, C, S)
DROT
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 dlalsd(UPLO, SMLSIZ, N, NRHS, D, E, B, LDB, RCOND, RANK, WORK, IWORK, INFO)
DLALSD uses the singular value decomposition of A to solve the least squares problem.
subroutine dgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
DGEMM
subroutine xerbla(SRNAME, INFO)
XERBLA
subroutine dlalsa(ICOMPQ, SMLSIZ, N, NRHS, B, LDB, BX, LDBX, U, LDU, VT, K, DIFL, DIFR, Z, POLES, GIVPTR, GIVCOL, LDGCOL, PERM, GIVNUM, C, S, WORK, IWORK, INFO)
DLALSA computes the SVD of the coefficient matrix in compact form. Used by sgelsd.
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 dlartg(F, G, CS, SN, R)
DLARTG generates a plane rotation with real cosine and real sine.