df1690d020d282c25d12d575f94a22ccef3cd3a6
[platform/upstream/lapack.git] / SRC / dlasq2.f
1 *> \brief \b DLASQ2 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 DLASQ2 + dependencies 
10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dlasq2.f"> 
11 *> [TGZ]</a> 
12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dlasq2.f"> 
13 *> [ZIP]</a> 
14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dlasq2.f"> 
15 *> [TXT]</a>
16 *> \endhtmlonly 
17 *
18 *  Definition:
19 *  ===========
20 *
21 *       SUBROUTINE DLASQ2( N, Z, INFO )
22
23 *       .. Scalar Arguments ..
24 *       INTEGER            INFO, N
25 *       ..
26 *       .. Array Arguments ..
27 *       DOUBLE PRECISION   Z( * )
28 *       ..
29 *  
30 *
31 *> \par Purpose:
32 *  =============
33 *>
34 *> \verbatim
35 *>
36 *> DLASQ2 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 : DLASQ2 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 DLASQ3.
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 DOUBLE PRECISION 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 DLASQ2( 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       DOUBLE PRECISION   Z( * )
125 *     ..
126 *
127 *  =====================================================================
128 *
129 *     .. Parameters ..
130       DOUBLE PRECISION   CBIAS
131       PARAMETER          ( CBIAS = 1.50D0 )
132       DOUBLE PRECISION   ZERO, HALF, ONE, TWO, FOUR, HUNDRD
133       PARAMETER          ( ZERO = 0.0D0, HALF = 0.5D0, ONE = 1.0D0,
134      $                     TWO = 2.0D0, FOUR = 4.0D0, HUNDRD = 100.0D0 )
135 *     ..
136 *     .. Local Scalars ..
137       LOGICAL            IEEE
138       INTEGER            I0, I1, I4, IINFO, IPN4, ITER, IWHILA, IWHILB,
139      $                   K, KMIN, N0, N1, NBIG, NDIV, NFAIL, PP, SPLT, 
140      $                   TTYPE
141       DOUBLE PRECISION   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           DLASQ3, DLASRT, XERBLA
148 *     ..
149 *     .. External Functions ..
150       INTEGER            ILAENV
151       DOUBLE PRECISION   DLAMCH
152       EXTERNAL           DLAMCH, ILAENV
153 *     ..
154 *     .. Intrinsic Functions ..
155       INTRINSIC          ABS, DBLE, MAX, MIN, SQRT
156 *     ..
157 *     .. Executable Statements ..
158 *      
159 *     Test the input arguments.
160 *     (in case DLASQ2 is not called by DLASQ1)
161 *
162       INFO = 0
163       EPS = DLAMCH( 'Precision' )
164       SAFMIN = DLAMCH( 'Safe minimum' )
165       TOL = EPS*HUNDRD
166       TOL2 = TOL**2
167 *
168       IF( N.LT.0 ) THEN
169          INFO = -1
170          CALL XERBLA( 'DLASQ2', 1 )
171          RETURN
172       ELSE IF( N.EQ.0 ) THEN
173          RETURN
174       ELSE IF( N.EQ.1 ) THEN
175 *
176 *        1-by-1 case.
177 *
178          IF( Z( 1 ).LT.ZERO ) THEN
179             INFO = -201
180             CALL XERBLA( 'DLASQ2', 2 )
181          END IF
182          RETURN
183       ELSE IF( N.EQ.2 ) THEN
184 *
185 *        2-by-2 case.
186 *
187          IF( Z( 2 ).LT.ZERO .OR. Z( 3 ).LT.ZERO ) THEN
188             INFO = -2
189             CALL XERBLA( 'DLASQ2', 2 )
190             RETURN
191          ELSE IF( Z( 3 ).GT.Z( 1 ) ) THEN
192             D = Z( 3 )
193             Z( 3 ) = Z( 1 )
194             Z( 1 ) = D
195          END IF
196          Z( 5 ) = Z( 1 ) + Z( 2 ) + Z( 3 )
197          IF( Z( 2 ).GT.Z( 3 )*TOL2 ) THEN
198             T = HALF*( ( Z( 1 )-Z( 3 ) )+Z( 2 ) ) 
199             S = Z( 3 )*( Z( 2 ) / T )
200             IF( S.LE.T ) THEN
201                S = Z( 3 )*( Z( 2 ) / ( T*( ONE+SQRT( ONE+S / T ) ) ) )
202             ELSE
203                S = Z( 3 )*( Z( 2 ) / ( T+SQRT( T )*SQRT( T+S ) ) )
204             END IF
205             T = Z( 1 ) + ( S+Z( 2 ) )
206             Z( 3 ) = Z( 3 )*( Z( 1 ) / T )
207             Z( 1 ) = T
208          END IF
209          Z( 2 ) = Z( 3 )
210          Z( 6 ) = Z( 2 ) + Z( 1 )
211          RETURN
212       END IF
213 *
214 *     Check for negative data and compute sums of q's and e's.
215 *
216       Z( 2*N ) = ZERO
217       EMIN = Z( 2 )
218       QMAX = ZERO
219       ZMAX = ZERO
220       D = ZERO
221       E = ZERO
222 *
223       DO 10 K = 1, 2*( N-1 ), 2
224          IF( Z( K ).LT.ZERO ) THEN
225             INFO = -( 200+K )
226             CALL XERBLA( 'DLASQ2', 2 )
227             RETURN
228          ELSE IF( Z( K+1 ).LT.ZERO ) THEN
229             INFO = -( 200+K+1 )
230             CALL XERBLA( 'DLASQ2', 2 )
231             RETURN
232          END IF
233          D = D + Z( K )
234          E = E + Z( K+1 )
235          QMAX = MAX( QMAX, Z( K ) )
236          EMIN = MIN( EMIN, Z( K+1 ) )
237          ZMAX = MAX( QMAX, ZMAX, Z( K+1 ) )
238    10 CONTINUE
239       IF( Z( 2*N-1 ).LT.ZERO ) THEN
240          INFO = -( 200+2*N-1 )
241          CALL XERBLA( 'DLASQ2', 2 )
242          RETURN
243       END IF
244       D = D + Z( 2*N-1 )
245       QMAX = MAX( QMAX, Z( 2*N-1 ) )
246       ZMAX = MAX( QMAX, ZMAX )
247 *
248 *     Check for diagonality.
249 *
250       IF( E.EQ.ZERO ) THEN
251          DO 20 K = 2, N
252             Z( K ) = Z( 2*K-1 )
253    20    CONTINUE
254          CALL DLASRT( 'D', N, Z, IINFO )
255          Z( 2*N-1 ) = D
256          RETURN
257       END IF
258 *
259       TRACE = D + E
260 *
261 *     Check for zero data.
262 *
263       IF( TRACE.EQ.ZERO ) THEN
264          Z( 2*N-1 ) = ZERO
265          RETURN
266       END IF
267 *         
268 *     Check whether the machine is IEEE conformable.
269 *         
270       IEEE = ILAENV( 10, 'DLASQ2', 'N', 1, 2, 3, 4 ).EQ.1 .AND.
271      $       ILAENV( 11, 'DLASQ2', 'N', 1, 2, 3, 4 ).EQ.1      
272 *         
273 *     Rearrange data for locality: Z=(q1,qq1,e1,ee1,q2,qq2,e2,ee2,...).
274 *
275       DO 30 K = 2*N, 2, -2
276          Z( 2*K ) = ZERO 
277          Z( 2*K-1 ) = Z( K ) 
278          Z( 2*K-2 ) = ZERO 
279          Z( 2*K-3 ) = Z( K-1 ) 
280    30 CONTINUE
281 *
282       I0 = 1
283       N0 = N
284 *
285 *     Reverse the qd-array, if warranted.
286 *
287       IF( CBIAS*Z( 4*I0-3 ).LT.Z( 4*N0-3 ) ) THEN
288          IPN4 = 4*( I0+N0 )
289          DO 40 I4 = 4*I0, 2*( I0+N0-1 ), 4
290             TEMP = Z( I4-3 )
291             Z( I4-3 ) = Z( IPN4-I4-3 )
292             Z( IPN4-I4-3 ) = TEMP
293             TEMP = Z( I4-1 )
294             Z( I4-1 ) = Z( IPN4-I4-5 )
295             Z( IPN4-I4-5 ) = TEMP
296    40    CONTINUE
297       END IF
298 *
299 *     Initial split checking via dqd and Li's test.
300 *
301       PP = 0
302 *
303       DO 80 K = 1, 2
304 *
305          D = Z( 4*N0+PP-3 )
306          DO 50 I4 = 4*( N0-1 ) + PP, 4*I0 + PP, -4
307             IF( Z( I4-1 ).LE.TOL2*D ) THEN
308                Z( I4-1 ) = -ZERO
309                D = Z( I4-3 )
310             ELSE
311                D = Z( I4-3 )*( D / ( D+Z( I4-1 ) ) )
312             END IF
313    50    CONTINUE
314 *
315 *        dqd maps Z to ZZ plus Li's test.
316 *
317          EMIN = Z( 4*I0+PP+1 )
318          D = Z( 4*I0+PP-3 )
319          DO 60 I4 = 4*I0 + PP, 4*( N0-1 ) + PP, 4
320             Z( I4-2*PP-2 ) = D + Z( I4-1 )
321             IF( Z( I4-1 ).LE.TOL2*D ) THEN
322                Z( I4-1 ) = -ZERO
323                Z( I4-2*PP-2 ) = D
324                Z( I4-2*PP ) = ZERO
325                D = Z( I4+1 )
326             ELSE IF( SAFMIN*Z( I4+1 ).LT.Z( I4-2*PP-2 ) .AND.
327      $               SAFMIN*Z( I4-2*PP-2 ).LT.Z( I4+1 ) ) THEN
328                TEMP = Z( I4+1 ) / Z( I4-2*PP-2 )
329                Z( I4-2*PP ) = Z( I4-1 )*TEMP
330                D = D*TEMP
331             ELSE
332                Z( I4-2*PP ) = Z( I4+1 )*( Z( I4-1 ) / Z( I4-2*PP-2 ) )
333                D = Z( I4+1 )*( D / Z( I4-2*PP-2 ) )
334             END IF
335             EMIN = MIN( EMIN, Z( I4-2*PP ) )
336    60    CONTINUE 
337          Z( 4*N0-PP-2 ) = D
338 *
339 *        Now find qmax.
340 *
341          QMAX = Z( 4*I0-PP-2 )
342          DO 70 I4 = 4*I0 - PP + 2, 4*N0 - PP - 2, 4
343             QMAX = MAX( QMAX, Z( I4 ) )
344    70    CONTINUE
345 *
346 *        Prepare for the next iteration on K.
347 *
348          PP = 1 - PP
349    80 CONTINUE
350 *
351 *     Initialise variables to pass to DLASQ3.
352 *
353       TTYPE = 0
354       DMIN1 = ZERO
355       DMIN2 = ZERO
356       DN    = ZERO
357       DN1   = ZERO
358       DN2   = ZERO
359       G     = ZERO
360       TAU   = ZERO
361 *
362       ITER = 2
363       NFAIL = 0
364       NDIV = 2*( N0-I0 )
365 *
366       DO 160 IWHILA = 1, N + 1
367          IF( N0.LT.1 ) 
368      $      GO TO 170
369 *
370 *        While array unfinished do 
371 *
372 *        E(N0) holds the value of SIGMA when submatrix in I0:N0
373 *        splits from the rest of the array, but is negated.
374 *      
375          DESIG = ZERO
376          IF( N0.EQ.N ) THEN
377             SIGMA = ZERO
378          ELSE
379             SIGMA = -Z( 4*N0-1 )
380          END IF
381          IF( SIGMA.LT.ZERO ) THEN
382             INFO = 1
383             RETURN
384          END IF
385 *
386 *        Find last unreduced submatrix's top index I0, find QMAX and
387 *        EMIN. Find Gershgorin-type bound if Q's much greater than E's.
388 *
389          EMAX = ZERO 
390          IF( N0.GT.I0 ) THEN
391             EMIN = ABS( Z( 4*N0-5 ) )
392          ELSE
393             EMIN = ZERO
394          END IF
395          QMIN = Z( 4*N0-3 )
396          QMAX = QMIN
397          DO 90 I4 = 4*N0, 8, -4
398             IF( Z( I4-5 ).LE.ZERO )
399      $         GO TO 100
400             IF( QMIN.GE.FOUR*EMAX ) THEN
401                QMIN = MIN( QMIN, Z( I4-3 ) )
402                EMAX = MAX( EMAX, Z( I4-5 ) )
403             END IF
404             QMAX = MAX( QMAX, Z( I4-7 )+Z( I4-5 ) )
405             EMIN = MIN( EMIN, Z( I4-5 ) )
406    90    CONTINUE
407          I4 = 4 
408 *
409   100    CONTINUE
410          I0 = I4 / 4
411          PP = 0
412 *
413          IF( N0-I0.GT.1 ) THEN
414             DEE = Z( 4*I0-3 )
415             DEEMIN = DEE
416             KMIN = I0
417             DO 110 I4 = 4*I0+1, 4*N0-3, 4
418                DEE = Z( I4 )*( DEE /( DEE+Z( I4-2 ) ) )
419                IF( DEE.LE.DEEMIN ) THEN
420                   DEEMIN = DEE
421                   KMIN = ( I4+3 )/4
422                END IF
423   110       CONTINUE
424             IF( (KMIN-I0)*2.LT.N0-KMIN .AND. 
425      $         DEEMIN.LE.HALF*Z(4*N0-3) ) THEN
426                IPN4 = 4*( I0+N0 )
427                PP = 2
428                DO 120 I4 = 4*I0, 2*( I0+N0-1 ), 4
429                   TEMP = Z( I4-3 )
430                   Z( I4-3 ) = Z( IPN4-I4-3 )
431                   Z( IPN4-I4-3 ) = TEMP
432                   TEMP = Z( I4-2 )
433                   Z( I4-2 ) = Z( IPN4-I4-2 )
434                   Z( IPN4-I4-2 ) = TEMP
435                   TEMP = Z( I4-1 )
436                   Z( I4-1 ) = Z( IPN4-I4-5 )
437                   Z( IPN4-I4-5 ) = TEMP
438                   TEMP = Z( I4 )
439                   Z( I4 ) = Z( IPN4-I4-4 )
440                   Z( IPN4-I4-4 ) = TEMP
441   120          CONTINUE
442             END IF
443          END IF
444 *
445 *        Put -(initial shift) into DMIN.
446 *
447          DMIN = -MAX( ZERO, QMIN-TWO*SQRT( QMIN )*SQRT( EMAX ) )
448 *
449 *        Now I0:N0 is unreduced. 
450 *        PP = 0 for ping, PP = 1 for pong.
451 *        PP = 2 indicates that flipping was applied to the Z array and
452 *               and that the tests for deflation upon entry in DLASQ3 
453 *               should not be performed.
454 *
455          NBIG = 100*( N0-I0+1 )
456          DO 140 IWHILB = 1, NBIG
457             IF( I0.GT.N0 ) 
458      $         GO TO 150
459 *
460 *           While submatrix unfinished take a good dqds step.
461 *
462             CALL DLASQ3( I0, N0, Z, PP, DMIN, SIGMA, DESIG, QMAX, NFAIL,
463      $                   ITER, NDIV, IEEE, TTYPE, DMIN1, DMIN2, DN, DN1,
464      $                   DN2, G, TAU )
465 *
466             PP = 1 - PP
467 *
468 *           When EMIN is very small check for splits.
469 *
470             IF( PP.EQ.0 .AND. N0-I0.GE.3 ) THEN
471                IF( Z( 4*N0 ).LE.TOL2*QMAX .OR.
472      $             Z( 4*N0-1 ).LE.TOL2*SIGMA ) THEN
473                   SPLT = I0 - 1
474                   QMAX = Z( 4*I0-3 )
475                   EMIN = Z( 4*I0-1 )
476                   OLDEMN = Z( 4*I0 )
477                   DO 130 I4 = 4*I0, 4*( N0-3 ), 4
478                      IF( Z( I4 ).LE.TOL2*Z( I4-3 ) .OR.
479      $                   Z( I4-1 ).LE.TOL2*SIGMA ) THEN
480                         Z( I4-1 ) = -SIGMA
481                         SPLT = I4 / 4
482                         QMAX = ZERO
483                         EMIN = Z( I4+3 )
484                         OLDEMN = Z( I4+4 )
485                      ELSE
486                         QMAX = MAX( QMAX, Z( I4+1 ) )
487                         EMIN = MIN( EMIN, Z( I4-1 ) )
488                         OLDEMN = MIN( OLDEMN, Z( I4 ) )
489                      END IF
490   130             CONTINUE
491                   Z( 4*N0-1 ) = EMIN
492                   Z( 4*N0 ) = OLDEMN
493                   I0 = SPLT + 1
494                END IF
495             END IF
496 *
497   140    CONTINUE
498 *
499          INFO = 2
500 *       
501 *        Maximum number of iterations exceeded, restore the shift 
502 *        SIGMA and place the new d's and e's in a qd array.
503 *        This might need to be done for several blocks
504 *
505          I1 = I0
506          N1 = N0
507  145     CONTINUE
508          TEMPQ = Z( 4*I0-3 )
509          Z( 4*I0-3 ) = Z( 4*I0-3 ) + SIGMA
510          DO K = I0+1, N0
511             TEMPE = Z( 4*K-5 )
512             Z( 4*K-5 ) = Z( 4*K-5 ) * (TEMPQ / Z( 4*K-7 ))
513             TEMPQ = Z( 4*K-3 )
514             Z( 4*K-3 ) = Z( 4*K-3 ) + SIGMA + TEMPE - Z( 4*K-5 )
515          END DO
516 *
517 *        Prepare to do this on the previous block if there is one
518 *
519          IF( I1.GT.1 ) THEN
520             N1 = I1-1
521             DO WHILE( ( I1.GE.2 ) .AND. ( Z(4*I1-5).GE.ZERO ) )
522                I1 = I1 - 1
523             END DO
524             SIGMA = -Z(4*N1-1)
525             GO TO 145
526          END IF
527
528          DO K = 1, N
529             Z( 2*K-1 ) = Z( 4*K-3 )
530 *
531 *        Only the block 1..N0 is unfinished.  The rest of the e's
532 *        must be essentially zero, although sometimes other data
533 *        has been stored in them.
534 *
535             IF( K.LT.N0 ) THEN
536                Z( 2*K ) = Z( 4*K-1 )
537             ELSE
538                Z( 2*K ) = 0
539             END IF
540          END DO
541          RETURN
542 *
543 *        end IWHILB
544 *
545   150    CONTINUE
546 *
547   160 CONTINUE
548 *
549       INFO = 3
550       RETURN
551 *
552 *     end IWHILA   
553 *
554   170 CONTINUE
555 *      
556 *     Move q's to the front.
557 *      
558       DO 180 K = 2, N
559          Z( K ) = Z( 4*K-3 )
560   180 CONTINUE
561 *      
562 *     Sort and compute sum of eigenvalues.
563 *
564       CALL DLASRT( 'D', N, Z, IINFO )
565 *
566       E = ZERO
567       DO 190 K = N, 1, -1
568          E = E + Z( K )
569   190 CONTINUE
570 *
571 *     Store trace, sum(eigenvalues) and information on performance.
572 *
573       Z( 2*N+1 ) = TRACE 
574       Z( 2*N+2 ) = E
575       Z( 2*N+3 ) = DBLE( ITER )
576       Z( 2*N+4 ) = DBLE( NDIV ) / DBLE( N**2 )
577       Z( 2*N+5 ) = HUNDRD*NFAIL / DBLE( ITER )
578       RETURN
579 *
580 *     End of DLASQ2
581 *
582       END