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

   1 *> \brief \b DLASQ6 computes one dqd 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 DLASQ6 + dependencies
  10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dlasq6.f">
  11 *> [TGZ]</a>
  12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dlasq6.f">
  13 *> [ZIP]</a>
  14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dlasq6.f">
  15 *> [TXT]</a>
  16 *> \endhtmlonly
  17 *
  18 *  Definition:
  19 *  ===========
  20 *
  21 *       SUBROUTINE DLASQ6( I0, N0, Z, PP, DMIN, DMIN1, DMIN2, DN,
  22 *                          DNM1, DNM2 )
  23 *
  24 *       .. Scalar Arguments ..
  25 *       INTEGER            I0, N0, PP
  26 *       DOUBLE PRECISION   DMIN, DMIN1, DMIN2, DN, DNM1, DNM2
  27 *       ..
  28 *       .. Array Arguments ..
  29 *       DOUBLE PRECISION   Z( * )
  30 *       ..
  31 *
  32 *
  33 *> \par Purpose:
  34 *  =============
  35 *>
  36 *> \verbatim
  37 *>
  38 *> DLASQ6 computes one dqd (shift equal to zero) transform in
  39 *> ping-pong form, with protection against underflow and overflow.
  40 *> \endverbatim
  41 *
  42 *  Arguments:
  43 *  ==========
  44 *
  45 *> \param[in] I0
  46 *> \verbatim
  47 *>          I0 is INTEGER
  48 *>        First index.
  49 *> \endverbatim
  50 *>
  51 *> \param[in] N0
  52 *> \verbatim
  53 *>          N0 is INTEGER
  54 *>        Last index.
  55 *> \endverbatim
  56 *>
  57 *> \param[in] Z
  58 *> \verbatim
  59 *>          Z is DOUBLE PRECISION array, dimension ( 4*N )
  60 *>        Z holds the qd array. EMIN is stored in Z(4*N0) to avoid
  61 *>        an extra argument.
  62 *> \endverbatim
  63 *>
  64 *> \param[in] PP
  65 *> \verbatim
  66 *>          PP is INTEGER
  67 *>        PP=0 for ping, PP=1 for pong.
  68 *> \endverbatim
  69 *>
  70 *> \param[out] DMIN
  71 *> \verbatim
  72 *>          DMIN is DOUBLE PRECISION
  73 *>        Minimum value of d.
  74 *> \endverbatim
  75 *>
  76 *> \param[out] DMIN1
  77 *> \verbatim
  78 *>          DMIN1 is DOUBLE PRECISION
  79 *>        Minimum value of d, excluding D( N0 ).
  80 *> \endverbatim
  81 *>
  82 *> \param[out] DMIN2
  83 *> \verbatim
  84 *>          DMIN2 is DOUBLE PRECISION
  85 *>        Minimum value of d, excluding D( N0 ) and D( N0-1 ).
  86 *> \endverbatim
  87 *>
  88 *> \param[out] DN
  89 *> \verbatim
  90 *>          DN is DOUBLE PRECISION
  91 *>        d(N0), the last value of d.
  92 *> \endverbatim
  93 *>
  94 *> \param[out] DNM1
  95 *> \verbatim
  96 *>          DNM1 is DOUBLE PRECISION
  97 *>        d(N0-1).
  98 *> \endverbatim
  99 *>
 100 *> \param[out] DNM2
 101 *> \verbatim
 102 *>          DNM2 is DOUBLE PRECISION
 103 *>        d(N0-2).
 104 *> \endverbatim
 105 *
 106 *  Authors:
 107 *  ========
 108 *
 109 *> \author Univ. of Tennessee
 110 *> \author Univ. of California Berkeley
 111 *> \author Univ. of Colorado Denver
 112 *> \author NAG Ltd.
 113 *
 114 *> \date September 2012
 115 *
 116 *> \ingroup auxOTHERcomputational
 117 *
 118 *  =====================================================================
 119       SUBROUTINE DLASQ6( I0, N0, Z, PP, DMIN, DMIN1, DMIN2, DN,
 120      $                   DNM1, DNM2 )
 121 *
 122 *  -- LAPACK computational routine (version 3.4.2) --
 123 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
 124 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
 125 *     September 2012
 126 *
 127 *     .. Scalar Arguments ..
 128       INTEGER            I0, N0, PP
 129       DOUBLE PRECISION   DMIN, DMIN1, DMIN2, DN, DNM1, DNM2
 130 *     ..
 131 *     .. Array Arguments ..
 132       DOUBLE PRECISION   Z( * )
 133 *     ..
 134 *
 135 *  =====================================================================
 136 *
 137 *     .. Parameter ..
 138       DOUBLE PRECISION   ZERO
 139       PARAMETER          ( ZERO = 0.0D0 )
 140 *     ..
 141 *     .. Local Scalars ..
 142       INTEGER            J4, J4P2
 143       DOUBLE PRECISION   D, EMIN, SAFMIN, TEMP
 144 *     ..
 145 *     .. External Function ..
 146       DOUBLE PRECISION   DLAMCH
 147       EXTERNAL           DLAMCH
 148 *     ..
 149 *     .. Intrinsic Functions ..
 150       INTRINSIC          MIN
 151 *     ..
 152 *     .. Executable Statements ..
 153 *
 154       IF( ( N0-I0-1 ).LE.0 )
 155      $   RETURN
 156 *
 157
 158       print *, "In dlasq6"
 159       STOP
 160
 161       SAFMIN = DLAMCH( 'Safe minimum' )
 162       J4 = 4*I0 + PP - 3
 163       EMIN = Z( J4+4 )
 164       D = Z( J4 )
 165       DMIN = D
 166 *
 167       IF( PP.EQ.0 ) THEN
 168          DO 10 J4 = 4*I0, 4*( N0-3 ), 4
 169             Z( J4-2 ) = D + Z( J4-1 )
 170             IF( Z( J4-2 ).EQ.ZERO ) THEN
 171                Z( J4 ) = ZERO
 172                D = Z( J4+1 )
 173                DMIN = D
 174                EMIN = ZERO
 175             ELSE IF( SAFMIN*Z( J4+1 ).LT.Z( J4-2 ) .AND.
 176      $               SAFMIN*Z( J4-2 ).LT.Z( J4+1 ) ) THEN
 177                TEMP = Z( J4+1 ) / Z( J4-2 )
 178                Z( J4 ) = Z( J4-1 )*TEMP
 179                D = D*TEMP
 180             ELSE
 181                Z( J4 ) = Z( J4+1 )*( Z( J4-1 ) / Z( J4-2 ) )
 182                D = Z( J4+1 )*( D / Z( J4-2 ) )
 183             END IF
 184             DMIN = MIN( DMIN, D )
 185             EMIN = MIN( EMIN, Z( J4 ) )
 186    10    CONTINUE
 187       ELSE
 188          DO 20 J4 = 4*I0, 4*( N0-3 ), 4
 189             Z( J4-3 ) = D + Z( J4 )
 190             IF( Z( J4-3 ).EQ.ZERO ) THEN
 191                Z( J4-1 ) = ZERO
 192                D = Z( J4+2 )
 193                DMIN = D
 194                EMIN = ZERO
 195             ELSE IF( SAFMIN*Z( J4+2 ).LT.Z( J4-3 ) .AND.
 196      $               SAFMIN*Z( J4-3 ).LT.Z( J4+2 ) ) THEN
 197                TEMP = Z( J4+2 ) / Z( J4-3 )
 198                Z( J4-1 ) = Z( J4 )*TEMP
 199                D = D*TEMP
 200             ELSE
 201                Z( J4-1 ) = Z( J4+2 )*( Z( J4 ) / Z( J4-3 ) )
 202                D = Z( J4+2 )*( D / Z( J4-3 ) )
 203             END IF
 204             DMIN = MIN( DMIN, D )
 205             EMIN = MIN( EMIN, Z( J4-1 ) )
 206    20    CONTINUE
 207       END IF
 208 *
 209 *     Unroll last two steps.
 210 *
 211       DNM2 = D
 212       DMIN2 = DMIN
 213       J4 = 4*( N0-2 ) - PP
 214       J4P2 = J4 + 2*PP - 1
 215       Z( J4-2 ) = DNM2 + Z( J4P2 )
 216       IF( Z( J4-2 ).EQ.ZERO ) THEN
 217          Z( J4 ) = ZERO
 218          DNM1 = Z( J4P2+2 )
 219          DMIN = DNM1
 220          EMIN = ZERO
 221       ELSE IF( SAFMIN*Z( J4P2+2 ).LT.Z( J4-2 ) .AND.
 222      $         SAFMIN*Z( J4-2 ).LT.Z( J4P2+2 ) ) THEN
 223          TEMP = Z( J4P2+2 ) / Z( J4-2 )
 224          Z( J4 ) = Z( J4P2 )*TEMP
 225          DNM1 = DNM2*TEMP
 226       ELSE
 227          Z( J4 ) = Z( J4P2+2 )*( Z( J4P2 ) / Z( J4-2 ) )
 228          DNM1 = Z( J4P2+2 )*( DNM2 / Z( J4-2 ) )
 229       END IF
 230       DMIN = MIN( DMIN, DNM1 )
 231 *
 232       DMIN1 = DMIN
 233       J4 = J4 + 4
 234       J4P2 = J4 + 2*PP - 1
 235       Z( J4-2 ) = DNM1 + Z( J4P2 )
 236       IF( Z( J4-2 ).EQ.ZERO ) THEN
 237          Z( J4 ) = ZERO
 238          DN = Z( J4P2+2 )
 239          DMIN = DN
 240          EMIN = ZERO
 241       ELSE IF( SAFMIN*Z( J4P2+2 ).LT.Z( J4-2 ) .AND.
 242      $         SAFMIN*Z( J4-2 ).LT.Z( J4P2+2 ) ) THEN
 243          TEMP = Z( J4P2+2 ) / Z( J4-2 )
 244          Z( J4 ) = Z( J4P2 )*TEMP
 245          DN = DNM1*TEMP
 246       ELSE
 247          Z( J4 ) = Z( J4P2+2 )*( Z( J4P2 ) / Z( J4-2 ) )
 248          DN = Z( J4P2+2 )*( DNM1 / Z( J4-2 ) )
 249       END IF
 250       DMIN = MIN( DMIN, DN )
 251 *
 252       Z( J4+2 ) = DN
 253       Z( 4*N0-PP ) = EMIN
 254       RETURN
 255 *
 256 *     End of DLASQ6
 257 *
 258       END