Annotation of rpl/lapack/lapack/ztrevc3.f, revision 1.7
1.1 bertrand 1: *> \brief \b ZTREVC3
2: *
3: * =========== DOCUMENTATION ===========
4: *
1.3 bertrand 5: * Online html documentation available at
6: * http://www.netlib.org/lapack/explore-html/
1.1 bertrand 7: *
8: *> \htmlonly
1.3 bertrand 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">
1.1 bertrand 15: *> [TXT]</a>
1.3 bertrand 16: *> \endhtmlonly
1.1 bertrand 17: *
18: * Definition:
19: * ===========
20: *
1.3 bertrand 21: * SUBROUTINE ZTREVC3( SIDE, HOWMNY, SELECT, N, T, LDT, VL, LDVL, VR,
22: * $ LDVR, MM, M, WORK, LWORK, RWORK, LRWORK, INFO)
1.1 bertrand 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: *> \ingroup complex16OTHERcomputational
226: *
227: *> \par Further Details:
228: * =====================
229: *>
230: *> \verbatim
231: *>
232: *> The algorithm used in this program is basically backward (forward)
233: *> substitution, with scaling to make the the code robust against
234: *> possible overflow.
235: *>
236: *> Each eigenvector is normalized so that the element of largest
237: *> magnitude has magnitude 1; here the magnitude of a complex number
238: *> (x,y) is taken to be |x| + |y|.
239: *> \endverbatim
240: *>
241: * =====================================================================
242: SUBROUTINE ZTREVC3( SIDE, HOWMNY, SELECT, N, T, LDT, VL, LDVL, VR,
243: $ LDVR, MM, M, WORK, LWORK, RWORK, LRWORK, INFO)
244: IMPLICIT NONE
245: *
1.7 ! bertrand 246: * -- LAPACK computational routine --
1.1 bertrand 247: * -- LAPACK is a software package provided by Univ. of Tennessee, --
248: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
249: *
250: * .. Scalar Arguments ..
251: CHARACTER HOWMNY, SIDE
252: INTEGER INFO, LDT, LDVL, LDVR, LWORK, LRWORK, M, MM, N
253: * ..
254: * .. Array Arguments ..
255: LOGICAL SELECT( * )
256: DOUBLE PRECISION RWORK( * )
257: COMPLEX*16 T( LDT, * ), VL( LDVL, * ), VR( LDVR, * ),
258: $ WORK( * )
259: * ..
260: *
261: * =====================================================================
262: *
263: * .. Parameters ..
264: DOUBLE PRECISION ZERO, ONE
265: PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 )
266: COMPLEX*16 CZERO, CONE
267: PARAMETER ( CZERO = ( 0.0D+0, 0.0D+0 ),
268: $ CONE = ( 1.0D+0, 0.0D+0 ) )
269: INTEGER NBMIN, NBMAX
270: PARAMETER ( NBMIN = 8, NBMAX = 128 )
271: * ..
272: * .. Local Scalars ..
273: LOGICAL ALLV, BOTHV, LEFTV, LQUERY, OVER, RIGHTV, SOMEV
274: INTEGER I, II, IS, J, K, KI, IV, MAXWRK, NB
275: DOUBLE PRECISION OVFL, REMAX, SCALE, SMIN, SMLNUM, ULP, UNFL
276: COMPLEX*16 CDUM
277: * ..
278: * .. External Functions ..
279: LOGICAL LSAME
280: INTEGER ILAENV, IZAMAX
281: DOUBLE PRECISION DLAMCH, DZASUM
282: EXTERNAL LSAME, ILAENV, IZAMAX, DLAMCH, DZASUM
283: * ..
284: * .. External Subroutines ..
1.5 bertrand 285: EXTERNAL XERBLA, ZCOPY, ZDSCAL, ZGEMV, ZLATRS,
286: $ ZGEMM, DLABAD, ZLASET, ZLACPY
1.1 bertrand 287: * ..
288: * .. Intrinsic Functions ..
1.7 ! bertrand 289: INTRINSIC ABS, DBLE, DCMPLX, CONJG, DIMAG, MAX
1.1 bertrand 290: * ..
291: * .. Statement Functions ..
292: DOUBLE PRECISION CABS1
293: * ..
294: * .. Statement Function definitions ..
1.7 ! bertrand 295: CABS1( CDUM ) = ABS( DBLE( CDUM ) ) + ABS( DIMAG( CDUM ) )
1.1 bertrand 296: * ..
297: * .. Executable Statements ..
298: *
299: * Decode and test the input parameters
300: *
301: BOTHV = LSAME( SIDE, 'B' )
302: RIGHTV = LSAME( SIDE, 'R' ) .OR. BOTHV
303: LEFTV = LSAME( SIDE, 'L' ) .OR. BOTHV
304: *
305: ALLV = LSAME( HOWMNY, 'A' )
306: OVER = LSAME( HOWMNY, 'B' )
307: SOMEV = LSAME( HOWMNY, 'S' )
308: *
309: * Set M to the number of columns required to store the selected
310: * eigenvectors.
311: *
312: IF( SOMEV ) THEN
313: M = 0
314: DO 10 J = 1, N
315: IF( SELECT( J ) )
316: $ M = M + 1
317: 10 CONTINUE
318: ELSE
319: M = N
320: END IF
321: *
322: INFO = 0
323: NB = ILAENV( 1, 'ZTREVC', SIDE // HOWMNY, N, -1, -1, -1 )
324: MAXWRK = N + 2*N*NB
325: WORK(1) = MAXWRK
326: RWORK(1) = N
327: LQUERY = ( LWORK.EQ.-1 .OR. LRWORK.EQ.-1 )
328: IF( .NOT.RIGHTV .AND. .NOT.LEFTV ) THEN
329: INFO = -1
330: ELSE IF( .NOT.ALLV .AND. .NOT.OVER .AND. .NOT.SOMEV ) THEN
331: INFO = -2
332: ELSE IF( N.LT.0 ) THEN
333: INFO = -4
334: ELSE IF( LDT.LT.MAX( 1, N ) ) THEN
335: INFO = -6
336: ELSE IF( LDVL.LT.1 .OR. ( LEFTV .AND. LDVL.LT.N ) ) THEN
337: INFO = -8
338: ELSE IF( LDVR.LT.1 .OR. ( RIGHTV .AND. LDVR.LT.N ) ) THEN
339: INFO = -10
340: ELSE IF( MM.LT.M ) THEN
341: INFO = -11
342: ELSE IF( LWORK.LT.MAX( 1, 2*N ) .AND. .NOT.LQUERY ) THEN
343: INFO = -14
344: ELSE IF ( LRWORK.LT.MAX( 1, N ) .AND. .NOT.LQUERY ) THEN
345: INFO = -16
346: END IF
347: IF( INFO.NE.0 ) THEN
348: CALL XERBLA( 'ZTREVC3', -INFO )
349: RETURN
350: ELSE IF( LQUERY ) THEN
351: RETURN
352: END IF
353: *
354: * Quick return if possible.
355: *
356: IF( N.EQ.0 )
357: $ RETURN
358: *
359: * Use blocked version of back-transformation if sufficient workspace.
360: * Zero-out the workspace to avoid potential NaN propagation.
361: *
362: IF( OVER .AND. LWORK .GE. N + 2*N*NBMIN ) THEN
363: NB = (LWORK - N) / (2*N)
364: NB = MIN( NB, NBMAX )
365: CALL ZLASET( 'F', N, 1+2*NB, CZERO, CZERO, WORK, N )
366: ELSE
367: NB = 1
368: END IF
369: *
370: * Set the constants to control overflow.
371: *
372: UNFL = DLAMCH( 'Safe minimum' )
373: OVFL = ONE / UNFL
374: CALL DLABAD( UNFL, OVFL )
375: ULP = DLAMCH( 'Precision' )
376: SMLNUM = UNFL*( N / ULP )
377: *
378: * Store the diagonal elements of T in working array WORK.
379: *
380: DO 20 I = 1, N
381: WORK( I ) = T( I, I )
382: 20 CONTINUE
383: *
384: * Compute 1-norm of each column of strictly upper triangular
385: * part of T to control overflow in triangular solver.
386: *
387: RWORK( 1 ) = ZERO
388: DO 30 J = 2, N
389: RWORK( J ) = DZASUM( J-1, T( 1, J ), 1 )
390: 30 CONTINUE
391: *
392: IF( RIGHTV ) THEN
393: *
394: * ============================================================
395: * Compute right eigenvectors.
396: *
397: * IV is index of column in current block.
398: * Non-blocked version always uses IV=NB=1;
399: * blocked version starts with IV=NB, goes down to 1.
400: * (Note the "0-th" column is used to store the original diagonal.)
401: IV = NB
402: IS = M
403: DO 80 KI = N, 1, -1
404: IF( SOMEV ) THEN
405: IF( .NOT.SELECT( KI ) )
406: $ GO TO 80
407: END IF
408: SMIN = MAX( ULP*( CABS1( T( KI, KI ) ) ), SMLNUM )
409: *
410: * --------------------------------------------------------
411: * Complex right eigenvector
412: *
413: WORK( KI + IV*N ) = CONE
414: *
415: * Form right-hand side.
416: *
417: DO 40 K = 1, KI - 1
418: WORK( K + IV*N ) = -T( K, KI )
419: 40 CONTINUE
420: *
421: * Solve upper triangular system:
422: * [ T(1:KI-1,1:KI-1) - T(KI,KI) ]*X = SCALE*WORK.
423: *
424: DO 50 K = 1, KI - 1
425: T( K, K ) = T( K, K ) - T( KI, KI )
426: IF( CABS1( T( K, K ) ).LT.SMIN )
427: $ T( K, K ) = SMIN
428: 50 CONTINUE
429: *
430: IF( KI.GT.1 ) THEN
431: CALL ZLATRS( 'Upper', 'No transpose', 'Non-unit', 'Y',
432: $ KI-1, T, LDT, WORK( 1 + IV*N ), SCALE,
433: $ RWORK, INFO )
434: WORK( KI + IV*N ) = SCALE
435: END IF
436: *
437: * Copy the vector x or Q*x to VR and normalize.
438: *
439: IF( .NOT.OVER ) THEN
440: * ------------------------------
441: * no back-transform: copy x to VR and normalize.
442: CALL ZCOPY( KI, WORK( 1 + IV*N ), 1, VR( 1, IS ), 1 )
443: *
444: II = IZAMAX( KI, VR( 1, IS ), 1 )
445: REMAX = ONE / CABS1( VR( II, IS ) )
446: CALL ZDSCAL( KI, REMAX, VR( 1, IS ), 1 )
447: *
448: DO 60 K = KI + 1, N
449: VR( K, IS ) = CZERO
450: 60 CONTINUE
451: *
452: ELSE IF( NB.EQ.1 ) THEN
453: * ------------------------------
454: * version 1: back-transform each vector with GEMV, Q*x.
455: IF( KI.GT.1 )
456: $ CALL ZGEMV( 'N', N, KI-1, CONE, VR, LDVR,
457: $ WORK( 1 + IV*N ), 1, DCMPLX( SCALE ),
458: $ VR( 1, KI ), 1 )
459: *
460: II = IZAMAX( N, VR( 1, KI ), 1 )
461: REMAX = ONE / CABS1( VR( II, KI ) )
462: CALL ZDSCAL( N, REMAX, VR( 1, KI ), 1 )
463: *
464: ELSE
465: * ------------------------------
466: * version 2: back-transform block of vectors with GEMM
467: * zero out below vector
468: DO K = KI + 1, N
469: WORK( K + IV*N ) = CZERO
470: END DO
471: *
472: * Columns IV:NB of work are valid vectors.
473: * When the number of vectors stored reaches NB,
474: * or if this was last vector, do the GEMM
475: IF( (IV.EQ.1) .OR. (KI.EQ.1) ) THEN
476: CALL ZGEMM( 'N', 'N', N, NB-IV+1, KI+NB-IV, CONE,
477: $ VR, LDVR,
478: $ WORK( 1 + (IV)*N ), N,
479: $ CZERO,
480: $ WORK( 1 + (NB+IV)*N ), N )
481: * normalize vectors
482: DO K = IV, NB
483: II = IZAMAX( N, WORK( 1 + (NB+K)*N ), 1 )
484: REMAX = ONE / CABS1( WORK( II + (NB+K)*N ) )
485: CALL ZDSCAL( N, REMAX, WORK( 1 + (NB+K)*N ), 1 )
486: END DO
487: CALL ZLACPY( 'F', N, NB-IV+1,
488: $ WORK( 1 + (NB+IV)*N ), N,
489: $ VR( 1, KI ), LDVR )
490: IV = NB
491: ELSE
492: IV = IV - 1
493: END IF
494: END IF
495: *
496: * Restore the original diagonal elements of T.
497: *
498: DO 70 K = 1, KI - 1
499: T( K, K ) = WORK( K )
500: 70 CONTINUE
501: *
502: IS = IS - 1
503: 80 CONTINUE
504: END IF
505: *
506: IF( LEFTV ) THEN
507: *
508: * ============================================================
509: * Compute left eigenvectors.
510: *
511: * IV is index of column in current block.
512: * Non-blocked version always uses IV=1;
513: * blocked version starts with IV=1, goes up to NB.
514: * (Note the "0-th" column is used to store the original diagonal.)
515: IV = 1
516: IS = 1
517: DO 130 KI = 1, N
518: *
519: IF( SOMEV ) THEN
520: IF( .NOT.SELECT( KI ) )
521: $ GO TO 130
522: END IF
523: SMIN = MAX( ULP*( CABS1( T( KI, KI ) ) ), SMLNUM )
524: *
525: * --------------------------------------------------------
526: * Complex left eigenvector
527: *
528: WORK( KI + IV*N ) = CONE
529: *
530: * Form right-hand side.
531: *
532: DO 90 K = KI + 1, N
533: WORK( K + IV*N ) = -CONJG( T( KI, K ) )
534: 90 CONTINUE
535: *
536: * Solve conjugate-transposed triangular system:
537: * [ T(KI+1:N,KI+1:N) - T(KI,KI) ]**H * X = SCALE*WORK.
538: *
539: DO 100 K = KI + 1, N
540: T( K, K ) = T( K, K ) - T( KI, KI )
541: IF( CABS1( T( K, K ) ).LT.SMIN )
542: $ T( K, K ) = SMIN
543: 100 CONTINUE
544: *
545: IF( KI.LT.N ) THEN
546: CALL ZLATRS( 'Upper', 'Conjugate transpose', 'Non-unit',
547: $ 'Y', N-KI, T( KI+1, KI+1 ), LDT,
548: $ WORK( KI+1 + IV*N ), SCALE, RWORK, INFO )
549: WORK( KI + IV*N ) = SCALE
550: END IF
551: *
552: * Copy the vector x or Q*x to VL and normalize.
553: *
554: IF( .NOT.OVER ) THEN
555: * ------------------------------
556: * no back-transform: copy x to VL and normalize.
557: CALL ZCOPY( N-KI+1, WORK( KI + IV*N ), 1, VL(KI,IS), 1 )
558: *
559: II = IZAMAX( N-KI+1, VL( KI, IS ), 1 ) + KI - 1
560: REMAX = ONE / CABS1( VL( II, IS ) )
561: CALL ZDSCAL( N-KI+1, REMAX, VL( KI, IS ), 1 )
562: *
563: DO 110 K = 1, KI - 1
564: VL( K, IS ) = CZERO
565: 110 CONTINUE
566: *
567: ELSE IF( NB.EQ.1 ) THEN
568: * ------------------------------
569: * version 1: back-transform each vector with GEMV, Q*x.
570: IF( KI.LT.N )
571: $ CALL ZGEMV( 'N', N, N-KI, CONE, VL( 1, KI+1 ), LDVL,
572: $ WORK( KI+1 + IV*N ), 1, DCMPLX( SCALE ),
573: $ VL( 1, KI ), 1 )
574: *
575: II = IZAMAX( N, VL( 1, KI ), 1 )
576: REMAX = ONE / CABS1( VL( II, KI ) )
577: CALL ZDSCAL( N, REMAX, VL( 1, KI ), 1 )
578: *
579: ELSE
580: * ------------------------------
581: * version 2: back-transform block of vectors with GEMM
582: * zero out above vector
583: * could go from KI-NV+1 to KI-1
584: DO K = 1, KI - 1
585: WORK( K + IV*N ) = CZERO
586: END DO
587: *
588: * Columns 1:IV of work are valid vectors.
589: * When the number of vectors stored reaches NB,
590: * or if this was last vector, do the GEMM
591: IF( (IV.EQ.NB) .OR. (KI.EQ.N) ) THEN
1.3 bertrand 592: CALL ZGEMM( 'N', 'N', N, IV, N-KI+IV, CONE,
1.1 bertrand 593: $ VL( 1, KI-IV+1 ), LDVL,
594: $ WORK( KI-IV+1 + (1)*N ), N,
595: $ CZERO,
596: $ WORK( 1 + (NB+1)*N ), N )
597: * normalize vectors
598: DO K = 1, IV
599: II = IZAMAX( N, WORK( 1 + (NB+K)*N ), 1 )
600: REMAX = ONE / CABS1( WORK( II + (NB+K)*N ) )
601: CALL ZDSCAL( N, REMAX, WORK( 1 + (NB+K)*N ), 1 )
602: END DO
603: CALL ZLACPY( 'F', N, IV,
604: $ WORK( 1 + (NB+1)*N ), N,
605: $ VL( 1, KI-IV+1 ), LDVL )
606: IV = 1
607: ELSE
608: IV = IV + 1
609: END IF
610: END IF
611: *
612: * Restore the original diagonal elements of T.
613: *
614: DO 120 K = KI + 1, N
615: T( K, K ) = WORK( K )
616: 120 CONTINUE
617: *
618: IS = IS + 1
619: 130 CONTINUE
620: END IF
621: *
622: RETURN
623: *
624: * End of ZTREVC3
625: *
626: 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