LAPACK  3.4.2
LAPACK: Linear Algebra PACKage
 All Files Functions Groups
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 *> \htmlonly
9 *> Download ZTREVC + dependencies
10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/ztrevc.f">
11 *> [TGZ]</a>
12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/ztrevc.f">
13 *> [ZIP]</a>
14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/ztrevc.f">
15 *> [TXT]</a>
16 *> \endhtmlonly
17 *
18 * Definition:
19 * ===========
20 *
21 * SUBROUTINE ZTREVC( 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 * DOUBLE PRECISION RWORK( * )
31 * COMPLEX*16 T( LDT, * ), VL( LDVL, * ), VR( LDVR, * ),
32 * $ WORK( * )
33 * ..
34 *
35 *
36 *> \par Purpose:
37 * =============
38 *>
39 *> \verbatim
40 *>
41 *> ZTREVC 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 ZHSEQR.
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*16 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*16 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 ZHSEQR).
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*16 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 ZHSEQR).
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*16 array, dimension (2*N)
177 *> \endverbatim
178 *>
179 *> \param[out] RWORK
180 *> \verbatim
181 *> RWORK is DOUBLE PRECISION 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 *> \date November 2011
200 *
201 *> \ingroup complex16OTHERcomputational
202 *
203 *> \par Further Details:
204 * =====================
205 *>
206 *> \verbatim
207 *>
208 *> The algorithm used in this program is basically backward (forward)
209 *> substitution, with scaling to make the the code robust against
210 *> possible overflow.
211 *>
212 *> Each eigenvector is normalized so that the element of largest
213 *> magnitude has magnitude 1; here the magnitude of a complex number
214 *> (x,y) is taken to be |x| + |y|.
215 *> \endverbatim
216 *>
217 * =====================================================================
218  SUBROUTINE ztrevc( SIDE, HOWMNY, SELECT, N, T, LDT, VL, LDVL, VR,
219  $ ldvr, mm, m, work, rwork, info )
220 *
221 * -- LAPACK computational routine (version 3.4.0) --
222 * -- LAPACK is a software package provided by Univ. of Tennessee, --
223 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
224 * November 2011
225 *
226 * .. Scalar Arguments ..
227  CHARACTER howmny, side
228  INTEGER info, ldt, ldvl, ldvr, m, mm, n
229 * ..
230 * .. Array Arguments ..
231  LOGICAL select( * )
232  DOUBLE PRECISION rwork( * )
233  COMPLEX*16 t( ldt, * ), vl( ldvl, * ), vr( ldvr, * ),
234  $ work( * )
235 * ..
236 *
237 * =====================================================================
238 *
239 * .. Parameters ..
240  DOUBLE PRECISION zero, one
241  parameter( zero = 0.0d+0, one = 1.0d+0 )
242  COMPLEX*16 cmzero, cmone
243  parameter( cmzero = ( 0.0d+0, 0.0d+0 ),
244  $ cmone = ( 1.0d+0, 0.0d+0 ) )
245 * ..
246 * .. Local Scalars ..
247  LOGICAL allv, bothv, leftv, over, rightv, somev
248  INTEGER i, ii, is, j, k, ki
249  DOUBLE PRECISION ovfl, remax, scale, smin, smlnum, ulp, unfl
250  COMPLEX*16 cdum
251 * ..
252 * .. External Functions ..
253  LOGICAL lsame
254  INTEGER izamax
255  DOUBLE PRECISION dlamch, dzasum
256  EXTERNAL lsame, izamax, dlamch, dzasum
257 * ..
258 * .. External Subroutines ..
259  EXTERNAL xerbla, zcopy, zdscal, zgemv, zlatrs
260 * ..
261 * .. Intrinsic Functions ..
262  INTRINSIC abs, dble, dcmplx, dconjg, dimag, max
263 * ..
264 * .. Statement Functions ..
265  DOUBLE PRECISION cabs1
266 * ..
267 * .. Statement Function definitions ..
268  cabs1( cdum ) = abs( dble( cdum ) ) + abs( dimag( cdum ) )
269 * ..
270 * .. Executable Statements ..
271 *
272 * Decode and test the input parameters
273 *
274  bothv = lsame( side, 'B' )
275  rightv = lsame( side, 'R' ) .OR. bothv
276  leftv = lsame( side, 'L' ) .OR. bothv
277 *
278  allv = lsame( howmny, 'A' )
279  over = lsame( howmny, 'B' )
280  somev = lsame( howmny, 'S' )
281 *
282 * Set M to the number of columns required to store the selected
283 * eigenvectors.
284 *
285  IF( somev ) THEN
286  m = 0
287  DO 10 j = 1, n
288  IF( SELECT( j ) )
289  $ m = m + 1
290  10 continue
291  ELSE
292  m = n
293  END IF
294 *
295  info = 0
296  IF( .NOT.rightv .AND. .NOT.leftv ) THEN
297  info = -1
298  ELSE IF( .NOT.allv .AND. .NOT.over .AND. .NOT.somev ) THEN
299  info = -2
300  ELSE IF( n.LT.0 ) THEN
301  info = -4
302  ELSE IF( ldt.LT.max( 1, n ) ) THEN
303  info = -6
304  ELSE IF( ldvl.LT.1 .OR. ( leftv .AND. ldvl.LT.n ) ) THEN
305  info = -8
306  ELSE IF( ldvr.LT.1 .OR. ( rightv .AND. ldvr.LT.n ) ) THEN
307  info = -10
308  ELSE IF( mm.LT.m ) THEN
309  info = -11
310  END IF
311  IF( info.NE.0 ) THEN
312  CALL xerbla( 'ZTREVC', -info )
313  return
314  END IF
315 *
316 * Quick return if possible.
317 *
318  IF( n.EQ.0 )
319  $ return
320 *
321 * Set the constants to control overflow.
322 *
323  unfl = dlamch( 'Safe minimum' )
324  ovfl = one / unfl
325  CALL dlabad( unfl, ovfl )
326  ulp = dlamch( 'Precision' )
327  smlnum = unfl*( n / ulp )
328 *
329 * Store the diagonal elements of T in working array WORK.
330 *
331  DO 20 i = 1, n
332  work( i+n ) = t( i, i )
333  20 continue
334 *
335 * Compute 1-norm of each column of strictly upper triangular
336 * part of T to control overflow in triangular solver.
337 *
338  rwork( 1 ) = zero
339  DO 30 j = 2, n
340  rwork( j ) = dzasum( j-1, t( 1, j ), 1 )
341  30 continue
342 *
343  IF( rightv ) THEN
344 *
345 * Compute right eigenvectors.
346 *
347  is = m
348  DO 80 ki = n, 1, -1
349 *
350  IF( somev ) THEN
351  IF( .NOT.SELECT( ki ) )
352  $ go to 80
353  END IF
354  smin = max( ulp*( cabs1( t( ki, ki ) ) ), smlnum )
355 *
356  work( 1 ) = cmone
357 *
358 * Form right-hand side.
359 *
360  DO 40 k = 1, ki - 1
361  work( k ) = -t( k, ki )
362  40 continue
363 *
364 * Solve the triangular system:
365 * (T(1:KI-1,1:KI-1) - T(KI,KI))*X = SCALE*WORK.
366 *
367  DO 50 k = 1, ki - 1
368  t( k, k ) = t( k, k ) - t( ki, ki )
369  IF( cabs1( t( k, k ) ).LT.smin )
370  $ t( k, k ) = smin
371  50 continue
372 *
373  IF( ki.GT.1 ) THEN
374  CALL zlatrs( 'Upper', 'No transpose', 'Non-unit', 'Y',
375  $ ki-1, t, ldt, work( 1 ), scale, rwork,
376  $ info )
377  work( ki ) = scale
378  END IF
379 *
380 * Copy the vector x or Q*x to VR and normalize.
381 *
382  IF( .NOT.over ) THEN
383  CALL zcopy( ki, work( 1 ), 1, vr( 1, is ), 1 )
384 *
385  ii = izamax( ki, vr( 1, is ), 1 )
386  remax = one / cabs1( vr( ii, is ) )
387  CALL zdscal( ki, remax, vr( 1, is ), 1 )
388 *
389  DO 60 k = ki + 1, n
390  vr( k, is ) = cmzero
391  60 continue
392  ELSE
393  IF( ki.GT.1 )
394  $ CALL zgemv( 'N', n, ki-1, cmone, vr, ldvr, work( 1 ),
395  $ 1, dcmplx( scale ), vr( 1, ki ), 1 )
396 *
397  ii = izamax( n, vr( 1, ki ), 1 )
398  remax = one / cabs1( vr( ii, ki ) )
399  CALL zdscal( n, remax, vr( 1, ki ), 1 )
400  END IF
401 *
402 * Set back the original diagonal elements of T.
403 *
404  DO 70 k = 1, ki - 1
405  t( k, k ) = work( k+n )
406  70 continue
407 *
408  is = is - 1
409  80 continue
410  END IF
411 *
412  IF( leftv ) THEN
413 *
414 * Compute left eigenvectors.
415 *
416  is = 1
417  DO 130 ki = 1, n
418 *
419  IF( somev ) THEN
420  IF( .NOT.SELECT( ki ) )
421  $ go to 130
422  END IF
423  smin = max( ulp*( cabs1( t( ki, ki ) ) ), smlnum )
424 *
425  work( n ) = cmone
426 *
427 * Form right-hand side.
428 *
429  DO 90 k = ki + 1, n
430  work( k ) = -dconjg( t( ki, k ) )
431  90 continue
432 *
433 * Solve the triangular system:
434 * (T(KI+1:N,KI+1:N) - T(KI,KI))**H * X = SCALE*WORK.
435 *
436  DO 100 k = ki + 1, n
437  t( k, k ) = t( k, k ) - t( ki, ki )
438  IF( cabs1( t( k, k ) ).LT.smin )
439  $ t( k, k ) = smin
440  100 continue
441 *
442  IF( ki.LT.n ) THEN
443  CALL zlatrs( 'Upper', 'Conjugate transpose', 'Non-unit',
444  $ 'Y', n-ki, t( ki+1, ki+1 ), ldt,
445  $ work( ki+1 ), scale, rwork, info )
446  work( ki ) = scale
447  END IF
448 *
449 * Copy the vector x or Q*x to VL and normalize.
450 *
451  IF( .NOT.over ) THEN
452  CALL zcopy( n-ki+1, work( ki ), 1, vl( ki, is ), 1 )
453 *
454  ii = izamax( n-ki+1, vl( ki, is ), 1 ) + ki - 1
455  remax = one / cabs1( vl( ii, is ) )
456  CALL zdscal( n-ki+1, remax, vl( ki, is ), 1 )
457 *
458  DO 110 k = 1, ki - 1
459  vl( k, is ) = cmzero
460  110 continue
461  ELSE
462  IF( ki.LT.n )
463  $ CALL zgemv( 'N', n, n-ki, cmone, vl( 1, ki+1 ), ldvl,
464  $ work( ki+1 ), 1, dcmplx( scale ),
465  $ vl( 1, ki ), 1 )
466 *
467  ii = izamax( n, vl( 1, ki ), 1 )
468  remax = one / cabs1( vl( ii, ki ) )
469  CALL zdscal( n, remax, vl( 1, ki ), 1 )
470  END IF
471 *
472 * Set back the original diagonal elements of T.
473 *
474  DO 120 k = ki + 1, n
475  t( k, k ) = work( k+n )
476  120 continue
477 *
478  is = is + 1
479  130 continue
480  END IF
481 *
482  return
483 *
484 * End of ZTREVC
485 *
486  END