vendor/gonum.org/v1/gonum/lapack/internal/testdata/dlasqtest/dlasq4.f

   1 *> \brief \b DLASQ4 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 DLASQ4 + dependencies
  10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dlasq4.f">
  11 *> [TGZ]</a>
  12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dlasq4.f">
  13 *> [ZIP]</a>
  14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dlasq4.f">
  15 *> [TXT]</a>
  16 *> \endhtmlonly
  17 *
  18 *  Definition:
  19 *  ===========
  20 *
  21 *       SUBROUTINE DLASQ4( 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 *       DOUBLE PRECISION   DMIN, DMIN1, DMIN2, DN, DN1, DN2, G, TAU
  27 *       ..
  28 *       .. Array Arguments ..
  29 *       DOUBLE PRECISION   Z( * )
  30 *       ..
  31 *
  32 *
  33 *> \par Purpose:
  34 *  =============
  35 *>
  36 *> \verbatim
  37 *>
  38 *> DLASQ4 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 DOUBLE PRECISION array, dimension ( 4*N )
  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 DOUBLE PRECISION
  78 *>        Minimum value of d.
  79 *> \endverbatim
  80 *>
  81 *> \param[in] DMIN1
  82 *> \verbatim
  83 *>          DMIN1 is DOUBLE PRECISION
  84 *>        Minimum value of d, excluding D( N0 ).
  85 *> \endverbatim
  86 *>
  87 *> \param[in] DMIN2
  88 *> \verbatim
  89 *>          DMIN2 is DOUBLE PRECISION
  90 *>        Minimum value of d, excluding D( N0 ) and D( N0-1 ).
  91 *> \endverbatim
  92 *>
  93 *> \param[in] DN
  94 *> \verbatim
  95 *>          DN is DOUBLE PRECISION
  96 *>        d(N)
  97 *> \endverbatim
  98 *>
  99 *> \param[in] DN1
 100 *> \verbatim
 101 *>          DN1 is DOUBLE PRECISION
 102 *>        d(N-1)
 103 *> \endverbatim
 104 *>
 105 *> \param[in] DN2
 106 *> \verbatim
 107 *>          DN2 is DOUBLE PRECISION
 108 *>        d(N-2)
 109 *> \endverbatim
 110 *>
 111 *> \param[out] TAU
 112 *> \verbatim
 113 *>          TAU is DOUBLE PRECISION
 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 DLASQ4.
 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 September 2012
 139 *
 140 *> \ingroup auxOTHERcomputational
 141 *
 142 *> \par Further Details:
 143 *  =====================
 144 *>
 145 *> \verbatim
 146 *>
 147 *>  CNST1 = 9/16
 148 *> \endverbatim
 149 *>
 150 *  =====================================================================
 151       SUBROUTINE DLASQ4( I0, N0, Z, PP, N0IN, DMIN, DMIN1, DMIN2, DN,
 152      $                   DN1, DN2, TAU, TTYPE, G )
 153 *
 154 *  -- LAPACK computational routine (version 3.4.2) --
 155 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
 156 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
 157 *     September 2012
 158 *
 159 *     .. Scalar Arguments ..
 160       INTEGER            I0, N0, N0IN, PP, TTYPE
 161       DOUBLE PRECISION   DMIN, DMIN1, DMIN2, DN, DN1, DN2, G, TAU
 162 *     ..
 163 *     .. Array Arguments ..
 164       DOUBLE PRECISION   Z( * )
 165 *     ..
 166 *
 167 *  =====================================================================
 168 *
 169 *     .. Parameters ..
 170       DOUBLE PRECISION   CNST1, CNST2, CNST3
 171       PARAMETER          ( CNST1 = 0.5630D0, CNST2 = 1.010D0,
 172      $                   CNST3 = 1.050D0 )
 173       DOUBLE PRECISION   QURTR, THIRD, HALF, ZERO, ONE, TWO, HUNDRD
 174       PARAMETER          ( QURTR = 0.250D0, THIRD = 0.3330D0,
 175      $                   HALF = 0.50D0, ZERO = 0.0D0, ONE = 1.0D0,
 176      $                   TWO = 2.0D0, HUNDRD = 100.0D0 )
 177 *     ..
 178 *     .. Local Scalars ..
 179       INTEGER            I4, NN, NP
 180       DOUBLE PRECISION   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
 191       IF( DMIN.LE.ZERO ) THEN
 192          TAU = -DMIN
 193          TTYPE = -1
 194          RETURN
 195       END IF
 196 *
 197       NN = 4*N0 + PP
 198       IF( N0IN.EQ.N0 ) THEN
 199 *
 200 *        No eigenvalues deflated.
 201 *
 202          IF( DMIN.EQ.DN .OR. DMIN.EQ.DN1 ) THEN
 203 *
 204             B1 = SQRT( Z( NN-3 ) )*SQRT( Z( NN-5 ) )
 205             B2 = SQRT( Z( NN-7 ) )*SQRT( Z( NN-9 ) )
 206             A2 = Z( NN-7 ) + Z( NN-5 )
 207 *
 208 *           Cases 2 and 3.
 209 *
 210             IF( DMIN.EQ.DN .AND. DMIN1.EQ.DN1 ) THEN
 211
 212                GAP2 = DMIN2 - A2 - DMIN2*QURTR
 213                IF( GAP2.GT.ZERO .AND. GAP2.GT.B2 ) THEN
 214                   GAP1 = A2 - DN - ( B2 / GAP2 )*B2
 215                ELSE
 216                   GAP1 = A2 - DN - ( B1+B2 )
 217                END IF
 218                IF( GAP1.GT.ZERO .AND. GAP1.GT.B1 ) THEN
 219                   S = MAX( DN-( B1 / GAP1 )*B1, HALF*DMIN )
 220                   TTYPE = -2
 221                ELSE
 222                   S = ZERO
 223                   IF( DN.GT.B1 )
 224      $               S = DN - B1
 225                   IF( A2.GT.( B1+B2 ) )
 226      $               S = MIN( S, A2-( B1+B2 ) )
 227                   S = MAX( S, THIRD*DMIN )
 228                   TTYPE = -3
 229                END IF
 230             ELSE
 231 *
 232 *              Case 4.
 233 *
 234                TTYPE = -4
 235                S = QURTR*DMIN
 236                IF( DMIN.EQ.DN ) THEN
 237                   GAM = DN
 238                   A2 = ZERO
 239                   IF( Z( NN-5 ) .GT. Z( NN-7 ) )
 240      $               RETURN
 241                   B2 = Z( NN-5 ) / Z( NN-7 )
 242                   NP = NN - 9
 243                ELSE
 244                   NP = NN - 2*PP
 245                   B2 = Z( NP-2 )
 246                   GAM = DN1
 247                   IF( Z( NP-4 ) .GT. Z( NP-2 ) )
 248      $               RETURN
 249                   A2 = Z( NP-4 ) / Z( NP-2 )
 250                   IF( Z( NN-9 ) .GT. Z( NN-11 ) )
 251      $               RETURN
 252                   B2 = Z( NN-9 ) / Z( NN-11 )
 253                   NP = NN - 13
 254                END IF
 255 *
 256 *              Approximate contribution to norm squared from I < NN-1.
 257 *
 258                A2 = A2 + B2
 259                DO 10 I4 = NP, 4*I0 - 1 + PP, -4
 260                   IF( B2.EQ.ZERO )
 261      $               GO TO 20
 262                   B1 = B2
 263                   IF( Z( I4 ) .GT. Z( I4-2 ) )
 264      $               RETURN
 265                   B2 = B2*( Z( I4 ) / Z( I4-2 ) )
 266                   A2 = A2 + B2
 267                   IF( HUNDRD*MAX( B2, B1 ).LT.A2 .OR. CNST1.LT.A2 )
 268      $               GO TO 20
 269    10          CONTINUE
 270    20          CONTINUE
 271                A2 = CNST3*A2
 272 *
 273 *              Rayleigh quotient residual bound.
 274 *
 275                IF( A2.LT.CNST1 )
 276      $            S = GAM*( ONE-SQRT( A2 ) ) / ( ONE+A2 )
 277             END IF
 278          ELSE IF( DMIN.EQ.DN2 ) THEN
 279 *
 280 *           Case 5.
 281 *
 282             TTYPE = -5
 283             S = QURTR*DMIN
 284 *
 285 *           Compute contribution to norm squared from I > NN-2.
 286 *
 287             NP = NN - 2*PP
 288             B1 = Z( NP-2 )
 289             B2 = Z( NP-6 )
 290             GAM = DN2
 291             IF( Z( NP-8 ).GT.B2 .OR. Z( NP-4 ).GT.B1 )
 292      $         RETURN
 293             A2 = ( Z( NP-8 ) / B2 )*( ONE+Z( NP-4 ) / B1 )
 294 *
 295 *           Approximate contribution to norm squared from I < NN-2.
 296 *
 297             IF( N0-I0.GT.2 ) THEN
 298                B2 = Z( NN-13 ) / Z( NN-15 )
 299                A2 = A2 + B2
 300                DO 30 I4 = NN - 17, 4*I0 - 1 + PP, -4
 301                   IF( B2.EQ.ZERO )
 302      $               GO TO 40
 303                   B1 = B2
 304                   IF( Z( I4 ) .GT. Z( I4-2 ) )
 305      $               RETURN
 306                   B2 = B2*( Z( I4 ) / Z( I4-2 ) )
 307                   A2 = A2 + B2
 308                   IF( HUNDRD*MAX( B2, B1 ).LT.A2 .OR. CNST1.LT.A2 )
 309      $               GO TO 40
 310    30          CONTINUE
 311    40          CONTINUE
 312                A2 = CNST3*A2
 313             END IF
 314 *
 315             IF( A2.LT.CNST1 )
 316      $         S = GAM*( ONE-SQRT( A2 ) ) / ( ONE+A2 )
 317          ELSE
 318 *
 319 *           Case 6, no information to guide us.
 320 *
 321             IF( TTYPE.EQ.-6 ) THEN
 322                G = G + THIRD*( ONE-G )
 323             ELSE IF( TTYPE.EQ.-18 ) THEN
 324                G = QURTR*THIRD
 325             ELSE
 326                G = QURTR
 327             END IF
 328             S = G*DMIN
 329             TTYPE = -6
 330          END IF
 331 *
 332       ELSE IF( N0IN.EQ.( N0+1 ) ) THEN
 333 *
 334 *        One eigenvalue just deflated. Use DMIN1, DN1 for DMIN and DN.
 335 *
 336          IF( DMIN1.EQ.DN1 .AND. DMIN2.EQ.DN2 ) THEN
 337 *
 338 *           Cases 7 and 8.
 339 *
 340             TTYPE = -7
 341             S = THIRD*DMIN1
 342             IF( Z( NN-5 ).GT.Z( NN-7 ) )
 343      $         RETURN
 344             B1 = Z( NN-5 ) / Z( NN-7 )
 345             B2 = B1
 346             IF( B2.EQ.ZERO )
 347      $         GO TO 60
 348             DO 50 I4 = 4*N0 - 9 + PP, 4*I0 - 1 + PP, -4
 349                A2 = B1
 350                IF( Z( I4 ).GT.Z( I4-2 ) )
 351      $            RETURN
 352                B1 = B1*( Z( I4 ) / Z( I4-2 ) )
 353                B2 = B2 + B1
 354                IF( HUNDRD*MAX( B1, A2 ).LT.B2 )
 355      $            GO TO 60
 356    50       CONTINUE
 357    60       CONTINUE
 358             B2 = SQRT( CNST3*B2 )
 359             A2 = DMIN1 / ( ONE+B2**2 )
 360             GAP2 = HALF*DMIN2 - A2
 361             IF( GAP2.GT.ZERO .AND. GAP2.GT.B2*A2 ) THEN
 362                S = MAX( S, A2*( ONE-CNST2*A2*( B2 / GAP2 )*B2 ) )
 363             ELSE
 364                S = MAX( S, A2*( ONE-CNST2*B2 ) )
 365                TTYPE = -8
 366             END IF
 367          ELSE
 368 *
 369 *           Case 9.
 370 *
 371             S = QURTR*DMIN1
 372             IF( DMIN1.EQ.DN1 )
 373      $         S = HALF*DMIN1
 374             TTYPE = -9
 375          END IF
 376 *
 377       ELSE IF( N0IN.EQ.( N0+2 ) ) THEN
 378 *
 379 *        Two eigenvalues deflated. Use DMIN2, DN2 for DMIN and DN.
 380 *
 381 *        Cases 10 and 11.
 382 *
 383          IF( DMIN2.EQ.DN2 .AND. TWO*Z( NN-5 ).LT.Z( NN-7 ) ) THEN
 384             TTYPE = -10
 385             S = THIRD*DMIN2
 386             IF( Z( NN-5 ).GT.Z( NN-7 ) )
 387      $         RETURN
 388             B1 = Z( NN-5 ) / Z( NN-7 )
 389             B2 = B1
 390             IF( B2.EQ.ZERO )
 391      $         GO TO 80
 392             DO 70 I4 = 4*N0 - 9 + PP, 4*I0 - 1 + PP, -4
 393                IF( Z( I4 ).GT.Z( I4-2 ) )
 394      $            RETURN
 395                B1 = B1*( Z( I4 ) / Z( I4-2 ) )
 396                B2 = B2 + B1
 397                IF( HUNDRD*B1.LT.B2 )
 398      $            GO TO 80
 399    70       CONTINUE
 400    80       CONTINUE
 401             B2 = SQRT( CNST3*B2 )
 402             A2 = DMIN2 / ( ONE+B2**2 )
 403             GAP2 = Z( NN-7 ) + Z( NN-9 ) -
 404      $             SQRT( Z( NN-11 ) )*SQRT( Z( NN-9 ) ) - A2
 405             IF( GAP2.GT.ZERO .AND. GAP2.GT.B2*A2 ) THEN
 406                S = MAX( S, A2*( ONE-CNST2*A2*( B2 / GAP2 )*B2 ) )
 407             ELSE
 408                S = MAX( S, A2*( ONE-CNST2*B2 ) )
 409             END IF
 410          ELSE
 411             S = QURTR*DMIN2
 412             TTYPE = -11
 413          END IF
 414       ELSE IF( N0IN.GT.( N0+2 ) ) THEN
 415 *
 416 *        Case 12, more than two eigenvalues deflated. No information.
 417 *
 418          S = ZERO
 419          TTYPE = -12
 420       END IF
 421 *
 422       TAU = S
 423       RETURN
 424 *
 425 *     End of DLASQ4
 426 *
 427       END