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