aboutsummaryrefslogtreecommitdiff
path: root/SRC/ssptrd.f
diff options
context:
space:
mode:
authorjason <jason@8a072113-8704-0410-8d35-dd094bca7971>2008-10-28 01:38:50 +0000
committerjason <jason@8a072113-8704-0410-8d35-dd094bca7971>2008-10-28 01:38:50 +0000
commitbaba851215b44ac3b60b9248eb02bcce7eb76247 (patch)
tree8c0f5c006875532a30d4409f5e94b0f310ff00a7 /SRC/ssptrd.f
Move LAPACK trunk into position.
Diffstat (limited to 'SRC/ssptrd.f')
-rw-r--r--SRC/ssptrd.f227
1 files changed, 227 insertions, 0 deletions
diff --git a/SRC/ssptrd.f b/SRC/ssptrd.f
new file mode 100644
index 00000000..68d3d459
--- /dev/null
+++ b/SRC/ssptrd.f
@@ -0,0 +1,227 @@
+ SUBROUTINE SSPTRD( UPLO, N, AP, D, E, TAU, INFO )
+*
+* -- LAPACK routine (version 3.1) --
+* Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
+* November 2006
+*
+* .. Scalar Arguments ..
+ CHARACTER UPLO
+ INTEGER INFO, N
+* ..
+* .. Array Arguments ..
+ REAL AP( * ), D( * ), E( * ), TAU( * )
+* ..
+*
+* Purpose
+* =======
+*
+* SSPTRD 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.
+*
+* Arguments
+* =========
+*
+* UPLO (input) CHARACTER*1
+* = 'U': Upper triangle of A is stored;
+* = 'L': Lower triangle of A is stored.
+*
+* N (input) INTEGER
+* The order of the matrix A. N >= 0.
+*
+* AP (input/output) REAL 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.
+*
+* D (output) REAL array, dimension (N)
+* The diagonal elements of the tridiagonal matrix T:
+* D(i) = A(i,i).
+*
+* E (output) 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'.
+*
+* TAU (output) REAL array, dimension (N-1)
+* The scalar factors of the elementary reflectors (see Further
+* Details).
+*
+* INFO (output) INTEGER
+* = 0: successful exit
+* < 0: if INFO = -i, the i-th argument had an illegal value
+*
+* 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'
+*
+* 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'
+*
+* 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).
+*
+* =====================================================================
+*
+* .. Parameters ..
+ REAL ONE, ZERO, HALF
+ PARAMETER ( ONE = 1.0, ZERO = 0.0, HALF = 1.0 / 2.0 )
+* ..
+* .. Local Scalars ..
+ LOGICAL UPPER
+ INTEGER I, I1, I1I1, II
+ REAL ALPHA, TAUI
+* ..
+* .. External Subroutines ..
+ EXTERNAL SAXPY, SLARFG, SSPMV, SSPR2, XERBLA
+* ..
+* .. External Functions ..
+ LOGICAL LSAME
+ REAL SDOT
+ EXTERNAL LSAME, SDOT
+* ..
+* .. Executable Statements ..
+*
+* Test the input parameters
+*
+ INFO = 0
+ UPPER = LSAME( UPLO, 'U' )
+ IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
+ INFO = -1
+ ELSE IF( N.LT.0 ) THEN
+ INFO = -2
+ END IF
+ IF( INFO.NE.0 ) THEN
+ CALL XERBLA( 'SSPTRD', -INFO )
+ RETURN
+ END IF
+*
+* Quick return if possible
+*
+ IF( N.LE.0 )
+ $ RETURN
+*
+ IF( UPPER ) THEN
+*
+* Reduce the upper triangle of A.
+* I1 is the index in AP of A(1,I+1).
+*
+ I1 = N*( N-1 ) / 2 + 1
+ DO 10 I = N - 1, 1, -1
+*
+* Generate elementary reflector H(i) = I - tau * v * v'
+* to annihilate A(1:i-1,i+1)
+*
+ CALL SLARFG( I, AP( I1+I-1 ), AP( I1 ), 1, TAUI )
+ E( I ) = AP( I1+I-1 )
+*
+ IF( TAUI.NE.ZERO ) THEN
+*
+* Apply H(i) from both sides to A(1:i,1:i)
+*
+ AP( I1+I-1 ) = ONE
+*
+* Compute y := tau * A * v storing y in TAU(1:i)
+*
+ CALL SSPMV( UPLO, I, TAUI, AP, AP( I1 ), 1, ZERO, TAU,
+ $ 1 )
+*
+* Compute w := y - 1/2 * tau * (y'*v) * v
+*
+ ALPHA = -HALF*TAUI*SDOT( I, TAU, 1, AP( I1 ), 1 )
+ CALL SAXPY( I, ALPHA, AP( I1 ), 1, TAU, 1 )
+*
+* Apply the transformation as a rank-2 update:
+* A := A - v * w' - w * v'
+*
+ CALL SSPR2( UPLO, I, -ONE, AP( I1 ), 1, TAU, 1, AP )
+*
+ AP( I1+I-1 ) = E( I )
+ END IF
+ D( I+1 ) = AP( I1+I )
+ TAU( I ) = TAUI
+ I1 = I1 - I
+ 10 CONTINUE
+ D( 1 ) = AP( 1 )
+ ELSE
+*
+* Reduce the lower triangle of A. II is the index in AP of
+* A(i,i) and I1I1 is the index of A(i+1,i+1).
+*
+ II = 1
+ DO 20 I = 1, N - 1
+ I1I1 = II + N - I + 1
+*
+* Generate elementary reflector H(i) = I - tau * v * v'
+* to annihilate A(i+2:n,i)
+*
+ CALL SLARFG( N-I, AP( II+1 ), AP( II+2 ), 1, TAUI )
+ E( I ) = AP( II+1 )
+*
+ IF( TAUI.NE.ZERO ) THEN
+*
+* Apply H(i) from both sides to A(i+1:n,i+1:n)
+*
+ AP( II+1 ) = ONE
+*
+* Compute y := tau * A * v storing y in TAU(i:n-1)
+*
+ CALL SSPMV( UPLO, N-I, TAUI, AP( I1I1 ), AP( II+1 ), 1,
+ $ ZERO, TAU( I ), 1 )
+*
+* Compute w := y - 1/2 * tau * (y'*v) * v
+*
+ ALPHA = -HALF*TAUI*SDOT( N-I, TAU( I ), 1, AP( II+1 ),
+ $ 1 )
+ CALL SAXPY( N-I, ALPHA, AP( II+1 ), 1, TAU( I ), 1 )
+*
+* Apply the transformation as a rank-2 update:
+* A := A - v * w' - w * v'
+*
+ CALL SSPR2( UPLO, N-I, -ONE, AP( II+1 ), 1, TAU( I ), 1,
+ $ AP( I1I1 ) )
+*
+ AP( II+1 ) = E( I )
+ END IF
+ D( I ) = AP( II )
+ TAU( I ) = TAUI
+ II = I1I1
+ 20 CONTINUE
+ D( N ) = AP( II )
+ END IF
+*
+ RETURN
+*
+* End of SSPTRD
+*
+ END