dlalsd()

subroutine dlalsd	(	character	uplo,
		integer	smlsiz,
		integer	n,
		integer	nrhs,
		double precision, dimension( * )	d,
		double precision, dimension( * )	e,
		double precision, dimension( ldb, * )	b,
		integer	ldb,
		double precision	rcond,
		integer	rank,
		double precision, dimension( * )	work,
		integer, dimension( * )	iwork,
		integer	info )

DLALSD uses the singular value decomposition of A to solve the least squares problem.

Download DLALSD + dependencies [TGZ] [ZIP] [TXT]

Purpose:

!>
!> DLALSD uses the singular value decomposition of A to solve the least
!> squares problem of finding X to minimize the Euclidean norm of each
!> column of A*X-B, where A is N-by-N upper bidiagonal, and X and B
!> are N-by-NRHS. The solution X overwrites B.
!>
!> The singular values of A smaller than RCOND times the largest
!> singular value are treated as zero in solving the least squares
!> problem; in this case a minimum norm solution is returned.
!> The actual singular values are returned in D in ascending order.
!>
!> 

Parameters

[in]	UPLO	!> UPLO is CHARACTER*1 !> = 'U': D and E define an upper bidiagonal matrix. !> = 'L': D and E define a lower bidiagonal matrix. !>
[in]	SMLSIZ	!> SMLSIZ is INTEGER !> The maximum size of the subproblems at the bottom of the !> computation tree. !>
[in]	N	!> N is INTEGER !> The dimension of the bidiagonal matrix. N >= 0. !>
[in]	NRHS	!> NRHS is INTEGER !> The number of columns of B. NRHS must be at least 1. !>
[in,out]	D	!> D is DOUBLE PRECISION array, dimension (N) !> On entry D contains the main diagonal of the bidiagonal !> matrix. On exit, if INFO = 0, D contains its singular values. !>
[in,out]	E	!> E is DOUBLE PRECISION array, dimension (N-1) !> Contains the super-diagonal entries of the bidiagonal matrix. !> On exit, E has been destroyed. !>
[in,out]	B	!> B is DOUBLE PRECISION array, dimension (LDB,NRHS) !> On input, B contains the right hand sides of the least !> squares problem. On output, B contains the solution X. !>
[in]	LDB	!> LDB is INTEGER !> The leading dimension of B in the calling subprogram. !> LDB must be at least max(1,N). !>
[in]	RCOND	!> RCOND is DOUBLE PRECISION !> The singular values of A less than or equal to RCOND times !> the largest singular value are treated as zero in solving !> the least squares problem. If RCOND is negative, !> machine precision is used instead. !> For example, if diag(S)*X=B were the least squares problem, !> where diag(S) is a diagonal matrix of singular values, the !> solution would be X(i) = B(i) / S(i) if S(i) is greater than !> RCOND*max(S), and X(i) = 0 if S(i) is less than or equal to !> RCOND*max(S). !>
[out]	RANK	!> RANK is INTEGER !> The number of singular values of A greater than RCOND times !> the largest singular value. !>
[out]	WORK	!> WORK is DOUBLE PRECISION array, dimension at least !> (9*N + 2*N*SMLSIZ + 8*N*NLVL + N*NRHS + (SMLSIZ+1)**2), !> where NLVL = max(0, INT(log_2 (N/(SMLSIZ+1))) + 1). !>
[out]	IWORK	!> IWORK is INTEGER array, dimension at least !> (3*N*NLVL + 11*N) !>
[out]	INFO	!> INFO is INTEGER !> = 0: successful exit. !> < 0: if INFO = -i, the i-th argument had an illegal value. !> > 0: The algorithm failed to compute a singular value while !> working on the submatrix lying in rows and columns !> INFO/(N+1) through MOD(INFO,N+1). !>

Author

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

Contributors:

Ming Gu and Ren-Cang Li, Computer Science Division, University of California at Berkeley, USA
Osni Marques, LBNL/NERSC, USA

Definition at line 169 of file dlalsd.f.

*
*  -- LAPACK computational routine --
*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
*
*     .. Scalar Arguments ..
      CHARACTER          UPLO
      INTEGER            INFO, LDB, N, NRHS, RANK, SMLSIZ
      DOUBLE PRECISION   RCOND
*     ..
*     .. Array Arguments ..
      INTEGER            IWORK( * )
      DOUBLE PRECISION   B( LDB, * ), D( * ), E( * ), WORK( * )
*     ..
*
*  =====================================================================
*
*     .. Parameters ..
      DOUBLE PRECISION   ZERO, ONE, TWO
      parameter( zero = 0.0d0, one = 1.0d0, two = 2.0d0 )
*     ..
*     .. Local Scalars ..
      INTEGER            BX, BXST, C, DIFL, DIFR, GIVCOL, GIVNUM,
     $                   GIVPTR, I, ICMPQ1, ICMPQ2, IWK, J, K, NLVL,
     $                   NM1, NSIZE, NSUB, NWORK, PERM, POLES, S, SIZEI,
     $                   SMLSZP, SQRE, ST, ST1, U, VT, Z
      DOUBLE PRECISION   CS, EPS, ORGNRM, R, RCND, SN, TOL
*     ..
*     .. External Functions ..
      INTEGER            IDAMAX
      DOUBLE PRECISION   DLAMCH, DLANST
      EXTERNAL           idamax, dlamch, dlanst
*     ..
*     .. External Subroutines ..
      EXTERNAL           dcopy, dgemm, dlacpy, dlalsa, dlartg,
     $                   dlascl,
     $                   dlasda, dlasdq, dlaset, dlasrt, drot, xerbla
*     ..
*     .. Intrinsic Functions ..
      INTRINSIC          abs, dble, int, log, sign
*     ..
*     .. Executable Statements ..
*
*     Test the input parameters.
*
      info = 0
*
      IF( n.LT.0 ) THEN
         info = -3
      ELSE IF( nrhs.LT.1 ) THEN
         info = -4
      ELSE IF( ( ldb.LT.1 ) .OR. ( ldb.LT.n ) ) THEN
         info = -8
      END IF
      IF( info.NE.0 ) THEN
         CALL xerbla( 'DLALSD', -info )
         RETURN
      END IF
*
      eps = dlamch( 'Epsilon' )
*
*     Set up the tolerance.
*
      IF( ( rcond.LE.zero ) .OR. ( rcond.GE.one ) ) THEN
         rcnd = eps
      ELSE
         rcnd = rcond
      END IF
*
      rank = 0
*
*     Quick return if possible.
*
      IF( n.EQ.0 ) THEN
         RETURN
      ELSE IF( n.EQ.1 ) THEN
         IF( d( 1 ).EQ.zero ) THEN
            CALL dlaset( 'A', 1, nrhs, zero, zero, b, ldb )
         ELSE
            rank = 1
            CALL dlascl( 'G', 0, 0, d( 1 ), one, 1, nrhs, b, ldb,
     $                   info )
            d( 1 ) = abs( d( 1 ) )
         END IF
         RETURN
      END IF
*
*     Rotate the matrix if it is lower bidiagonal.
*
      IF( uplo.EQ.'L' ) THEN
         DO 10 i = 1, n - 1
            CALL dlartg( d( i ), e( i ), cs, sn, r )
            d( i ) = r
            e( i ) = sn*d( i+1 )
            d( i+1 ) = cs*d( i+1 )
            IF( nrhs.EQ.1 ) THEN
               CALL drot( 1, b( i, 1 ), 1, b( i+1, 1 ), 1, cs, sn )
            ELSE
               work( i*2-1 ) = cs
               work( i*2 ) = sn
            END IF
   10    CONTINUE
         IF( nrhs.GT.1 ) THEN
            DO 30 i = 1, nrhs
               DO 20 j = 1, n - 1
                  cs = work( j*2-1 )
                  sn = work( j*2 )
                  CALL drot( 1, b( j, i ), 1, b( j+1, i ), 1, cs,
     $                       sn )
   20          CONTINUE
   30       CONTINUE
         END IF
      END IF
*
*     Scale.
*
      nm1 = n - 1
      orgnrm = dlanst( 'M', n, d, e )
      IF( orgnrm.EQ.zero ) THEN
         CALL dlaset( 'A', n, nrhs, zero, zero, b, ldb )
         RETURN
      END IF
*
      CALL dlascl( 'G', 0, 0, orgnrm, one, n, 1, d, n, info )
      CALL dlascl( 'G', 0, 0, orgnrm, one, nm1, 1, e, nm1, info )
