001:       SUBROUTINE SLASQ2( N, Z, INFO )
002: *
003: *  -- LAPACK routine (version 3.2)                                    --
004: *
005: *  -- Contributed by Osni Marques of the Lawrence Berkeley National   --
006: *  -- Laboratory and Beresford Parlett of the Univ. of California at  --
007: *  -- Berkeley                                                        --
008: *  -- November 2008                                                   --
009: *
010: *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
011: *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
012: *
013: *     .. Scalar Arguments ..
014:       INTEGER            INFO, N
015: *     ..
016: *     .. Array Arguments ..
017:       REAL               Z( * )
018: *     ..
019: *
020: *  Purpose
021: *  =======
022: *
023: *  SLASQ2 computes all the eigenvalues of the symmetric positive 
024: *  definite tridiagonal matrix associated with the qd array Z to high
025: *  relative accuracy are computed to high relative accuracy, in the
026: *  absence of denormalization, underflow and overflow.
027: *
028: *  To see the relation of Z to the tridiagonal matrix, let L be a
029: *  unit lower bidiagonal matrix with subdiagonals Z(2,4,6,,..) and
030: *  let U be an upper bidiagonal matrix with 1's above and diagonal
031: *  Z(1,3,5,,..). The tridiagonal is L*U or, if you prefer, the
032: *  symmetric tridiagonal to which it is similar.
033: *
034: *  Note : SLASQ2 defines a logical variable, IEEE, which is true
035: *  on machines which follow ieee-754 floating-point standard in their
036: *  handling of infinities and NaNs, and false otherwise. This variable
037: *  is passed to SLASQ3.
038: *
039: *  Arguments
040: *  =========
041: *
042: *  N     (input) INTEGER
043: *        The number of rows and columns in the matrix. N >= 0.
044: *
045: *  Z     (input/output) REAL array, dimension ( 4*N )
046: *        On entry Z holds the qd array. On exit, entries 1 to N hold
047: *        the eigenvalues in decreasing order, Z( 2*N+1 ) holds the
048: *        trace, and Z( 2*N+2 ) holds the sum of the eigenvalues. If
049: *        N > 2, then Z( 2*N+3 ) holds the iteration count, Z( 2*N+4 )
050: *        holds NDIVS/NIN^2, and Z( 2*N+5 ) holds the percentage of
051: *        shifts that failed.
052: *
053: *  INFO  (output) INTEGER
054: *        = 0: successful exit
055: *        < 0: if the i-th argument is a scalar and had an illegal
056: *             value, then INFO = -i, if the i-th argument is an
057: *             array and the j-entry had an illegal value, then
058: *             INFO = -(i*100+j)
059: *        > 0: the algorithm failed
060: *              = 1, a split was marked by a positive value in E
061: *              = 2, current block of Z not diagonalized after 30*N
062: *                   iterations (in inner while loop)
063: *              = 3, termination criterion of outer while loop not met 
064: *                   (program created more than N unreduced blocks)
065: *
066: *  Further Details
067: *  ===============
068: *  Local Variables: I0:N0 defines a current unreduced segment of Z.
069: *  The shifts are accumulated in SIGMA. Iteration count is in ITER.
070: *  Ping-pong is controlled by PP (alternates between 0 and 1).
071: *
072: *  =====================================================================
073: *
074: *     .. Parameters ..
075:       REAL               CBIAS
076:       PARAMETER          ( CBIAS = 1.50E0 )
077:       REAL               ZERO, HALF, ONE, TWO, FOUR, HUNDRD
078:       PARAMETER          ( ZERO = 0.0E0, HALF = 0.5E0, ONE = 1.0E0,
079:      $                     TWO = 2.0E0, FOUR = 4.0E0, HUNDRD = 100.0E0 )
080: *     ..
081: *     .. Local Scalars ..
082:       LOGICAL            IEEE
083:       INTEGER            I0, I4, IINFO, IPN4, ITER, IWHILA, IWHILB, K,
084:      $                   KMIN, N0, NBIG, NDIV, NFAIL, PP, SPLT, TTYPE
085:       REAL               D, DEE, DEEMIN, DESIG, DMIN, DMIN1, DMIN2, DN,
086:      $                   DN1, DN2, E, EMAX, EMIN, EPS, G, OLDEMN, QMAX,
087:      $                   QMIN, S, SAFMIN, SIGMA, T, TAU, TEMP, TOL,
088:      $                   TOL2, TRACE, ZMAX
089: *     ..
090: *     .. External Subroutines ..
091:       EXTERNAL           SLASQ3, SLASRT, XERBLA
092: *     ..
093: *     .. External Functions ..
094:       INTEGER            ILAENV
095:       REAL               SLAMCH
096:       EXTERNAL           ILAENV, SLAMCH
097: *     ..
098: *     .. Intrinsic Functions ..
099:       INTRINSIC          ABS, MAX, MIN, REAL, SQRT
100: *     ..
101: *     .. Executable Statements ..
102: *      
103: *     Test the input arguments.
104: *     (in case SLASQ2 is not called by SLASQ1)
105: *
106:       INFO = 0
107:       EPS = SLAMCH( 'Precision' )
108:       SAFMIN = SLAMCH( 'Safe minimum' )
109:       TOL = EPS*HUNDRD
110:       TOL2 = TOL**2
111: *
112:       IF( N.LT.0 ) THEN
113:          INFO = -1
114:          CALL XERBLA( 'SLASQ2', 1 )
115:          RETURN
116:       ELSE IF( N.EQ.0 ) THEN
117:          RETURN
118:       ELSE IF( N.EQ.1 ) THEN
119: *
120: *        1-by-1 case.
121: *
122:          IF( Z( 1 ).LT.ZERO ) THEN
123:             INFO = -201
124:             CALL XERBLA( 'SLASQ2', 2 )
125:          END IF
126:          RETURN
127:       ELSE IF( N.EQ.2 ) THEN
128: *
129: *        2-by-2 case.
130: *
131:          IF( Z( 2 ).LT.ZERO .OR. Z( 3 ).LT.ZERO ) THEN
132:             INFO = -2
133:             CALL XERBLA( 'SLASQ2', 2 )
134:             RETURN
135:          ELSE IF( Z( 3 ).GT.Z( 1 ) ) THEN
136:             D = Z( 3 )
137:             Z( 3 ) = Z( 1 )
138:             Z( 1 ) = D
139:          END IF
140:          Z( 5 ) = Z( 1 ) + Z( 2 ) + Z( 3 )
141:          IF( Z( 2 ).GT.Z( 3 )*TOL2 ) THEN
142:             T = HALF*( ( Z( 1 )-Z( 3 ) )+Z( 2 ) ) 
143:             S = Z( 3 )*( Z( 2 ) / T )
144:             IF( S.LE.T ) THEN
145:                S = Z( 3 )*( Z( 2 ) / ( T*( ONE+SQRT( ONE+S / T ) ) ) )
146:             ELSE
147:                S = Z( 3 )*( Z( 2 ) / ( T+SQRT( T )*SQRT( T+S ) ) )
148:             END IF
149:             T = Z( 1 ) + ( S+Z( 2 ) )
150:             Z( 3 ) = Z( 3 )*( Z( 1 ) / T )
151:             Z( 1 ) = T
152:          END IF
153:          Z( 2 ) = Z( 3 )
154:          Z( 6 ) = Z( 2 ) + Z( 1 )
155:          RETURN
156:       END IF
157: *
158: *     Check for negative data and compute sums of q's and e's.
159: *
160:       Z( 2*N ) = ZERO
161:       EMIN = Z( 2 )
162:       QMAX = ZERO
163:       ZMAX = ZERO
164:       D = ZERO
165:       E = ZERO
166: *
167:       DO 10 K = 1, 2*( N-1 ), 2
168:          IF( Z( K ).LT.ZERO ) THEN
169:             INFO = -( 200+K )
170:             CALL XERBLA( 'SLASQ2', 2 )
171:             RETURN
172:          ELSE IF( Z( K+1 ).LT.ZERO ) THEN
173:             INFO = -( 200+K+1 )
174:             CALL XERBLA( 'SLASQ2', 2 )
175:             RETURN
176:          END IF
177:          D = D + Z( K )
178:          E = E + Z( K+1 )
179:          QMAX = MAX( QMAX, Z( K ) )
180:          EMIN = MIN( EMIN, Z( K+1 ) )
181:          ZMAX = MAX( QMAX, ZMAX, Z( K+1 ) )
182:    10 CONTINUE
183:       IF( Z( 2*N-1 ).LT.ZERO ) THEN
184:          INFO = -( 200+2*N-1 )
185:          CALL XERBLA( 'SLASQ2', 2 )
186:          RETURN
187:       END IF
188:       D = D + Z( 2*N-1 )
189:       QMAX = MAX( QMAX, Z( 2*N-1 ) )
190:       ZMAX = MAX( QMAX, ZMAX )
191: *
192: *     Check for diagonality.
193: *
194:       IF( E.EQ.ZERO ) THEN
195:          DO 20 K = 2, N
196:             Z( K ) = Z( 2*K-1 )
197:    20    CONTINUE
198:          CALL SLASRT( 'D', N, Z, IINFO )
199:          Z( 2*N-1 ) = D
200:          RETURN
201:       END IF
202: *
203:       TRACE = D + E
204: *
205: *     Check for zero data.
206: *
207:       IF( TRACE.EQ.ZERO ) THEN
208:          Z( 2*N-1 ) = ZERO
209:          RETURN
210:       END IF
211: *         
212: *     Check whether the machine is IEEE conformable.
213: *         
214: *     IEEE = ILAENV( 10, 'SLASQ2', 'N', 1, 2, 3, 4 ).EQ.1 .AND.
215: *    $       ILAENV( 11, 'SLASQ2', 'N', 1, 2, 3, 4 ).EQ.1      
216: *
217: *     [11/15/2008] The case IEEE=.TRUE. has a problem in single precision with
218: *     some the test matrices of type 16. The double precision code is fine.
219: *
220:       IEEE = .FALSE.
221: *         
222: *     Rearrange data for locality: Z=(q1,qq1,e1,ee1,q2,qq2,e2,ee2,...).
223: *
224:       DO 30 K = 2*N, 2, -2
225:          Z( 2*K ) = ZERO 
226:          Z( 2*K-1 ) = Z( K ) 
227:          Z( 2*K-2 ) = ZERO 
228:          Z( 2*K-3 ) = Z( K-1 ) 
229:    30 CONTINUE
230: *
231:       I0 = 1
232:       N0 = N
233: *
234: *     Reverse the qd-array, if warranted.
235: *
236:       IF( CBIAS*Z( 4*I0-3 ).LT.Z( 4*N0-3 ) ) THEN
237:          IPN4 = 4*( I0+N0 )
238:          DO 40 I4 = 4*I0, 2*( I0+N0-1 ), 4
239:             TEMP = Z( I4-3 )
240:             Z( I4-3 ) = Z( IPN4-I4-3 )
241:             Z( IPN4-I4-3 ) = TEMP
242:             TEMP = Z( I4-1 )
243:             Z( I4-1 ) = Z( IPN4-I4-5 )
244:             Z( IPN4-I4-5 ) = TEMP
245:    40    CONTINUE
246:       END IF
247: *
248: *     Initial split checking via dqd and Li's test.
249: *
250:       PP = 0
251: *
252:       DO 80 K = 1, 2
253: *
254:          D = Z( 4*N0+PP-3 )
255:          DO 50 I4 = 4*( N0-1 ) + PP, 4*I0 + PP, -4
256:             IF( Z( I4-1 ).LE.TOL2*D ) THEN
257:                Z( I4-1 ) = -ZERO
258:                D = Z( I4-3 )
259:             ELSE
260:                D = Z( I4-3 )*( D / ( D+Z( I4-1 ) ) )
261:             END IF
262:    50    CONTINUE
263: *
264: *        dqd maps Z to ZZ plus Li's test.
265: *
266:          EMIN = Z( 4*I0+PP+1 )
267:          D = Z( 4*I0+PP-3 )
268:          DO 60 I4 = 4*I0 + PP, 4*( N0-1 ) + PP, 4
269:             Z( I4-2*PP-2 ) = D + Z( I4-1 )
270:             IF( Z( I4-1 ).LE.TOL2*D ) THEN
271:                Z( I4-1 ) = -ZERO
272:                Z( I4-2*PP-2 ) = D
273:                Z( I4-2*PP ) = ZERO
274:                D = Z( I4+1 )
275:             ELSE IF( SAFMIN*Z( I4+1 ).LT.Z( I4-2*PP-2 ) .AND.
276:      $               SAFMIN*Z( I4-2*PP-2 ).LT.Z( I4+1 ) ) THEN
277:                TEMP = Z( I4+1 ) / Z( I4-2*PP-2 )
278:                Z( I4-2*PP ) = Z( I4-1 )*TEMP
279:                D = D*TEMP
280:             ELSE
281:                Z( I4-2*PP ) = Z( I4+1 )*( Z( I4-1 ) / Z( I4-2*PP-2 ) )
282:                D = Z( I4+1 )*( D / Z( I4-2*PP-2 ) )
283:             END IF
284:             EMIN = MIN( EMIN, Z( I4-2*PP ) )
285:    60    CONTINUE 
286:          Z( 4*N0-PP-2 ) = D
287: *
288: *        Now find qmax.
289: *
290:          QMAX = Z( 4*I0-PP-2 )
291:          DO 70 I4 = 4*I0 - PP + 2, 4*N0 - PP - 2, 4
292:             QMAX = MAX( QMAX, Z( I4 ) )
293:    70    CONTINUE
294: *
295: *        Prepare for the next iteration on K.
296: *
297:          PP = 1 - PP
298:    80 CONTINUE
299: *
300: *     Initialise variables to pass to SLASQ3.
301: *
302:       TTYPE = 0
303:       DMIN1 = ZERO
304:       DMIN2 = ZERO
305:       DN    = ZERO
306:       DN1   = ZERO
307:       DN2   = ZERO
308:       G     = ZERO
309:       TAU   = ZERO
310: *
311:       ITER = 2
312:       NFAIL = 0
313:       NDIV = 2*( N0-I0 )
314: *
315:       DO 160 IWHILA = 1, N + 1
316:          IF( N0.LT.1 ) 
317:      $      GO TO 170
318: *
319: *        While array unfinished do 
320: *
321: *        E(N0) holds the value of SIGMA when submatrix in I0:N0
322: *        splits from the rest of the array, but is negated.
323: *      
324:          DESIG = ZERO
325:          IF( N0.EQ.N ) THEN
326:             SIGMA = ZERO
327:          ELSE
328:             SIGMA = -Z( 4*N0-1 )
329:          END IF
330:          IF( SIGMA.LT.ZERO ) THEN
331:             INFO = 1
332:             RETURN
333:          END IF
334: *
335: *        Find last unreduced submatrix's top index I0, find QMAX and
336: *        EMIN. Find Gershgorin-type bound if Q's much greater than E's.
337: *
338:          EMAX = ZERO 
339:          IF( N0.GT.I0 ) THEN
340:             EMIN = ABS( Z( 4*N0-5 ) )
341:          ELSE
342:             EMIN = ZERO
343:          END IF
344:          QMIN = Z( 4*N0-3 )
345:          QMAX = QMIN
346:          DO 90 I4 = 4*N0, 8, -4
347:             IF( Z( I4-5 ).LE.ZERO )
348:      $         GO TO 100
349:             IF( QMIN.GE.FOUR*EMAX ) THEN
350:                QMIN = MIN( QMIN, Z( I4-3 ) )
351:                EMAX = MAX( EMAX, Z( I4-5 ) )
352:             END IF
353:             QMAX = MAX( QMAX, Z( I4-7 )+Z( I4-5 ) )
354:             EMIN = MIN( EMIN, Z( I4-5 ) )
355:    90    CONTINUE
356:          I4 = 4 
357: *
358:   100    CONTINUE
359:          I0 = I4 / 4
360:          PP = 0
361: *
362:          IF( N0-I0.GT.1 ) THEN
363:             DEE = Z( 4*I0-3 )
364:             DEEMIN = DEE
365:             KMIN = I0
366:             DO 110 I4 = 4*I0+1, 4*N0-3, 4
367:                DEE = Z( I4 )*( DEE /( DEE+Z( I4-2 ) ) )
368:                IF( DEE.LE.DEEMIN ) THEN
369:                   DEEMIN = DEE
370:                   KMIN = ( I4+3 )/4
371:                END IF
372:   110       CONTINUE
373:             IF( (KMIN-I0)*2.LT.N0-KMIN .AND. 
374:      $         DEEMIN.LE.HALF*Z(4*N0-3) ) THEN
375:                IPN4 = 4*( I0+N0 )
376:                PP = 2
377:                DO 120 I4 = 4*I0, 2*( I0+N0-1 ), 4
378:                   TEMP = Z( I4-3 )
379:                   Z( I4-3 ) = Z( IPN4-I4-3 )
380:                   Z( IPN4-I4-3 ) = TEMP
381:                   TEMP = Z( I4-2 )
382:                   Z( I4-2 ) = Z( IPN4-I4-2 )
383:                   Z( IPN4-I4-2 ) = TEMP
384:                   TEMP = Z( I4-1 )
385:                   Z( I4-1 ) = Z( IPN4-I4-5 )
386:                   Z( IPN4-I4-5 ) = TEMP
387:                   TEMP = Z( I4 )
388:                   Z( I4 ) = Z( IPN4-I4-4 )
389:                   Z( IPN4-I4-4 ) = TEMP
390:   120          CONTINUE
391:             END IF
392:          END IF
393: *
394: *        Put -(initial shift) into DMIN.
395: *
396:          DMIN = -MAX( ZERO, QMIN-TWO*SQRT( QMIN )*SQRT( EMAX ) )
397: *
398: *        Now I0:N0 is unreduced. 
399: *        PP = 0 for ping, PP = 1 for pong.
400: *        PP = 2 indicates that flipping was applied to the Z array and
401: *               and that the tests for deflation upon entry in SLASQ3 
402: *               should not be performed.
403: *
404:          NBIG = 30*( N0-I0+1 )
405:          DO 140 IWHILB = 1, NBIG
406:             IF( I0.GT.N0 ) 
407:      $         GO TO 150
408: *
409: *           While submatrix unfinished take a good dqds step.
410: *
411:             CALL SLASQ3( I0, N0, Z, PP, DMIN, SIGMA, DESIG, QMAX, NFAIL,
412:      $                   ITER, NDIV, IEEE, TTYPE, DMIN1, DMIN2, DN, DN1,
413:      $                   DN2, G, TAU )
414: *
415:             PP = 1 - PP
416: *
417: *           When EMIN is very small check for splits.
418: *
419:             IF( PP.EQ.0 .AND. N0-I0.GE.3 ) THEN
420:                IF( Z( 4*N0 ).LE.TOL2*QMAX .OR.
421:      $             Z( 4*N0-1 ).LE.TOL2*SIGMA ) THEN
422:                   SPLT = I0 - 1
423:                   QMAX = Z( 4*I0-3 )
424:                   EMIN = Z( 4*I0-1 )
425:                   OLDEMN = Z( 4*I0 )
426:                   DO 130 I4 = 4*I0, 4*( N0-3 ), 4
427:                      IF( Z( I4 ).LE.TOL2*Z( I4-3 ) .OR.
428:      $                   Z( I4-1 ).LE.TOL2*SIGMA ) THEN
429:                         Z( I4-1 ) = -SIGMA
430:                         SPLT = I4 / 4
431:                         QMAX = ZERO
432:                         EMIN = Z( I4+3 )
433:                         OLDEMN = Z( I4+4 )
434:                      ELSE
435:                         QMAX = MAX( QMAX, Z( I4+1 ) )
436:                         EMIN = MIN( EMIN, Z( I4-1 ) )
437:                         OLDEMN = MIN( OLDEMN, Z( I4 ) )
438:                      END IF
439:   130             CONTINUE
440:                   Z( 4*N0-1 ) = EMIN
441:                   Z( 4*N0 ) = OLDEMN
442:                   I0 = SPLT + 1
443:                END IF
444:             END IF
445: *
446:   140    CONTINUE
447: *
448:          INFO = 2
449:          RETURN
450: *
451: *        end IWHILB
452: *
453:   150    CONTINUE
454: *
455:   160 CONTINUE
456: *
457:       INFO = 3
458:       RETURN
459: *
460: *     end IWHILA   
461: *
462:   170 CONTINUE
463: *      
464: *     Move q's to the front.
465: *      
466:       DO 180 K = 2, N
467:          Z( K ) = Z( 4*K-3 )
468:   180 CONTINUE
469: *      
470: *     Sort and compute sum of eigenvalues.
471: *
472:       CALL SLASRT( 'D', N, Z, IINFO )
473: *
474:       E = ZERO
475:       DO 190 K = N, 1, -1
476:          E = E + Z( K )
477:   190 CONTINUE
478: *
479: *     Store trace, sum(eigenvalues) and information on performance.
480: *
481:       Z( 2*N+1 ) = TRACE 
482:       Z( 2*N+2 ) = E
483:       Z( 2*N+3 ) = REAL( ITER )
484:       Z( 2*N+4 ) = REAL( NDIV ) / REAL( N**2 )
485:       Z( 2*N+5 ) = HUNDRD*NFAIL / REAL( ITER )
486:       RETURN
487: *
488: *     End of SLASQ2
489: *
490:       END
491: