001:       SUBROUTINE CLATRD( UPLO, N, NB, A, LDA, E, TAU, W, LDW )
002: *
003: *  -- LAPACK auxiliary routine (version 3.2) --
004: *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
005: *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
006: *     November 2006
007: *
008: *     .. Scalar Arguments ..
009:       CHARACTER          UPLO
010:       INTEGER            LDA, LDW, N, NB
011: *     ..
012: *     .. Array Arguments ..
013:       REAL               E( * )
014:       COMPLEX            A( LDA, * ), TAU( * ), W( LDW, * )
015: *     ..
016: *
017: *  Purpose
018: *  =======
019: *
020: *  CLATRD reduces NB rows and columns of a complex Hermitian matrix A to
021: *  Hermitian tridiagonal form by a unitary similarity
022: *  transformation Q' * A * Q, and returns the matrices V and W which are
023: *  needed to apply the transformation to the unreduced part of A.
024: *
025: *  If UPLO = 'U', CLATRD reduces the last NB rows and columns of a
026: *  matrix, of which the upper triangle is supplied;
027: *  if UPLO = 'L', CLATRD reduces the first NB rows and columns of a
028: *  matrix, of which the lower triangle is supplied.
029: *
030: *  This is an auxiliary routine called by CHETRD.
031: *
032: *  Arguments
033: *  =========
034: *
035: *  UPLO    (input) CHARACTER*1
036: *          Specifies whether the upper or lower triangular part of the
037: *          Hermitian matrix A is stored:
038: *          = 'U': Upper triangular
039: *          = 'L': Lower triangular
040: *
041: *  N       (input) INTEGER
042: *          The order of the matrix A.
043: *
044: *  NB      (input) INTEGER
045: *          The number of rows and columns to be reduced.
046: *
047: *  A       (input/output) COMPLEX array, dimension (LDA,N)
048: *          On entry, the Hermitian matrix A.  If UPLO = 'U', the leading
049: *          n-by-n upper triangular part of A contains the upper
050: *          triangular part of the matrix A, and the strictly lower
051: *          triangular part of A is not referenced.  If UPLO = 'L', the
052: *          leading n-by-n lower triangular part of A contains the lower
053: *          triangular part of the matrix A, and the strictly upper
054: *          triangular part of A is not referenced.
055: *          On exit:
056: *          if UPLO = 'U', the last NB columns have been reduced to
057: *            tridiagonal form, with the diagonal elements overwriting
058: *            the diagonal elements of A; the elements above the diagonal
059: *            with the array TAU, represent the unitary matrix Q as a
060: *            product of elementary reflectors;
061: *          if UPLO = 'L', the first NB columns have been reduced to
062: *            tridiagonal form, with the diagonal elements overwriting
063: *            the diagonal elements of A; the elements below the diagonal
064: *            with the array TAU, represent the  unitary matrix Q as a
065: *            product of elementary reflectors.
066: *          See Further Details.
067: *
068: *  LDA     (input) INTEGER
069: *          The leading dimension of the array A.  LDA >= max(1,N).
070: *
071: *  E       (output) REAL array, dimension (N-1)
072: *          If UPLO = 'U', E(n-nb:n-1) contains the superdiagonal
073: *          elements of the last NB columns of the reduced matrix;
074: *          if UPLO = 'L', E(1:nb) contains the subdiagonal elements of
075: *          the first NB columns of the reduced matrix.
076: *
077: *  TAU     (output) COMPLEX array, dimension (N-1)
078: *          The scalar factors of the elementary reflectors, stored in
079: *          TAU(n-nb:n-1) if UPLO = 'U', and in TAU(1:nb) if UPLO = 'L'.
080: *          See Further Details.
081: *
082: *  W       (output) COMPLEX array, dimension (LDW,NB)
083: *          The n-by-nb matrix W required to update the unreduced part
084: *          of A.
085: *
086: *  LDW     (input) INTEGER
087: *          The leading dimension of the array W. LDW >= max(1,N).
088: *
089: *  Further Details
090: *  ===============
091: *
092: *  If UPLO = 'U', the matrix Q is represented as a product of elementary
093: *  reflectors
094: *
095: *     Q = H(n) H(n-1) . . . H(n-nb+1).
096: *
097: *  Each H(i) has the form
098: *
099: *     H(i) = I - tau * v * v'
100: *
101: *  where tau is a complex scalar, and v is a complex vector with
102: *  v(i:n) = 0 and v(i-1) = 1; v(1:i-1) is stored on exit in A(1:i-1,i),
103: *  and tau in TAU(i-1).
104: *
105: *  If UPLO = 'L', the matrix Q is represented as a product of elementary
106: *  reflectors
107: *
108: *     Q = H(1) H(2) . . . H(nb).
109: *
110: *  Each H(i) has the form
111: *
112: *     H(i) = I - tau * v * v'
113: *
114: *  where tau is a complex scalar, and v is a complex vector with
115: *  v(1:i) = 0 and v(i+1) = 1; v(i+1:n) is stored on exit in A(i+1:n,i),
116: *  and tau in TAU(i).
117: *
118: *  The elements of the vectors v together form the n-by-nb matrix V
119: *  which is needed, with W, to apply the transformation to the unreduced
120: *  part of the matrix, using a Hermitian rank-2k update of the form:
121: *  A := A - V*W' - W*V'.
122: *
123: *  The contents of A on exit are illustrated by the following examples
124: *  with n = 5 and nb = 2:
125: *
126: *  if UPLO = 'U':                       if UPLO = 'L':
127: *
128: *    (  a   a   a   v4  v5 )              (  d                  )
129: *    (      a   a   v4  v5 )              (  1   d              )
130: *    (          a   1   v5 )              (  v1  1   a          )
131: *    (              d   1  )              (  v1  v2  a   a      )
132: *    (                  d  )              (  v1  v2  a   a   a  )
133: *
134: *  where d denotes a diagonal element of the reduced matrix, a denotes
135: *  an element of the original matrix that is unchanged, and vi denotes
136: *  an element of the vector defining H(i).
137: *
138: *  =====================================================================
139: *
140: *     .. Parameters ..
141:       COMPLEX            ZERO, ONE, HALF
142:       PARAMETER          ( ZERO = ( 0.0E+0, 0.0E+0 ),
143:      $                   ONE = ( 1.0E+0, 0.0E+0 ),
144:      $                   HALF = ( 0.5E+0, 0.0E+0 ) )
145: *     ..
146: *     .. Local Scalars ..
147:       INTEGER            I, IW
148:       COMPLEX            ALPHA
149: *     ..
150: *     .. External Subroutines ..
151:       EXTERNAL           CAXPY, CGEMV, CHEMV, CLACGV, CLARFG, CSCAL
152: *     ..
153: *     .. External Functions ..
154:       LOGICAL            LSAME
155:       COMPLEX            CDOTC
156:       EXTERNAL           LSAME, CDOTC
157: *     ..
158: *     .. Intrinsic Functions ..
159:       INTRINSIC          MIN, REAL
160: *     ..
161: *     .. Executable Statements ..
162: *
163: *     Quick return if possible
164: *
165:       IF( N.LE.0 )
166:      $   RETURN
167: *
168:       IF( LSAME( UPLO, 'U' ) ) THEN
169: *
170: *        Reduce last NB columns of upper triangle
171: *
172:          DO 10 I = N, N - NB + 1, -1
173:             IW = I - N + NB
174:             IF( I.LT.N ) THEN
175: *
176: *              Update A(1:i,i)
177: *
178:                A( I, I ) = REAL( A( I, I ) )
179:                CALL CLACGV( N-I, W( I, IW+1 ), LDW )
180:                CALL CGEMV( 'No transpose', I, N-I, -ONE, A( 1, I+1 ),
181:      $                     LDA, W( I, IW+1 ), LDW, ONE, A( 1, I ), 1 )
182:                CALL CLACGV( N-I, W( I, IW+1 ), LDW )
183:                CALL CLACGV( N-I, A( I, I+1 ), LDA )
184:                CALL CGEMV( 'No transpose', I, N-I, -ONE, W( 1, IW+1 ),
185:      $                     LDW, A( I, I+1 ), LDA, ONE, A( 1, I ), 1 )
186:                CALL CLACGV( N-I, A( I, I+1 ), LDA )
187:                A( I, I ) = REAL( A( I, I ) )
188:             END IF
189:             IF( I.GT.1 ) THEN
190: *
191: *              Generate elementary reflector H(i) to annihilate
192: *              A(1:i-2,i)
193: *
194:                ALPHA = A( I-1, I )
195:                CALL CLARFG( I-1, ALPHA, A( 1, I ), 1, TAU( I-1 ) )
196:                E( I-1 ) = ALPHA
197:                A( I-1, I ) = ONE
198: *
199: *              Compute W(1:i-1,i)
200: *
201:                CALL CHEMV( 'Upper', I-1, ONE, A, LDA, A( 1, I ), 1,
202:      $                     ZERO, W( 1, IW ), 1 )
203:                IF( I.LT.N ) THEN
204:                   CALL CGEMV( 'Conjugate transpose', I-1, N-I, ONE,
205:      $                        W( 1, IW+1 ), LDW, A( 1, I ), 1, ZERO,
206:      $                        W( I+1, IW ), 1 )
207:                   CALL CGEMV( 'No transpose', I-1, N-I, -ONE,
208:      $                        A( 1, I+1 ), LDA, W( I+1, IW ), 1, ONE,
209:      $                        W( 1, IW ), 1 )
210:                   CALL CGEMV( 'Conjugate transpose', I-1, N-I, ONE,
211:      $                        A( 1, I+1 ), LDA, A( 1, I ), 1, ZERO,
212:      $                        W( I+1, IW ), 1 )
213:                   CALL CGEMV( 'No transpose', I-1, N-I, -ONE,
214:      $                        W( 1, IW+1 ), LDW, W( I+1, IW ), 1, ONE,
215:      $                        W( 1, IW ), 1 )
216:                END IF
217:                CALL CSCAL( I-1, TAU( I-1 ), W( 1, IW ), 1 )
218:                ALPHA = -HALF*TAU( I-1 )*CDOTC( I-1, W( 1, IW ), 1,
219:      $                 A( 1, I ), 1 )
220:                CALL CAXPY( I-1, ALPHA, A( 1, I ), 1, W( 1, IW ), 1 )
221:             END IF
222: *
223:    10    CONTINUE
224:       ELSE
225: *
226: *        Reduce first NB columns of lower triangle
227: *
228:          DO 20 I = 1, NB
229: *
230: *           Update A(i:n,i)
231: *
232:             A( I, I ) = REAL( A( I, I ) )
233:             CALL CLACGV( I-1, W( I, 1 ), LDW )
234:             CALL CGEMV( 'No transpose', N-I+1, I-1, -ONE, A( I, 1 ),
235:      $                  LDA, W( I, 1 ), LDW, ONE, A( I, I ), 1 )
236:             CALL CLACGV( I-1, W( I, 1 ), LDW )
237:             CALL CLACGV( I-1, A( I, 1 ), LDA )
238:             CALL CGEMV( 'No transpose', N-I+1, I-1, -ONE, W( I, 1 ),
239:      $                  LDW, A( I, 1 ), LDA, ONE, A( I, I ), 1 )
240:             CALL CLACGV( I-1, A( I, 1 ), LDA )
241:             A( I, I ) = REAL( A( I, I ) )
242:             IF( I.LT.N ) THEN
243: *
244: *              Generate elementary reflector H(i) to annihilate
245: *              A(i+2:n,i)
246: *
247:                ALPHA = A( I+1, I )
248:                CALL CLARFG( N-I, ALPHA, A( MIN( I+2, N ), I ), 1,
249:      $                      TAU( I ) )
250:                E( I ) = ALPHA
251:                A( I+1, I ) = ONE
252: *
253: *              Compute W(i+1:n,i)
254: *
255:                CALL CHEMV( 'Lower', N-I, ONE, A( I+1, I+1 ), LDA,
256:      $                     A( I+1, I ), 1, ZERO, W( I+1, I ), 1 )
257:                CALL CGEMV( 'Conjugate transpose', N-I, I-1, ONE,
258:      $                     W( I+1, 1 ), LDW, A( I+1, I ), 1, ZERO,
259:      $                     W( 1, I ), 1 )
260:                CALL CGEMV( 'No transpose', N-I, I-1, -ONE, A( I+1, 1 ),
261:      $                     LDA, W( 1, I ), 1, ONE, W( I+1, I ), 1 )
262:                CALL CGEMV( 'Conjugate transpose', N-I, I-1, ONE,
263:      $                     A( I+1, 1 ), LDA, A( I+1, I ), 1, ZERO,
264:      $                     W( 1, I ), 1 )
265:                CALL CGEMV( 'No transpose', N-I, I-1, -ONE, W( I+1, 1 ),
266:      $                     LDW, W( 1, I ), 1, ONE, W( I+1, I ), 1 )
267:                CALL CSCAL( N-I, TAU( I ), W( I+1, I ), 1 )
268:                ALPHA = -HALF*TAU( I )*CDOTC( N-I, W( I+1, I ), 1,
269:      $                 A( I+1, I ), 1 )
270:                CALL CAXPY( N-I, ALPHA, A( I+1, I ), 1, W( I+1, I ), 1 )
271:             END IF
272: *
273:    20    CONTINUE
274:       END IF
275: *
276:       RETURN
277: *
278: *     End of CLATRD
279: *
280:       END
281: