dlasd1()

subroutine dlasd1	(	integer	nl,
		integer	nr,
		integer	sqre,
		double precision, dimension( * )	d,
		double precision	alpha,
		double precision	beta,
		double precision, dimension( ldu, * )	u,
		integer	ldu,
		double precision, dimension( ldvt, * )	vt,
		integer	ldvt,
		integer, dimension( * )	idxq,
		integer, dimension( * )	iwork,
		double precision, dimension( * )	work,
		integer	info )

DLASD1 computes the SVD of an upper bidiagonal matrix B of the specified size. Used by sbdsdc.

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

Purpose:

!>
!> DLASD1 computes the SVD of an upper bidiagonal N-by-M matrix B,
!> where N = NL + NR + 1 and M = N + SQRE. DLASD1 is called from DLASD0.
!>
!> A related subroutine DLASD7 handles the case in which the singular
!> values (and the singular vectors in factored form) are desired.
!>
!> DLASD1 computes the SVD as follows:
!>
!>               ( D1(in)    0    0       0 )
!>   B = U(in) * (   Z1**T   a   Z2**T    b ) * VT(in)
!>               (   0       0   D2(in)   0 )
!>
!>     = U(out) * ( D(out) 0) * VT(out)
!>
!> where Z**T = (Z1**T a Z2**T b) = u**T VT**T, and u is a vector of dimension M
!> with ALPHA and BETA in the NL+1 and NL+2 th entries and zeros
!> elsewhere; and the entry b is empty if SQRE = 0.
!>
!> The left singular vectors of the original matrix are stored in U, and
!> the transpose of the right singular vectors are stored in VT, and the
!> singular values are in D.  The algorithm consists of three stages:
!>
!>    The first stage consists of deflating the size of the problem
!>    when there are multiple singular values or when there are zeros in
!>    the Z vector.  For each such occurrence the dimension of the
!>    secular equation problem is reduced by one.  This stage is
!>    performed by the routine DLASD2.
!>
!>    The second stage consists of calculating the updated
!>    singular values. This is done by finding the square roots of the
!>    roots of the secular equation via the routine DLASD4 (as called
!>    by DLASD3). This routine also calculates the singular vectors of
!>    the current problem.
!>
!>    The final stage consists of computing the updated singular vectors
!>    directly using the updated singular values.  The singular vectors
!>    for the current problem are multiplied with the singular vectors
!>    from the overall problem.
!> 

Parameters

[in]	NL	!> NL is INTEGER !> The row dimension of the upper block. NL >= 1. !>
[in]	NR	!> NR is INTEGER !> The row dimension of the lower block. NR >= 1. !>
[in]	SQRE	!> SQRE is INTEGER !> = 0: the lower block is an NR-by-NR square matrix. !> = 1: the lower block is an NR-by-(NR+1) rectangular matrix. !> !> The bidiagonal matrix has row dimension N = NL + NR + 1, !> and column dimension M = N + SQRE. !>
[in,out]	D	!> D is DOUBLE PRECISION array, !> dimension (N = NL+NR+1). !> On entry D(1:NL,1:NL) contains the singular values of the !> upper block; and D(NL+2:N) contains the singular values of !> the lower block. On exit D(1:N) contains the singular values !> of the modified matrix. !>
[in,out]	ALPHA	!> ALPHA is DOUBLE PRECISION !> Contains the diagonal element associated with the added row. !>
[in,out]	BETA	!> BETA is DOUBLE PRECISION !> Contains the off-diagonal element associated with the added !> row. !>
[in,out]	U	!> U is DOUBLE PRECISION array, dimension(LDU,N) !> On entry U(1:NL, 1:NL) contains the left singular vectors of !> the upper block; U(NL+2:N, NL+2:N) contains the left singular !> vectors of the lower block. On exit U contains the left !> singular vectors of the bidiagonal matrix. !>
[in]	LDU	!> LDU is INTEGER !> The leading dimension of the array U. LDU >= max( 1, N ). !>
[in,out]	VT	!> VT is DOUBLE PRECISION array, dimension(LDVT,M) !> where M = N + SQRE. !> On entry VT(1:NL+1, 1:NL+1)**T contains the right singular !> vectors of the upper block; VT(NL+2:M, NL+2:M)**T contains !> the right singular vectors of the lower block. On exit !> VT**T contains the right singular vectors of the !> bidiagonal matrix. !>
[in]	LDVT	!> LDVT is INTEGER !> The leading dimension of the array VT. LDVT >= max( 1, M ). !>
[in,out]	IDXQ	!> IDXQ is INTEGER array, dimension(N) !> This contains the permutation which will reintegrate the !> subproblem just solved back into sorted order, i.e. !> D( IDXQ( I = 1, N ) ) will be in ascending order. !>
[out]	IWORK	!> IWORK is INTEGER array, dimension( 4 * N ) !>
[out]	WORK	!> WORK is DOUBLE PRECISION array, dimension( 3*M**2 + 2*M ) !>
[out]	INFO	!> INFO is INTEGER !> = 0: successful exit. !> < 0: if INFO = -i, the i-th argument had an illegal value. !> > 0: if INFO = 1, a singular value did not converge !>

