SRC/slasq4.f

   1 *> \brief \b SLASQ4 computes an approximation to the smallest eigenvalue using values of d from the previous transform. 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 SLASQ4 + dependencies
  10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/slasq4.f">
  11 *> [TGZ]</a>
  12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/slasq4.f">
  13 *> [ZIP]</a>
  14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/slasq4.f">
  15 *> [TXT]</a>
  16 *> \endhtmlonly
  17 *
  18 *  Definition:
  19 *  ===========
  20 *
  21 *       SUBROUTINE SLASQ4( I0, N0, Z, PP, N0IN, DMIN, DMIN1, DMIN2, DN,
  22 *                          DN1, DN2, TAU, TTYPE, G )
  23 *
  24 *       .. Scalar Arguments ..
  25 *       INTEGER            I0, N0, N0IN, PP, TTYPE
  26 *       REAL               DMIN, DMIN1, DMIN2, DN, DN1, DN2, G, TAU
  27 *       ..
  28 *       .. Array Arguments ..
  29 *       REAL               Z( * )
  30 *       ..
  31 *
  32 *
  33 *> \par Purpose:
  34 *  =============
  35 *>
  36 *> \verbatim
  37 *>
  38 *> SLASQ4 computes an approximation TAU to the smallest eigenvalue
  39 *> using values of d from the previous transform.
  40 *> \endverbatim
  41 *
  42 *  Arguments:
  43 *  ==========
  44 *
  45 *> \param[in] I0
  46 *> \verbatim
  47 *>          I0 is INTEGER
  48 *>        First index.
  49 *> \endverbatim
  50 *>
  51 *> \param[in] N0
  52 *> \verbatim
  53 *>          N0 is INTEGER
  54 *>        Last index.
  55 *> \endverbatim
  56 *>
  57 *> \param[in] Z
  58 *> \verbatim
  59 *>          Z is REAL array, dimension ( 4*N0 )
  60 *>        Z holds the qd array.
  61 *> \endverbatim
  62 *>
  63 *> \param[in] PP
  64 *> \verbatim
  65 *>          PP is INTEGER
  66 *>        PP=0 for ping, PP=1 for pong.
  67 *> \endverbatim
  68 *>
  69 *> \param[in] N0IN
  70 *> \verbatim
  71 *>          N0IN is INTEGER
  72 *>        The value of N0 at start of EIGTEST.
  73 *> \endverbatim
  74 *>
  75 *> \param[in] DMIN
  76 *> \verbatim
  77 *>          DMIN is REAL
  78 *>        Minimum value of d.
  79 *> \endverbatim
  80 *>
  81 *> \param[in] DMIN1
  82 *> \verbatim
  83 *>          DMIN1 is REAL
  84 *>        Minimum value of d, excluding D( N0 ).
  85 *> \endverbatim
  86 *>
  87 *> \param[in] DMIN2
  88 *> \verbatim
  89 *>          DMIN2 is REAL
  90 *>        Minimum value of d, excluding D( N0 ) and D( N0-1 ).
  91 *> \endverbatim
  92 *>
  93 *> \param[in] DN
  94 *> \verbatim
  95 *>          DN is REAL
  96 *>        d(N)
  97 *> \endverbatim
  98 *>
  99 *> \param[in] DN1
 100 *> \verbatim
 101 *>          DN1 is REAL
 102 *>        d(N-1)
 103 *> \endverbatim
 104 *>
 105 *> \param[in] DN2
 106 *> \verbatim
 107 *>          DN2 is REAL
 108 *>        d(N-2)
 109 *> \endverbatim
 110 *>
 111 *> \param[out] TAU
 112 *> \verbatim
 113 *>          TAU is REAL
 114 *>        This is the shift.
 115 *> \endverbatim
 116 *>
 117 *> \param[out] TTYPE
 118 *> \verbatim
 119 *>          TTYPE is INTEGER
 120 *>        Shift type.
 121 *> \endverbatim
 122 *>
 123 *> \param[in,out] G
 124 *> \verbatim
 125 *>          G is REAL
 126 *>        G is passed as an argument in order to save its value between
 127 *>        calls to SLASQ4.
 128 *> \endverbatim
 129 *
 130 *  Authors:
 131 *  ========
 132 *
 133 *> \author Univ. of Tennessee
 134 *> \author Univ. of California Berkeley
 135 *> \author Univ. of Colorado Denver
 136 *> \author NAG Ltd.
 137 *
 138 *> \date June 2016
 139 *
 140 *> \ingroup auxOTHERcomputational
 141 *
 142 *> \par Further Details:
 143 *  =====================
 144 *>
 145 *> \verbatim
 146 *>
 147 *>  CNST1 = 9/16
 148 *> \endverbatim
 149 *>
 150 *  =====================================================================
 151       SUBROUTINE SLASQ4( I0, N0, Z, PP, N0IN, DMIN, DMIN1, DMIN2, DN,
 152      $                   DN1, DN2, TAU, TTYPE, G )
 153 *
 154 *  -- LAPACK computational routine (version 3.6.1) --
 155 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
 156 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
 157 *     June 2016
 158 *
 159 *     .. Scalar Arguments ..
 160       INTEGER            I0, N0, N0IN, PP, TTYPE
 161       REAL               DMIN, DMIN1, DMIN2, DN, DN1, DN2, G, TAU
 162 *     ..
 163 *     .. Array Arguments ..
 164       REAL               Z( * )
 165 *     ..
 166 *
 167 *  =====================================================================
 168 *
 169 *     .. Parameters ..
 170       REAL               CNST1, CNST2, CNST3
 171       PARAMETER          ( CNST1 = 0.5630E0, CNST2 = 1.010E0,
 172      $                   CNST3 = 1.050E0 )
 173       REAL               QURTR, THIRD, HALF, ZERO, ONE, TWO, HUNDRD
 174       PARAMETER          ( QURTR = 0.250E0, THIRD = 0.3330E0,
 175      $                   HALF = 0.50E0, ZERO = 0.0E0, ONE = 1.0E0,
 176      $                   TWO = 2.0E0, HUNDRD = 100.0E0 )
 177 *     ..
 178 *     .. Local Scalars ..
 179       INTEGER            I4, NN, NP
 180       REAL               A2, B1, B2, GAM, GAP1, GAP2, S
 181 *     ..
 182 *     .. Intrinsic Functions ..
 183       INTRINSIC          MAX, MIN, SQRT
 184 *     ..
 185 *     .. Executable Statements ..
 186 *
 187 *     A negative DMIN forces the shift to take that absolute value
 188 *     TTYPE records the type of shift.
 189 *
 190       IF( DMIN.LE.ZERO ) THEN
 191          TAU = -DMIN
 192          TTYPE = -1
 193          RETURN
 194       END IF
 195 *
 196       NN = 4*N0 + PP
 197       IF( N0IN.EQ.N0 ) THEN
 198 *
 199 *        No eigenvalues deflated.
 200 *
 201          IF( DMIN.EQ.DN .OR. DMIN.EQ.DN1 ) THEN
 202 *
 203             B1 = SQRT( Z( NN-3 ) )*SQRT( Z( NN-5 ) )
 204             B2 = SQRT( Z( NN-7 ) )*SQRT( Z( NN-9 ) )
 205             A2 = Z( NN-7 ) + Z( NN-5 )
 206 *
 207 *           Cases 2 and 3.
 208 *
 209             IF( DMIN.EQ.DN .AND. DMIN1.EQ.DN1 ) THEN
 210                GAP2 = DMIN2 - A2 - DMIN2*QURTR
 211                IF( GAP2.GT.ZERO .AND. GAP2.GT.B2 ) THEN
 212                   GAP1 = A2 - DN - ( B2 / GAP2 )*B2
 213                ELSE
 214                   GAP1 = A2 - DN - ( B1+B2 )
 215                END IF
 216                IF( GAP1.GT.ZERO .AND. GAP1.GT.B1 ) THEN
 217                   S = MAX( DN-( B1 / GAP1 )*B1, HALF*DMIN )
 218                   TTYPE = -2
 219                ELSE
 220                   S = ZERO
 221                   IF( DN.GT.B1 )
 222      $               S = DN - B1
 223                   IF( A2.GT.( B1+B2 ) )
 224      $               S = MIN( S, A2-( B1+B2 ) )
 225                   S = MAX( S, THIRD*DMIN )
 226                   TTYPE = -3
 227                END IF
 228             ELSE
 229 *
 230 *              Case 4.
 231 *
 232                TTYPE = -4
 233                S = QURTR*DMIN
 234                IF( DMIN.EQ.DN ) THEN
 235                   GAM = DN
 236                   A2 = ZERO
 237                   IF( Z( NN-5 ) .GT. Z( NN-7 ) )
 238      $               RETURN
 239                   B2 = Z( NN-5 ) / Z( NN-7 )
 240                   NP = NN - 9
 241                ELSE
 242                   NP = NN - 2*PP
 243                   B2 = Z( NP-2 )
 244                   GAM = DN1
 245                   IF( Z( NP-4 ) .GT. Z( NP-2 ) )
 246      $               RETURN
 247                   A2 = Z( NP-4 ) / Z( NP-2 )
 248                   IF( Z( NN-9 ) .GT. Z( NN-11 ) )
 249      $               RETURN
 250                   B2 = Z( NN-9 ) / Z( NN-11 )
 251                   NP = NN - 13
 252                END IF
 253 *
 254 *              Approximate contribution to norm squared from I < NN-1.
 255 *
 256                A2 = A2 + B2
 257                DO 10 I4 = NP, 4*I0 - 1 + PP, -4
 258                   IF( B2.EQ.ZERO )
 259      $               GO TO 20
 260                   B1 = B2
 261                   IF( Z( I4 ) .GT. Z( I4-2 ) )
 262      $               RETURN
 263                   B2 = B2*( Z( I4 ) / Z( I4-2 ) )
 264                   A2 = A2 + B2
 265                   IF( HUNDRD*MAX( B2, B1 ).LT.A2 .OR. CNST1.LT.A2 )
 266      $               GO TO 20
 267    10          CONTINUE
 268    20          CONTINUE
 269                A2 = CNST3*A2
 270 *
 271 *              Rayleigh quotient residual bound.
 272 *
 273                IF( A2.LT.CNST1 )
 274      $            S = GAM*( ONE-SQRT( A2 ) ) / ( ONE+A2 )
 275             END IF
 276          ELSE IF( DMIN.EQ.DN2 ) THEN
 277 *
 278 *           Case 5.
 279 *
 280             TTYPE = -5
 281             S = QURTR*DMIN
 282 *
 283 *           Compute contribution to norm squared from I > NN-2.
 284 *
 285             NP = NN - 2*PP
 286             B1 = Z( NP-2 )
 287             B2 = Z( NP-6 )
 288             GAM = DN2
 289             IF( Z( NP-8 ).GT.B2 .OR. Z( NP-4 ).GT.B1 )
 290      $         RETURN
 291             A2 = ( Z( NP-8 ) / B2 )*( ONE+Z( NP-4 ) / B1 )
 292 *
 293 *           Approximate contribution to norm squared from I < NN-2.
 294 *
 295             IF( N0-I0.GT.2 ) THEN
 296                B2 = Z( NN-13 ) / Z( NN-15 )
 297                A2 = A2 + B2
 298                DO 30 I4 = NN - 17, 4*I0 - 1 + PP, -4
 299                   IF( B2.EQ.ZERO )
 300      $               GO TO 40
 301                   B1 = B2
 302                   IF( Z( I4 ) .GT. Z( I4-2 ) )
 303      $               RETURN
 304                   B2 = B2*( Z( I4 ) / Z( I4-2 ) )
 305                   A2 = A2 + B2
 306                   IF( HUNDRD*MAX( B2, B1 ).LT.A2 .OR. CNST1.LT.A2 )
 307      $               GO TO 40
 308    30          CONTINUE
 309    40          CONTINUE
 310                A2 = CNST3*A2
 311             END IF
 312 *
 313             IF( A2.LT.CNST1 )
 314      $         S = GAM*( ONE-SQRT( A2 ) ) / ( ONE+A2 )
 315          ELSE
 316 *
 317 *           Case 6, no information to guide us.
 318 *
 319             IF( TTYPE.EQ.-6 ) THEN
 320                G = G + THIRD*( ONE-G )
 321             ELSE IF( TTYPE.EQ.-18 ) THEN
 322                G = QURTR*THIRD
 323             ELSE
 324                G = QURTR
 325             END IF
 326             S = G*DMIN
 327             TTYPE = -6
 328          END IF
 329 *
 330       ELSE IF( N0IN.EQ.( N0+1 ) ) THEN
 331 *
 332 *        One eigenvalue just deflated. Use DMIN1, DN1 for DMIN and DN.
 333 *
 334          IF( DMIN1.EQ.DN1 .AND. DMIN2.EQ.DN2 ) THEN
 335 *
 336 *           Cases 7 and 8.
 337 *
 338             TTYPE = -7
 339             S = THIRD*DMIN1
 340             IF( Z( NN-5 ).GT.Z( NN-7 ) )
 341      $         RETURN
 342             B1 = Z( NN-5 ) / Z( NN-7 )
 343             B2 = B1
 344             IF( B2.EQ.ZERO )
 345      $         GO TO 60
 346             DO 50 I4 = 4*N0 - 9 + PP, 4*I0 - 1 + PP, -4
 347                A2 = B1
 348                IF( Z( I4 ).GT.Z( I4-2 ) )
 349      $            RETURN
 350                B1 = B1*( Z( I4 ) / Z( I4-2 ) )
 351                B2 = B2 + B1
 352                IF( HUNDRD*MAX( B1, A2 ).LT.B2 )
 353      $            GO TO 60
 354    50       CONTINUE
 355    60       CONTINUE
 356             B2 = SQRT( CNST3*B2 )
 357             A2 = DMIN1 / ( ONE+B2**2 )
 358             GAP2 = HALF*DMIN2 - A2
 359             IF( GAP2.GT.ZERO .AND. GAP2.GT.B2*A2 ) THEN
 360                S = MAX( S, A2*( ONE-CNST2*A2*( B2 / GAP2 )*B2 ) )
 361             ELSE
 362                S = MAX( S, A2*( ONE-CNST2*B2 ) )
 363                TTYPE = -8
 364             END IF
 365          ELSE
 366 *
 367 *           Case 9.
 368 *
 369             S = QURTR*DMIN1
 370             IF( DMIN1.EQ.DN1 )
 371      $         S = HALF*DMIN1
 372             TTYPE = -9
 373          END IF
 374 *
 375       ELSE IF( N0IN.EQ.( N0+2 ) ) THEN
 376 *
 377 *        Two eigenvalues deflated. Use DMIN2, DN2 for DMIN and DN.
 378 *
 379 *        Cases 10 and 11.
 380 *
 381          IF( DMIN2.EQ.DN2 .AND. TWO*Z( NN-5 ).LT.Z( NN-7 ) ) THEN
 382             TTYPE = -10
 383             S = THIRD*DMIN2
 384             IF( Z( NN-5 ).GT.Z( NN-7 ) )
 385      $         RETURN
 386             B1 = Z( NN-5 ) / Z( NN-7 )
 387             B2 = B1
 388             IF( B2.EQ.ZERO )
 389      $         GO TO 80
 390             DO 70 I4 = 4*N0 - 9 + PP, 4*I0 - 1 + PP, -4
 391                IF( Z( I4 ).GT.Z( I4-2 ) )
 392      $            RETURN
 393                B1 = B1*( Z( I4 ) / Z( I4-2 ) )
 394                B2 = B2 + B1
 395                IF( HUNDRD*B1.LT.B2 )
 396      $            GO TO 80
 397    70       CONTINUE
 398    80       CONTINUE
 399             B2 = SQRT( CNST3*B2 )
 400             A2 = DMIN2 / ( ONE+B2**2 )
 401             GAP2 = Z( NN-7 ) + Z( NN-9 ) -
 402      $             SQRT( Z( NN-11 ) )*SQRT( Z( NN-9 ) ) - A2
 403             IF( GAP2.GT.ZERO .AND. GAP2.GT.B2*A2 ) THEN
 404                S = MAX( S, A2*( ONE-CNST2*A2*( B2 / GAP2 )*B2 ) )
 405             ELSE
 406                S = MAX( S, A2*( ONE-CNST2*B2 ) )
 407             END IF
 408          ELSE
 409             S = QURTR*DMIN2
 410             TTYPE = -11
 411          END IF
 412       ELSE IF( N0IN.GT.( N0+2 ) ) THEN
 413 *
 414 *        Case 12, more than two eigenvalues deflated. No information.
 415 *
 416          S = ZERO
 417          TTYPE = -12
 418       END IF
 419 *
 420       TAU = S
 421       RETURN
 422 *
 423 *     End of SLASQ4
 424 *
 425       END