1*> \brief \b DSYTRD 2* 3* =========== DOCUMENTATION =========== 4* 5* Online html documentation available at 6* http://www.netlib.org/lapack/explore-html/ 7* 8*> \htmlonly 9*> Download DSYTRD + dependencies 10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dsytrd.f"> 11*> [TGZ]</a> 12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dsytrd.f"> 13*> [ZIP]</a> 14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dsytrd.f"> 15*> [TXT]</a> 16*> \endhtmlonly 17* 18* Definition: 19* =========== 20* 21* SUBROUTINE DSYTRD( UPLO, N, A, LDA, D, E, TAU, WORK, LWORK, INFO ) 22* 23* .. Scalar Arguments .. 24* CHARACTER UPLO 25* INTEGER INFO, LDA, LWORK, N 26* .. 27* .. Array Arguments .. 28* DOUBLE PRECISION A( LDA, * ), D( * ), E( * ), TAU( * ), 29* $ WORK( * ) 30* .. 31* 32* 33*> \par Purpose: 34* ============= 35*> 36*> \verbatim 37*> 38*> DSYTRD reduces a real symmetric matrix A to real symmetric 39*> tridiagonal form T by an orthogonal similarity transformation: 40*> Q**T * A * Q = T. 41*> \endverbatim 42* 43* Arguments: 44* ========== 45* 46*> \param[in] UPLO 47*> \verbatim 48*> UPLO is CHARACTER*1 49*> = 'U': Upper triangle of A is stored; 50*> = 'L': Lower triangle of A is stored. 51*> \endverbatim 52*> 53*> \param[in] N 54*> \verbatim 55*> N is INTEGER 56*> The order of the matrix A. N >= 0. 57*> \endverbatim 58*> 59*> \param[in,out] A 60*> \verbatim 61*> A is DOUBLE PRECISION array, dimension (LDA,N) 62*> On entry, the symmetric matrix A. If UPLO = 'U', the leading 63*> N-by-N upper triangular part of A contains the upper 64*> triangular part of the matrix A, and the strictly lower 65*> triangular part of A is not referenced. If UPLO = 'L', the 66*> leading N-by-N lower triangular part of A contains the lower 67*> triangular part of the matrix A, and the strictly upper 68*> triangular part of A is not referenced. 69*> On exit, if UPLO = 'U', the diagonal and first superdiagonal 70*> of A are overwritten by the corresponding elements of the 71*> tridiagonal matrix T, and the elements above the first 72*> superdiagonal, with the array TAU, represent the orthogonal 73*> matrix Q as a product of elementary reflectors; if UPLO 74*> = 'L', the diagonal and first subdiagonal of A are over- 75*> written by the corresponding elements of the tridiagonal 76*> matrix T, and the elements below the first subdiagonal, with 77*> the array TAU, represent the orthogonal matrix Q as a product 78*> of elementary reflectors. See Further Details. 79*> \endverbatim 80*> 81*> \param[in] LDA 82*> \verbatim 83*> LDA is INTEGER 84*> The leading dimension of the array A. LDA >= max(1,N). 85*> \endverbatim 86*> 87*> \param[out] D 88*> \verbatim 89*> D is DOUBLE PRECISION array, dimension (N) 90*> The diagonal elements of the tridiagonal matrix T: 91*> D(i) = A(i,i). 92*> \endverbatim 93*> 94*> \param[out] E 95*> \verbatim 96*> E is DOUBLE PRECISION array, dimension (N-1) 97*> The off-diagonal elements of the tridiagonal matrix T: 98*> E(i) = A(i,i+1) if UPLO = 'U', E(i) = A(i+1,i) if UPLO = 'L'. 99*> \endverbatim 100*> 101*> \param[out] TAU 102*> \verbatim 103*> TAU is DOUBLE PRECISION array, dimension (N-1) 104*> The scalar factors of the elementary reflectors (see Further 105*> Details). 106*> \endverbatim 107*> 108*> \param[out] WORK 109*> \verbatim 110*> WORK is DOUBLE PRECISION array, dimension (MAX(1,LWORK)) 111*> On exit, if INFO = 0, WORK(1) returns the optimal LWORK. 112*> \endverbatim 113*> 114*> \param[in] LWORK 115*> \verbatim 116*> LWORK is INTEGER 117*> The dimension of the array WORK. LWORK >= 1. 118*> For optimum performance LWORK >= N*NB, where NB is the 119*> optimal blocksize. 120*> 121*> If LWORK = -1, then a workspace query is assumed; the routine 122*> only calculates the optimal size of the WORK array, returns 123*> this value as the first entry of the WORK array, and no error 124*> message related to LWORK is issued by XERBLA. 125*> \endverbatim 126*> 127*> \param[out] INFO 128*> \verbatim 129*> INFO is INTEGER 130*> = 0: successful exit 131*> < 0: if INFO = -i, the i-th argument had an illegal value 132*> \endverbatim 133* 134* Authors: 135* ======== 136* 137*> \author Univ. of Tennessee 138*> \author Univ. of California Berkeley 139*> \author Univ. of Colorado Denver 140*> \author NAG Ltd. 141* 142*> \date December 2016 143* 144*> \ingroup doubleSYcomputational 145* 146*> \par Further Details: 147* ===================== 148*> 149*> \verbatim 150*> 151*> If UPLO = 'U', the matrix Q is represented as a product of elementary 152*> reflectors 153*> 154*> Q = H(n-1) . . . H(2) H(1). 155*> 156*> Each H(i) has the form 157*> 158*> H(i) = I - tau * v * v**T 159*> 160*> where tau is a real scalar, and v is a real vector with 161*> v(i+1:n) = 0 and v(i) = 1; v(1:i-1) is stored on exit in 162*> A(1:i-1,i+1), and tau in TAU(i). 163*> 164*> If UPLO = 'L', the matrix Q is represented as a product of elementary 165*> reflectors 166*> 167*> Q = H(1) H(2) . . . H(n-1). 168*> 169*> Each H(i) has the form 170*> 171*> H(i) = I - tau * v * v**T 172*> 173*> where tau is a real scalar, and v is a real vector with 174*> v(1:i) = 0 and v(i+1) = 1; v(i+2:n) is stored on exit in A(i+2:n,i), 175*> and tau in TAU(i). 176*> 177*> The contents of A on exit are illustrated by the following examples 178*> with n = 5: 179*> 180*> if UPLO = 'U': if UPLO = 'L': 181*> 182*> ( d e v2 v3 v4 ) ( d ) 183*> ( d e v3 v4 ) ( e d ) 184*> ( d e v4 ) ( v1 e d ) 185*> ( d e ) ( v1 v2 e d ) 186*> ( d ) ( v1 v2 v3 e d ) 187*> 188*> where d and e denote diagonal and off-diagonal elements of T, and vi 189*> denotes an element of the vector defining H(i). 190*> \endverbatim 191*> 192* ===================================================================== 193 SUBROUTINE DSYTRD( UPLO, N, A, LDA, D, E, TAU, WORK, LWORK, INFO ) 194* 195* -- LAPACK computational routine (version 3.7.0) -- 196* -- LAPACK is a software package provided by Univ. of Tennessee, -- 197* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 198* December 2016 199* 200* .. Scalar Arguments .. 201 CHARACTER UPLO 202 INTEGER INFO, LDA, LWORK, N 203* .. 204* .. Array Arguments .. 205 DOUBLE PRECISION A( LDA, * ), D( * ), E( * ), TAU( * ), 206 $ WORK( * ) 207* .. 208* 209* ===================================================================== 210* 211* .. Parameters .. 212 DOUBLE PRECISION ONE 213 PARAMETER ( ONE = 1.0D+0 ) 214* .. 215* .. Local Scalars .. 216 LOGICAL LQUERY, UPPER 217 INTEGER I, IINFO, IWS, J, KK, LDWORK, LWKOPT, NB, 218 $ NBMIN, NX 219* .. 220* .. External Subroutines .. 221 EXTERNAL DLATRD, DSYR2K, DSYTD2, XERBLA 222* .. 223* .. Intrinsic Functions .. 224 INTRINSIC MAX 225* .. 226* .. External Functions .. 227 LOGICAL LSAME 228 INTEGER ILAENV 229 EXTERNAL LSAME, ILAENV 230* .. 231* .. Executable Statements .. 232* 233* Test the input parameters 234* 235 INFO = 0 236 UPPER = LSAME( UPLO, 'U' ) 237 LQUERY = ( LWORK.EQ.-1 ) 238 IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN 239 INFO = -1 240 ELSE IF( N.LT.0 ) THEN 241 INFO = -2 242 ELSE IF( LDA.LT.MAX( 1, N ) ) THEN 243 INFO = -4 244 ELSE IF( LWORK.LT.1 .AND. .NOT.LQUERY ) THEN 245 INFO = -9 246 END IF 247* 248 IF( INFO.EQ.0 ) THEN 249* 250* Determine the block size. 251* 252 NB = ILAENV( 1, 'DSYTRD', UPLO, N, -1, -1, -1 ) 253 LWKOPT = N*NB 254 WORK( 1 ) = LWKOPT 255 END IF 256* 257 IF( INFO.NE.0 ) THEN 258 CALL XERBLA( 'DSYTRD', -INFO ) 259 RETURN 260 ELSE IF( LQUERY ) THEN 261 RETURN 262 END IF 263* 264* Quick return if possible 265* 266 IF( N.EQ.0 ) THEN 267 WORK( 1 ) = 1 268 RETURN 269 END IF 270* 271 NX = N 272 IWS = 1 273 IF( NB.GT.1 .AND. NB.LT.N ) THEN 274* 275* Determine when to cross over from blocked to unblocked code 276* (last block is always handled by unblocked code). 277* 278 NX = MAX( NB, ILAENV( 3, 'DSYTRD', UPLO, N, -1, -1, -1 ) ) 279 IF( NX.LT.N ) THEN 280* 281* Determine if workspace is large enough for blocked code. 282* 283 LDWORK = N 284 IWS = LDWORK*NB 285 IF( LWORK.LT.IWS ) THEN 286* 287* Not enough workspace to use optimal NB: determine the 288* minimum value of NB, and reduce NB or force use of 289* unblocked code by setting NX = N. 290* 291 NB = MAX( LWORK / LDWORK, 1 ) 292 NBMIN = ILAENV( 2, 'DSYTRD', UPLO, N, -1, -1, -1 ) 293 IF( NB.LT.NBMIN ) 294 $ NX = N 295 END IF 296 ELSE 297 NX = N 298 END IF 299 ELSE 300 NB = 1 301 END IF 302* 303 IF( UPPER ) THEN 304* 305* Reduce the upper triangle of A. 306* Columns 1:kk are handled by the unblocked method. 307* 308 KK = N - ( ( N-NX+NB-1 ) / NB )*NB 309 DO 20 I = N - NB + 1, KK + 1, -NB 310* 311* Reduce columns i:i+nb-1 to tridiagonal form and form the 312* matrix W which is needed to update the unreduced part of 313* the matrix 314* 315 CALL DLATRD( UPLO, I+NB-1, NB, A, LDA, E, TAU, WORK, 316 $ LDWORK ) 317* 318* Update the unreduced submatrix A(1:i-1,1:i-1), using an 319* update of the form: A := A - V*W**T - W*V**T 320* 321 CALL DSYR2K( UPLO, 'No transpose', I-1, NB, -ONE, A( 1, I ), 322 $ LDA, WORK, LDWORK, ONE, A, LDA ) 323* 324* Copy superdiagonal elements back into A, and diagonal 325* elements into D 326* 327 DO 10 J = I, I + NB - 1 328 A( J-1, J ) = E( J-1 ) 329 D( J ) = A( J, J ) 330 10 CONTINUE 331 20 CONTINUE 332* 333* Use unblocked code to reduce the last or only block 334* 335 CALL DSYTD2( UPLO, KK, A, LDA, D, E, TAU, IINFO ) 336 ELSE 337* 338* Reduce the lower triangle of A 339* 340 DO 40 I = 1, N - NX, NB 341* 342* Reduce columns i:i+nb-1 to tridiagonal form and form the 343* matrix W which is needed to update the unreduced part of 344* the matrix 345* 346 CALL DLATRD( UPLO, N-I+1, NB, A( I, I ), LDA, E( I ), 347 $ TAU( I ), WORK, LDWORK ) 348* 349* Update the unreduced submatrix A(i+ib:n,i+ib:n), using 350* an update of the form: A := A - V*W**T - W*V**T 351* 352 CALL DSYR2K( UPLO, 'No transpose', N-I-NB+1, NB, -ONE, 353 $ A( I+NB, I ), LDA, WORK( NB+1 ), LDWORK, ONE, 354 $ A( I+NB, I+NB ), LDA ) 355* 356* Copy subdiagonal elements back into A, and diagonal 357* elements into D 358* 359 DO 30 J = I, I + NB - 1 360 A( J+1, J ) = E( J ) 361 D( J ) = A( J, J ) 362 30 CONTINUE 363 40 CONTINUE 364* 365* Use unblocked code to reduce the last or only block 366* 367 CALL DSYTD2( UPLO, N-I+1, A( I, I ), LDA, D( I ), E( I ), 368 $ TAU( I ), IINFO ) 369 END IF 370* 371 WORK( 1 ) = LWKOPT 372 RETURN 373* 374* End of DSYTRD 375* 376 END 377