SRC/dlasd7.f

   1 *> \brief \b DLASD7 merges the two sets of singular values together into a single sorted set. Then it tries to deflate the size of the problem. Used by sbdsdc.
   2 *
   3 *  =========== DOCUMENTATION ===========
   4 *
   5 * Online html documentation available at
   6 *            http://www.netlib.org/lapack/explore-html/
   7 *
   8 *> \htmlonly
   9 *> Download DLASD7 + dependencies
  10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dlasd7.f">
  11 *> [TGZ]</a>
  12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dlasd7.f">
  13 *> [ZIP]</a>
  14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dlasd7.f">
  15 *> [TXT]</a>
  16 *> \endhtmlonly
  17 *
  18 *  Definition:
  19 *  ===========
  20 *
  21 *       SUBROUTINE DLASD7( ICOMPQ, NL, NR, SQRE, K, D, Z, ZW, VF, VFW, VL,
  22 *                          VLW, ALPHA, BETA, DSIGMA, IDX, IDXP, IDXQ,
  23 *                          PERM, GIVPTR, GIVCOL, LDGCOL, GIVNUM, LDGNUM,
  24 *                          C, S, INFO )
  25 *
  26 *       .. Scalar Arguments ..
  27 *       INTEGER            GIVPTR, ICOMPQ, INFO, K, LDGCOL, LDGNUM, NL,
  28 *      $                   NR, SQRE
  29 *       DOUBLE PRECISION   ALPHA, BETA, C, S
  30 *       ..
  31 *       .. Array Arguments ..
  32 *       INTEGER            GIVCOL( LDGCOL, * ), IDX( * ), IDXP( * ),
  33 *      $                   IDXQ( * ), PERM( * )
  34 *       DOUBLE PRECISION   D( * ), DSIGMA( * ), GIVNUM( LDGNUM, * ),
  35 *      $                   VF( * ), VFW( * ), VL( * ), VLW( * ), Z( * ),
  36 *      $                   ZW( * )
  37 *       ..
  38 *
  39 *
  40 *> \par Purpose:
  41 *  =============
  42 *>
  43 *> \verbatim
  44 *>
  45 *> DLASD7 merges the two sets of singular values together into a single
  46 *> sorted set. Then it tries to deflate the size of the problem. There
  47 *> are two ways in which deflation can occur:  when two or more singular
  48 *> values are close together or if there is a tiny entry in the Z
  49 *> vector. For each such occurrence the order of the related
  50 *> secular equation problem is reduced by one.
  51 *>
  52 *> DLASD7 is called from DLASD6.
  53 *> \endverbatim
  54 *
  55 *  Arguments:
  56 *  ==========
  57 *
  58 *> \param[in] ICOMPQ
  59 *> \verbatim
  60 *>          ICOMPQ is INTEGER
  61 *>          Specifies whether singular vectors are to be computed
  62 *>          in compact form, as follows:
  63 *>          = 0: Compute singular values only.
  64 *>          = 1: Compute singular vectors of upper
  65 *>               bidiagonal matrix in compact form.
  66 *> \endverbatim
  67 *>
  68 *> \param[in] NL
  69 *> \verbatim
  70 *>          NL is INTEGER
  71 *>         The row dimension of the upper block. NL >= 1.
  72 *> \endverbatim
  73 *>
  74 *> \param[in] NR
  75 *> \verbatim
  76 *>          NR is INTEGER
  77 *>         The row dimension of the lower block. NR >= 1.
  78 *> \endverbatim
  79 *>
  80 *> \param[in] SQRE
  81 *> \verbatim
  82 *>          SQRE is INTEGER
  83 *>         = 0: the lower block is an NR-by-NR square matrix.
  84 *>         = 1: the lower block is an NR-by-(NR+1) rectangular matrix.
  85 *>
  86 *>         The bidiagonal matrix has
  87 *>         N = NL + NR + 1 rows and
  88 *>         M = N + SQRE >= N columns.
  89 *> \endverbatim
  90 *>
  91 *> \param[out] K
  92 *> \verbatim
  93 *>          K is INTEGER
  94 *>         Contains the dimension of the non-deflated matrix, this is
  95 *>         the order of the related secular equation. 1 <= K <=N.
  96 *> \endverbatim
  97 *>
  98 *> \param[in,out] D
  99 *> \verbatim
 100 *>          D is DOUBLE PRECISION array, dimension ( N )
 101 *>         On entry D contains the singular values of the two submatrices
 102 *>         to be combined. On exit D contains the trailing (N-K) updated
 103 *>         singular values (those which were deflated) sorted into
 104 *>         increasing order.
 105 *> \endverbatim
 106 *>
 107 *> \param[out] Z
 108 *> \verbatim
 109 *>          Z is DOUBLE PRECISION array, dimension ( M )
 110 *>         On exit Z contains the updating row vector in the secular
 111 *>         equation.
 112 *> \endverbatim
 113 *>
 114 *> \param[out] ZW
 115 *> \verbatim
 116 *>          ZW is DOUBLE PRECISION array, dimension ( M )
 117 *>         Workspace for Z.
 118 *> \endverbatim
 119 *>
 120 *> \param[in,out] VF
 121 *> \verbatim
 122 *>          VF is DOUBLE PRECISION array, dimension ( M )
 123 *>         On entry, VF(1:NL+1) contains the first components of all
 124 *>         right singular vectors of the upper block; and VF(NL+2:M)
 125 *>         contains the first components of all right singular vectors
 126 *>         of the lower block. On exit, VF contains the first components
 127 *>         of all right singular vectors of the bidiagonal matrix.
 128 *> \endverbatim
 129 *>
 130 *> \param[out] VFW
 131 *> \verbatim
 132 *>          VFW is DOUBLE PRECISION array, dimension ( M )
 133 *>         Workspace for VF.
 134 *> \endverbatim
 135 *>
 136 *> \param[in,out] VL
 137 *> \verbatim
 138 *>          VL is DOUBLE PRECISION array, dimension ( M )
 139 *>         On entry, VL(1:NL+1) contains the  last components of all
 140 *>         right singular vectors of the upper block; and VL(NL+2:M)
 141 *>         contains the last components of all right singular vectors
 142 *>         of the lower block. On exit, VL contains the last components
 143 *>         of all right singular vectors of the bidiagonal matrix.
 144 *> \endverbatim
 145 *>
 146 *> \param[out] VLW
 147 *> \verbatim
 148 *>          VLW is DOUBLE PRECISION array, dimension ( M )
 149 *>         Workspace for VL.
 150 *> \endverbatim
 151 *>
 152 *> \param[in] ALPHA
 153 *> \verbatim
 154 *>          ALPHA is DOUBLE PRECISION
 155 *>         Contains the diagonal element associated with the added row.
 156 *> \endverbatim
 157 *>
 158 *> \param[in] BETA
 159 *> \verbatim
 160 *>          BETA is DOUBLE PRECISION
 161 *>         Contains the off-diagonal element associated with the added
 162 *>         row.
 163 *> \endverbatim
 164 *>
 165 *> \param[out] DSIGMA
 166 *> \verbatim
 167 *>          DSIGMA is DOUBLE PRECISION array, dimension ( N )
 168 *>         Contains a copy of the diagonal elements (K-1 singular values
 169 *>         and one zero) in the secular equation.
 170 *> \endverbatim
 171 *>
 172 *> \param[out] IDX
 173 *> \verbatim
 174 *>          IDX is INTEGER array, dimension ( N )
 175 *>         This will contain the permutation used to sort the contents of
 176 *>         D into ascending order.
 177 *> \endverbatim
 178 *>
 179 *> \param[out] IDXP
 180 *> \verbatim
 181 *>          IDXP is INTEGER array, dimension ( N )
 182 *>         This will contain the permutation used to place deflated
 183 *>         values of D at the end of the array. On output IDXP(2:K)
 184 *>         points to the nondeflated D-values and IDXP(K+1:N)
 185 *>         points to the deflated singular values.
 186 *> \endverbatim
 187 *>
 188 *> \param[in] IDXQ
 189 *> \verbatim
 190 *>          IDXQ is INTEGER array, dimension ( N )
 191 *>         This contains the permutation which separately sorts the two
 192 *>         sub-problems in D into ascending order.  Note that entries in
 193 *>         the first half of this permutation must first be moved one
 194 *>         position backward; and entries in the second half
 195 *>         must first have NL+1 added to their values.
 196 *> \endverbatim
 197 *>
 198 *> \param[out] PERM
 199 *> \verbatim
 200 *>          PERM is INTEGER array, dimension ( N )
 201 *>         The permutations (from deflation and sorting) to be applied
 202 *>         to each singular block. Not referenced if ICOMPQ = 0.
 203 *> \endverbatim
 204 *>
 205 *> \param[out] GIVPTR
 206 *> \verbatim
 207 *>          GIVPTR is INTEGER
 208 *>         The number of Givens rotations which took place in this
 209 *>         subproblem. Not referenced if ICOMPQ = 0.
 210 *> \endverbatim
 211 *>
 212 *> \param[out] GIVCOL
 213 *> \verbatim
 214 *>          GIVCOL is INTEGER array, dimension ( LDGCOL, 2 )
 215 *>         Each pair of numbers indicates a pair of columns to take place
 216 *>         in a Givens rotation. Not referenced if ICOMPQ = 0.
 217 *> \endverbatim
 218 *>
 219 *> \param[in] LDGCOL
 220 *> \verbatim
 221 *>          LDGCOL is INTEGER
 222 *>         The leading dimension of GIVCOL, must be at least N.
 223 *> \endverbatim
 224 *>
 225 *> \param[out] GIVNUM
 226 *> \verbatim
 227 *>          GIVNUM is DOUBLE PRECISION array, dimension ( LDGNUM, 2 )
 228 *>         Each number indicates the C or S value to be used in the
 229 *>         corresponding Givens rotation. Not referenced if ICOMPQ = 0.
 230 *> \endverbatim
 231 *>
 232 *> \param[in] LDGNUM
 233 *> \verbatim
 234 *>          LDGNUM is INTEGER
 235 *>         The leading dimension of GIVNUM, must be at least N.
 236 *> \endverbatim
 237 *>
 238 *> \param[out] C
 239 *> \verbatim
 240 *>          C is DOUBLE PRECISION
 241 *>         C contains garbage if SQRE =0 and the C-value of a Givens
 242 *>         rotation related to the right null space if SQRE = 1.
 243 *> \endverbatim
 244 *>
 245 *> \param[out] S
 246 *> \verbatim
 247 *>          S is DOUBLE PRECISION
 248 *>         S contains garbage if SQRE =0 and the S-value of a Givens
 249 *>         rotation related to the right null space if SQRE = 1.
 250 *> \endverbatim
 251 *>
 252 *> \param[out] INFO
 253 *> \verbatim
 254 *>          INFO is INTEGER
 255 *>         = 0:  successful exit.
 256 *>         < 0:  if INFO = -i, the i-th argument had an illegal value.
 257 *> \endverbatim
 258 *
 259 *  Authors:
 260 *  ========
 261 *
 262 *> \author Univ. of Tennessee
 263 *> \author Univ. of California Berkeley
 264 *> \author Univ. of Colorado Denver
 265 *> \author NAG Ltd.
 266 *
 267 *> \date September 2012
 268 *
 269 *> \ingroup auxOTHERauxiliary
 270 *
 271 *> \par Contributors:
 272 *  ==================
 273 *>
 274 *>     Ming Gu and Huan Ren, Computer Science Division, University of
 275 *>     California at Berkeley, USA
 276 *>
 277 *  =====================================================================
 278       SUBROUTINE DLASD7( ICOMPQ, NL, NR, SQRE, K, D, Z, ZW, VF, VFW, VL,
 279      $                   VLW, ALPHA, BETA, DSIGMA, IDX, IDXP, IDXQ,
 280      $                   PERM, GIVPTR, GIVCOL, LDGCOL, GIVNUM, LDGNUM,
 281      $                   C, S, INFO )
 282 *
 283 *  -- LAPACK auxiliary routine (version 3.4.2) --
 284 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
 285 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
 286 *     September 2012
 287 *
 288 *     .. Scalar Arguments ..
 289       INTEGER            GIVPTR, ICOMPQ, INFO, K, LDGCOL, LDGNUM, NL,
 290      $                   NR, SQRE
 291       DOUBLE PRECISION   ALPHA, BETA, C, S
 292 *     ..
 293 *     .. Array Arguments ..
 294       INTEGER            GIVCOL( LDGCOL, * ), IDX( * ), IDXP( * ),
 295      $                   IDXQ( * ), PERM( * )
 296       DOUBLE PRECISION   D( * ), DSIGMA( * ), GIVNUM( LDGNUM, * ),
 297      $                   VF( * ), VFW( * ), VL( * ), VLW( * ), Z( * ),
 298      $                   ZW( * )
 299 *     ..
 300 *
 301 *  =====================================================================
 302 *
 303 *     .. Parameters ..
 304       DOUBLE PRECISION   ZERO, ONE, TWO, EIGHT
 305       PARAMETER          ( ZERO = 0.0D+0, ONE = 1.0D+0, TWO = 2.0D+0,
 306      $                   EIGHT = 8.0D+0 )
 307 *     ..
 308 *     .. Local Scalars ..
 309 *
 310       INTEGER            I, IDXI, IDXJ, IDXJP, J, JP, JPREV, K2, M, N,
 311      $                   NLP1, NLP2
 312       DOUBLE PRECISION   EPS, HLFTOL, TAU, TOL, Z1
 313 *     ..
 314 *     .. External Subroutines ..
 315       EXTERNAL           DCOPY, DLAMRG, DROT, XERBLA
 316 *     ..
 317 *     .. External Functions ..
 318       DOUBLE PRECISION   DLAMCH, DLAPY2
 319       EXTERNAL           DLAMCH, DLAPY2
 320 *     ..
 321 *     .. Intrinsic Functions ..
 322       INTRINSIC          ABS, MAX
 323 *     ..
 324 *     .. Executable Statements ..
 325 *
 326 *     Test the input parameters.
 327 *
 328       INFO = 0
 329       N = NL + NR + 1
 330       M = N + SQRE
 331 *
 332       IF( ( ICOMPQ.LT.0 ) .OR. ( ICOMPQ.GT.1 ) ) THEN
 333          INFO = -1
 334       ELSE IF( NL.LT.1 ) THEN
 335          INFO = -2
 336       ELSE IF( NR.LT.1 ) THEN
 337          INFO = -3
 338       ELSE IF( ( SQRE.LT.0 ) .OR. ( SQRE.GT.1 ) ) THEN
 339          INFO = -4
 340       ELSE IF( LDGCOL.LT.N ) THEN
 341          INFO = -22
 342       ELSE IF( LDGNUM.LT.N ) THEN
 343          INFO = -24
 344       END IF
 345       IF( INFO.NE.0 ) THEN
 346          CALL XERBLA( 'DLASD7', -INFO )
 347          RETURN
 348       END IF
 349 *
 350       NLP1 = NL + 1
 351       NLP2 = NL + 2
 352       IF( ICOMPQ.EQ.1 ) THEN
 353          GIVPTR = 0
 354       END IF
 355 *
 356 *     Generate the first part of the vector Z and move the singular
 357 *     values in the first part of D one position backward.
 358 *
 359       Z1 = ALPHA*VL( NLP1 )
 360       VL( NLP1 ) = ZERO
 361       TAU = VF( NLP1 )
 362       DO 10 I = NL, 1, -1
 363          Z( I+1 ) = ALPHA*VL( I )
 364          VL( I ) = ZERO
 365          VF( I+1 ) = VF( I )
 366          D( I+1 ) = D( I )
 367          IDXQ( I+1 ) = IDXQ( I ) + 1
 368    10 CONTINUE
 369       VF( 1 ) = TAU
 370 *
 371 *     Generate the second part of the vector Z.
 372 *
 373       DO 20 I = NLP2, M
 374          Z( I ) = BETA*VF( I )
 375          VF( I ) = ZERO
 376    20 CONTINUE
 377 *
 378 *     Sort the singular values into increasing order
 379 *
 380       DO 30 I = NLP2, N
 381          IDXQ( I ) = IDXQ( I ) + NLP1
 382    30 CONTINUE
 383 *
 384 *     DSIGMA, IDXC, IDXC, and ZW are used as storage space.
 385 *
 386       DO 40 I = 2, N
 387          DSIGMA( I ) = D( IDXQ( I ) )
 388          ZW( I ) = Z( IDXQ( I ) )
 389          VFW( I ) = VF( IDXQ( I ) )
 390          VLW( I ) = VL( IDXQ( I ) )
 391    40 CONTINUE
 392 *
 393       CALL DLAMRG( NL, NR, DSIGMA( 2 ), 1, 1, IDX( 2 ) )
 394 *
 395       DO 50 I = 2, N
 396          IDXI = 1 + IDX( I )
 397          D( I ) = DSIGMA( IDXI )
 398          Z( I ) = ZW( IDXI )
 399          VF( I ) = VFW( IDXI )
 400          VL( I ) = VLW( IDXI )
 401    50 CONTINUE
 402 *
 403 *     Calculate the allowable deflation tolerence
 404 *
 405       EPS = DLAMCH( 'Epsilon' )
 406       TOL = MAX( ABS( ALPHA ), ABS( BETA ) )
 407       TOL = EIGHT*EIGHT*EPS*MAX( ABS( D( N ) ), TOL )
 408 *
 409 *     There are 2 kinds of deflation -- first a value in the z-vector
 410 *     is small, second two (or more) singular values are very close
 411 *     together (their difference is small).
 412 *
 413 *     If the value in the z-vector is small, we simply permute the
 414 *     array so that the corresponding singular value is moved to the
 415 *     end.
 416 *
 417 *     If two values in the D-vector are close, we perform a two-sided
 418 *     rotation designed to make one of the corresponding z-vector
 419 *     entries zero, and then permute the array so that the deflated
 420 *     singular value is moved to the end.
 421 *
 422 *     If there are multiple singular values then the problem deflates.
 423 *     Here the number of equal singular values are found.  As each equal
 424 *     singular value is found, an elementary reflector is computed to
 425 *     rotate the corresponding singular subspace so that the
 426 *     corresponding components of Z are zero in this new basis.
 427 *
 428       K = 1
 429       K2 = N + 1
 430       DO 60 J = 2, N
 431          IF( ABS( Z( J ) ).LE.TOL ) THEN
 432 *
 433 *           Deflate due to small z component.
 434 *
 435             K2 = K2 - 1
 436             IDXP( K2 ) = J
 437             IF( J.EQ.N )
 438      $         GO TO 100
 439          ELSE
 440             JPREV = J
 441             GO TO 70
 442          END IF
 443    60 CONTINUE
 444    70 CONTINUE
 445       J = JPREV
 446    80 CONTINUE
 447       J = J + 1
 448       IF( J.GT.N )
 449      $   GO TO 90
 450       IF( ABS( Z( J ) ).LE.TOL ) THEN
 451 *
 452 *        Deflate due to small z component.
 453 *
 454          K2 = K2 - 1
 455          IDXP( K2 ) = J
 456       ELSE
 457 *
 458 *        Check if singular values are close enough to allow deflation.
 459 *
 460          IF( ABS( D( J )-D( JPREV ) ).LE.TOL ) THEN
 461 *
 462 *           Deflation is possible.
 463 *
 464             S = Z( JPREV )
 465             C = Z( J )
 466 *
 467 *           Find sqrt(a**2+b**2) without overflow or
 468 *           destructive underflow.
 469 *
 470             TAU = DLAPY2( C, S )
 471             Z( J ) = TAU
 472             Z( JPREV ) = ZERO
 473             C = C / TAU
 474             S = -S / TAU
 475 *
 476 *           Record the appropriate Givens rotation
 477 *
 478             IF( ICOMPQ.EQ.1 ) THEN
 479                GIVPTR = GIVPTR + 1
 480                IDXJP = IDXQ( IDX( JPREV )+1 )
 481                IDXJ = IDXQ( IDX( J )+1 )
 482                IF( IDXJP.LE.NLP1 ) THEN
 483                   IDXJP = IDXJP - 1
 484                END IF
 485                IF( IDXJ.LE.NLP1 ) THEN
 486                   IDXJ = IDXJ - 1
 487                END IF
 488                GIVCOL( GIVPTR, 2 ) = IDXJP
 489                GIVCOL( GIVPTR, 1 ) = IDXJ
 490                GIVNUM( GIVPTR, 2 ) = C
 491                GIVNUM( GIVPTR, 1 ) = S
 492             END IF
 493             CALL DROT( 1, VF( JPREV ), 1, VF( J ), 1, C, S )
 494             CALL DROT( 1, VL( JPREV ), 1, VL( J ), 1, C, S )
 495             K2 = K2 - 1
 496             IDXP( K2 ) = JPREV
 497             JPREV = J
 498          ELSE
 499             K = K + 1
 500             ZW( K ) = Z( JPREV )
 501             DSIGMA( K ) = D( JPREV )
 502             IDXP( K ) = JPREV
 503             JPREV = J
 504          END IF
 505       END IF
 506       GO TO 80
 507    90 CONTINUE
 508 *
 509 *     Record the last singular value.
 510 *
 511       K = K + 1
 512       ZW( K ) = Z( JPREV )
 513       DSIGMA( K ) = D( JPREV )
 514       IDXP( K ) = JPREV
 515 *
 516   100 CONTINUE
 517 *
 518 *     Sort the singular values into DSIGMA. The singular values which
 519 *     were not deflated go into the first K slots of DSIGMA, except
 520 *     that DSIGMA(1) is treated separately.
 521 *
 522       DO 110 J = 2, N
 523          JP = IDXP( J )
 524          DSIGMA( J ) = D( JP )
 525          VFW( J ) = VF( JP )
 526          VLW( J ) = VL( JP )
 527   110 CONTINUE
 528       IF( ICOMPQ.EQ.1 ) THEN
 529          DO 120 J = 2, N
 530             JP = IDXP( J )
 531             PERM( J ) = IDXQ( IDX( JP )+1 )
 532             IF( PERM( J ).LE.NLP1 ) THEN
 533                PERM( J ) = PERM( J ) - 1
 534             END IF
 535   120    CONTINUE
 536       END IF
 537 *
 538 *     The deflated singular values go back into the last N - K slots of
 539 *     D.
 540 *
 541       CALL DCOPY( N-K, DSIGMA( K+1 ), 1, D( K+1 ), 1 )
 542 *
 543 *     Determine DSIGMA(1), DSIGMA(2), Z(1), VF(1), VL(1), VF(M), and
 544 *     VL(M).
 545 *
 546       DSIGMA( 1 ) = ZERO
 547       HLFTOL = TOL / TWO
 548       IF( ABS( DSIGMA( 2 ) ).LE.HLFTOL )
 549      $   DSIGMA( 2 ) = HLFTOL
 550       IF( M.GT.N ) THEN
 551          Z( 1 ) = DLAPY2( Z1, Z( M ) )
 552          IF( Z( 1 ).LE.TOL ) THEN
 553             C = ONE
 554             S = ZERO
 555             Z( 1 ) = TOL
 556          ELSE
 557             C = Z1 / Z( 1 )
 558             S = -Z( M ) / Z( 1 )
 559          END IF
 560          CALL DROT( 1, VF( M ), 1, VF( 1 ), 1, C, S )
 561          CALL DROT( 1, VL( M ), 1, VL( 1 ), 1, C, S )
 562       ELSE
 563          IF( ABS( Z1 ).LE.TOL ) THEN
 564             Z( 1 ) = TOL
 565          ELSE
 566             Z( 1 ) = Z1
 567          END IF
 568       END IF
 569 *
 570 *     Restore Z, VF, and VL.
 571 *
 572       CALL DCOPY( K-1, ZW( 2 ), 1, Z( 2 ), 1 )
 573       CALL DCOPY( N-1, VFW( 2 ), 1, VF( 2 ), 1 )
 574       CALL DCOPY( N-1, VLW( 2 ), 1, VL( 2 ), 1 )
 575 *
 576       RETURN
 577 *
 578 *     End of DLASD7
 579 *
 580       END