LAPACK  3.4.1
LAPACK: Linear Algebra PACKage
dlasd6.f
Go to the documentation of this file.
00001 *> \brief \b DLASD6
00002 *
00003 *  =========== DOCUMENTATION ===========
00004 *
00005 * Online html documentation available at 
00006 *            http://www.netlib.org/lapack/explore-html/ 
00007 *
00008 *> \htmlonly
00009 *> Download DLASD6 + dependencies 
00010 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dlasd6.f"> 
00011 *> [TGZ]</a> 
00012 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dlasd6.f"> 
00013 *> [ZIP]</a> 
00014 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dlasd6.f"> 
00015 *> [TXT]</a>
00016 *> \endhtmlonly 
00017 *
00018 *  Definition:
00019 *  ===========
00020 *
00021 *       SUBROUTINE DLASD6( ICOMPQ, NL, NR, SQRE, D, VF, VL, ALPHA, BETA,
00022 *                          IDXQ, PERM, GIVPTR, GIVCOL, LDGCOL, GIVNUM,
00023 *                          LDGNUM, POLES, DIFL, DIFR, Z, K, C, S, WORK,
00024 *                          IWORK, INFO )
00025 * 
00026 *       .. Scalar Arguments ..
00027 *       INTEGER            GIVPTR, ICOMPQ, INFO, K, LDGCOL, LDGNUM, NL,
00028 *      $                   NR, SQRE
00029 *       DOUBLE PRECISION   ALPHA, BETA, C, S
00030 *       ..
00031 *       .. Array Arguments ..
00032 *       INTEGER            GIVCOL( LDGCOL, * ), IDXQ( * ), IWORK( * ),
00033 *      $                   PERM( * )
00034 *       DOUBLE PRECISION   D( * ), DIFL( * ), DIFR( * ),
00035 *      $                   GIVNUM( LDGNUM, * ), POLES( LDGNUM, * ),
00036 *      $                   VF( * ), VL( * ), WORK( * ), Z( * )
00037 *       ..
00038 *  
00039 *
00040 *> \par Purpose:
00041 *  =============
00042 *>
00043 *> \verbatim
00044 *>
00045 *> DLASD6 computes the SVD of an updated upper bidiagonal matrix B
00046 *> obtained by merging two smaller ones by appending a row. This
00047 *> routine is used only for the problem which requires all singular
00048 *> values and optionally singular vector matrices in factored form.
00049 *> B is an N-by-M matrix with N = NL + NR + 1 and M = N + SQRE.
00050 *> A related subroutine, DLASD1, handles the case in which all singular
00051 *> values and singular vectors of the bidiagonal matrix are desired.
00052 *>
00053 *> DLASD6 computes the SVD as follows:
00054 *>
00055 *>               ( D1(in)    0    0       0 )
00056 *>   B = U(in) * (   Z1**T   a   Z2**T    b ) * VT(in)
00057 *>               (   0       0   D2(in)   0 )
00058 *>
00059 *>     = U(out) * ( D(out) 0) * VT(out)
00060 *>
00061 *> where Z**T = (Z1**T a Z2**T b) = u**T VT**T, and u is a vector of dimension M
00062 *> with ALPHA and BETA in the NL+1 and NL+2 th entries and zeros
00063 *> elsewhere; and the entry b is empty if SQRE = 0.
00064 *>
00065 *> The singular values of B can be computed using D1, D2, the first
00066 *> components of all the right singular vectors of the lower block, and
00067 *> the last components of all the right singular vectors of the upper
00068 *> block. These components are stored and updated in VF and VL,
00069 *> respectively, in DLASD6. Hence U and VT are not explicitly
00070 *> referenced.
00071 *>
00072 *> The singular values are stored in D. The algorithm consists of two
00073 *> stages:
00074 *>
00075 *>       The first stage consists of deflating the size of the problem
00076 *>       when there are multiple singular values or if there is a zero
00077 *>       in the Z vector. For each such occurence the dimension of the
00078 *>       secular equation problem is reduced by one. This stage is
00079 *>       performed by the routine DLASD7.
00080 *>
00081 *>       The second stage consists of calculating the updated
00082 *>       singular values. This is done by finding the roots of the
00083 *>       secular equation via the routine DLASD4 (as called by DLASD8).
00084 *>       This routine also updates VF and VL and computes the distances
00085 *>       between the updated singular values and the old singular
00086 *>       values.
00087 *>
00088 *> DLASD6 is called from DLASDA.
00089 *> \endverbatim
00090 *
00091 *  Arguments:
00092 *  ==========
00093 *
00094 *> \param[in] ICOMPQ
00095 *> \verbatim
00096 *>          ICOMPQ is INTEGER
00097 *>         Specifies whether singular vectors are to be computed in
00098 *>         factored form:
00099 *>         = 0: Compute singular values only.
00100 *>         = 1: Compute singular vectors in factored form as well.
00101 *> \endverbatim
00102 *>
00103 *> \param[in] NL
00104 *> \verbatim
00105 *>          NL is INTEGER
00106 *>         The row dimension of the upper block.  NL >= 1.
00107 *> \endverbatim
00108 *>
00109 *> \param[in] NR
00110 *> \verbatim
00111 *>          NR is INTEGER
00112 *>         The row dimension of the lower block.  NR >= 1.
00113 *> \endverbatim
00114 *>
00115 *> \param[in] SQRE
00116 *> \verbatim
00117 *>          SQRE is INTEGER
00118 *>         = 0: the lower block is an NR-by-NR square matrix.
00119 *>         = 1: the lower block is an NR-by-(NR+1) rectangular matrix.
00120 *>
00121 *>         The bidiagonal matrix has row dimension N = NL + NR + 1,
00122 *>         and column dimension M = N + SQRE.
00123 *> \endverbatim
00124 *>
00125 *> \param[in,out] D
00126 *> \verbatim
00127 *>          D is DOUBLE PRECISION array, dimension ( NL+NR+1 ).
00128 *>         On entry D(1:NL,1:NL) contains the singular values of the
00129 *>         upper block, and D(NL+2:N) contains the singular values
00130 *>         of the lower block. On exit D(1:N) contains the singular
00131 *>         values of the modified matrix.
00132 *> \endverbatim
00133 *>
00134 *> \param[in,out] VF
00135 *> \verbatim
00136 *>          VF is DOUBLE PRECISION array, dimension ( M )
00137 *>         On entry, VF(1:NL+1) contains the first components of all
00138 *>         right singular vectors of the upper block; and VF(NL+2:M)
00139 *>         contains the first components of all right singular vectors
00140 *>         of the lower block. On exit, VF contains the first components
00141 *>         of all right singular vectors of the bidiagonal matrix.
00142 *> \endverbatim
00143 *>
00144 *> \param[in,out] VL
00145 *> \verbatim
00146 *>          VL is DOUBLE PRECISION array, dimension ( M )
00147 *>         On entry, VL(1:NL+1) contains the  last components of all
00148 *>         right singular vectors of the upper block; and VL(NL+2:M)
00149 *>         contains the last components of all right singular vectors of
00150 *>         the lower block. On exit, VL contains the last components of
00151 *>         all right singular vectors of the bidiagonal matrix.
00152 *> \endverbatim
00153 *>
00154 *> \param[in,out] ALPHA
00155 *> \verbatim
00156 *>          ALPHA is DOUBLE PRECISION
00157 *>         Contains the diagonal element associated with the added row.
00158 *> \endverbatim
00159 *>
00160 *> \param[in,out] BETA
00161 *> \verbatim
00162 *>          BETA is DOUBLE PRECISION
00163 *>         Contains the off-diagonal element associated with the added
00164 *>         row.
00165 *> \endverbatim
00166 *>
00167 *> \param[out] IDXQ
00168 *> \verbatim
00169 *>          IDXQ is INTEGER array, dimension ( N )
00170 *>         This contains the permutation which will reintegrate the
00171 *>         subproblem just solved back into sorted order, i.e.
00172 *>         D( IDXQ( I = 1, N ) ) will be in ascending order.
00173 *> \endverbatim
00174 *>
00175 *> \param[out] PERM
00176 *> \verbatim
00177 *>          PERM is INTEGER array, dimension ( N )
00178 *>         The permutations (from deflation and sorting) to be applied
00179 *>         to each block. Not referenced if ICOMPQ = 0.
00180 *> \endverbatim
00181 *>
00182 *> \param[out] GIVPTR
00183 *> \verbatim
00184 *>          GIVPTR is INTEGER
00185 *>         The number of Givens rotations which took place in this
00186 *>         subproblem. Not referenced if ICOMPQ = 0.
00187 *> \endverbatim
00188 *>
00189 *> \param[out] GIVCOL
00190 *> \verbatim
00191 *>          GIVCOL is INTEGER array, dimension ( LDGCOL, 2 )
00192 *>         Each pair of numbers indicates a pair of columns to take place
00193 *>         in a Givens rotation. Not referenced if ICOMPQ = 0.
00194 *> \endverbatim
00195 *>
00196 *> \param[in] LDGCOL
00197 *> \verbatim
00198 *>          LDGCOL is INTEGER
00199 *>         leading dimension of GIVCOL, must be at least N.
00200 *> \endverbatim
00201 *>
00202 *> \param[out] GIVNUM
00203 *> \verbatim
00204 *>          GIVNUM is DOUBLE PRECISION array, dimension ( LDGNUM, 2 )
00205 *>         Each number indicates the C or S value to be used in the
00206 *>         corresponding Givens rotation. Not referenced if ICOMPQ = 0.
00207 *> \endverbatim
00208 *>
00209 *> \param[in] LDGNUM
00210 *> \verbatim
00211 *>          LDGNUM is INTEGER
00212 *>         The leading dimension of GIVNUM and POLES, must be at least N.
00213 *> \endverbatim
00214 *>
00215 *> \param[out] POLES
00216 *> \verbatim
00217 *>          POLES is DOUBLE PRECISION array, dimension ( LDGNUM, 2 )
00218 *>         On exit, POLES(1,*) is an array containing the new singular
00219 *>         values obtained from solving the secular equation, and
00220 *>         POLES(2,*) is an array containing the poles in the secular
00221 *>         equation. Not referenced if ICOMPQ = 0.
00222 *> \endverbatim
00223 *>
00224 *> \param[out] DIFL
00225 *> \verbatim
00226 *>          DIFL is DOUBLE PRECISION array, dimension ( N )
00227 *>         On exit, DIFL(I) is the distance between I-th updated
00228 *>         (undeflated) singular value and the I-th (undeflated) old
00229 *>         singular value.
00230 *> \endverbatim
00231 *>
00232 *> \param[out] DIFR
00233 *> \verbatim
00234 *>          DIFR is DOUBLE PRECISION array,
00235 *>                  dimension ( LDGNUM, 2 ) if ICOMPQ = 1 and
00236 *>                  dimension ( N ) if ICOMPQ = 0.
00237 *>         On exit, DIFR(I, 1) is the distance between I-th updated
00238 *>         (undeflated) singular value and the I+1-th (undeflated) old
00239 *>         singular value.
00240 *>
00241 *>         If ICOMPQ = 1, DIFR(1:K,2) is an array containing the
00242 *>         normalizing factors for the right singular vector matrix.
00243 *>
00244 *>         See DLASD8 for details on DIFL and DIFR.
00245 *> \endverbatim
00246 *>
00247 *> \param[out] Z
00248 *> \verbatim
00249 *>          Z is DOUBLE PRECISION array, dimension ( M )
00250 *>         The first elements of this array contain the components
00251 *>         of the deflation-adjusted updating row vector.
00252 *> \endverbatim
00253 *>
00254 *> \param[out] K
00255 *> \verbatim
00256 *>          K is INTEGER
00257 *>         Contains the dimension of the non-deflated matrix,
00258 *>         This is the order of the related secular equation. 1 <= K <=N.
00259 *> \endverbatim
00260 *>
00261 *> \param[out] C
00262 *> \verbatim
00263 *>          C is DOUBLE PRECISION
00264 *>         C contains garbage if SQRE =0 and the C-value of a Givens
00265 *>         rotation related to the right null space if SQRE = 1.
00266 *> \endverbatim
00267 *>
00268 *> \param[out] S
00269 *> \verbatim
00270 *>          S is DOUBLE PRECISION
00271 *>         S contains garbage if SQRE =0 and the S-value of a Givens
00272 *>         rotation related to the right null space if SQRE = 1.
00273 *> \endverbatim
00274 *>
00275 *> \param[out] WORK
00276 *> \verbatim
00277 *>          WORK is DOUBLE PRECISION array, dimension ( 4 * M )
00278 *> \endverbatim
00279 *>
00280 *> \param[out] IWORK
00281 *> \verbatim
00282 *>          IWORK is INTEGER array, dimension ( 3 * N )
00283 *> \endverbatim
00284 *>
00285 *> \param[out] INFO
00286 *> \verbatim
00287 *>          INFO is INTEGER
00288 *>          = 0:  successful exit.
00289 *>          < 0:  if INFO = -i, the i-th argument had an illegal value.
00290 *>          > 0:  if INFO = 1, a singular value did not converge
00291 *> \endverbatim
00292 *
00293 *  Authors:
00294 *  ========
00295 *
00296 *> \author Univ. of Tennessee 
00297 *> \author Univ. of California Berkeley 
00298 *> \author Univ. of Colorado Denver 
00299 *> \author NAG Ltd. 
00300 *
00301 *> \date November 2011
00302 *
00303 *> \ingroup auxOTHERauxiliary
00304 *
00305 *> \par Contributors:
00306 *  ==================
00307 *>
00308 *>     Ming Gu and Huan Ren, Computer Science Division, University of
00309 *>     California at Berkeley, USA
00310 *>
00311 *  =====================================================================
00312       SUBROUTINE DLASD6( ICOMPQ, NL, NR, SQRE, D, VF, VL, ALPHA, BETA,
00313      $                   IDXQ, PERM, GIVPTR, GIVCOL, LDGCOL, GIVNUM,
00314      $                   LDGNUM, POLES, DIFL, DIFR, Z, K, C, S, WORK,
00315      $                   IWORK, INFO )
00316 *
00317 *  -- LAPACK auxiliary routine (version 3.4.0) --
00318 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
00319 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
00320 *     November 2011
00321 *
00322 *     .. Scalar Arguments ..
00323       INTEGER            GIVPTR, ICOMPQ, INFO, K, LDGCOL, LDGNUM, NL,
00324      $                   NR, SQRE
00325       DOUBLE PRECISION   ALPHA, BETA, C, S
00326 *     ..
00327 *     .. Array Arguments ..
00328       INTEGER            GIVCOL( LDGCOL, * ), IDXQ( * ), IWORK( * ),
00329      $                   PERM( * )
00330       DOUBLE PRECISION   D( * ), DIFL( * ), DIFR( * ),
00331      $                   GIVNUM( LDGNUM, * ), POLES( LDGNUM, * ),
00332      $                   VF( * ), VL( * ), WORK( * ), Z( * )
00333 *     ..
00334 *
00335 *  =====================================================================
00336 *
00337 *     .. Parameters ..
00338       DOUBLE PRECISION   ONE, ZERO
00339       PARAMETER          ( ONE = 1.0D+0, ZERO = 0.0D+0 )
00340 *     ..
00341 *     .. Local Scalars ..
00342       INTEGER            I, IDX, IDXC, IDXP, ISIGMA, IVFW, IVLW, IW, M,
00343      $                   N, N1, N2
00344       DOUBLE PRECISION   ORGNRM
00345 *     ..
00346 *     .. External Subroutines ..
00347       EXTERNAL           DCOPY, DLAMRG, DLASCL, DLASD7, DLASD8, XERBLA
00348 *     ..
00349 *     .. Intrinsic Functions ..
00350       INTRINSIC          ABS, MAX
00351 *     ..
00352 *     .. Executable Statements ..
00353 *
00354 *     Test the input parameters.
00355 *
00356       INFO = 0
00357       N = NL + NR + 1
00358       M = N + SQRE
00359 *
00360       IF( ( ICOMPQ.LT.0 ) .OR. ( ICOMPQ.GT.1 ) ) THEN
00361          INFO = -1
00362       ELSE IF( NL.LT.1 ) THEN
00363          INFO = -2
00364       ELSE IF( NR.LT.1 ) THEN
00365          INFO = -3
00366       ELSE IF( ( SQRE.LT.0 ) .OR. ( SQRE.GT.1 ) ) THEN
00367          INFO = -4
00368       ELSE IF( LDGCOL.LT.N ) THEN
00369          INFO = -14
00370       ELSE IF( LDGNUM.LT.N ) THEN
00371          INFO = -16
00372       END IF
00373       IF( INFO.NE.0 ) THEN
00374          CALL XERBLA( 'DLASD6', -INFO )
00375          RETURN
00376       END IF
00377 *
00378 *     The following values are for bookkeeping purposes only.  They are
00379 *     integer pointers which indicate the portion of the workspace
00380 *     used by a particular array in DLASD7 and DLASD8.
00381 *
00382       ISIGMA = 1
00383       IW = ISIGMA + N
00384       IVFW = IW + M
00385       IVLW = IVFW + M
00386 *
00387       IDX = 1
00388       IDXC = IDX + N
00389       IDXP = IDXC + N
00390 *
00391 *     Scale.
00392 *
00393       ORGNRM = MAX( ABS( ALPHA ), ABS( BETA ) )
00394       D( NL+1 ) = ZERO
00395       DO 10 I = 1, N
00396          IF( ABS( D( I ) ).GT.ORGNRM ) THEN
00397             ORGNRM = ABS( D( I ) )
00398          END IF
00399    10 CONTINUE
00400       CALL DLASCL( 'G', 0, 0, ORGNRM, ONE, N, 1, D, N, INFO )
00401       ALPHA = ALPHA / ORGNRM
00402       BETA = BETA / ORGNRM
00403 *
00404 *     Sort and Deflate singular values.
00405 *
00406       CALL DLASD7( ICOMPQ, NL, NR, SQRE, K, D, Z, WORK( IW ), VF,
00407      $             WORK( IVFW ), VL, WORK( IVLW ), ALPHA, BETA,
00408      $             WORK( ISIGMA ), IWORK( IDX ), IWORK( IDXP ), IDXQ,
00409      $             PERM, GIVPTR, GIVCOL, LDGCOL, GIVNUM, LDGNUM, C, S,
00410      $             INFO )
00411 *
00412 *     Solve Secular Equation, compute DIFL, DIFR, and update VF, VL.
00413 *
00414       CALL DLASD8( ICOMPQ, K, D, Z, VF, VL, DIFL, DIFR, LDGNUM,
00415      $             WORK( ISIGMA ), WORK( IW ), INFO )
00416 *
00417 *     Handle error returned
00418 *
00419       IF( INFO.NE.0 ) THEN
00420          CALL XERBLA( 'DLASD8', -INFO )
00421          RETURN
00422       END IF
00423 *
00424 *     Save the poles if ICOMPQ = 1.
00425 *
00426       IF( ICOMPQ.EQ.1 ) THEN
00427          CALL DCOPY( K, D, 1, POLES( 1, 1 ), 1 )
00428          CALL DCOPY( K, WORK( ISIGMA ), 1, POLES( 1, 2 ), 1 )
00429       END IF
00430 *
00431 *     Unscale.
00432 *
00433       CALL DLASCL( 'G', 0, 0, ONE, ORGNRM, N, 1, D, N, INFO )
00434 *
00435 *     Prepare the IDXQ sorting permutation.
00436 *
00437       N1 = K
00438       N2 = N - K
00439       CALL DLAMRG( N1, N2, D, 1, -1, IDXQ )
00440 *
00441       RETURN
00442 *
00443 *     End of DLASD6
00444 *
00445       END
 All Files Functions