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, VR,
22: * $ LDVR, MM, M, WORK, LWORK, RWORK, LRWORK, 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 December 2016
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.7.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: * December 2016
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: $ ZGEMM, DLABAD, ZLASET
292: * ..
293: * .. Intrinsic Functions ..
294: INTRINSIC ABS, DBLE, DCMPLX, CONJG, AIMAG, MAX
295: * ..
296: * .. Statement Functions ..
297: DOUBLE PRECISION CABS1
298: * ..
299: * .. Statement Function definitions ..
300: CABS1( CDUM ) = ABS( DBLE( CDUM ) ) + ABS( AIMAG( CDUM ) )
301: * ..
302: * .. Executable Statements ..
303: *
304: * Decode and test the input parameters
305: *
306: BOTHV = LSAME( SIDE, 'B' )
307: RIGHTV = LSAME( SIDE, 'R' ) .OR. BOTHV
308: LEFTV = LSAME( SIDE, 'L' ) .OR. BOTHV
309: *
310: ALLV = LSAME( HOWMNY, 'A' )
311: OVER = LSAME( HOWMNY, 'B' )
312: SOMEV = LSAME( HOWMNY, 'S' )
313: *
314: * Set M to the number of columns required to store the selected
315: * eigenvectors.
316: *
317: IF( SOMEV ) THEN
318: M = 0
319: DO 10 J = 1, N
320: IF( SELECT( J ) )
321: $ M = M + 1
322: 10 CONTINUE
323: ELSE
324: M = N
325: END IF
326: *
327: INFO = 0
328: NB = ILAENV( 1, 'ZTREVC', SIDE // HOWMNY, N, -1, -1, -1 )
329: MAXWRK = N + 2*N*NB
330: WORK(1) = MAXWRK
331: RWORK(1) = N
332: LQUERY = ( LWORK.EQ.-1 .OR. LRWORK.EQ.-1 )
333: IF( .NOT.RIGHTV .AND. .NOT.LEFTV ) THEN
334: INFO = -1
335: ELSE IF( .NOT.ALLV .AND. .NOT.OVER .AND. .NOT.SOMEV ) THEN
336: INFO = -2
337: ELSE IF( N.LT.0 ) THEN
338: INFO = -4
339: ELSE IF( LDT.LT.MAX( 1, N ) ) THEN
340: INFO = -6
341: ELSE IF( LDVL.LT.1 .OR. ( LEFTV .AND. LDVL.LT.N ) ) THEN
342: INFO = -8
343: ELSE IF( LDVR.LT.1 .OR. ( RIGHTV .AND. LDVR.LT.N ) ) THEN
344: INFO = -10
345: ELSE IF( MM.LT.M ) THEN
346: INFO = -11
347: ELSE IF( LWORK.LT.MAX( 1, 2*N ) .AND. .NOT.LQUERY ) THEN
348: INFO = -14
349: ELSE IF ( LRWORK.LT.MAX( 1, N ) .AND. .NOT.LQUERY ) THEN
350: INFO = -16
351: END IF
352: IF( INFO.NE.0 ) THEN
353: CALL XERBLA( 'ZTREVC3', -INFO )
354: RETURN
355: ELSE IF( LQUERY ) THEN
356: RETURN
357: END IF
358: *
359: * Quick return if possible.
360: *
361: IF( N.EQ.0 )
362: $ RETURN
363: *
364: * Use blocked version of back-transformation if sufficient workspace.
365: * Zero-out the workspace to avoid potential NaN propagation.
366: *
367: IF( OVER .AND. LWORK .GE. N + 2*N*NBMIN ) THEN
368: NB = (LWORK - N) / (2*N)
369: NB = MIN( NB, NBMAX )
370: CALL ZLASET( 'F', N, 1+2*NB, CZERO, CZERO, WORK, N )
371: ELSE
372: NB = 1
373: END IF
374: *
375: * Set the constants to control overflow.
376: *
377: UNFL = DLAMCH( 'Safe minimum' )
378: OVFL = ONE / UNFL
379: CALL DLABAD( UNFL, OVFL )
380: ULP = DLAMCH( 'Precision' )
381: SMLNUM = UNFL*( N / ULP )
382: *
383: * Store the diagonal elements of T in working array WORK.
384: *
385: DO 20 I = 1, N
386: WORK( I ) = T( I, I )
387: 20 CONTINUE
388: *
389: * Compute 1-norm of each column of strictly upper triangular
390: * part of T to control overflow in triangular solver.
391: *
392: RWORK( 1 ) = ZERO
393: DO 30 J = 2, N
394: RWORK( J ) = DZASUM( J-1, T( 1, J ), 1 )
395: 30 CONTINUE
396: *
397: IF( RIGHTV ) THEN
398: *
399: * ============================================================
400: * Compute right eigenvectors.
401: *
402: * IV is index of column in current block.
403: * Non-blocked version always uses IV=NB=1;
404: * blocked version starts with IV=NB, goes down to 1.
405: * (Note the "0-th" column is used to store the original diagonal.)
406: IV = NB
407: IS = M
408: DO 80 KI = N, 1, -1
409: IF( SOMEV ) THEN
410: IF( .NOT.SELECT( KI ) )
411: $ GO TO 80
412: END IF
413: SMIN = MAX( ULP*( CABS1( T( KI, KI ) ) ), SMLNUM )
414: *
415: * --------------------------------------------------------
416: * Complex right eigenvector
417: *
418: WORK( KI + IV*N ) = CONE
419: *
420: * Form right-hand side.
421: *
422: DO 40 K = 1, KI - 1
423: WORK( K + IV*N ) = -T( K, KI )
424: 40 CONTINUE
425: *
426: * Solve upper triangular system:
427: * [ T(1:KI-1,1:KI-1) - T(KI,KI) ]*X = SCALE*WORK.
428: *
429: DO 50 K = 1, KI - 1
430: T( K, K ) = T( K, K ) - T( KI, KI )
431: IF( CABS1( T( K, K ) ).LT.SMIN )
432: $ T( K, K ) = SMIN
433: 50 CONTINUE
434: *
435: IF( KI.GT.1 ) THEN
436: CALL ZLATRS( 'Upper', 'No transpose', 'Non-unit', 'Y',
437: $ KI-1, T, LDT, WORK( 1 + IV*N ), SCALE,
438: $ RWORK, INFO )
439: WORK( KI + IV*N ) = SCALE
440: END IF
441: *
442: * Copy the vector x or Q*x to VR and normalize.
443: *
444: IF( .NOT.OVER ) THEN
445: * ------------------------------
446: * no back-transform: copy x to VR and normalize.
447: CALL ZCOPY( KI, WORK( 1 + IV*N ), 1, VR( 1, IS ), 1 )
448: *
449: II = IZAMAX( KI, VR( 1, IS ), 1 )
450: REMAX = ONE / CABS1( VR( II, IS ) )
451: CALL ZDSCAL( KI, REMAX, VR( 1, IS ), 1 )
452: *
453: DO 60 K = KI + 1, N
454: VR( K, IS ) = CZERO
455: 60 CONTINUE
456: *
457: ELSE IF( NB.EQ.1 ) THEN
458: * ------------------------------
459: * version 1: back-transform each vector with GEMV, Q*x.
460: IF( KI.GT.1 )
461: $ CALL ZGEMV( 'N', N, KI-1, CONE, VR, LDVR,
462: $ WORK( 1 + IV*N ), 1, DCMPLX( SCALE ),
463: $ VR( 1, KI ), 1 )
464: *
465: II = IZAMAX( N, VR( 1, KI ), 1 )
466: REMAX = ONE / CABS1( VR( II, KI ) )
467: CALL ZDSCAL( N, REMAX, VR( 1, KI ), 1 )
468: *
469: ELSE
470: * ------------------------------
471: * version 2: back-transform block of vectors with GEMM
472: * zero out below vector
473: DO K = KI + 1, N
474: WORK( K + IV*N ) = CZERO
475: END DO
476: *
477: * Columns IV:NB of work are valid vectors.
478: * When the number of vectors stored reaches NB,
479: * or if this was last vector, do the GEMM
480: IF( (IV.EQ.1) .OR. (KI.EQ.1) ) THEN
481: CALL ZGEMM( 'N', 'N', N, NB-IV+1, KI+NB-IV, CONE,
482: $ VR, LDVR,
483: $ WORK( 1 + (IV)*N ), N,
484: $ CZERO,
485: $ WORK( 1 + (NB+IV)*N ), N )
486: * normalize vectors
487: DO K = IV, NB
488: II = IZAMAX( N, WORK( 1 + (NB+K)*N ), 1 )
489: REMAX = ONE / CABS1( WORK( II + (NB+K)*N ) )
490: CALL ZDSCAL( N, REMAX, WORK( 1 + (NB+K)*N ), 1 )
491: END DO
492: CALL ZLACPY( 'F', N, NB-IV+1,
493: $ WORK( 1 + (NB+IV)*N ), N,
494: $ VR( 1, KI ), LDVR )
495: IV = NB
496: ELSE
497: IV = IV - 1
498: END IF
499: END IF
500: *
501: * Restore the original diagonal elements of T.
502: *
503: DO 70 K = 1, KI - 1
504: T( K, K ) = WORK( K )
505: 70 CONTINUE
506: *
507: IS = IS - 1
508: 80 CONTINUE
509: END IF
510: *
511: IF( LEFTV ) THEN
512: *
513: * ============================================================
514: * Compute left eigenvectors.
515: *
516: * IV is index of column in current block.
517: * Non-blocked version always uses IV=1;
518: * blocked version starts with IV=1, goes up to NB.
519: * (Note the "0-th" column is used to store the original diagonal.)
520: IV = 1
521: IS = 1
522: DO 130 KI = 1, N
523: *
524: IF( SOMEV ) THEN
525: IF( .NOT.SELECT( KI ) )
526: $ GO TO 130
527: END IF
528: SMIN = MAX( ULP*( CABS1( T( KI, KI ) ) ), SMLNUM )
529: *
530: * --------------------------------------------------------
531: * Complex left eigenvector
532: *
533: WORK( KI + IV*N ) = CONE
534: *
535: * Form right-hand side.
536: *
537: DO 90 K = KI + 1, N
538: WORK( K + IV*N ) = -CONJG( T( KI, K ) )
539: 90 CONTINUE
540: *
541: * Solve conjugate-transposed triangular system:
542: * [ T(KI+1:N,KI+1:N) - T(KI,KI) ]**H * X = SCALE*WORK.
543: *
544: DO 100 K = KI + 1, N
545: T( K, K ) = T( K, K ) - T( KI, KI )
546: IF( CABS1( T( K, K ) ).LT.SMIN )
547: $ T( K, K ) = SMIN
548: 100 CONTINUE
549: *
550: IF( KI.LT.N ) THEN
551: CALL ZLATRS( 'Upper', 'Conjugate transpose', 'Non-unit',
552: $ 'Y', N-KI, T( KI+1, KI+1 ), LDT,
553: $ WORK( KI+1 + IV*N ), SCALE, RWORK, INFO )
554: WORK( KI + IV*N ) = SCALE
555: END IF
556: *
557: * Copy the vector x or Q*x to VL and normalize.
558: *
559: IF( .NOT.OVER ) THEN
560: * ------------------------------
561: * no back-transform: copy x to VL and normalize.
562: CALL ZCOPY( N-KI+1, WORK( KI + IV*N ), 1, VL(KI,IS), 1 )
563: *
564: II = IZAMAX( N-KI+1, VL( KI, IS ), 1 ) + KI - 1
565: REMAX = ONE / CABS1( VL( II, IS ) )
566: CALL ZDSCAL( N-KI+1, REMAX, VL( KI, IS ), 1 )
567: *
568: DO 110 K = 1, KI - 1
569: VL( K, IS ) = CZERO
570: 110 CONTINUE
571: *
572: ELSE IF( NB.EQ.1 ) THEN
573: * ------------------------------
574: * version 1: back-transform each vector with GEMV, Q*x.
575: IF( KI.LT.N )
576: $ CALL ZGEMV( 'N', N, N-KI, CONE, VL( 1, KI+1 ), LDVL,
577: $ WORK( KI+1 + IV*N ), 1, DCMPLX( SCALE ),
578: $ VL( 1, KI ), 1 )
579: *
580: II = IZAMAX( N, VL( 1, KI ), 1 )
581: REMAX = ONE / CABS1( VL( II, KI ) )
582: CALL ZDSCAL( N, REMAX, VL( 1, KI ), 1 )
583: *
584: ELSE
585: * ------------------------------
586: * version 2: back-transform block of vectors with GEMM
587: * zero out above vector
588: * could go from KI-NV+1 to KI-1
589: DO K = 1, KI - 1
590: WORK( K + IV*N ) = CZERO
591: END DO
592: *
593: * Columns 1:IV of work are valid vectors.
594: * When the number of vectors stored reaches NB,
595: * or if this was last vector, do the GEMM
596: IF( (IV.EQ.NB) .OR. (KI.EQ.N) ) THEN
597: CALL ZGEMM( 'N', 'N', N, IV, N-KI+IV, CONE,
598: $ VL( 1, KI-IV+1 ), LDVL,
599: $ WORK( KI-IV+1 + (1)*N ), N,
600: $ CZERO,
601: $ WORK( 1 + (NB+1)*N ), N )
602: * normalize vectors
603: DO K = 1, IV
604: II = IZAMAX( N, WORK( 1 + (NB+K)*N ), 1 )
605: REMAX = ONE / CABS1( WORK( II + (NB+K)*N ) )
606: CALL ZDSCAL( N, REMAX, WORK( 1 + (NB+K)*N ), 1 )
607: END DO
608: CALL ZLACPY( 'F', N, IV,
609: $ WORK( 1 + (NB+1)*N ), N,
610: $ VL( 1, KI-IV+1 ), LDVL )
611: IV = 1
612: ELSE
613: IV = IV + 1
614: END IF
615: END IF
616: *
617: * Restore the original diagonal elements of T.
618: *
619: DO 120 K = KI + 1, N
620: T( K, K ) = WORK( K )
621: 120 CONTINUE
622: *
623: IS = IS + 1
624: 130 CONTINUE
625: END IF
626: *
627: RETURN
628: *
629: * End of ZTREVC3
630: *
631: 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