Return to dlasq5.f CVS log

Up to [local] / rpl / lapack / lapack

Cohérence

    1:       SUBROUTINE DLASQ5( I0, N0, Z, PP, TAU, DMIN, DMIN1, DMIN2, DN,
    2:      $                   DNM1, DNM2, IEEE )
    3: *
    4: *  -- LAPACK routine (version 3.2)                                    --
    5: *
    6: *  -- Contributed by Osni Marques of the Lawrence Berkeley National   --
    7: *  -- Laboratory and Beresford Parlett of the Univ. of California at  --
    8: *  -- Berkeley                                                        --
    9: *  -- November 2008                                                   --
   10: *
   11: *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
   12: *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
   13: *
   14: *     .. Scalar Arguments ..
   15:       LOGICAL            IEEE
   16:       INTEGER            I0, N0, PP
   17:       DOUBLE PRECISION   DMIN, DMIN1, DMIN2, DN, DNM1, DNM2, TAU
   18: *     ..
   19: *     .. Array Arguments ..
   20:       DOUBLE PRECISION   Z( * )
   21: *     ..
   22: *
   23: *  Purpose
   24: *  =======
   25: *
   26: *  DLASQ5 computes one dqds transform in ping-pong form, one
   27: *  version for IEEE machines another for non IEEE machines.
   28: *
   29: *  Arguments
   30: *  =========
   31: *
   32: *  I0    (input) INTEGER
   33: *        First index.
   34: *
   35: *  N0    (input) INTEGER
   36: *        Last index.
   37: *
   38: *  Z     (input) DOUBLE PRECISION array, dimension ( 4*N )
   39: *        Z holds the qd array. EMIN is stored in Z(4*N0) to avoid
   40: *        an extra argument.
   41: *
   42: *  PP    (input) INTEGER
   43: *        PP=0 for ping, PP=1 for pong.
   44: *
   45: *  TAU   (input) DOUBLE PRECISION
   46: *        This is the shift.
   47: *
   48: *  DMIN  (output) DOUBLE PRECISION
   49: *        Minimum value of d.
   50: *
   51: *  DMIN1 (output) DOUBLE PRECISION
   52: *        Minimum value of d, excluding D( N0 ).
   53: *
   54: *  DMIN2 (output) DOUBLE PRECISION
   55: *        Minimum value of d, excluding D( N0 ) and D( N0-1 ).
   56: *
   57: *  DN    (output) DOUBLE PRECISION
   58: *        d(N0), the last value of d.
   59: *
   60: *  DNM1  (output) DOUBLE PRECISION
   61: *        d(N0-1).
   62: *
   63: *  DNM2  (output) DOUBLE PRECISION
   64: *        d(N0-2).
   65: *
   66: *  IEEE  (input) LOGICAL
   67: *        Flag for IEEE or non IEEE arithmetic.
   68: *
   69: *  =====================================================================
   70: *
   71: *     .. Parameter ..
   72:       DOUBLE PRECISION   ZERO
   73:       PARAMETER          ( ZERO = 0.0D0 )
   74: *     ..
   75: *     .. Local Scalars ..
   76:       INTEGER            J4, J4P2
   77:       DOUBLE PRECISION   D, EMIN, TEMP
   78: *     ..
   79: *     .. Intrinsic Functions ..
   80:       INTRINSIC          MIN
   81: *     ..
   82: *     .. Executable Statements ..
   83: *
   84:       IF( ( N0-I0-1 ).LE.0 )
   85:      $   RETURN
   86: *
   87:       J4 = 4*I0 + PP - 3
   88:       EMIN = Z( J4+4 ) 
   89:       D = Z( J4 ) - TAU
   90:       DMIN = D
   91:       DMIN1 = -Z( J4 )
   92: *
   93:       IF( IEEE ) THEN
   94: *
   95: *        Code for IEEE arithmetic.
   96: *
   97:          IF( PP.EQ.0 ) THEN
   98:             DO 10 J4 = 4*I0, 4*( N0-3 ), 4
   99:                Z( J4-2 ) = D + Z( J4-1 ) 
  100:                TEMP = Z( J4+1 ) / Z( J4-2 )
  101:                D = D*TEMP - TAU
  102:                DMIN = MIN( DMIN, D )
  103:                Z( J4 ) = Z( J4-1 )*TEMP
  104:                EMIN = MIN( Z( J4 ), EMIN )
  105:    10       CONTINUE
  106:          ELSE
  107:             DO 20 J4 = 4*I0, 4*( N0-3 ), 4
  108:                Z( J4-3 ) = D + Z( J4 ) 
  109:                TEMP = Z( J4+2 ) / Z( J4-3 )
  110:                D = D*TEMP - TAU
  111:                DMIN = MIN( DMIN, D )
  112:                Z( J4-1 ) = Z( J4 )*TEMP
  113:                EMIN = MIN( Z( J4-1 ), EMIN )
  114:    20       CONTINUE
  115:          END IF
  116: *
  117: *        Unroll last two steps. 
  118: *
  119:          DNM2 = D
  120:          DMIN2 = DMIN
  121:          J4 = 4*( N0-2 ) - PP
  122:          J4P2 = J4 + 2*PP - 1
  123:          Z( J4-2 ) = DNM2 + Z( J4P2 )
  124:          Z( J4 ) = Z( J4P2+2 )*( Z( J4P2 ) / Z( J4-2 ) )
  125:          DNM1 = Z( J4P2+2 )*( DNM2 / Z( J4-2 ) ) - TAU
  126:          DMIN = MIN( DMIN, DNM1 )
  127: *
  128:          DMIN1 = DMIN
  129:          J4 = J4 + 4
  130:          J4P2 = J4 + 2*PP - 1
  131:          Z( J4-2 ) = DNM1 + Z( J4P2 )
  132:          Z( J4 ) = Z( J4P2+2 )*( Z( J4P2 ) / Z( J4-2 ) )
  133:          DN = Z( J4P2+2 )*( DNM1 / Z( J4-2 ) ) - TAU
  134:          DMIN = MIN( DMIN, DN )
  135: *
  136:       ELSE
  137: *
  138: *        Code for non IEEE arithmetic.
  139: *
  140:          IF( PP.EQ.0 ) THEN
  141:             DO 30 J4 = 4*I0, 4*( N0-3 ), 4
  142:                Z( J4-2 ) = D + Z( J4-1 ) 
  143:                IF( D.LT.ZERO ) THEN
  144:                   RETURN
  145:                ELSE 
  146:                   Z( J4 ) = Z( J4+1 )*( Z( J4-1 ) / Z( J4-2 ) )
  147:                   D = Z( J4+1 )*( D / Z( J4-2 ) ) - TAU
  148:                END IF
  149:                DMIN = MIN( DMIN, D )
  150:                EMIN = MIN( EMIN, Z( J4 ) )
  151:    30       CONTINUE
  152:          ELSE
  153:             DO 40 J4 = 4*I0, 4*( N0-3 ), 4
  154:                Z( J4-3 ) = D + Z( J4 ) 
  155:                IF( D.LT.ZERO ) THEN
  156:                   RETURN
  157:                ELSE 
  158:                   Z( J4-1 ) = Z( J4+2 )*( Z( J4 ) / Z( J4-3 ) )
  159:                   D = Z( J4+2 )*( D / Z( J4-3 ) ) - TAU
  160:                END IF
  161:                DMIN = MIN( DMIN, D )
  162:                EMIN = MIN( EMIN, Z( J4-1 ) )
  163:    40       CONTINUE
  164:          END IF
  165: *
  166: *        Unroll last two steps. 
  167: *
  168:          DNM2 = D
  169:          DMIN2 = DMIN
  170:          J4 = 4*( N0-2 ) - PP
  171:          J4P2 = J4 + 2*PP - 1
  172:          Z( J4-2 ) = DNM2 + Z( J4P2 )
  173:          IF( DNM2.LT.ZERO ) THEN
  174:             RETURN
  175:          ELSE
  176:             Z( J4 ) = Z( J4P2+2 )*( Z( J4P2 ) / Z( J4-2 ) )
  177:             DNM1 = Z( J4P2+2 )*( DNM2 / Z( J4-2 ) ) - TAU
  178:          END IF
  179:          DMIN = MIN( DMIN, DNM1 )
  180: *
  181:          DMIN1 = DMIN
  182:          J4 = J4 + 4
  183:          J4P2 = J4 + 2*PP - 1
  184:          Z( J4-2 ) = DNM1 + Z( J4P2 )
  185:          IF( DNM1.LT.ZERO ) THEN
  186:             RETURN
  187:          ELSE
  188:             Z( J4 ) = Z( J4P2+2 )*( Z( J4P2 ) / Z( J4-2 ) )
  189:             DN = Z( J4P2+2 )*( DNM1 / Z( J4-2 ) ) - TAU
  190:          END IF
  191:          DMIN = MIN( DMIN, DN )
  192: *
  193:       END IF
  194: *
  195:       Z( J4+2 ) = DN
  196:       Z( 4*N0-PP ) = EMIN
  197:       RETURN
  198: *
  199: *     End of DLASQ5
  200: *
  201:       END