1*> \brief \b CSTT21 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 CSTT21( N, KBAND, AD, AE, SD, SE, U, LDU, WORK, RWORK, 12* RESULT ) 13* 14* .. Scalar Arguments .. 15* INTEGER KBAND, LDU, N 16* .. 17* .. Array Arguments .. 18* REAL AD( * ), AE( * ), RESULT( 2 ), RWORK( * ), 19* $ SD( * ), SE( * ) 20* COMPLEX U( LDU, * ), WORK( * ) 21* .. 22* 23* 24*> \par Purpose: 25* ============= 26*> 27*> \verbatim 28*> 29*> CSTT21 checks a decomposition of the form 30*> 31*> A = U S U**H 32*> 33*> where **H means conjugate transpose, A is real symmetric tridiagonal, 34*> U is unitary, and S is real and diagonal (if KBAND=0) or symmetric 35*> tridiagonal (if KBAND=1). Two tests are performed: 36*> 37*> RESULT(1) = | A - U S U**H | / ( |A| n ulp ) 38*> 39*> RESULT(2) = | I - U U**H | / ( n ulp ) 40*> \endverbatim 41* 42* Arguments: 43* ========== 44* 45*> \param[in] N 46*> \verbatim 47*> N is INTEGER 48*> The size of the matrix. If it is zero, CSTT21 does nothing. 49*> It must be at least zero. 50*> \endverbatim 51*> 52*> \param[in] KBAND 53*> \verbatim 54*> KBAND is INTEGER 55*> The bandwidth of the matrix S. It may only be zero or one. 56*> If zero, then S is diagonal, and SE is not referenced. If 57*> one, then S is symmetric tri-diagonal. 58*> \endverbatim 59*> 60*> \param[in] AD 61*> \verbatim 62*> AD is REAL array, dimension (N) 63*> The diagonal of the original (unfactored) matrix A. A is 64*> assumed to be real symmetric tridiagonal. 65*> \endverbatim 66*> 67*> \param[in] AE 68*> \verbatim 69*> AE is REAL array, dimension (N-1) 70*> The off-diagonal of the original (unfactored) matrix A. A 71*> is assumed to be symmetric tridiagonal. AE(1) is the (1,2) 72*> and (2,1) element, AE(2) is the (2,3) and (3,2) element, etc. 73*> \endverbatim 74*> 75*> \param[in] SD 76*> \verbatim 77*> SD is REAL array, dimension (N) 78*> The diagonal of the real (symmetric tri-) diagonal matrix S. 79*> \endverbatim 80*> 81*> \param[in] SE 82*> \verbatim 83*> SE is REAL array, dimension (N-1) 84*> The off-diagonal of the (symmetric tri-) diagonal matrix S. 85*> Not referenced if KBSND=0. If KBAND=1, then AE(1) is the 86*> (1,2) and (2,1) element, SE(2) is the (2,3) and (3,2) 87*> element, etc. 88*> \endverbatim 89*> 90*> \param[in] U 91*> \verbatim 92*> U is COMPLEX array, dimension (LDU, N) 93*> The unitary matrix in the decomposition. 94*> \endverbatim 95*> 96*> \param[in] LDU 97*> \verbatim 98*> LDU is INTEGER 99*> The leading dimension of U. LDU must be at least N. 100*> \endverbatim 101*> 102*> \param[out] WORK 103*> \verbatim 104*> WORK is COMPLEX array, dimension (N**2) 105*> \endverbatim 106*> 107*> \param[out] RWORK 108*> \verbatim 109*> RWORK is REAL array, dimension (N) 110*> \endverbatim 111*> 112*> \param[out] RESULT 113*> \verbatim 114*> RESULT is REAL array, dimension (2) 115*> The values computed by the two tests described above. The 116*> values are currently limited to 1/ulp, to avoid overflow. 117*> RESULT(1) is always modified. 118*> \endverbatim 119* 120* Authors: 121* ======== 122* 123*> \author Univ. of Tennessee 124*> \author Univ. of California Berkeley 125*> \author Univ. of Colorado Denver 126*> \author NAG Ltd. 127* 128*> \ingroup complex_eig 129* 130* ===================================================================== 131 SUBROUTINE CSTT21( N, KBAND, AD, AE, SD, SE, U, LDU, WORK, RWORK, 132 $ RESULT ) 133* 134* -- LAPACK test routine -- 135* -- LAPACK is a software package provided by Univ. of Tennessee, -- 136* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 137* 138* .. Scalar Arguments .. 139 INTEGER KBAND, LDU, N 140* .. 141* .. Array Arguments .. 142 REAL AD( * ), AE( * ), RESULT( 2 ), RWORK( * ), 143 $ SD( * ), SE( * ) 144 COMPLEX U( LDU, * ), WORK( * ) 145* .. 146* 147* ===================================================================== 148* 149* .. Parameters .. 150 REAL ZERO, ONE 151 PARAMETER ( ZERO = 0.0E+0, ONE = 1.0E+0 ) 152 COMPLEX CZERO, CONE 153 PARAMETER ( CZERO = ( 0.0E+0, 0.0E+0 ), 154 $ CONE = ( 1.0E+0, 0.0E+0 ) ) 155* .. 156* .. Local Scalars .. 157 INTEGER J 158 REAL ANORM, TEMP1, TEMP2, ULP, UNFL, WNORM 159* .. 160* .. External Functions .. 161 REAL CLANGE, CLANHE, SLAMCH 162 EXTERNAL CLANGE, CLANHE, SLAMCH 163* .. 164* .. External Subroutines .. 165 EXTERNAL CGEMM, CHER, CHER2, CLASET 166* .. 167* .. Intrinsic Functions .. 168 INTRINSIC ABS, CMPLX, MAX, MIN, REAL 169* .. 170* .. Executable Statements .. 171* 172* 1) Constants 173* 174 RESULT( 1 ) = ZERO 175 RESULT( 2 ) = ZERO 176 IF( N.LE.0 ) 177 $ RETURN 178* 179 UNFL = SLAMCH( 'Safe minimum' ) 180 ULP = SLAMCH( 'Precision' ) 181* 182* Do Test 1 183* 184* Copy A & Compute its 1-Norm: 185* 186 CALL CLASET( 'Full', N, N, CZERO, CZERO, WORK, N ) 187* 188 ANORM = ZERO 189 TEMP1 = ZERO 190* 191 DO 10 J = 1, N - 1 192 WORK( ( N+1 )*( J-1 )+1 ) = AD( J ) 193 WORK( ( N+1 )*( J-1 )+2 ) = AE( J ) 194 TEMP2 = ABS( AE( J ) ) 195 ANORM = MAX( ANORM, ABS( AD( J ) )+TEMP1+TEMP2 ) 196 TEMP1 = TEMP2 197 10 CONTINUE 198* 199 WORK( N**2 ) = AD( N ) 200 ANORM = MAX( ANORM, ABS( AD( N ) )+TEMP1, UNFL ) 201* 202* Norm of A - U S U**H 203* 204 DO 20 J = 1, N 205 CALL CHER( 'L', N, -SD( J ), U( 1, J ), 1, WORK, N ) 206 20 CONTINUE 207* 208 IF( N.GT.1 .AND. KBAND.EQ.1 ) THEN 209 DO 30 J = 1, N - 1 210 CALL CHER2( 'L', N, -CMPLX( SE( J ) ), U( 1, J ), 1, 211 $ U( 1, J+1 ), 1, WORK, N ) 212 30 CONTINUE 213 END IF 214* 215 WNORM = CLANHE( '1', 'L', N, WORK, N, RWORK ) 216* 217 IF( ANORM.GT.WNORM ) THEN 218 RESULT( 1 ) = ( WNORM / ANORM ) / ( N*ULP ) 219 ELSE 220 IF( ANORM.LT.ONE ) THEN 221 RESULT( 1 ) = ( MIN( WNORM, N*ANORM ) / ANORM ) / ( N*ULP ) 222 ELSE 223 RESULT( 1 ) = MIN( WNORM / ANORM, REAL( N ) ) / ( N*ULP ) 224 END IF 225 END IF 226* 227* Do Test 2 228* 229* Compute U U**H - I 230* 231 CALL CGEMM( 'N', 'C', N, N, N, CONE, U, LDU, U, LDU, CZERO, WORK, 232 $ N ) 233* 234 DO 40 J = 1, N 235 WORK( ( N+1 )*( J-1 )+1 ) = WORK( ( N+1 )*( J-1 )+1 ) - CONE 236 40 CONTINUE 237* 238 RESULT( 2 ) = MIN( REAL( N ), CLANGE( '1', N, N, WORK, N, 239 $ RWORK ) ) / ( N*ULP ) 240* 241 RETURN 242* 243* End of CSTT21 244* 245 END 246