![]() |
LAPACK
3.4.1
LAPACK: Linear Algebra PACKage
|
00001 *> \brief \b SSBT21 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 SSBT21( UPLO, N, KA, KS, A, LDA, D, E, U, LDU, WORK, 00012 * RESULT ) 00013 * 00014 * .. Scalar Arguments .. 00015 * CHARACTER UPLO 00016 * INTEGER KA, KS, LDA, LDU, N 00017 * .. 00018 * .. Array Arguments .. 00019 * REAL A( LDA, * ), D( * ), E( * ), RESULT( 2 ), 00020 * $ U( LDU, * ), WORK( * ) 00021 * .. 00022 * 00023 * 00024 *> \par Purpose: 00025 * ============= 00026 *> 00027 *> \verbatim 00028 *> 00029 *> SSBT21 generally checks a decomposition of the form 00030 *> 00031 *> A = U S U' 00032 *> 00033 *> where ' means transpose, A is symmetric banded, U is 00034 *> orthogonal, and S is diagonal (if KS=0) or symmetric 00035 *> tridiagonal (if KS=1). 00036 *> 00037 *> Specifically: 00038 *> 00039 *> RESULT(1) = | A - U S U' | / ( |A| n ulp ) *andC> RESULT(2) = | I - UU' | / ( n ulp ) 00040 *> \endverbatim 00041 * 00042 * Arguments: 00043 * ========== 00044 * 00045 *> \param[in] UPLO 00046 *> \verbatim 00047 *> UPLO is CHARACTER 00048 *> If UPLO='U', the upper triangle of A and V will be used and 00049 *> the (strictly) lower triangle will not be referenced. 00050 *> If UPLO='L', the lower triangle of A and V will be used and 00051 *> the (strictly) upper triangle will not be referenced. 00052 *> \endverbatim 00053 *> 00054 *> \param[in] N 00055 *> \verbatim 00056 *> N is INTEGER 00057 *> The size of the matrix. If it is zero, SSBT21 does nothing. 00058 *> It must be at least zero. 00059 *> \endverbatim 00060 *> 00061 *> \param[in] KA 00062 *> \verbatim 00063 *> KA is INTEGER 00064 *> The bandwidth of the matrix A. It must be at least zero. If 00065 *> it is larger than N-1, then max( 0, N-1 ) will be used. 00066 *> \endverbatim 00067 *> 00068 *> \param[in] KS 00069 *> \verbatim 00070 *> KS is INTEGER 00071 *> The bandwidth of the matrix S. It may only be zero or one. 00072 *> If zero, then S is diagonal, and E is not referenced. If 00073 *> one, then S is symmetric tri-diagonal. 00074 *> \endverbatim 00075 *> 00076 *> \param[in] A 00077 *> \verbatim 00078 *> A is REAL array, dimension (LDA, N) 00079 *> The original (unfactored) matrix. It is assumed to be 00080 *> symmetric, and only the upper (UPLO='U') or only the lower 00081 *> (UPLO='L') will be referenced. 00082 *> \endverbatim 00083 *> 00084 *> \param[in] LDA 00085 *> \verbatim 00086 *> LDA is INTEGER 00087 *> The leading dimension of A. It must be at least 1 00088 *> and at least min( KA, N-1 ). 00089 *> \endverbatim 00090 *> 00091 *> \param[in] D 00092 *> \verbatim 00093 *> D is REAL array, dimension (N) 00094 *> The diagonal of the (symmetric tri-) diagonal matrix S. 00095 *> \endverbatim 00096 *> 00097 *> \param[in] E 00098 *> \verbatim 00099 *> E is REAL array, dimension (N-1) 00100 *> The off-diagonal of the (symmetric tri-) diagonal matrix S. 00101 *> E(1) is the (1,2) and (2,1) element, E(2) is the (2,3) and 00102 *> (3,2) element, etc. 00103 *> Not referenced if KS=0. 00104 *> \endverbatim 00105 *> 00106 *> \param[in] U 00107 *> \verbatim 00108 *> U is REAL array, dimension (LDU, N) 00109 *> The orthogonal matrix in the decomposition, expressed as a 00110 *> dense matrix (i.e., not as a product of Householder 00111 *> transformations, Givens transformations, etc.) 00112 *> \endverbatim 00113 *> 00114 *> \param[in] LDU 00115 *> \verbatim 00116 *> LDU is INTEGER 00117 *> The leading dimension of U. LDU must be at least N and 00118 *> at least 1. 00119 *> \endverbatim 00120 *> 00121 *> \param[out] WORK 00122 *> \verbatim 00123 *> WORK is REAL array, dimension (N**2+N) 00124 *> \endverbatim 00125 *> 00126 *> \param[out] RESULT 00127 *> \verbatim 00128 *> RESULT is REAL array, dimension (2) 00129 *> The values computed by the two tests described above. The 00130 *> values are currently limited to 1/ulp, to avoid overflow. 00131 *> \endverbatim 00132 * 00133 * Authors: 00134 * ======== 00135 * 00136 *> \author Univ. of Tennessee 00137 *> \author Univ. of California Berkeley 00138 *> \author Univ. of Colorado Denver 00139 *> \author NAG Ltd. 00140 * 00141 *> \date November 2011 00142 * 00143 *> \ingroup single_eig 00144 * 00145 * ===================================================================== 00146 SUBROUTINE SSBT21( UPLO, N, KA, KS, A, LDA, D, E, U, LDU, WORK, 00147 $ RESULT ) 00148 * 00149 * -- LAPACK test routine (version 3.4.0) -- 00150 * -- LAPACK is a software package provided by Univ. of Tennessee, -- 00151 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 00152 * November 2011 00153 * 00154 * .. Scalar Arguments .. 00155 CHARACTER UPLO 00156 INTEGER KA, KS, LDA, LDU, N 00157 * .. 00158 * .. Array Arguments .. 00159 REAL A( LDA, * ), D( * ), E( * ), RESULT( 2 ), 00160 $ U( LDU, * ), WORK( * ) 00161 * .. 00162 * 00163 * ===================================================================== 00164 * 00165 * .. Parameters .. 00166 REAL ZERO, ONE 00167 PARAMETER ( ZERO = 0.0E0, ONE = 1.0E0 ) 00168 * .. 00169 * .. Local Scalars .. 00170 LOGICAL LOWER 00171 CHARACTER CUPLO 00172 INTEGER IKA, J, JC, JR, LW 00173 REAL ANORM, ULP, UNFL, WNORM 00174 * .. 00175 * .. External Functions .. 00176 LOGICAL LSAME 00177 REAL SLAMCH, SLANGE, SLANSB, SLANSP 00178 EXTERNAL LSAME, SLAMCH, SLANGE, SLANSB, SLANSP 00179 * .. 00180 * .. External Subroutines .. 00181 EXTERNAL SGEMM, SSPR, SSPR2 00182 * .. 00183 * .. Intrinsic Functions .. 00184 INTRINSIC MAX, MIN, REAL 00185 * .. 00186 * .. Executable Statements .. 00187 * 00188 * Constants 00189 * 00190 RESULT( 1 ) = ZERO 00191 RESULT( 2 ) = ZERO 00192 IF( N.LE.0 ) 00193 $ RETURN 00194 * 00195 IKA = MAX( 0, MIN( N-1, KA ) ) 00196 LW = ( N*( N+1 ) ) / 2 00197 * 00198 IF( LSAME( UPLO, 'U' ) ) THEN 00199 LOWER = .FALSE. 00200 CUPLO = 'U' 00201 ELSE 00202 LOWER = .TRUE. 00203 CUPLO = 'L' 00204 END IF 00205 * 00206 UNFL = SLAMCH( 'Safe minimum' ) 00207 ULP = SLAMCH( 'Epsilon' )*SLAMCH( 'Base' ) 00208 * 00209 * Some Error Checks 00210 * 00211 * Do Test 1 00212 * 00213 * Norm of A: 00214 * 00215 ANORM = MAX( SLANSB( '1', CUPLO, N, IKA, A, LDA, WORK ), UNFL ) 00216 * 00217 * Compute error matrix: Error = A - U S U' 00218 * 00219 * Copy A from SB to SP storage format. 00220 * 00221 J = 0 00222 DO 50 JC = 1, N 00223 IF( LOWER ) THEN 00224 DO 10 JR = 1, MIN( IKA+1, N+1-JC ) 00225 J = J + 1 00226 WORK( J ) = A( JR, JC ) 00227 10 CONTINUE 00228 DO 20 JR = IKA + 2, N + 1 - JC 00229 J = J + 1 00230 WORK( J ) = ZERO 00231 20 CONTINUE 00232 ELSE 00233 DO 30 JR = IKA + 2, JC 00234 J = J + 1 00235 WORK( J ) = ZERO 00236 30 CONTINUE 00237 DO 40 JR = MIN( IKA, JC-1 ), 0, -1 00238 J = J + 1 00239 WORK( J ) = A( IKA+1-JR, JC ) 00240 40 CONTINUE 00241 END IF 00242 50 CONTINUE 00243 * 00244 DO 60 J = 1, N 00245 CALL SSPR( CUPLO, N, -D( J ), U( 1, J ), 1, WORK ) 00246 60 CONTINUE 00247 * 00248 IF( N.GT.1 .AND. KS.EQ.1 ) THEN 00249 DO 70 J = 1, N - 1 00250 CALL SSPR2( CUPLO, N, -E( J ), U( 1, J ), 1, U( 1, J+1 ), 1, 00251 $ WORK ) 00252 70 CONTINUE 00253 END IF 00254 WNORM = SLANSP( '1', CUPLO, N, WORK, WORK( LW+1 ) ) 00255 * 00256 IF( ANORM.GT.WNORM ) THEN 00257 RESULT( 1 ) = ( WNORM / ANORM ) / ( N*ULP ) 00258 ELSE 00259 IF( ANORM.LT.ONE ) THEN 00260 RESULT( 1 ) = ( MIN( WNORM, N*ANORM ) / ANORM ) / ( N*ULP ) 00261 ELSE 00262 RESULT( 1 ) = MIN( WNORM / ANORM, REAL( N ) ) / ( N*ULP ) 00263 END IF 00264 END IF 00265 * 00266 * Do Test 2 00267 * 00268 * Compute UU' - I 00269 * 00270 CALL SGEMM( 'N', 'C', N, N, N, ONE, U, LDU, U, LDU, ZERO, WORK, 00271 $ N ) 00272 * 00273 DO 80 J = 1, N 00274 WORK( ( N+1 )*( J-1 )+1 ) = WORK( ( N+1 )*( J-1 )+1 ) - ONE 00275 80 CONTINUE 00276 * 00277 RESULT( 2 ) = MIN( SLANGE( '1', N, N, WORK, N, WORK( N**2+1 ) ), 00278 $ REAL( N ) ) / ( N*ULP ) 00279 * 00280 RETURN 00281 * 00282 * End of SSBT21 00283 * 00284 END