SRC/slasq5.f

   1 *> \brief \b SLASQ5 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 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 *>          SIGMA is REAL
  79 *>        This is the accumulated shift up to this step.
  80 *> \endverbatim
  81 *>
  82 *> \param[out] DMIN
  83 *> \verbatim
  84 *>          DMIN is REAL
  85 *>        Minimum value of d.
  86 *> \endverbatim
  87 *>
  88 *> \param[out] DMIN1
  89 *> \verbatim
  90 *>          DMIN1 is REAL
  91 *>        Minimum value of d, excluding D( N0 ).
  92 *> \endverbatim
  93 *>
  94 *> \param[out] DMIN2
  95 *> \verbatim
  96 *>          DMIN2 is REAL
  97 *>        Minimum value of d, excluding D( N0 ) and D( N0-1 ).
  98 *> \endverbatim
  99 *>
 100 *> \param[out] DN
 101 *> \verbatim
 102 *>          DN is REAL
 103 *>        d(N0), the last value of d.
 104 *> \endverbatim
 105 *>
 106 *> \param[out] DNM1
 107 *> \verbatim
 108 *>          DNM1 is REAL
 109 *>        d(N0-1).
 110 *> \endverbatim
 111 *>
 112 *> \param[out] DNM2
 113 *> \verbatim
 114 *>          DNM2 is REAL
 115 *>        d(N0-2).
 116 *> \endverbatim
 117 *>
 118 *> \param[in] IEEE
 119 *> \verbatim
 120 *>          IEEE is LOGICAL
 121 *>        Flag for IEEE or non IEEE arithmetic.
 122 *> \endverbatim
 123 *>
 124 *> \param[in] EPS
 125 *> \verbatim
 126 *>         EPS is REAL
 127 *>        This is the value of epsilon used.
 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 November 2015
 139 *
 140 *> \ingroup auxOTHERcomputational
 141 *
 142 *  =====================================================================
 143       SUBROUTINE SLASQ5( I0, N0, Z, PP, TAU, SIGMA, DMIN, DMIN1, DMIN2,
 144      $                   DN, DNM1, DNM2, IEEE, EPS )
 145 *
 146 *  -- LAPACK computational routine (version 3.6.0) --
 147 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
 148 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
 149 *     November 2015
 150 *
 151 *     .. Scalar Arguments ..
 152       LOGICAL            IEEE
 153       INTEGER            I0, N0, PP
 154       REAL               DMIN, DMIN1, DMIN2, DN, DNM1, DNM2, TAU,
 155      $                   SIGMA, EPS
 156 *     ..
 157 *     .. Array Arguments ..
 158       REAL               Z( * )
 159 *     ..
 160 *
 161 *  =====================================================================
 162 *
 163 *     .. Parameter ..
 164       REAL               ZERO, HALF
 165       PARAMETER          ( ZERO = 0.0E0, HALF = 0.5 )
 166 *     ..
 167 *     .. Local Scalars ..
 168       INTEGER            J4, J4P2
 169       REAL               D, EMIN, TEMP, DTHRESH
 170 *     ..
 171 *     .. Intrinsic Functions ..
 172       INTRINSIC          MIN
 173 *     ..
 174 *     .. Executable Statements ..
 175 *
 176       IF( ( N0-I0-1 ).LE.0 )
 177      $   RETURN
 178 *
 179       DTHRESH = EPS*(SIGMA+TAU)
 180       IF( TAU.LT.DTHRESH*HALF ) TAU = ZERO
 181       IF( TAU.NE.ZERO ) THEN
 182          J4 = 4*I0 + PP - 3
 183          EMIN = Z( J4+4 )
 184          D = Z( J4 ) - TAU
 185          DMIN = D
 186          DMIN1 = -Z( J4 )
 187 *
 188          IF( IEEE ) THEN
 189 *
 190 *     Code for IEEE arithmetic.
 191 *
 192             IF( PP.EQ.0 ) THEN
 193                DO 10 J4 = 4*I0, 4*( N0-3 ), 4
 194                   Z( J4-2 ) = D + Z( J4-1 )
 195                   TEMP = Z( J4+1 ) / Z( J4-2 )
 196                   D = D*TEMP - TAU
 197                   DMIN = MIN( DMIN, D )
 198                   Z( J4 ) = Z( J4-1 )*TEMP
 199                   EMIN = MIN( Z( J4 ), EMIN )
 200  10            CONTINUE
 201             ELSE
 202                DO 20 J4 = 4*I0, 4*( N0-3 ), 4
 203                   Z( J4-3 ) = D + Z( J4 )
 204                   TEMP = Z( J4+2 ) / Z( J4-3 )
 205                   D = D*TEMP - TAU
 206                   DMIN = MIN( DMIN, D )
 207                   Z( J4-1 ) = Z( J4 )*TEMP
 208                   EMIN = MIN( Z( J4-1 ), EMIN )
 209  20            CONTINUE
 210             END IF
 211 *
 212 *     Unroll last two steps.
 213 *
 214             DNM2 = D
 215             DMIN2 = DMIN
 216             J4 = 4*( N0-2 ) - PP
 217             J4P2 = J4 + 2*PP - 1
 218             Z( J4-2 ) = DNM2 + Z( J4P2 )
 219             Z( J4 ) = Z( J4P2+2 )*( Z( J4P2 ) / Z( J4-2 ) )
 220             DNM1 = Z( J4P2+2 )*( DNM2 / Z( J4-2 ) ) - TAU
 221             DMIN = MIN( DMIN, DNM1 )
 222 *
 223             DMIN1 = DMIN
 224             J4 = J4 + 4
 225             J4P2 = J4 + 2*PP - 1
 226             Z( J4-2 ) = DNM1 + Z( J4P2 )
 227             Z( J4 ) = Z( J4P2+2 )*( Z( J4P2 ) / Z( J4-2 ) )
 228             DN = Z( J4P2+2 )*( DNM1 / Z( J4-2 ) ) - TAU
 229             DMIN = MIN( DMIN, DN )
 230 *
 231          ELSE
 232 *
 233 *     Code for non IEEE arithmetic.
 234 *
 235             IF( PP.EQ.0 ) THEN
 236                DO 30 J4 = 4*I0, 4*( N0-3 ), 4
 237                   Z( J4-2 ) = D + Z( J4-1 )
 238                   IF( D.LT.ZERO ) THEN
 239                      RETURN
 240                   ELSE
 241                      Z( J4 ) = Z( J4+1 )*( Z( J4-1 ) / Z( J4-2 ) )
 242                      D = Z( J4+1 )*( D / Z( J4-2 ) ) - TAU
 243                   END IF
 244                   DMIN = MIN( DMIN, D )
 245                   EMIN = MIN( EMIN, Z( J4 ) )
 246  30            CONTINUE
 247             ELSE
 248                DO 40 J4 = 4*I0, 4*( N0-3 ), 4
 249                   Z( J4-3 ) = D + Z( J4 )
 250                   IF( D.LT.ZERO ) THEN
 251                      RETURN
 252                   ELSE
 253                      Z( J4-1 ) = Z( J4+2 )*( Z( J4 ) / Z( J4-3 ) )
 254                      D = Z( J4+2 )*( D / Z( J4-3 ) ) - TAU
 255                   END IF
 256                   DMIN = MIN( DMIN, D )
 257                   EMIN = MIN( EMIN, Z( J4-1 ) )
 258  40            CONTINUE
 259             END IF
 260 *
 261 *     Unroll last two steps.
 262 *
 263             DNM2 = D
 264             DMIN2 = DMIN
 265             J4 = 4*( N0-2 ) - PP
 266             J4P2 = J4 + 2*PP - 1
 267             Z( J4-2 ) = DNM2 + Z( J4P2 )
 268             IF( DNM2.LT.ZERO ) THEN
 269                RETURN
 270             ELSE
 271                Z( J4 ) = Z( J4P2+2 )*( Z( J4P2 ) / Z( J4-2 ) )
 272                DNM1 = Z( J4P2+2 )*( DNM2 / Z( J4-2 ) ) - TAU
 273             END IF
 274             DMIN = MIN( DMIN, DNM1 )
 275 *
 276             DMIN1 = DMIN
 277             J4 = J4 + 4
 278             J4P2 = J4 + 2*PP - 1
 279             Z( J4-2 ) = DNM1 + Z( J4P2 )
 280             IF( DNM1.LT.ZERO ) THEN
 281                RETURN
 282             ELSE
 283                Z( J4 ) = Z( J4P2+2 )*( Z( J4P2 ) / Z( J4-2 ) )
 284                DN = Z( J4P2+2 )*( DNM1 / Z( J4-2 ) ) - TAU
 285             END IF
 286             DMIN = MIN( DMIN, DN )
 287 *
 288          END IF
 289 *
 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 *
 403       END IF
 404       Z( J4+2 ) = DN
 405       Z( 4*N0-PP ) = EMIN
 406       RETURN
 407 *
 408 *     End of SLASQ5
 409 *
 410       END