![]() |
LAPACK
3.4.1
LAPACK: Linear Algebra PACKage
|
00001 *> \brief \b DSTT21 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 DSTT21( N, KBAND, AD, AE, SD, SE, U, LDU, WORK, 00012 * RESULT ) 00013 * 00014 * .. Scalar Arguments .. 00015 * INTEGER KBAND, LDU, N 00016 * .. 00017 * .. Array Arguments .. 00018 * DOUBLE PRECISION AD( * ), AE( * ), RESULT( 2 ), SD( * ), 00019 * $ SE( * ), U( LDU, * ), WORK( * ) 00020 * .. 00021 * 00022 * 00023 *> \par Purpose: 00024 * ============= 00025 *> 00026 *> \verbatim 00027 *> 00028 *> DSTT21 checks a decomposition of the form 00029 *> 00030 *> A = U S U' 00031 *> 00032 *> where ' means transpose, A is symmetric tridiagonal, U is orthogonal, 00033 *> and S is diagonal (if KBAND=0) or symmetric tridiagonal (if KBAND=1). 00034 *> 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, DSTT21 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 DOUBLE PRECISION array, dimension (N) 00062 *> The diagonal of the original (unfactored) matrix A. A is 00063 *> assumed to be symmetric tridiagonal. 00064 *> \endverbatim 00065 *> 00066 *> \param[in] AE 00067 *> \verbatim 00068 *> AE is DOUBLE PRECISION 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 DOUBLE PRECISION array, dimension (N) 00077 *> The diagonal of the (symmetric tri-) diagonal matrix S. 00078 *> \endverbatim 00079 *> 00080 *> \param[in] SE 00081 *> \verbatim 00082 *> SE is DOUBLE PRECISION 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 DOUBLE PRECISION array, dimension (LDU, N) 00092 *> The orthogonal 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 DOUBLE PRECISION array, dimension (N*(N+1)) 00104 *> \endverbatim 00105 *> 00106 *> \param[out] RESULT 00107 *> \verbatim 00108 *> RESULT is DOUBLE PRECISION array, dimension (2) 00109 *> The values computed by the two tests described above. The 00110 *> values are currently limited to 1/ulp, to avoid overflow. 00111 *> RESULT(1) is always modified. 00112 *> \endverbatim 00113 * 00114 * Authors: 00115 * ======== 00116 * 00117 *> \author Univ. of Tennessee 00118 *> \author Univ. of California Berkeley 00119 *> \author Univ. of Colorado Denver 00120 *> \author NAG Ltd. 00121 * 00122 *> \date November 2011 00123 * 00124 *> \ingroup double_eig 00125 * 00126 * ===================================================================== 00127 SUBROUTINE DSTT21( N, KBAND, AD, AE, SD, SE, U, LDU, WORK, 00128 $ RESULT ) 00129 * 00130 * -- LAPACK test routine (version 3.4.0) -- 00131 * -- LAPACK is a software package provided by Univ. of Tennessee, -- 00132 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 00133 * November 2011 00134 * 00135 * .. Scalar Arguments .. 00136 INTEGER KBAND, LDU, N 00137 * .. 00138 * .. Array Arguments .. 00139 DOUBLE PRECISION AD( * ), AE( * ), RESULT( 2 ), SD( * ), 00140 $ SE( * ), U( LDU, * ), WORK( * ) 00141 * .. 00142 * 00143 * ===================================================================== 00144 * 00145 * .. Parameters .. 00146 DOUBLE PRECISION ZERO, ONE 00147 PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0 ) 00148 * .. 00149 * .. Local Scalars .. 00150 INTEGER J 00151 DOUBLE PRECISION ANORM, TEMP1, TEMP2, ULP, UNFL, WNORM 00152 * .. 00153 * .. External Functions .. 00154 DOUBLE PRECISION DLAMCH, DLANGE, DLANSY 00155 EXTERNAL DLAMCH, DLANGE, DLANSY 00156 * .. 00157 * .. External Subroutines .. 00158 EXTERNAL DGEMM, DLASET, DSYR, DSYR2 00159 * .. 00160 * .. Intrinsic Functions .. 00161 INTRINSIC ABS, DBLE, MAX, MIN 00162 * .. 00163 * .. Executable Statements .. 00164 * 00165 * 1) Constants 00166 * 00167 RESULT( 1 ) = ZERO 00168 RESULT( 2 ) = ZERO 00169 IF( N.LE.0 ) 00170 $ RETURN 00171 * 00172 UNFL = DLAMCH( 'Safe minimum' ) 00173 ULP = DLAMCH( 'Precision' ) 00174 * 00175 * Do Test 1 00176 * 00177 * Copy A & Compute its 1-Norm: 00178 * 00179 CALL DLASET( 'Full', N, N, ZERO, ZERO, WORK, N ) 00180 * 00181 ANORM = ZERO 00182 TEMP1 = ZERO 00183 * 00184 DO 10 J = 1, N - 1 00185 WORK( ( N+1 )*( J-1 )+1 ) = AD( J ) 00186 WORK( ( N+1 )*( J-1 )+2 ) = AE( J ) 00187 TEMP2 = ABS( AE( J ) ) 00188 ANORM = MAX( ANORM, ABS( AD( J ) )+TEMP1+TEMP2 ) 00189 TEMP1 = TEMP2 00190 10 CONTINUE 00191 * 00192 WORK( N**2 ) = AD( N ) 00193 ANORM = MAX( ANORM, ABS( AD( N ) )+TEMP1, UNFL ) 00194 * 00195 * Norm of A - USU' 00196 * 00197 DO 20 J = 1, N 00198 CALL DSYR( 'L', N, -SD( J ), U( 1, J ), 1, WORK, N ) 00199 20 CONTINUE 00200 * 00201 IF( N.GT.1 .AND. KBAND.EQ.1 ) THEN 00202 DO 30 J = 1, N - 1 00203 CALL DSYR2( 'L', N, -SE( J ), U( 1, J ), 1, U( 1, J+1 ), 1, 00204 $ WORK, N ) 00205 30 CONTINUE 00206 END IF 00207 * 00208 WNORM = DLANSY( '1', 'L', N, WORK, N, WORK( N**2+1 ) ) 00209 * 00210 IF( ANORM.GT.WNORM ) THEN 00211 RESULT( 1 ) = ( WNORM / ANORM ) / ( N*ULP ) 00212 ELSE 00213 IF( ANORM.LT.ONE ) THEN 00214 RESULT( 1 ) = ( MIN( WNORM, N*ANORM ) / ANORM ) / ( N*ULP ) 00215 ELSE 00216 RESULT( 1 ) = MIN( WNORM / ANORM, DBLE( N ) ) / ( N*ULP ) 00217 END IF 00218 END IF 00219 * 00220 * Do Test 2 00221 * 00222 * Compute UU' - I 00223 * 00224 CALL DGEMM( 'N', 'C', N, N, N, ONE, U, LDU, U, LDU, ZERO, WORK, 00225 $ N ) 00226 * 00227 DO 40 J = 1, N 00228 WORK( ( N+1 )*( J-1 )+1 ) = WORK( ( N+1 )*( J-1 )+1 ) - ONE 00229 40 CONTINUE 00230 * 00231 RESULT( 2 ) = MIN( DBLE( N ), DLANGE( '1', N, N, WORK, N, 00232 $ WORK( N**2+1 ) ) ) / ( N*ULP ) 00233 * 00234 RETURN 00235 * 00236 * End of DSTT21 00237 * 00238 END