1*> \brief \b DSYT21 2* 3* =========== DOCUMENTATION =========== 4* 5* Online html documentation available at 6* http://www.netlib.org/lapack/explore-html/ 7* 8* Definition: 9* =========== 10* 11* SUBROUTINE DSYT21( ITYPE, UPLO, N, KBAND, A, LDA, D, E, U, LDU, V, 12* LDV, TAU, WORK, RESULT ) 13* 14* .. Scalar Arguments .. 15* CHARACTER UPLO 16* INTEGER ITYPE, KBAND, LDA, LDU, LDV, N 17* .. 18* .. Array Arguments .. 19* DOUBLE PRECISION A( LDA, * ), D( * ), E( * ), RESULT( 2 ), 20* $ TAU( * ), U( LDU, * ), V( LDV, * ), WORK( * ) 21* .. 22* 23* 24*> \par Purpose: 25* ============= 26*> 27*> \verbatim 28*> 29*> DSYT21 generally checks a decomposition of the form 30*> 31*> A = U S U**T 32*> 33*> where **T means transpose, A is symmetric, U is orthogonal, and S is 34*> diagonal (if KBAND=0) or symmetric tridiagonal (if KBAND=1). 35*> 36*> If ITYPE=1, then U is represented as a dense matrix; otherwise U is 37*> expressed as a product of Householder transformations, whose vectors 38*> are stored in the array "V" and whose scaling constants are in "TAU". 39*> We shall use the letter "V" to refer to the product of Householder 40*> transformations (which should be equal to U). 41*> 42*> Specifically, if ITYPE=1, then: 43*> 44*> RESULT(1) = | A - U S U**T | / ( |A| n ulp ) and 45*> RESULT(2) = | I - U U**T | / ( n ulp ) 46*> 47*> If ITYPE=2, then: 48*> 49*> RESULT(1) = | A - V S V**T | / ( |A| n ulp ) 50*> 51*> If ITYPE=3, then: 52*> 53*> RESULT(1) = | I - V U**T | / ( n ulp ) 54*> 55*> For ITYPE > 1, the transformation U is expressed as a product 56*> V = H(1)...H(n-2), where H(j) = I - tau(j) v(j) v(j)**T and each 57*> vector v(j) has its first j elements 0 and the remaining n-j elements 58*> stored in V(j+1:n,j). 59*> \endverbatim 60* 61* Arguments: 62* ========== 63* 64*> \param[in] ITYPE 65*> \verbatim 66*> ITYPE is INTEGER 67*> Specifies the type of tests to be performed. 68*> 1: U expressed as a dense orthogonal matrix: 69*> RESULT(1) = | A - U S U**T | / ( |A| n ulp ) and 70*> RESULT(2) = | I - U U**T | / ( n ulp ) 71*> 72*> 2: U expressed as a product V of Housholder transformations: 73*> RESULT(1) = | A - V S V**T | / ( |A| n ulp ) 74*> 75*> 3: U expressed both as a dense orthogonal matrix and 76*> as a product of Housholder transformations: 77*> RESULT(1) = | I - V U**T | / ( n ulp ) 78*> \endverbatim 79*> 80*> \param[in] UPLO 81*> \verbatim 82*> UPLO is CHARACTER 83*> If UPLO='U', the upper triangle of A and V will be used and 84*> the (strictly) lower triangle will not be referenced. 85*> If UPLO='L', the lower triangle of A and V will be used and 86*> the (strictly) upper triangle will not be referenced. 87*> \endverbatim 88*> 89*> \param[in] N 90*> \verbatim 91*> N is INTEGER 92*> The size of the matrix. If it is zero, DSYT21 does nothing. 93*> It must be at least zero. 94*> \endverbatim 95*> 96*> \param[in] KBAND 97*> \verbatim 98*> KBAND is INTEGER 99*> The bandwidth of the matrix. It may only be zero or one. 100*> If zero, then S is diagonal, and E is not referenced. If 101*> one, then S is symmetric tri-diagonal. 102*> \endverbatim 103*> 104*> \param[in] A 105*> \verbatim 106*> A is DOUBLE PRECISION array, dimension (LDA, N) 107*> The original (unfactored) matrix. It is assumed to be 108*> symmetric, and only the upper (UPLO='U') or only the lower 109*> (UPLO='L') will be referenced. 110*> \endverbatim 111*> 112*> \param[in] LDA 113*> \verbatim 114*> LDA is INTEGER 115*> The leading dimension of A. It must be at least 1 116*> and at least N. 117*> \endverbatim 118*> 119*> \param[in] D 120*> \verbatim 121*> D is DOUBLE PRECISION array, dimension (N) 122*> The diagonal of the (symmetric tri-) diagonal matrix. 123*> \endverbatim 124*> 125*> \param[in] E 126*> \verbatim 127*> E is DOUBLE PRECISION array, dimension (N-1) 128*> The off-diagonal of the (symmetric tri-) diagonal matrix. 129*> E(1) is the (1,2) and (2,1) element, E(2) is the (2,3) and 130*> (3,2) element, etc. 131*> Not referenced if KBAND=0. 132*> \endverbatim 133*> 134*> \param[in] U 135*> \verbatim 136*> U is DOUBLE PRECISION array, dimension (LDU, N) 137*> If ITYPE=1 or 3, this contains the orthogonal matrix in 138*> the decomposition, expressed as a dense matrix. If ITYPE=2, 139*> then it is not referenced. 140*> \endverbatim 141*> 142*> \param[in] LDU 143*> \verbatim 144*> LDU is INTEGER 145*> The leading dimension of U. LDU must be at least N and 146*> at least 1. 147*> \endverbatim 148*> 149*> \param[in] V 150*> \verbatim 151*> V is DOUBLE PRECISION array, dimension (LDV, N) 152*> If ITYPE=2 or 3, the columns of this array contain the 153*> Householder vectors used to describe the orthogonal matrix 154*> in the decomposition. If UPLO='L', then the vectors are in 155*> the lower triangle, if UPLO='U', then in the upper 156*> triangle. 157*> *NOTE* If ITYPE=2 or 3, V is modified and restored. The 158*> subdiagonal (if UPLO='L') or the superdiagonal (if UPLO='U') 159*> is set to one, and later reset to its original value, during 160*> the course of the calculation. 161*> If ITYPE=1, then it is neither referenced nor modified. 162*> \endverbatim 163*> 164*> \param[in] LDV 165*> \verbatim 166*> LDV is INTEGER 167*> The leading dimension of V. LDV must be at least N and 168*> at least 1. 169*> \endverbatim 170*> 171*> \param[in] TAU 172*> \verbatim 173*> TAU is DOUBLE PRECISION array, dimension (N) 174*> If ITYPE >= 2, then TAU(j) is the scalar factor of 175*> v(j) v(j)**T in the Householder transformation H(j) of 176*> the product U = H(1)...H(n-2) 177*> If ITYPE < 2, then TAU is not referenced. 178*> \endverbatim 179*> 180*> \param[out] WORK 181*> \verbatim 182*> WORK is DOUBLE PRECISION array, dimension (2*N**2) 183*> \endverbatim 184*> 185*> \param[out] RESULT 186*> \verbatim 187*> RESULT is DOUBLE PRECISION array, dimension (2) 188*> The values computed by the two tests described above. The 189*> values are currently limited to 1/ulp, to avoid overflow. 190*> RESULT(1) is always modified. RESULT(2) is modified only 191*> if ITYPE=1. 192*> \endverbatim 193* 194* Authors: 195* ======== 196* 197*> \author Univ. of Tennessee 198*> \author Univ. of California Berkeley 199*> \author Univ. of Colorado Denver 200*> \author NAG Ltd. 201* 202*> \ingroup double_eig 203* 204* ===================================================================== 205 SUBROUTINE DSYT21( ITYPE, UPLO, N, KBAND, A, LDA, D, E, U, LDU, V, 206 $ LDV, TAU, WORK, RESULT ) 207* 208* -- LAPACK test routine -- 209* -- LAPACK is a software package provided by Univ. of Tennessee, -- 210* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 211* 212* .. Scalar Arguments .. 213 CHARACTER UPLO 214 INTEGER ITYPE, KBAND, LDA, LDU, LDV, N 215* .. 216* .. Array Arguments .. 217 DOUBLE PRECISION A( LDA, * ), D( * ), E( * ), RESULT( 2 ), 218 $ TAU( * ), U( LDU, * ), V( LDV, * ), WORK( * ) 219* .. 220* 221* ===================================================================== 222* 223* .. Parameters .. 224 DOUBLE PRECISION ZERO, ONE, TEN 225 PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0, TEN = 10.0D0 ) 226* .. 227* .. Local Scalars .. 228 LOGICAL LOWER 229 CHARACTER CUPLO 230 INTEGER IINFO, J, JCOL, JR, JROW 231 DOUBLE PRECISION ANORM, ULP, UNFL, VSAVE, WNORM 232* .. 233* .. External Functions .. 234 LOGICAL LSAME 235 DOUBLE PRECISION DLAMCH, DLANGE, DLANSY 236 EXTERNAL LSAME, DLAMCH, DLANGE, DLANSY 237* .. 238* .. External Subroutines .. 239 EXTERNAL DGEMM, DLACPY, DLARFY, DLASET, DORM2L, DORM2R, 240 $ DSYR, DSYR2 241* .. 242* .. Intrinsic Functions .. 243 INTRINSIC DBLE, MAX, MIN 244* .. 245* .. Executable Statements .. 246* 247 RESULT( 1 ) = ZERO 248 IF( ITYPE.EQ.1 ) 249 $ RESULT( 2 ) = ZERO 250 IF( N.LE.0 ) 251 $ RETURN 252* 253 IF( LSAME( UPLO, 'U' ) ) THEN 254 LOWER = .FALSE. 255 CUPLO = 'U' 256 ELSE 257 LOWER = .TRUE. 258 CUPLO = 'L' 259 END IF 260* 261 UNFL = DLAMCH( 'Safe minimum' ) 262 ULP = DLAMCH( 'Epsilon' )*DLAMCH( 'Base' ) 263* 264* Some Error Checks 265* 266 IF( ITYPE.LT.1 .OR. ITYPE.GT.3 ) THEN 267 RESULT( 1 ) = TEN / ULP 268 RETURN 269 END IF 270* 271* Do Test 1 272* 273* Norm of A: 274* 275 IF( ITYPE.EQ.3 ) THEN 276 ANORM = ONE 277 ELSE 278 ANORM = MAX( DLANSY( '1', CUPLO, N, A, LDA, WORK ), UNFL ) 279 END IF 280* 281* Compute error matrix: 282* 283 IF( ITYPE.EQ.1 ) THEN 284* 285* ITYPE=1: error = A - U S U**T 286* 287 CALL DLASET( 'Full', N, N, ZERO, ZERO, WORK, N ) 288 CALL DLACPY( CUPLO, N, N, A, LDA, WORK, N ) 289* 290 DO 10 J = 1, N 291 CALL DSYR( CUPLO, N, -D( J ), U( 1, J ), 1, WORK, N ) 292 10 CONTINUE 293* 294 IF( N.GT.1 .AND. KBAND.EQ.1 ) THEN 295 DO 20 J = 1, N - 1 296 CALL DSYR2( CUPLO, N, -E( J ), U( 1, J ), 1, U( 1, J+1 ), 297 $ 1, WORK, N ) 298 20 CONTINUE 299 END IF 300 WNORM = DLANSY( '1', CUPLO, N, WORK, N, WORK( N**2+1 ) ) 301* 302 ELSE IF( ITYPE.EQ.2 ) THEN 303* 304* ITYPE=2: error = V S V**T - A 305* 306 CALL DLASET( 'Full', N, N, ZERO, ZERO, WORK, N ) 307* 308 IF( LOWER ) THEN 309 WORK( N**2 ) = D( N ) 310 DO 40 J = N - 1, 1, -1 311 IF( KBAND.EQ.1 ) THEN 312 WORK( ( N+1 )*( J-1 )+2 ) = ( ONE-TAU( J ) )*E( J ) 313 DO 30 JR = J + 2, N 314 WORK( ( J-1 )*N+JR ) = -TAU( J )*E( J )*V( JR, J ) 315 30 CONTINUE 316 END IF 317* 318 VSAVE = V( J+1, J ) 319 V( J+1, J ) = ONE 320 CALL DLARFY( 'L', N-J, V( J+1, J ), 1, TAU( J ), 321 $ WORK( ( N+1 )*J+1 ), N, WORK( N**2+1 ) ) 322 V( J+1, J ) = VSAVE 323 WORK( ( N+1 )*( J-1 )+1 ) = D( J ) 324 40 CONTINUE 325 ELSE 326 WORK( 1 ) = D( 1 ) 327 DO 60 J = 1, N - 1 328 IF( KBAND.EQ.1 ) THEN 329 WORK( ( N+1 )*J ) = ( ONE-TAU( J ) )*E( J ) 330 DO 50 JR = 1, J - 1 331 WORK( J*N+JR ) = -TAU( J )*E( J )*V( JR, J+1 ) 332 50 CONTINUE 333 END IF 334* 335 VSAVE = V( J, J+1 ) 336 V( J, J+1 ) = ONE 337 CALL DLARFY( 'U', J, V( 1, J+1 ), 1, TAU( J ), WORK, N, 338 $ WORK( N**2+1 ) ) 339 V( J, J+1 ) = VSAVE 340 WORK( ( N+1 )*J+1 ) = D( J+1 ) 341 60 CONTINUE 342 END IF 343* 344 DO 90 JCOL = 1, N 345 IF( LOWER ) THEN 346 DO 70 JROW = JCOL, N 347 WORK( JROW+N*( JCOL-1 ) ) = WORK( JROW+N*( JCOL-1 ) ) 348 $ - A( JROW, JCOL ) 349 70 CONTINUE 350 ELSE 351 DO 80 JROW = 1, JCOL 352 WORK( JROW+N*( JCOL-1 ) ) = WORK( JROW+N*( JCOL-1 ) ) 353 $ - A( JROW, JCOL ) 354 80 CONTINUE 355 END IF 356 90 CONTINUE 357 WNORM = DLANSY( '1', CUPLO, N, WORK, N, WORK( N**2+1 ) ) 358* 359 ELSE IF( ITYPE.EQ.3 ) THEN 360* 361* ITYPE=3: error = U V**T - I 362* 363 IF( N.LT.2 ) 364 $ RETURN 365 CALL DLACPY( ' ', N, N, U, LDU, WORK, N ) 366 IF( LOWER ) THEN 367 CALL DORM2R( 'R', 'T', N, N-1, N-1, V( 2, 1 ), LDV, TAU, 368 $ WORK( N+1 ), N, WORK( N**2+1 ), IINFO ) 369 ELSE 370 CALL DORM2L( 'R', 'T', N, N-1, N-1, V( 1, 2 ), LDV, TAU, 371 $ WORK, N, WORK( N**2+1 ), IINFO ) 372 END IF 373 IF( IINFO.NE.0 ) THEN 374 RESULT( 1 ) = TEN / ULP 375 RETURN 376 END IF 377* 378 DO 100 J = 1, N 379 WORK( ( N+1 )*( J-1 )+1 ) = WORK( ( N+1 )*( J-1 )+1 ) - ONE 380 100 CONTINUE 381* 382 WNORM = DLANGE( '1', N, N, WORK, N, WORK( N**2+1 ) ) 383 END IF 384* 385 IF( ANORM.GT.WNORM ) THEN 386 RESULT( 1 ) = ( WNORM / ANORM ) / ( N*ULP ) 387 ELSE 388 IF( ANORM.LT.ONE ) THEN 389 RESULT( 1 ) = ( MIN( WNORM, N*ANORM ) / ANORM ) / ( N*ULP ) 390 ELSE 391 RESULT( 1 ) = MIN( WNORM / ANORM, DBLE( N ) ) / ( N*ULP ) 392 END IF 393 END IF 394* 395* Do Test 2 396* 397* Compute U U**T - I 398* 399 IF( ITYPE.EQ.1 ) THEN 400 CALL DGEMM( 'N', 'C', N, N, N, ONE, U, LDU, U, LDU, ZERO, WORK, 401 $ N ) 402* 403 DO 110 J = 1, N 404 WORK( ( N+1 )*( J-1 )+1 ) = WORK( ( N+1 )*( J-1 )+1 ) - ONE 405 110 CONTINUE 406* 407 RESULT( 2 ) = MIN( DLANGE( '1', N, N, WORK, N, 408 $ WORK( N**2+1 ) ), DBLE( N ) ) / ( N*ULP ) 409 END IF 410* 411 RETURN 412* 413* End of DSYT21 414* 415 END 416