*
*     If N is smaller than the minimum divide size SMLSIZ, then solve
*     the problem with another solver.
*
      IF( n.LE.smlsiz ) THEN
         nwork = 1 + n*n
         CALL dlaset( 'A', n, n, zero, one, work, n )
         CALL dlasdq( 'U', 0, n, n, 0, nrhs, d, e, work, n, work, n,
     $                b,
     $                ldb, work( nwork ), info )
         IF( info.NE.0 ) THEN
            RETURN
         END IF
         tol = rcnd*abs( d( idamax( n, d, 1 ) ) )
         DO 40 i = 1, n
            IF( d( i ).LE.tol ) THEN
               CALL dlaset( 'A', 1, nrhs, zero, zero, b( i, 1 ),
     $                      ldb )
            ELSE
               CALL dlascl( 'G', 0, 0, d( i ), one, 1, nrhs, b( i,
     $                      1 ),
     $                      ldb, info )
               rank = rank + 1
            END IF
   40    CONTINUE
         CALL dgemm( 'T', 'N', n, nrhs, n, one, work, n, b, ldb,
     $               zero,
     $               work( nwork ), n )
         CALL dlacpy( 'A', n, nrhs, work( nwork ), n, b, ldb )
*
*        Unscale.
*
         CALL dlascl( 'G', 0, 0, one, orgnrm, n, 1, d, n, info )
         CALL dlasrt( 'D', n, d, info )
         CALL dlascl( 'G', 0, 0, orgnrm, one, n, nrhs, b, ldb, info )
*
         RETURN
      END IF
*
*     Book-keeping and setting up some constants.
*
      nlvl = int( log( dble( n ) / dble( smlsiz+1 ) ) / log( two ) ) + 1
*
      smlszp = smlsiz + 1
*
      u = 1
      vt = 1 + smlsiz*n
      difl = vt + smlszp*n
      difr = difl + nlvl*n
      z = difr + nlvl*n*2
      c = z + nlvl*n
      s = c + n
      poles = s + n
      givnum = poles + 2*nlvl*n
      bx = givnum + 2*nlvl*n
      nwork = bx + n*nrhs
*
      sizei = 1 + n
      k = sizei + n
      givptr = k + n
      perm = givptr + n
      givcol = perm + nlvl*n
      iwk = givcol + nlvl*n*2
*
      st = 1
      sqre = 0
      icmpq1 = 1
      icmpq2 = 0
      nsub = 0
*
      DO 50 i = 1, n
         IF( abs( d( i ) ).LT.eps ) THEN
            d( i ) = sign( eps, d( i ) )
         END IF
   50 CONTINUE
*
      DO 60 i = 1, nm1
         IF( ( abs( e( i ) ).LT.eps ) .OR. ( i.EQ.nm1 ) ) THEN
            nsub = nsub + 1
            iwork( nsub ) = st
*
*           Subproblem found. First determine its size and then
*           apply divide and conquer on it.
*
            IF( i.LT.nm1 ) THEN
*
*              A subproblem with E(I) small for I < NM1.
*
               nsize = i - st + 1
               iwork( sizei+nsub-1 ) = nsize
            ELSE IF( abs( e( i ) ).GE.eps ) THEN
*
*              A subproblem with E(NM1) not too small but I = NM1.
*
               nsize = n - st + 1
               iwork( sizei+nsub-1 ) = nsize
            ELSE
*
*              A subproblem with E(NM1) small. This implies an
*              1-by-1 subproblem at D(N), which is not solved
*              explicitly.
*
               nsize = i - st + 1
               iwork( sizei+nsub-1 ) = nsize
               nsub = nsub + 1
               iwork( nsub ) = n
               iwork( sizei+nsub-1 ) = 1
               CALL dcopy( nrhs, b( n, 1 ), ldb, work( bx+nm1 ), n )
            END IF
            st1 = st - 1
            IF( nsize.EQ.1 ) THEN
*
*              This is a 1-by-1 subproblem and is not solved
*              explicitly.
*
               CALL dcopy( nrhs, b( st, 1 ), ldb, work( bx+st1 ), n )
            ELSE IF( nsize.LE.smlsiz ) THEN
*
*              This is a small subproblem and is solved by DLASDQ.
*
               CALL dlaset( 'A', nsize, nsize, zero, one,
     $                      work( vt+st1 ), n )
               CALL dlasdq( 'U', 0, nsize, nsize, 0, nrhs, d( st ),
     $                      e( st ), work( vt+st1 ), n, work( nwork ),
     $                      n, b( st, 1 ), ldb, work( nwork ), info )
               IF( info.NE.0 ) THEN
                  RETURN
               END IF
               CALL dlacpy( 'A', nsize, nrhs, b( st, 1 ), ldb,
     $                      work( bx+st1 ), n )
            ELSE
