![]() |
LAPACK
3.4.1
LAPACK: Linear Algebra PACKage
|
00001 *> \brief \b DSPTRD 00002 * 00003 * =========== DOCUMENTATION =========== 00004 * 00005 * Online html documentation available at 00006 * http://www.netlib.org/lapack/explore-html/ 00007 * 00008 *> \htmlonly 00009 *> Download DSPTRD + dependencies 00010 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dsptrd.f"> 00011 *> [TGZ]</a> 00012 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dsptrd.f"> 00013 *> [ZIP]</a> 00014 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dsptrd.f"> 00015 *> [TXT]</a> 00016 *> \endhtmlonly 00017 * 00018 * Definition: 00019 * =========== 00020 * 00021 * SUBROUTINE DSPTRD( UPLO, N, AP, D, E, TAU, INFO ) 00022 * 00023 * .. Scalar Arguments .. 00024 * CHARACTER UPLO 00025 * INTEGER INFO, N 00026 * .. 00027 * .. Array Arguments .. 00028 * DOUBLE PRECISION AP( * ), D( * ), E( * ), TAU( * ) 00029 * .. 00030 * 00031 * 00032 *> \par Purpose: 00033 * ============= 00034 *> 00035 *> \verbatim 00036 *> 00037 *> DSPTRD 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 doubleOTHERcomputational 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 DSPTRD( 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 DOUBLE PRECISION AP( * ), D( * ), E( * ), TAU( * ) 00164 * .. 00165 * 00166 * ===================================================================== 00167 * 00168 * .. Parameters .. 00169 DOUBLE PRECISION ONE, ZERO, HALF 00170 PARAMETER ( ONE = 1.0D0, ZERO = 0.0D0, 00171 $ HALF = 1.0D0 / 2.0D0 ) 00172 * .. 00173 * .. Local Scalars .. 00174 LOGICAL UPPER 00175 INTEGER I, I1, I1I1, II 00176 DOUBLE PRECISION ALPHA, TAUI 00177 * .. 00178 * .. External Subroutines .. 00179 EXTERNAL DAXPY, DLARFG, DSPMV, DSPR2, XERBLA 00180 * .. 00181 * .. External Functions .. 00182 LOGICAL LSAME 00183 DOUBLE PRECISION DDOT 00184 EXTERNAL LSAME, DDOT 00185 * .. 00186 * .. Executable Statements .. 00187 * 00188 * Test the input parameters 00189 * 00190 INFO = 0 00191 UPPER = LSAME( UPLO, 'U' ) 00192 IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN 00193 INFO = -1 00194 ELSE IF( N.LT.0 ) THEN 00195 INFO = -2 00196 END IF 00197 IF( INFO.NE.0 ) THEN 00198 CALL XERBLA( 'DSPTRD', -INFO ) 00199 RETURN 00200 END IF 00201 * 00202 * Quick return if possible 00203 * 00204 IF( N.LE.0 ) 00205 $ RETURN 00206 * 00207 IF( UPPER ) THEN 00208 * 00209 * Reduce the upper triangle of A. 00210 * I1 is the index in AP of A(1,I+1). 00211 * 00212 I1 = N*( N-1 ) / 2 + 1 00213 DO 10 I = N - 1, 1, -1 00214 * 00215 * Generate elementary reflector H(i) = I - tau * v * v**T 00216 * to annihilate A(1:i-1,i+1) 00217 * 00218 CALL DLARFG( I, AP( I1+I-1 ), AP( I1 ), 1, TAUI ) 00219 E( I ) = AP( I1+I-1 ) 00220 * 00221 IF( TAUI.NE.ZERO ) THEN 00222 * 00223 * Apply H(i) from both sides to A(1:i,1:i) 00224 * 00225 AP( I1+I-1 ) = ONE 00226 * 00227 * Compute y := tau * A * v storing y in TAU(1:i) 00228 * 00229 CALL DSPMV( UPLO, I, TAUI, AP, AP( I1 ), 1, ZERO, TAU, 00230 $ 1 ) 00231 * 00232 * Compute w := y - 1/2 * tau * (y**T *v) * v 00233 * 00234 ALPHA = -HALF*TAUI*DDOT( I, TAU, 1, AP( I1 ), 1 ) 00235 CALL DAXPY( I, ALPHA, AP( I1 ), 1, TAU, 1 ) 00236 * 00237 * Apply the transformation as a rank-2 update: 00238 * A := A - v * w**T - w * v**T 00239 * 00240 CALL DSPR2( UPLO, I, -ONE, AP( I1 ), 1, TAU, 1, AP ) 00241 * 00242 AP( I1+I-1 ) = E( I ) 00243 END IF 00244 D( I+1 ) = AP( I1+I ) 00245 TAU( I ) = TAUI 00246 I1 = I1 - I 00247 10 CONTINUE 00248 D( 1 ) = AP( 1 ) 00249 ELSE 00250 * 00251 * Reduce the lower triangle of A. II is the index in AP of 00252 * A(i,i) and I1I1 is the index of A(i+1,i+1). 00253 * 00254 II = 1 00255 DO 20 I = 1, N - 1 00256 I1I1 = II + N - I + 1 00257 * 00258 * Generate elementary reflector H(i) = I - tau * v * v**T 00259 * to annihilate A(i+2:n,i) 00260 * 00261 CALL DLARFG( N-I, AP( II+1 ), AP( II+2 ), 1, TAUI ) 00262 E( I ) = AP( II+1 ) 00263 * 00264 IF( TAUI.NE.ZERO ) THEN 00265 * 00266 * Apply H(i) from both sides to A(i+1:n,i+1:n) 00267 * 00268 AP( II+1 ) = ONE 00269 * 00270 * Compute y := tau * A * v storing y in TAU(i:n-1) 00271 * 00272 CALL DSPMV( UPLO, N-I, TAUI, AP( I1I1 ), AP( II+1 ), 1, 00273 $ ZERO, TAU( I ), 1 ) 00274 * 00275 * Compute w := y - 1/2 * tau * (y**T *v) * v 00276 * 00277 ALPHA = -HALF*TAUI*DDOT( N-I, TAU( I ), 1, AP( II+1 ), 00278 $ 1 ) 00279 CALL DAXPY( N-I, ALPHA, AP( II+1 ), 1, TAU( I ), 1 ) 00280 * 00281 * Apply the transformation as a rank-2 update: 00282 * A := A - v * w**T - w * v**T 00283 * 00284 CALL DSPR2( UPLO, N-I, -ONE, AP( II+1 ), 1, TAU( I ), 1, 00285 $ AP( I1I1 ) ) 00286 * 00287 AP( II+1 ) = E( I ) 00288 END IF 00289 D( I ) = AP( II ) 00290 TAU( I ) = TAUI 00291 II = I1I1 00292 20 CONTINUE 00293 D( N ) = AP( II ) 00294 END IF 00295 * 00296 RETURN 00297 * 00298 * End of DSPTRD 00299 * 00300 END