SCALAPACK 2.2.2
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches
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
min
#define min(A, B)
Definition pcgemr.c:181
zlahqr2
subroutine zlahqr2(wantt, wantz, n, ilo, ihi, h, ldh, w, iloz, ihiz, z, ldz, info)
Definition zlahqr2.f:3
zlanv2
subroutine zlanv2(a, b, c, d, rt1, rt2, cs, sn)
Definition zlanv2.f:2