SRC/slasd8.f

   1 *> \brief \b SLASD8 finds the square roots of the roots of the secular equation, and stores, for each element in D, the distance to its two nearest poles. 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 SLASD8 + dependencies
  10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/slasd8.f">
  11 *> [TGZ]</a>
  12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/slasd8.f">
  13 *> [ZIP]</a>
  14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/slasd8.f">
  15 *> [TXT]</a>
  16 *> \endhtmlonly
  17 *
  18 *  Definition:
  19 *  ===========
  20 *
  21 *       SUBROUTINE SLASD8( ICOMPQ, K, D, Z, VF, VL, DIFL, DIFR, LDDIFR,
  22 *                          DSIGMA, WORK, INFO )
  23 *
  24 *       .. Scalar Arguments ..
  25 *       INTEGER            ICOMPQ, INFO, K, LDDIFR
  26 *       ..
  27 *       .. Array Arguments ..
  28 *       REAL               D( * ), DIFL( * ), DIFR( LDDIFR, * ),
  29 *      $                   DSIGMA( * ), VF( * ), VL( * ), WORK( * ),
  30 *      $                   Z( * )
  31 *       ..
  32 *
  33 *
  34 *> \par Purpose:
  35 *  =============
  36 *>
  37 *> \verbatim
  38 *>
  39 *> SLASD8 finds the square roots of the roots of the secular equation,
  40 *> as defined by the values in DSIGMA and Z. It makes the appropriate
  41 *> calls to SLASD4, and stores, for each  element in D, the distance
  42 *> to its two nearest poles (elements in DSIGMA). It also updates
  43 *> the arrays VF and VL, the first and last components of all the
  44 *> right singular vectors of the original bidiagonal matrix.
  45 *>
  46 *> SLASD8 is called from SLASD6.
  47 *> \endverbatim
  48 *
  49 *  Arguments:
  50 *  ==========
  51 *
  52 *> \param[in] ICOMPQ
  53 *> \verbatim
  54 *>          ICOMPQ is INTEGER
  55 *>          Specifies whether singular vectors are to be computed in
  56 *>          factored form in the calling routine:
  57 *>          = 0: Compute singular values only.
  58 *>          = 1: Compute singular vectors in factored form as well.
  59 *> \endverbatim
  60 *>
  61 *> \param[in] K
  62 *> \verbatim
  63 *>          K is INTEGER
  64 *>          The number of terms in the rational function to be solved
  65 *>          by SLASD4.  K >= 1.
  66 *> \endverbatim
  67 *>
  68 *> \param[out] D
  69 *> \verbatim
  70 *>          D is REAL array, dimension ( K )
  71 *>          On output, D contains the updated singular values.
  72 *> \endverbatim
  73 *>
  74 *> \param[in,out] Z
  75 *> \verbatim
  76 *>          Z is REAL array, dimension ( K )
  77 *>          On entry, the first K elements of this array contain the
  78 *>          components of the deflation-adjusted updating row vector.
  79 *>          On exit, Z is updated.
  80 *> \endverbatim
  81 *>
  82 *> \param[in,out] VF
  83 *> \verbatim
  84 *>          VF is REAL array, dimension ( K )
  85 *>          On entry, VF contains  information passed through DBEDE8.
  86 *>          On exit, VF contains the first K components of the first
  87 *>          components of all right singular vectors of the bidiagonal
  88 *>          matrix.
  89 *> \endverbatim
  90 *>
  91 *> \param[in,out] VL
  92 *> \verbatim
  93 *>          VL is REAL array, dimension ( K )
  94 *>          On entry, VL contains  information passed through DBEDE8.
  95 *>          On exit, VL contains the first K components of the last
  96 *>          components of all right singular vectors of the bidiagonal
  97 *>          matrix.
  98 *> \endverbatim
  99 *>
 100 *> \param[out] DIFL
 101 *> \verbatim
 102 *>          DIFL is REAL array, dimension ( K )
 103 *>          On exit, DIFL(I) = D(I) - DSIGMA(I).
 104 *> \endverbatim
 105 *>
 106 *> \param[out] DIFR
 107 *> \verbatim
 108 *>          DIFR is REAL array,
 109 *>                   dimension ( LDDIFR, 2 ) if ICOMPQ = 1 and
 110 *>                   dimension ( K ) if ICOMPQ = 0.
 111 *>          On exit, DIFR(I,1) = D(I) - DSIGMA(I+1), DIFR(K,1) is not
 112 *>          defined and will not be referenced.
 113 *>
 114 *>          If ICOMPQ = 1, DIFR(1:K,2) is an array containing the
 115 *>          normalizing factors for the right singular vector matrix.
 116 *> \endverbatim
 117 *>
 118 *> \param[in] LDDIFR
 119 *> \verbatim
 120 *>          LDDIFR is INTEGER
 121 *>          The leading dimension of DIFR, must be at least K.
 122 *> \endverbatim
 123 *>
 124 *> \param[in,out] DSIGMA
 125 *> \verbatim
 126 *>          DSIGMA is REAL array, dimension ( K )
 127 *>          On entry, the first K elements of this array contain the old
 128 *>          roots of the deflated updating problem.  These are the poles
 129 *>          of the secular equation.
 130 *>          On exit, the elements of DSIGMA may be very slightly altered
 131 *>          in value.
 132 *> \endverbatim
 133 *>
 134 *> \param[out] WORK
 135 *> \verbatim
 136 *>          WORK is REAL array, dimension at least 3 * K
 137 *> \endverbatim
 138 *>
 139 *> \param[out] INFO
 140 *> \verbatim
 141 *>          INFO is INTEGER
 142 *>          = 0:  successful exit.
 143 *>          < 0:  if INFO = -i, the i-th argument had an illegal value.
 144 *>          > 0:  if INFO = 1, a singular value did not converge
 145 *> \endverbatim
 146 *
 147 *  Authors:
 148 *  ========
 149 *
 150 *> \author Univ. of Tennessee
 151 *> \author Univ. of California Berkeley
 152 *> \author Univ. of Colorado Denver
 153 *> \author NAG Ltd.
 154 *
 155 *> \date November 2015
 156 *
 157 *> \ingroup auxOTHERauxiliary
 158 *
 159 *> \par Contributors:
 160 *  ==================
 161 *>
 162 *>     Ming Gu and Huan Ren, Computer Science Division, University of
 163 *>     California at Berkeley, USA
 164 *>
 165 *  =====================================================================
 166       SUBROUTINE SLASD8( ICOMPQ, K, D, Z, VF, VL, DIFL, DIFR, LDDIFR,
 167      $                   DSIGMA, WORK, INFO )
 168 *
 169 *  -- LAPACK auxiliary routine (version 3.6.0) --
 170 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
 171 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
 172 *     November 2015
 173 *
 174 *     .. Scalar Arguments ..
 175       INTEGER            ICOMPQ, INFO, K, LDDIFR
 176 *     ..
 177 *     .. Array Arguments ..
 178       REAL               D( * ), DIFL( * ), DIFR( LDDIFR, * ),
 179      $                   DSIGMA( * ), VF( * ), VL( * ), WORK( * ),
 180      $                   Z( * )
 181 *     ..
 182 *
 183 *  =====================================================================
 184 *
 185 *     .. Parameters ..
 186       REAL               ONE
 187       PARAMETER          ( ONE = 1.0E+0 )
 188 *     ..
 189 *     .. Local Scalars ..
 190       INTEGER            I, IWK1, IWK2, IWK2I, IWK3, IWK3I, J
 191       REAL               DIFLJ, DIFRJ, DJ, DSIGJ, DSIGJP, RHO, TEMP
 192 *     ..
 193 *     .. External Subroutines ..
 194       EXTERNAL           SCOPY, SLASCL, SLASD4, SLASET, XERBLA
 195 *     ..
 196 *     .. External Functions ..
 197       REAL               SDOT, SLAMC3, SNRM2
 198       EXTERNAL           SDOT, SLAMC3, SNRM2
 199 *     ..
 200 *     .. Intrinsic Functions ..
 201       INTRINSIC          ABS, SIGN, SQRT
 202 *     ..
 203 *     .. Executable Statements ..
 204 *
 205 *     Test the input parameters.
 206 *
 207       INFO = 0
 208 *
 209       IF( ( ICOMPQ.LT.0 ) .OR. ( ICOMPQ.GT.1 ) ) THEN
 210          INFO = -1
 211       ELSE IF( K.LT.1 ) THEN
 212          INFO = -2
 213       ELSE IF( LDDIFR.LT.K ) THEN
 214          INFO = -9
 215       END IF
 216       IF( INFO.NE.0 ) THEN
 217          CALL XERBLA( 'SLASD8', -INFO )
 218          RETURN
 219       END IF
 220 *
 221 *     Quick return if possible
 222 *
 223       IF( K.EQ.1 ) THEN
 224          D( 1 ) = ABS( Z( 1 ) )
 225          DIFL( 1 ) = D( 1 )
 226          IF( ICOMPQ.EQ.1 ) THEN
 227             DIFL( 2 ) = ONE
 228             DIFR( 1, 2 ) = ONE
 229          END IF
 230          RETURN
 231       END IF
 232 *
 233 *     Modify values DSIGMA(i) to make sure all DSIGMA(i)-DSIGMA(j) can
 234 *     be computed with high relative accuracy (barring over/underflow).
 235 *     This is a problem on machines without a guard digit in
 236 *     add/subtract (Cray XMP, Cray YMP, Cray C 90 and Cray 2).
 237 *     The following code replaces DSIGMA(I) by 2*DSIGMA(I)-DSIGMA(I),
 238 *     which on any of these machines zeros out the bottommost
 239 *     bit of DSIGMA(I) if it is 1; this makes the subsequent
 240 *     subtractions DSIGMA(I)-DSIGMA(J) unproblematic when cancellation
 241 *     occurs. On binary machines with a guard digit (almost all
 242 *     machines) it does not change DSIGMA(I) at all. On hexadecimal
 243 *     and decimal machines with a guard digit, it slightly
 244 *     changes the bottommost bits of DSIGMA(I). It does not account
 245 *     for hexadecimal or decimal machines without guard digits
 246 *     (we know of none). We use a subroutine call to compute
 247 *     2*DLAMBDA(I) to prevent optimizing compilers from eliminating
 248 *     this code.
 249 *
 250       DO 10 I = 1, K
 251          DSIGMA( I ) = SLAMC3( DSIGMA( I ), DSIGMA( I ) ) - DSIGMA( I )
 252    10 CONTINUE
 253 *
 254 *     Book keeping.
 255 *
 256       IWK1 = 1
 257       IWK2 = IWK1 + K
 258       IWK3 = IWK2 + K
 259       IWK2I = IWK2 - 1
 260       IWK3I = IWK3 - 1
 261 *
 262 *     Normalize Z.
 263 *
 264       RHO = SNRM2( K, Z, 1 )
 265       CALL SLASCL( 'G', 0, 0, RHO, ONE, K, 1, Z, K, INFO )
 266       RHO = RHO*RHO
 267 *
 268 *     Initialize WORK(IWK3).
 269 *
 270       CALL SLASET( 'A', K, 1, ONE, ONE, WORK( IWK3 ), K )
 271 *
 272 *     Compute the updated singular values, the arrays DIFL, DIFR,
 273 *     and the updated Z.
 274 *
 275       DO 40 J = 1, K
 276          CALL SLASD4( K, J, DSIGMA, Z, WORK( IWK1 ), RHO, D( J ),
 277      $                WORK( IWK2 ), INFO )
 278 *
 279 *        If the root finder fails, report the convergence failure.
 280 *
 281          IF( INFO.NE.0 ) THEN
 282             RETURN
 283          END IF
 284          WORK( IWK3I+J ) = WORK( IWK3I+J )*WORK( J )*WORK( IWK2I+J )
 285          DIFL( J ) = -WORK( J )
 286          DIFR( J, 1 ) = -WORK( J+1 )
 287          DO 20 I = 1, J - 1
 288             WORK( IWK3I+I ) = WORK( IWK3I+I )*WORK( I )*
 289      $                        WORK( IWK2I+I ) / ( DSIGMA( I )-
 290      $                        DSIGMA( J ) ) / ( DSIGMA( I )+
 291      $                        DSIGMA( J ) )
 292    20    CONTINUE
 293          DO 30 I = J + 1, K
 294             WORK( IWK3I+I ) = WORK( IWK3I+I )*WORK( I )*
 295      $                        WORK( IWK2I+I ) / ( DSIGMA( I )-
 296      $                        DSIGMA( J ) ) / ( DSIGMA( I )+
 297      $                        DSIGMA( J ) )
 298    30    CONTINUE
 299    40 CONTINUE
 300 *
 301 *     Compute updated Z.
 302 *
 303       DO 50 I = 1, K
 304          Z( I ) = SIGN( SQRT( ABS( WORK( IWK3I+I ) ) ), Z( I ) )
 305    50 CONTINUE
 306 *
 307 *     Update VF and VL.
 308 *
 309       DO 80 J = 1, K
 310          DIFLJ = DIFL( J )
 311          DJ = D( J )
 312          DSIGJ = -DSIGMA( J )
 313          IF( J.LT.K ) THEN
 314             DIFRJ = -DIFR( J, 1 )
 315             DSIGJP = -DSIGMA( J+1 )
 316          END IF
 317          WORK( J ) = -Z( J ) / DIFLJ / ( DSIGMA( J )+DJ )
 318          DO 60 I = 1, J - 1
 319             WORK( I ) = Z( I ) / ( SLAMC3( DSIGMA( I ), DSIGJ )-DIFLJ )
 320      $                   / ( DSIGMA( I )+DJ )
 321    60    CONTINUE
 322          DO 70 I = J + 1, K
 323             WORK( I ) = Z( I ) / ( SLAMC3( DSIGMA( I ), DSIGJP )+DIFRJ )
 324      $                   / ( DSIGMA( I )+DJ )
 325    70    CONTINUE
 326          TEMP = SNRM2( K, WORK, 1 )
 327          WORK( IWK2I+J ) = SDOT( K, WORK, 1, VF, 1 ) / TEMP
 328          WORK( IWK3I+J ) = SDOT( K, WORK, 1, VL, 1 ) / TEMP
 329          IF( ICOMPQ.EQ.1 ) THEN
 330             DIFR( J, 2 ) = TEMP
 331          END IF
 332    80 CONTINUE
 333 *
 334       CALL SCOPY( K, WORK( IWK2 ), 1, VF, 1 )
 335       CALL SCOPY( K, WORK( IWK3 ), 1, VL, 1 )
 336 *
 337       RETURN
 338 *
 339 *     End of SLASD8
 340 *
 341       END
 342