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
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 )
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 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 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
subroutine dcopy(N, DX, INCX, DY, INCY)
DCOPY
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 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 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.
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.
double precision function dlamch(CMACH)
DLAMCH
set ue cd $ADTTMP cat<< EOF > tmp f Program LinearEquations Implicit none Real j
subroutine drot(N, DX, INCX, DY, INCY, C, S)
DROT
integer function idamax(N, DX, INCX)
IDAMAX
subroutine dlacpy(UPLO, M, N, A, LDA, B, LDB)
DLACPY copies all or part of one two-dimensional array to another.
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...