001:       SUBROUTINE SLASQ3( I0, N0, Z, PP, DMIN, SIGMA, DESIG, QMAX, NFAIL,
002:      $                   ITER, NDIV, IEEE, TTYPE, DMIN1, DMIN2, DN, DN1,
003:      $                   DN2, G, TAU )
004: *
005: *  -- LAPACK routine (version 3.2)                                    --
006: *
007: *  -- Contributed by Osni Marques of the Lawrence Berkeley National   --
008: *  -- Laboratory and Beresford Parlett of the Univ. of California at  --
009: *  -- Berkeley                                                        --
010: *  -- November 2008                                                   --
011: *
012: *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
013: *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
014: *
015: *     .. Scalar Arguments ..
016:       LOGICAL            IEEE
017:       INTEGER            I0, ITER, N0, NDIV, NFAIL, PP
018:       REAL               DESIG, DMIN, DMIN1, DMIN2, DN, DN1, DN2, G,
019:      $                   QMAX, SIGMA, TAU
020: *     ..
021: *     .. Array Arguments ..
022:       REAL               Z( * )
023: *     ..
024: *
025: *  Purpose
026: *  =======
027: *
028: *  SLASQ3 checks for deflation, computes a shift (TAU) and calls dqds.
029: *  In case of failure it changes shifts, and tries again until output
030: *  is positive.
031: *
032: *  Arguments
033: *  =========
034: *
035: *  I0     (input) INTEGER
036: *         First index.
037: *
038: *  N0     (input) INTEGER
039: *         Last index.
040: *
041: *  Z      (input) REAL array, dimension ( 4*N )
042: *         Z holds the qd array.
043: *
044: *  PP     (input/output) INTEGER
045: *         PP=0 for ping, PP=1 for pong.
046: *         PP=2 indicates that flipping was applied to the Z array   
047: *         and that the initial tests for deflation should not be 
048: *         performed.
049: *
050: *  DMIN   (output) REAL
051: *         Minimum value of d.
052: *
053: *  SIGMA  (output) REAL
054: *         Sum of shifts used in current segment.
055: *
056: *  DESIG  (input/output) REAL
057: *         Lower order part of SIGMA
058: *
059: *  QMAX   (input) REAL
060: *         Maximum value of q.
061: *
062: *  NFAIL  (output) INTEGER
063: *         Number of times shift was too big.
064: *
065: *  ITER   (output) INTEGER
066: *         Number of iterations.
067: *
068: *  NDIV   (output) INTEGER
069: *         Number of divisions.
070: *
071: *  IEEE   (input) LOGICAL
072: *         Flag for IEEE or non IEEE arithmetic (passed to SLASQ5).
073: *
074: *  TTYPE  (input/output) INTEGER
075: *         Shift type.
076: *
077: *  DMIN1, DMIN2, DN, DN1, DN2, G, TAU (input/output) REAL
078: *         These are passed as arguments in order to save their values
079: *         between calls to SLASQ3.
080: *
081: *  =====================================================================
082: *
083: *     .. Parameters ..
084:       REAL               CBIAS
085:       PARAMETER          ( CBIAS = 1.50E0 )
086:       REAL               ZERO, QURTR, HALF, ONE, TWO, HUNDRD
087:       PARAMETER          ( ZERO = 0.0E0, QURTR = 0.250E0, HALF = 0.5E0,
088:      $                     ONE = 1.0E0, TWO = 2.0E0, HUNDRD = 100.0E0 )
089: *     ..
090: *     .. Local Scalars ..
091:       INTEGER            IPN4, J4, N0IN, NN, TTYPE
092:       REAL               EPS, S, T, TEMP, TOL, TOL2
093: *     ..
094: *     .. External Subroutines ..
095:       EXTERNAL           SLASQ4, SLASQ5, SLASQ6
096: *     ..
097: *     .. External Function ..
098:       REAL               SLAMCH
099:       LOGICAL            SISNAN
100:       EXTERNAL           SISNAN, SLAMCH
101: *     ..
102: *     .. Intrinsic Functions ..
103:       INTRINSIC          ABS, MAX, MIN, SQRT
104: *     ..
105: *     .. Executable Statements ..
106: *
107:       N0IN = N0
108:       EPS = SLAMCH( 'Precision' )
109:       TOL = EPS*HUNDRD
110:       TOL2 = TOL**2
111: *
112: *     Check for deflation.
113: *
114:    10 CONTINUE
115: *
116:       IF( N0.LT.I0 )
117:      $   RETURN
118:       IF( N0.EQ.I0 )
119:      $   GO TO 20
120:       NN = 4*N0 + PP
121:       IF( N0.EQ.( I0+1 ) )
122:      $   GO TO 40
123: *
124: *     Check whether E(N0-1) is negligible, 1 eigenvalue.
125: *
126:       IF( Z( NN-5 ).GT.TOL2*( SIGMA+Z( NN-3 ) ) .AND.
127:      $    Z( NN-2*PP-4 ).GT.TOL2*Z( NN-7 ) )
128:      $   GO TO 30
129: *
130:    20 CONTINUE
131: *
132:       Z( 4*N0-3 ) = Z( 4*N0+PP-3 ) + SIGMA
133:       N0 = N0 - 1
134:       GO TO 10
135: *
136: *     Check  whether E(N0-2) is negligible, 2 eigenvalues.
137: *
138:    30 CONTINUE
139: *
140:       IF( Z( NN-9 ).GT.TOL2*SIGMA .AND.
141:      $    Z( NN-2*PP-8 ).GT.TOL2*Z( NN-11 ) )
142:      $   GO TO 50
143: *
144:    40 CONTINUE
145: *
146:       IF( Z( NN-3 ).GT.Z( NN-7 ) ) THEN
147:          S = Z( NN-3 )
148:          Z( NN-3 ) = Z( NN-7 )
149:          Z( NN-7 ) = S
150:       END IF
151:       IF( Z( NN-5 ).GT.Z( NN-3 )*TOL2 ) THEN
152:          T = HALF*( ( Z( NN-7 )-Z( NN-3 ) )+Z( NN-5 ) )
153:          S = Z( NN-3 )*( Z( NN-5 ) / T )
154:          IF( S.LE.T ) THEN
155:             S = Z( NN-3 )*( Z( NN-5 ) /
156:      $          ( T*( ONE+SQRT( ONE+S / T ) ) ) )
157:          ELSE
158:             S = Z( NN-3 )*( Z( NN-5 ) / ( T+SQRT( T )*SQRT( T+S ) ) )
159:          END IF
160:          T = Z( NN-7 ) + ( S+Z( NN-5 ) )
161:          Z( NN-3 ) = Z( NN-3 )*( Z( NN-7 ) / T )
162:          Z( NN-7 ) = T
163:       END IF
164:       Z( 4*N0-7 ) = Z( NN-7 ) + SIGMA
165:       Z( 4*N0-3 ) = Z( NN-3 ) + SIGMA
166:       N0 = N0 - 2
167:       GO TO 10
168: *
169:    50 CONTINUE
170:       IF( PP.EQ.2 ) 
171:      $   PP = 0
172: *
173: *     Reverse the qd-array, if warranted.
174: *
175:       IF( DMIN.LE.ZERO .OR. N0.LT.N0IN ) THEN
176:          IF( CBIAS*Z( 4*I0+PP-3 ).LT.Z( 4*N0+PP-3 ) ) THEN
177:             IPN4 = 4*( I0+N0 )
178:             DO 60 J4 = 4*I0, 2*( I0+N0-1 ), 4
179:                TEMP = Z( J4-3 )
180:                Z( J4-3 ) = Z( IPN4-J4-3 )
181:                Z( IPN4-J4-3 ) = TEMP
182:                TEMP = Z( J4-2 )
183:                Z( J4-2 ) = Z( IPN4-J4-2 )
184:                Z( IPN4-J4-2 ) = TEMP
185:                TEMP = Z( J4-1 )
186:                Z( J4-1 ) = Z( IPN4-J4-5 )
187:                Z( IPN4-J4-5 ) = TEMP
188:                TEMP = Z( J4 )
189:                Z( J4 ) = Z( IPN4-J4-4 )
190:                Z( IPN4-J4-4 ) = TEMP
191:    60       CONTINUE
192:             IF( N0-I0.LE.4 ) THEN
193:                Z( 4*N0+PP-1 ) = Z( 4*I0+PP-1 )
194:                Z( 4*N0-PP ) = Z( 4*I0-PP )
195:             END IF
196:             DMIN2 = MIN( DMIN2, Z( 4*N0+PP-1 ) )
197:             Z( 4*N0+PP-1 ) = MIN( Z( 4*N0+PP-1 ), Z( 4*I0+PP-1 ),
198:      $                            Z( 4*I0+PP+3 ) )
199:             Z( 4*N0-PP ) = MIN( Z( 4*N0-PP ), Z( 4*I0-PP ),
200:      $                          Z( 4*I0-PP+4 ) )
201:             QMAX = MAX( QMAX, Z( 4*I0+PP-3 ), Z( 4*I0+PP+1 ) )
202:             DMIN = -ZERO
203:          END IF
204:       END IF
205: *
206: *     Choose a shift.
207: *
208:       CALL SLASQ4( I0, N0, Z, PP, N0IN, DMIN, DMIN1, DMIN2, DN, DN1,
209:      $             DN2, TAU, TTYPE, G )
210: *
211: *     Call dqds until DMIN > 0.
212: *
213:    70 CONTINUE
214: *
215:       CALL SLASQ5( I0, N0, Z, PP, TAU, DMIN, DMIN1, DMIN2, DN,
216:      $             DN1, DN2, IEEE )
217: *
218:       NDIV = NDIV + ( N0-I0+2 )
219:       ITER = ITER + 1
220: *
221: *     Check status.
222: *
223:       IF( DMIN.GE.ZERO .AND. DMIN1.GT.ZERO ) THEN
224: *
225: *        Success.
226: *
227:          GO TO 90
228: *
229:       ELSE IF( DMIN.LT.ZERO .AND. DMIN1.GT.ZERO .AND. 
230:      $         Z( 4*( N0-1 )-PP ).LT.TOL*( SIGMA+DN1 ) .AND.
231:      $         ABS( DN ).LT.TOL*SIGMA ) THEN
232: *
233: *        Convergence hidden by negative DN.
234: *
235:          Z( 4*( N0-1 )-PP+2 ) = ZERO
236:          DMIN = ZERO
237:          GO TO 90
238:       ELSE IF( DMIN.LT.ZERO ) THEN
239: *
240: *        TAU too big. Select new TAU and try again.
241: *
242:          NFAIL = NFAIL + 1
243:          IF( TTYPE.LT.-22 ) THEN
244: *
245: *           Failed twice. Play it safe.
246: *
247:             TAU = ZERO
248:          ELSE IF( DMIN1.GT.ZERO ) THEN
249: *
250: *           Late failure. Gives excellent shift.
251: *
252:             TAU = ( TAU+DMIN )*( ONE-TWO*EPS )
253:             TTYPE = TTYPE - 11
254:          ELSE
255: *
256: *           Early failure. Divide by 4.
257: *
258:             TAU = QURTR*TAU
259:             TTYPE = TTYPE - 12
260:          END IF
261:          GO TO 70
262:       ELSE IF( SISNAN( DMIN ) ) THEN
263: *
264: *        NaN.
265: *
266:          IF( TAU.EQ.ZERO ) THEN
267:             GO TO 80
268:          ELSE
269:             TAU = ZERO
270:             GO TO 70
271:          END IF
272:       ELSE
273: *            
274: *        Possible underflow. Play it safe.
275: *
276:          GO TO 80
277:       END IF
278: *
279: *     Risk of underflow.
280: *
281:    80 CONTINUE
282:       CALL SLASQ6( I0, N0, Z, PP, DMIN, DMIN1, DMIN2, DN, DN1, DN2 )
283:       NDIV = NDIV + ( N0-I0+2 )
284:       ITER = ITER + 1
285:       TAU = ZERO
286: *
287:    90 CONTINUE
288:       IF( TAU.LT.SIGMA ) THEN
289:          DESIG = DESIG + TAU
290:          T = SIGMA + DESIG
291:          DESIG = DESIG - ( T-SIGMA )
292:       ELSE
293:          T = SIGMA + TAU
294:          DESIG = SIGMA - ( T-TAU ) + DESIG
295:       END IF
296:       SIGMA = T
297: *
298:       RETURN
299: *
300: *     End of SLASQ3
301: *
302:       END
303: