LAPACK  3.10.0
LAPACK: Linear Algebra PACKage

◆ chptrd()

subroutine chptrd ( character  UPLO,
integer  N,
complex, dimension( * )  AP,
real, dimension( * )  D,
real, dimension( * )  E,
complex, dimension( * )  TAU,
integer  INFO 
)

CHPTRD

Download CHPTRD + dependencies [TGZ] [ZIP] [TXT]

Purpose:
 CHPTRD reduces a complex Hermitian matrix A stored in packed form to
 real symmetric tridiagonal form T by a unitary similarity
 transformation: Q**H * A * Q = T.
Parameters
[in]UPLO
          UPLO is CHARACTER*1
          = 'U':  Upper triangle of A is stored;
          = 'L':  Lower triangle of A is stored.
[in]N
          N is INTEGER
          The order of the matrix A.  N >= 0.
[in,out]AP
          AP is COMPLEX array, dimension (N*(N+1)/2)
          On entry, the upper or lower triangle of the Hermitian matrix
          A, packed columnwise in a linear array.  The j-th column of A
          is stored in the array AP as follows:
          if UPLO = 'U', AP(i + (j-1)*j/2) = A(i,j) for 1<=i<=j;
          if UPLO = 'L', AP(i + (j-1)*(2*n-j)/2) = A(i,j) for j<=i<=n.
          On exit, if UPLO = 'U', the diagonal and first superdiagonal
          of A are overwritten by the corresponding elements of the
          tridiagonal matrix T, and the elements above the first
          superdiagonal, with the array TAU, represent the unitary
          matrix Q as a product of elementary reflectors; if UPLO
          = 'L', the diagonal and first subdiagonal of A are over-
          written by the corresponding elements of the tridiagonal
          matrix T, and the elements below the first subdiagonal, with
          the array TAU, represent the unitary matrix Q as a product
          of elementary reflectors. See Further Details.
[out]D
          D is REAL array, dimension (N)
          The diagonal elements of the tridiagonal matrix T:
          D(i) = A(i,i).
[out]E
          E is REAL array, dimension (N-1)
          The off-diagonal elements of the tridiagonal matrix T:
          E(i) = A(i,i+1) if UPLO = 'U', E(i) = A(i+1,i) if UPLO = 'L'.
[out]TAU
          TAU is COMPLEX array, dimension (N-1)
          The scalar factors of the elementary reflectors (see Further
          Details).
[out]INFO
          INFO is INTEGER
          = 0:  successful exit
          < 0:  if INFO = -i, the i-th argument had an illegal value
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Further Details:
  If UPLO = 'U', the matrix Q is represented as a product of elementary
  reflectors

     Q = H(n-1) . . . H(2) H(1).

  Each H(i) has the form

     H(i) = I - tau * v * v**H

  where tau is a complex scalar, and v is a complex vector with
  v(i+1:n) = 0 and v(i) = 1; v(1:i-1) is stored on exit in AP,
  overwriting A(1:i-1,i+1), and tau is stored in TAU(i).

  If UPLO = 'L', the matrix Q is represented as a product of elementary
  reflectors

     Q = H(1) H(2) . . . H(n-1).

  Each H(i) has the form

     H(i) = I - tau * v * v**H

  where tau is a complex scalar, and v is a complex vector with
  v(1:i) = 0 and v(i+1) = 1; v(i+2:n) is stored on exit in AP,
  overwriting A(i+2:n,i), and tau is stored in TAU(i).

Definition at line 150 of file chptrd.f.

