1: *> \brief \b ZTREVC3
2: *
3: * =========== DOCUMENTATION ===========
4: *
5: * Online html documentation available at
6: * http://www.netlib.org/lapack/explore-html/
7: *
8: *> \htmlonly
9: *> Download ZTREVC3 + dependencies
10: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/ztrevc3.f">
11: *> [TGZ]</a>
12: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/ztrevc3.f">
13: *> [ZIP]</a>
14: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/ztrevc3.f">
15: *> [TXT]</a>
16: *> \endhtmlonly
17: *
18: * Definition:
19: * ===========
20: *
21: * SUBROUTINE ZTREVC3( SIDE, HOWMNY, SELECT, N, T, LDT, VL, LDVL,
22: * VR, LDVR, MM, M, WORK, LWORK, RWORK, INFO )
23: *
24: * .. Scalar Arguments ..
25: * CHARACTER HOWMNY, SIDE
26: * INTEGER INFO, LDT, LDVL, LDVR, LWORK, 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: *> ZTREVC3 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: *>
61: *> This uses a Level 3 BLAS version of the back transformation.
62: *> \endverbatim
63: *
64: * Arguments:
65: * ==========
66: *
67: *> \param[in] SIDE
68: *> \verbatim
69: *> SIDE is CHARACTER*1
70: *> = 'R': compute right eigenvectors only;
71: *> = 'L': compute left eigenvectors only;
72: *> = 'B': compute both right and left eigenvectors.
73: *> \endverbatim
74: *>
75: *> \param[in] HOWMNY
76: *> \verbatim
77: *> HOWMNY is CHARACTER*1
78: *> = 'A': compute all right and/or left eigenvectors;
79: *> = 'B': compute all right and/or left eigenvectors,
80: *> backtransformed using the matrices supplied in
81: *> VR and/or VL;
82: *> = 'S': compute selected right and/or left eigenvectors,
83: *> as indicated by the logical array SELECT.
84: *> \endverbatim
85: *>
86: *> \param[in] SELECT
87: *> \verbatim
88: *> SELECT is LOGICAL array, dimension (N)
89: *> If HOWMNY = 'S', SELECT specifies the eigenvectors to be
90: *> computed.
91: *> The eigenvector corresponding to the j-th eigenvalue is
92: *> computed if SELECT(j) = .TRUE..
93: *> Not referenced if HOWMNY = 'A' or 'B'.
94: *> \endverbatim
95: *>
96: *> \param[in] N
97: *> \verbatim
98: *> N is INTEGER
99: *> The order of the matrix T. N >= 0.
100: *> \endverbatim
101: *>
102: *> \param[in,out] T
103: *> \verbatim
104: *> T is COMPLEX*16 array, dimension (LDT,N)
105: *> The upper triangular matrix T. T is modified, but restored
106: *> on exit.
107: *> \endverbatim
108: *>
109: *> \param[in] LDT
110: *> \verbatim
111: *> LDT is INTEGER
112: *> The leading dimension of the array T. LDT >= max(1,N).
113: *> \endverbatim
114: *>
115: *> \param[in,out] VL
116: *> \verbatim
117: *> VL is COMPLEX*16 array, dimension (LDVL,MM)
118: *> On entry, if SIDE = 'L' or 'B' and HOWMNY = 'B', VL must
119: *> contain an N-by-N matrix Q (usually the unitary matrix Q of
120: *> Schur vectors returned by ZHSEQR).
121: *> On exit, if SIDE = 'L' or 'B', VL contains:
122: *> if HOWMNY = 'A', the matrix Y of left eigenvectors of T;
123: *> if HOWMNY = 'B', the matrix Q*Y;
124: *> if HOWMNY = 'S', the left eigenvectors of T specified by
125: *> SELECT, stored consecutively in the columns
126: *> of VL, in the same order as their
127: *> eigenvalues.
128: *> Not referenced if SIDE = 'R'.
129: *> \endverbatim
130: *>
131: *> \param[in] LDVL
132: *> \verbatim
133: *> LDVL is INTEGER
134: *> The leading dimension of the array VL.
135: *> LDVL >= 1, and if SIDE = 'L' or 'B', LDVL >= N.
136: *> \endverbatim
137: *>
138: *> \param[in,out] VR
139: *> \verbatim
140: *> VR is COMPLEX*16 array, dimension (LDVR,MM)
141: *> On entry, if SIDE = 'R' or 'B' and HOWMNY = 'B', VR must
142: *> contain an N-by-N matrix Q (usually the unitary matrix Q of
143: *> Schur vectors returned by ZHSEQR).
144: *> On exit, if SIDE = 'R' or 'B', VR contains:
145: *> if HOWMNY = 'A', the matrix X of right eigenvectors of T;
146: *> if HOWMNY = 'B', the matrix Q*X;
147: *> if HOWMNY = 'S', the right eigenvectors of T specified by
148: *> SELECT, stored consecutively in the columns
149: *> of VR, in the same order as their
150: *> eigenvalues.
151: *> Not referenced if SIDE = 'L'.
152: *> \endverbatim
153: *>
154: *> \param[in] LDVR
155: *> \verbatim
156: *> LDVR is INTEGER
157: *> The leading dimension of the array VR.
158: *> LDVR >= 1, and if SIDE = 'R' or 'B', LDVR >= N.
159: *> \endverbatim
160: *>
161: *> \param[in] MM
162: *> \verbatim
163: *> MM is INTEGER
164: *> The number of columns in the arrays VL and/or VR. MM >= M.
165: *> \endverbatim
166: *>
167: *> \param[out] M
168: *> \verbatim
169: *> M is INTEGER
170: *> The number of columns in the arrays VL and/or VR actually
171: *> used to store the eigenvectors.
172: *> If HOWMNY = 'A' or 'B', M is set to N.
173: *> Each selected eigenvector occupies one column.
174: *> \endverbatim
175: *>
176: *> \param[out] WORK
177: *> \verbatim
178: *> WORK is COMPLEX*16 array, dimension (MAX(1,LWORK))
179: *> \endverbatim
180: *>
181: *> \param[in] LWORK
182: *> \verbatim
183: *> LWORK is INTEGER
184: *> The dimension of array WORK. LWORK >= max(1,2*N).
185: *> For optimum performance, LWORK >= N + 2*N*NB, where NB is
186: *> the optimal blocksize.
187: *>
188: *> If LWORK = -1, then a workspace query is assumed; the routine
189: *> only calculates the optimal size of the WORK array, returns
190: *> this value as the first entry of the WORK array, and no error
191: *> message related to LWORK is issued by XERBLA.
192: *> \endverbatim
193: *>
194: *> \param[out] RWORK
195: *> \verbatim
196: *> RWORK is DOUBLE PRECISION array, dimension (LRWORK)
197: *> \endverbatim
198: *>
199: *> \param[in] LRWORK
200: *> \verbatim
201: *> LRWORK is INTEGER
202: *> The dimension of array RWORK. LRWORK >= max(1,N).
203: *>
204: *> If LRWORK = -1, then a workspace query is assumed; the routine
205: *> only calculates the optimal size of the RWORK array, returns
206: *> this value as the first entry of the RWORK array, and no error
207: *> message related to LRWORK is issued by XERBLA.
208: *> \endverbatim
209: *>
210: *> \param[out] INFO
211: *> \verbatim
212: *> INFO is INTEGER
213: *> = 0: successful exit
214: *> < 0: if INFO = -i, the i-th argument had an illegal value
215: *> \endverbatim
216: *
217: * Authors:
218: * ========
219: *
220: *> \author Univ. of Tennessee
221: *> \author Univ. of California Berkeley
222: *> \author Univ. of Colorado Denver
223: *> \author NAG Ltd.
224: *
225: *> \date November 2011
226: *
227: * @precisions fortran z -> c
228: *
229: *> \ingroup complex16OTHERcomputational
230: *
231: *> \par Further Details:
232: * =====================
233: *>
234: *> \verbatim
235: *>
236: *> The algorithm used in this program is basically backward (forward)
237: *> substitution, with scaling to make the the code robust against
238: *> possible overflow.
239: *>
240: *> Each eigenvector is normalized so that the element of largest
241: *> magnitude has magnitude 1; here the magnitude of a complex number
242: *> (x,y) is taken to be |x| + |y|.
243: *> \endverbatim
244: *>
245: * =====================================================================
246: SUBROUTINE ZTREVC3( SIDE, HOWMNY, SELECT, N, T, LDT, VL, LDVL, VR,
247: $ LDVR, MM, M, WORK, LWORK, RWORK, LRWORK, INFO)
248: IMPLICIT NONE
249: *
250: * -- LAPACK computational routine (version 3.4.0) --
251: * -- LAPACK is a software package provided by Univ. of Tennessee, --
252: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
253: * November 2011
254: *
255: * .. Scalar Arguments ..
256: CHARACTER HOWMNY, SIDE
257: INTEGER INFO, LDT, LDVL, LDVR, LWORK, LRWORK, M, MM, N
258: * ..
259: * .. Array Arguments ..
260: LOGICAL SELECT( * )
261: DOUBLE PRECISION RWORK( * )
262: COMPLEX*16 T( LDT, * ), VL( LDVL, * ), VR( LDVR, * ),
263: $ WORK( * )
264: * ..
265: *
266: * =====================================================================
267: *
268: * .. Parameters ..
269: DOUBLE PRECISION ZERO, ONE
270: PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 )
271: COMPLEX*16 CZERO, CONE
272: PARAMETER ( CZERO = ( 0.0D+0, 0.0D+0 ),
273: $ CONE = ( 1.0D+0, 0.0D+0 ) )
274: INTEGER NBMIN, NBMAX
275: PARAMETER ( NBMIN = 8, NBMAX = 128 )
276: * ..
277: * .. Local Scalars ..
278: LOGICAL ALLV, BOTHV, LEFTV, LQUERY, OVER, RIGHTV, SOMEV
279: INTEGER I, II, IS, J, K, KI, IV, MAXWRK, NB
280: DOUBLE PRECISION OVFL, REMAX, SCALE, SMIN, SMLNUM, ULP, UNFL
281: COMPLEX*16 CDUM
282: * ..
283: * .. External Functions ..
284: LOGICAL LSAME
285: INTEGER ILAENV, IZAMAX
286: DOUBLE PRECISION DLAMCH, DZASUM
287: EXTERNAL LSAME, ILAENV, IZAMAX, DLAMCH, DZASUM
288: * ..
289: * .. External Subroutines ..
290: EXTERNAL XERBLA, ZCOPY, ZDSCAL, ZGEMV, ZLATRS
291: * ..
292: * .. Intrinsic Functions ..
293: INTRINSIC ABS, DBLE, DCMPLX, CONJG, AIMAG, MAX
294: * ..
295: * .. Statement Functions ..
296: DOUBLE PRECISION CABS1
297: * ..
298: * .. Statement Function definitions ..
299: CABS1( CDUM ) = ABS( DBLE( CDUM ) ) + ABS( AIMAG( CDUM ) )
300: * ..
301: * .. Executable Statements ..
302: *
303: * Decode and test the input parameters
304: *
305: BOTHV = LSAME( SIDE, 'B' )
306: RIGHTV = LSAME( SIDE, 'R' ) .OR. BOTHV
307: LEFTV = LSAME( SIDE, 'L' ) .OR. BOTHV
308: *
309: ALLV = LSAME( HOWMNY, 'A' )
310: OVER = LSAME( HOWMNY, 'B' )
311: SOMEV = LSAME( HOWMNY, 'S' )
312: *
313: * Set M to the number of columns required to store the selected
314: * eigenvectors.
315: *
316: IF( SOMEV ) THEN
317: M = 0
318: DO 10 J = 1, N
319: IF( SELECT( J ) )
320: $ M = M + 1
321: 10 CONTINUE
322: ELSE
323: M = N
324: END IF
325: *
326: INFO = 0
327: NB = ILAENV( 1, 'ZTREVC', SIDE // HOWMNY, N, -1, -1, -1 )
328: MAXWRK = N + 2*N*NB
329: WORK(1) = MAXWRK
330: RWORK(1) = N
331: LQUERY = ( LWORK.EQ.-1 .OR. LRWORK.EQ.-1 )
332: IF( .NOT.RIGHTV .AND. .NOT.LEFTV ) THEN
333: INFO = -1
334: ELSE IF( .NOT.ALLV .AND. .NOT.OVER .AND. .NOT.SOMEV ) THEN
335: INFO = -2
336: ELSE IF( N.LT.0 ) THEN
337: INFO = -4
338: ELSE IF( LDT.LT.MAX( 1, N ) ) THEN
339: INFO = -6
340: ELSE IF( LDVL.LT.1 .OR. ( LEFTV .AND. LDVL.LT.N ) ) THEN
341: INFO = -8
342: ELSE IF( LDVR.LT.1 .OR. ( RIGHTV .AND. LDVR.LT.N ) ) THEN
343: INFO = -10
344: ELSE IF( MM.LT.M ) THEN
345: INFO = -11
346: ELSE IF( LWORK.LT.MAX( 1, 2*N ) .AND. .NOT.LQUERY ) THEN
347: INFO = -14
348: ELSE IF ( LRWORK.LT.MAX( 1, N ) .AND. .NOT.LQUERY ) THEN
349: INFO = -16
350: END IF
351: IF( INFO.NE.0 ) THEN
352: CALL XERBLA( 'ZTREVC3', -INFO )
353: RETURN
354: ELSE IF( LQUERY ) THEN
355: RETURN
356: END IF
357: *
358: * Quick return if possible.
359: *
360: IF( N.EQ.0 )
361: $ RETURN
362: *
363: * Use blocked version of back-transformation if sufficient workspace.
364: * Zero-out the workspace to avoid potential NaN propagation.
365: *
366: IF( OVER .AND. LWORK .GE. N + 2*N*NBMIN ) THEN
367: NB = (LWORK - N) / (2*N)
368: NB = MIN( NB, NBMAX )
369: CALL ZLASET( 'F', N, 1+2*NB, CZERO, CZERO, WORK, N )
370: ELSE
371: NB = 1
372: END IF
373: *
374: * Set the constants to control overflow.
375: *
376: UNFL = DLAMCH( 'Safe minimum' )
377: OVFL = ONE / UNFL
378: CALL DLABAD( UNFL, OVFL )
379: ULP = DLAMCH( 'Precision' )
380: SMLNUM = UNFL*( N / ULP )
381: *
382: * Store the diagonal elements of T in working array WORK.
383: *
384: DO 20 I = 1, N
385: WORK( I ) = T( I, I )
386: 20 CONTINUE
387: *
388: * Compute 1-norm of each column of strictly upper triangular
389: * part of T to control overflow in triangular solver.
390: *
391: RWORK( 1 ) = ZERO
392: DO 30 J = 2, N
393: RWORK( J ) = DZASUM( J-1, T( 1, J ), 1 )
394: 30 CONTINUE
395: *
396: IF( RIGHTV ) THEN
397: *
398: * ============================================================
399: * Compute right eigenvectors.
400: *
401: * IV is index of column in current block.
402: * Non-blocked version always uses IV=NB=1;
403: * blocked version starts with IV=NB, goes down to 1.
404: * (Note the "0-th" column is used to store the original diagonal.)
405: IV = NB
406: IS = M
407: DO 80 KI = N, 1, -1
408: IF( SOMEV ) THEN
409: IF( .NOT.SELECT( KI ) )
410: $ GO TO 80
411: END IF
412: SMIN = MAX( ULP*( CABS1( T( KI, KI ) ) ), SMLNUM )
413: *
414: * --------------------------------------------------------
415: * Complex right eigenvector
416: *
417: WORK( KI + IV*N ) = CONE
418: *
419: * Form right-hand side.
420: *
421: DO 40 K = 1, KI - 1
422: WORK( K + IV*N ) = -T( K, KI )
423: 40 CONTINUE
424: *
425: * Solve upper triangular system:
426: * [ T(1:KI-1,1:KI-1) - T(KI,KI) ]*X = SCALE*WORK.
427: *
428: DO 50 K = 1, KI - 1
429: T( K, K ) = T( K, K ) - T( KI, KI )
430: IF( CABS1( T( K, K ) ).LT.SMIN )
431: $ T( K, K ) = SMIN
432: 50 CONTINUE
433: *
434: IF( KI.GT.1 ) THEN
435: CALL ZLATRS( 'Upper', 'No transpose', 'Non-unit', 'Y',
436: $ KI-1, T, LDT, WORK( 1 + IV*N ), SCALE,
437: $ RWORK, INFO )
438: WORK( KI + IV*N ) = SCALE
439: END IF
440: *
441: * Copy the vector x or Q*x to VR and normalize.
442: *
443: IF( .NOT.OVER ) THEN
444: * ------------------------------
445: * no back-transform: copy x to VR and normalize.
446: CALL ZCOPY( KI, WORK( 1 + IV*N ), 1, VR( 1, IS ), 1 )
447: *
448: II = IZAMAX( KI, VR( 1, IS ), 1 )
449: REMAX = ONE / CABS1( VR( II, IS ) )
450: CALL ZDSCAL( KI, REMAX, VR( 1, IS ), 1 )
451: *
452: DO 60 K = KI + 1, N
453: VR( K, IS ) = CZERO
454: 60 CONTINUE
455: *
456: ELSE IF( NB.EQ.1 ) THEN
457: * ------------------------------
458: * version 1: back-transform each vector with GEMV, Q*x.
459: IF( KI.GT.1 )
460: $ CALL ZGEMV( 'N', N, KI-1, CONE, VR, LDVR,
461: $ WORK( 1 + IV*N ), 1, DCMPLX( SCALE ),
462: $ VR( 1, KI ), 1 )
463: *
464: II = IZAMAX( N, VR( 1, KI ), 1 )
465: REMAX = ONE / CABS1( VR( II, KI ) )
466: CALL ZDSCAL( N, REMAX, VR( 1, KI ), 1 )
467: *
468: ELSE
469: * ------------------------------
470: * version 2: back-transform block of vectors with GEMM
471: * zero out below vector
472: DO K = KI + 1, N
473: WORK( K + IV*N ) = CZERO
474: END DO
475: *
476: * Columns IV:NB of work are valid vectors.
477: * When the number of vectors stored reaches NB,
478: * or if this was last vector, do the GEMM
479: IF( (IV.EQ.1) .OR. (KI.EQ.1) ) THEN
480: CALL ZGEMM( 'N', 'N', N, NB-IV+1, KI+NB-IV, CONE,
481: $ VR, LDVR,
482: $ WORK( 1 + (IV)*N ), N,
483: $ CZERO,
484: $ WORK( 1 + (NB+IV)*N ), N )
485: * normalize vectors
486: DO K = IV, NB
487: II = IZAMAX( N, WORK( 1 + (NB+K)*N ), 1 )
488: REMAX = ONE / CABS1( WORK( II + (NB+K)*N ) )
489: CALL ZDSCAL( N, REMAX, WORK( 1 + (NB+K)*N ), 1 )
490: END DO
491: CALL ZLACPY( 'F', N, NB-IV+1,
492: $ WORK( 1 + (NB+IV)*N ), N,
493: $ VR( 1, KI ), LDVR )
494: IV = NB
495: ELSE
496: IV = IV - 1
497: END IF
498: END IF
499: *
500: * Restore the original diagonal elements of T.
501: *
502: DO 70 K = 1, KI - 1
503: T( K, K ) = WORK( K )
504: 70 CONTINUE
505: *
506: IS = IS - 1
507: 80 CONTINUE
508: END IF
509: *
510: IF( LEFTV ) THEN
511: *
512: * ============================================================
513: * Compute left eigenvectors.
514: *
515: * IV is index of column in current block.
516: * Non-blocked version always uses IV=1;
517: * blocked version starts with IV=1, goes up to NB.
518: * (Note the "0-th" column is used to store the original diagonal.)
519: IV = 1
520: IS = 1
521: DO 130 KI = 1, N
522: *
523: IF( SOMEV ) THEN
524: IF( .NOT.SELECT( KI ) )
525: $ GO TO 130
526: END IF
527: SMIN = MAX( ULP*( CABS1( T( KI, KI ) ) ), SMLNUM )
528: *
529: * --------------------------------------------------------
530: * Complex left eigenvector
531: *
532: WORK( KI + IV*N ) = CONE
533: *
534: * Form right-hand side.
535: *
536: DO 90 K = KI + 1, N
537: WORK( K + IV*N ) = -CONJG( T( KI, K ) )
538: 90 CONTINUE
539: *
540: * Solve conjugate-transposed triangular system:
541: * [ T(KI+1:N,KI+1:N) - T(KI,KI) ]**H * X = SCALE*WORK.
542: *
543: DO 100 K = KI + 1, N
544: T( K, K ) = T( K, K ) - T( KI, KI )
545: IF( CABS1( T( K, K ) ).LT.SMIN )
546: $ T( K, K ) = SMIN
547: 100 CONTINUE
548: *
549: IF( KI.LT.N ) THEN
550: CALL ZLATRS( 'Upper', 'Conjugate transpose', 'Non-unit',
551: $ 'Y', N-KI, T( KI+1, KI+1 ), LDT,
552: $ WORK( KI+1 + IV*N ), SCALE, RWORK, INFO )
553: WORK( KI + IV*N ) = SCALE
554: END IF
555: *
556: * Copy the vector x or Q*x to VL and normalize.
557: *
558: IF( .NOT.OVER ) THEN
559: * ------------------------------
560: * no back-transform: copy x to VL and normalize.
561: CALL ZCOPY( N-KI+1, WORK( KI + IV*N ), 1, VL(KI,IS), 1 )
562: *
563: II = IZAMAX( N-KI+1, VL( KI, IS ), 1 ) + KI - 1
564: REMAX = ONE / CABS1( VL( II, IS ) )
565: CALL ZDSCAL( N-KI+1, REMAX, VL( KI, IS ), 1 )
566: *
567: DO 110 K = 1, KI - 1
568: VL( K, IS ) = CZERO
569: 110 CONTINUE
570: *
571: ELSE IF( NB.EQ.1 ) THEN
572: * ------------------------------
573: * version 1: back-transform each vector with GEMV, Q*x.
574: IF( KI.LT.N )
575: $ CALL ZGEMV( 'N', N, N-KI, CONE, VL( 1, KI+1 ), LDVL,
576: $ WORK( KI+1 + IV*N ), 1, DCMPLX( SCALE ),
577: $ VL( 1, KI ), 1 )
578: *
579: II = IZAMAX( N, VL( 1, KI ), 1 )
580: REMAX = ONE / CABS1( VL( II, KI ) )
581: CALL ZDSCAL( N, REMAX, VL( 1, KI ), 1 )
582: *
583: ELSE
584: * ------------------------------
585: * version 2: back-transform block of vectors with GEMM
586: * zero out above vector
587: * could go from KI-NV+1 to KI-1
588: DO K = 1, KI - 1
589: WORK( K + IV*N ) = CZERO
590: END DO
591: *
592: * Columns 1:IV of work are valid vectors.
593: * When the number of vectors stored reaches NB,
594: * or if this was last vector, do the GEMM
595: IF( (IV.EQ.NB) .OR. (KI.EQ.N) ) THEN
596: CALL ZGEMM( 'N', 'N', N, IV, N-KI+IV, ONE,
597: $ VL( 1, KI-IV+1 ), LDVL,
598: $ WORK( KI-IV+1 + (1)*N ), N,
599: $ CZERO,
600: $ WORK( 1 + (NB+1)*N ), N )
601: * normalize vectors
602: DO K = 1, IV
603: II = IZAMAX( N, WORK( 1 + (NB+K)*N ), 1 )
604: REMAX = ONE / CABS1( WORK( II + (NB+K)*N ) )
605: CALL ZDSCAL( N, REMAX, WORK( 1 + (NB+K)*N ), 1 )
606: END DO
607: CALL ZLACPY( 'F', N, IV,
608: $ WORK( 1 + (NB+1)*N ), N,
609: $ VL( 1, KI-IV+1 ), LDVL )
610: IV = 1
611: ELSE
612: IV = IV + 1
613: END IF
614: END IF
615: *
616: * Restore the original diagonal elements of T.
617: *
618: DO 120 K = KI + 1, N
619: T( K, K ) = WORK( K )
620: 120 CONTINUE
621: *
622: IS = IS + 1
623: 130 CONTINUE
624: END IF
625: *
626: RETURN
627: *
628: * End of ZTREVC3
629: *
630: END
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to ztrevc3.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack