Lots of trailing whitespaces in the files of Syd. Cleaning this. No big deal.
[platform/upstream/lapack.git] / SRC / slasq2.f
1 *> \brief \b SLASQ2 computes all the eigenvalues of the symmetric positive definite tridiagonal matrix associated with the qd Array Z to high relative accuracy. Used by sbdsqr and sstegr.
2 *
3 *  =========== DOCUMENTATION ===========
4 *
5 * Online html documentation available at
6 *            http://www.netlib.org/lapack/explore-html/
7 *
8 *> \htmlonly
9 *> Download SLASQ2 + dependencies
10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/slasq2.f">
11 *> [TGZ]</a>
12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/slasq2.f">
13 *> [ZIP]</a>
14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/slasq2.f">
15 *> [TXT]</a>
16 *> \endhtmlonly
17 *
18 *  Definition:
19 *  ===========
20 *
21 *       SUBROUTINE SLASQ2( N, Z, INFO )
22 *
23 *       .. Scalar Arguments ..
24 *       INTEGER            INFO, N
25 *       ..
26 *       .. Array Arguments ..
27 *       REAL               Z( * )
28 *       ..
29 *
30 *
31 *> \par Purpose:
32 *  =============
33 *>
34 *> \verbatim
35 *>
36 *> SLASQ2 computes all the eigenvalues of the symmetric positive
37 *> definite tridiagonal matrix associated with the qd array Z to high
38 *> relative accuracy are computed to high relative accuracy, in the
39 *> absence of denormalization, underflow and overflow.
40 *>
41 *> To see the relation of Z to the tridiagonal matrix, let L be a
42 *> unit lower bidiagonal matrix with subdiagonals Z(2,4,6,,..) and
43 *> let U be an upper bidiagonal matrix with 1's above and diagonal
44 *> Z(1,3,5,,..). The tridiagonal is L*U or, if you prefer, the
45 *> symmetric tridiagonal to which it is similar.
46 *>
47 *> Note : SLASQ2 defines a logical variable, IEEE, which is true
48 *> on machines which follow ieee-754 floating-point standard in their
49 *> handling of infinities and NaNs, and false otherwise. This variable
50 *> is passed to SLASQ3.
51 *> \endverbatim
52 *
53 *  Arguments:
54 *  ==========
55 *
56 *> \param[in] N
57 *> \verbatim
58 *>          N is INTEGER
59 *>        The number of rows and columns in the matrix. N >= 0.
60 *> \endverbatim
61 *>
62 *> \param[in,out] Z
63 *> \verbatim
64 *>          Z is REAL array, dimension ( 4*N )
65 *>        On entry Z holds the qd array. On exit, entries 1 to N hold
66 *>        the eigenvalues in decreasing order, Z( 2*N+1 ) holds the
67 *>        trace, and Z( 2*N+2 ) holds the sum of the eigenvalues. If
68 *>        N > 2, then Z( 2*N+3 ) holds the iteration count, Z( 2*N+4 )
69 *>        holds NDIVS/NIN^2, and Z( 2*N+5 ) holds the percentage of
70 *>        shifts that failed.
71 *> \endverbatim
72 *>
73 *> \param[out] INFO
74 *> \verbatim
75 *>          INFO is INTEGER
76 *>        = 0: successful exit
77 *>        < 0: if the i-th argument is a scalar and had an illegal
78 *>             value, then INFO = -i, if the i-th argument is an
79 *>             array and the j-entry had an illegal value, then
80 *>             INFO = -(i*100+j)
81 *>        > 0: the algorithm failed
82 *>              = 1, a split was marked by a positive value in E
83 *>              = 2, current block of Z not diagonalized after 100*N
84 *>                   iterations (in inner while loop).  On exit Z holds
85 *>                   a qd array with the same eigenvalues as the given Z.
86 *>              = 3, termination criterion of outer while loop not met
87 *>                   (program created more than N unreduced blocks)
88 *> \endverbatim
89 *
90 *  Authors:
91 *  ========
92 *
93 *> \author Univ. of Tennessee
94 *> \author Univ. of California Berkeley
95 *> \author Univ. of Colorado Denver
96 *> \author NAG Ltd.
97 *
98 *> \date September 2012
99 *
100 *> \ingroup auxOTHERcomputational
101 *
102 *> \par Further Details:
103 *  =====================
104 *>
105 *> \verbatim
106 *>
107 *>  Local Variables: I0:N0 defines a current unreduced segment of Z.
108 *>  The shifts are accumulated in SIGMA. Iteration count is in ITER.
109 *>  Ping-pong is controlled by PP (alternates between 0 and 1).
110 *> \endverbatim
111 *>
112 *  =====================================================================
113       SUBROUTINE SLASQ2( N, Z, INFO )
114 *
115 *  -- LAPACK computational routine (version 3.4.2) --
116 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
117 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
118 *     September 2012
119 *
120 *     .. Scalar Arguments ..
121       INTEGER            INFO, N
122 *     ..
123 *     .. Array Arguments ..
124       REAL               Z( * )
125 *     ..
126 *
127 *  =====================================================================
128 *
129 *     .. Parameters ..
130       REAL               CBIAS
131       PARAMETER          ( CBIAS = 1.50E0 )
132       REAL               ZERO, HALF, ONE, TWO, FOUR, HUNDRD
133       PARAMETER          ( ZERO = 0.0E0, HALF = 0.5E0, ONE = 1.0E0,
134      $                     TWO = 2.0E0, FOUR = 4.0E0, HUNDRD = 100.0E0 )
135 *     ..
136 *     .. Local Scalars ..
137       LOGICAL            IEEE
138       INTEGER            I0, I4, IINFO, IPN4, ITER, IWHILA, IWHILB, K,
139      $                   KMIN, N0, NBIG, NDIV, NFAIL, PP, SPLT, TTYPE,
140      $                   I1, N1
141       REAL               D, DEE, DEEMIN, DESIG, DMIN, DMIN1, DMIN2, DN,
142      $                   DN1, DN2, E, EMAX, EMIN, EPS, G, OLDEMN, QMAX,
143      $                   QMIN, S, SAFMIN, SIGMA, T, TAU, TEMP, TOL,
144      $                   TOL2, TRACE, ZMAX, TEMPE, TEMPQ
145 *     ..
146 *     .. External Subroutines ..
147       EXTERNAL           SLASQ3, SLASRT, XERBLA
148 *     ..
149 *     .. External Functions ..
150       REAL               SLAMCH
151       EXTERNAL           SLAMCH
152 *     ..
153 *     .. Intrinsic Functions ..
154       INTRINSIC          ABS, MAX, MIN, REAL, SQRT
155 *     ..
156 *     .. Executable Statements ..
157 *
158 *     Test the input arguments.
159 *     (in case SLASQ2 is not called by SLASQ1)
160 *
161       INFO = 0
162       EPS = SLAMCH( 'Precision' )
163       SAFMIN = SLAMCH( 'Safe minimum' )
164       TOL = EPS*HUNDRD
165       TOL2 = TOL**2
166 *
167       IF( N.LT.0 ) THEN
168          INFO = -1
169          CALL XERBLA( 'SLASQ2', 1 )
170          RETURN
171       ELSE IF( N.EQ.0 ) THEN
172          RETURN
173       ELSE IF( N.EQ.1 ) THEN
174 *
175 *        1-by-1 case.
176 *
177          IF( Z( 1 ).LT.ZERO ) THEN
178             INFO = -201
179             CALL XERBLA( 'SLASQ2', 2 )
180          END IF
181          RETURN
182       ELSE IF( N.EQ.2 ) THEN
183 *
184 *        2-by-2 case.
185 *
186          IF( Z( 2 ).LT.ZERO .OR. Z( 3 ).LT.ZERO ) THEN
187             INFO = -2
188             CALL XERBLA( 'SLASQ2', 2 )
189             RETURN
190          ELSE IF( Z( 3 ).GT.Z( 1 ) ) THEN
191             D = Z( 3 )
192             Z( 3 ) = Z( 1 )
193             Z( 1 ) = D
194          END IF
195          Z( 5 ) = Z( 1 ) + Z( 2 ) + Z( 3 )
196          IF( Z( 2 ).GT.Z( 3 )*TOL2 ) THEN
197             T = HALF*( ( Z( 1 )-Z( 3 ) )+Z( 2 ) )
198             S = Z( 3 )*( Z( 2 ) / T )
199             IF( S.LE.T ) THEN
200                S = Z( 3 )*( Z( 2 ) / ( T*( ONE+SQRT( ONE+S / T ) ) ) )
201             ELSE
202                S = Z( 3 )*( Z( 2 ) / ( T+SQRT( T )*SQRT( T+S ) ) )
203             END IF
204             T = Z( 1 ) + ( S+Z( 2 ) )
205             Z( 3 ) = Z( 3 )*( Z( 1 ) / T )
206             Z( 1 ) = T
207          END IF
208          Z( 2 ) = Z( 3 )
209          Z( 6 ) = Z( 2 ) + Z( 1 )
210          RETURN
211       END IF
212 *
213 *     Check for negative data and compute sums of q's and e's.
214 *
215       Z( 2*N ) = ZERO
216       EMIN = Z( 2 )
217       QMAX = ZERO
218       ZMAX = ZERO
219       D = ZERO
220       E = ZERO
221 *
222       DO 10 K = 1, 2*( N-1 ), 2
223          IF( Z( K ).LT.ZERO ) THEN
224             INFO = -( 200+K )
225             CALL XERBLA( 'SLASQ2', 2 )
226             RETURN
227          ELSE IF( Z( K+1 ).LT.ZERO ) THEN
228             INFO = -( 200+K+1 )
229             CALL XERBLA( 'SLASQ2', 2 )
230             RETURN
231          END IF
232          D = D + Z( K )
233          E = E + Z( K+1 )
234          QMAX = MAX( QMAX, Z( K ) )
235          EMIN = MIN( EMIN, Z( K+1 ) )
236          ZMAX = MAX( QMAX, ZMAX, Z( K+1 ) )
237    10 CONTINUE
238       IF( Z( 2*N-1 ).LT.ZERO ) THEN
239          INFO = -( 200+2*N-1 )
240          CALL XERBLA( 'SLASQ2', 2 )
241          RETURN
242       END IF
243       D = D + Z( 2*N-1 )
244       QMAX = MAX( QMAX, Z( 2*N-1 ) )
245       ZMAX = MAX( QMAX, ZMAX )
246 *
247 *     Check for diagonality.
248 *
249       IF( E.EQ.ZERO ) THEN
250          DO 20 K = 2, N
251             Z( K ) = Z( 2*K-1 )
252    20    CONTINUE
253          CALL SLASRT( 'D', N, Z, IINFO )
254          Z( 2*N-1 ) = D
255          RETURN
256       END IF
257 *
258       TRACE = D + E
259 *
260 *     Check for zero data.
261 *
262       IF( TRACE.EQ.ZERO ) THEN
263          Z( 2*N-1 ) = ZERO
264          RETURN
265       END IF
266 *
267 *     Check whether the machine is IEEE conformable.
268 *
269 *     IEEE = ILAENV( 10, 'SLASQ2', 'N', 1, 2, 3, 4 ).EQ.1 .AND.
270 *    $       ILAENV( 11, 'SLASQ2', 'N', 1, 2, 3, 4 ).EQ.1
271 *
272 *     [11/15/2008] The case IEEE=.TRUE. has a problem in single precision with
273 *     some the test matrices of type 16. The double precision code is fine.
274 *
275       IEEE = .FALSE.
276 *
277 *     Rearrange data for locality: Z=(q1,qq1,e1,ee1,q2,qq2,e2,ee2,...).
278 *
279       DO 30 K = 2*N, 2, -2
280          Z( 2*K ) = ZERO
281          Z( 2*K-1 ) = Z( K )
282          Z( 2*K-2 ) = ZERO
283          Z( 2*K-3 ) = Z( K-1 )
284    30 CONTINUE
285 *
286       I0 = 1
287       N0 = N
288 *
289 *     Reverse the qd-array, if warranted.
290 *
291       IF( CBIAS*Z( 4*I0-3 ).LT.Z( 4*N0-3 ) ) THEN
292          IPN4 = 4*( I0+N0 )
293          DO 40 I4 = 4*I0, 2*( I0+N0-1 ), 4
294             TEMP = Z( I4-3 )
295             Z( I4-3 ) = Z( IPN4-I4-3 )
296             Z( IPN4-I4-3 ) = TEMP
297             TEMP = Z( I4-1 )
298             Z( I4-1 ) = Z( IPN4-I4-5 )
299             Z( IPN4-I4-5 ) = TEMP
300    40    CONTINUE
301       END IF
302 *
303 *     Initial split checking via dqd and Li's test.
304 *
305       PP = 0
306 *
307       DO 80 K = 1, 2
308 *
309          D = Z( 4*N0+PP-3 )
310          DO 50 I4 = 4*( N0-1 ) + PP, 4*I0 + PP, -4
311             IF( Z( I4-1 ).LE.TOL2*D ) THEN
312                Z( I4-1 ) = -ZERO
313                D = Z( I4-3 )
314             ELSE
315                D = Z( I4-3 )*( D / ( D+Z( I4-1 ) ) )
316             END IF
317    50    CONTINUE
318 *
319 *        dqd maps Z to ZZ plus Li's test.
320 *
321          EMIN = Z( 4*I0+PP+1 )
322          D = Z( 4*I0+PP-3 )
323          DO 60 I4 = 4*I0 + PP, 4*( N0-1 ) + PP, 4
324             Z( I4-2*PP-2 ) = D + Z( I4-1 )
325             IF( Z( I4-1 ).LE.TOL2*D ) THEN
326                Z( I4-1 ) = -ZERO
327                Z( I4-2*PP-2 ) = D
328                Z( I4-2*PP ) = ZERO
329                D = Z( I4+1 )
330             ELSE IF( SAFMIN*Z( I4+1 ).LT.Z( I4-2*PP-2 ) .AND.
331      $               SAFMIN*Z( I4-2*PP-2 ).LT.Z( I4+1 ) ) THEN
332                TEMP = Z( I4+1 ) / Z( I4-2*PP-2 )
333                Z( I4-2*PP ) = Z( I4-1 )*TEMP
334                D = D*TEMP
335             ELSE
336                Z( I4-2*PP ) = Z( I4+1 )*( Z( I4-1 ) / Z( I4-2*PP-2 ) )
337                D = Z( I4+1 )*( D / Z( I4-2*PP-2 ) )
338             END IF
339             EMIN = MIN( EMIN, Z( I4-2*PP ) )
340    60    CONTINUE
341          Z( 4*N0-PP-2 ) = D
342 *
343 *        Now find qmax.
344 *
345          QMAX = Z( 4*I0-PP-2 )
346          DO 70 I4 = 4*I0 - PP + 2, 4*N0 - PP - 2, 4
347             QMAX = MAX( QMAX, Z( I4 ) )
348    70    CONTINUE
349 *
350 *        Prepare for the next iteration on K.
351 *
352          PP = 1 - PP
353    80 CONTINUE
354 *
355 *     Initialise variables to pass to SLASQ3.
356 *
357       TTYPE = 0
358       DMIN1 = ZERO
359       DMIN2 = ZERO
360       DN    = ZERO
361       DN1   = ZERO
362       DN2   = ZERO
363       G     = ZERO
364       TAU   = ZERO
365 *
366       ITER = 2
367       NFAIL = 0
368       NDIV = 2*( N0-I0 )
369 *
370       DO 160 IWHILA = 1, N + 1
371          IF( N0.LT.1 )
372      $      GO TO 170
373 *
374 *        While array unfinished do
375 *
376 *        E(N0) holds the value of SIGMA when submatrix in I0:N0
377 *        splits from the rest of the array, but is negated.
378 *
379          DESIG = ZERO
380          IF( N0.EQ.N ) THEN
381             SIGMA = ZERO
382          ELSE
383             SIGMA = -Z( 4*N0-1 )
384          END IF
385          IF( SIGMA.LT.ZERO ) THEN
386             INFO = 1
387             RETURN
388          END IF
389 *
390 *        Find last unreduced submatrix's top index I0, find QMAX and
391 *        EMIN. Find Gershgorin-type bound if Q's much greater than E's.
392 *
393          EMAX = ZERO
394          IF( N0.GT.I0 ) THEN
395             EMIN = ABS( Z( 4*N0-5 ) )
396          ELSE
397             EMIN = ZERO
398          END IF
399          QMIN = Z( 4*N0-3 )
400          QMAX = QMIN
401          DO 90 I4 = 4*N0, 8, -4
402             IF( Z( I4-5 ).LE.ZERO )
403      $         GO TO 100
404             IF( QMIN.GE.FOUR*EMAX ) THEN
405                QMIN = MIN( QMIN, Z( I4-3 ) )
406                EMAX = MAX( EMAX, Z( I4-5 ) )
407             END IF
408             QMAX = MAX( QMAX, Z( I4-7 )+Z( I4-5 ) )
409             EMIN = MIN( EMIN, Z( I4-5 ) )
410    90    CONTINUE
411          I4 = 4
412 *
413   100    CONTINUE
414          I0 = I4 / 4
415          PP = 0
416 *
417          IF( N0-I0.GT.1 ) THEN
418             DEE = Z( 4*I0-3 )
419             DEEMIN = DEE
420             KMIN = I0
421             DO 110 I4 = 4*I0+1, 4*N0-3, 4
422                DEE = Z( I4 )*( DEE /( DEE+Z( I4-2 ) ) )
423                IF( DEE.LE.DEEMIN ) THEN
424                   DEEMIN = DEE
425                   KMIN = ( I4+3 )/4
426                END IF
427   110       CONTINUE
428             IF( (KMIN-I0)*2.LT.N0-KMIN .AND.
429      $         DEEMIN.LE.HALF*Z(4*N0-3) ) THEN
430                IPN4 = 4*( I0+N0 )
431                PP = 2
432                DO 120 I4 = 4*I0, 2*( I0+N0-1 ), 4
433                   TEMP = Z( I4-3 )
434                   Z( I4-3 ) = Z( IPN4-I4-3 )
435                   Z( IPN4-I4-3 ) = TEMP
436                   TEMP = Z( I4-2 )
437                   Z( I4-2 ) = Z( IPN4-I4-2 )
438                   Z( IPN4-I4-2 ) = TEMP
439                   TEMP = Z( I4-1 )
440                   Z( I4-1 ) = Z( IPN4-I4-5 )
441                   Z( IPN4-I4-5 ) = TEMP
442                   TEMP = Z( I4 )
443                   Z( I4 ) = Z( IPN4-I4-4 )
444                   Z( IPN4-I4-4 ) = TEMP
445   120          CONTINUE
446             END IF
447          END IF
448 *
449 *        Put -(initial shift) into DMIN.
450 *
451          DMIN = -MAX( ZERO, QMIN-TWO*SQRT( QMIN )*SQRT( EMAX ) )
452 *
453 *        Now I0:N0 is unreduced.
454 *        PP = 0 for ping, PP = 1 for pong.
455 *        PP = 2 indicates that flipping was applied to the Z array and
456 *               and that the tests for deflation upon entry in SLASQ3
457 *               should not be performed.
458 *
459          NBIG = 100*( N0-I0+1 )
460          DO 140 IWHILB = 1, NBIG
461             IF( I0.GT.N0 )
462      $         GO TO 150
463 *
464 *           While submatrix unfinished take a good dqds step.
465 *
466             CALL SLASQ3( I0, N0, Z, PP, DMIN, SIGMA, DESIG, QMAX, NFAIL,
467      $                   ITER, NDIV, IEEE, TTYPE, DMIN1, DMIN2, DN, DN1,
468      $                   DN2, G, TAU )
469 *
470             PP = 1 - PP
471 *
472 *           When EMIN is very small check for splits.
473 *
474             IF( PP.EQ.0 .AND. N0-I0.GE.3 ) THEN
475                IF( Z( 4*N0 ).LE.TOL2*QMAX .OR.
476      $             Z( 4*N0-1 ).LE.TOL2*SIGMA ) THEN
477                   SPLT = I0 - 1
478                   QMAX = Z( 4*I0-3 )
479                   EMIN = Z( 4*I0-1 )
480                   OLDEMN = Z( 4*I0 )
481                   DO 130 I4 = 4*I0, 4*( N0-3 ), 4
482                      IF( Z( I4 ).LE.TOL2*Z( I4-3 ) .OR.
483      $                   Z( I4-1 ).LE.TOL2*SIGMA ) THEN
484                         Z( I4-1 ) = -SIGMA
485                         SPLT = I4 / 4
486                         QMAX = ZERO
487                         EMIN = Z( I4+3 )
488                         OLDEMN = Z( I4+4 )
489                      ELSE
490                         QMAX = MAX( QMAX, Z( I4+1 ) )
491                         EMIN = MIN( EMIN, Z( I4-1 ) )
492                         OLDEMN = MIN( OLDEMN, Z( I4 ) )
493                      END IF
494   130             CONTINUE
495                   Z( 4*N0-1 ) = EMIN
496                   Z( 4*N0 ) = OLDEMN
497                   I0 = SPLT + 1
498                END IF
499             END IF
500 *
501   140    CONTINUE
502 *
503          INFO = 2
504 *
505 *        Maximum number of iterations exceeded, restore the shift
506 *        SIGMA and place the new d's and e's in a qd array.
507 *        This might need to be done for several blocks
508 *
509          I1 = I0
510          N1 = N0
511  145     CONTINUE
512          TEMPQ = Z( 4*I0-3 )
513          Z( 4*I0-3 ) = Z( 4*I0-3 ) + SIGMA
514          DO K = I0+1, N0
515             TEMPE = Z( 4*K-5 )
516             Z( 4*K-5 ) = Z( 4*K-5 ) * (TEMPQ / Z( 4*K-7 ))
517             TEMPQ = Z( 4*K-3 )
518             Z( 4*K-3 ) = Z( 4*K-3 ) + SIGMA + TEMPE - Z( 4*K-5 )
519          END DO
520 *
521 *        Prepare to do this on the previous block if there is one
522 *
523          IF( I1.GT.1 ) THEN
524             N1 = I1-1
525             DO WHILE( ( I1.GE.2 ) .AND. ( Z(4*I1-5).GE.ZERO ) )
526                I1 = I1 - 1
527             END DO
528             IF( I1.GE.1 ) THEN
529                SIGMA = -Z(4*N1-1)
530                GO TO 145
531             END IF
532          END IF
533
534          DO K = 1, N
535             Z( 2*K-1 ) = Z( 4*K-3 )
536 *
537 *        Only the block 1..N0 is unfinished.  The rest of the e's
538 *        must be essentially zero, although sometimes other data
539 *        has been stored in them.
540 *
541             IF( K.LT.N0 ) THEN
542                Z( 2*K ) = Z( 4*K-1 )
543             ELSE
544                Z( 2*K ) = 0
545             END IF
546          END DO
547          RETURN
548 *
549 *        end IWHILB
550 *
551   150    CONTINUE
552 *
553   160 CONTINUE
554 *
555       INFO = 3
556       RETURN
557 *
558 *     end IWHILA
559 *
560   170 CONTINUE
561 *
562 *     Move q's to the front.
563 *
564       DO 180 K = 2, N
565          Z( K ) = Z( 4*K-3 )
566   180 CONTINUE
567 *
568 *     Sort and compute sum of eigenvalues.
569 *
570       CALL SLASRT( 'D', N, Z, IINFO )
571 *
572       E = ZERO
573       DO 190 K = N, 1, -1
574          E = E + Z( K )
575   190 CONTINUE
576 *
577 *     Store trace, sum(eigenvalues) and information on performance.
578 *
579       Z( 2*N+1 ) = TRACE
580       Z( 2*N+2 ) = E
581       Z( 2*N+3 ) = REAL( ITER )
582       Z( 2*N+4 ) = REAL( NDIV ) / REAL( N**2 )
583       Z( 2*N+5 ) = HUNDRD*NFAIL / REAL( ITER )
584       RETURN
585 *
586 *     End of SLASQ2
587 *
588       END