SRC/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       SAFMIN = DLAMCH( 'Safe minimum' )
 158       J4 = 4*I0 + PP - 3
 159       EMIN = Z( J4+4 )
 160       D = Z( J4 )
 161       DMIN = D
 162 *
 163       IF( PP.EQ.0 ) THEN
 164          DO 10 J4 = 4*I0, 4*( N0-3 ), 4
 165             Z( J4-2 ) = D + Z( J4-1 )
 166             IF( Z( J4-2 ).EQ.ZERO ) THEN
 167                Z( J4 ) = ZERO
 168                D = Z( J4+1 )
 169                DMIN = D
 170                EMIN = ZERO
 171             ELSE IF( SAFMIN*Z( J4+1 ).LT.Z( J4-2 ) .AND.
 172      $               SAFMIN*Z( J4-2 ).LT.Z( J4+1 ) ) THEN
 173                TEMP = Z( J4+1 ) / Z( J4-2 )
 174                Z( J4 ) = Z( J4-1 )*TEMP
 175                D = D*TEMP
 176             ELSE
 177                Z( J4 ) = Z( J4+1 )*( Z( J4-1 ) / Z( J4-2 ) )
 178                D = Z( J4+1 )*( D / Z( J4-2 ) )
 179             END IF
 180             DMIN = MIN( DMIN, D )
 181             EMIN = MIN( EMIN, Z( J4 ) )
 182    10    CONTINUE
 183       ELSE
 184          DO 20 J4 = 4*I0, 4*( N0-3 ), 4
 185             Z( J4-3 ) = D + Z( J4 )
 186             IF( Z( J4-3 ).EQ.ZERO ) THEN
 187                Z( J4-1 ) = ZERO
 188                D = Z( J4+2 )
 189                DMIN = D
 190                EMIN = ZERO
 191             ELSE IF( SAFMIN*Z( J4+2 ).LT.Z( J4-3 ) .AND.
 192      $               SAFMIN*Z( J4-3 ).LT.Z( J4+2 ) ) THEN
 193                TEMP = Z( J4+2 ) / Z( J4-3 )
 194                Z( J4-1 ) = Z( J4 )*TEMP
 195                D = D*TEMP
 196             ELSE
 197                Z( J4-1 ) = Z( J4+2 )*( Z( J4 ) / Z( J4-3 ) )
 198                D = Z( J4+2 )*( D / Z( J4-3 ) )
 199             END IF
 200             DMIN = MIN( DMIN, D )
 201             EMIN = MIN( EMIN, Z( J4-1 ) )
 202    20    CONTINUE
 203       END IF
 204 *
 205 *     Unroll last two steps.
 206 *
 207       DNM2 = D
 208       DMIN2 = DMIN
 209       J4 = 4*( N0-2 ) - PP
 210       J4P2 = J4 + 2*PP - 1
 211       Z( J4-2 ) = DNM2 + Z( J4P2 )
 212       IF( Z( J4-2 ).EQ.ZERO ) THEN
 213          Z( J4 ) = ZERO
 214          DNM1 = Z( J4P2+2 )
 215          DMIN = DNM1
 216          EMIN = ZERO
 217       ELSE IF( SAFMIN*Z( J4P2+2 ).LT.Z( J4-2 ) .AND.
 218      $         SAFMIN*Z( J4-2 ).LT.Z( J4P2+2 ) ) THEN
 219          TEMP = Z( J4P2+2 ) / Z( J4-2 )
 220          Z( J4 ) = Z( J4P2 )*TEMP
 221          DNM1 = DNM2*TEMP
 222       ELSE
 223          Z( J4 ) = Z( J4P2+2 )*( Z( J4P2 ) / Z( J4-2 ) )
 224          DNM1 = Z( J4P2+2 )*( DNM2 / Z( J4-2 ) )
 225       END IF
 226       DMIN = MIN( DMIN, DNM1 )
 227 *
 228       DMIN1 = DMIN
 229       J4 = J4 + 4
 230       J4P2 = J4 + 2*PP - 1
 231       Z( J4-2 ) = DNM1 + Z( J4P2 )
 232       IF( Z( J4-2 ).EQ.ZERO ) THEN
 233          Z( J4 ) = ZERO
 234          DN = Z( J4P2+2 )
 235          DMIN = DN
 236          EMIN = ZERO
 237       ELSE IF( SAFMIN*Z( J4P2+2 ).LT.Z( J4-2 ) .AND.
 238      $         SAFMIN*Z( J4-2 ).LT.Z( J4P2+2 ) ) THEN
 239          TEMP = Z( J4P2+2 ) / Z( J4-2 )
 240          Z( J4 ) = Z( J4P2 )*TEMP
 241          DN = DNM1*TEMP
 242       ELSE
 243          Z( J4 ) = Z( J4P2+2 )*( Z( J4P2 ) / Z( J4-2 ) )
 244          DN = Z( J4P2+2 )*( DNM1 / Z( J4-2 ) )
 245       END IF
 246       DMIN = MIN( DMIN, DN )
 247 *
 248       Z( J4+2 ) = DN
 249       Z( 4*N0-PP ) = EMIN
 250       RETURN
 251 *
 252 *     End of DLASQ6
 253 *
 254       END