SRC/slasq5.f

   1 *> \brief <b> SLASQ5 computes one dqds transform in ping-pong form. Used by sbdsqr and sstegr. </b>
   2 *
   3 *  =========== DOCUMENTATION ===========
   4 *
   5 * Online html documentation available at
   6 *            http://www.netlib.org/lapack/explore-html/
   7 *
   8 *> \htmlonly
   9 *> Download SLASQ5 + dependencies
  10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/slasq5.f">
  11 *> [TGZ]</a>
  12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/slasq5.f">
  13 *> [ZIP]</a>
  14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/slasq5.f">
  15 *> [TXT]</a>
  16 *> \endhtmlonly
  17 *
  18 *  Definition:
  19 *  ===========
  20 *
  21 *       SUBROUTINE SLASQ5( 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 *       REAL               EPS, DMIN, DMIN1, DMIN2, DN, DNM1, DNM2, SIGMA, TAU
  28 *       ..
  29 *       .. Array Arguments ..
  30 *       REAL               Z( * )
  31 *       ..
  32 *
  33 *
  34 *> \par Purpose:
  35 *  =============
  36 *>
  37 *> \verbatim
  38 *>
  39 *> SLASQ5 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 REAL 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 REAL
  74 *>        This is the shift.
  75 *> \endverbatim
  76 *>
  77 *> \param[in] SIGMA
  78 *> \verbatim
  79 *>          SIGMA is REAL
  80 *>        This is the accumulated shift up to this step.
  81 *> \endverbatim
  82 *>
  83 *> \param[out] DMIN
  84 *> \verbatim
  85 *>          DMIN is REAL
  86 *>        Minimum value of d.
  87 *> \endverbatim
  88 *>
  89 *> \param[out] DMIN1
  90 *> \verbatim
  91 *>          DMIN1 is REAL
  92 *>        Minimum value of d, excluding D( N0 ).
  93 *> \endverbatim
  94 *>
  95 *> \param[out] DMIN2
  96 *> \verbatim
  97 *>          DMIN2 is REAL
  98 *>        Minimum value of d, excluding D( N0 ) and D( N0-1 ).
  99 *> \endverbatim
 100 *>
 101 *> \param[out] DN
 102 *> \verbatim
 103 *>          DN is REAL
 104 *>        d(N0), the last value of d.
 105 *> \endverbatim
 106 *>
 107 *> \param[out] DNM1
 108 *> \verbatim
 109 *>          DNM1 is REAL
 110 *>        d(N0-1).
 111 *> \endverbatim
 112 *>
 113 *> \param[out] DNM2
 114 *> \verbatim
 115 *>          DNM2 is REAL
 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 REAL
 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 November 2015
 140 *
 141 *> \ingroup auxOTHERcomputational
 142 *
 143 *  =====================================================================
 144       SUBROUTINE SLASQ5( I0, N0, Z, PP, TAU, SIGMA, DMIN, DMIN1, DMIN2,
 145      $                   DN, DNM1, DNM2, IEEE, EPS )
 146 *
 147 *  -- LAPACK computational routine (version 3.6.0) --
 148 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
 149 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
 150 *     November 2015
 151 *
 152 *     .. Scalar Arguments ..
 153       LOGICAL            IEEE
 154       INTEGER            I0, N0, PP
 155       REAL               DMIN, DMIN1, DMIN2, DN, DNM1, DNM2, TAU,
 156      $                   SIGMA, EPS
 157 *     ..
 158 *     .. Array Arguments ..
 159       REAL               Z( * )
 160 *     ..
 161 *
 162 *  =====================================================================
 163 *
 164 *     .. Parameter ..
 165       REAL               ZERO, HALF
 166       PARAMETER          ( ZERO = 0.0E0, HALF = 0.5 )
 167 *     ..
 168 *     .. Local Scalars ..
 169       INTEGER            J4, J4P2
 170       REAL               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 *
 291       ELSE
 292 *     This is the version that sets d's to zero if they are small enough
 293          J4 = 4*I0 + PP - 3
 294          EMIN = Z( J4+4 )
 295          D = Z( J4 ) - TAU
 296          DMIN = D
 297          DMIN1 = -Z( J4 )
 298          IF( IEEE ) THEN
 299 *
 300 *     Code for IEEE arithmetic.
 301 *
 302             IF( PP.EQ.0 ) THEN
 303                DO 50 J4 = 4*I0, 4*( N0-3 ), 4
 304                   Z( J4-2 ) = D + Z( J4-1 )
 305                   TEMP = Z( J4+1 ) / Z( J4-2 )
 306                   D = D*TEMP - TAU
 307                   IF( D.LT.DTHRESH ) D = ZERO
 308                   DMIN = MIN( DMIN, D )
 309                   Z( J4 ) = Z( J4-1 )*TEMP
 310                   EMIN = MIN( Z( J4 ), EMIN )
 311  50            CONTINUE
 312             ELSE
 313                DO 60 J4 = 4*I0, 4*( N0-3 ), 4
 314                   Z( J4-3 ) = D + Z( J4 )
 315                   TEMP = Z( J4+2 ) / Z( J4-3 )
 316                   D = D*TEMP - TAU
 317                   IF( D.LT.DTHRESH ) D = ZERO
 318                   DMIN = MIN( DMIN, D )
 319                   Z( J4-1 ) = Z( J4 )*TEMP
 320                   EMIN = MIN( Z( J4-1 ), EMIN )
 321  60            CONTINUE
 322             END IF
 323 *
 324 *     Unroll last two steps.
 325 *
 326             DNM2 = D
 327             DMIN2 = DMIN
 328             J4 = 4*( N0-2 ) - PP
 329             J4P2 = J4 + 2*PP - 1
 330             Z( J4-2 ) = DNM2 + Z( J4P2 )
 331             Z( J4 ) = Z( J4P2+2 )*( Z( J4P2 ) / Z( J4-2 ) )
 332             DNM1 = Z( J4P2+2 )*( DNM2 / Z( J4-2 ) ) - TAU
 333             DMIN = MIN( DMIN, DNM1 )
 334 *
 335             DMIN1 = DMIN
 336             J4 = J4 + 4
 337             J4P2 = J4 + 2*PP - 1
 338             Z( J4-2 ) = DNM1 + Z( J4P2 )
 339             Z( J4 ) = Z( J4P2+2 )*( Z( J4P2 ) / Z( J4-2 ) )
 340             DN = Z( J4P2+2 )*( DNM1 / Z( J4-2 ) ) - TAU
 341             DMIN = MIN( DMIN, DN )
 342 *
 343          ELSE
 344 *
 345 *     Code for non IEEE arithmetic.
 346 *
 347             IF( PP.EQ.0 ) THEN
 348                DO 70 J4 = 4*I0, 4*( N0-3 ), 4
 349                   Z( J4-2 ) = D + Z( J4-1 )
 350                   IF( D.LT.ZERO ) THEN
 351                      RETURN
 352                   ELSE
 353                      Z( J4 ) = Z( J4+1 )*( Z( J4-1 ) / Z( J4-2 ) )
 354                      D = Z( J4+1 )*( D / Z( J4-2 ) ) - TAU
 355                   END IF
 356                   IF( D.LT.DTHRESH ) D = ZERO
 357                   DMIN = MIN( DMIN, D )
 358                   EMIN = MIN( EMIN, Z( J4 ) )
 359  70            CONTINUE
 360             ELSE
 361                DO 80 J4 = 4*I0, 4*( N0-3 ), 4
 362                   Z( J4-3 ) = D + Z( J4 )
 363                   IF( D.LT.ZERO ) THEN
 364                      RETURN
 365                   ELSE
 366                      Z( J4-1 ) = Z( J4+2 )*( Z( J4 ) / Z( J4-3 ) )
 367                      D = Z( J4+2 )*( D / Z( J4-3 ) ) - TAU
 368                   END IF
 369                   IF( D.LT.DTHRESH ) D = ZERO
 370                   DMIN = MIN( DMIN, D )
 371                   EMIN = MIN( EMIN, Z( J4-1 ) )
 372  80            CONTINUE
 373             END IF
 374 *
 375 *     Unroll last two steps.
 376 *
 377             DNM2 = D
 378             DMIN2 = DMIN
 379             J4 = 4*( N0-2 ) - PP
 380             J4P2 = J4 + 2*PP - 1
 381             Z( J4-2 ) = DNM2 + Z( J4P2 )
 382             IF( DNM2.LT.ZERO ) THEN
 383                RETURN
 384             ELSE
 385                Z( J4 ) = Z( J4P2+2 )*( Z( J4P2 ) / Z( J4-2 ) )
 386                DNM1 = Z( J4P2+2 )*( DNM2 / Z( J4-2 ) ) - TAU
 387             END IF
 388             DMIN = MIN( DMIN, DNM1 )
 389 *
 390             DMIN1 = DMIN
 391             J4 = J4 + 4
 392             J4P2 = J4 + 2*PP - 1
 393             Z( J4-2 ) = DNM1 + Z( J4P2 )
 394             IF( DNM1.LT.ZERO ) THEN
 395                RETURN
 396             ELSE
 397                Z( J4 ) = Z( J4P2+2 )*( Z( J4P2 ) / Z( J4-2 ) )
 398                DN = Z( J4P2+2 )*( DNM1 / Z( J4-2 ) ) - TAU
 399             END IF
 400             DMIN = MIN( DMIN, DN )
 401 *
 402          END IF
 403 *
 404       END IF
 405       Z( J4+2 ) = DN
 406       Z( 4*N0-PP ) = EMIN
 407       RETURN
 408 *
 409 *     End of SLASQ5
 410 *
 411       END