*
*              A large problem. Solve it using divide and conquer.
*
               CALL dlasda( icmpq1, smlsiz, nsize, sqre, d( st ),
     $                      e( st ), work( u+st1 ), n, work( vt+st1 ),
     $                      iwork( k+st1 ), work( difl+st1 ),
     $                      work( difr+st1 ), work( z+st1 ),
     $                      work( poles+st1 ), iwork( givptr+st1 ),
     $                      iwork( givcol+st1 ), n, iwork( perm+st1 ),
     $                      work( givnum+st1 ), work( c+st1 ),
     $                      work( s+st1 ), work( nwork ), iwork( iwk ),
     $                      info )
               IF( info.NE.0 ) THEN
                  RETURN
               END IF
               bxst = bx + st1
               CALL dlalsa( icmpq2, smlsiz, nsize, nrhs, b( st, 1 ),
     $                      ldb, work( bxst ), n, work( u+st1 ), n,
     $                      work( vt+st1 ), iwork( k+st1 ),
     $                      work( difl+st1 ), work( difr+st1 ),
     $                      work( z+st1 ), work( poles+st1 ),
     $                      iwork( givptr+st1 ), iwork( givcol+st1 ), n,
     $                      iwork( perm+st1 ), work( givnum+st1 ),
     $                      work( c+st1 ), work( s+st1 ), work( nwork ),
     $                      iwork( iwk ), info )
               IF( info.NE.0 ) THEN
                  RETURN
               END IF
            END IF
            st = i + 1
         END IF
   60 CONTINUE
*
*     Apply the singular values and treat the tiny ones as zero.
*
      tol = rcnd*abs( d( idamax( n, d, 1 ) ) )
*
      DO 70 i = 1, n
*
*        Some of the elements in D can be negative because 1-by-1
*        subproblems were not solved explicitly.
*
         IF( abs( d( i ) ).LE.tol ) THEN
            CALL dlaset( 'A', 1, nrhs, zero, zero, work( bx+i-1 ),
     $                   n )
         ELSE
            rank = rank + 1
            CALL dlascl( 'G', 0, 0, d( i ), one, 1, nrhs,
     $                   work( bx+i-1 ), n, info )
         END IF
         d( i ) = abs( d( i ) )
   70 CONTINUE
*
*     Now apply back the right singular vectors.
*
      icmpq2 = 1
      DO 80 i = 1, nsub
         st = iwork( i )
         st1 = st - 1
         nsize = iwork( sizei+i-1 )
         bxst = bx + st1
         IF( nsize.EQ.1 ) THEN
            CALL dcopy( nrhs, work( bxst ), n, b( st, 1 ), ldb )
         ELSE IF( nsize.LE.smlsiz ) THEN
            CALL dgemm( 'T', 'N', nsize, nrhs, nsize, one,
     $                  work( vt+st1 ), n, work( bxst ), n, zero,
     $                  b( st, 1 ), ldb )
         ELSE
            CALL dlalsa( icmpq2, smlsiz, nsize, nrhs, work( bxst ),
     $                   n,
     $                   b( st, 1 ), ldb, work( u+st1 ), n,
     $                   work( vt+st1 ), iwork( k+st1 ),
     $                   work( difl+st1 ), work( difr+st1 ),
     $                   work( z+st1 ), work( poles+st1 ),
     $                   iwork( givptr+st1 ), iwork( givcol+st1 ), n,
     $                   iwork( perm+st1 ), work( givnum+st1 ),
     $                   work( c+st1 ), work( s+st1 ), work( nwork ),
     $                   iwork( iwk ), info )
            IF( info.NE.0 ) THEN
               RETURN
            END IF
         END IF
   80 CONTINUE
*
*     Unscale and sort the singular values.
*
      CALL dlascl( 'G', 0, 0, one, orgnrm, n, 1, d, n, info )
      CALL dlasrt( 'D', n, d, info )
      CALL dlascl( 'G', 0, 0, orgnrm, one, n, nrhs, b, ldb, info )
*
      RETURN
*
*     End of DLALSD
*

Here is the call graph for this function:

Here is the caller graph for this function: