1       SUBROUTINE ZLATRD( UPLO, N, NB, A, LDA, E, TAU, W, LDW )
  2 *
  3 *  -- LAPACK auxiliary routine (version 3.3.1) --
  4 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
  5 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
  6 *  -- April 2011                                                      --
  7 *
  8 *     .. Scalar Arguments ..
  9       CHARACTER          UPLO
 10       INTEGER            LDA, LDW, N, NB
 11 *     ..
 12 *     .. Array Arguments ..
 13       DOUBLE PRECISION   E( * )
 14       COMPLEX*16         A( LDA, * ), TAU( * ), W( LDW, * )
 15 *     ..
 16 *
 17 *  Purpose
 18 *  =======
 19 *
 20 *  ZLATRD reduces NB rows and columns of a complex Hermitian matrix A to
 21 *  Hermitian tridiagonal form by a unitary similarity
 22 *  transformation Q**H * A * Q, and returns the matrices V and W which are
 23 *  needed to apply the transformation to the unreduced part of A.
 24 *
 25 *  If UPLO = 'U', ZLATRD reduces the last NB rows and columns of a
 26 *  matrix, of which the upper triangle is supplied;
 27 *  if UPLO = 'L', ZLATRD reduces the first NB rows and columns of a
 28 *  matrix, of which the lower triangle is supplied.
 29 *
 30 *  This is an auxiliary routine called by ZHETRD.
 31 *
 32 *  Arguments
 33 *  =========
 34 *
 35 *  UPLO    (input) CHARACTER*1
 36 *          Specifies whether the upper or lower triangular part of the
 37 *          Hermitian matrix A is stored:
 38 *          = 'U': Upper triangular
 39 *          = 'L': Lower triangular
 40 *
 41 *  N       (input) INTEGER
 42 *          The order of the matrix A.
 43 *
 44 *  NB      (input) INTEGER
 45 *          The number of rows and columns to be reduced.
 46 *
 47 *  A       (input/output) COMPLEX*16 array, dimension (LDA,N)
 48 *          On entry, the Hermitian matrix A.  If UPLO = 'U', the leading
 49 *          n-by-n upper triangular part of A contains the upper
 50 *          triangular part of the matrix A, and the strictly lower
 51 *          triangular part of A is not referenced.  If UPLO = 'L', the
 52 *          leading n-by-n lower triangular part of A contains the lower
 53 *          triangular part of the matrix A, and the strictly upper
 54 *          triangular part of A is not referenced.
 55 *          On exit:
 56 *          if UPLO = 'U', the last NB columns have been reduced to
 57 *            tridiagonal form, with the diagonal elements overwriting
 58 *            the diagonal elements of A; the elements above the diagonal
 59 *            with the array TAU, represent the unitary matrix Q as a
 60 *            product of elementary reflectors;
 61 *          if UPLO = 'L', the first NB columns have been reduced to
 62 *            tridiagonal form, with the diagonal elements overwriting
 63 *            the diagonal elements of A; the elements below the diagonal
 64 *            with the array TAU, represent the  unitary matrix Q as a
 65 *            product of elementary reflectors.
 66 *          See Further Details.
 67 *
 68 *  LDA     (input) INTEGER
 69 *          The leading dimension of the array A.  LDA >= max(1,N).
 70 *
 71 *  E       (output) DOUBLE PRECISION array, dimension (N-1)
 72 *          If UPLO = 'U', E(n-nb:n-1) contains the superdiagonal
 73 *          elements of the last NB columns of the reduced matrix;
 74 *          if UPLO = 'L', E(1:nb) contains the subdiagonal elements of
 75 *          the first NB columns of the reduced matrix.
 76 *
 77 *  TAU     (output) COMPLEX*16 array, dimension (N-1)
 78 *          The scalar factors of the elementary reflectors, stored in
 79 *          TAU(n-nb:n-1) if UPLO = 'U', and in TAU(1:nb) if UPLO = 'L'.
 80 *          See Further Details.
 81 *
 82 *  W       (output) COMPLEX*16 array, dimension (LDW,NB)
 83 *          The n-by-nb matrix W required to update the unreduced part
 84 *          of A.
 85 *
 86 *  LDW     (input) INTEGER
 87 *          The leading dimension of the array W. LDW >= max(1,N).
 88 *
 89 *  Further Details
 90 *  ===============
 91 *
 92 *  If UPLO = 'U', the matrix Q is represented as a product of elementary
 93 *  reflectors
 94 *
 95 *     Q = H(n) H(n-1) . . . H(n-nb+1).
 96 *
 97 *  Each H(i) has the form
 98 *
 99 *     H(i) = I - tau * v * v**H
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**H
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**H - W*V**H.
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*16         ZERO, ONE, HALF
142       PARAMETER          ( ZERO = ( 0.0D+00.0D+0 ),
143      $                   ONE = ( 1.0D+00.0D+0 ),
144      $                   HALF = ( 0.5D+00.0D+0 ) )
145 *     ..
146 *     .. Local Scalars ..
147       INTEGER            I, IW
148       COMPLEX*16         ALPHA
149 *     ..
150 *     .. External Subroutines ..
151       EXTERNAL           ZAXPY, ZGEMV, ZHEMV, ZLACGV, ZLARFG, ZSCAL
152 *     ..
153 *     .. External Functions ..
154       LOGICAL            LSAME
155       COMPLEX*16         ZDOTC
156       EXTERNAL           LSAME, ZDOTC
157 *     ..
158 *     .. Intrinsic Functions ..
159       INTRINSIC          DBLEMIN
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 ) = DBLE( A( I, I ) )
179                CALL ZLACGV( N-I, W( I, IW+1 ), LDW )
180                CALL ZGEMV( 'No transpose', I, N-I, -ONE, A( 1, I+1 ),
181      $                     LDA, W( I, IW+1 ), LDW, ONE, A( 1, I ), 1 )
182                CALL ZLACGV( N-I, W( I, IW+1 ), LDW )
183                CALL ZLACGV( N-I, A( I, I+1 ), LDA )
184                CALL ZGEMV( 'No transpose', I, N-I, -ONE, W( 1, IW+1 ),
185      $                     LDW, A( I, I+1 ), LDA, ONE, A( 1, I ), 1 )
186                CALL ZLACGV( N-I, A( I, I+1 ), LDA )
187                A( I, I ) = DBLE( 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 ZLARFG( 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 ZHEMV( 'Upper', I-1, ONE, A, LDA, A( 1, I ), 1,
202      $                     ZERO, W( 1, IW ), 1 )
203                IF( I.LT.N ) THEN
204                   CALL ZGEMV( '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 ZGEMV( '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 ZGEMV( '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 ZGEMV( '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 ZSCAL( I-1, TAU( I-1 ), W( 1, IW ), 1 )
218                ALPHA = -HALF*TAU( I-1 )*ZDOTC( I-1, W( 1, IW ), 1,
219      $                 A( 1, I ), 1 )
220                CALL ZAXPY( 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 ) = DBLE( A( I, I ) )
233             CALL ZLACGV( I-1, W( I, 1 ), LDW )
234             CALL ZGEMV( 'No transpose', N-I+1, I-1-ONE, A( I, 1 ),
235      $                  LDA, W( I, 1 ), LDW, ONE, A( I, I ), 1 )
236             CALL ZLACGV( I-1, W( I, 1 ), LDW )
237             CALL ZLACGV( I-1, A( I, 1 ), LDA )
238             CALL ZGEMV( 'No transpose', N-I+1, I-1-ONE, W( I, 1 ),
239      $                  LDW, A( I, 1 ), LDA, ONE, A( I, I ), 1 )
240             CALL ZLACGV( I-1, A( I, 1 ), LDA )
241             A( I, I ) = DBLE( 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 ZLARFG( 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 ZHEMV( 'Lower', N-I, ONE, A( I+1, I+1 ), LDA,
256      $                     A( I+1, I ), 1, ZERO, W( I+1, I ), 1 )
257                CALL ZGEMV( 'Conjugate transpose', N-I, I-1, ONE,
258      $                     W( I+11 ), LDW, A( I+1, I ), 1, ZERO,
259      $                     W( 1, I ), 1 )
260                CALL ZGEMV( 'No transpose', N-I, I-1-ONE, A( I+11 ),
261      $                     LDA, W( 1, I ), 1, ONE, W( I+1, I ), 1 )
262                CALL ZGEMV( 'Conjugate transpose', N-I, I-1, ONE,
263      $                     A( I+11 ), LDA, A( I+1, I ), 1, ZERO,
264      $                     W( 1, I ), 1 )
265                CALL ZGEMV( 'No transpose', N-I, I-1-ONE, W( I+11 ),
266      $                     LDW, W( 1, I ), 1, ONE, W( I+1, I ), 1 )
267                CALL ZSCAL( N-I, TAU( I ), W( I+1, I ), 1 )
268                ALPHA = -HALF*TAU( I )*ZDOTC( N-I, W( I+1, I ), 1,
269      $                 A( I+1, I ), 1 )
270                CALL ZAXPY( 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 ZLATRD
279 *
280       END