LAPACK  3.6.1
LAPACK: Linear Algebra PACKage
subroutine dsptrd ( character  UPLO,
integer  N,
double precision, dimension( * )  AP,
double precision, dimension( * )  D,
double precision, dimension( * )  E,
double precision, dimension( * )  TAU,
integer  INFO 
)

DSPTRD

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

Purpose: 
 DSPTRD reduces a real symmetric matrix A stored in packed form to
 symmetric tridiagonal form T by an orthogonal similarity
 transformation: Q**T * 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 DOUBLE PRECISION array, dimension (N*(N+1)/2)
          On entry, the upper or lower triangle of the symmetric 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 orthogonal
          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 orthogonal matrix Q as a product
          of elementary reflectors. See Further Details.
[out]D
          D is DOUBLE PRECISION array, dimension (N)
          The diagonal elements of the tridiagonal matrix T:
          D(i) = A(i,i).
[out]E
          E is DOUBLE PRECISION 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 DOUBLE PRECISION 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.
Date
November 2011
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**T

  where tau is a real scalar, and v is a real 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**T

  where tau is a real scalar, and v is a real 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 152 of file dsptrd.f.

152 *
153 * -- LAPACK computational routine (version 3.4.0) --
154 * -- LAPACK is a software package provided by Univ. of Tennessee, --
155 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
156 * November 2011
157 *
158 * .. Scalar Arguments ..
159  CHARACTER uplo
160  INTEGER info, n
161 * ..
162 * .. Array Arguments ..
163  DOUBLE PRECISION ap( * ), d( * ), e( * ), tau( * )
164 * ..
165 *
166 * =====================================================================
167 *
168 * .. Parameters ..
169  DOUBLE PRECISION one, zero, half
170  parameter ( one = 1.0d0, zero = 0.0d0,
171  $ half = 1.0d0 / 2.0d0 )
172 * ..
173 * .. Local Scalars ..
174  LOGICAL upper
175  INTEGER i, i1, i1i1, ii
176  DOUBLE PRECISION alpha, taui
177 * ..
178 * .. External Subroutines ..
179  EXTERNAL daxpy, dlarfg, dspmv, dspr2, xerbla
180 * ..
181 * .. External Functions ..
182  LOGICAL lsame
183  DOUBLE PRECISION ddot
184  EXTERNAL lsame, ddot
185 * ..
186 * .. Executable Statements ..
187 *
188 * Test the input parameters
189 *
190  info = 0
191  upper = lsame( uplo, 'U' )
192  IF( .NOT.upper .AND. .NOT.lsame( uplo, 'L' ) ) THEN
193  info = -1
194  ELSE IF( n.LT.0 ) THEN
195  info = -2
196  END IF
197  IF( info.NE.0 ) THEN
198  CALL xerbla( 'DSPTRD', -info )
199  RETURN
200  END IF
201 *
202 * Quick return if possible
203 *
204  IF( n.LE.0 )
205  $ RETURN
206 *
207  IF( upper ) THEN
208 *
209 * Reduce the upper triangle of A.
210 * I1 is the index in AP of A(1,I+1).
211 *
212  i1 = n*( n-1 ) / 2 + 1
213  DO 10 i = n - 1, 1, -1
214 *
215 * Generate elementary reflector H(i) = I - tau * v * v**T
216 * to annihilate A(1:i-1,i+1)
217 *
218  CALL dlarfg( i, ap( i1+i-1 ), ap( i1 ), 1, taui )
219  e( i ) = ap( i1+i-1 )
220 *
221  IF( taui.NE.zero ) THEN
222 *
223 * Apply H(i) from both sides to A(1:i,1:i)
224 *
225  ap( i1+i-1 ) = one
226 *
227 * Compute y := tau * A * v storing y in TAU(1:i)
228 *
229  CALL dspmv( uplo, i, taui, ap, ap( i1 ), 1, zero, tau,
230  $ 1 )
231 *
232 * Compute w := y - 1/2 * tau * (y**T *v) * v
233 *
234  alpha = -half*taui*ddot( i, tau, 1, ap( i1 ), 1 )
235  CALL daxpy( i, alpha, ap( i1 ), 1, tau, 1 )
236 *
237 * Apply the transformation as a rank-2 update:
238 * A := A - v * w**T - w * v**T
239 *
240  CALL dspr2( uplo, i, -one, ap( i1 ), 1, tau, 1, ap )
241 *
242  ap( i1+i-1 ) = e( i )
243  END IF
244  d( i+1 ) = ap( i1+i )
245  tau( i ) = taui
246  i1 = i1 - i
247  10 CONTINUE
248  d( 1 ) = ap( 1 )
249  ELSE
250 *
251 * Reduce the lower triangle of A. II is the index in AP of
252 * A(i,i) and I1I1 is the index of A(i+1,i+1).
253 *
254  ii = 1
255  DO 20 i = 1, n - 1
256  i1i1 = ii + n - i + 1
257 *
258 * Generate elementary reflector H(i) = I - tau * v * v**T
259 * to annihilate A(i+2:n,i)
260 *
261  CALL dlarfg( n-i, ap( ii+1 ), ap( ii+2 ), 1, taui )
262  e( i ) = ap( ii+1 )
263 *
264  IF( taui.NE.zero ) THEN
265 *
266 * Apply H(i) from both sides to A(i+1:n,i+1:n)
267 *
268  ap( ii+1 ) = one
269 *
270 * Compute y := tau * A * v storing y in TAU(i:n-1)
271 *
272  CALL dspmv( uplo, n-i, taui, ap( i1i1 ), ap( ii+1 ), 1,
273  $ zero, tau( i ), 1 )
274 *
275 * Compute w := y - 1/2 * tau * (y**T *v) * v
276 *
277  alpha = -half*taui*ddot( n-i, tau( i ), 1, ap( ii+1 ),
278  $ 1 )
279  CALL daxpy( n-i, alpha, ap( ii+1 ), 1, tau( i ), 1 )
280 *
281 * Apply the transformation as a rank-2 update:
282 * A := A - v * w**T - w * v**T
283 *
284  CALL dspr2( uplo, n-i, -one, ap( ii+1 ), 1, tau( i ), 1,
285  $ ap( i1i1 ) )
286 *
287  ap( ii+1 ) = e( i )
288  END IF
289  d( i ) = ap( ii )
290  tau( i ) = taui
291  ii = i1i1
292  20 CONTINUE
293  d( n ) = ap( ii )
294  END IF
295 *
296  RETURN
297 *
298 * End of DSPTRD
299 *
ddot
double precision function ddot(N, DX, INCX, DY, INCY)
DDOT
Definition: ddot.f:53
dspmv
subroutine dspmv(UPLO, N, ALPHA, AP, X, INCX, BETA, Y, INCY)
DSPMV
Definition: dspmv.f:149
daxpy
subroutine daxpy(N, DA, DX, INCX, DY, INCY)
DAXPY
Definition: daxpy.f:54
xerbla
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
dlarfg
subroutine dlarfg(N, ALPHA, X, INCX, TAU)
DLARFG generates an elementary reflector (Householder matrix).
Definition: dlarfg.f:108
dspr2
subroutine dspr2(UPLO, N, ALPHA, X, INCX, Y, INCY, AP)
DSPR2
Definition: dspr2.f:144
lsame
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:55

Here is the call graph for this function:

Here is the caller graph for this function: