LAPACK  3.4.1
LAPACK: Linear Algebra PACKage
dlasd7.f
Go to the documentation of this file.
00001 *> \brief \b DLASD7
00002 *
00003 *  =========== DOCUMENTATION ===========
00004 *
00005 * Online html documentation available at 
00006 *            http://www.netlib.org/lapack/explore-html/ 
00007 *
00008 *> \htmlonly
00009 *> Download DLASD7 + dependencies 
00010 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dlasd7.f"> 
00011 *> [TGZ]</a> 
00012 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dlasd7.f"> 
00013 *> [ZIP]</a> 
00014 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dlasd7.f"> 
00015 *> [TXT]</a>
00016 *> \endhtmlonly 
00017 *
00018 *  Definition:
00019 *  ===========
00020 *
00021 *       SUBROUTINE DLASD7( ICOMPQ, NL, NR, SQRE, K, D, Z, ZW, VF, VFW, VL,
00022 *                          VLW, ALPHA, BETA, DSIGMA, IDX, IDXP, IDXQ,
00023 *                          PERM, GIVPTR, GIVCOL, LDGCOL, GIVNUM, LDGNUM,
00024 *                          C, S, 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, * ), IDX( * ), IDXP( * ),
00033 *      $                   IDXQ( * ), PERM( * )
00034 *       DOUBLE PRECISION   D( * ), DSIGMA( * ), GIVNUM( LDGNUM, * ),
00035 *      $                   VF( * ), VFW( * ), VL( * ), VLW( * ), Z( * ),
00036 *      $                   ZW( * )
00037 *       ..
00038 *  
00039 *
00040 *> \par Purpose:
00041 *  =============
00042 *>
00043 *> \verbatim
00044 *>
00045 *> DLASD7 merges the two sets of singular values together into a single
00046 *> sorted set. Then it tries to deflate the size of the problem. There
00047 *> are two ways in which deflation can occur:  when two or more singular
00048 *> values are close together or if there is a tiny entry in the Z
00049 *> vector. For each such occurrence the order of the related
00050 *> secular equation problem is reduced by one.
00051 *>
00052 *> DLASD7 is called from DLASD6.
00053 *> \endverbatim
00054 *
00055 *  Arguments:
00056 *  ==========
00057 *
00058 *> \param[in] ICOMPQ
00059 *> \verbatim
00060 *>          ICOMPQ is INTEGER
00061 *>          Specifies whether singular vectors are to be computed
00062 *>          in compact form, as follows:
00063 *>          = 0: Compute singular values only.
00064 *>          = 1: Compute singular vectors of upper
00065 *>               bidiagonal matrix in compact form.
00066 *> \endverbatim
00067 *>
00068 *> \param[in] NL
00069 *> \verbatim
00070 *>          NL is INTEGER
00071 *>         The row dimension of the upper block. NL >= 1.
00072 *> \endverbatim
00073 *>
00074 *> \param[in] NR
00075 *> \verbatim
00076 *>          NR is INTEGER
00077 *>         The row dimension of the lower block. NR >= 1.
00078 *> \endverbatim
00079 *>
00080 *> \param[in] SQRE
00081 *> \verbatim
00082 *>          SQRE is INTEGER
00083 *>         = 0: the lower block is an NR-by-NR square matrix.
00084 *>         = 1: the lower block is an NR-by-(NR+1) rectangular matrix.
00085 *>
00086 *>         The bidiagonal matrix has
00087 *>         N = NL + NR + 1 rows and
00088 *>         M = N + SQRE >= N columns.
00089 *> \endverbatim
00090 *>
00091 *> \param[out] K
00092 *> \verbatim
00093 *>          K is INTEGER
00094 *>         Contains the dimension of the non-deflated matrix, this is
00095 *>         the order of the related secular equation. 1 <= K <=N.
00096 *> \endverbatim
00097 *>
00098 *> \param[in,out] D
00099 *> \verbatim
00100 *>          D is DOUBLE PRECISION array, dimension ( N )
00101 *>         On entry D contains the singular values of the two submatrices
00102 *>         to be combined. On exit D contains the trailing (N-K) updated
00103 *>         singular values (those which were deflated) sorted into
00104 *>         increasing order.
00105 *> \endverbatim
00106 *>
00107 *> \param[out] Z
00108 *> \verbatim
00109 *>          Z is DOUBLE PRECISION array, dimension ( M )
00110 *>         On exit Z contains the updating row vector in the secular
00111 *>         equation.
00112 *> \endverbatim
00113 *>
00114 *> \param[out] ZW
00115 *> \verbatim
00116 *>          ZW is DOUBLE PRECISION array, dimension ( M )
00117 *>         Workspace for Z.
00118 *> \endverbatim
00119 *>
00120 *> \param[in,out] VF
00121 *> \verbatim
00122 *>          VF is DOUBLE PRECISION array, dimension ( M )
00123 *>         On entry, VF(1:NL+1) contains the first components of all
00124 *>         right singular vectors of the upper block; and VF(NL+2:M)
00125 *>         contains the first components of all right singular vectors
00126 *>         of the lower block. On exit, VF contains the first components
00127 *>         of all right singular vectors of the bidiagonal matrix.
00128 *> \endverbatim
00129 *>
00130 *> \param[out] VFW
00131 *> \verbatim
00132 *>          VFW is DOUBLE PRECISION array, dimension ( M )
00133 *>         Workspace for VF.
00134 *> \endverbatim
00135 *>
00136 *> \param[in,out] VL
00137 *> \verbatim
00138 *>          VL is DOUBLE PRECISION array, dimension ( M )
00139 *>         On entry, VL(1:NL+1) contains the  last components of all
00140 *>         right singular vectors of the upper block; and VL(NL+2:M)
00141 *>         contains the last components of all right singular vectors
00142 *>         of the lower block. On exit, VL contains the last components
00143 *>         of all right singular vectors of the bidiagonal matrix.
00144 *> \endverbatim
00145 *>
00146 *> \param[out] VLW
00147 *> \verbatim
00148 *>          VLW is DOUBLE PRECISION array, dimension ( M )
00149 *>         Workspace for VL.
00150 *> \endverbatim
00151 *>
00152 *> \param[in] ALPHA
00153 *> \verbatim
00154 *>          ALPHA is DOUBLE PRECISION
00155 *>         Contains the diagonal element associated with the added row.
00156 *> \endverbatim
00157 *>
00158 *> \param[in] BETA
00159 *> \verbatim
00160 *>          BETA is DOUBLE PRECISION
00161 *>         Contains the off-diagonal element associated with the added
00162 *>         row.
00163 *> \endverbatim
00164 *>
00165 *> \param[out] DSIGMA
00166 *> \verbatim
00167 *>          DSIGMA is DOUBLE PRECISION array, dimension ( N )
00168 *>         Contains a copy of the diagonal elements (K-1 singular values
00169 *>         and one zero) in the secular equation.
00170 *> \endverbatim
00171 *>
00172 *> \param[out] IDX
00173 *> \verbatim
00174 *>          IDX is INTEGER array, dimension ( N )
00175 *>         This will contain the permutation used to sort the contents of
00176 *>         D into ascending order.
00177 *> \endverbatim
00178 *>
00179 *> \param[out] IDXP
00180 *> \verbatim
00181 *>          IDXP is INTEGER array, dimension ( N )
00182 *>         This will contain the permutation used to place deflated
00183 *>         values of D at the end of the array. On output IDXP(2:K)
00184 *>         points to the nondeflated D-values and IDXP(K+1:N)
00185 *>         points to the deflated singular values.
00186 *> \endverbatim
00187 *>
00188 *> \param[in] IDXQ
00189 *> \verbatim
00190 *>          IDXQ is INTEGER array, dimension ( N )
00191 *>         This contains the permutation which separately sorts the two
00192 *>         sub-problems in D into ascending order.  Note that entries in
00193 *>         the first half of this permutation must first be moved one
00194 *>         position backward; and entries in the second half
00195 *>         must first have NL+1 added to their values.
00196 *> \endverbatim
00197 *>
00198 *> \param[out] PERM
00199 *> \verbatim
00200 *>          PERM is INTEGER array, dimension ( N )
00201 *>         The permutations (from deflation and sorting) to be applied
00202 *>         to each singular block. Not referenced if ICOMPQ = 0.
00203 *> \endverbatim
00204 *>
00205 *> \param[out] GIVPTR
00206 *> \verbatim
00207 *>          GIVPTR is INTEGER
00208 *>         The number of Givens rotations which took place in this
00209 *>         subproblem. Not referenced if ICOMPQ = 0.
00210 *> \endverbatim
00211 *>
00212 *> \param[out] GIVCOL
00213 *> \verbatim
00214 *>          GIVCOL is INTEGER array, dimension ( LDGCOL, 2 )
00215 *>         Each pair of numbers indicates a pair of columns to take place
00216 *>         in a Givens rotation. Not referenced if ICOMPQ = 0.
00217 *> \endverbatim
00218 *>
00219 *> \param[in] LDGCOL
00220 *> \verbatim
00221 *>          LDGCOL is INTEGER
00222 *>         The leading dimension of GIVCOL, must be at least N.
00223 *> \endverbatim
00224 *>
00225 *> \param[out] GIVNUM
00226 *> \verbatim
00227 *>          GIVNUM is DOUBLE PRECISION array, dimension ( LDGNUM, 2 )
00228 *>         Each number indicates the C or S value to be used in the
00229 *>         corresponding Givens rotation. Not referenced if ICOMPQ = 0.
00230 *> \endverbatim
00231 *>
00232 *> \param[in] LDGNUM
00233 *> \verbatim
00234 *>          LDGNUM is INTEGER
00235 *>         The leading dimension of GIVNUM, must be at least N.
00236 *> \endverbatim
00237 *>
00238 *> \param[out] C
00239 *> \verbatim
00240 *>          C is DOUBLE PRECISION
00241 *>         C contains garbage if SQRE =0 and the C-value of a Givens
00242 *>         rotation related to the right null space if SQRE = 1.
00243 *> \endverbatim
00244 *>
00245 *> \param[out] S
00246 *> \verbatim
00247 *>          S is DOUBLE PRECISION
00248 *>         S contains garbage if SQRE =0 and the S-value of a Givens
00249 *>         rotation related to the right null space if SQRE = 1.
00250 *> \endverbatim
00251 *>
00252 *> \param[out] INFO
00253 *> \verbatim
00254 *>          INFO is INTEGER
00255 *>         = 0:  successful exit.
00256 *>         < 0:  if INFO = -i, the i-th argument had an illegal value.
00257 *> \endverbatim
00258 *
00259 *  Authors:
00260 *  ========
00261 *
00262 *> \author Univ. of Tennessee 
00263 *> \author Univ. of California Berkeley 
00264 *> \author Univ. of Colorado Denver 
00265 *> \author NAG Ltd. 
00266 *
00267 *> \date November 2011
00268 *
00269 *> \ingroup auxOTHERauxiliary
00270 *
00271 *> \par Contributors:
00272 *  ==================
00273 *>
00274 *>     Ming Gu and Huan Ren, Computer Science Division, University of
00275 *>     California at Berkeley, USA
00276 *>
00277 *  =====================================================================
00278       SUBROUTINE DLASD7( ICOMPQ, NL, NR, SQRE, K, D, Z, ZW, VF, VFW, VL,
00279      $                   VLW, ALPHA, BETA, DSIGMA, IDX, IDXP, IDXQ,
00280      $                   PERM, GIVPTR, GIVCOL, LDGCOL, GIVNUM, LDGNUM,
00281      $                   C, S, INFO )
00282 *
00283 *  -- LAPACK auxiliary routine (version 3.4.0) --
00284 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
00285 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
00286 *     November 2011
00287 *
00288 *     .. Scalar Arguments ..
00289       INTEGER            GIVPTR, ICOMPQ, INFO, K, LDGCOL, LDGNUM, NL,
00290      $                   NR, SQRE
00291       DOUBLE PRECISION   ALPHA, BETA, C, S
00292 *     ..
00293 *     .. Array Arguments ..
00294       INTEGER            GIVCOL( LDGCOL, * ), IDX( * ), IDXP( * ),
00295      $                   IDXQ( * ), PERM( * )
00296       DOUBLE PRECISION   D( * ), DSIGMA( * ), GIVNUM( LDGNUM, * ),
00297      $                   VF( * ), VFW( * ), VL( * ), VLW( * ), Z( * ),
00298      $                   ZW( * )
00299 *     ..
00300 *
00301 *  =====================================================================
00302 *
00303 *     .. Parameters ..
00304       DOUBLE PRECISION   ZERO, ONE, TWO, EIGHT
00305       PARAMETER          ( ZERO = 0.0D+0, ONE = 1.0D+0, TWO = 2.0D+0,
00306      $                   EIGHT = 8.0D+0 )
00307 *     ..
00308 *     .. Local Scalars ..
00309 *
00310       INTEGER            I, IDXI, IDXJ, IDXJP, J, JP, JPREV, K2, M, N,
00311      $                   NLP1, NLP2
00312       DOUBLE PRECISION   EPS, HLFTOL, TAU, TOL, Z1
00313 *     ..
00314 *     .. External Subroutines ..
00315       EXTERNAL           DCOPY, DLAMRG, DROT, XERBLA
00316 *     ..
00317 *     .. External Functions ..
00318       DOUBLE PRECISION   DLAMCH, DLAPY2
00319       EXTERNAL           DLAMCH, DLAPY2
00320 *     ..
00321 *     .. Intrinsic Functions ..
00322       INTRINSIC          ABS, MAX
00323 *     ..
00324 *     .. Executable Statements ..
00325 *
00326 *     Test the input parameters.
00327 *
00328       INFO = 0
00329       N = NL + NR + 1
00330       M = N + SQRE
00331 *
00332       IF( ( ICOMPQ.LT.0 ) .OR. ( ICOMPQ.GT.1 ) ) THEN
00333          INFO = -1
00334       ELSE IF( NL.LT.1 ) THEN
00335          INFO = -2
00336       ELSE IF( NR.LT.1 ) THEN
00337          INFO = -3
00338       ELSE IF( ( SQRE.LT.0 ) .OR. ( SQRE.GT.1 ) ) THEN
00339          INFO = -4
00340       ELSE IF( LDGCOL.LT.N ) THEN
00341          INFO = -22
00342       ELSE IF( LDGNUM.LT.N ) THEN
00343          INFO = -24
00344       END IF
00345       IF( INFO.NE.0 ) THEN
00346          CALL XERBLA( 'DLASD7', -INFO )
00347          RETURN
00348       END IF
00349 *
00350       NLP1 = NL + 1
00351       NLP2 = NL + 2
00352       IF( ICOMPQ.EQ.1 ) THEN
00353          GIVPTR = 0
00354       END IF
00355 *
00356 *     Generate the first part of the vector Z and move the singular
00357 *     values in the first part of D one position backward.
00358 *
00359       Z1 = ALPHA*VL( NLP1 )
00360       VL( NLP1 ) = ZERO
00361       TAU = VF( NLP1 )
00362       DO 10 I = NL, 1, -1
00363          Z( I+1 ) = ALPHA*VL( I )
00364          VL( I ) = ZERO
00365          VF( I+1 ) = VF( I )
00366          D( I+1 ) = D( I )
00367          IDXQ( I+1 ) = IDXQ( I ) + 1
00368    10 CONTINUE
00369       VF( 1 ) = TAU
00370 *
00371 *     Generate the second part of the vector Z.
00372 *
00373       DO 20 I = NLP2, M
00374          Z( I ) = BETA*VF( I )
00375          VF( I ) = ZERO
00376    20 CONTINUE
00377 *
00378 *     Sort the singular values into increasing order
00379 *
00380       DO 30 I = NLP2, N
00381          IDXQ( I ) = IDXQ( I ) + NLP1
00382    30 CONTINUE
00383 *
00384 *     DSIGMA, IDXC, IDXC, and ZW are used as storage space.
00385 *
00386       DO 40 I = 2, N
00387          DSIGMA( I ) = D( IDXQ( I ) )
00388          ZW( I ) = Z( IDXQ( I ) )
00389          VFW( I ) = VF( IDXQ( I ) )
00390          VLW( I ) = VL( IDXQ( I ) )
00391    40 CONTINUE
00392 *
00393       CALL DLAMRG( NL, NR, DSIGMA( 2 ), 1, 1, IDX( 2 ) )
00394 *
00395       DO 50 I = 2, N
00396          IDXI = 1 + IDX( I )
00397          D( I ) = DSIGMA( IDXI )
00398          Z( I ) = ZW( IDXI )
00399          VF( I ) = VFW( IDXI )
00400          VL( I ) = VLW( IDXI )
00401    50 CONTINUE
00402 *
00403 *     Calculate the allowable deflation tolerence
00404 *
00405       EPS = DLAMCH( 'Epsilon' )
00406       TOL = MAX( ABS( ALPHA ), ABS( BETA ) )
00407       TOL = EIGHT*EIGHT*EPS*MAX( ABS( D( N ) ), TOL )
00408 *
00409 *     There are 2 kinds of deflation -- first a value in the z-vector
00410 *     is small, second two (or more) singular values are very close
00411 *     together (their difference is small).
00412 *
00413 *     If the value in the z-vector is small, we simply permute the
00414 *     array so that the corresponding singular value is moved to the
00415 *     end.
00416 *
00417 *     If two values in the D-vector are close, we perform a two-sided
00418 *     rotation designed to make one of the corresponding z-vector
00419 *     entries zero, and then permute the array so that the deflated
00420 *     singular value is moved to the end.
00421 *
00422 *     If there are multiple singular values then the problem deflates.
00423 *     Here the number of equal singular values are found.  As each equal
00424 *     singular value is found, an elementary reflector is computed to
00425 *     rotate the corresponding singular subspace so that the
00426 *     corresponding components of Z are zero in this new basis.
00427 *
00428       K = 1
00429       K2 = N + 1
00430       DO 60 J = 2, N
00431          IF( ABS( Z( J ) ).LE.TOL ) THEN
00432 *
00433 *           Deflate due to small z component.
00434 *
00435             K2 = K2 - 1
00436             IDXP( K2 ) = J
00437             IF( J.EQ.N )
00438      $         GO TO 100
00439          ELSE
00440             JPREV = J
00441             GO TO 70
00442          END IF
00443    60 CONTINUE
00444    70 CONTINUE
00445       J = JPREV
00446    80 CONTINUE
00447       J = J + 1
00448       IF( J.GT.N )
00449      $   GO TO 90
00450       IF( ABS( Z( J ) ).LE.TOL ) THEN
00451 *
00452 *        Deflate due to small z component.
00453 *
00454          K2 = K2 - 1
00455          IDXP( K2 ) = J
00456       ELSE
00457 *
00458 *        Check if singular values are close enough to allow deflation.
00459 *
00460          IF( ABS( D( J )-D( JPREV ) ).LE.TOL ) THEN
00461 *
00462 *           Deflation is possible.
00463 *
00464             S = Z( JPREV )
00465             C = Z( J )
00466 *
00467 *           Find sqrt(a**2+b**2) without overflow or
00468 *           destructive underflow.
00469 *
00470             TAU = DLAPY2( C, S )
00471             Z( J ) = TAU
00472             Z( JPREV ) = ZERO
00473             C = C / TAU
00474             S = -S / TAU
00475 *
00476 *           Record the appropriate Givens rotation
00477 *
00478             IF( ICOMPQ.EQ.1 ) THEN
00479                GIVPTR = GIVPTR + 1
00480                IDXJP = IDXQ( IDX( JPREV )+1 )
00481                IDXJ = IDXQ( IDX( J )+1 )
00482                IF( IDXJP.LE.NLP1 ) THEN
00483                   IDXJP = IDXJP - 1
00484                END IF
00485                IF( IDXJ.LE.NLP1 ) THEN
00486                   IDXJ = IDXJ - 1
00487                END IF
00488                GIVCOL( GIVPTR, 2 ) = IDXJP
00489                GIVCOL( GIVPTR, 1 ) = IDXJ
00490                GIVNUM( GIVPTR, 2 ) = C
00491                GIVNUM( GIVPTR, 1 ) = S
00492             END IF
00493             CALL DROT( 1, VF( JPREV ), 1, VF( J ), 1, C, S )
00494             CALL DROT( 1, VL( JPREV ), 1, VL( J ), 1, C, S )
00495             K2 = K2 - 1
00496             IDXP( K2 ) = JPREV
00497             JPREV = J
00498          ELSE
00499             K = K + 1
00500             ZW( K ) = Z( JPREV )
00501             DSIGMA( K ) = D( JPREV )
00502             IDXP( K ) = JPREV
00503             JPREV = J
00504          END IF
00505       END IF
00506       GO TO 80
00507    90 CONTINUE
00508 *
00509 *     Record the last singular value.
00510 *
00511       K = K + 1
00512       ZW( K ) = Z( JPREV )
00513       DSIGMA( K ) = D( JPREV )
00514       IDXP( K ) = JPREV
00515 *
00516   100 CONTINUE
00517 *
00518 *     Sort the singular values into DSIGMA. The singular values which
00519 *     were not deflated go into the first K slots of DSIGMA, except
00520 *     that DSIGMA(1) is treated separately.
00521 *
00522       DO 110 J = 2, N
00523          JP = IDXP( J )
00524          DSIGMA( J ) = D( JP )
00525          VFW( J ) = VF( JP )
00526          VLW( J ) = VL( JP )
00527   110 CONTINUE
00528       IF( ICOMPQ.EQ.1 ) THEN
00529          DO 120 J = 2, N
00530             JP = IDXP( J )
00531             PERM( J ) = IDXQ( IDX( JP )+1 )
00532             IF( PERM( J ).LE.NLP1 ) THEN
00533                PERM( J ) = PERM( J ) - 1
00534             END IF
00535   120    CONTINUE
00536       END IF
00537 *
00538 *     The deflated singular values go back into the last N - K slots of
00539 *     D.
00540 *
00541       CALL DCOPY( N-K, DSIGMA( K+1 ), 1, D( K+1 ), 1 )
00542 *
00543 *     Determine DSIGMA(1), DSIGMA(2), Z(1), VF(1), VL(1), VF(M), and
00544 *     VL(M).
00545 *
00546       DSIGMA( 1 ) = ZERO
00547       HLFTOL = TOL / TWO
00548       IF( ABS( DSIGMA( 2 ) ).LE.HLFTOL )
00549      $   DSIGMA( 2 ) = HLFTOL
00550       IF( M.GT.N ) THEN
00551          Z( 1 ) = DLAPY2( Z1, Z( M ) )
00552          IF( Z( 1 ).LE.TOL ) THEN
00553             C = ONE
00554             S = ZERO
00555             Z( 1 ) = TOL
00556          ELSE
00557             C = Z1 / Z( 1 )
00558             S = -Z( M ) / Z( 1 )
00559          END IF
00560          CALL DROT( 1, VF( M ), 1, VF( 1 ), 1, C, S )
00561          CALL DROT( 1, VL( M ), 1, VL( 1 ), 1, C, S )
00562       ELSE
00563          IF( ABS( Z1 ).LE.TOL ) THEN
00564             Z( 1 ) = TOL
00565          ELSE
00566             Z( 1 ) = Z1
00567          END IF
00568       END IF
00569 *
00570 *     Restore Z, VF, and VL.
00571 *
00572       CALL DCOPY( K-1, ZW( 2 ), 1, Z( 2 ), 1 )
00573       CALL DCOPY( N-1, VFW( 2 ), 1, VF( 2 ), 1 )
00574       CALL DCOPY( N-1, VLW( 2 ), 1, VL( 2 ), 1 )
00575 *
00576       RETURN
00577 *
00578 *     End of DLASD7
00579 *
00580       END
 All Files Functions