151 *
152 * -- LAPACK computational routine --
153 * -- LAPACK is a software package provided by Univ. of Tennessee, --
154 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
155 *
156 * .. Scalar Arguments ..
157  CHARACTER UPLO
158  INTEGER INFO, N
159 * ..
160 * .. Array Arguments ..
161  REAL D( * ), E( * )
162  COMPLEX AP( * ), TAU( * )
163 * ..
164 *
165 * =====================================================================
166 *
167 * .. Parameters ..
168  COMPLEX ONE, ZERO, HALF
169  parameter( one = ( 1.0e+0, 0.0e+0 ),
170  $ zero = ( 0.0e+0, 0.0e+0 ),
171  $ half = ( 0.5e+0, 0.0e+0 ) )
172 * ..
173 * .. Local Scalars ..
174  LOGICAL UPPER
175  INTEGER I, I1, I1I1, II
176  COMPLEX ALPHA, TAUI
177 * ..
178 * .. External Subroutines ..
179  EXTERNAL caxpy, chpmv, chpr2, clarfg, xerbla
180 * ..
181 * .. External Functions ..
182  LOGICAL LSAME
183  COMPLEX CDOTC
184  EXTERNAL lsame, cdotc
185 * ..
186 * .. Intrinsic Functions ..
187  INTRINSIC real
188 * ..
189 * .. Executable Statements ..
190 *
191 * Test the input parameters
192 *
193  info = 0
194  upper = lsame( uplo, 'U' )
195  IF( .NOT.upper .AND. .NOT.lsame( uplo, 'L' ) ) THEN
196  info = -1
197  ELSE IF( n.LT.0 ) THEN
198  info = -2
199  END IF
200  IF( info.NE.0 ) THEN
201  CALL xerbla( 'CHPTRD', -info )
202  RETURN
203  END IF
204 *
205 * Quick return if possible
206 *
207  IF( n.LE.0 )
208  $ RETURN
209 *
210  IF( upper ) THEN
211 *
212 * Reduce the upper triangle of A.
213 * I1 is the index in AP of A(1,I+1).
214 *
215  i1 = n*( n-1 ) / 2 + 1
216  ap( i1+n-1 ) = real( ap( i1+n-1 ) )
217  DO 10 i = n - 1, 1, -1
218 *
219 * Generate elementary reflector H(i) = I - tau * v * v**H
220 * to annihilate A(1:i-1,i+1)
221 *
222  alpha = ap( i1+i-1 )
223  CALL clarfg( i, alpha, ap( i1 ), 1, taui )
224  e( i ) = real( alpha )
225 *
226  IF( taui.NE.zero ) THEN
227 *
228 * Apply H(i) from both sides to A(1:i,1:i)
229 *
230  ap( i1+i-1 ) = one
231 *
232 * Compute y := tau * A * v storing y in TAU(1:i)
233 *
234  CALL chpmv( uplo, i, taui, ap, ap( i1 ), 1, zero, tau,
235  $ 1 )
236 *
237 * Compute w := y - 1/2 * tau * (y**H *v) * v
238 *
239  alpha = -half*taui*cdotc( i, tau, 1, ap( i1 ), 1 )
240  CALL caxpy( i, alpha, ap( i1 ), 1, tau, 1 )
241 *
242 * Apply the transformation as a rank-2 update:
243 * A := A - v * w**H - w * v**H
244 *
245  CALL chpr2( uplo, i, -one, ap( i1 ), 1, tau, 1, ap )
246 *
247  END IF
248  ap( i1+i-1 ) = e( i )
249  d( i+1 ) = real( ap( i1+i ) )
250  tau( i ) = taui
251  i1 = i1 - i
252  10 CONTINUE
253  d( 1 ) = real( ap( 1 ) )
254  ELSE
255 *
256 * Reduce the lower triangle of A. II is the index in AP of
257 * A(i,i) and I1I1 is the index of A(i+1,i+1).
258 *
259  ii = 1
260  ap( 1 ) = real( ap( 1 ) )
261  DO 20 i = 1, n - 1
262  i1i1 = ii + n - i + 1
263 *
264 * Generate elementary reflector H(i) = I - tau * v * v**H
265 * to annihilate A(i+2:n,i)
266 *
267  alpha = ap( ii+1 )
268  CALL clarfg( n-i, alpha, ap( ii+2 ), 1, taui )
269  e( i ) = real( alpha )
270 *
271  IF( taui.NE.zero ) THEN
272 *
273 * Apply H(i) from both sides to A(i+1:n,i+1:n)
274 *
275  ap( ii+1 ) = one
276 *
277 * Compute y := tau * A * v storing y in TAU(i:n-1)
278 *
279  CALL chpmv( uplo, n-i, taui, ap( i1i1 ), ap( ii+1 ), 1,
280  $ zero, tau( i ), 1 )
281 *
282 * Compute w := y - 1/2 * tau * (y**H *v) * v
283 *
284  alpha = -half*taui*cdotc( n-i, tau( i ), 1, ap( ii+1 ),
285  $ 1 )
286  CALL caxpy( n-i, alpha, ap( ii+1 ), 1, tau( i ), 1 )
287 *
288 * Apply the transformation as a rank-2 update:
289 * A := A - v * w**H - w * v**H
290 *
291  CALL chpr2( uplo, n-i, -one, ap( ii+1 ), 1, tau( i ), 1,
292  $ ap( i1i1 ) )
293 *
294  END IF
295  ap( ii+1 ) = e( i )
296  d( i ) = real( ap( ii ) )
297  tau( i ) = taui
298  ii = i1i1
299  20 CONTINUE
300  d( n ) = real( ap( ii ) )
301  END IF
302 *
303  RETURN
304 *
305 * End of CHPTRD
306 *
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:60
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:53
complex function cdotc(N, CX, INCX, CY, INCY)
CDOTC
Definition: cdotc.f:83
subroutine caxpy(N, CA, CX, INCX, CY, INCY)
CAXPY
Definition: caxpy.f:88
subroutine chpr2(UPLO, N, ALPHA, X, INCX, Y, INCY, AP)
CHPR2
Definition: chpr2.f:145
subroutine chpmv(UPLO, N, ALPHA, AP, X, INCX, BETA, Y, INCY)
CHPMV
Definition: chpmv.f:149
subroutine clarfg(N, ALPHA, X, INCX, TAU)
CLARFG generates an elementary reflector (Householder matrix).
Definition: clarfg.f:106
Here is the call graph for this function:
Here is the caller graph for this function: