LAPACK 3.12.1
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches
ztrevc.f
Go to the documentation of this file.
1*> \brief \b ZTREVC
2*
3* =========== DOCUMENTATION ===========
4*
5* Online html documentation available at
6* http://www.netlib.org/lapack/explore-html/
7*
8*> Download ZTREVC + dependencies
9*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/ztrevc.f">
10*> [TGZ]</a>
11*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/ztrevc.f">
12*> [ZIP]</a>
13*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/ztrevc.f">
14*> [TXT]</a>
15*
16* Definition:
17* ===========
18*
19* SUBROUTINE ZTREVC( SIDE, HOWMNY, SELECT, N, T, LDT, VL, LDVL, VR,
20* LDVR, MM, M, WORK, RWORK, INFO )
21*
22* .. Scalar Arguments ..
23* CHARACTER HOWMNY, SIDE
24* INTEGER INFO, LDT, LDVL, LDVR, M, MM, N
25* ..
26* .. Array Arguments ..
27* LOGICAL SELECT( * )
28* DOUBLE PRECISION RWORK( * )
29* COMPLEX*16 T( LDT, * ), VL( LDVL, * ), VR( LDVR, * ),
30* $ WORK( * )
31* ..
32*
33*
34*> \par Purpose:
35* =============
36*>
37*> \verbatim
38*>
39*> ZTREVC computes some or all of the right and/or left eigenvectors of
40*> a complex upper triangular matrix T.
41*> Matrices of this type are produced by the Schur factorization of
42*> a complex general matrix: A = Q*T*Q**H, as computed by ZHSEQR.
43*>
44*> The right eigenvector x and the left eigenvector y of T corresponding
45*> to an eigenvalue w are defined by:
46*>
47*> T*x = w*x, (y**H)*T = w*(y**H)
48*>
49*> where y**H denotes the conjugate transpose of the vector y.
50*> The eigenvalues are not input to this routine, but are read directly
51*> from the diagonal of T.
52*>
53*> This routine returns the matrices X and/or Y of right and left
54*> eigenvectors of T, or the products Q*X and/or Q*Y, where Q is an
55*> input matrix. If Q is the unitary factor that reduces a matrix A to
56*> Schur form T, then Q*X and Q*Y are the matrices of right and left
57*> eigenvectors of A.
58*> \endverbatim
59*
60* Arguments:
61* ==========
62*
63*> \param[in] SIDE
64*> \verbatim
65*> SIDE is CHARACTER*1
66*> = 'R': compute right eigenvectors only;
67*> = 'L': compute left eigenvectors only;
68*> = 'B': compute both right and left eigenvectors.
69*> \endverbatim
70*>
71*> \param[in] HOWMNY
72*> \verbatim
73*> HOWMNY is CHARACTER*1
74*> = 'A': compute all right and/or left eigenvectors;
75*> = 'B': compute all right and/or left eigenvectors,
76*> backtransformed using the matrices supplied in
77*> VR and/or VL;
78*> = 'S': compute selected right and/or left eigenvectors,
79*> as indicated by the logical array SELECT.
80*> \endverbatim
81*>
82*> \param[in] SELECT
83*> \verbatim
84*> SELECT is LOGICAL array, dimension (N)
85*> If HOWMNY = 'S', SELECT specifies the eigenvectors to be
86*> computed.
87*> The eigenvector corresponding to the j-th eigenvalue is
88*> computed if SELECT(j) = .TRUE..
89*> Not referenced if HOWMNY = 'A' or 'B'.
90*> \endverbatim
91*>
92*> \param[in] N
93*> \verbatim
94*> N is INTEGER
95*> The order of the matrix T. N >= 0.
96*> \endverbatim
97*>
98*> \param[in,out] T
99*> \verbatim
100*> T is COMPLEX*16 array, dimension (LDT,N)
101*> The upper triangular matrix T. T is modified, but restored
102*> on exit.
103*> \endverbatim
104*>
105*> \param[in] LDT
106*> \verbatim
107*> LDT is INTEGER
108*> The leading dimension of the array T. LDT >= max(1,N).
109*> \endverbatim
110*>
111*> \param[in,out] VL
112*> \verbatim
113*> VL is COMPLEX*16 array, dimension (LDVL,MM)
114*> On entry, if SIDE = 'L' or 'B' and HOWMNY = 'B', VL must
115*> contain an N-by-N matrix Q (usually the unitary matrix Q of
116*> Schur vectors returned by ZHSEQR).
117*> On exit, if SIDE = 'L' or 'B', VL contains:
118*> if HOWMNY = 'A', the matrix Y of left eigenvectors of T;
119*> if HOWMNY = 'B', the matrix Q*Y;
120*> if HOWMNY = 'S', the left eigenvectors of T specified by
121*> SELECT, stored consecutively in the columns
122*> of VL, in the same order as their
123*> eigenvalues.
124*> Not referenced if SIDE = 'R'.
125*> \endverbatim
126*>
127*> \param[in] LDVL
128*> \verbatim
129*> LDVL is INTEGER
130*> The leading dimension of the array VL. LDVL >= 1, and if
131*> SIDE = 'L' or 'B', LDVL >= N.
132*> \endverbatim
133*>
134*> \param[in,out] VR
135*> \verbatim
136*> VR is COMPLEX*16 array, dimension (LDVR,MM)
137*> On entry, if SIDE = 'R' or 'B' and HOWMNY = 'B', VR must
138*> contain an N-by-N matrix Q (usually the unitary matrix Q of
139*> Schur vectors returned by ZHSEQR).
140*> On exit, if SIDE = 'R' or 'B', VR contains:
141*> if HOWMNY = 'A', the matrix X of right eigenvectors of T;
142*> if HOWMNY = 'B', the matrix Q*X;
143*> if HOWMNY = 'S', the right eigenvectors of T specified by
144*> SELECT, stored consecutively in the columns
145*> of VR, in the same order as their
146*> eigenvalues.
147*> Not referenced if SIDE = 'L'.
148*> \endverbatim
149*>
150*> \param[in] LDVR
151*> \verbatim
152*> LDVR is INTEGER
153*> The leading dimension of the array VR. LDVR >= 1, and if
154*> SIDE = 'R' or 'B'; LDVR >= N.
155*> \endverbatim
156*>
157*> \param[in] MM
158*> \verbatim
159*> MM is INTEGER
160*> The number of columns in the arrays VL and/or VR. MM >= M.
161*> \endverbatim
162*>
163*> \param[out] M
164*> \verbatim
165*> M is INTEGER
166*> The number of columns in the arrays VL and/or VR actually
167*> used to store the eigenvectors. If HOWMNY = 'A' or 'B', M
168*> is set to N. Each selected eigenvector occupies one
169*> column.
170*> \endverbatim
171*>
172*> \param[out] WORK
173*> \verbatim
174*> WORK is COMPLEX*16 array, dimension (2*N)
175*> \endverbatim
176*>
177*> \param[out] RWORK
178*> \verbatim
179*> RWORK is DOUBLE PRECISION array, dimension (N)
180*> \endverbatim
181*>
182*> \param[out] INFO
183*> \verbatim
184*> INFO is INTEGER
185*> = 0: successful exit
186*> < 0: if INFO = -i, the i-th argument had an illegal value
187*> \endverbatim
188*
189* Authors:
190* ========
191*
192*> \author Univ. of Tennessee
193*> \author Univ. of California Berkeley
194*> \author Univ. of Colorado Denver
195*> \author NAG Ltd.
196*
197*> \ingroup trevc
198*
199*> \par Further Details:
200* =====================
201*>
202*> \verbatim
203*>
204*> The algorithm used in this program is basically backward (forward)
205*> substitution, with scaling to make the the code robust against
206*> possible overflow.
207*>
208*> Each eigenvector is normalized so that the element of largest
209*> magnitude has magnitude 1; here the magnitude of a complex number
210*> (x,y) is taken to be |x| + |y|.
211*> \endverbatim
212*>
213* =====================================================================
214 SUBROUTINE ztrevc( SIDE, HOWMNY, SELECT, N, T, LDT, VL, LDVL,
215 $ VR,
216 $ LDVR, MM, M, WORK, RWORK, INFO )
217*
218* -- LAPACK computational routine --
219* -- LAPACK is a software package provided by Univ. of Tennessee, --
220* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
221*
222* .. Scalar Arguments ..
223 CHARACTER HOWMNY, SIDE
224 INTEGER INFO, LDT, LDVL, LDVR, M, MM, N
225* ..
226* .. Array Arguments ..
227 LOGICAL SELECT( * )
228 DOUBLE PRECISION RWORK( * )
229 COMPLEX*16 T( LDT, * ), VL( LDVL, * ), VR( LDVR, * ),
230 $ work( * )
231* ..
232*
233* =====================================================================
234*
235* .. Parameters ..
236 DOUBLE PRECISION ZERO, ONE
237 PARAMETER ( ZERO = 0.0d+0, one = 1.0d+0 )
238 COMPLEX*16 CMZERO, CMONE
239 parameter( cmzero = ( 0.0d+0, 0.0d+0 ),
240 $ cmone = ( 1.0d+0, 0.0d+0 ) )
241* ..
242* .. Local Scalars ..
243 LOGICAL ALLV, BOTHV, LEFTV, OVER, RIGHTV, SOMEV
244 INTEGER I, II, IS, J, K, KI
245 DOUBLE PRECISION OVFL, REMAX, SCALE, SMIN, SMLNUM, ULP, UNFL
246 COMPLEX*16 CDUM
247* ..
248* .. External Functions ..
249 LOGICAL LSAME
250 INTEGER IZAMAX
251 DOUBLE PRECISION DLAMCH, DZASUM
252 EXTERNAL lsame, izamax, dlamch, dzasum
253* ..
254* .. External Subroutines ..
255 EXTERNAL xerbla, zcopy, zdscal, zgemv,
256 $ zlatrs
257* ..
258* .. Intrinsic Functions ..
259 INTRINSIC abs, dble, dcmplx, dconjg, dimag, max
260* ..
261* .. Statement Functions ..
262 DOUBLE PRECISION CABS1
263* ..
264* .. Statement Function definitions ..
265 cabs1( cdum ) = abs( dble( cdum ) ) + abs( dimag( cdum ) )
266* ..
267* .. Executable Statements ..
268*
269* Decode and test the input parameters
270*
271 bothv = lsame( side, 'B' )
272 rightv = lsame( side, 'R' ) .OR. bothv
273 leftv = lsame( side, 'L' ) .OR. bothv
274*
275 allv = lsame( howmny, 'A' )
276 over = lsame( howmny, 'B' )
277 somev = lsame( howmny, 'S' )
278*
279* Set M to the number of columns required to store the selected
280* eigenvectors.
281*
282 IF( somev ) THEN
283 m = 0
284 DO 10 j = 1, n
285 IF( SELECT( j ) )
286 $ m = m + 1
287 10 CONTINUE
288 ELSE
289 m = n
290 END IF
291*
292 info = 0
293 IF( .NOT.rightv .AND. .NOT.leftv ) THEN
294 info = -1
295 ELSE IF( .NOT.allv .AND. .NOT.over .AND. .NOT.somev ) THEN
296 info = -2
297 ELSE IF( n.LT.0 ) THEN
298 info = -4
299 ELSE IF( ldt.LT.max( 1, n ) ) THEN
300 info = -6
301 ELSE IF( ldvl.LT.1 .OR. ( leftv .AND. ldvl.LT.n ) ) THEN
302 info = -8
303 ELSE IF( ldvr.LT.1 .OR. ( rightv .AND. ldvr.LT.n ) ) THEN
304 info = -10
305 ELSE IF( mm.LT.m ) THEN
306 info = -11
307 END IF
308 IF( info.NE.0 ) THEN
309 CALL xerbla( 'ZTREVC', -info )
310 RETURN
311 END IF
312*
313* Quick return if possible.
314*
315 IF( n.EQ.0 )
316 $ RETURN
317*
318* Set the constants to control overflow.
319*
320 unfl = dlamch( 'Safe minimum' )
321 ovfl = one / unfl
322 ulp = dlamch( 'Precision' )
323 smlnum = unfl*( n / ulp )
324*
325* Store the diagonal elements of T in working array WORK.
326*
327 DO 20 i = 1, n
328 work( i+n ) = t( i, i )
329 20 CONTINUE
330*
331* Compute 1-norm of each column of strictly upper triangular
332* part of T to control overflow in triangular solver.
333*
334 rwork( 1 ) = zero
335 DO 30 j = 2, n
336 rwork( j ) = dzasum( j-1, t( 1, j ), 1 )
337 30 CONTINUE
338*
339 IF( rightv ) THEN
340*
341* Compute right eigenvectors.
342*
343 is = m
344 DO 80 ki = n, 1, -1
345*
346 IF( somev ) THEN
347 IF( .NOT.SELECT( ki ) )
348 $ GO TO 80
349 END IF
350 smin = max( ulp*( cabs1( t( ki, ki ) ) ), smlnum )
351*
352 work( 1 ) = cmone
353*
354* Form right-hand side.
355*
356 DO 40 k = 1, ki - 1
357 work( k ) = -t( k, ki )
358 40 CONTINUE
359*
360* Solve the triangular system:
361* (T(1:KI-1,1:KI-1) - T(KI,KI))*X = SCALE*WORK.
362*
363 DO 50 k = 1, ki - 1
364 t( k, k ) = t( k, k ) - t( ki, ki )
365 IF( cabs1( t( k, k ) ).LT.smin )
366 $ t( k, k ) = smin
367 50 CONTINUE
368*
369 IF( ki.GT.1 ) THEN
370 CALL zlatrs( 'Upper', 'No transpose', 'Non-unit', 'Y',
371 $ ki-1, t, ldt, work( 1 ), scale, rwork,
372 $ info )
373 work( ki ) = scale
374 END IF
375*
376* Copy the vector x or Q*x to VR and normalize.
377*
378 IF( .NOT.over ) THEN
379 CALL zcopy( ki, work( 1 ), 1, vr( 1, is ), 1 )
380*
381 ii = izamax( ki, vr( 1, is ), 1 )
382 remax = one / cabs1( vr( ii, is ) )
383 CALL zdscal( ki, remax, vr( 1, is ), 1 )
384*
385 DO 60 k = ki + 1, n
386 vr( k, is ) = cmzero
387 60 CONTINUE
388 ELSE
389 IF( ki.GT.1 )
390 $ CALL zgemv( 'N', n, ki-1, cmone, vr, ldvr,
391 $ work( 1 ),
392 $ 1, dcmplx( scale ), vr( 1, ki ), 1 )
393*
394 ii = izamax( n, vr( 1, ki ), 1 )
395 remax = one / cabs1( vr( ii, ki ) )
396 CALL zdscal( n, remax, vr( 1, ki ), 1 )
397 END IF
398*
399* Set back the original diagonal elements of T.
400*
401 DO 70 k = 1, ki - 1
402 t( k, k ) = work( k+n )
403 70 CONTINUE
404*
405 is = is - 1
406 80 CONTINUE
407 END IF
408*
409 IF( leftv ) THEN
410*
411* Compute left eigenvectors.
412*
413 is = 1
414 DO 130 ki = 1, n
415*
416 IF( somev ) THEN
417 IF( .NOT.SELECT( ki ) )
418 $ GO TO 130
419 END IF
420 smin = max( ulp*( cabs1( t( ki, ki ) ) ), smlnum )
421*
422 work( n ) = cmone
423*
424* Form right-hand side.
425*
426 DO 90 k = ki + 1, n
427 work( k ) = -dconjg( t( ki, k ) )
428 90 CONTINUE
429*
430* Solve the triangular system:
431* (T(KI+1:N,KI+1:N) - T(KI,KI))**H * X = SCALE*WORK.
432*
433 DO 100 k = ki + 1, n
434 t( k, k ) = t( k, k ) - t( ki, ki )
435 IF( cabs1( t( k, k ) ).LT.smin )
436 $ t( k, k ) = smin
437 100 CONTINUE
438*
439 IF( ki.LT.n ) THEN
440 CALL zlatrs( 'Upper', 'Conjugate transpose',
441 $ 'Non-unit',
442 $ 'Y', n-ki, t( ki+1, ki+1 ), ldt,
443 $ work( ki+1 ), scale, rwork, info )
444 work( ki ) = scale
445 END IF
446*
447* Copy the vector x or Q*x to VL and normalize.
448*
449 IF( .NOT.over ) THEN
450 CALL zcopy( n-ki+1, work( ki ), 1, vl( ki, is ), 1 )
451*
452 ii = izamax( n-ki+1, vl( ki, is ), 1 ) + ki - 1
453 remax = one / cabs1( vl( ii, is ) )
454 CALL zdscal( n-ki+1, remax, vl( ki, is ), 1 )
455*
456 DO 110 k = 1, ki - 1
457 vl( k, is ) = cmzero
458 110 CONTINUE
459 ELSE
460 IF( ki.LT.n )
461 $ CALL zgemv( 'N', n, n-ki, cmone, vl( 1, ki+1 ),
462 $ ldvl,
463 $ work( ki+1 ), 1, dcmplx( scale ),
464 $ vl( 1, ki ), 1 )
465*
466 ii = izamax( n, vl( 1, ki ), 1 )
467 remax = one / cabs1( vl( ii, ki ) )
468 CALL zdscal( n, remax, vl( 1, ki ), 1 )
469 END IF
470*
471* Set back the original diagonal elements of T.
472*
473 DO 120 k = ki + 1, n
474 t( k, k ) = work( k+n )
475 120 CONTINUE
476*
477 is = is + 1
478 130 CONTINUE
479 END IF
480*
481 RETURN
482*
483* End of ZTREVC
484*
485 END
subroutine xerbla(srname, info)
Definition cblat2.f:3285
subroutine zcopy(n, zx, incx, zy, incy)
ZCOPY
Definition zcopy.f:81
subroutine zgemv(trans, m, n, alpha, a, lda, x, incx, beta, y, incy)
ZGEMV
Definition zgemv.f:160
subroutine zlatrs(uplo, trans, diag, normin, n, a, lda, x, scale, cnorm, info)
ZLATRS solves a triangular system of equations with the scale factor set to prevent overflow.
Definition zlatrs.f:238
subroutine zdscal(n, da, zx, incx)
ZDSCAL
Definition zdscal.f:78
subroutine ztrevc(side, howmny, select, n, t, ldt, vl, ldvl, vr, ldvr, mm, m, work, rwork, info)
ZTREVC
Definition ztrevc.f:217