*> \brief \b DLASD1 computes the SVD of an upper bidiagonal matrix B of the specified size. Used by sbdsdc. * * =========== DOCUMENTATION =========== * * Online html documentation available at * http://www.netlib.org/lapack/explore-html/ * *> \htmlonly *> Download DLASD1 + dependencies *> *> [TGZ] *> *> [ZIP] *> *> [TXT] *> \endhtmlonly * * Definition: * =========== * * SUBROUTINE DLASD1( NL, NR, SQRE, D, ALPHA, BETA, U, LDU, VT, LDVT, * IDXQ, IWORK, WORK, INFO ) * * .. 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( * ) * .. * * *> \par Purpose: * ============= *> *> \verbatim *> *> 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. *> \endverbatim * * Arguments: * ========== * *> \param[in] NL *> \verbatim *> NL is INTEGER *> The row dimension of the upper block. NL >= 1. *> \endverbatim *> *> \param[in] NR *> \verbatim *> NR is INTEGER *> The row dimension of the lower block. NR >= 1. *> \endverbatim *> *> \param[in] SQRE *> \verbatim *> 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. *> \endverbatim *> *> \param[in,out] D *> \verbatim *> 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. *> \endverbatim *> *> \param[in,out] ALPHA *> \verbatim *> ALPHA is DOUBLE PRECISION *> Contains the diagonal element associated with the added row. *> \endverbatim *> *> \param[in,out] BETA *> \verbatim *> BETA is DOUBLE PRECISION *> Contains the off-diagonal element associated with the added *> row. *> \endverbatim *> *> \param[in,out] U *> \verbatim *> 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. *> \endverbatim *> *> \param[in] LDU *> \verbatim *> LDU is INTEGER *> The leading dimension of the array U. LDU >= max( 1, N ). *> \endverbatim *> *> \param[in,out] VT *> \verbatim *> 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. *> \endverbatim *> *> \param[in] LDVT *> \verbatim *> LDVT is INTEGER *> The leading dimension of the array VT. LDVT >= max( 1, M ). *> \endverbatim *> *> \param[in,out] IDXQ *> \verbatim *> 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. *> \endverbatim *> *> \param[out] IWORK *> \verbatim *> IWORK is INTEGER array, dimension( 4 * N ) *> \endverbatim *> *> \param[out] WORK *> \verbatim *> WORK is DOUBLE PRECISION array, dimension( 3*M**2 + 2*M ) *> \endverbatim *> *> \param[out] INFO *> \verbatim *> 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 *> \endverbatim * * Authors: * ======== * *> \author Univ. of Tennessee *> \author Univ. of California Berkeley *> \author Univ. of Colorado Denver *> \author NAG Ltd. * *> \ingroup OTHERauxiliary * *> \par Contributors: * ================== *> *> Ming Gu and Huan Ren, Computer Science Division, University of *> California at Berkeley, USA *> * ===================================================================== SUBROUTINE DLASD1( NL, NR, SQRE, D, ALPHA, BETA, U, LDU, VT, LDVT, \$ IDXQ, IWORK, WORK, INFO ) * * -- 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 * END