![]() |
LAPACK
3.4.1
LAPACK: Linear Algebra PACKage
|
00001 *> \brief \b SSPTRD 00002 * 00003 * =========== DOCUMENTATION =========== 00004 * 00005 * Online html documentation available at 00006 * http://www.netlib.org/lapack/explore-html/ 00007 * 00008 *> \htmlonly 00009 *> Download SSPTRD + dependencies 00010 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/ssptrd.f"> 00011 *> [TGZ]</a> 00012 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/ssptrd.f"> 00013 *> [ZIP]</a> 00014 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/ssptrd.f"> 00015 *> [TXT]</a> 00016 *> \endhtmlonly 00017 * 00018 * Definition: 00019 * =========== 00020 * 00021 * SUBROUTINE SSPTRD( UPLO, N, AP, D, E, TAU, INFO ) 00022 * 00023 * .. Scalar Arguments .. 00024 * CHARACTER UPLO 00025 * INTEGER INFO, N 00026 * .. 00027 * .. Array Arguments .. 00028 * REAL AP( * ), D( * ), E( * ), TAU( * ) 00029 * .. 00030 * 00031 * 00032 *> \par Purpose: 00033 * ============= 00034 *> 00035 *> \verbatim 00036 *> 00037 *> SSPTRD reduces a real symmetric matrix A stored in packed form to 00038 *> symmetric tridiagonal form T by an orthogonal similarity 00039 *> transformation: Q**T * A * Q = T. 00040 *> \endverbatim 00041 * 00042 * Arguments: 00043 * ========== 00044 * 00045 *> \param[in] UPLO 00046 *> \verbatim 00047 *> UPLO is CHARACTER*1 00048 *> = 'U': Upper triangle of A is stored; 00049 *> = 'L': Lower triangle of A is stored. 00050 *> \endverbatim 00051 *> 00052 *> \param[in] N 00053 *> \verbatim 00054 *> N is INTEGER 00055 *> The order of the matrix A. N >= 0. 00056 *> \endverbatim 00057 *> 00058 *> \param[in,out] AP 00059 *> \verbatim 00060 *> AP is REAL array, dimension (N*(N+1)/2) 00061 *> On entry, the upper or lower triangle of the symmetric matrix 00062 *> A, packed columnwise in a linear array. The j-th column of A 00063 *> is stored in the array AP as follows: 00064 *> if UPLO = 'U', AP(i + (j-1)*j/2) = A(i,j) for 1<=i<=j; 00065 *> if UPLO = 'L', AP(i + (j-1)*(2*n-j)/2) = A(i,j) for j<=i<=n. 00066 *> On exit, if UPLO = 'U', the diagonal and first superdiagonal 00067 *> of A are overwritten by the corresponding elements of the 00068 *> tridiagonal matrix T, and the elements above the first 00069 *> superdiagonal, with the array TAU, represent the orthogonal 00070 *> matrix Q as a product of elementary reflectors; if UPLO 00071 *> = 'L', the diagonal and first subdiagonal of A are over- 00072 *> written by the corresponding elements of the tridiagonal 00073 *> matrix T, and the elements below the first subdiagonal, with 00074 *> the array TAU, represent the orthogonal matrix Q as a product 00075 *> of elementary reflectors. See Further Details. 00076 *> \endverbatim 00077 *> 00078 *> \param[out] D 00079 *> \verbatim 00080 *> D is REAL array, dimension (N) 00081 *> The diagonal elements of the tridiagonal matrix T: 00082 *> D(i) = A(i,i). 00083 *> \endverbatim 00084 *> 00085 *> \param[out] E 00086 *> \verbatim 00087 *> E is REAL array, dimension (N-1) 00088 *> The off-diagonal elements of the tridiagonal matrix T: 00089 *> E(i) = A(i,i+1) if UPLO = 'U', E(i) = A(i+1,i) if UPLO = 'L'. 00090 *> \endverbatim 00091 *> 00092 *> \param[out] TAU 00093 *> \verbatim 00094 *> TAU is REAL array, dimension (N-1) 00095 *> The scalar factors of the elementary reflectors (see Further 00096 *> Details). 00097 *> \endverbatim 00098 *> 00099 *> \param[out] INFO 00100 *> \verbatim 00101 *> INFO is INTEGER 00102 *> = 0: successful exit 00103 *> < 0: if INFO = -i, the i-th argument had an illegal value 00104 *> \endverbatim 00105 * 00106 * Authors: 00107 * ======== 00108 * 00109 *> \author Univ. of Tennessee 00110 *> \author Univ. of California Berkeley 00111 *> \author Univ. of Colorado Denver 00112 *> \author NAG Ltd. 00113 * 00114 *> \date November 2011 00115 * 00116 *> \ingroup realOTHERcomputational 00117 * 00118 *> \par Further Details: 00119 * ===================== 00120 *> 00121 *> \verbatim 00122 *> 00123 *> If UPLO = 'U', the matrix Q is represented as a product of elementary 00124 *> reflectors 00125 *> 00126 *> Q = H(n-1) . . . H(2) H(1). 00127 *> 00128 *> Each H(i) has the form 00129 *> 00130 *> H(i) = I - tau * v * v**T 00131 *> 00132 *> where tau is a real scalar, and v is a real vector with 00133 *> v(i+1:n) = 0 and v(i) = 1; v(1:i-1) is stored on exit in AP, 00134 *> overwriting A(1:i-1,i+1), and tau is stored in TAU(i). 00135 *> 00136 *> If UPLO = 'L', the matrix Q is represented as a product of elementary 00137 *> reflectors 00138 *> 00139 *> Q = H(1) H(2) . . . H(n-1). 00140 *> 00141 *> Each H(i) has the form 00142 *> 00143 *> H(i) = I - tau * v * v**T 00144 *> 00145 *> where tau is a real scalar, and v is a real vector with 00146 *> v(1:i) = 0 and v(i+1) = 1; v(i+2:n) is stored on exit in AP, 00147 *> overwriting A(i+2:n,i), and tau is stored in TAU(i). 00148 *> \endverbatim 00149 *> 00150 * ===================================================================== 00151 SUBROUTINE SSPTRD( UPLO, N, AP, D, E, TAU, INFO ) 00152 * 00153 * -- LAPACK computational routine (version 3.4.0) -- 00154 * -- LAPACK is a software package provided by Univ. of Tennessee, -- 00155 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 00156 * November 2011 00157 * 00158 * .. Scalar Arguments .. 00159 CHARACTER UPLO 00160 INTEGER INFO, N 00161 * .. 00162 * .. Array Arguments .. 00163 REAL AP( * ), D( * ), E( * ), TAU( * ) 00164 * .. 00165 * 00166 * ===================================================================== 00167 * 00168 * .. Parameters .. 00169 REAL ONE, ZERO, HALF 00170 PARAMETER ( ONE = 1.0, ZERO = 0.0, HALF = 1.0 / 2.0 ) 00171 * .. 00172 * .. Local Scalars .. 00173 LOGICAL UPPER 00174 INTEGER I, I1, I1I1, II 00175 REAL ALPHA, TAUI 00176 * .. 00177 * .. External Subroutines .. 00178 EXTERNAL SAXPY, SLARFG, SSPMV, SSPR2, XERBLA 00179 * .. 00180 * .. External Functions .. 00181 LOGICAL LSAME 00182 REAL SDOT 00183 EXTERNAL LSAME, SDOT 00184 * .. 00185 * .. Executable Statements .. 00186 * 00187 * Test the input parameters 00188 * 00189 INFO = 0 00190 UPPER = LSAME( UPLO, 'U' ) 00191 IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN 00192 INFO = -1 00193 ELSE IF( N.LT.0 ) THEN 00194 INFO = -2 00195 END IF 00196 IF( INFO.NE.0 ) THEN 00197 CALL XERBLA( 'SSPTRD', -INFO ) 00198 RETURN 00199 END IF 00200 * 00201 * Quick return if possible 00202 * 00203 IF( N.LE.0 ) 00204 $ RETURN 00205 * 00206 IF( UPPER ) THEN 00207 * 00208 * Reduce the upper triangle of A. 00209 * I1 is the index in AP of A(1,I+1). 00210 * 00211 I1 = N*( N-1 ) / 2 + 1 00212 DO 10 I = N - 1, 1, -1 00213 * 00214 * Generate elementary reflector H(i) = I - tau * v * v**T 00215 * to annihilate A(1:i-1,i+1) 00216 * 00217 CALL SLARFG( I, AP( I1+I-1 ), AP( I1 ), 1, TAUI ) 00218 E( I ) = AP( I1+I-1 ) 00219 * 00220 IF( TAUI.NE.ZERO ) THEN 00221 * 00222 * Apply H(i) from both sides to A(1:i,1:i) 00223 * 00224 AP( I1+I-1 ) = ONE 00225 * 00226 * Compute y := tau * A * v storing y in TAU(1:i) 00227 * 00228 CALL SSPMV( UPLO, I, TAUI, AP, AP( I1 ), 1, ZERO, TAU, 00229 $ 1 ) 00230 * 00231 * Compute w := y - 1/2 * tau * (y**T *v) * v 00232 * 00233 ALPHA = -HALF*TAUI*SDOT( I, TAU, 1, AP( I1 ), 1 ) 00234 CALL SAXPY( I, ALPHA, AP( I1 ), 1, TAU, 1 ) 00235 * 00236 * Apply the transformation as a rank-2 update: 00237 * A := A - v * w**T - w * v**T 00238 * 00239 CALL SSPR2( UPLO, I, -ONE, AP( I1 ), 1, TAU, 1, AP ) 00240 * 00241 AP( I1+I-1 ) = E( I ) 00242 END IF 00243 D( I+1 ) = AP( I1+I ) 00244 TAU( I ) = TAUI 00245 I1 = I1 - I 00246 10 CONTINUE 00247 D( 1 ) = AP( 1 ) 00248 ELSE 00249 * 00250 * Reduce the lower triangle of A. II is the index in AP of 00251 * A(i,i) and I1I1 is the index of A(i+1,i+1). 00252 * 00253 II = 1 00254 DO 20 I = 1, N - 1 00255 I1I1 = II + N - I + 1 00256 * 00257 * Generate elementary reflector H(i) = I - tau * v * v**T 00258 * to annihilate A(i+2:n,i) 00259 * 00260 CALL SLARFG( N-I, AP( II+1 ), AP( II+2 ), 1, TAUI ) 00261 E( I ) = AP( II+1 ) 00262 * 00263 IF( TAUI.NE.ZERO ) THEN 00264 * 00265 * Apply H(i) from both sides to A(i+1:n,i+1:n) 00266 * 00267 AP( II+1 ) = ONE 00268 * 00269 * Compute y := tau * A * v storing y in TAU(i:n-1) 00270 * 00271 CALL SSPMV( UPLO, N-I, TAUI, AP( I1I1 ), AP( II+1 ), 1, 00272 $ ZERO, TAU( I ), 1 ) 00273 * 00274 * Compute w := y - 1/2 * tau * (y**T *v) * v 00275 * 00276 ALPHA = -HALF*TAUI*SDOT( N-I, TAU( I ), 1, AP( II+1 ), 00277 $ 1 ) 00278 CALL SAXPY( N-I, ALPHA, AP( II+1 ), 1, TAU( I ), 1 ) 00279 * 00280 * Apply the transformation as a rank-2 update: 00281 * A := A - v * w**T - w * v**T 00282 * 00283 CALL SSPR2( UPLO, N-I, -ONE, AP( II+1 ), 1, TAU( I ), 1, 00284 $ AP( I1I1 ) ) 00285 * 00286 AP( II+1 ) = E( I ) 00287 END IF 00288 D( I ) = AP( II ) 00289 TAU( I ) = TAUI 00290 II = I1I1 00291 20 CONTINUE 00292 D( N ) = AP( II ) 00293 END IF 00294 * 00295 RETURN 00296 * 00297 * End of SSPTRD 00298 * 00299 END