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