SRC/dlasq5.f

   1 *> \brief \b DLASQ5 computes one dqds transform in ping-pong form. Used by sbdsqr and sstegr.
   2 *
   3 *  =========== DOCUMENTATION ===========
   4 *
   5 * Online html documentation available at
   6 *            http://www.netlib.org/lapack/explore-html/
   7 *
   8 *> \htmlonly
   9 *> Download DLASQ5 + dependencies
  10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dlasq5.f">
  11 *> [TGZ]</a>
  12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dlasq5.f">
  13 *> [ZIP]</a>
  14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dlasq5.f">
  15 *> [TXT]</a>
  16 *> \endhtmlonly
  17 *
  18 *  Definition:
  19 *  ===========
  20 *
  21 *       SUBROUTINE DLASQ5( I0, N0, Z, PP, TAU, SIGMA, DMIN, DMIN1, DMIN2, DN,
  22 *                          DNM1, DNM2, IEEE, EPS )
  23 *
  24 *       .. Scalar Arguments ..
  25 *       LOGICAL            IEEE
  26 *       INTEGER            I0, N0, PP
  27 *       DOUBLE PRECISION   DMIN, DMIN1, DMIN2, DN, DNM1, DNM2, TAU, SIGMA, EPS
  28 *       ..
  29 *       .. Array Arguments ..
  30 *       DOUBLE PRECISION   Z( * )
  31 *       ..
  32 *
  33 *
  34 *> \par Purpose:
  35 *  =============
  36 *>
  37 *> \verbatim
  38 *>
  39 *> DLASQ5 computes one dqds transform in ping-pong form, one
  40 *> version for IEEE machines another for non IEEE machines.
  41 *> \endverbatim
  42 *
  43 *  Arguments:
  44 *  ==========
  45 *
  46 *> \param[in] I0
  47 *> \verbatim
  48 *>          I0 is INTEGER
  49 *>        First index.
  50 *> \endverbatim
  51 *>
  52 *> \param[in] N0
  53 *> \verbatim
  54 *>          N0 is INTEGER
  55 *>        Last index.
  56 *> \endverbatim
  57 *>
  58 *> \param[in] Z
  59 *> \verbatim
  60 *>          Z is DOUBLE PRECISION array, dimension ( 4*N )
  61 *>        Z holds the qd array. EMIN is stored in Z(4*N0) to avoid
  62 *>        an extra argument.
  63 *> \endverbatim
  64 *>
  65 *> \param[in] PP
  66 *> \verbatim
  67 *>          PP is INTEGER
  68 *>        PP=0 for ping, PP=1 for pong.
  69 *> \endverbatim
  70 *>
  71 *> \param[in] TAU
  72 *> \verbatim
  73 *>          TAU is DOUBLE PRECISION
  74 *>        This is the shift.
  75 *> \endverbatim
  76 *>
  77 *> \param[in] SIGMA
  78 *> \verbatim
  79 *>          SIGMA is DOUBLE PRECISION
  80 *>        This is the accumulated shift up to this step.
  81 *> \endverbatim
  82 *>
  83 *> \param[out] DMIN
  84 *> \verbatim
  85 *>          DMIN is DOUBLE PRECISION
  86 *>        Minimum value of d.
  87 *> \endverbatim
  88 *>
  89 *> \param[out] DMIN1
  90 *> \verbatim
  91 *>          DMIN1 is DOUBLE PRECISION
  92 *>        Minimum value of d, excluding D( N0 ).
  93 *> \endverbatim
  94 *>
  95 *> \param[out] DMIN2
  96 *> \verbatim
  97 *>          DMIN2 is DOUBLE PRECISION
  98 *>        Minimum value of d, excluding D( N0 ) and D( N0-1 ).
  99 *> \endverbatim
 100 *>
 101 *> \param[out] DN
 102 *> \verbatim
 103 *>          DN is DOUBLE PRECISION
 104 *>        d(N0), the last value of d.
 105 *> \endverbatim
 106 *>
 107 *> \param[out] DNM1
 108 *> \verbatim
 109 *>          DNM1 is DOUBLE PRECISION
 110 *>        d(N0-1).
 111 *> \endverbatim
 112 *>
 113 *> \param[out] DNM2
 114 *> \verbatim
 115 *>          DNM2 is DOUBLE PRECISION
 116 *>        d(N0-2).
 117 *> \endverbatim
 118 *>
 119 *> \param[in] IEEE
 120 *> \verbatim
 121 *>          IEEE is LOGICAL
 122 *>        Flag for IEEE or non IEEE arithmetic.
 123 *> \endverbatim
 124 *
 125 *> \param[in] EPS
 126 *> \verbatim
 127 *>          EPS is DOUBLE PRECISION
 128 *>        This is the value of epsilon used.
 129 *> \endverbatim
 130 *>
 131 *  Authors:
 132 *  ========
 133 *
 134 *> \author Univ. of Tennessee
 135 *> \author Univ. of California Berkeley
 136 *> \author Univ. of Colorado Denver
 137 *> \author NAG Ltd.
 138 *
 139 *> \date September 2012
 140 *
 141 *> \ingroup auxOTHERcomputational
 142 *
 143 *  =====================================================================
 144       SUBROUTINE DLASQ5( I0, N0, Z, PP, TAU, SIGMA, DMIN, DMIN1, DMIN2,
 145      $                   DN, DNM1, DNM2, IEEE, EPS )
 146 *
 147 *  -- LAPACK computational routine (version 3.4.2) --
 148 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
 149 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
 150 *     September 2012
 151 *
 152 *     .. Scalar Arguments ..
 153       LOGICAL            IEEE
 154       INTEGER            I0, N0, PP
 155       DOUBLE PRECISION   DMIN, DMIN1, DMIN2, DN, DNM1, DNM2, TAU,
 156      $                   SIGMA, EPS
 157 *     ..
 158 *     .. Array Arguments ..
 159       DOUBLE PRECISION   Z( * )
 160 *     ..
 161 *
 162 *  =====================================================================
 163 *
 164 *     .. Parameter ..
 165       DOUBLE PRECISION   ZERO, HALF
 166       PARAMETER          ( ZERO = 0.0D0, HALF = 0.5 )
 167 *     ..
 168 *     .. Local Scalars ..
 169       INTEGER            J4, J4P2
 170       DOUBLE PRECISION   D, EMIN, TEMP, DTHRESH
 171 *     ..
 172 *     .. Intrinsic Functions ..
 173       INTRINSIC          MIN
 174 *     ..
 175 *     .. Executable Statements ..
 176 *
 177       IF( ( N0-I0-1 ).LE.0 )
 178      $   RETURN
 179 *
 180       DTHRESH = EPS*(SIGMA+TAU)
 181       IF( TAU.LT.DTHRESH*HALF ) TAU = ZERO
 182       IF( TAU.NE.ZERO ) THEN
 183       J4 = 4*I0 + PP - 3
 184       EMIN = Z( J4+4 )
 185       D = Z( J4 ) - TAU
 186       DMIN = D
 187       DMIN1 = -Z( J4 )
 188 *
 189       IF( IEEE ) THEN
 190 *
 191 *        Code for IEEE arithmetic.
 192 *
 193          IF( PP.EQ.0 ) THEN
 194             DO 10 J4 = 4*I0, 4*( N0-3 ), 4
 195                Z( J4-2 ) = D + Z( J4-1 )
 196                TEMP = Z( J4+1 ) / Z( J4-2 )
 197                D = D*TEMP - TAU
 198                DMIN = MIN( DMIN, D )
 199                Z( J4 ) = Z( J4-1 )*TEMP
 200                EMIN = MIN( Z( J4 ), EMIN )
 201    10       CONTINUE
 202          ELSE
 203             DO 20 J4 = 4*I0, 4*( N0-3 ), 4
 204                Z( J4-3 ) = D + Z( J4 )
 205                TEMP = Z( J4+2 ) / Z( J4-3 )
 206                D = D*TEMP - TAU
 207                DMIN = MIN( DMIN, D )
 208                Z( J4-1 ) = Z( J4 )*TEMP
 209                EMIN = MIN( Z( J4-1 ), EMIN )
 210    20       CONTINUE
 211          END IF
 212 *
 213 *        Unroll last two steps.
 214 *
 215          DNM2 = D
 216          DMIN2 = DMIN
 217          J4 = 4*( N0-2 ) - PP
 218          J4P2 = J4 + 2*PP - 1
 219          Z( J4-2 ) = DNM2 + Z( J4P2 )
 220          Z( J4 ) = Z( J4P2+2 )*( Z( J4P2 ) / Z( J4-2 ) )
 221          DNM1 = Z( J4P2+2 )*( DNM2 / Z( J4-2 ) ) - TAU
 222          DMIN = MIN( DMIN, DNM1 )
 223 *
 224          DMIN1 = DMIN
 225          J4 = J4 + 4
 226          J4P2 = J4 + 2*PP - 1
 227          Z( J4-2 ) = DNM1 + Z( J4P2 )
 228          Z( J4 ) = Z( J4P2+2 )*( Z( J4P2 ) / Z( J4-2 ) )
 229          DN = Z( J4P2+2 )*( DNM1 / Z( J4-2 ) ) - TAU
 230          DMIN = MIN( DMIN, DN )
 231 *
 232       ELSE
 233 *
 234 *        Code for non IEEE arithmetic.
 235 *
 236          IF( PP.EQ.0 ) THEN
 237             DO 30 J4 = 4*I0, 4*( N0-3 ), 4
 238                Z( J4-2 ) = D + Z( J4-1 )
 239                IF( D.LT.ZERO ) THEN
 240                   RETURN
 241                ELSE
 242                   Z( J4 ) = Z( J4+1 )*( Z( J4-1 ) / Z( J4-2 ) )
 243                   D = Z( J4+1 )*( D / Z( J4-2 ) ) - TAU
 244                END IF
 245                DMIN = MIN( DMIN, D )
 246                EMIN = MIN( EMIN, Z( J4 ) )
 247    30       CONTINUE
 248          ELSE
 249             DO 40 J4 = 4*I0, 4*( N0-3 ), 4
 250                Z( J4-3 ) = D + Z( J4 )
 251                IF( D.LT.ZERO ) THEN
 252                   RETURN
 253                ELSE
 254                   Z( J4-1 ) = Z( J4+2 )*( Z( J4 ) / Z( J4-3 ) )
 255                   D = Z( J4+2 )*( D / Z( J4-3 ) ) - TAU
 256                END IF
 257                DMIN = MIN( DMIN, D )
 258                EMIN = MIN( EMIN, Z( J4-1 ) )
 259    40       CONTINUE
 260          END IF
 261 *
 262 *        Unroll last two steps.
 263 *
 264          DNM2 = D
 265          DMIN2 = DMIN
 266          J4 = 4*( N0-2 ) - PP
 267          J4P2 = J4 + 2*PP - 1
 268          Z( J4-2 ) = DNM2 + Z( J4P2 )
 269          IF( DNM2.LT.ZERO ) THEN
 270             RETURN
 271          ELSE
 272             Z( J4 ) = Z( J4P2+2 )*( Z( J4P2 ) / Z( J4-2 ) )
 273             DNM1 = Z( J4P2+2 )*( DNM2 / Z( J4-2 ) ) - TAU
 274          END IF
 275          DMIN = MIN( DMIN, DNM1 )
 276 *
 277          DMIN1 = DMIN
 278          J4 = J4 + 4
 279          J4P2 = J4 + 2*PP - 1
 280          Z( J4-2 ) = DNM1 + Z( J4P2 )
 281          IF( DNM1.LT.ZERO ) THEN
 282             RETURN
 283          ELSE
 284             Z( J4 ) = Z( J4P2+2 )*( Z( J4P2 ) / Z( J4-2 ) )
 285             DN = Z( J4P2+2 )*( DNM1 / Z( J4-2 ) ) - TAU
 286          END IF
 287          DMIN = MIN( DMIN, DN )
 288 *
 289       END IF
 290       ELSE
 291 *     This is the version that sets d's to zero if they are small enough
 292          J4 = 4*I0 + PP - 3
 293          EMIN = Z( J4+4 )
 294          D = Z( J4 ) - TAU
 295          DMIN = D
 296          DMIN1 = -Z( J4 )
 297          IF( IEEE ) THEN
 298 *
 299 *     Code for IEEE arithmetic.
 300 *
 301             IF( PP.EQ.0 ) THEN
 302                DO 50 J4 = 4*I0, 4*( N0-3 ), 4
 303                   Z( J4-2 ) = D + Z( J4-1 )
 304                   TEMP = Z( J4+1 ) / Z( J4-2 )
 305                   D = D*TEMP - TAU
 306                   IF( D.LT.DTHRESH ) D = ZERO
 307                   DMIN = MIN( DMIN, D )
 308                   Z( J4 ) = Z( J4-1 )*TEMP
 309                   EMIN = MIN( Z( J4 ), EMIN )
 310  50            CONTINUE
 311             ELSE
 312                DO 60 J4 = 4*I0, 4*( N0-3 ), 4
 313                   Z( J4-3 ) = D + Z( J4 )
 314                   TEMP = Z( J4+2 ) / Z( J4-3 )
 315                   D = D*TEMP - TAU
 316                   IF( D.LT.DTHRESH ) D = ZERO
 317                   DMIN = MIN( DMIN, D )
 318                   Z( J4-1 ) = Z( J4 )*TEMP
 319                   EMIN = MIN( Z( J4-1 ), EMIN )
 320  60            CONTINUE
 321             END IF
 322 *
 323 *     Unroll last two steps.
 324 *
 325             DNM2 = D
 326             DMIN2 = DMIN
 327             J4 = 4*( N0-2 ) - PP
 328             J4P2 = J4 + 2*PP - 1
 329             Z( J4-2 ) = DNM2 + Z( J4P2 )
 330             Z( J4 ) = Z( J4P2+2 )*( Z( J4P2 ) / Z( J4-2 ) )
 331             DNM1 = Z( J4P2+2 )*( DNM2 / Z( J4-2 ) ) - TAU
 332             DMIN = MIN( DMIN, DNM1 )
 333 *
 334             DMIN1 = DMIN
 335             J4 = J4 + 4
 336             J4P2 = J4 + 2*PP - 1
 337             Z( J4-2 ) = DNM1 + Z( J4P2 )
 338             Z( J4 ) = Z( J4P2+2 )*( Z( J4P2 ) / Z( J4-2 ) )
 339             DN = Z( J4P2+2 )*( DNM1 / Z( J4-2 ) ) - TAU
 340             DMIN = MIN( DMIN, DN )
 341 *
 342          ELSE
 343 *
 344 *     Code for non IEEE arithmetic.
 345 *
 346             IF( PP.EQ.0 ) THEN
 347                DO 70 J4 = 4*I0, 4*( N0-3 ), 4
 348                   Z( J4-2 ) = D + Z( J4-1 )
 349                   IF( D.LT.ZERO ) THEN
 350                      RETURN
 351                   ELSE
 352                      Z( J4 ) = Z( J4+1 )*( Z( J4-1 ) / Z( J4-2 ) )
 353                      D = Z( J4+1 )*( D / Z( J4-2 ) ) - TAU
 354                   END IF
 355                   IF( D.LT.DTHRESH) D = ZERO
 356                   DMIN = MIN( DMIN, D )
 357                   EMIN = MIN( EMIN, Z( J4 ) )
 358  70            CONTINUE
 359             ELSE
 360                DO 80 J4 = 4*I0, 4*( N0-3 ), 4
 361                   Z( J4-3 ) = D + Z( J4 )
 362                   IF( D.LT.ZERO ) THEN
 363                      RETURN
 364                   ELSE
 365                      Z( J4-1 ) = Z( J4+2 )*( Z( J4 ) / Z( J4-3 ) )
 366                      D = Z( J4+2 )*( D / Z( J4-3 ) ) - TAU
 367                   END IF
 368                   IF( D.LT.DTHRESH) D = ZERO
 369                   DMIN = MIN( DMIN, D )
 370                   EMIN = MIN( EMIN, Z( J4-1 ) )
 371  80            CONTINUE
 372             END IF
 373 *
 374 *     Unroll last two steps.
 375 *
 376             DNM2 = D
 377             DMIN2 = DMIN
 378             J4 = 4*( N0-2 ) - PP
 379             J4P2 = J4 + 2*PP - 1
 380             Z( J4-2 ) = DNM2 + Z( J4P2 )
 381             IF( DNM2.LT.ZERO ) THEN
 382                RETURN
 383             ELSE
 384                Z( J4 ) = Z( J4P2+2 )*( Z( J4P2 ) / Z( J4-2 ) )
 385                DNM1 = Z( J4P2+2 )*( DNM2 / Z( J4-2 ) ) - TAU
 386             END IF
 387             DMIN = MIN( DMIN, DNM1 )
 388 *
 389             DMIN1 = DMIN
 390             J4 = J4 + 4
 391             J4P2 = J4 + 2*PP - 1
 392             Z( J4-2 ) = DNM1 + Z( J4P2 )
 393             IF( DNM1.LT.ZERO ) THEN
 394                RETURN
 395             ELSE
 396                Z( J4 ) = Z( J4P2+2 )*( Z( J4P2 ) / Z( J4-2 ) )
 397                DN = Z( J4P2+2 )*( DNM1 / Z( J4-2 ) ) - TAU
 398             END IF
 399             DMIN = MIN( DMIN, DN )
 400 *
 401          END IF
 402       END IF
 403 *
 404       Z( J4+2 ) = DN
 405       Z( 4*N0-PP ) = EMIN
 406       RETURN
 407 *
 408 *     End of DLASQ5
 409 *
 410       END