ScaLAPACK 2.1  2.1
ScaLAPACK: Scalable Linear Algebra PACKage
zlahqr2.f
Go to the documentation of this file.
1  SUBROUTINE zlahqr2( WANTT, WANTZ, N, ILO, IHI, H, LDH, W, ILOZ,
2  $ IHIZ, Z, LDZ, INFO )
3 *
4 * -- ScaLAPACK routine (version 1.7) --
5 * Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
6 * Courant Institute, Argonne National Lab, and Rice University
7 * June 22, 2000
8 *
9 * .. Scalar Arguments ..
10  LOGICAL WANTT, WANTZ
11  INTEGER IHI, IHIZ, ILO, ILOZ, INFO, LDH, LDZ, N
12 * ..
13 * .. Array Arguments ..
14  COMPLEX*16 H( LDH, * ), W( * ), Z( LDZ, * )
15 * ..
16 *
17 * Purpose
18 * =======
19 *
20 * ZLAHQR2 is an auxiliary routine called by ZHSEQR to update the
21 * eigenvalues and Schur decomposition already computed by ZHSEQR, by
22 * dealing with the Hessenberg submatrix in rows and columns ILO to IHI.
23 * This version of ZLAHQR (not the standard LAPACK version) uses a
24 * double-shift algorithm (like LAPACK's DLAHQR).
25 * Unlike the standard LAPACK convention, this does not assume the
26 * subdiagonal is real, nor does it work to preserve this quality if
27 * given.
28 *
29 * Arguments
30 * =========
31 *
32 * WANTT (input) LOGICAL
33 * = .TRUE. : the full Schur form T is required;
34 * = .FALSE.: only eigenvalues are required.
35 *
36 * WANTZ (input) LOGICAL
37 * = .TRUE. : the matrix of Schur vectors Z is required;
38 * = .FALSE.: Schur vectors are not required.
39 *
40 * N (input) INTEGER
41 * The order of the matrix H. N >= 0.
42 *
43 * ILO (input) INTEGER
44 * IHI (input) INTEGER
45 * It is assumed that H is already upper triangular in rows and
46 * columns IHI+1:N, and that H(ILO,ILO-1) = 0 (unless ILO = 1).
47 * ZLAHQR works primarily with the Hessenberg submatrix in rows
48 * and columns ILO to IHI, but applies transformations to all of
49 * H if WANTT is .TRUE..
50 * 1 <= ILO <= max(1,IHI); IHI <= N.
51 *
52 * H (input/output) COMPLEX*16 array, dimension (LDH,N)
53 * On entry, the upper Hessenberg matrix H.
54 * On exit, if WANTT is .TRUE., H is upper triangular in rows
55 * and columns ILO:IHI. If WANTT is .FALSE., the contents of H
56 * are unspecified on exit.
57 *
58 * LDH (input) INTEGER
59 * The leading dimension of the array H. LDH >= max(1,N).
60 *
61 * W (output) COMPLEX*16 array, dimension (N)
62 * The computed eigenvalues ILO to IHI are stored in the
63 * corresponding elements of W. If WANTT is .TRUE., the
64 * eigenvalues are stored in the same order as on the diagonal
65 * of the Schur form returned in H, with W(i) = H(i,i).
66 *
67 * ILOZ (input) INTEGER
68 * IHIZ (input) INTEGER
69 * Specify the rows of Z to which transformations must be
70 * applied if WANTZ is .TRUE..
71 * 1 <= ILOZ <= ILO; IHI <= IHIZ <= N.
72 *
73 * Z (input/output) COMPLEX*16 array, dimension (LDZ,N)
74 * If WANTZ is .TRUE., on entry Z must contain the current
75 * matrix Z of transformations, and on exit Z has been updated;
76 * transformations are applied only to the submatrix
77 * Z(ILOZ:IHIZ,ILO:IHI). If WANTZ is .FALSE., Z is not
78 * referenced.
79 *
80 * LDZ (input) INTEGER
81 * The leading dimension of the array Z. LDZ >= max(1,N).
82 *
83 * INFO (output) INTEGER
84 * = 0: successful exit
85 * > 0: if INFO = i, ZLAHQR failed to compute all the
86 * eigenvalues ILO to IHI in a total of 30*(IHI-ILO+1)
87 * iterations; elements i+1:ihi of W contain those
88 * eigenvalues which have been successfully computed.
89 *
90 * Further Details
91 * ===============
92 *
93 * Modified by Mark R. Fahey, June, 2000
94 *
95 * =====================================================================
96 *
97 * .. Parameters ..
98  COMPLEX*16 ZERO
99  parameter( zero = ( 0.0d+0, 0.0d+0 ) )
100  DOUBLE PRECISION RZERO, RONE
101  parameter( rzero = 0.0d+0, rone = 1.0d+0 )
102  DOUBLE PRECISION DAT1, DAT2
103  parameter( dat1 = 0.75d+0, dat2 = -0.4375d+0 )
104 * ..
105 * .. Local Scalars ..
106  INTEGER I, I1, I2, ITN, ITS, J, K, L, M, NH, NR, NZ
107  DOUBLE PRECISION CS, OVFL, S, SMLNUM, TST1, ULP, UNFL
108  COMPLEX*16 CDUM, H00, H10, H11, H12, H21, H22, H33, H33S,
109  $ h43h34, h44, h44s, sn, sum, t1, t2, t3, v1, v2,
110  $ v3
111 * ..
112 * .. Local Arrays ..
113  DOUBLE PRECISION RWORK( 1 )
114  COMPLEX*16 V( 3 )
115 * ..
116 * .. External Functions ..
117  DOUBLE PRECISION DLAMCH, ZLANHS
118  EXTERNAL dlamch, zlanhs
119 * ..
120 * .. External Subroutines ..
121  EXTERNAL dlabad, zcopy, zlanv2, zlarfg, zrot
122 * ..
123 * .. Intrinsic Functions ..
124  INTRINSIC abs, dble, dconjg, dimag, max, min
125 * ..
126 * .. Statement Functions ..
127  DOUBLE PRECISION CABS1
128 * ..
129 * .. Statement Function definitions ..
130  cabs1( cdum ) = abs( dble( cdum ) ) + abs( dimag( cdum ) )
131 * ..
132 * .. Executable Statements ..
133 *
134  info = 0
135 *
136 * Quick return if possible
137 *
138  IF( n.EQ.0 )
139  $ RETURN
140  IF( ilo.EQ.ihi ) THEN
141  w( ilo ) = h( ilo, ilo )
142  RETURN
143  END IF
144 *
145  nh = ihi - ilo + 1
146  nz = ihiz - iloz + 1
147 *
148 * Set machine-dependent constants for the stopping criterion.
149 * If norm(H) <= sqrt(OVFL), overflow should not occur.
150 *
151  unfl = dlamch( 'Safe minimum' )
152  ovfl = rone / unfl
153  CALL dlabad( unfl, ovfl )
154  ulp = dlamch( 'Precision' )
155  smlnum = unfl*( nh / ulp )
156 *
157 * I1 and I2 are the indices of the first row and last column of H
158 * to which transformations must be applied. If eigenvalues only are
159 * being computed, I1 and I2 are set inside the main loop.
160 *
161  IF( wantt ) THEN
162  i1 = 1
163  i2 = n
164  END IF
165 *
166 * ITN is the total number of QR iterations allowed.
167 *
168  itn = 30*nh
169 *
170 * The main loop begins here. I is the loop index and decreases from
171 * IHI to ILO in steps of 1 or 2. Each iteration of the loop works
172 * with the active submatrix in rows and columns L to I.
173 * Eigenvalues I+1 to IHI have already converged. Either L = ILO, or
174 * H(L,L-1) is negligible so that the matrix splits.
175 *
176  i = ihi
177  10 CONTINUE
178  l = ilo
179  IF( i.LT.ilo )
180  $ GO TO 150
181 *
182 * Perform QR iterations on rows and columns ILO to I until a
183 * submatrix of order 1 or 2 splits off at the bottom because a
184 * subdiagonal element has become negligible.
185 *
186  DO 130 its = 0, itn
187 *
188 * Look for a single small subdiagonal element.
189 *
190  DO 20 k = i, l + 1, -1
191  tst1 = cabs1( h( k-1, k-1 ) ) + cabs1( h( k, k ) )
192  IF( tst1.EQ.rzero )
193  $ tst1 = zlanhs( '1', i-l+1, h( l, l ), ldh, rwork )
194  IF( cabs1( h( k, k-1 ) ).LE.max( ulp*tst1, smlnum ) )
195  $ GO TO 30
196  20 CONTINUE
197  30 CONTINUE
198  l = k
199  IF( l.GT.ilo ) THEN
200 *
201 * H(L,L-1) is negligible
202 *
203  h( l, l-1 ) = zero
204  END IF
205 *
206 * Exit from loop if a submatrix of order 1 or 2 has split off.
207 *
208  IF( l.GE.i-1 )
209  $ GO TO 140
210 *
211 * Now the active submatrix is in rows and columns L to I. If
212 * eigenvalues only are being computed, only the active submatrix
213 * need be transformed.
214 *
215  IF( .NOT.wantt ) THEN
216  i1 = l
217  i2 = i
218  END IF
219 *
220  IF( its.EQ.10 .OR. its.EQ.20 ) THEN
221 *
222 * Exceptional shift.
223 *
224 * S = ABS( DBLE( H( I,I-1 ) ) ) + ABS( DBLE( H( I-1,I-2 ) ) )
225  s = cabs1( h( i, i-1 ) ) + cabs1( h( i-1, i-2 ) )
226  h44 = dat1*s
227  h33 = h44
228  h43h34 = dat2*s*s
229  ELSE
230 *
231 * Prepare to use Wilkinson's shift.
232 *
233  h44 = h( i, i )
234  h33 = h( i-1, i-1 )
235  h43h34 = h( i, i-1 )*h( i-1, i )
236  END IF
237 *
238 * Look for two consecutive small subdiagonal elements.
239 *
240  DO 40 m = i - 2, l, -1
241 *
242 * Determine the effect of starting the double-shift QR
243 * iteration at row M, and see if this would make H(M,M-1)
244 * negligible.
245 *
246  h11 = h( m, m )
247  h22 = h( m+1, m+1 )
248  h21 = h( m+1, m )
249  h12 = h( m, m+1 )
250  h44s = h44 - h11
251  h33s = h33 - h11
252  v1 = ( h33s*h44s-h43h34 ) / h21 + h12
253  v2 = h22 - h11 - h33s - h44s
254  v3 = h( m+2, m+1 )
255  s = cabs1( v1 ) + cabs1( v2 ) + abs( v3 )
256  v1 = v1 / s
257  v2 = v2 / s
258  v3 = v3 / s
259  v( 1 ) = v1
260  v( 2 ) = v2
261  v( 3 ) = v3
262  IF( m.EQ.l )
263  $ GO TO 50
264  h00 = h( m-1, m-1 )
265  h10 = h( m, m-1 )
266  tst1 = cabs1( v1 )*( cabs1( h00 )+cabs1( h11 )+
267  $ cabs1( h22 ) )
268  IF( cabs1( h10 )*( cabs1( v2 )+cabs1( v3 ) ).LE.ulp*tst1 )
269  $ GO TO 50
270  40 CONTINUE
271  50 CONTINUE
272 *
273 * Double-shift QR step
274 *
275  DO 120 k = m, i - 1
276 *
277 * The first iteration of this loop determines a reflection G
278 * from the vector V and applies it from left and right to H,
279 * thus creating a nonzero bulge below the subdiagonal.
280 *
281 * Each subsequent iteration determines a reflection G to
282 * restore the Hessenberg form in the (K-1)th column, and thus
283 * chases the bulge one step toward the bottom of the active
284 * submatrix. NR is the order of G
285 *
286  nr = min( 3, i-k+1 )
287  IF( k.GT.m )
288  $ CALL zcopy( nr, h( k, k-1 ), 1, v, 1 )
289  CALL zlarfg( nr, v( 1 ), v( 2 ), 1, t1 )
290  IF( k.GT.m ) THEN
291  h( k, k-1 ) = v( 1 )
292  h( k+1, k-1 ) = zero
293  IF( k.LT.i-1 )
294  $ h( k+2, k-1 ) = zero
295  ELSE IF( m.GT.l ) THEN
296 * The real double-shift code uses H( K, K-1 ) = -H( K, K-1 )
297 * instead of the following.
298  h( k, k-1 ) = h( k, k-1 ) - dconjg( t1 )*h( k, k-1 )
299  END IF
300  v2 = v( 2 )
301  t2 = t1*v2
302  IF( nr.EQ.3 ) THEN
303  v3 = v( 3 )
304  t3 = t1*v3
305 *
306 * Apply G from the left to transform the rows of the matrix
307 * in columns K to I2.
308 *
309  DO 60 j = k, i2
310  sum = dconjg( t1 )*h( k, j ) +
311  $ dconjg( t2 )*h( k+1, j ) +
312  $ dconjg( t3 )*h( k+2, j )
313  h( k, j ) = h( k, j ) - sum
314  h( k+1, j ) = h( k+1, j ) - sum*v2
315  h( k+2, j ) = h( k+2, j ) - sum*v3
316  60 CONTINUE
317 *
318 * Apply G from the right to transform the columns of the
319 * matrix in rows I1 to min(K+3,I).
320 *
321  DO 70 j = i1, min( k+3, i )
322  sum = t1*h( j, k ) + t2*h( j, k+1 ) + t3*h( j, k+2 )
323  h( j, k ) = h( j, k ) - sum
324  h( j, k+1 ) = h( j, k+1 ) - sum*dconjg( v2 )
325  h( j, k+2 ) = h( j, k+2 ) - sum*dconjg( v3 )
326  70 CONTINUE
327 *
328  IF( wantz ) THEN
329 *
330 * Accumulate transformations in the matrix Z
331 *
332  DO 80 j = iloz, ihiz
333  sum = t1*z( j, k ) + t2*z( j, k+1 ) +
334  $ t3*z( j, k+2 )
335  z( j, k ) = z( j, k ) - sum
336  z( j, k+1 ) = z( j, k+1 ) - sum*dconjg( v2 )
337  z( j, k+2 ) = z( j, k+2 ) - sum*dconjg( v3 )
338  80 CONTINUE
339  END IF
340  ELSE IF( nr.EQ.2 ) THEN
341 *
342 * Apply G from the left to transform the rows of the matrix
343 * in columns K to I2.
344 *
345  DO 90 j = k, i2
346  sum = dconjg( t1 )*h( k, j ) +
347  $ dconjg( t2 )*h( k+1, j )
348  h( k, j ) = h( k, j ) - sum
349  h( k+1, j ) = h( k+1, j ) - sum*v2
350  90 CONTINUE
351 *
352 * Apply G from the right to transform the columns of the
353 * matrix in rows I1 to min(K+2,I).
354 *
355  DO 100 j = i1, min( k+2, i )
356  sum = t1*h( j, k ) + t2*h( j, k+1 )
357  h( j, k ) = h( j, k ) - sum
358  h( j, k+1 ) = h( j, k+1 ) - sum*dconjg( v2 )
359  100 CONTINUE
360 *
361  IF( wantz ) THEN
362 *
363 * Accumulate transformations in the matrix Z
364 *
365  DO 110 j = iloz, ihiz
366  sum = t1*z( j, k ) + t2*z( j, k+1 )
367  z( j, k ) = z( j, k ) - sum
368  z( j, k+1 ) = z( j, k+1 ) - sum*dconjg( v2 )
369  110 CONTINUE
370  END IF
371  END IF
372 *
373 * Since at the start of the QR step we have for M > L
374 * H( K, K-1 ) = H( K, K-1 ) - DCONJG( T1 )*H( K, K-1 )
375 * then we don't need to do the following
376 * IF( K.EQ.M .AND. M.GT.L ) THEN
377 * If the QR step was started at row M > L because two
378 * consecutive small subdiagonals were found, then H(M,M-1)
379 * must also be updated by a factor of (1-T1).
380 * TEMP = ONE - T1
381 * H( m, m-1 ) = H( m, m-1 )*DCONJG( TEMP )
382 * END IF
383  120 CONTINUE
384 *
385 * Ensure that H(I,I-1) is real.
386 *
387  130 CONTINUE
388 *
389 * Failure to converge in remaining number of iterations
390 *
391  info = i
392  RETURN
393 *
394  140 CONTINUE
395 *
396  IF( l.EQ.i ) THEN
397 *
398 * H(I,I-1) is negligible: one eigenvalue has converged.
399 *
400  w( i ) = h( i, i )
401 *
402  ELSE IF( l.EQ.i-1 ) THEN
403 *
404 * H(I-1,I-2) is negligible: a pair of eigenvalues have converged.
405 *
406 * Transform the 2-by-2 submatrix to standard Schur form,
407 * and compute and store the eigenvalues.
408 *
409  CALL zlanv2( h( i-1, i-1 ), h( i-1, i ), h( i, i-1 ),
410  $ h( i, i ), w( i-1 ), w( i ), cs, sn )
411 *
412  IF( wantt ) THEN
413 *
414 * Apply the transformation to the rest of H.
415 *
416  IF( i2.GT.i )
417  $ CALL zrot( i2-i, h( i-1, i+1 ), ldh, h( i, i+1 ), ldh,
418  $ cs, sn )
419  CALL zrot( i-i1-1, h( i1, i-1 ), 1, h( i1, i ), 1, cs,
420  $ dconjg( sn ) )
421  END IF
422  IF( wantz ) THEN
423 *
424 * Apply the transformation to Z.
425 *
426  CALL zrot( nz, z( iloz, i-1 ), 1, z( iloz, i ), 1, cs,
427  $ dconjg( sn ) )
428  END IF
429 *
430  END IF
431 *
432 * Decrement number of remaining iterations, and return to start of
433 * the main loop with new value of I.
434 *
435  itn = itn - its
436  i = l - 1
437  GO TO 10
438 *
439  150 CONTINUE
440  RETURN
441 *
442 * End of ZLAHQR2
443 *
444  END
max
#define max(A, B)
Definition: pcgemr.c:180
zlanv2
subroutine zlanv2(A, B, C, D, RT1, RT2, CS, SN)
Definition: zlanv2.f:2
zlahqr2
subroutine zlahqr2(WANTT, WANTZ, N, ILO, IHI, H, LDH, W, ILOZ, IHIZ, Z, LDZ, INFO)
Definition: zlahqr2.f:3
min
#define min(A, B)
Definition: pcgemr.c:181