e554f8a0376d464f4d05cf708d5cc8f3720a80fd
[platform/upstream/lapack.git] / SRC / zlahqr.f
1 *> \brief \b ZLAHQR computes the eigenvalues and Schur factorization of an upper Hessenberg matrix, using the double-shift/single-shift QR algorithm.
2 *
3 *  =========== DOCUMENTATION ===========
4 *
5 * Online html documentation available at 
6 *            http://www.netlib.org/lapack/explore-html/ 
7 *
8 *> \htmlonly
9 *> Download ZLAHQR + dependencies 
10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zlahqr.f"> 
11 *> [TGZ]</a> 
12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zlahqr.f"> 
13 *> [ZIP]</a> 
14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zlahqr.f"> 
15 *> [TXT]</a>
16 *> \endhtmlonly 
17 *
18 *  Definition:
19 *  ===========
20 *
21 *       SUBROUTINE ZLAHQR( WANTT, WANTZ, N, ILO, IHI, H, LDH, W, ILOZ,
22 *                          IHIZ, Z, LDZ, INFO )
23
24 *       .. Scalar Arguments ..
25 *       INTEGER            IHI, IHIZ, ILO, ILOZ, INFO, LDH, LDZ, N
26 *       LOGICAL            WANTT, WANTZ
27 *       ..
28 *       .. Array Arguments ..
29 *       COMPLEX*16         H( LDH, * ), W( * ), Z( LDZ, * )
30 *       ..
31 *  
32 *
33 *> \par Purpose:
34 *  =============
35 *>
36 *> \verbatim
37 *>
38 *>    ZLAHQR is an auxiliary routine called by CHSEQR to update the
39 *>    eigenvalues and Schur decomposition already computed by CHSEQR, by
40 *>    dealing with the Hessenberg submatrix in rows and columns ILO to
41 *>    IHI.
42 *> \endverbatim
43 *
44 *  Arguments:
45 *  ==========
46 *
47 *> \param[in] WANTT
48 *> \verbatim
49 *>          WANTT is LOGICAL
50 *>          = .TRUE. : the full Schur form T is required;
51 *>          = .FALSE.: only eigenvalues are required.
52 *> \endverbatim
53 *>
54 *> \param[in] WANTZ
55 *> \verbatim
56 *>          WANTZ is LOGICAL
57 *>          = .TRUE. : the matrix of Schur vectors Z is required;
58 *>          = .FALSE.: Schur vectors are not required.
59 *> \endverbatim
60 *>
61 *> \param[in] N
62 *> \verbatim
63 *>          N is INTEGER
64 *>          The order of the matrix H.  N >= 0.
65 *> \endverbatim
66 *>
67 *> \param[in] ILO
68 *> \verbatim
69 *>          ILO is INTEGER
70 *> \endverbatim
71 *>
72 *> \param[in] IHI
73 *> \verbatim
74 *>          IHI is INTEGER
75 *>          It is assumed that H is already upper triangular in rows and
76 *>          columns IHI+1:N, and that H(ILO,ILO-1) = 0 (unless ILO = 1).
77 *>          ZLAHQR works primarily with the Hessenberg submatrix in rows
78 *>          and columns ILO to IHI, but applies transformations to all of
79 *>          H if WANTT is .TRUE..
80 *>          1 <= ILO <= max(1,IHI); IHI <= N.
81 *> \endverbatim
82 *>
83 *> \param[in,out] H
84 *> \verbatim
85 *>          H is COMPLEX*16 array, dimension (LDH,N)
86 *>          On entry, the upper Hessenberg matrix H.
87 *>          On exit, if INFO is zero and if WANTT is .TRUE., then H
88 *>          is upper triangular in rows and columns ILO:IHI.  If INFO
89 *>          is zero and if WANTT is .FALSE., then the contents of H
90 *>          are unspecified on exit.  The output state of H in case
91 *>          INF is positive is below under the description of INFO.
92 *> \endverbatim
93 *>
94 *> \param[in] LDH
95 *> \verbatim
96 *>          LDH is INTEGER
97 *>          The leading dimension of the array H. LDH >= max(1,N).
98 *> \endverbatim
99 *>
100 *> \param[out] W
101 *> \verbatim
102 *>          W is COMPLEX*16 array, dimension (N)
103 *>          The computed eigenvalues ILO to IHI are stored in the
104 *>          corresponding elements of W. If WANTT is .TRUE., the
105 *>          eigenvalues are stored in the same order as on the diagonal
106 *>          of the Schur form returned in H, with W(i) = H(i,i).
107 *> \endverbatim
108 *>
109 *> \param[in] ILOZ
110 *> \verbatim
111 *>          ILOZ is INTEGER
112 *> \endverbatim
113 *>
114 *> \param[in] IHIZ
115 *> \verbatim
116 *>          IHIZ is INTEGER
117 *>          Specify the rows of Z to which transformations must be
118 *>          applied if WANTZ is .TRUE..
119 *>          1 <= ILOZ <= ILO; IHI <= IHIZ <= N.
120 *> \endverbatim
121 *>
122 *> \param[in,out] Z
123 *> \verbatim
124 *>          Z is COMPLEX*16 array, dimension (LDZ,N)
125 *>          If WANTZ is .TRUE., on entry Z must contain the current
126 *>          matrix Z of transformations accumulated by CHSEQR, and on
127 *>          exit Z has been updated; transformations are applied only to
128 *>          the submatrix Z(ILOZ:IHIZ,ILO:IHI).
129 *>          If WANTZ is .FALSE., Z is not referenced.
130 *> \endverbatim
131 *>
132 *> \param[in] LDZ
133 *> \verbatim
134 *>          LDZ is INTEGER
135 *>          The leading dimension of the array Z. LDZ >= max(1,N).
136 *> \endverbatim
137 *>
138 *> \param[out] INFO
139 *> \verbatim
140 *>          INFO is INTEGER
141 *>           =   0: successful exit
142 *>          .GT. 0: if INFO = i, ZLAHQR failed to compute all the
143 *>                  eigenvalues ILO to IHI in a total of 30 iterations
144 *>                  per eigenvalue; elements i+1:ihi of W contain
145 *>                  those eigenvalues which have been successfully
146 *>                  computed.
147 *>
148 *>                  If INFO .GT. 0 and WANTT is .FALSE., then on exit,
149 *>                  the remaining unconverged eigenvalues are the
150 *>                  eigenvalues of the upper Hessenberg matrix
151 *>                  rows and columns ILO thorugh INFO of the final,
152 *>                  output value of H.
153 *>
154 *>                  If INFO .GT. 0 and WANTT is .TRUE., then on exit
155 *>          (*)       (initial value of H)*U  = U*(final value of H)
156 *>                  where U is an orthognal matrix.    The final
157 *>                  value of H is upper Hessenberg and triangular in
158 *>                  rows and columns INFO+1 through IHI.
159 *>
160 *>                  If INFO .GT. 0 and WANTZ is .TRUE., then on exit
161 *>                      (final value of Z)  = (initial value of Z)*U
162 *>                  where U is the orthogonal matrix in (*)
163 *>                  (regardless of the value of WANTT.)
164 *> \endverbatim
165 *
166 *  Authors:
167 *  ========
168 *
169 *> \author Univ. of Tennessee 
170 *> \author Univ. of California Berkeley 
171 *> \author Univ. of Colorado Denver 
172 *> \author NAG Ltd. 
173 *
174 *> \date November 2015
175 *
176 *> \ingroup complex16OTHERauxiliary
177 *
178 *> \par Contributors:
179 *  ==================
180 *>
181 *> \verbatim
182 *>
183 *>     02-96 Based on modifications by
184 *>     David Day, Sandia National Laboratory, USA
185 *>
186 *>     12-04 Further modifications by
187 *>     Ralph Byers, University of Kansas, USA
188 *>     This is a modified version of ZLAHQR from LAPACK version 3.0.
189 *>     It is (1) more robust against overflow and underflow and
190 *>     (2) adopts the more conservative Ahues & Tisseur stopping
191 *>     criterion (LAWN 122, 1997).
192 *> \endverbatim
193 *>
194 *  =====================================================================
195       SUBROUTINE ZLAHQR( WANTT, WANTZ, N, ILO, IHI, H, LDH, W, ILOZ,
196      $                   IHIZ, Z, LDZ, INFO )
197 *
198 *  -- LAPACK auxiliary routine (version 3.6.0) --
199 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
200 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
201 *     November 2015
202 *
203 *     .. Scalar Arguments ..
204       INTEGER            IHI, IHIZ, ILO, ILOZ, INFO, LDH, LDZ, N
205       LOGICAL            WANTT, WANTZ
206 *     ..
207 *     .. Array Arguments ..
208       COMPLEX*16         H( LDH, * ), W( * ), Z( LDZ, * )
209 *     ..
210 *
211 *  =========================================================
212 *
213 *     .. Parameters ..
214       COMPLEX*16         ZERO, ONE
215       PARAMETER          ( ZERO = ( 0.0d0, 0.0d0 ),
216      $                   ONE = ( 1.0d0, 0.0d0 ) )
217       DOUBLE PRECISION   RZERO, RONE, HALF
218       PARAMETER          ( RZERO = 0.0d0, RONE = 1.0d0, HALF = 0.5d0 )
219       DOUBLE PRECISION   DAT1
220       PARAMETER          ( DAT1 = 3.0d0 / 4.0d0 )
221 *     ..
222 *     .. Local Scalars ..
223       COMPLEX*16         CDUM, H11, H11S, H22, SC, SUM, T, T1, TEMP, U,
224      $                   V2, X, Y
225       DOUBLE PRECISION   AA, AB, BA, BB, H10, H21, RTEMP, S, SAFMAX,
226      $                   SAFMIN, SMLNUM, SX, T2, TST, ULP
227       INTEGER            I, I1, I2, ITS, ITMAX, J, JHI, JLO, K, L, M,
228      $                   NH, NZ
229 *     ..
230 *     .. Local Arrays ..
231       COMPLEX*16         V( 2 )
232 *     ..
233 *     .. External Functions ..
234       COMPLEX*16         ZLADIV
235       DOUBLE PRECISION   DLAMCH
236       EXTERNAL           ZLADIV, DLAMCH
237 *     ..
238 *     .. External Subroutines ..
239       EXTERNAL           DLABAD, ZCOPY, ZLARFG, ZSCAL
240 *     ..
241 *     .. Statement Functions ..
242       DOUBLE PRECISION   CABS1
243 *     ..
244 *     .. Intrinsic Functions ..
245       INTRINSIC          ABS, DBLE, DCONJG, DIMAG, MAX, MIN, SQRT
246 *     ..
247 *     .. Statement Function definitions ..
248       CABS1( CDUM ) = ABS( DBLE( CDUM ) ) + ABS( DIMAG( CDUM ) )
249 *     ..
250 *     .. Executable Statements ..
251 *
252       INFO = 0
253 *
254 *     Quick return if possible
255 *
256       IF( N.EQ.0 )
257      $   RETURN
258       IF( ILO.EQ.IHI ) THEN
259          W( ILO ) = H( ILO, ILO )
260          RETURN
261       END IF
262 *
263 *     ==== clear out the trash ====
264       DO 10 J = ILO, IHI - 3
265          H( J+2, J ) = ZERO
266          H( J+3, J ) = ZERO
267    10 CONTINUE
268       IF( ILO.LE.IHI-2 )
269      $   H( IHI, IHI-2 ) = ZERO
270 *     ==== ensure that subdiagonal entries are real ====
271       IF( WANTT ) THEN
272          JLO = 1
273          JHI = N
274       ELSE
275          JLO = ILO
276          JHI = IHI
277       END IF
278       DO 20 I = ILO + 1, IHI
279          IF( DIMAG( H( I, I-1 ) ).NE.RZERO ) THEN
280 *           ==== The following redundant normalization
281 *           .    avoids problems with both gradual and
282 *           .    sudden underflow in ABS(H(I,I-1)) ====
283             SC = H( I, I-1 ) / CABS1( H( I, I-1 ) )
284             SC = DCONJG( SC ) / ABS( SC )
285             H( I, I-1 ) = ABS( H( I, I-1 ) )
286             CALL ZSCAL( JHI-I+1, SC, H( I, I ), LDH )
287             CALL ZSCAL( MIN( JHI, I+1 )-JLO+1, DCONJG( SC ),
288      $                  H( JLO, I ), 1 )
289             IF( WANTZ )
290      $         CALL ZSCAL( IHIZ-ILOZ+1, DCONJG( SC ), Z( ILOZ, I ), 1 )
291          END IF
292    20 CONTINUE
293 *
294       NH = IHI - ILO + 1
295       NZ = IHIZ - ILOZ + 1
296 *
297 *     Set machine-dependent constants for the stopping criterion.
298 *
299       SAFMIN = DLAMCH( 'SAFE MINIMUM' )
300       SAFMAX = RONE / SAFMIN
301       CALL DLABAD( SAFMIN, SAFMAX )
302       ULP = DLAMCH( 'PRECISION' )
303       SMLNUM = SAFMIN*( DBLE( NH ) / ULP )
304 *
305 *     I1 and I2 are the indices of the first row and last column of H
306 *     to which transformations must be applied. If eigenvalues only are
307 *     being computed, I1 and I2 are set inside the main loop.
308 *
309       IF( WANTT ) THEN
310          I1 = 1
311          I2 = N
312       END IF
313 *
314 *     ITMAX is the total number of QR iterations allowed.
315 *
316       ITMAX = 30 * MAX( 10, NH ) 
317 *
318 *     The main loop begins here. I is the loop index and decreases from
319 *     IHI to ILO in steps of 1. Each iteration of the loop works
320 *     with the active submatrix in rows and columns L to I.
321 *     Eigenvalues I+1 to IHI have already converged. Either L = ILO, or
322 *     H(L,L-1) is negligible so that the matrix splits.
323 *
324       I = IHI
325    30 CONTINUE
326       IF( I.LT.ILO )
327      $   GO TO 150
328 *
329 *     Perform QR iterations on rows and columns ILO to I until a
330 *     submatrix of order 1 splits off at the bottom because a
331 *     subdiagonal element has become negligible.
332 *
333       L = ILO
334       DO 130 ITS = 0, ITMAX
335 *
336 *        Look for a single small subdiagonal element.
337 *
338          DO 40 K = I, L + 1, -1
339             IF( CABS1( H( K, K-1 ) ).LE.SMLNUM )
340      $         GO TO 50
341             TST = CABS1( H( K-1, K-1 ) ) + CABS1( H( K, K ) )
342             IF( TST.EQ.ZERO ) THEN
343                IF( K-2.GE.ILO )
344      $            TST = TST + ABS( DBLE( H( K-1, K-2 ) ) )
345                IF( K+1.LE.IHI )
346      $            TST = TST + ABS( DBLE( H( K+1, K ) ) )
347             END IF
348 *           ==== The following is a conservative small subdiagonal
349 *           .    deflation criterion due to Ahues & Tisseur (LAWN 122,
350 *           .    1997). It has better mathematical foundation and
351 *           .    improves accuracy in some examples.  ====
352             IF( ABS( DBLE( H( K, K-1 ) ) ).LE.ULP*TST ) THEN
353                AB = MAX( CABS1( H( K, K-1 ) ), CABS1( H( K-1, K ) ) )
354                BA = MIN( CABS1( H( K, K-1 ) ), CABS1( H( K-1, K ) ) )
355                AA = MAX( CABS1( H( K, K ) ),
356      $              CABS1( H( K-1, K-1 )-H( K, K ) ) )
357                BB = MIN( CABS1( H( K, K ) ),
358      $              CABS1( H( K-1, K-1 )-H( K, K ) ) )
359                S = AA + AB
360                IF( BA*( AB / S ).LE.MAX( SMLNUM,
361      $             ULP*( BB*( AA / S ) ) ) )GO TO 50
362             END IF
363    40    CONTINUE
364    50    CONTINUE
365          L = K
366          IF( L.GT.ILO ) THEN
367 *
368 *           H(L,L-1) is negligible
369 *
370             H( L, L-1 ) = ZERO
371          END IF
372 *
373 *        Exit from loop if a submatrix of order 1 has split off.
374 *
375          IF( L.GE.I )
376      $      GO TO 140
377 *
378 *        Now the active submatrix is in rows and columns L to I. If
379 *        eigenvalues only are being computed, only the active submatrix
380 *        need be transformed.
381 *
382          IF( .NOT.WANTT ) THEN
383             I1 = L
384             I2 = I
385          END IF
386 *
387          IF( ITS.EQ.10 ) THEN
388 *
389 *           Exceptional shift.
390 *
391             S = DAT1*ABS( DBLE( H( L+1, L ) ) )
392             T = S + H( L, L )
393          ELSE IF( ITS.EQ.20 ) THEN
394 *
395 *           Exceptional shift.
396 *
397             S = DAT1*ABS( DBLE( H( I, I-1 ) ) )
398             T = S + H( I, I )
399          ELSE
400 *
401 *           Wilkinson's shift.
402 *
403             T = H( I, I )
404             U = SQRT( H( I-1, I ) )*SQRT( H( I, I-1 ) )
405             S = CABS1( U )
406             IF( S.NE.RZERO ) THEN
407                X = HALF*( H( I-1, I-1 )-T )
408                SX = CABS1( X )
409                S = MAX( S, CABS1( X ) )
410                Y = S*SQRT( ( X / S )**2+( U / S )**2 )
411                IF( SX.GT.RZERO ) THEN
412                   IF( DBLE( X / SX )*DBLE( Y )+DIMAG( X / SX )*
413      $                DIMAG( Y ).LT.RZERO )Y = -Y
414                END IF
415                T = T - U*ZLADIV( U, ( X+Y ) )
416             END IF
417          END IF
418 *
419 *        Look for two consecutive small subdiagonal elements.
420 *
421          DO 60 M = I - 1, L + 1, -1
422 *
423 *           Determine the effect of starting the single-shift QR
424 *           iteration at row M, and see if this would make H(M,M-1)
425 *           negligible.
426 *
427             H11 = H( M, M )
428             H22 = H( M+1, M+1 )
429             H11S = H11 - T
430             H21 = DBLE( H( M+1, M ) )
431             S = CABS1( H11S ) + ABS( H21 )
432             H11S = H11S / S
433             H21 = H21 / S
434             V( 1 ) = H11S
435             V( 2 ) = H21
436             H10 = DBLE( H( M, M-1 ) )
437             IF( ABS( H10 )*ABS( H21 ).LE.ULP*
438      $          ( CABS1( H11S )*( CABS1( H11 )+CABS1( H22 ) ) ) )
439      $          GO TO 70
440    60    CONTINUE
441          H11 = H( L, L )
442          H22 = H( L+1, L+1 )
443          H11S = H11 - T
444          H21 = DBLE( H( L+1, L ) )
445          S = CABS1( H11S ) + ABS( H21 )
446          H11S = H11S / S
447          H21 = H21 / S
448          V( 1 ) = H11S
449          V( 2 ) = H21
450    70    CONTINUE
451 *
452 *        Single-shift QR step
453 *
454          DO 120 K = M, I - 1
455 *
456 *           The first iteration of this loop determines a reflection G
457 *           from the vector V and applies it from left and right to H,
458 *           thus creating a nonzero bulge below the subdiagonal.
459 *
460 *           Each subsequent iteration determines a reflection G to
461 *           restore the Hessenberg form in the (K-1)th column, and thus
462 *           chases the bulge one step toward the bottom of the active
463 *           submatrix.
464 *
465 *           V(2) is always real before the call to ZLARFG, and hence
466 *           after the call T2 ( = T1*V(2) ) is also real.
467 *
468             IF( K.GT.M )
469      $         CALL ZCOPY( 2, H( K, K-1 ), 1, V, 1 )
470             CALL ZLARFG( 2, V( 1 ), V( 2 ), 1, T1 )
471             IF( K.GT.M ) THEN
472                H( K, K-1 ) = V( 1 )
473                H( K+1, K-1 ) = ZERO
474             END IF
475             V2 = V( 2 )
476             T2 = DBLE( T1*V2 )
477 *
478 *           Apply G from the left to transform the rows of the matrix
479 *           in columns K to I2.
480 *
481             DO 80 J = K, I2
482                SUM = DCONJG( T1 )*H( K, J ) + T2*H( K+1, J )
483                H( K, J ) = H( K, J ) - SUM
484                H( K+1, J ) = H( K+1, J ) - SUM*V2
485    80       CONTINUE
486 *
487 *           Apply G from the right to transform the columns of the
488 *           matrix in rows I1 to min(K+2,I).
489 *
490             DO 90 J = I1, MIN( K+2, I )
491                SUM = T1*H( J, K ) + T2*H( J, K+1 )
492                H( J, K ) = H( J, K ) - SUM
493                H( J, K+1 ) = H( J, K+1 ) - SUM*DCONJG( V2 )
494    90       CONTINUE
495 *
496             IF( WANTZ ) THEN
497 *
498 *              Accumulate transformations in the matrix Z
499 *
500                DO 100 J = ILOZ, IHIZ
501                   SUM = T1*Z( J, K ) + T2*Z( J, K+1 )
502                   Z( J, K ) = Z( J, K ) - SUM
503                   Z( J, K+1 ) = Z( J, K+1 ) - SUM*DCONJG( V2 )
504   100          CONTINUE
505             END IF
506 *
507             IF( K.EQ.M .AND. M.GT.L ) THEN
508 *
509 *              If the QR step was started at row M > L because two
510 *              consecutive small subdiagonals were found, then extra
511 *              scaling must be performed to ensure that H(M,M-1) remains
512 *              real.
513 *
514                TEMP = ONE - T1
515                TEMP = TEMP / ABS( TEMP )
516                H( M+1, M ) = H( M+1, M )*DCONJG( TEMP )
517                IF( M+2.LE.I )
518      $            H( M+2, M+1 ) = H( M+2, M+1 )*TEMP
519                DO 110 J = M, I
520                   IF( J.NE.M+1 ) THEN
521                      IF( I2.GT.J )
522      $                  CALL ZSCAL( I2-J, TEMP, H( J, J+1 ), LDH )
523                      CALL ZSCAL( J-I1, DCONJG( TEMP ), H( I1, J ), 1 )
524                      IF( WANTZ ) THEN
525                         CALL ZSCAL( NZ, DCONJG( TEMP ), Z( ILOZ, J ),
526      $                              1 )
527                      END IF
528                   END IF
529   110          CONTINUE
530             END IF
531   120    CONTINUE
532 *
533 *        Ensure that H(I,I-1) is real.
534 *
535          TEMP = H( I, I-1 )
536          IF( DIMAG( TEMP ).NE.RZERO ) THEN
537             RTEMP = ABS( TEMP )
538             H( I, I-1 ) = RTEMP
539             TEMP = TEMP / RTEMP
540             IF( I2.GT.I )
541      $         CALL ZSCAL( I2-I, DCONJG( TEMP ), H( I, I+1 ), LDH )
542             CALL ZSCAL( I-I1, TEMP, H( I1, I ), 1 )
543             IF( WANTZ ) THEN
544                CALL ZSCAL( NZ, TEMP, Z( ILOZ, I ), 1 )
545             END IF
546          END IF
547 *
548   130 CONTINUE
549 *
550 *     Failure to converge in remaining number of iterations
551 *
552       INFO = I
553       RETURN
554 *
555   140 CONTINUE
556 *
557 *     H(I,I-1) is negligible: one eigenvalue has converged.
558 *
559       W( I ) = H( I, I )
560 *
561 *     return to start of the main loop with new value of I.
562 *
563       I = L - 1
564       GO TO 30
565 *
566   150 CONTINUE
567       RETURN
568 *
569 *     End of ZLAHQR
570 *
571       END