001:       SUBROUTINE DLASQ6( I0, N0, Z, PP, DMIN, DMIN1, DMIN2, DN,
002:      $                   DNM1, DNM2 )
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:       INTEGER            I0, N0, PP
016:       DOUBLE PRECISION   DMIN, DMIN1, DMIN2, DN, DNM1, DNM2
017: *     ..
018: *     .. Array Arguments ..
019:       DOUBLE PRECISION   Z( * )
020: *     ..
021: *
022: *  Purpose
023: *  =======
024: *
025: *  DLASQ6 computes one dqd (shift equal to zero) transform in
026: *  ping-pong form, with protection against underflow and overflow.
027: *
028: *  Arguments
029: *  =========
030: *
031: *  I0    (input) INTEGER
032: *        First index.
033: *
034: *  N0    (input) INTEGER
035: *        Last index.
036: *
037: *  Z     (input) DOUBLE PRECISION array, dimension ( 4*N )
038: *        Z holds the qd array. EMIN is stored in Z(4*N0) to avoid
039: *        an extra argument.
040: *
041: *  PP    (input) INTEGER
042: *        PP=0 for ping, PP=1 for pong.
043: *
044: *  DMIN  (output) DOUBLE PRECISION
045: *        Minimum value of d.
046: *
047: *  DMIN1 (output) DOUBLE PRECISION
048: *        Minimum value of d, excluding D( N0 ).
049: *
050: *  DMIN2 (output) DOUBLE PRECISION
051: *        Minimum value of d, excluding D( N0 ) and D( N0-1 ).
052: *
053: *  DN    (output) DOUBLE PRECISION
054: *        d(N0), the last value of d.
055: *
056: *  DNM1  (output) DOUBLE PRECISION
057: *        d(N0-1).
058: *
059: *  DNM2  (output) DOUBLE PRECISION
060: *        d(N0-2).
061: *
062: *  =====================================================================
063: *
064: *     .. Parameter ..
065:       DOUBLE PRECISION   ZERO
066:       PARAMETER          ( ZERO = 0.0D0 )
067: *     ..
068: *     .. Local Scalars ..
069:       INTEGER            J4, J4P2
070:       DOUBLE PRECISION   D, EMIN, SAFMIN, TEMP
071: *     ..
072: *     .. External Function ..
073:       DOUBLE PRECISION   DLAMCH
074:       EXTERNAL           DLAMCH
075: *     ..
076: *     .. Intrinsic Functions ..
077:       INTRINSIC          MIN
078: *     ..
079: *     .. Executable Statements ..
080: *
081:       IF( ( N0-I0-1 ).LE.0 )
082:      $   RETURN
083: *
084:       SAFMIN = DLAMCH( 'Safe minimum' )
085:       J4 = 4*I0 + PP - 3
086:       EMIN = Z( J4+4 ) 
087:       D = Z( J4 )
088:       DMIN = D
089: *
090:       IF( PP.EQ.0 ) THEN
091:          DO 10 J4 = 4*I0, 4*( N0-3 ), 4
092:             Z( J4-2 ) = D + Z( J4-1 ) 
093:             IF( Z( J4-2 ).EQ.ZERO ) THEN
094:                Z( J4 ) = ZERO
095:                D = Z( J4+1 )
096:                DMIN = D
097:                EMIN = ZERO
098:             ELSE IF( SAFMIN*Z( J4+1 ).LT.Z( J4-2 ) .AND.
099:      $               SAFMIN*Z( J4-2 ).LT.Z( J4+1 ) ) THEN
100:                TEMP = Z( J4+1 ) / Z( J4-2 )
101:                Z( J4 ) = Z( J4-1 )*TEMP
102:                D = D*TEMP
103:             ELSE 
104:                Z( J4 ) = Z( J4+1 )*( Z( J4-1 ) / Z( J4-2 ) )
105:                D = Z( J4+1 )*( D / Z( J4-2 ) )
106:             END IF
107:             DMIN = MIN( DMIN, D )
108:             EMIN = MIN( EMIN, Z( J4 ) )
109:    10    CONTINUE
110:       ELSE
111:          DO 20 J4 = 4*I0, 4*( N0-3 ), 4
112:             Z( J4-3 ) = D + Z( J4 ) 
113:             IF( Z( J4-3 ).EQ.ZERO ) THEN
114:                Z( J4-1 ) = ZERO
115:                D = Z( J4+2 )
116:                DMIN = D
117:                EMIN = ZERO
118:             ELSE IF( SAFMIN*Z( J4+2 ).LT.Z( J4-3 ) .AND.
119:      $               SAFMIN*Z( J4-3 ).LT.Z( J4+2 ) ) THEN
120:                TEMP = Z( J4+2 ) / Z( J4-3 )
121:                Z( J4-1 ) = Z( J4 )*TEMP
122:                D = D*TEMP
123:             ELSE 
124:                Z( J4-1 ) = Z( J4+2 )*( Z( J4 ) / Z( J4-3 ) )
125:                D = Z( J4+2 )*( D / Z( J4-3 ) )
126:             END IF
127:             DMIN = MIN( DMIN, D )
128:             EMIN = MIN( EMIN, Z( J4-1 ) )
129:    20    CONTINUE
130:       END IF
131: *
132: *     Unroll last two steps. 
133: *
134:       DNM2 = D
135:       DMIN2 = DMIN
136:       J4 = 4*( N0-2 ) - PP
137:       J4P2 = J4 + 2*PP - 1
138:       Z( J4-2 ) = DNM2 + Z( J4P2 )
139:       IF( Z( J4-2 ).EQ.ZERO ) THEN
140:          Z( J4 ) = ZERO
141:          DNM1 = Z( J4P2+2 )
142:          DMIN = DNM1
143:          EMIN = ZERO
144:       ELSE IF( SAFMIN*Z( J4P2+2 ).LT.Z( J4-2 ) .AND.
145:      $         SAFMIN*Z( J4-2 ).LT.Z( J4P2+2 ) ) THEN
146:          TEMP = Z( J4P2+2 ) / Z( J4-2 )
147:          Z( J4 ) = Z( J4P2 )*TEMP
148:          DNM1 = DNM2*TEMP
149:       ELSE
150:          Z( J4 ) = Z( J4P2+2 )*( Z( J4P2 ) / Z( J4-2 ) )
151:          DNM1 = Z( J4P2+2 )*( DNM2 / Z( J4-2 ) )
152:       END IF
153:       DMIN = MIN( DMIN, DNM1 )
154: *
155:       DMIN1 = DMIN
156:       J4 = J4 + 4
157:       J4P2 = J4 + 2*PP - 1
158:       Z( J4-2 ) = DNM1 + Z( J4P2 )
159:       IF( Z( J4-2 ).EQ.ZERO ) THEN
160:          Z( J4 ) = ZERO
161:          DN = Z( J4P2+2 )
162:          DMIN = DN
163:          EMIN = ZERO
164:       ELSE IF( SAFMIN*Z( J4P2+2 ).LT.Z( J4-2 ) .AND.
165:      $         SAFMIN*Z( J4-2 ).LT.Z( J4P2+2 ) ) THEN
166:          TEMP = Z( J4P2+2 ) / Z( J4-2 )
167:          Z( J4 ) = Z( J4P2 )*TEMP
168:          DN = DNM1*TEMP
169:       ELSE
170:          Z( J4 ) = Z( J4P2+2 )*( Z( J4P2 ) / Z( J4-2 ) )
171:          DN = Z( J4P2+2 )*( DNM1 / Z( J4-2 ) )
172:       END IF
173:       DMIN = MIN( DMIN, DN )
174: *
175:       Z( J4+2 ) = DN
176:       Z( 4*N0-PP ) = EMIN
177:       RETURN
178: *
179: *     End of DLASQ6
180: *
181:       END
182: