![]() |
LAPACK
3.4.1
LAPACK: Linear Algebra PACKage
|
00001 *> \brief \b CSTT21 00002 * 00003 * =========== DOCUMENTATION =========== 00004 * 00005 * Online html documentation available at 00006 * http://www.netlib.org/lapack/explore-html/ 00007 * 00008 * Definition: 00009 * =========== 00010 * 00011 * SUBROUTINE CSTT21( N, KBAND, AD, AE, SD, SE, U, LDU, WORK, RWORK, 00012 * RESULT ) 00013 * 00014 * .. Scalar Arguments .. 00015 * INTEGER KBAND, LDU, N 00016 * .. 00017 * .. Array Arguments .. 00018 * REAL AD( * ), AE( * ), RESULT( 2 ), RWORK( * ), 00019 * $ SD( * ), SE( * ) 00020 * COMPLEX U( LDU, * ), WORK( * ) 00021 * .. 00022 * 00023 * 00024 *> \par Purpose: 00025 * ============= 00026 *> 00027 *> \verbatim 00028 *> 00029 *> CSTT21 checks a decomposition of the form 00030 *> 00031 *> A = U S UC> 00032 *> where * means conjugate transpose, A is real symmetric tridiagonal, 00033 *> U is unitary, and S is real and diagonal (if KBAND=0) or symmetric 00034 *> tridiagonal (if KBAND=1). Two tests are performed: 00035 *> 00036 *> RESULT(1) = | A - U S U* | / ( |A| n ulp ) 00037 *> 00038 *> RESULT(2) = | I - UU* | / ( n ulp ) 00039 *> \endverbatim 00040 * 00041 * Arguments: 00042 * ========== 00043 * 00044 *> \param[in] N 00045 *> \verbatim 00046 *> N is INTEGER 00047 *> The size of the matrix. If it is zero, CSTT21 does nothing. 00048 *> It must be at least zero. 00049 *> \endverbatim 00050 *> 00051 *> \param[in] KBAND 00052 *> \verbatim 00053 *> KBAND is INTEGER 00054 *> The bandwidth of the matrix S. It may only be zero or one. 00055 *> If zero, then S is diagonal, and SE is not referenced. If 00056 *> one, then S is symmetric tri-diagonal. 00057 *> \endverbatim 00058 *> 00059 *> \param[in] AD 00060 *> \verbatim 00061 *> AD is REAL array, dimension (N) 00062 *> The diagonal of the original (unfactored) matrix A. A is 00063 *> assumed to be real symmetric tridiagonal. 00064 *> \endverbatim 00065 *> 00066 *> \param[in] AE 00067 *> \verbatim 00068 *> AE is REAL array, dimension (N-1) 00069 *> The off-diagonal of the original (unfactored) matrix A. A 00070 *> is assumed to be symmetric tridiagonal. AE(1) is the (1,2) 00071 *> and (2,1) element, AE(2) is the (2,3) and (3,2) element, etc. 00072 *> \endverbatim 00073 *> 00074 *> \param[in] SD 00075 *> \verbatim 00076 *> SD is REAL array, dimension (N) 00077 *> The diagonal of the real (symmetric tri-) diagonal matrix S. 00078 *> \endverbatim 00079 *> 00080 *> \param[in] SE 00081 *> \verbatim 00082 *> SE is REAL array, dimension (N-1) 00083 *> The off-diagonal of the (symmetric tri-) diagonal matrix S. 00084 *> Not referenced if KBSND=0. If KBAND=1, then AE(1) is the 00085 *> (1,2) and (2,1) element, SE(2) is the (2,3) and (3,2) 00086 *> element, etc. 00087 *> \endverbatim 00088 *> 00089 *> \param[in] U 00090 *> \verbatim 00091 *> U is COMPLEX array, dimension (LDU, N) 00092 *> The unitary matrix in the decomposition. 00093 *> \endverbatim 00094 *> 00095 *> \param[in] LDU 00096 *> \verbatim 00097 *> LDU is INTEGER 00098 *> The leading dimension of U. LDU must be at least N. 00099 *> \endverbatim 00100 *> 00101 *> \param[out] WORK 00102 *> \verbatim 00103 *> WORK is COMPLEX array, dimension (N**2) 00104 *> \endverbatim 00105 *> 00106 *> \param[out] RWORK 00107 *> \verbatim 00108 *> RWORK is REAL array, dimension (N) 00109 *> \endverbatim 00110 *> 00111 *> \param[out] RESULT 00112 *> \verbatim 00113 *> RESULT is REAL array, dimension (2) 00114 *> The values computed by the two tests described above. The 00115 *> values are currently limited to 1/ulp, to avoid overflow. 00116 *> RESULT(1) is always modified. 00117 *> \endverbatim 00118 * 00119 * Authors: 00120 * ======== 00121 * 00122 *> \author Univ. of Tennessee 00123 *> \author Univ. of California Berkeley 00124 *> \author Univ. of Colorado Denver 00125 *> \author NAG Ltd. 00126 * 00127 *> \date November 2011 00128 * 00129 *> \ingroup complex_eig 00130 * 00131 * ===================================================================== 00132 SUBROUTINE CSTT21( N, KBAND, AD, AE, SD, SE, U, LDU, WORK, RWORK, 00133 $ RESULT ) 00134 * 00135 * -- LAPACK test routine (version 3.4.0) -- 00136 * -- LAPACK is a software package provided by Univ. of Tennessee, -- 00137 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 00138 * November 2011 00139 * 00140 * .. Scalar Arguments .. 00141 INTEGER KBAND, LDU, N 00142 * .. 00143 * .. Array Arguments .. 00144 REAL AD( * ), AE( * ), RESULT( 2 ), RWORK( * ), 00145 $ SD( * ), SE( * ) 00146 COMPLEX U( LDU, * ), WORK( * ) 00147 * .. 00148 * 00149 * ===================================================================== 00150 * 00151 * .. Parameters .. 00152 REAL ZERO, ONE 00153 PARAMETER ( ZERO = 0.0E+0, ONE = 1.0E+0 ) 00154 COMPLEX CZERO, CONE 00155 PARAMETER ( CZERO = ( 0.0E+0, 0.0E+0 ), 00156 $ CONE = ( 1.0E+0, 0.0E+0 ) ) 00157 * .. 00158 * .. Local Scalars .. 00159 INTEGER J 00160 REAL ANORM, TEMP1, TEMP2, ULP, UNFL, WNORM 00161 * .. 00162 * .. External Functions .. 00163 REAL CLANGE, CLANHE, SLAMCH 00164 EXTERNAL CLANGE, CLANHE, SLAMCH 00165 * .. 00166 * .. External Subroutines .. 00167 EXTERNAL CGEMM, CHER, CHER2, CLASET 00168 * .. 00169 * .. Intrinsic Functions .. 00170 INTRINSIC ABS, CMPLX, MAX, MIN, REAL 00171 * .. 00172 * .. Executable Statements .. 00173 * 00174 * 1) Constants 00175 * 00176 RESULT( 1 ) = ZERO 00177 RESULT( 2 ) = ZERO 00178 IF( N.LE.0 ) 00179 $ RETURN 00180 * 00181 UNFL = SLAMCH( 'Safe minimum' ) 00182 ULP = SLAMCH( 'Precision' ) 00183 * 00184 * Do Test 1 00185 * 00186 * Copy A & Compute its 1-Norm: 00187 * 00188 CALL CLASET( 'Full', N, N, CZERO, CZERO, WORK, N ) 00189 * 00190 ANORM = ZERO 00191 TEMP1 = ZERO 00192 * 00193 DO 10 J = 1, N - 1 00194 WORK( ( N+1 )*( J-1 )+1 ) = AD( J ) 00195 WORK( ( N+1 )*( J-1 )+2 ) = AE( J ) 00196 TEMP2 = ABS( AE( J ) ) 00197 ANORM = MAX( ANORM, ABS( AD( J ) )+TEMP1+TEMP2 ) 00198 TEMP1 = TEMP2 00199 10 CONTINUE 00200 * 00201 WORK( N**2 ) = AD( N ) 00202 ANORM = MAX( ANORM, ABS( AD( N ) )+TEMP1, UNFL ) 00203 * 00204 * Norm of A - USU* 00205 * 00206 DO 20 J = 1, N 00207 CALL CHER( 'L', N, -SD( J ), U( 1, J ), 1, WORK, N ) 00208 20 CONTINUE 00209 * 00210 IF( N.GT.1 .AND. KBAND.EQ.1 ) THEN 00211 DO 30 J = 1, N - 1 00212 CALL CHER2( 'L', N, -CMPLX( SE( J ) ), U( 1, J ), 1, 00213 $ U( 1, J+1 ), 1, WORK, N ) 00214 30 CONTINUE 00215 END IF 00216 * 00217 WNORM = CLANHE( '1', 'L', N, WORK, N, RWORK ) 00218 * 00219 IF( ANORM.GT.WNORM ) THEN 00220 RESULT( 1 ) = ( WNORM / ANORM ) / ( N*ULP ) 00221 ELSE 00222 IF( ANORM.LT.ONE ) THEN 00223 RESULT( 1 ) = ( MIN( WNORM, N*ANORM ) / ANORM ) / ( N*ULP ) 00224 ELSE 00225 RESULT( 1 ) = MIN( WNORM / ANORM, REAL( N ) ) / ( N*ULP ) 00226 END IF 00227 END IF 00228 * 00229 * Do Test 2 00230 * 00231 * Compute UU* - I 00232 * 00233 CALL CGEMM( 'N', 'C', N, N, N, CONE, U, LDU, U, LDU, CZERO, WORK, 00234 $ N ) 00235 * 00236 DO 40 J = 1, N 00237 WORK( ( N+1 )*( J-1 )+1 ) = WORK( ( N+1 )*( J-1 )+1 ) - CONE 00238 40 CONTINUE 00239 * 00240 RESULT( 2 ) = MIN( REAL( N ), CLANGE( '1', N, N, WORK, N, 00241 $ RWORK ) ) / ( N*ULP ) 00242 * 00243 RETURN 00244 * 00245 * End of CSTT21 00246 * 00247 END