Author

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

Contributors:

Ming Gu and Huan Ren, Computer Science Division, University of California at Berkeley, USA

Definition at line 200 of file dlasd1.f.

*
*  -- LAPACK auxiliary routine --
*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
*
*     .. Scalar Arguments ..
      INTEGER            INFO, LDU, LDVT, NL, NR, SQRE
      DOUBLE PRECISION   ALPHA, BETA
*     ..
*     .. Array Arguments ..
      INTEGER            IDXQ( * ), IWORK( * )
      DOUBLE PRECISION   D( * ), U( LDU, * ), VT( LDVT, * ), WORK( * )
*     ..
*
*  =====================================================================
*
*     .. Parameters ..
*
      DOUBLE PRECISION   ONE, ZERO
      parameter( one = 1.0d+0, zero = 0.0d+0 )
*     ..
*     .. Local Scalars ..
      INTEGER            COLTYP, I, IDX, IDXC, IDXP, IQ, ISIGMA, IU2,
     $                   IVT2, IZ, K, LDQ, LDU2, LDVT2, M, N, N1, N2
      DOUBLE PRECISION   ORGNRM
*     ..
*     .. External Subroutines ..
      EXTERNAL           dlamrg, dlascl, dlasd2, dlasd3,
     $                   xerbla
*     ..
*     .. Intrinsic Functions ..
      INTRINSIC          abs, max
*     ..
*     .. Executable Statements ..
*
*     Test the input parameters.
*
      info = 0
*
      IF( nl.LT.1 ) THEN
         info = -1
      ELSE IF( nr.LT.1 ) THEN
         info = -2
      ELSE IF( ( sqre.LT.0 ) .OR. ( sqre.GT.1 ) ) THEN
         info = -3
      END IF
      IF( info.NE.0 ) THEN
         CALL xerbla( 'DLASD1', -info )
         RETURN
      END IF
*
      n = nl + nr + 1
      m = n + sqre
*
*     The following values are for bookkeeping purposes only.  They are
*     integer pointers which indicate the portion of the workspace
*     used by a particular array in DLASD2 and DLASD3.
*
      ldu2 = n
      ldvt2 = m
*
      iz = 1
      isigma = iz + m
      iu2 = isigma + n
      ivt2 = iu2 + ldu2*n
      iq = ivt2 + ldvt2*m
*
      idx = 1
      idxc = idx + n
      coltyp = idxc + n
      idxp = coltyp + n
*
*     Scale.
*
      orgnrm = max( abs( alpha ), abs( beta ) )
      d( nl+1 ) = zero
      DO 10 i = 1, n
         IF( abs( d( i ) ).GT.orgnrm ) THEN
            orgnrm = abs( d( i ) )
         END IF
   10 CONTINUE
      CALL dlascl( 'G', 0, 0, orgnrm, one, n, 1, d, n, info )
      alpha = alpha / orgnrm
      beta = beta / orgnrm
*
*     Deflate singular values.
*
      CALL dlasd2( nl, nr, sqre, k, d, work( iz ), alpha, beta, u,
     $             ldu,
     $             vt, ldvt, work( isigma ), work( iu2 ), ldu2,
     $             work( ivt2 ), ldvt2, iwork( idxp ), iwork( idx ),
     $             iwork( idxc ), idxq, iwork( coltyp ), info )
*
*     Solve Secular Equation and update singular vectors.
*
      ldq = k
      CALL dlasd3( nl, nr, sqre, k, d, work( iq ), ldq,
     $             work( isigma ),
     $             u, ldu, work( iu2 ), ldu2, vt, ldvt, work( ivt2 ),
     $             ldvt2, iwork( idxc ), iwork( coltyp ), work( iz ),
     $             info )
*
*     Report the convergence failure.
*
      IF( info.NE.0 ) THEN
         RETURN
      END IF
*
*     Unscale.
*
      CALL dlascl( 'G', 0, 0, one, orgnrm, n, 1, d, n, info )
*
*     Prepare the IDXQ sorting permutation.
*
      n1 = k
      n2 = n - k
      CALL dlamrg( n1, n2, d, 1, -1, idxq )
*
      RETURN
*
*     End of DLASD1
*

Here is the call graph for this function:

Here is the caller graph for this function: