001:       SUBROUTINE ZHETD2( UPLO, N, A, LDA, D, E, TAU, INFO )
002: *
003: *  -- LAPACK 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            INFO, LDA, N
011: *     ..
012: *     .. Array Arguments ..
013:       DOUBLE PRECISION   D( * ), E( * )
014:       COMPLEX*16         A( LDA, * ), TAU( * )
015: *     ..
016: *
017: *  Purpose
018: *  =======
019: *
020: *  ZHETD2 reduces a complex Hermitian matrix A to real symmetric
021: *  tridiagonal form T by a unitary similarity transformation:
022: *  Q' * A * Q = T.
023: *
024: *  Arguments
025: *  =========
026: *
027: *  UPLO    (input) CHARACTER*1
028: *          Specifies whether the upper or lower triangular part of the
029: *          Hermitian matrix A is stored:
030: *          = 'U':  Upper triangular
031: *          = 'L':  Lower triangular
032: *
033: *  N       (input) INTEGER
034: *          The order of the matrix A.  N >= 0.
035: *
036: *  A       (input/output) COMPLEX*16 array, dimension (LDA,N)
037: *          On entry, the Hermitian matrix A.  If UPLO = 'U', the leading
038: *          n-by-n upper triangular part of A contains the upper
039: *          triangular part of the matrix A, and the strictly lower
040: *          triangular part of A is not referenced.  If UPLO = 'L', the
041: *          leading n-by-n lower triangular part of A contains the lower
042: *          triangular part of the matrix A, and the strictly upper
043: *          triangular part of A is not referenced.
044: *          On exit, if UPLO = 'U', the diagonal and first superdiagonal
045: *          of A are overwritten by the corresponding elements of the
046: *          tridiagonal matrix T, and the elements above the first
047: *          superdiagonal, with the array TAU, represent the unitary
048: *          matrix Q as a product of elementary reflectors; if UPLO
049: *          = 'L', the diagonal and first subdiagonal of A are over-
050: *          written by the corresponding elements of the tridiagonal
051: *          matrix T, and the elements below the first subdiagonal, with
052: *          the array TAU, represent the unitary matrix Q as a product
053: *          of elementary reflectors. See Further Details.
054: *
055: *  LDA     (input) INTEGER
056: *          The leading dimension of the array A.  LDA >= max(1,N).
057: *
058: *  D       (output) DOUBLE PRECISION array, dimension (N)
059: *          The diagonal elements of the tridiagonal matrix T:
060: *          D(i) = A(i,i).
061: *
062: *  E       (output) DOUBLE PRECISION array, dimension (N-1)
063: *          The off-diagonal elements of the tridiagonal matrix T:
064: *          E(i) = A(i,i+1) if UPLO = 'U', E(i) = A(i+1,i) if UPLO = 'L'.
065: *
066: *  TAU     (output) COMPLEX*16 array, dimension (N-1)
067: *          The scalar factors of the elementary reflectors (see Further
068: *          Details).
069: *
070: *  INFO    (output) INTEGER
071: *          = 0:  successful exit
072: *          < 0:  if INFO = -i, the i-th argument had an illegal value.
073: *
074: *  Further Details
075: *  ===============
076: *
077: *  If UPLO = 'U', the matrix Q is represented as a product of elementary
078: *  reflectors
079: *
080: *     Q = H(n-1) . . . H(2) H(1).
081: *
082: *  Each H(i) has the form
083: *
084: *     H(i) = I - tau * v * v'
085: *
086: *  where tau is a complex scalar, and v is a complex vector with
087: *  v(i+1:n) = 0 and v(i) = 1; v(1:i-1) is stored on exit in
088: *  A(1:i-1,i+1), and tau in TAU(i).
089: *
090: *  If UPLO = 'L', the matrix Q is represented as a product of elementary
091: *  reflectors
092: *
093: *     Q = H(1) H(2) . . . H(n-1).
094: *
095: *  Each H(i) has the form
096: *
097: *     H(i) = I - tau * v * v'
098: *
099: *  where tau is a complex scalar, and v is a complex vector with
100: *  v(1:i) = 0 and v(i+1) = 1; v(i+2:n) is stored on exit in A(i+2:n,i),
101: *  and tau in TAU(i).
102: *
103: *  The contents of A on exit are illustrated by the following examples
104: *  with n = 5:
105: *
106: *  if UPLO = 'U':                       if UPLO = 'L':
107: *
108: *    (  d   e   v2  v3  v4 )              (  d                  )
109: *    (      d   e   v3  v4 )              (  e   d              )
110: *    (          d   e   v4 )              (  v1  e   d          )
111: *    (              d   e  )              (  v1  v2  e   d      )
112: *    (                  d  )              (  v1  v2  v3  e   d  )
113: *
114: *  where d and e denote diagonal and off-diagonal elements of T, and vi
115: *  denotes an element of the vector defining H(i).
116: *
117: *  =====================================================================
118: *
119: *     .. Parameters ..
120:       COMPLEX*16         ONE, ZERO, HALF
121:       PARAMETER          ( ONE = ( 1.0D+0, 0.0D+0 ),
122:      $                   ZERO = ( 0.0D+0, 0.0D+0 ),
123:      $                   HALF = ( 0.5D+0, 0.0D+0 ) )
124: *     ..
125: *     .. Local Scalars ..
126:       LOGICAL            UPPER
127:       INTEGER            I
128:       COMPLEX*16         ALPHA, TAUI
129: *     ..
130: *     .. External Subroutines ..
131:       EXTERNAL           XERBLA, ZAXPY, ZHEMV, ZHER2, ZLARFG
132: *     ..
133: *     .. External Functions ..
134:       LOGICAL            LSAME
135:       COMPLEX*16         ZDOTC
136:       EXTERNAL           LSAME, ZDOTC
137: *     ..
138: *     .. Intrinsic Functions ..
139:       INTRINSIC          DBLE, MAX, MIN
140: *     ..
141: *     .. Executable Statements ..
142: *
143: *     Test the input parameters
144: *
145:       INFO = 0
146:       UPPER = LSAME( UPLO, 'U' )
147:       IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
148:          INFO = -1
149:       ELSE IF( N.LT.0 ) THEN
150:          INFO = -2
151:       ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
152:          INFO = -4
153:       END IF
154:       IF( INFO.NE.0 ) THEN
155:          CALL XERBLA( 'ZHETD2', -INFO )
156:          RETURN
157:       END IF
158: *
159: *     Quick return if possible
160: *
161:       IF( N.LE.0 )
162:      $   RETURN
163: *
164:       IF( UPPER ) THEN
165: *
166: *        Reduce the upper triangle of A
167: *
168:          A( N, N ) = DBLE( A( N, N ) )
169:          DO 10 I = N - 1, 1, -1
170: *
171: *           Generate elementary reflector H(i) = I - tau * v * v'
172: *           to annihilate A(1:i-1,i+1)
173: *
174:             ALPHA = A( I, I+1 )
175:             CALL ZLARFG( I, ALPHA, A( 1, I+1 ), 1, TAUI )
176:             E( I ) = ALPHA
177: *
178:             IF( TAUI.NE.ZERO ) THEN
179: *
180: *              Apply H(i) from both sides to A(1:i,1:i)
181: *
182:                A( I, I+1 ) = ONE
183: *
184: *              Compute  x := tau * A * v  storing x in TAU(1:i)
185: *
186:                CALL ZHEMV( UPLO, I, TAUI, A, LDA, A( 1, I+1 ), 1, ZERO,
187:      $                     TAU, 1 )
188: *
189: *              Compute  w := x - 1/2 * tau * (x'*v) * v
190: *
191:                ALPHA = -HALF*TAUI*ZDOTC( I, TAU, 1, A( 1, I+1 ), 1 )
192:                CALL ZAXPY( I, ALPHA, A( 1, I+1 ), 1, TAU, 1 )
193: *
194: *              Apply the transformation as a rank-2 update:
195: *                 A := A - v * w' - w * v'
196: *
197:                CALL ZHER2( UPLO, I, -ONE, A( 1, I+1 ), 1, TAU, 1, A,
198:      $                     LDA )
199: *
200:             ELSE
201:                A( I, I ) = DBLE( A( I, I ) )
202:             END IF
203:             A( I, I+1 ) = E( I )
204:             D( I+1 ) = A( I+1, I+1 )
205:             TAU( I ) = TAUI
206:    10    CONTINUE
207:          D( 1 ) = A( 1, 1 )
208:       ELSE
209: *
210: *        Reduce the lower triangle of A
211: *
212:          A( 1, 1 ) = DBLE( A( 1, 1 ) )
213:          DO 20 I = 1, N - 1
214: *
215: *           Generate elementary reflector H(i) = I - tau * v * v'
216: *           to annihilate A(i+2:n,i)
217: *
218:             ALPHA = A( I+1, I )
219:             CALL ZLARFG( N-I, ALPHA, A( MIN( I+2, N ), I ), 1, TAUI )
220:             E( I ) = ALPHA
221: *
222:             IF( TAUI.NE.ZERO ) THEN
223: *
224: *              Apply H(i) from both sides to A(i+1:n,i+1:n)
225: *
226:                A( I+1, I ) = ONE
227: *
228: *              Compute  x := tau * A * v  storing y in TAU(i:n-1)
229: *
230:                CALL ZHEMV( UPLO, N-I, TAUI, A( I+1, I+1 ), LDA,
231:      $                     A( I+1, I ), 1, ZERO, TAU( I ), 1 )
232: *
233: *              Compute  w := x - 1/2 * tau * (x'*v) * v
234: *
235:                ALPHA = -HALF*TAUI*ZDOTC( N-I, TAU( I ), 1, A( I+1, I ),
236:      $                 1 )
237:                CALL ZAXPY( N-I, ALPHA, A( I+1, I ), 1, TAU( I ), 1 )
238: *
239: *              Apply the transformation as a rank-2 update:
240: *                 A := A - v * w' - w * v'
241: *
242:                CALL ZHER2( UPLO, N-I, -ONE, A( I+1, I ), 1, TAU( I ), 1,
243:      $                     A( I+1, I+1 ), LDA )
244: *
245:             ELSE
246:                A( I+1, I+1 ) = DBLE( A( I+1, I+1 ) )
247:             END IF
248:             A( I+1, I ) = E( I )
249:             D( I ) = A( I, I )
250:             TAU( I ) = TAUI
251:    20    CONTINUE
252:          D( N ) = A( N, N )
253:       END IF
254: *
255:       RETURN
256: *
257: *     End of ZHETD2
258: *
259:       END
260: