OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
zlahqr2.f File Reference

Go to the source code of this file.

Functions/Subroutines

subroutine zlahqr2 (wantt, wantz, n, ilo, ihi, h, ldh, w, iloz, ihiz, z, ldz, info)

Function/Subroutine Documentation

◆ zlahqr2()

subroutine zlahqr2 ( logical wantt,
logical wantz,
integer n,
integer ilo,
integer ihi,
complex*16, dimension( ldh, * ) h,
integer ldh,
complex*16, dimension( * ) w,
integer iloz,
integer ihiz,
complex*16, dimension( ldz, * ) z,
integer ldz,
integer info )

Definition at line 1 of file zlahqr2.f.

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.LT. IF( IILO )
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.EQ. IF( TST1RZERO )
193 $ TST1 = ZLANHS( '1', I-L+1, H( L, L ), LDH, RWORK )
194.LE. IF( CABS1( H( K, K-1 ) )MAX( ULP*TST1, SMLNUM ) )
195 $ GO TO 30
196 20 CONTINUE
197 30 CONTINUE
198 L = K
199.GT. IF( LILO ) 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.GE. IF( LI-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.NOT. IF( WANTT ) THEN
216 I1 = L
217 I2 = I
218 END IF
219*
220.EQ..OR..EQ. IF( ITS10 ITS20 ) 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.EQ. IF( ML )
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.LE. IF( CABS1( H10 )*( CABS1( V2 )+CABS1( V3 ) )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.GT. IF( KM )
288 $ CALL ZCOPY( NR, H( K, K-1 ), 1, V, 1 )
289 CALL ZLARFG( NR, V( 1 ), V( 2 ), 1, T1 )
290.GT. IF( KM ) THEN
291 H( K, K-1 ) = V( 1 )
292 H( K+1, K-1 ) = ZERO
293.LT. IF( KI-1 )
294 $ H( K+2, K-1 ) = ZERO
295.GT. ELSE IF( ML ) 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.EQ. IF( NR3 ) 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.EQ. ELSE IF( NR2 ) 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.EQ. IF( LI ) THEN
397*
398* H(I,I-1) is negligible: one eigenvalue has converged.
399*
400 W( I ) = H( I, I )
401*
402.EQ. ELSE IF( LI-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.GT. IF( I2I )
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*
dlabad
subroutine dlabad(small, large)
DLABAD
Definition dlabad.f:74
zlanhs
double precision function zlanhs(norm, n, a, lda, work)
ZLANHS returns the value of the 1-norm, Frobenius norm, infinity-norm, or the largest absolute value ...
Definition zlanhs.f:109
zrot
subroutine zrot(n, cx, incx, cy, incy, c, s)
ZROT applies a plane rotation with real cosine and complex sine to a pair of complex vectors.
Definition zrot.f:103
zlarfg
subroutine zlarfg(n, alpha, x, incx, tau)
ZLARFG generates an elementary reflector (Householder matrix).
Definition zlarfg.f:106
zcopy
subroutine zcopy(n, zx, incx, zy, incy)
ZCOPY
Definition zcopy.f:81
dlamch
double precision function dlamch(cmach)
DLAMCH
Definition dlamch.f:69
min
#define min(a, b)
Definition macros.h:20
max
#define max(a, b)
Definition macros.h:21
zlanv2
subroutine zlanv2(a, b, c, d, rt1, rt2, cs, sn)
Definition zlanv2.f:2