vendor/gonum.org/v1/gonum/lapack/internal/testdata/dlasqtest/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
 178       IF( ( N0-I0-1 ).LE.0 )
 179      $   RETURN
 180 *
 181       DTHRESH = EPS*(SIGMA+TAU)
 182       IF( TAU.LT.DTHRESH*HALF ) TAU = ZERO
 183       IF( TAU.NE.ZERO ) THEN
 184       J4 = 4*I0 + PP - 3
 185       EMIN = Z( J4+4 )
 186       D = Z( J4 ) - TAU
 187       DMIN = D
 188       DMIN1 = -Z( J4 )
 189 *
 190       IF( IEEE ) THEN
 191 *
 192 *        Code for IEEE arithmetic.
 193 *
 194          IF( PP.EQ.0 ) THEN
 195             DO 10 J4 = 4*I0, 4*( N0-3 ), 4
 196                Z( J4-2 ) = D + Z( J4-1 )
 197                TEMP = Z( J4+1 ) / Z( J4-2 )
 198                D = D*TEMP - TAU
 199                DMIN = MIN( DMIN, D )
 200                Z( J4 ) = Z( J4-1 )*TEMP
 201                EMIN = MIN( Z( J4 ), EMIN )
 202    10       CONTINUE
 203          ELSE
 204             DO 20 J4 = 4*I0, 4*( N0-3 ), 4
 205                Z( J4-3 ) = D + Z( J4 )
 206                TEMP = Z( J4+2 ) / Z( J4-3 )
 207                D = D*TEMP - TAU
 208                DMIN = MIN( DMIN, D )
 209                Z( J4-1 ) = Z( J4 )*TEMP
 210                EMIN = MIN( Z( J4-1 ), EMIN )
 211    20       CONTINUE
 212          END IF
 213
 214 *
 215 *        Unroll last two steps.
 216 *
 217          DNM2 = D
 218          DMIN2 = DMIN
 219          J4 = 4*( N0-2 ) - PP
 220          J4P2 = J4 + 2*PP - 1
 221          Z( J4-2 ) = DNM2 + Z( J4P2 )
 222          Z( J4 ) = Z( J4P2+2 )*( Z( J4P2 ) / Z( J4-2 ) )
 223          DNM1 = Z( J4P2+2 )*( DNM2 / Z( J4-2 ) ) - TAU
 224          DMIN = MIN( DMIN, DNM1 )
 225 *
 226          DMIN1 = DMIN
 227          J4 = J4 + 4
 228          J4P2 = J4 + 2*PP - 1
 229          Z( J4-2 ) = DNM1 + Z( J4P2 )
 230          Z( J4 ) = Z( J4P2+2 )*( Z( J4P2 ) / Z( J4-2 ) )
 231          DN = Z( J4P2+2 )*( DNM1 / Z( J4-2 ) ) - TAU
 232          DMIN = MIN( DMIN, DN )
 233 *
 234       ELSE
 235 *
 236 *        Code for non IEEE arithmetic.
 237 *
 238          IF( PP.EQ.0 ) THEN
 239             DO 30 J4 = 4*I0, 4*( N0-3 ), 4
 240                Z( J4-2 ) = D + Z( J4-1 )
 241                IF( D.LT.ZERO ) THEN
 242                   RETURN
 243                ELSE
 244                   Z( J4 ) = Z( J4+1 )*( Z( J4-1 ) / Z( J4-2 ) )
 245                   D = Z( J4+1 )*( D / Z( J4-2 ) ) - TAU
 246                END IF
 247                DMIN = MIN( DMIN, D )
 248                EMIN = MIN( EMIN, Z( J4 ) )
 249    30       CONTINUE
 250          ELSE
 251             DO 40 J4 = 4*I0, 4*( N0-3 ), 4
 252                Z( J4-3 ) = D + Z( J4 )
 253                IF( D.LT.ZERO ) THEN
 254                   RETURN
 255                ELSE
 256                   Z( J4-1 ) = Z( J4+2 )*( Z( J4 ) / Z( J4-3 ) )
 257                   D = Z( J4+2 )*( D / Z( J4-3 ) ) - TAU
 258                END IF
 259                DMIN = MIN( DMIN, D )
 260                EMIN = MIN( EMIN, Z( J4-1 ) )
 261    40       CONTINUE
 262          END IF
 263 *
 264 *        Unroll last two steps.
 265 *
 266          DNM2 = D
 267          DMIN2 = DMIN
 268          J4 = 4*( N0-2 ) - PP
 269          J4P2 = J4 + 2*PP - 1
 270          Z( J4-2 ) = DNM2 + Z( J4P2 )
 271          IF( DNM2.LT.ZERO ) THEN
 272             RETURN
 273          ELSE
 274             Z( J4 ) = Z( J4P2+2 )*( Z( J4P2 ) / Z( J4-2 ) )
 275             DNM1 = Z( J4P2+2 )*( DNM2 / Z( J4-2 ) ) - TAU
 276          END IF
 277          DMIN = MIN( DMIN, DNM1 )
 278 *
 279          DMIN1 = DMIN
 280          J4 = J4 + 4
 281          J4P2 = J4 + 2*PP - 1
 282          Z( J4-2 ) = DNM1 + Z( J4P2 )
 283          IF( DNM1.LT.ZERO ) THEN
 284             RETURN
 285          ELSE
 286             Z( J4 ) = Z( J4P2+2 )*( Z( J4P2 ) / Z( J4-2 ) )
 287             DN = Z( J4P2+2 )*( DNM1 / Z( J4-2 ) ) - TAU
 288          END IF
 289          DMIN = MIN( DMIN, DN )
 290 *
 291       END IF
 292       ELSE
 293 *     This is the version that sets d's to zero if they are small enough
 294          J4 = 4*I0 + PP - 3
 295          EMIN = Z( J4+4 )
 296          D = Z( J4 ) - TAU
 297          DMIN = D
 298          DMIN1 = -Z( J4 )
 299          IF( IEEE ) THEN
 300 *
 301 *     Code for IEEE arithmetic.
 302 *
 303
 304             IF( PP.EQ.0 ) THEN
 305                DO 50 J4 = 4*I0, 4*( N0-3 ), 4
 306                   Z( J4-2 ) = D + Z( J4-1 )
 307                   TEMP = Z( J4+1 ) / Z( J4-2 )
 308                   D = D*TEMP - TAU
 309                   IF( D.LT.DTHRESH ) D = ZERO
 310                   DMIN = MIN( DMIN, D )
 311                   Z( J4 ) = Z( J4-1 )*TEMP
 312                   EMIN = MIN( Z( J4 ), EMIN )
 313  50            CONTINUE
 314             ELSE
 315                DO 60 J4 = 4*I0, 4*( N0-3 ), 4
 316                   Z( J4-3 ) = D + Z( J4 )
 317                   TEMP = Z( J4+2 ) / Z( J4-3 )
 318                   D = D*TEMP - TAU
 319                   IF( D.LT.DTHRESH ) D = ZERO
 320                   DMIN = MIN( DMIN, D )
 321                   Z( J4-1 ) = Z( J4 )*TEMP
 322                   EMIN = MIN( Z( J4-1 ), EMIN )
 323  60            CONTINUE
 324             END IF
 325 *
 326 *     Unroll last two steps.
 327 *
 328             DNM2 = D
 329             DMIN2 = DMIN
 330             J4 = 4*( N0-2 ) - PP
 331             J4P2 = J4 + 2*PP - 1
 332             Z( J4-2 ) = DNM2 + Z( J4P2 )
 333             Z( J4 ) = Z( J4P2+2 )*( Z( J4P2 ) / Z( J4-2 ) )
 334             DNM1 = Z( J4P2+2 )*( DNM2 / Z( J4-2 ) ) - TAU
 335             DMIN = MIN( DMIN, DNM1 )
 336 *
 337             DMIN1 = DMIN
 338             J4 = J4 + 4
 339             J4P2 = J4 + 2*PP - 1
 340             Z( J4-2 ) = DNM1 + Z( J4P2 )
 341             Z( J4 ) = Z( J4P2+2 )*( Z( J4P2 ) / Z( J4-2 ) )
 342             DN = Z( J4P2+2 )*( DNM1 / Z( J4-2 ) ) - TAU
 343             DMIN = MIN( DMIN, DN )
 344 *
 345          ELSE
 346 *
 347 *     Code for non IEEE arithmetic.
 348 *
 349             IF( PP.EQ.0 ) THEN
 350                DO 70 J4 = 4*I0, 4*( N0-3 ), 4
 351                   Z( J4-2 ) = D + Z( J4-1 )
 352                   IF( D.LT.ZERO ) THEN
 353                      RETURN
 354                   ELSE
 355                      Z( J4 ) = Z( J4+1 )*( Z( J4-1 ) / Z( J4-2 ) )
 356                      D = Z( J4+1 )*( D / Z( J4-2 ) ) - TAU
 357                   END IF
 358                   IF( D.LT.DTHRESH) D = ZERO
 359                   DMIN = MIN( DMIN, D )
 360                   EMIN = MIN( EMIN, Z( J4 ) )
 361  70            CONTINUE
 362             ELSE
 363                DO 80 J4 = 4*I0, 4*( N0-3 ), 4
 364                   Z( J4-3 ) = D + Z( J4 )
 365                   IF( D.LT.ZERO ) THEN
 366                      RETURN
 367                   ELSE
 368                      Z( J4-1 ) = Z( J4+2 )*( Z( J4 ) / Z( J4-3 ) )
 369                      D = Z( J4+2 )*( D / Z( J4-3 ) ) - TAU
 370                   END IF
 371                   IF( D.LT.DTHRESH) D = ZERO
 372                   DMIN = MIN( DMIN, D )
 373                   EMIN = MIN( EMIN, Z( J4-1 ) )
 374  80            CONTINUE
 375             END IF
 376 *
 377 *     Unroll last two steps.
 378 *
 379             DNM2 = D
 380             DMIN2 = DMIN
 381             J4 = 4*( N0-2 ) - PP
 382             J4P2 = J4 + 2*PP - 1
 383             Z( J4-2 ) = DNM2 + Z( J4P2 )
 384             IF( DNM2.LT.ZERO ) THEN
 385                RETURN
 386             ELSE
 387                Z( J4 ) = Z( J4P2+2 )*( Z( J4P2 ) / Z( J4-2 ) )
 388                DNM1 = Z( J4P2+2 )*( DNM2 / Z( J4-2 ) ) - TAU
 389             END IF
 390             DMIN = MIN( DMIN, DNM1 )
 391 *
 392             DMIN1 = DMIN
 393             J4 = J4 + 4
 394             J4P2 = J4 + 2*PP - 1
 395             Z( J4-2 ) = DNM1 + Z( J4P2 )
 396             IF( DNM1.LT.ZERO ) THEN
 397                RETURN
 398             ELSE
 399                Z( J4 ) = Z( J4P2+2 )*( Z( J4P2 ) / Z( J4-2 ) )
 400                DN = Z( J4P2+2 )*( DNM1 / Z( J4-2 ) ) - TAU
 401             END IF
 402             DMIN = MIN( DMIN, DN )
 403 *
 404          END IF
 405       END IF
 406 *
 407       Z( J4+2 ) = DN
 408       Z( 4*N0-PP ) = EMIN
 409       RETURN
 410 *
 411 *     End of DLASQ5
 412 *
 413       END