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

   1 *> \brief \b DLASQ3 checks for deflation, computes a shift and calls dqds. 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 DLASQ3 + dependencies
  10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dlasq3.f">
  11 *> [TGZ]</a>
  12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dlasq3.f">
  13 *> [ZIP]</a>
  14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dlasq3.f">
  15 *> [TXT]</a>
  16 *> \endhtmlonly
  17 *
  18 *  Definition:
  19 *  ===========
  20 *
  21 *       SUBROUTINE DLASQ3( I0, N0, Z, PP, DMIN, SIGMA, DESIG, QMAX, NFAIL,
  22 *                          ITER, NDIV, IEEE, TTYPE, DMIN1, DMIN2, DN, DN1,
  23 *                          DN2, G, TAU )
  24 *
  25 *       .. Scalar Arguments ..
  26 *       LOGICAL            IEEE
  27 *       INTEGER            I0, ITER, N0, NDIV, NFAIL, PP
  28 *       DOUBLE PRECISION   DESIG, DMIN, DMIN1, DMIN2, DN, DN1, DN2, G,
  29 *      $                   QMAX, SIGMA, TAU
  30 *       ..
  31 *       .. Array Arguments ..
  32 *       DOUBLE PRECISION   Z( * )
  33 *       ..
  34 *
  35 *
  36 *> \par Purpose:
  37 *  =============
  38 *>
  39 *> \verbatim
  40 *>
  41 *> DLASQ3 checks for deflation, computes a shift (TAU) and calls dqds.
  42 *> In case of failure it changes shifts, and tries again until output
  43 *> is positive.
  44 *> \endverbatim
  45 *
  46 *  Arguments:
  47 *  ==========
  48 *
  49 *> \param[in] I0
  50 *> \verbatim
  51 *>          I0 is INTEGER
  52 *>         First index.
  53 *> \endverbatim
  54 *>
  55 *> \param[in,out] N0
  56 *> \verbatim
  57 *>          N0 is INTEGER
  58 *>         Last index.
  59 *> \endverbatim
  60 *>
  61 *> \param[in] Z
  62 *> \verbatim
  63 *>          Z is DOUBLE PRECISION array, dimension ( 4*N )
  64 *>         Z holds the qd array.
  65 *> \endverbatim
  66 *>
  67 *> \param[in,out] PP
  68 *> \verbatim
  69 *>          PP is INTEGER
  70 *>         PP=0 for ping, PP=1 for pong.
  71 *>         PP=2 indicates that flipping was applied to the Z array
  72 *>         and that the initial tests for deflation should not be
  73 *>         performed.
  74 *> \endverbatim
  75 *>
  76 *> \param[out] DMIN
  77 *> \verbatim
  78 *>          DMIN is DOUBLE PRECISION
  79 *>         Minimum value of d.
  80 *> \endverbatim
  81 *>
  82 *> \param[out] SIGMA
  83 *> \verbatim
  84 *>          SIGMA is DOUBLE PRECISION
  85 *>         Sum of shifts used in current segment.
  86 *> \endverbatim
  87 *>
  88 *> \param[in,out] DESIG
  89 *> \verbatim
  90 *>          DESIG is DOUBLE PRECISION
  91 *>         Lower order part of SIGMA
  92 *> \endverbatim
  93 *>
  94 *> \param[in] QMAX
  95 *> \verbatim
  96 *>          QMAX is DOUBLE PRECISION
  97 *>         Maximum value of q.
  98 *> \endverbatim
  99 *>
 100 *> \param[out] NFAIL
 101 *> \verbatim
 102 *>          NFAIL is INTEGER
 103 *>         Number of times shift was too big.
 104 *> \endverbatim
 105 *>
 106 *> \param[out] ITER
 107 *> \verbatim
 108 *>          ITER is INTEGER
 109 *>         Number of iterations.
 110 *> \endverbatim
 111 *>
 112 *> \param[out] NDIV
 113 *> \verbatim
 114 *>          NDIV is INTEGER
 115 *>         Number of divisions.
 116 *> \endverbatim
 117 *>
 118 *> \param[in] IEEE
 119 *> \verbatim
 120 *>          IEEE is LOGICAL
 121 *>         Flag for IEEE or non IEEE arithmetic (passed to DLASQ5).
 122 *> \endverbatim
 123 *>
 124 *> \param[in,out] TTYPE
 125 *> \verbatim
 126 *>          TTYPE is INTEGER
 127 *>         Shift type.
 128 *> \endverbatim
 129 *>
 130 *> \param[in,out] DMIN1
 131 *> \verbatim
 132 *>          DMIN1 is DOUBLE PRECISION
 133 *> \endverbatim
 134 *>
 135 *> \param[in,out] DMIN2
 136 *> \verbatim
 137 *>          DMIN2 is DOUBLE PRECISION
 138 *> \endverbatim
 139 *>
 140 *> \param[in,out] DN
 141 *> \verbatim
 142 *>          DN is DOUBLE PRECISION
 143 *> \endverbatim
 144 *>
 145 *> \param[in,out] DN1
 146 *> \verbatim
 147 *>          DN1 is DOUBLE PRECISION
 148 *> \endverbatim
 149 *>
 150 *> \param[in,out] DN2
 151 *> \verbatim
 152 *>          DN2 is DOUBLE PRECISION
 153 *> \endverbatim
 154 *>
 155 *> \param[in,out] G
 156 *> \verbatim
 157 *>          G is DOUBLE PRECISION
 158 *> \endverbatim
 159 *>
 160 *> \param[in,out] TAU
 161 *> \verbatim
 162 *>          TAU is DOUBLE PRECISION
 163 *>
 164 *>         These are passed as arguments in order to save their values
 165 *>         between calls to DLASQ3.
 166 *> \endverbatim
 167 *
 168 *  Authors:
 169 *  ========
 170 *
 171 *> \author Univ. of Tennessee
 172 *> \author Univ. of California Berkeley
 173 *> \author Univ. of Colorado Denver
 174 *> \author NAG Ltd.
 175 *
 176 *> \date September 2012
 177 *
 178 *> \ingroup auxOTHERcomputational
 179 *
 180 *  =====================================================================
 181       SUBROUTINE DLASQ3( I0, N0, Z, PP, DMIN, SIGMA, DESIG, QMAX, NFAIL,
 182      $                   ITER, NDIV, IEEE, TTYPE, DMIN1, DMIN2, DN, DN1,
 183      $                   DN2, G, TAU )
 184 *
 185 *  -- LAPACK computational routine (version 3.4.2) --
 186 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
 187 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
 188 *     September 2012
 189 *
 190 *     .. Scalar Arguments ..
 191       LOGICAL            IEEE
 192       INTEGER            I0, ITER, N0, NDIV, NFAIL, PP
 193       DOUBLE PRECISION   DESIG, DMIN, DMIN1, DMIN2, DN, DN1, DN2, G,
 194      $                   QMAX, SIGMA, TAU
 195 *     ..
 196 *     .. Array Arguments ..
 197       DOUBLE PRECISION   Z( * )
 198 *     ..
 199 *
 200 *  =====================================================================
 201 *
 202 *     .. Parameters ..
 203       DOUBLE PRECISION   CBIAS
 204       PARAMETER          ( CBIAS = 1.50D0 )
 205       DOUBLE PRECISION   ZERO, QURTR, HALF, ONE, TWO, HUNDRD
 206       PARAMETER          ( ZERO = 0.0D0, QURTR = 0.250D0, HALF = 0.5D0,
 207      $                     ONE = 1.0D0, TWO = 2.0D0, HUNDRD = 100.0D0 )
 208 *     ..
 209 *     .. Local Scalars ..
 210       INTEGER            IPN4, J4, N0IN, NN, TTYPE
 211       DOUBLE PRECISION   EPS, S, T, TEMP, TOL, TOL2
 212 *     ..
 213 *     .. External Subroutines ..
 214       EXTERNAL           DLASQ4, DLASQ5, DLASQ6
 215 *     ..
 216 *     .. External Function ..
 217       DOUBLE PRECISION   DLAMCH
 218       LOGICAL            DISNAN
 219       EXTERNAL           DISNAN, DLAMCH
 220 *     ..
 221 *     .. Intrinsic Functions ..
 222       INTRINSIC          ABS, MAX, MIN, SQRT
 223 *     ..
 224 *     .. Executable Statements ..
 225 *
 226
 227       N0IN = N0
 228       EPS = DLAMCH( 'Precision' )
 229       TOL = EPS*HUNDRD
 230       TOL2 = TOL**2
 231 *
 232 *     Check for deflation.
 233 *
 234    10 CONTINUE
 235 *
 236       IF( N0.LT.I0 )
 237      $   RETURN
 238       IF( N0.EQ.I0 )
 239      $   GO TO 20
 240       NN = 4*N0 + PP
 241       IF( N0.EQ.( I0+1 ) )
 242      $   GO TO 40
 243 *
 244 *     Check whether E(N0-1) is negligible, 1 eigenvalue.
 245 *
 246       IF( Z( NN-5 ).GT.TOL2*( SIGMA+Z( NN-3 ) ) .AND.
 247      $    Z( NN-2*PP-4 ).GT.TOL2*Z( NN-7 ) )
 248      $   GO TO 30
 249 *
 250    20 CONTINUE
 251 *
 252       Z( 4*N0-3 ) = Z( 4*N0+PP-3 ) + SIGMA
 253       N0 = N0 - 1
 254       GO TO 10
 255 *
 256 *     Check  whether E(N0-2) is negligible, 2 eigenvalues.
 257 *
 258    30 CONTINUE
 259 *
 260       IF( Z( NN-9 ).GT.TOL2*SIGMA .AND.
 261      $    Z( NN-2*PP-8 ).GT.TOL2*Z( NN-11 ) )
 262      $   GO TO 50
 263 *
 264    40 CONTINUE
 265 *
 266       IF( Z( NN-3 ).GT.Z( NN-7 ) ) THEN
 267          S = Z( NN-3 )
 268          Z( NN-3 ) = Z( NN-7 )
 269          Z( NN-7 ) = S
 270       END IF
 271       T = HALF*( ( Z( NN-7 )-Z( NN-3 ) )+Z( NN-5 ) )
 272       IF( Z( NN-5 ).GT.Z( NN-3 )*TOL2.AND.T.NE.ZERO ) THEN
 273          S = Z( NN-3 )*( Z( NN-5 ) / T )
 274          IF( S.LE.T ) THEN
 275             S = Z( NN-3 )*( Z( NN-5 ) /
 276      $          ( T*( ONE+SQRT( ONE+S / T ) ) ) )
 277          ELSE
 278             S = Z( NN-3 )*( Z( NN-5 ) / ( T+SQRT( T )*SQRT( T+S ) ) )
 279          END IF
 280          T = Z( NN-7 ) + ( S+Z( NN-5 ) )
 281          Z( NN-3 ) = Z( NN-3 )*( Z( NN-7 ) / T )
 282          Z( NN-7 ) = T
 283       END IF
 284       Z( 4*N0-7 ) = Z( NN-7 ) + SIGMA
 285       Z( 4*N0-3 ) = Z( NN-3 ) + SIGMA
 286       N0 = N0 - 2
 287       GO TO 10
 288 *
 289    50 CONTINUE
 290       IF( PP.EQ.2 )
 291      $   PP = 0
 292 *
 293 *     Reverse the qd-array, if warranted.
 294 *
 295
 296       IF( DMIN.LE.ZERO .OR. N0.LT.N0IN ) THEN
 297          IF( CBIAS*Z( 4*I0+PP-3 ).LT.Z( 4*N0+PP-3 ) ) THEN
 298             IPN4 = 4*( I0+N0 )
 299             DO 60 J4 = 4*I0, 2*( I0+N0-1 ), 4
 300                TEMP = Z( J4-3 )
 301                Z( J4-3 ) = Z( IPN4-J4-3 )
 302                Z( IPN4-J4-3 ) = TEMP
 303                TEMP = Z( J4-2 )
 304                Z( J4-2 ) = Z( IPN4-J4-2 )
 305                Z( IPN4-J4-2 ) = TEMP
 306                TEMP = Z( J4-1 )
 307                Z( J4-1 ) = Z( IPN4-J4-5 )
 308                Z( IPN4-J4-5 ) = TEMP
 309                TEMP = Z( J4 )
 310                Z( J4 ) = Z( IPN4-J4-4 )
 311                Z( IPN4-J4-4 ) = TEMP
 312    60       CONTINUE
 313             IF( N0-I0.LE.4 ) THEN
 314                Z( 4*N0+PP-1 ) = Z( 4*I0+PP-1 )
 315                Z( 4*N0-PP ) = Z( 4*I0-PP )
 316             END IF
 317             DMIN2 = MIN( DMIN2, Z( 4*N0+PP-1 ) )
 318             Z( 4*N0+PP-1 ) = MIN( Z( 4*N0+PP-1 ), Z( 4*I0+PP-1 ),
 319      $                            Z( 4*I0+PP+3 ) )
 320             Z( 4*N0-PP ) = MIN( Z( 4*N0-PP ), Z( 4*I0-PP ),
 321      $                          Z( 4*I0-PP+4 ) )
 322             QMAX = MAX( QMAX, Z( 4*I0+PP-3 ), Z( 4*I0+PP+1 ) )
 323             DMIN = -ZERO
 324          END IF
 325       END IF
 326 *
 327 *     Choose a shift.
 328 *
 329       ! Print out DLASQ4 test cases
 330       write(4,*) "{"
 331       write(4,'(9999(g0))',advance="no") "z: []float64{"
 332       do i = 1, NN
 333         write (4,'(99999(e24.16,a))',advance="no") z(i), ","
 334       end do
 335       write (4,*) "},"
 336       write (4,*) "i0: ", I0, ","
 337       write (4,*) "n0: ", N0, ","
 338       write (4,*) "pp:   ", PP, ","
 339       write (4,*) "n0in: ", N0IN, ","
 340       write (4,*) "dmin: ", DMIN, ","
 341       write (4,*) "dmin1:", DMIN1, ","
 342       write (4,*) "dmin2:", DMIN2, ","
 343       write (4,*) "dn:   ", DN, ","
 344       write (4,*) "dn1:  ", DN1, ","
 345       write (4,*) "dn2:  ", DN2, ","
 346       write (4,*) "tau: ", TAU, ","
 347       write (4,*) "ttype: ", TTYPE, ","
 348       write (4,*) "g:    ", G, ","
 349       CALL DLASQ4( I0, N0, Z, PP, N0IN, DMIN, DMIN1, DMIN2, DN, DN1,
 350      $             DN2, TAU, TTYPE, G )
 351
 352       write(4,'(9999(g0))',advance="no") "zOut: []float64{"
 353       do i = 1, NN
 354         write (4,'(99999(e24.16,a))',advance="no") z(i), ","
 355       end do
 356       write (4,*) "},"
 357       write (4,*) "tauOut: ", TAU, ","
 358       write (4,*) "ttypeOut: ", TTYPE, ","
 359       write (4,*) "gOut:   ", G, ","
 360       write (4,*) "},"
 361
 362 *
 363 *     Call dqds until DMIN > 0.
 364 *
 365    70 CONTINUE
 366 *
 367
 368       write(5,*) "{"
 369       write(5,'(9999(g0))',advance="no") "z: []float64{"
 370       do i = 1, NN
 371         write (5,'(99999(e24.16,a))',advance="no") z(i), ","
 372       end do
 373       write (5,*) "},"
 374       write (5,*) "i0: ", I0, ","
 375       write (5,*) "n0: ", N0, ","
 376       write (5,*) "pp:   ", PP, ","
 377       write (5,*) "tau: ", TAU, ","
 378       write (5,*) "sigma: ", SIGMA, ","
 379       write (5,*) "dmin: ", DMIN, ","
 380       write (5,*) "dmin1:", DMIN1, ","
 381       write (5,*) "dmin2:", DMIN2, ","
 382       write (5,*) "dn:   ", DN, ","
 383       write (5,*) "dnm1:  ", DN1, ","
 384       write (5,*) "dnm2:  ", DN2, ","
 385
 386
 387       CALL DLASQ5( I0, N0, Z, PP, TAU, SIGMA, DMIN, DMIN1, DMIN2, DN,
 388      $             DN1, DN2, IEEE, EPS )
 389
 390
 391
 392       write (5,*) "i0Out: ", I0, ","
 393       write (5,*) "n0Out: ", N0, ","
 394       write (5,*) "ppOut:   ", PP, ","
 395       write (5,*) "tauOut: ", TAU, ","
 396       write (5,*) "sigmaOut: ", SIGMA, ","
 397       write (5,*) "dminOut: ", DMIN, ","
 398       write (5,*) "dmin1Out:", DMIN1, ","
 399       write (5,*) "dmin2Out:", DMIN2, ","
 400       write (5,*) "dnOut:   ", DN, ","
 401       write (5,*) "dnm1Out:  ", DN1, ","
 402       write (5,*) "dnm2Out:  ", DN2, ","
 403       write (5,*) "},"
 404
 405 *
 406       NDIV = NDIV + ( N0-I0+2 )
 407
 408       ITER = ITER + 1
 409 *
 410 *     Check status.
 411 *
 412
 413       IF( DMIN.GE.ZERO .AND. DMIN1.GE.ZERO ) THEN
 414 *
 415 *        Success.
 416 *
 417          GO TO 90
 418 *
 419       ELSE IF( DMIN.LT.ZERO .AND. DMIN1.GT.ZERO .AND.
 420      $         Z( 4*( N0-1 )-PP ).LT.TOL*( SIGMA+DN1 ) .AND.
 421      $         ABS( DN ).LT.TOL*SIGMA ) THEN
 422
 423 *
 424 *        Convergence hidden by negative DN.
 425 *
 426          Z( 4*( N0-1 )-PP+2 ) = ZERO
 427          DMIN = ZERO
 428          GO TO 90
 429       ELSE IF( DMIN.LT.ZERO ) THEN
 430
 431 *
 432 *        TAU too big. Select new TAU and try again.
 433 *
 434          NFAIL = NFAIL + 1
 435          IF( TTYPE.LT.-22 ) THEN
 436 *
 437 *           Failed twice. Play it safe.
 438 *
 439             TAU = ZERO
 440          ELSE IF( DMIN1.GT.ZERO ) THEN
 441 *
 442 *           Late failure. Gives excellent shift.
 443 *
 444             TAU = ( TAU+DMIN )*( ONE-TWO*EPS )
 445             TTYPE = TTYPE - 11
 446          ELSE
 447 *
 448 *           Early failure. Divide by 4.
 449 *
 450             TAU = QURTR*TAU
 451             TTYPE = TTYPE - 12
 452          END IF
 453          GO TO 70
 454       ELSE IF( DISNAN( DMIN ) ) THEN
 455 *
 456 *        NaN.
 457 *
 458          IF( TAU.EQ.ZERO ) THEN
 459             GO TO 80
 460          ELSE
 461             TAU = ZERO
 462             GO TO 70
 463          END IF
 464       ELSE
 465 *
 466 *        Possible underflow. Play it safe.
 467 *
 468          GO TO 80
 469       END IF
 470 *
 471 *     Risk of underflow.
 472 *
 473    80 CONTINUE
 474
 475       CALL DLASQ6( I0, N0, Z, PP, DMIN, DMIN1, DMIN2, DN, DN1, DN2 )
 476
 477
 478       NDIV = NDIV + ( N0-I0+2 )
 479       ITER = ITER + 1
 480       TAU = ZERO
 481 *
 482    90 CONTINUE
 483
 484       IF( TAU.LT.SIGMA ) THEN
 485          DESIG = DESIG + TAU
 486          T = SIGMA + DESIG
 487          DESIG = DESIG - ( T-SIGMA )
 488       ELSE
 489          T = SIGMA + TAU
 490          DESIG = SIGMA - ( T-TAU ) + DESIG
 491       END IF
 492       SIGMA = T
 493 *
 494       RETURN
 495 *
 496 *     End of DLASQ3
 497 *
 498       END