SRC/slasd6.f

   1 *> \brief \b SLASD6 computes the SVD of an updated upper bidiagonal matrix obtained by merging two smaller ones by appending a row. 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 SLASD6 + dependencies
  10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/slasd6.f">
  11 *> [TGZ]</a>
  12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/slasd6.f">
  13 *> [ZIP]</a>
  14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/slasd6.f">
  15 *> [TXT]</a>
  16 *> \endhtmlonly
  17 *
  18 *  Definition:
  19 *  ===========
  20 *
  21 *       SUBROUTINE SLASD6( ICOMPQ, NL, NR, SQRE, D, VF, VL, ALPHA, BETA,
  22 *                          IDXQ, PERM, GIVPTR, GIVCOL, LDGCOL, GIVNUM,
  23 *                          LDGNUM, POLES, DIFL, DIFR, Z, K, C, S, WORK,
  24 *                          IWORK, INFO )
  25 *
  26 *       .. Scalar Arguments ..
  27 *       INTEGER            GIVPTR, ICOMPQ, INFO, K, LDGCOL, LDGNUM, NL,
  28 *      $                   NR, SQRE
  29 *       REAL               ALPHA, BETA, C, S
  30 *       ..
  31 *       .. Array Arguments ..
  32 *       INTEGER            GIVCOL( LDGCOL, * ), IDXQ( * ), IWORK( * ),
  33 *      $                   PERM( * )
  34 *       REAL               D( * ), DIFL( * ), DIFR( * ),
  35 *      $                   GIVNUM( LDGNUM, * ), POLES( LDGNUM, * ),
  36 *      $                   VF( * ), VL( * ), WORK( * ), Z( * )
  37 *       ..
  38 *
  39 *
  40 *> \par Purpose:
  41 *  =============
  42 *>
  43 *> \verbatim
  44 *>
  45 *> SLASD6 computes the SVD of an updated upper bidiagonal matrix B
  46 *> obtained by merging two smaller ones by appending a row. This
  47 *> routine is used only for the problem which requires all singular
  48 *> values and optionally singular vector matrices in factored form.
  49 *> B is an N-by-M matrix with N = NL + NR + 1 and M = N + SQRE.
  50 *> A related subroutine, SLASD1, handles the case in which all singular
  51 *> values and singular vectors of the bidiagonal matrix are desired.
  52 *>
  53 *> SLASD6 computes the SVD as follows:
  54 *>
  55 *>               ( D1(in)    0    0       0 )
  56 *>   B = U(in) * (   Z1**T   a   Z2**T    b ) * VT(in)
  57 *>               (   0       0   D2(in)   0 )
  58 *>
  59 *>     = U(out) * ( D(out) 0) * VT(out)
  60 *>
  61 *> where Z**T = (Z1**T a Z2**T b) = u**T VT**T, and u is a vector of dimension M
  62 *> with ALPHA and BETA in the NL+1 and NL+2 th entries and zeros
  63 *> elsewhere; and the entry b is empty if SQRE = 0.
  64 *>
  65 *> The singular values of B can be computed using D1, D2, the first
  66 *> components of all the right singular vectors of the lower block, and
  67 *> the last components of all the right singular vectors of the upper
  68 *> block. These components are stored and updated in VF and VL,
  69 *> respectively, in SLASD6. Hence U and VT are not explicitly
  70 *> referenced.
  71 *>
  72 *> The singular values are stored in D. The algorithm consists of two
  73 *> stages:
  74 *>
  75 *>       The first stage consists of deflating the size of the problem
  76 *>       when there are multiple singular values or if there is a zero
  77 *>       in the Z vector. For each such occurrence the dimension of the
  78 *>       secular equation problem is reduced by one. This stage is
  79 *>       performed by the routine SLASD7.
  80 *>
  81 *>       The second stage consists of calculating the updated
  82 *>       singular values. This is done by finding the roots of the
  83 *>       secular equation via the routine SLASD4 (as called by SLASD8).
  84 *>       This routine also updates VF and VL and computes the distances
  85 *>       between the updated singular values and the old singular
  86 *>       values.
  87 *>
  88 *> SLASD6 is called from SLASDA.
  89 *> \endverbatim
  90 *
  91 *  Arguments:
  92 *  ==========
  93 *
  94 *> \param[in] ICOMPQ
  95 *> \verbatim
  96 *>          ICOMPQ is INTEGER
  97 *>         Specifies whether singular vectors are to be computed in
  98 *>         factored form:
  99 *>         = 0: Compute singular values only.
 100 *>         = 1: Compute singular vectors in factored form as well.
 101 *> \endverbatim
 102 *>
 103 *> \param[in] NL
 104 *> \verbatim
 105 *>          NL is INTEGER
 106 *>         The row dimension of the upper block.  NL >= 1.
 107 *> \endverbatim
 108 *>
 109 *> \param[in] NR
 110 *> \verbatim
 111 *>          NR is INTEGER
 112 *>         The row dimension of the lower block.  NR >= 1.
 113 *> \endverbatim
 114 *>
 115 *> \param[in] SQRE
 116 *> \verbatim
 117 *>          SQRE is INTEGER
 118 *>         = 0: the lower block is an NR-by-NR square matrix.
 119 *>         = 1: the lower block is an NR-by-(NR+1) rectangular matrix.
 120 *>
 121 *>         The bidiagonal matrix has row dimension N = NL + NR + 1,
 122 *>         and column dimension M = N + SQRE.
 123 *> \endverbatim
 124 *>
 125 *> \param[in,out] D
 126 *> \verbatim
 127 *>          D is REAL array, dimension (NL+NR+1).
 128 *>         On entry D(1:NL,1:NL) contains the singular values of the
 129 *>         upper block, and D(NL+2:N) contains the singular values
 130 *>         of the lower block. On exit D(1:N) contains the singular
 131 *>         values of the modified matrix.
 132 *> \endverbatim
 133 *>
 134 *> \param[in,out] VF
 135 *> \verbatim
 136 *>          VF is REAL array, dimension (M)
 137 *>         On entry, VF(1:NL+1) contains the first components of all
 138 *>         right singular vectors of the upper block; and VF(NL+2:M)
 139 *>         contains the first components of all right singular vectors
 140 *>         of the lower block. On exit, VF contains the first components
 141 *>         of all right singular vectors of the bidiagonal matrix.
 142 *> \endverbatim
 143 *>
 144 *> \param[in,out] VL
 145 *> \verbatim
 146 *>          VL is REAL array, dimension (M)
 147 *>         On entry, VL(1:NL+1) contains the  last components of all
 148 *>         right singular vectors of the upper block; and VL(NL+2:M)
 149 *>         contains the last components of all right singular vectors of
 150 *>         the lower block. On exit, VL contains the last components of
 151 *>         all right singular vectors of the bidiagonal matrix.
 152 *> \endverbatim
 153 *>
 154 *> \param[in,out] ALPHA
 155 *> \verbatim
 156 *>          ALPHA is REAL
 157 *>         Contains the diagonal element associated with the added row.
 158 *> \endverbatim
 159 *>
 160 *> \param[in,out] BETA
 161 *> \verbatim
 162 *>          BETA is REAL
 163 *>         Contains the off-diagonal element associated with the added
 164 *>         row.
 165 *> \endverbatim
 166 *>
 167 *> \param[in,out] IDXQ
 168 *> \verbatim
 169 *>          IDXQ is INTEGER array, dimension (N)
 170 *>         This contains the permutation which will reintegrate the
 171 *>         subproblem just solved back into sorted order, i.e.
 172 *>         D( IDXQ( I = 1, N ) ) will be in ascending order.
 173 *> \endverbatim
 174 *>
 175 *> \param[out] PERM
 176 *> \verbatim
 177 *>          PERM is INTEGER array, dimension ( N )
 178 *>         The permutations (from deflation and sorting) to be applied
 179 *>         to each block. Not referenced if ICOMPQ = 0.
 180 *> \endverbatim
 181 *>
 182 *> \param[out] GIVPTR
 183 *> \verbatim
 184 *>          GIVPTR is INTEGER
 185 *>         The number of Givens rotations which took place in this
 186 *>         subproblem. Not referenced if ICOMPQ = 0.
 187 *> \endverbatim
 188 *>
 189 *> \param[out] GIVCOL
 190 *> \verbatim
 191 *>          GIVCOL is INTEGER array, dimension ( LDGCOL, 2 )
 192 *>         Each pair of numbers indicates a pair of columns to take place
 193 *>         in a Givens rotation. Not referenced if ICOMPQ = 0.
 194 *> \endverbatim
 195 *>
 196 *> \param[in] LDGCOL
 197 *> \verbatim
 198 *>          LDGCOL is INTEGER
 199 *>         leading dimension of GIVCOL, must be at least N.
 200 *> \endverbatim
 201 *>
 202 *> \param[out] GIVNUM
 203 *> \verbatim
 204 *>          GIVNUM is REAL array, dimension ( LDGNUM, 2 )
 205 *>         Each number indicates the C or S value to be used in the
 206 *>         corresponding Givens rotation. Not referenced if ICOMPQ = 0.
 207 *> \endverbatim
 208 *>
 209 *> \param[in] LDGNUM
 210 *> \verbatim
 211 *>          LDGNUM is INTEGER
 212 *>         The leading dimension of GIVNUM and POLES, must be at least N.
 213 *> \endverbatim
 214 *>
 215 *> \param[out] POLES
 216 *> \verbatim
 217 *>          POLES is REAL array, dimension ( LDGNUM, 2 )
 218 *>         On exit, POLES(1,*) is an array containing the new singular
 219 *>         values obtained from solving the secular equation, and
 220 *>         POLES(2,*) is an array containing the poles in the secular
 221 *>         equation. Not referenced if ICOMPQ = 0.
 222 *> \endverbatim
 223 *>
 224 *> \param[out] DIFL
 225 *> \verbatim
 226 *>          DIFL is REAL array, dimension ( N )
 227 *>         On exit, DIFL(I) is the distance between I-th updated
 228 *>         (undeflated) singular value and the I-th (undeflated) old
 229 *>         singular value.
 230 *> \endverbatim
 231 *>
 232 *> \param[out] DIFR
 233 *> \verbatim
 234 *>          DIFR is REAL array,
 235 *>                   dimension ( LDDIFR, 2 ) if ICOMPQ = 1 and
 236 *>                   dimension ( K ) if ICOMPQ = 0.
 237 *>          On exit, DIFR(I,1) = D(I) - DSIGMA(I+1), DIFR(K,1) is not
 238 *>          defined and will not be referenced.
 239 *>
 240 *>          If ICOMPQ = 1, DIFR(1:K,2) is an array containing the
 241 *>          normalizing factors for the right singular vector matrix.
 242 *>
 243 *>         See SLASD8 for details on DIFL and DIFR.
 244 *> \endverbatim
 245 *>
 246 *> \param[out] Z
 247 *> \verbatim
 248 *>          Z is REAL array, dimension ( M )
 249 *>         The first elements of this array contain the components
 250 *>         of the deflation-adjusted updating row vector.
 251 *> \endverbatim
 252 *>
 253 *> \param[out] K
 254 *> \verbatim
 255 *>          K is INTEGER
 256 *>         Contains the dimension of the non-deflated matrix,
 257 *>         This is the order of the related secular equation. 1 <= K <=N.
 258 *> \endverbatim
 259 *>
 260 *> \param[out] C
 261 *> \verbatim
 262 *>          C is REAL
 263 *>         C contains garbage if SQRE =0 and the C-value of a Givens
 264 *>         rotation related to the right null space if SQRE = 1.
 265 *> \endverbatim
 266 *>
 267 *> \param[out] S
 268 *> \verbatim
 269 *>          S is REAL
 270 *>         S contains garbage if SQRE =0 and the S-value of a Givens
 271 *>         rotation related to the right null space if SQRE = 1.
 272 *> \endverbatim
 273 *>
 274 *> \param[out] WORK
 275 *> \verbatim
 276 *>          WORK is REAL array, dimension ( 4 * M )
 277 *> \endverbatim
 278 *>
 279 *> \param[out] IWORK
 280 *> \verbatim
 281 *>          IWORK is INTEGER array, dimension ( 3 * N )
 282 *> \endverbatim
 283 *>
 284 *> \param[out] INFO
 285 *> \verbatim
 286 *>          INFO is INTEGER
 287 *>          = 0:  successful exit.
 288 *>          < 0:  if INFO = -i, the i-th argument had an illegal value.
 289 *>          > 0:  if INFO = 1, a singular value did not converge
 290 *> \endverbatim
 291 *
 292 *  Authors:
 293 *  ========
 294 *
 295 *> \author Univ. of Tennessee
 296 *> \author Univ. of California Berkeley
 297 *> \author Univ. of Colorado Denver
 298 *> \author NAG Ltd.
 299 *
 300 *> \date June 2016
 301 *
 302 *> \ingroup OTHERauxiliary
 303 *
 304 *> \par Contributors:
 305 *  ==================
 306 *>
 307 *>     Ming Gu and Huan Ren, Computer Science Division, University of
 308 *>     California at Berkeley, USA
 309 *>
 310 *  =====================================================================
 311       SUBROUTINE SLASD6( ICOMPQ, NL, NR, SQRE, D, VF, VL, ALPHA, BETA,
 312      $                   IDXQ, PERM, GIVPTR, GIVCOL, LDGCOL, GIVNUM,
 313      $                   LDGNUM, POLES, DIFL, DIFR, Z, K, C, S, WORK,
 314      $                   IWORK, INFO )
 315 *
 316 *  -- LAPACK auxiliary routine (version 3.6.1) --
 317 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
 318 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
 319 *     June 2016
 320 *
 321 *     .. Scalar Arguments ..
 322       INTEGER            GIVPTR, ICOMPQ, INFO, K, LDGCOL, LDGNUM, NL,
 323      $                   NR, SQRE
 324       REAL               ALPHA, BETA, C, S
 325 *     ..
 326 *     .. Array Arguments ..
 327       INTEGER            GIVCOL( LDGCOL, * ), IDXQ( * ), IWORK( * ),
 328      $                   PERM( * )
 329       REAL               D( * ), DIFL( * ), DIFR( * ),
 330      $                   GIVNUM( LDGNUM, * ), POLES( LDGNUM, * ),
 331      $                   VF( * ), VL( * ), WORK( * ), Z( * )
 332 *     ..
 333 *
 334 *  =====================================================================
 335 *
 336 *     .. Parameters ..
 337       REAL               ONE, ZERO
 338       PARAMETER          ( ONE = 1.0E+0, ZERO = 0.0E+0 )
 339 *     ..
 340 *     .. Local Scalars ..
 341       INTEGER            I, IDX, IDXC, IDXP, ISIGMA, IVFW, IVLW, IW, M,
 342      $                   N, N1, N2
 343       REAL               ORGNRM
 344 *     ..
 345 *     .. External Subroutines ..
 346       EXTERNAL           SCOPY, SLAMRG, SLASCL, SLASD7, SLASD8, XERBLA
 347 *     ..
 348 *     .. Intrinsic Functions ..
 349       INTRINSIC          ABS, MAX
 350 *     ..
 351 *     .. Executable Statements ..
 352 *
 353 *     Test the input parameters.
 354 *
 355       INFO = 0
 356       N = NL + NR + 1
 357       M = N + SQRE
 358 *
 359       IF( ( ICOMPQ.LT.0 ) .OR. ( ICOMPQ.GT.1 ) ) THEN
 360          INFO = -1
 361       ELSE IF( NL.LT.1 ) THEN
 362          INFO = -2
 363       ELSE IF( NR.LT.1 ) THEN
 364          INFO = -3
 365       ELSE IF( ( SQRE.LT.0 ) .OR. ( SQRE.GT.1 ) ) THEN
 366          INFO = -4
 367       ELSE IF( LDGCOL.LT.N ) THEN
 368          INFO = -14
 369       ELSE IF( LDGNUM.LT.N ) THEN
 370          INFO = -16
 371       END IF
 372       IF( INFO.NE.0 ) THEN
 373          CALL XERBLA( 'SLASD6', -INFO )
 374          RETURN
 375       END IF
 376 *
 377 *     The following values are for bookkeeping purposes only.  They are
 378 *     integer pointers which indicate the portion of the workspace
 379 *     used by a particular array in SLASD7 and SLASD8.
 380 *
 381       ISIGMA = 1
 382       IW = ISIGMA + N
 383       IVFW = IW + M
 384       IVLW = IVFW + M
 385 *
 386       IDX = 1
 387       IDXC = IDX + N
 388       IDXP = IDXC + N
 389 *
 390 *     Scale.
 391 *
 392       ORGNRM = MAX( ABS( ALPHA ), ABS( BETA ) )
 393       D( NL+1 ) = ZERO
 394       DO 10 I = 1, N
 395          IF( ABS( D( I ) ).GT.ORGNRM ) THEN
 396             ORGNRM = ABS( D( I ) )
 397          END IF
 398    10 CONTINUE
 399       CALL SLASCL( 'G', 0, 0, ORGNRM, ONE, N, 1, D, N, INFO )
 400       ALPHA = ALPHA / ORGNRM
 401       BETA = BETA / ORGNRM
 402 *
 403 *     Sort and Deflate singular values.
 404 *
 405       CALL SLASD7( ICOMPQ, NL, NR, SQRE, K, D, Z, WORK( IW ), VF,
 406      $             WORK( IVFW ), VL, WORK( IVLW ), ALPHA, BETA,
 407      $             WORK( ISIGMA ), IWORK( IDX ), IWORK( IDXP ), IDXQ,
 408      $             PERM, GIVPTR, GIVCOL, LDGCOL, GIVNUM, LDGNUM, C, S,
 409      $             INFO )
 410 *
 411 *     Solve Secular Equation, compute DIFL, DIFR, and update VF, VL.
 412 *
 413       CALL SLASD8( ICOMPQ, K, D, Z, VF, VL, DIFL, DIFR, LDGNUM,
 414      $             WORK( ISIGMA ), WORK( IW ), INFO )
 415 *
 416 *     Report the possible convergence failure.
 417 *
 418       IF( INFO.NE.0 ) THEN
 419          RETURN
 420       END IF
 421 *
 422 *     Save the poles if ICOMPQ = 1.
 423 *
 424       IF( ICOMPQ.EQ.1 ) THEN
 425          CALL SCOPY( K, D, 1, POLES( 1, 1 ), 1 )
 426          CALL SCOPY( K, WORK( ISIGMA ), 1, POLES( 1, 2 ), 1 )
 427       END IF
 428 *
 429 *     Unscale.
 430 *
 431       CALL SLASCL( 'G', 0, 0, ONE, ORGNRM, N, 1, D, N, INFO )
 432 *
 433 *     Prepare the IDXQ sorting permutation.
 434 *
 435       N1 = K
 436       N2 = N - K
 437       CALL SLAMRG( N1, N2, D, 1, -1, IDXQ )
 438 *
 439       RETURN
 440 *
 441 *     End of SLASD6
 442 *
 443       END