LAPACK 3.12.0
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches
ctrevc.f
Go to the documentation of this file.
1*> \brief \b CTREVC
2*
3* =========== DOCUMENTATION ===========
4*
5* Online html documentation available at
6* http://www.netlib.org/lapack/explore-html/
7*
8*> \htmlonly
9*> Download CTREVC + dependencies
10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/ctrevc.f">
11*> [TGZ]</a>
12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/ctrevc.f">
13*> [ZIP]</a>
14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/ctrevc.f">
15*> [TXT]</a>
16*> \endhtmlonly
17*
18* Definition:
19* ===========
20*
21* SUBROUTINE CTREVC( SIDE, HOWMNY, SELECT, N, T, LDT, VL, LDVL, VR,
22* LDVR, MM, M, WORK, RWORK, INFO )
23*
24* .. Scalar Arguments ..
25* CHARACTER HOWMNY, SIDE
26* INTEGER INFO, LDT, LDVL, LDVR, M, MM, N
27* ..
28* .. Array Arguments ..
29* LOGICAL SELECT( * )
30* REAL RWORK( * )
31* COMPLEX T( LDT, * ), VL( LDVL, * ), VR( LDVR, * ),
32* $ WORK( * )
33* ..
34*
35*
36*> \par Purpose:
37* =============
38*>
39*> \verbatim
40*>
41*> CTREVC computes some or all of the right and/or left eigenvectors of
42*> a complex upper triangular matrix T.
43*> Matrices of this type are produced by the Schur factorization of
44*> a complex general matrix: A = Q*T*Q**H, as computed by CHSEQR.
45*>
46*> The right eigenvector x and the left eigenvector y of T corresponding
47*> to an eigenvalue w are defined by:
48*>
49*> T*x = w*x, (y**H)*T = w*(y**H)
50*>
51*> where y**H denotes the conjugate transpose of the vector y.
52*> The eigenvalues are not input to this routine, but are read directly
53*> from the diagonal of T.
54*>
55*> This routine returns the matrices X and/or Y of right and left
56*> eigenvectors of T, or the products Q*X and/or Q*Y, where Q is an
57*> input matrix. If Q is the unitary factor that reduces a matrix A to
58*> Schur form T, then Q*X and Q*Y are the matrices of right and left
59*> eigenvectors of A.
60*> \endverbatim
61*
62* Arguments:
63* ==========
64*
65*> \param[in] SIDE
66*> \verbatim
67*> SIDE is CHARACTER*1
68*> = 'R': compute right eigenvectors only;
69*> = 'L': compute left eigenvectors only;
70*> = 'B': compute both right and left eigenvectors.
71*> \endverbatim
72*>
73*> \param[in] HOWMNY
74*> \verbatim
75*> HOWMNY is CHARACTER*1
76*> = 'A': compute all right and/or left eigenvectors;
77*> = 'B': compute all right and/or left eigenvectors,
78*> backtransformed using the matrices supplied in
79*> VR and/or VL;
80*> = 'S': compute selected right and/or left eigenvectors,
81*> as indicated by the logical array SELECT.
82*> \endverbatim
83*>
84*> \param[in] SELECT
85*> \verbatim
86*> SELECT is LOGICAL array, dimension (N)
87*> If HOWMNY = 'S', SELECT specifies the eigenvectors to be
88*> computed.
89*> The eigenvector corresponding to the j-th eigenvalue is
90*> computed if SELECT(j) = .TRUE..
91*> Not referenced if HOWMNY = 'A' or 'B'.
92*> \endverbatim
93*>
94*> \param[in] N
95*> \verbatim
96*> N is INTEGER
97*> The order of the matrix T. N >= 0.
98*> \endverbatim
99*>
100*> \param[in,out] T
101*> \verbatim
102*> T is COMPLEX array, dimension (LDT,N)
103*> The upper triangular matrix T. T is modified, but restored
104*> on exit.
105*> \endverbatim
106*>
107*> \param[in] LDT
108*> \verbatim
109*> LDT is INTEGER
110*> The leading dimension of the array T. LDT >= max(1,N).
111*> \endverbatim
112*>
113*> \param[in,out] VL
114*> \verbatim
115*> VL is COMPLEX array, dimension (LDVL,MM)
116*> On entry, if SIDE = 'L' or 'B' and HOWMNY = 'B', VL must
117*> contain an N-by-N matrix Q (usually the unitary matrix Q of
118*> Schur vectors returned by CHSEQR).
119*> On exit, if SIDE = 'L' or 'B', VL contains:
120*> if HOWMNY = 'A', the matrix Y of left eigenvectors of T;
121*> if HOWMNY = 'B', the matrix Q*Y;
122*> if HOWMNY = 'S', the left eigenvectors of T specified by
123*> SELECT, stored consecutively in the columns
124*> of VL, in the same order as their
125*> eigenvalues.
126*> Not referenced if SIDE = 'R'.
127*> \endverbatim
128*>
129*> \param[in] LDVL
130*> \verbatim
131*> LDVL is INTEGER
132*> The leading dimension of the array VL. LDVL >= 1, and if
133*> SIDE = 'L' or 'B', LDVL >= N.
134*> \endverbatim
135*>
136*> \param[in,out] VR
137*> \verbatim
138*> VR is COMPLEX array, dimension (LDVR,MM)
139*> On entry, if SIDE = 'R' or 'B' and HOWMNY = 'B', VR must
140*> contain an N-by-N matrix Q (usually the unitary matrix Q of
141*> Schur vectors returned by CHSEQR).
142*> On exit, if SIDE = 'R' or 'B', VR contains:
143*> if HOWMNY = 'A', the matrix X of right eigenvectors of T;
144*> if HOWMNY = 'B', the matrix Q*X;
145*> if HOWMNY = 'S', the right eigenvectors of T specified by
146*> SELECT, stored consecutively in the columns
147*> of VR, in the same order as their
148*> eigenvalues.
149*> Not referenced if SIDE = 'L'.
150*> \endverbatim
151*>
152*> \param[in] LDVR
153*> \verbatim
154*> LDVR is INTEGER
155*> The leading dimension of the array VR. LDVR >= 1, and if
156*> SIDE = 'R' or 'B'; LDVR >= N.
157*> \endverbatim
158*>
159*> \param[in] MM
160*> \verbatim
161*> MM is INTEGER
162*> The number of columns in the arrays VL and/or VR. MM >= M.
163*> \endverbatim
164*>
165*> \param[out] M
166*> \verbatim
167*> M is INTEGER
168*> The number of columns in the arrays VL and/or VR actually
169*> used to store the eigenvectors. If HOWMNY = 'A' or 'B', M
170*> is set to N. Each selected eigenvector occupies one
171*> column.
172*> \endverbatim
173*>
174*> \param[out] WORK
175*> \verbatim
176*> WORK is COMPLEX array, dimension (2*N)
177*> \endverbatim
178*>
179*> \param[out] RWORK
180*> \verbatim
181*> RWORK is REAL array, dimension (N)
182*> \endverbatim
183*>
184*> \param[out] INFO
185*> \verbatim
186*> INFO is INTEGER
187*> = 0: successful exit
188*> < 0: if INFO = -i, the i-th argument had an illegal value
189*> \endverbatim
190*
191* Authors:
192* ========
193*
194*> \author Univ. of Tennessee
195*> \author Univ. of California Berkeley
196*> \author Univ. of Colorado Denver
197*> \author NAG Ltd.
198*
199*> \ingroup trevc
200*
201*> \par Further Details:
202* =====================
203*>
204*> \verbatim
205*>
206*> The algorithm used in this program is basically backward (forward)
207*> substitution, with scaling to make the the code robust against
208*> possible overflow.
209*>
210*> Each eigenvector is normalized so that the element of largest
211*> magnitude has magnitude 1; here the magnitude of a complex number
212*> (x,y) is taken to be |x| + |y|.
213*> \endverbatim
214*>
215* =====================================================================
216 SUBROUTINE ctrevc( SIDE, HOWMNY, SELECT, N, T, LDT, VL, LDVL, VR,
217 $ LDVR, MM, M, WORK, RWORK, INFO )
218*
219* -- LAPACK computational routine --
220* -- LAPACK is a software package provided by Univ. of Tennessee, --
221* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
222*
223* .. Scalar Arguments ..
224 CHARACTER HOWMNY, SIDE
225 INTEGER INFO, LDT, LDVL, LDVR, M, MM, N
226* ..
227* .. Array Arguments ..
228 LOGICAL SELECT( * )
229 REAL RWORK( * )
230 COMPLEX T( LDT, * ), VL( LDVL, * ), VR( LDVR, * ),
231 $ work( * )
232* ..
233*
234* =====================================================================
235*
236* .. Parameters ..
237 REAL ZERO, ONE
238 parameter( zero = 0.0e+0, one = 1.0e+0 )
239 COMPLEX CMZERO, CMONE
240 parameter( cmzero = ( 0.0e+0, 0.0e+0 ),
241 $ cmone = ( 1.0e+0, 0.0e+0 ) )
242* ..
243* .. Local Scalars ..
244 LOGICAL ALLV, BOTHV, LEFTV, OVER, RIGHTV, SOMEV
245 INTEGER I, II, IS, J, K, KI
246 REAL OVFL, REMAX, SCALE, SMIN, SMLNUM, ULP, UNFL
247 COMPLEX CDUM
248* ..
249* .. External Functions ..
250 LOGICAL LSAME
251 INTEGER ICAMAX
252 REAL SCASUM, SLAMCH
253 EXTERNAL lsame, icamax, scasum, slamch
254* ..
255* .. External Subroutines ..
256 EXTERNAL ccopy, cgemv, clatrs, csscal, xerbla
257* ..
258* .. Intrinsic Functions ..
259 INTRINSIC abs, aimag, cmplx, conjg, max, real
260* ..
261* .. Statement Functions ..
262 REAL CABS1
263* ..
264* .. Statement Function definitions ..
265 cabs1( cdum ) = abs( real( cdum ) ) + abs( aimag( 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( 'CTREVC', -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 = slamch( 'Safe minimum' )
321 ovfl = one / unfl
322 ulp = slamch( '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 ) = scasum( 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 clatrs( '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 ccopy( ki, work( 1 ), 1, vr( 1, is ), 1 )
380*
381 ii = icamax( ki, vr( 1, is ), 1 )
382 remax = one / cabs1( vr( ii, is ) )
383 CALL csscal( 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 cgemv( 'N', n, ki-1, cmone, vr, ldvr, work( 1 ),
391 $ 1, cmplx( scale ), vr( 1, ki ), 1 )
392*
393 ii = icamax( n, vr( 1, ki ), 1 )
394 remax = one / cabs1( vr( ii, ki ) )
395 CALL csscal( n, remax, vr( 1, ki ), 1 )
396 END IF
397*
398* Set back the original diagonal elements of T.
399*
400 DO 70 k = 1, ki - 1
401 t( k, k ) = work( k+n )
402 70 CONTINUE
403*
404 is = is - 1
405 80 CONTINUE
406 END IF
407*
408 IF( leftv ) THEN
409*
410* Compute left eigenvectors.
411*
412 is = 1
413 DO 130 ki = 1, n
414*
415 IF( somev ) THEN
416 IF( .NOT.SELECT( ki ) )
417 $ GO TO 130
418 END IF
419 smin = max( ulp*( cabs1( t( ki, ki ) ) ), smlnum )
420*
421 work( n ) = cmone
422*
423* Form right-hand side.
424*
425 DO 90 k = ki + 1, n
426 work( k ) = -conjg( t( ki, k ) )
427 90 CONTINUE
428*
429* Solve the triangular system:
430* (T(KI+1:N,KI+1:N) - T(KI,KI))**H*X = SCALE*WORK.
431*
432 DO 100 k = ki + 1, n
433 t( k, k ) = t( k, k ) - t( ki, ki )
434 IF( cabs1( t( k, k ) ).LT.smin )
435 $ t( k, k ) = smin
436 100 CONTINUE
437*
438 IF( ki.LT.n ) THEN
439 CALL clatrs( 'Upper', 'Conjugate transpose', 'Non-unit',
440 $ 'Y', n-ki, t( ki+1, ki+1 ), ldt,
441 $ work( ki+1 ), scale, rwork, info )
442 work( ki ) = scale
443 END IF
444*
445* Copy the vector x or Q*x to VL and normalize.
446*
447 IF( .NOT.over ) THEN
448 CALL ccopy( n-ki+1, work( ki ), 1, vl( ki, is ), 1 )
449*
450 ii = icamax( n-ki+1, vl( ki, is ), 1 ) + ki - 1
451 remax = one / cabs1( vl( ii, is ) )
452 CALL csscal( n-ki+1, remax, vl( ki, is ), 1 )
453*
454 DO 110 k = 1, ki - 1
455 vl( k, is ) = cmzero
456 110 CONTINUE
457 ELSE
458 IF( ki.LT.n )
459 $ CALL cgemv( 'N', n, n-ki, cmone, vl( 1, ki+1 ), ldvl,
460 $ work( ki+1 ), 1, cmplx( scale ),
461 $ vl( 1, ki ), 1 )
462*
463 ii = icamax( n, vl( 1, ki ), 1 )
464 remax = one / cabs1( vl( ii, ki ) )
465 CALL csscal( n, remax, vl( 1, ki ), 1 )
466 END IF
467*
468* Set back the original diagonal elements of T.
469*
470 DO 120 k = ki + 1, n
471 t( k, k ) = work( k+n )
472 120 CONTINUE
473*
474 is = is + 1
475 130 CONTINUE
476 END IF
477*
478 RETURN
479*
480* End of CTREVC
481*
482 END
subroutine xerbla(srname, info)
Definition cblat2.f:3285
subroutine ccopy(n, cx, incx, cy, incy)
CCOPY
Definition ccopy.f:81
subroutine cgemv(trans, m, n, alpha, a, lda, x, incx, beta, y, incy)
CGEMV
Definition cgemv.f:160
subroutine clatrs(uplo, trans, diag, normin, n, a, lda, x, scale, cnorm, info)
CLATRS solves a triangular system of equations with the scale factor set to prevent overflow.
Definition clatrs.f:239
subroutine csscal(n, sa, cx, incx)
CSSCAL
Definition csscal.f:78
subroutine ctrevc(side, howmny, select, n, t, ldt, vl, ldvl, vr, ldvr, mm, m, work, rwork, info)
CTREVC
Definition ctrevc.f:218