SRC/slasq3.f

   1 *> \brief \b SLASQ3 checks for deflation, computes a shift and calls dqds. Used by sbdsqr.
   2 *
   3 *  =========== DOCUMENTATION ===========
   4 *
   5 * Online html documentation available at
   6 *            http://www.netlib.org/lapack/explore-html/
   7 *
   8 *> \htmlonly
   9 *> Download SLASQ3 + dependencies
  10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/slasq3.f">
  11 *> [TGZ]</a>
  12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/slasq3.f">
  13 *> [ZIP]</a>
  14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/slasq3.f">
  15 *> [TXT]</a>
  16 *> \endhtmlonly
  17 *
  18 *  Definition:
  19 *  ===========
  20 *
  21 *       SUBROUTINE SLASQ3( I0, N0, Z, PP, DMIN, SIGMA, DESIG, QMAX, NFAIL,
  22 *                          ITER, NDIV, IEEE, TTYPE, DMIN1, DMIN2, DN, DN1,
  23 *                          DN2, G, TAU )
  24 *
  25 *       .. Scalar Arguments ..
  26 *       LOGICAL            IEEE
  27 *       INTEGER            I0, ITER, N0, NDIV, NFAIL, PP
  28 *       REAL               DESIG, DMIN, DMIN1, DMIN2, DN, DN1, DN2, G,
  29 *      $                   QMAX, SIGMA, TAU
  30 *       ..
  31 *       .. Array Arguments ..
  32 *       REAL               Z( * )
  33 *       ..
  34 *
  35 *
  36 *> \par Purpose:
  37 *  =============
  38 *>
  39 *> \verbatim
  40 *>
  41 *> SLASQ3 checks for deflation, computes a shift (TAU) and calls dqds.
  42 *> In case of failure it changes shifts, and tries again until output
  43 *> is positive.
  44 *> \endverbatim
  45 *
  46 *  Arguments:
  47 *  ==========
  48 *
  49 *> \param[in] I0
  50 *> \verbatim
  51 *>          I0 is INTEGER
  52 *>         First index.
  53 *> \endverbatim
  54 *>
  55 *> \param[in,out] N0
  56 *> \verbatim
  57 *>          N0 is INTEGER
  58 *>         Last index.
  59 *> \endverbatim
  60 *>
  61 *> \param[in,out] Z
  62 *> \verbatim
  63 *>          Z is REAL array, dimension ( 4*N0 )
  64 *>         Z holds the qd array.
  65 *> \endverbatim
  66 *>
  67 *> \param[in,out] PP
  68 *> \verbatim
  69 *>          PP is INTEGER
  70 *>         PP=0 for ping, PP=1 for pong.
  71 *>         PP=2 indicates that flipping was applied to the Z array
  72 *>         and that the initial tests for deflation should not be
  73 *>         performed.
  74 *> \endverbatim
  75 *>
  76 *> \param[out] DMIN
  77 *> \verbatim
  78 *>          DMIN is REAL
  79 *>         Minimum value of d.
  80 *> \endverbatim
  81 *>
  82 *> \param[out] SIGMA
  83 *> \verbatim
  84 *>          SIGMA is REAL
  85 *>         Sum of shifts used in current segment.
  86 *> \endverbatim
  87 *>
  88 *> \param[in,out] DESIG
  89 *> \verbatim
  90 *>          DESIG is REAL
  91 *>         Lower order part of SIGMA
  92 *> \endverbatim
  93 *>
  94 *> \param[in] QMAX
  95 *> \verbatim
  96 *>          QMAX is REAL
  97 *>         Maximum value of q.
  98 *> \endverbatim
  99 *>
 100 *> \param[in,out] NFAIL
 101 *> \verbatim
 102 *>          NFAIL is INTEGER
 103 *>         Increment NFAIL by 1 each time the shift was too big.
 104 *> \endverbatim
 105 *>
 106 *> \param[in,out] ITER
 107 *> \verbatim
 108 *>          ITER is INTEGER
 109 *>         Increment ITER by 1 for each iteration.
 110 *> \endverbatim
 111 *>
 112 *> \param[in,out] NDIV
 113 *> \verbatim
 114 *>          NDIV is INTEGER
 115 *>         Increment NDIV by 1 for each division.
 116 *> \endverbatim
 117 *>
 118 *> \param[in] IEEE
 119 *> \verbatim
 120 *>          IEEE is LOGICAL
 121 *>         Flag for IEEE or non IEEE arithmetic (passed to SLASQ5).
 122 *> \endverbatim
 123 *>
 124 *> \param[in,out] TTYPE
 125 *> \verbatim
 126 *>          TTYPE is INTEGER
 127 *>         Shift type.
 128 *> \endverbatim
 129 *>
 130 *> \param[in,out] DMIN1
 131 *> \verbatim
 132 *>          DMIN1 is REAL
 133 *> \endverbatim
 134 *>
 135 *> \param[in,out] DMIN2
 136 *> \verbatim
 137 *>          DMIN2 is REAL
 138 *> \endverbatim
 139 *>
 140 *> \param[in,out] DN
 141 *> \verbatim
 142 *>          DN is REAL
 143 *> \endverbatim
 144 *>
 145 *> \param[in,out] DN1
 146 *> \verbatim
 147 *>          DN1 is REAL
 148 *> \endverbatim
 149 *>
 150 *> \param[in,out] DN2
 151 *> \verbatim
 152 *>          DN2 is REAL
 153 *> \endverbatim
 154 *>
 155 *> \param[in,out] G
 156 *> \verbatim
 157 *>          G is REAL
 158 *> \endverbatim
 159 *>
 160 *> \param[in,out] TAU
 161 *> \verbatim
 162 *>          TAU is REAL
 163 *>
 164 *>         These are passed as arguments in order to save their values
 165 *>         between calls to SLASQ3.
 166 *> \endverbatim
 167 *
 168 *  Authors:
 169 *  ========
 170 *
 171 *> \author Univ. of Tennessee
 172 *> \author Univ. of California Berkeley
 173 *> \author Univ. of Colorado Denver
 174 *> \author NAG Ltd.
 175 *
 176 *> \date June 2016
 177 *
 178 *> \ingroup auxOTHERcomputational
 179 *
 180 *  =====================================================================
 181       SUBROUTINE SLASQ3( I0, N0, Z, PP, DMIN, SIGMA, DESIG, QMAX, NFAIL,
 182      $                   ITER, NDIV, IEEE, TTYPE, DMIN1, DMIN2, DN, DN1,
 183      $                   DN2, G, TAU )
 184 *
 185 *  -- LAPACK computational routine (version 3.6.1) --
 186 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
 187 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
 188 *     June 2016
 189 *
 190 *     .. Scalar Arguments ..
 191       LOGICAL            IEEE
 192       INTEGER            I0, ITER, N0, NDIV, NFAIL, PP
 193       REAL               DESIG, DMIN, DMIN1, DMIN2, DN, DN1, DN2, G,
 194      $                   QMAX, SIGMA, TAU
 195 *     ..
 196 *     .. Array Arguments ..
 197       REAL               Z( * )
 198 *     ..
 199 *
 200 *  =====================================================================
 201 *
 202 *     .. Parameters ..
 203       REAL               CBIAS
 204       PARAMETER          ( CBIAS = 1.50E0 )
 205       REAL               ZERO, QURTR, HALF, ONE, TWO, HUNDRD
 206       PARAMETER          ( ZERO = 0.0E0, QURTR = 0.250E0, HALF = 0.5E0,
 207      $                     ONE = 1.0E0, TWO = 2.0E0, HUNDRD = 100.0E0 )
 208 *     ..
 209 *     .. Local Scalars ..
 210       INTEGER            IPN4, J4, N0IN, NN, TTYPE
 211       REAL               EPS, S, T, TEMP, TOL, TOL2
 212 *     ..
 213 *     .. External Subroutines ..
 214       EXTERNAL           SLASQ4, SLASQ5, SLASQ6
 215 *     ..
 216 *     .. External Function ..
 217       REAL               SLAMCH
 218       LOGICAL            SISNAN
 219       EXTERNAL           SISNAN, SLAMCH
 220 *     ..
 221 *     .. Intrinsic Functions ..
 222       INTRINSIC          ABS, MAX, MIN, SQRT
 223 *     ..
 224 *     .. Executable Statements ..
 225 *
 226       N0IN = N0
 227       EPS = SLAMCH( 'Precision' )
 228       TOL = EPS*HUNDRD
 229       TOL2 = TOL**2
 230 *
 231 *     Check for deflation.
 232 *
 233    10 CONTINUE
 234 *
 235       IF( N0.LT.I0 )
 236      $   RETURN
 237       IF( N0.EQ.I0 )
 238      $   GO TO 20
 239       NN = 4*N0 + PP
 240       IF( N0.EQ.( I0+1 ) )
 241      $   GO TO 40
 242 *
 243 *     Check whether E(N0-1) is negligible, 1 eigenvalue.
 244 *
 245       IF( Z( NN-5 ).GT.TOL2*( SIGMA+Z( NN-3 ) ) .AND.
 246      $    Z( NN-2*PP-4 ).GT.TOL2*Z( NN-7 ) )
 247      $   GO TO 30
 248 *
 249    20 CONTINUE
 250 *
 251       Z( 4*N0-3 ) = Z( 4*N0+PP-3 ) + SIGMA
 252       N0 = N0 - 1
 253       GO TO 10
 254 *
 255 *     Check  whether E(N0-2) is negligible, 2 eigenvalues.
 256 *
 257    30 CONTINUE
 258 *
 259       IF( Z( NN-9 ).GT.TOL2*SIGMA .AND.
 260      $    Z( NN-2*PP-8 ).GT.TOL2*Z( NN-11 ) )
 261      $   GO TO 50
 262 *
 263    40 CONTINUE
 264 *
 265       IF( Z( NN-3 ).GT.Z( NN-7 ) ) THEN
 266          S = Z( NN-3 )
 267          Z( NN-3 ) = Z( NN-7 )
 268          Z( NN-7 ) = S
 269       END IF
 270       T = HALF*( ( Z( NN-7 )-Z( NN-3 ) )+Z( NN-5 ) )
 271       IF( Z( NN-5 ).GT.Z( NN-3 )*TOL2.AND.T.NE.ZERO ) THEN
 272          S = Z( NN-3 )*( Z( NN-5 ) / T )
 273          IF( S.LE.T ) THEN
 274             S = Z( NN-3 )*( Z( NN-5 ) /
 275      $          ( T*( ONE+SQRT( ONE+S / T ) ) ) )
 276          ELSE
 277             S = Z( NN-3 )*( Z( NN-5 ) / ( T+SQRT( T )*SQRT( T+S ) ) )
 278          END IF
 279          T = Z( NN-7 ) + ( S+Z( NN-5 ) )
 280          Z( NN-3 ) = Z( NN-3 )*( Z( NN-7 ) / T )
 281          Z( NN-7 ) = T
 282       END IF
 283       Z( 4*N0-7 ) = Z( NN-7 ) + SIGMA
 284       Z( 4*N0-3 ) = Z( NN-3 ) + SIGMA
 285       N0 = N0 - 2
 286       GO TO 10
 287 *
 288    50 CONTINUE
 289       IF( PP.EQ.2 )
 290      $   PP = 0
 291 *
 292 *     Reverse the qd-array, if warranted.
 293 *
 294       IF( DMIN.LE.ZERO .OR. N0.LT.N0IN ) THEN
 295          IF( CBIAS*Z( 4*I0+PP-3 ).LT.Z( 4*N0+PP-3 ) ) THEN
 296             IPN4 = 4*( I0+N0 )
 297             DO 60 J4 = 4*I0, 2*( I0+N0-1 ), 4
 298                TEMP = Z( J4-3 )
 299                Z( J4-3 ) = Z( IPN4-J4-3 )
 300                Z( IPN4-J4-3 ) = TEMP
 301                TEMP = Z( J4-2 )
 302                Z( J4-2 ) = Z( IPN4-J4-2 )
 303                Z( IPN4-J4-2 ) = TEMP
 304                TEMP = Z( J4-1 )
 305                Z( J4-1 ) = Z( IPN4-J4-5 )
 306                Z( IPN4-J4-5 ) = TEMP
 307                TEMP = Z( J4 )
 308                Z( J4 ) = Z( IPN4-J4-4 )
 309                Z( IPN4-J4-4 ) = TEMP
 310    60       CONTINUE
 311             IF( N0-I0.LE.4 ) THEN
 312                Z( 4*N0+PP-1 ) = Z( 4*I0+PP-1 )
 313                Z( 4*N0-PP ) = Z( 4*I0-PP )
 314             END IF
 315             DMIN2 = MIN( DMIN2, Z( 4*N0+PP-1 ) )
 316             Z( 4*N0+PP-1 ) = MIN( Z( 4*N0+PP-1 ), Z( 4*I0+PP-1 ),
 317      $                            Z( 4*I0+PP+3 ) )
 318             Z( 4*N0-PP ) = MIN( Z( 4*N0-PP ), Z( 4*I0-PP ),
 319      $                          Z( 4*I0-PP+4 ) )
 320             QMAX = MAX( QMAX, Z( 4*I0+PP-3 ), Z( 4*I0+PP+1 ) )
 321             DMIN = -ZERO
 322          END IF
 323       END IF
 324 *
 325 *     Choose a shift.
 326 *
 327       CALL SLASQ4( I0, N0, Z, PP, N0IN, DMIN, DMIN1, DMIN2, DN, DN1,
 328      $             DN2, TAU, TTYPE, G )
 329 *
 330 *     Call dqds until DMIN > 0.
 331 *
 332    70 CONTINUE
 333 *
 334       CALL SLASQ5( I0, N0, Z, PP, TAU, SIGMA, DMIN, DMIN1, DMIN2, DN,
 335      $             DN1, DN2, IEEE, EPS )
 336 *
 337       NDIV = NDIV + ( N0-I0+2 )
 338       ITER = ITER + 1
 339 *
 340 *     Check status.
 341 *
 342       IF( DMIN.GE.ZERO .AND. DMIN1.GE.ZERO ) THEN
 343 *
 344 *        Success.
 345 *
 346          GO TO 90
 347 *
 348       ELSE IF( DMIN.LT.ZERO .AND. DMIN1.GT.ZERO .AND.
 349      $         Z( 4*( N0-1 )-PP ).LT.TOL*( SIGMA+DN1 ) .AND.
 350      $         ABS( DN ).LT.TOL*SIGMA ) THEN
 351 *
 352 *        Convergence hidden by negative DN.
 353 *
 354          Z( 4*( N0-1 )-PP+2 ) = ZERO
 355          DMIN = ZERO
 356          GO TO 90
 357       ELSE IF( DMIN.LT.ZERO ) THEN
 358 *
 359 *        TAU too big. Select new TAU and try again.
 360 *
 361          NFAIL = NFAIL + 1
 362          IF( TTYPE.LT.-22 ) THEN
 363 *
 364 *           Failed twice. Play it safe.
 365 *
 366             TAU = ZERO
 367          ELSE IF( DMIN1.GT.ZERO ) THEN
 368 *
 369 *           Late failure. Gives excellent shift.
 370 *
 371             TAU = ( TAU+DMIN )*( ONE-TWO*EPS )
 372             TTYPE = TTYPE - 11
 373          ELSE
 374 *
 375 *           Early failure. Divide by 4.
 376 *
 377             TAU = QURTR*TAU
 378             TTYPE = TTYPE - 12
 379          END IF
 380          GO TO 70
 381       ELSE IF( SISNAN( DMIN ) ) THEN
 382 *
 383 *        NaN.
 384 *
 385          IF( TAU.EQ.ZERO ) THEN
 386             GO TO 80
 387          ELSE
 388             TAU = ZERO
 389             GO TO 70
 390          END IF
 391       ELSE
 392 *
 393 *        Possible underflow. Play it safe.
 394 *
 395          GO TO 80
 396       END IF
 397 *
 398 *     Risk of underflow.
 399 *
 400    80 CONTINUE
 401       CALL SLASQ6( I0, N0, Z, PP, DMIN, DMIN1, DMIN2, DN, DN1, DN2 )
 402       NDIV = NDIV + ( N0-I0+2 )
 403       ITER = ITER + 1
 404       TAU = ZERO
 405 *
 406    90 CONTINUE
 407       IF( TAU.LT.SIGMA ) THEN
 408          DESIG = DESIG + TAU
 409          T = SIGMA + DESIG
 410          DESIG = DESIG - ( T-SIGMA )
 411       ELSE
 412          T = SIGMA + TAU
 413          DESIG = SIGMA - ( T-TAU ) + DESIG
 414       END IF
 415       SIGMA = T
 416 *
 417       RETURN
 418 *
 419 *     End of SLASQ3
 420 *
 421       END