![]() |
LAPACK
3.4.1
LAPACK: Linear Algebra PACKage
|
00001 *> \brief \b SSYT21 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 SSYT21( ITYPE, UPLO, N, KBAND, A, LDA, D, E, U, LDU, V, 00012 * LDV, TAU, WORK, RESULT ) 00013 * 00014 * .. Scalar Arguments .. 00015 * CHARACTER UPLO 00016 * INTEGER ITYPE, KBAND, LDA, LDU, LDV, N 00017 * .. 00018 * .. Array Arguments .. 00019 * REAL A( LDA, * ), D( * ), E( * ), RESULT( 2 ), 00020 * $ TAU( * ), U( LDU, * ), V( LDV, * ), WORK( * ) 00021 * .. 00022 * 00023 * 00024 *> \par Purpose: 00025 * ============= 00026 *> 00027 *> \verbatim 00028 *> 00029 *> SSYT21 generally checks a decomposition of the form 00030 *> 00031 *> A = U S U' 00032 *> 00033 *> where ' means transpose, A is symmetric, U is orthogonal, and S is 00034 *> diagonal (if KBAND=0) or symmetric tridiagonal (if KBAND=1). 00035 *> 00036 *> If ITYPE=1, then U is represented as a dense matrix; otherwise U is 00037 *> expressed as a product of Householder transformations, whose vectors 00038 *> are stored in the array "V" and whose scaling constants are in "TAU". 00039 *> We shall use the letter "V" to refer to the product of Householder 00040 *> transformations (which should be equal to U). 00041 *> 00042 *> Specifically, if ITYPE=1, then: 00043 *> 00044 *> RESULT(1) = | A - U S U' | / ( |A| n ulp ) *andC> RESULT(2) = | I - UU' | / ( n ulp ) 00045 *> 00046 *> If ITYPE=2, then: 00047 *> 00048 *> RESULT(1) = | A - V S V' | / ( |A| n ulp ) 00049 *> 00050 *> If ITYPE=3, then: 00051 *> 00052 *> RESULT(1) = | I - VU' | / ( n ulp ) 00053 *> 00054 *> For ITYPE > 1, the transformation U is expressed as a product 00055 *> V = H(1)...H(n-2), where H(j) = I - tau(j) v(j) v(j)' and each 00056 *> vector v(j) has its first j elements 0 and the remaining n-j elements 00057 *> stored in V(j+1:n,j). 00058 *> \endverbatim 00059 * 00060 * Arguments: 00061 * ========== 00062 * 00063 *> \param[in] ITYPE 00064 *> \verbatim 00065 *> ITYPE is INTEGER 00066 *> Specifies the type of tests to be performed. 00067 *> 1: U expressed as a dense orthogonal matrix: 00068 *> RESULT(1) = | A - U S U' | / ( |A| n ulp ) *andC> RESULT(2) = | I - UU' | / ( n ulp ) 00069 *> 00070 *> 2: U expressed as a product V of Housholder transformations: 00071 *> RESULT(1) = | A - V S V' | / ( |A| n ulp ) 00072 *> 00073 *> 3: U expressed both as a dense orthogonal matrix and 00074 *> as a product of Housholder transformations: 00075 *> RESULT(1) = | I - VU' | / ( n ulp ) 00076 *> \endverbatim 00077 *> 00078 *> \param[in] UPLO 00079 *> \verbatim 00080 *> UPLO is CHARACTER 00081 *> If UPLO='U', the upper triangle of A and V will be used and 00082 *> the (strictly) lower triangle will not be referenced. 00083 *> If UPLO='L', the lower triangle of A and V will be used and 00084 *> the (strictly) upper triangle will not be referenced. 00085 *> \endverbatim 00086 *> 00087 *> \param[in] N 00088 *> \verbatim 00089 *> N is INTEGER 00090 *> The size of the matrix. If it is zero, SSYT21 does nothing. 00091 *> It must be at least zero. 00092 *> \endverbatim 00093 *> 00094 *> \param[in] KBAND 00095 *> \verbatim 00096 *> KBAND is INTEGER 00097 *> The bandwidth of the matrix. It may only be zero or one. 00098 *> If zero, then S is diagonal, and E is not referenced. If 00099 *> one, then S is symmetric tri-diagonal. 00100 *> \endverbatim 00101 *> 00102 *> \param[in] A 00103 *> \verbatim 00104 *> A is REAL array, dimension (LDA, N) 00105 *> The original (unfactored) matrix. It is assumed to be 00106 *> symmetric, and only the upper (UPLO='U') or only the lower 00107 *> (UPLO='L') will be referenced. 00108 *> \endverbatim 00109 *> 00110 *> \param[in] LDA 00111 *> \verbatim 00112 *> LDA is INTEGER 00113 *> The leading dimension of A. It must be at least 1 00114 *> and at least N. 00115 *> \endverbatim 00116 *> 00117 *> \param[in] D 00118 *> \verbatim 00119 *> D is REAL array, dimension (N) 00120 *> The diagonal of the (symmetric tri-) diagonal matrix. 00121 *> \endverbatim 00122 *> 00123 *> \param[in] E 00124 *> \verbatim 00125 *> E is REAL array, dimension (N-1) 00126 *> The off-diagonal of the (symmetric tri-) diagonal matrix. 00127 *> E(1) is the (1,2) and (2,1) element, E(2) is the (2,3) and 00128 *> (3,2) element, etc. 00129 *> Not referenced if KBAND=0. 00130 *> \endverbatim 00131 *> 00132 *> \param[in] U 00133 *> \verbatim 00134 *> U is REAL array, dimension (LDU, N) 00135 *> If ITYPE=1 or 3, this contains the orthogonal matrix in 00136 *> the decomposition, expressed as a dense matrix. If ITYPE=2, 00137 *> then it is not referenced. 00138 *> \endverbatim 00139 *> 00140 *> \param[in] LDU 00141 *> \verbatim 00142 *> LDU is INTEGER 00143 *> The leading dimension of U. LDU must be at least N and 00144 *> at least 1. 00145 *> \endverbatim 00146 *> 00147 *> \param[in] V 00148 *> \verbatim 00149 *> V is REAL array, dimension (LDV, N) 00150 *> If ITYPE=2 or 3, the columns of this array contain the 00151 *> Householder vectors used to describe the orthogonal matrix 00152 *> in the decomposition. If UPLO='L', then the vectors are in 00153 *> the lower triangle, if UPLO='U', then in the upper 00154 *> triangle. 00155 *> *NOTE* If ITYPE=2 or 3, V is modified and restored. The 00156 *> subdiagonal (if UPLO='L') or the superdiagonal (if UPLO='U') 00157 *> is set to one, and later reset to its original value, during 00158 *> the course of the calculation. 00159 *> If ITYPE=1, then it is neither referenced nor modified. 00160 *> \endverbatim 00161 *> 00162 *> \param[in] LDV 00163 *> \verbatim 00164 *> LDV is INTEGER 00165 *> The leading dimension of V. LDV must be at least N and 00166 *> at least 1. 00167 *> \endverbatim 00168 *> 00169 *> \param[in] TAU 00170 *> \verbatim 00171 *> TAU is REAL array, dimension (N) 00172 *> If ITYPE >= 2, then TAU(j) is the scalar factor of 00173 *> v(j) v(j)' in the Householder transformation H(j) of 00174 *> the product U = H(1)...H(n-2) 00175 *> If ITYPE < 2, then TAU is not referenced. 00176 *> \endverbatim 00177 *> 00178 *> \param[out] WORK 00179 *> \verbatim 00180 *> WORK is REAL array, dimension (2*N**2) 00181 *> \endverbatim 00182 *> 00183 *> \param[out] RESULT 00184 *> \verbatim 00185 *> RESULT is REAL array, dimension (2) 00186 *> The values computed by the two tests described above. The 00187 *> values are currently limited to 1/ulp, to avoid overflow. 00188 *> RESULT(1) is always modified. RESULT(2) is modified only 00189 *> if ITYPE=1. 00190 *> \endverbatim 00191 * 00192 * Authors: 00193 * ======== 00194 * 00195 *> \author Univ. of Tennessee 00196 *> \author Univ. of California Berkeley 00197 *> \author Univ. of Colorado Denver 00198 *> \author NAG Ltd. 00199 * 00200 *> \date November 2011 00201 * 00202 *> \ingroup single_eig 00203 * 00204 * ===================================================================== 00205 SUBROUTINE SSYT21( ITYPE, UPLO, N, KBAND, A, LDA, D, E, U, LDU, V, 00206 $ LDV, TAU, WORK, RESULT ) 00207 * 00208 * -- LAPACK test routine (version 3.4.0) -- 00209 * -- LAPACK is a software package provided by Univ. of Tennessee, -- 00210 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 00211 * November 2011 00212 * 00213 * .. Scalar Arguments .. 00214 CHARACTER UPLO 00215 INTEGER ITYPE, KBAND, LDA, LDU, LDV, N 00216 * .. 00217 * .. Array Arguments .. 00218 REAL A( LDA, * ), D( * ), E( * ), RESULT( 2 ), 00219 $ TAU( * ), U( LDU, * ), V( LDV, * ), WORK( * ) 00220 * .. 00221 * 00222 * ===================================================================== 00223 * 00224 * .. Parameters .. 00225 REAL ZERO, ONE, TEN 00226 PARAMETER ( ZERO = 0.0E0, ONE = 1.0E0, TEN = 10.0E0 ) 00227 * .. 00228 * .. Local Scalars .. 00229 LOGICAL LOWER 00230 CHARACTER CUPLO 00231 INTEGER IINFO, J, JCOL, JR, JROW 00232 REAL ANORM, ULP, UNFL, VSAVE, WNORM 00233 * .. 00234 * .. External Functions .. 00235 LOGICAL LSAME 00236 REAL SLAMCH, SLANGE, SLANSY 00237 EXTERNAL LSAME, SLAMCH, SLANGE, SLANSY 00238 * .. 00239 * .. External Subroutines .. 00240 EXTERNAL SGEMM, SLACPY, SLARFY, SLASET, SORM2L, SORM2R, 00241 $ SSYR, SSYR2 00242 * .. 00243 * .. Intrinsic Functions .. 00244 INTRINSIC MAX, MIN, REAL 00245 * .. 00246 * .. Executable Statements .. 00247 * 00248 RESULT( 1 ) = ZERO 00249 IF( ITYPE.EQ.1 ) 00250 $ RESULT( 2 ) = ZERO 00251 IF( N.LE.0 ) 00252 $ RETURN 00253 * 00254 IF( LSAME( UPLO, 'U' ) ) THEN 00255 LOWER = .FALSE. 00256 CUPLO = 'U' 00257 ELSE 00258 LOWER = .TRUE. 00259 CUPLO = 'L' 00260 END IF 00261 * 00262 UNFL = SLAMCH( 'Safe minimum' ) 00263 ULP = SLAMCH( 'Epsilon' )*SLAMCH( 'Base' ) 00264 * 00265 * Some Error Checks 00266 * 00267 IF( ITYPE.LT.1 .OR. ITYPE.GT.3 ) THEN 00268 RESULT( 1 ) = TEN / ULP 00269 RETURN 00270 END IF 00271 * 00272 * Do Test 1 00273 * 00274 * Norm of A: 00275 * 00276 IF( ITYPE.EQ.3 ) THEN 00277 ANORM = ONE 00278 ELSE 00279 ANORM = MAX( SLANSY( '1', CUPLO, N, A, LDA, WORK ), UNFL ) 00280 END IF 00281 * 00282 * Compute error matrix: 00283 * 00284 IF( ITYPE.EQ.1 ) THEN 00285 * 00286 * ITYPE=1: error = A - U S U' 00287 * 00288 CALL SLASET( 'Full', N, N, ZERO, ZERO, WORK, N ) 00289 CALL SLACPY( CUPLO, N, N, A, LDA, WORK, N ) 00290 * 00291 DO 10 J = 1, N 00292 CALL SSYR( CUPLO, N, -D( J ), U( 1, J ), 1, WORK, N ) 00293 10 CONTINUE 00294 * 00295 IF( N.GT.1 .AND. KBAND.EQ.1 ) THEN 00296 DO 20 J = 1, N - 1 00297 CALL SSYR2( CUPLO, N, -E( J ), U( 1, J ), 1, U( 1, J+1 ), 00298 $ 1, WORK, N ) 00299 20 CONTINUE 00300 END IF 00301 WNORM = SLANSY( '1', CUPLO, N, WORK, N, WORK( N**2+1 ) ) 00302 * 00303 ELSE IF( ITYPE.EQ.2 ) THEN 00304 * 00305 * ITYPE=2: error = V S V' - A 00306 * 00307 CALL SLASET( 'Full', N, N, ZERO, ZERO, WORK, N ) 00308 * 00309 IF( LOWER ) THEN 00310 WORK( N**2 ) = D( N ) 00311 DO 40 J = N - 1, 1, -1 00312 IF( KBAND.EQ.1 ) THEN 00313 WORK( ( N+1 )*( J-1 )+2 ) = ( ONE-TAU( J ) )*E( J ) 00314 DO 30 JR = J + 2, N 00315 WORK( ( J-1 )*N+JR ) = -TAU( J )*E( J )*V( JR, J ) 00316 30 CONTINUE 00317 END IF 00318 * 00319 VSAVE = V( J+1, J ) 00320 V( J+1, J ) = ONE 00321 CALL SLARFY( 'L', N-J, V( J+1, J ), 1, TAU( J ), 00322 $ WORK( ( N+1 )*J+1 ), N, WORK( N**2+1 ) ) 00323 V( J+1, J ) = VSAVE 00324 WORK( ( N+1 )*( J-1 )+1 ) = D( J ) 00325 40 CONTINUE 00326 ELSE 00327 WORK( 1 ) = D( 1 ) 00328 DO 60 J = 1, N - 1 00329 IF( KBAND.EQ.1 ) THEN 00330 WORK( ( N+1 )*J ) = ( ONE-TAU( J ) )*E( J ) 00331 DO 50 JR = 1, J - 1 00332 WORK( J*N+JR ) = -TAU( J )*E( J )*V( JR, J+1 ) 00333 50 CONTINUE 00334 END IF 00335 * 00336 VSAVE = V( J, J+1 ) 00337 V( J, J+1 ) = ONE 00338 CALL SLARFY( 'U', J, V( 1, J+1 ), 1, TAU( J ), WORK, N, 00339 $ WORK( N**2+1 ) ) 00340 V( J, J+1 ) = VSAVE 00341 WORK( ( N+1 )*J+1 ) = D( J+1 ) 00342 60 CONTINUE 00343 END IF 00344 * 00345 DO 90 JCOL = 1, N 00346 IF( LOWER ) THEN 00347 DO 70 JROW = JCOL, N 00348 WORK( JROW+N*( JCOL-1 ) ) = WORK( JROW+N*( JCOL-1 ) ) 00349 $ - A( JROW, JCOL ) 00350 70 CONTINUE 00351 ELSE 00352 DO 80 JROW = 1, JCOL 00353 WORK( JROW+N*( JCOL-1 ) ) = WORK( JROW+N*( JCOL-1 ) ) 00354 $ - A( JROW, JCOL ) 00355 80 CONTINUE 00356 END IF 00357 90 CONTINUE 00358 WNORM = SLANSY( '1', CUPLO, N, WORK, N, WORK( N**2+1 ) ) 00359 * 00360 ELSE IF( ITYPE.EQ.3 ) THEN 00361 * 00362 * ITYPE=3: error = U V' - I 00363 * 00364 IF( N.LT.2 ) 00365 $ RETURN 00366 CALL SLACPY( ' ', N, N, U, LDU, WORK, N ) 00367 IF( LOWER ) THEN 00368 CALL SORM2R( 'R', 'T', N, N-1, N-1, V( 2, 1 ), LDV, TAU, 00369 $ WORK( N+1 ), N, WORK( N**2+1 ), IINFO ) 00370 ELSE 00371 CALL SORM2L( 'R', 'T', N, N-1, N-1, V( 1, 2 ), LDV, TAU, 00372 $ WORK, N, WORK( N**2+1 ), IINFO ) 00373 END IF 00374 IF( IINFO.NE.0 ) THEN 00375 RESULT( 1 ) = TEN / ULP 00376 RETURN 00377 END IF 00378 * 00379 DO 100 J = 1, N 00380 WORK( ( N+1 )*( J-1 )+1 ) = WORK( ( N+1 )*( J-1 )+1 ) - ONE 00381 100 CONTINUE 00382 * 00383 WNORM = SLANGE( '1', N, N, WORK, N, WORK( N**2+1 ) ) 00384 END IF 00385 * 00386 IF( ANORM.GT.WNORM ) THEN 00387 RESULT( 1 ) = ( WNORM / ANORM ) / ( N*ULP ) 00388 ELSE 00389 IF( ANORM.LT.ONE ) THEN 00390 RESULT( 1 ) = ( MIN( WNORM, N*ANORM ) / ANORM ) / ( N*ULP ) 00391 ELSE 00392 RESULT( 1 ) = MIN( WNORM / ANORM, REAL( N ) ) / ( N*ULP ) 00393 END IF 00394 END IF 00395 * 00396 * Do Test 2 00397 * 00398 * Compute UU' - I 00399 * 00400 IF( ITYPE.EQ.1 ) THEN 00401 CALL SGEMM( 'N', 'C', N, N, N, ONE, U, LDU, U, LDU, ZERO, WORK, 00402 $ N ) 00403 * 00404 DO 110 J = 1, N 00405 WORK( ( N+1 )*( J-1 )+1 ) = WORK( ( N+1 )*( J-1 )+1 ) - ONE 00406 110 CONTINUE 00407 * 00408 RESULT( 2 ) = MIN( SLANGE( '1', N, N, WORK, N, 00409 $ WORK( N**2+1 ) ), REAL( N ) ) / ( N*ULP ) 00410 END IF 00411 * 00412 RETURN 00413 * 00414 * End of SSYT21 00415 * 00416 END