LAPACK  3.4.1
LAPACK: Linear Algebra PACKage
dlasd1.f
Go to the documentation of this file.
00001 *> \brief \b DLASD1
00002 *
00003 *  =========== DOCUMENTATION ===========
00004 *
00005 * Online html documentation available at 
00006 *            http://www.netlib.org/lapack/explore-html/ 
00007 *
00008 *> \htmlonly
00009 *> Download DLASD1 + dependencies 
00010 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dlasd1.f"> 
00011 *> [TGZ]</a> 
00012 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dlasd1.f"> 
00013 *> [ZIP]</a> 
00014 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dlasd1.f"> 
00015 *> [TXT]</a>
00016 *> \endhtmlonly 
00017 *
00018 *  Definition:
00019 *  ===========
00020 *
00021 *       SUBROUTINE DLASD1( NL, NR, SQRE, D, ALPHA, BETA, U, LDU, VT, LDVT,
00022 *                          IDXQ, IWORK, WORK, INFO )
00023 * 
00024 *       .. Scalar Arguments ..
00025 *       INTEGER            INFO, LDU, LDVT, NL, NR, SQRE
00026 *       DOUBLE PRECISION   ALPHA, BETA
00027 *       ..
00028 *       .. Array Arguments ..
00029 *       INTEGER            IDXQ( * ), IWORK( * )
00030 *       DOUBLE PRECISION   D( * ), U( LDU, * ), VT( LDVT, * ), WORK( * )
00031 *       ..
00032 *  
00033 *
00034 *> \par Purpose:
00035 *  =============
00036 *>
00037 *> \verbatim
00038 *>
00039 *> DLASD1 computes the SVD of an upper bidiagonal N-by-M matrix B,
00040 *> where N = NL + NR + 1 and M = N + SQRE. DLASD1 is called from DLASD0.
00041 *>
00042 *> A related subroutine DLASD7 handles the case in which the singular
00043 *> values (and the singular vectors in factored form) are desired.
00044 *>
00045 *> DLASD1 computes the SVD as follows:
00046 *>
00047 *>               ( D1(in)    0    0       0 )
00048 *>   B = U(in) * (   Z1**T   a   Z2**T    b ) * VT(in)
00049 *>               (   0       0   D2(in)   0 )
00050 *>
00051 *>     = U(out) * ( D(out) 0) * VT(out)
00052 *>
00053 *> where Z**T = (Z1**T a Z2**T b) = u**T VT**T, and u is a vector of dimension M
00054 *> with ALPHA and BETA in the NL+1 and NL+2 th entries and zeros
00055 *> elsewhere; and the entry b is empty if SQRE = 0.
00056 *>
00057 *> The left singular vectors of the original matrix are stored in U, and
00058 *> the transpose of the right singular vectors are stored in VT, and the
00059 *> singular values are in D.  The algorithm consists of three stages:
00060 *>
00061 *>    The first stage consists of deflating the size of the problem
00062 *>    when there are multiple singular values or when there are zeros in
00063 *>    the Z vector.  For each such occurence the dimension of the
00064 *>    secular equation problem is reduced by one.  This stage is
00065 *>    performed by the routine DLASD2.
00066 *>
00067 *>    The second stage consists of calculating the updated
00068 *>    singular values. This is done by finding the square roots of the
00069 *>    roots of the secular equation via the routine DLASD4 (as called
00070 *>    by DLASD3). This routine also calculates the singular vectors of
00071 *>    the current problem.
00072 *>
00073 *>    The final stage consists of computing the updated singular vectors
00074 *>    directly using the updated singular values.  The singular vectors
00075 *>    for the current problem are multiplied with the singular vectors
00076 *>    from the overall problem.
00077 *> \endverbatim
00078 *
00079 *  Arguments:
00080 *  ==========
00081 *
00082 *> \param[in] NL
00083 *> \verbatim
00084 *>          NL is INTEGER
00085 *>         The row dimension of the upper block.  NL >= 1.
00086 *> \endverbatim
00087 *>
00088 *> \param[in] NR
00089 *> \verbatim
00090 *>          NR is INTEGER
00091 *>         The row dimension of the lower block.  NR >= 1.
00092 *> \endverbatim
00093 *>
00094 *> \param[in] SQRE
00095 *> \verbatim
00096 *>          SQRE is INTEGER
00097 *>         = 0: the lower block is an NR-by-NR square matrix.
00098 *>         = 1: the lower block is an NR-by-(NR+1) rectangular matrix.
00099 *>
00100 *>         The bidiagonal matrix has row dimension N = NL + NR + 1,
00101 *>         and column dimension M = N + SQRE.
00102 *> \endverbatim
00103 *>
00104 *> \param[in,out] D
00105 *> \verbatim
00106 *>          D is DOUBLE PRECISION array,
00107 *>                        dimension (N = NL+NR+1).
00108 *>         On entry D(1:NL,1:NL) contains the singular values of the
00109 *>         upper block; and D(NL+2:N) contains the singular values of
00110 *>         the lower block. On exit D(1:N) contains the singular values
00111 *>         of the modified matrix.
00112 *> \endverbatim
00113 *>
00114 *> \param[in,out] ALPHA
00115 *> \verbatim
00116 *>          ALPHA is DOUBLE PRECISION
00117 *>         Contains the diagonal element associated with the added row.
00118 *> \endverbatim
00119 *>
00120 *> \param[in,out] BETA
00121 *> \verbatim
00122 *>          BETA is DOUBLE PRECISION
00123 *>         Contains the off-diagonal element associated with the added
00124 *>         row.
00125 *> \endverbatim
00126 *>
00127 *> \param[in,out] U
00128 *> \verbatim
00129 *>          U is DOUBLE PRECISION array, dimension(LDU,N)
00130 *>         On entry U(1:NL, 1:NL) contains the left singular vectors of
00131 *>         the upper block; U(NL+2:N, NL+2:N) contains the left singular
00132 *>         vectors of the lower block. On exit U contains the left
00133 *>         singular vectors of the bidiagonal matrix.
00134 *> \endverbatim
00135 *>
00136 *> \param[in] LDU
00137 *> \verbatim
00138 *>          LDU is INTEGER
00139 *>         The leading dimension of the array U.  LDU >= max( 1, N ).
00140 *> \endverbatim
00141 *>
00142 *> \param[in,out] VT
00143 *> \verbatim
00144 *>          VT is DOUBLE PRECISION array, dimension(LDVT,M)
00145 *>         where M = N + SQRE.
00146 *>         On entry VT(1:NL+1, 1:NL+1)**T contains the right singular
00147 *>         vectors of the upper block; VT(NL+2:M, NL+2:M)**T contains
00148 *>         the right singular vectors of the lower block. On exit
00149 *>         VT**T contains the right singular vectors of the
00150 *>         bidiagonal matrix.
00151 *> \endverbatim
00152 *>
00153 *> \param[in] LDVT
00154 *> \verbatim
00155 *>          LDVT is INTEGER
00156 *>         The leading dimension of the array VT.  LDVT >= max( 1, M ).
00157 *> \endverbatim
00158 *>
00159 *> \param[out] IDXQ
00160 *> \verbatim
00161 *>          IDXQ is INTEGER array, dimension(N)
00162 *>         This contains the permutation which will reintegrate the
00163 *>         subproblem just solved back into sorted order, i.e.
00164 *>         D( IDXQ( I = 1, N ) ) will be in ascending order.
00165 *> \endverbatim
00166 *>
00167 *> \param[out] IWORK
00168 *> \verbatim
00169 *>          IWORK is INTEGER array, dimension( 4 * N )
00170 *> \endverbatim
00171 *>
00172 *> \param[out] WORK
00173 *> \verbatim
00174 *>          WORK is DOUBLE PRECISION array, dimension( 3*M**2 + 2*M )
00175 *> \endverbatim
00176 *>
00177 *> \param[out] INFO
00178 *> \verbatim
00179 *>          INFO is INTEGER
00180 *>          = 0:  successful exit.
00181 *>          < 0:  if INFO = -i, the i-th argument had an illegal value.
00182 *>          > 0:  if INFO = 1, a singular value did not converge
00183 *> \endverbatim
00184 *
00185 *  Authors:
00186 *  ========
00187 *
00188 *> \author Univ. of Tennessee 
00189 *> \author Univ. of California Berkeley 
00190 *> \author Univ. of Colorado Denver 
00191 *> \author NAG Ltd. 
00192 *
00193 *> \date November 2011
00194 *
00195 *> \ingroup auxOTHERauxiliary
00196 *
00197 *> \par Contributors:
00198 *  ==================
00199 *>
00200 *>     Ming Gu and Huan Ren, Computer Science Division, University of
00201 *>     California at Berkeley, USA
00202 *>
00203 *  =====================================================================
00204       SUBROUTINE DLASD1( NL, NR, SQRE, D, ALPHA, BETA, U, LDU, VT, LDVT,
00205      $                   IDXQ, IWORK, WORK, INFO )
00206 *
00207 *  -- LAPACK auxiliary routine (version 3.4.0) --
00208 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
00209 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
00210 *     November 2011
00211 *
00212 *     .. Scalar Arguments ..
00213       INTEGER            INFO, LDU, LDVT, NL, NR, SQRE
00214       DOUBLE PRECISION   ALPHA, BETA
00215 *     ..
00216 *     .. Array Arguments ..
00217       INTEGER            IDXQ( * ), IWORK( * )
00218       DOUBLE PRECISION   D( * ), U( LDU, * ), VT( LDVT, * ), WORK( * )
00219 *     ..
00220 *
00221 *  =====================================================================
00222 *
00223 *     .. Parameters ..
00224 *
00225       DOUBLE PRECISION   ONE, ZERO
00226       PARAMETER          ( ONE = 1.0D+0, ZERO = 0.0D+0 )
00227 *     ..
00228 *     .. Local Scalars ..
00229       INTEGER            COLTYP, I, IDX, IDXC, IDXP, IQ, ISIGMA, IU2,
00230      $                   IVT2, IZ, K, LDQ, LDU2, LDVT2, M, N, N1, N2
00231       DOUBLE PRECISION   ORGNRM
00232 *     ..
00233 *     .. External Subroutines ..
00234       EXTERNAL           DLAMRG, DLASCL, DLASD2, DLASD3, XERBLA
00235 *     ..
00236 *     .. Intrinsic Functions ..
00237       INTRINSIC          ABS, MAX
00238 *     ..
00239 *     .. Executable Statements ..
00240 *
00241 *     Test the input parameters.
00242 *
00243       INFO = 0
00244 *
00245       IF( NL.LT.1 ) THEN
00246          INFO = -1
00247       ELSE IF( NR.LT.1 ) THEN
00248          INFO = -2
00249       ELSE IF( ( SQRE.LT.0 ) .OR. ( SQRE.GT.1 ) ) THEN
00250          INFO = -3
00251       END IF
00252       IF( INFO.NE.0 ) THEN
00253          CALL XERBLA( 'DLASD1', -INFO )
00254          RETURN
00255       END IF
00256 *
00257       N = NL + NR + 1
00258       M = N + SQRE
00259 *
00260 *     The following values are for bookkeeping purposes only.  They are
00261 *     integer pointers which indicate the portion of the workspace
00262 *     used by a particular array in DLASD2 and DLASD3.
00263 *
00264       LDU2 = N
00265       LDVT2 = M
00266 *
00267       IZ = 1
00268       ISIGMA = IZ + M
00269       IU2 = ISIGMA + N
00270       IVT2 = IU2 + LDU2*N
00271       IQ = IVT2 + LDVT2*M
00272 *
00273       IDX = 1
00274       IDXC = IDX + N
00275       COLTYP = IDXC + N
00276       IDXP = COLTYP + N
00277 *
00278 *     Scale.
00279 *
00280       ORGNRM = MAX( ABS( ALPHA ), ABS( BETA ) )
00281       D( NL+1 ) = ZERO
00282       DO 10 I = 1, N
00283          IF( ABS( D( I ) ).GT.ORGNRM ) THEN
00284             ORGNRM = ABS( D( I ) )
00285          END IF
00286    10 CONTINUE
00287       CALL DLASCL( 'G', 0, 0, ORGNRM, ONE, N, 1, D, N, INFO )
00288       ALPHA = ALPHA / ORGNRM
00289       BETA = BETA / ORGNRM
00290 *
00291 *     Deflate singular values.
00292 *
00293       CALL DLASD2( NL, NR, SQRE, K, D, WORK( IZ ), ALPHA, BETA, U, LDU,
00294      $             VT, LDVT, WORK( ISIGMA ), WORK( IU2 ), LDU2,
00295      $             WORK( IVT2 ), LDVT2, IWORK( IDXP ), IWORK( IDX ),
00296      $             IWORK( IDXC ), IDXQ, IWORK( COLTYP ), INFO )
00297 *
00298 *     Solve Secular Equation and update singular vectors.
00299 *
00300       LDQ = K
00301       CALL DLASD3( NL, NR, SQRE, K, D, WORK( IQ ), LDQ, WORK( ISIGMA ),
00302      $             U, LDU, WORK( IU2 ), LDU2, VT, LDVT, WORK( IVT2 ),
00303      $             LDVT2, IWORK( IDXC ), IWORK( COLTYP ), WORK( IZ ),
00304      $             INFO )
00305       IF( INFO.NE.0 ) THEN
00306          RETURN
00307       END IF
00308 *
00309 *     Unscale.
00310 *
00311       CALL DLASCL( 'G', 0, 0, ONE, ORGNRM, N, 1, D, N, INFO )
00312 *
00313 *     Prepare the IDXQ sorting permutation.
00314 *
00315       N1 = K
00316       N2 = N - K
00317       CALL DLAMRG( N1, N2, D, 1, -1, IDXQ )
00318 *
00319       RETURN
00320 *
00321 *     End of DLASD1
00322 *
00323       END
 All Files Functions