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