ScaLAPACK 2.1  2.1 ScaLAPACK: Scalable Linear Algebra PACKage
clahqr2.f
Go to the documentation of this file.
1  SUBROUTINE clahqr2( 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 H( LDH, * ), W( * ), Z( LDZ, * )
15 * ..
16 *
17 * Purpose
18 * =======
19 *
20 * CLAHQR2 is an auxiliary routine called by CHSEQR to update the
21 * eigenvalues and Schur decomposition already computed by CHSEQR, by
22 * dealing with the Hessenberg submatrix in rows and columns ILO to IHI.
23 * This version of CLAHQR (not the standard LAPACK version) uses a
24 * double-shift algorithm (like LAPACK's SLAHQR).
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 * CLAHQR 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 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 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 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, CLAHQR 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 ZERO
99  parameter( zero = ( 0.0e+0, 0.0e+0 ) )
100  REAL RZERO, RONE
101  parameter( rzero = 0.0e+0, rone = 1.0e+0 )
102  REAL DAT1, DAT2
103  parameter( dat1 = 0.75e+0, dat2 = -0.4375e+0 )
104 * ..
105 * .. Local Scalars ..
106  INTEGER I, I1, I2, ITN, ITS, J, K, L, M, NH, NR, NZ
107  REAL CS, OVFL, S, SMLNUM, TST1, ULP, UNFL
108  COMPLEX 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  REAL RWORK( 1 )
114  COMPLEX V( 3 )
115 * ..
116 * .. External Functions ..
117  REAL SLAMCH, CLANHS
118  EXTERNAL slamch, clanhs
119 * ..
120 * .. External Subroutines ..
121  EXTERNAL slabad, ccopy, clanv2, clarfg, crot
122 * ..
123 * .. Intrinsic Functions ..
124  INTRINSIC abs, real, conjg, aimag, max, min
125 * ..
126 * .. Statement Functions ..
127  REAL CABS1
128 * ..
129 * .. Statement Function definitions ..
130  cabs1( cdum ) = abs( real( cdum ) ) + abs( aimag( 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 = slamch( 'Safe minimum' )
152  ovfl = rone / unfl
153  CALL slabad( unfl, ovfl )
154  ulp = slamch( '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 = clanhs( '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( REAL( H( I,I-1 ) ) ) + ABS( REAL( 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 ccopy( nr, h( k, k-1 ), 1, v, 1 )
289  CALL clarfg( 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 ) - conjg( 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 = conjg( t1 )*h( k, j ) +
311  \$ conjg( t2 )*h( k+1, j ) +
312  \$ conjg( 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*conjg( v2 )
325  h( j, k+2 ) = h( j, k+2 ) - sum*conjg( 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*conjg( v2 )
337  z( j, k+2 ) = z( j, k+2 ) - sum*conjg( 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 = conjg( t1 )*h( k, j ) +
347  \$ conjg( 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*conjg( 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*conjg( 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 ) - CONJG( 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 )*CONJG( 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 clanv2( 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 crot( i2-i, h( i-1, i+1 ), ldh, h( i, i+1 ), ldh,
418  \$ cs, sn )
419  CALL crot( i-i1-1, h( i1, i-1 ), 1, h( i1, i ), 1, cs,
420  \$ conjg( sn ) )
421  END IF
422  IF( wantz ) THEN
423 *
424 * Apply the transformation to Z.
425 *
426  CALL crot( nz, z( iloz, i-1 ), 1, z( iloz, i ), 1, cs,
427  \$ conjg( 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 CLAHQR2
443 *
444  END
max
#define max(A, B)
Definition: pcgemr.c:180
clanv2
subroutine clanv2(A, B, C, D, RT1, RT2, CS, SN)
Definition: clanv2.f:2
clahqr2
subroutine clahqr2(WANTT, WANTZ, N, ILO, IHI, H, LDH, W, ILOZ, IHIZ, Z, LDZ, INFO)
Definition: clahqr2.f:3
min
#define min(A, B)
Definition: pcgemr.c:181