001:       SUBROUTINE DLASQ5( I0, N0, Z, PP, TAU, DMIN, DMIN1, DMIN2, DN,
002:      $                   DNM1, DNM2, IEEE )
003: *
004: *  -- LAPACK routine (version 3.2)                                    --
005: *
006: *  -- Contributed by Osni Marques of the Lawrence Berkeley National   --
007: *  -- Laboratory and Beresford Parlett of the Univ. of California at  --
008: *  -- Berkeley                                                        --
009: *  -- November 2008                                                   --
010: *
011: *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
012: *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
013: *
014: *     .. Scalar Arguments ..
015:       LOGICAL            IEEE
016:       INTEGER            I0, N0, PP
017:       DOUBLE PRECISION   DMIN, DMIN1, DMIN2, DN, DNM1, DNM2, TAU
018: *     ..
019: *     .. Array Arguments ..
020:       DOUBLE PRECISION   Z( * )
021: *     ..
022: *
023: *  Purpose
024: *  =======
025: *
026: *  DLASQ5 computes one dqds transform in ping-pong form, one
027: *  version for IEEE machines another for non IEEE machines.
028: *
029: *  Arguments
030: *  =========
031: *
032: *  I0    (input) INTEGER
033: *        First index.
034: *
035: *  N0    (input) INTEGER
036: *        Last index.
037: *
038: *  Z     (input) DOUBLE PRECISION array, dimension ( 4*N )
039: *        Z holds the qd array. EMIN is stored in Z(4*N0) to avoid
040: *        an extra argument.
041: *
042: *  PP    (input) INTEGER
043: *        PP=0 for ping, PP=1 for pong.
044: *
045: *  TAU   (input) DOUBLE PRECISION
046: *        This is the shift.
047: *
048: *  DMIN  (output) DOUBLE PRECISION
049: *        Minimum value of d.
050: *
051: *  DMIN1 (output) DOUBLE PRECISION
052: *        Minimum value of d, excluding D( N0 ).
053: *
054: *  DMIN2 (output) DOUBLE PRECISION
055: *        Minimum value of d, excluding D( N0 ) and D( N0-1 ).
056: *
057: *  DN    (output) DOUBLE PRECISION
058: *        d(N0), the last value of d.
059: *
060: *  DNM1  (output) DOUBLE PRECISION
061: *        d(N0-1).
062: *
063: *  DNM2  (output) DOUBLE PRECISION
064: *        d(N0-2).
065: *
066: *  IEEE  (input) LOGICAL
067: *        Flag for IEEE or non IEEE arithmetic.
068: *
069: *  =====================================================================
070: *
071: *     .. Parameter ..
072:       DOUBLE PRECISION   ZERO
073:       PARAMETER          ( ZERO = 0.0D0 )
074: *     ..
075: *     .. Local Scalars ..
076:       INTEGER            J4, J4P2
077:       DOUBLE PRECISION   D, EMIN, TEMP
078: *     ..
079: *     .. Intrinsic Functions ..
080:       INTRINSIC          MIN
081: *     ..
082: *     .. Executable Statements ..
083: *
084:       IF( ( N0-I0-1 ).LE.0 )
085:      $   RETURN
086: *
087:       J4 = 4*I0 + PP - 3
088:       EMIN = Z( J4+4 ) 
089:       D = Z( J4 ) - TAU
090:       DMIN = D
091:       DMIN1 = -Z( J4 )
092: *
093:       IF( IEEE ) THEN
094: *
095: *        Code for IEEE arithmetic.
096: *
097:          IF( PP.EQ.0 ) THEN
098:             DO 10 J4 = 4*I0, 4*( N0-3 ), 4
099:                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
202: