 |
LAPACK 3.3.1
Linear Algebra PACKage
|
|
LAPACK 3.3.1
Linear Algebra PACKage
|
00001 SUBROUTINE DSYT21( ITYPE, UPLO, N, KBAND, A, LDA, D, E, U, LDU, V, 00002 $ LDV, TAU, WORK, RESULT ) 00003 * 00004 * -- LAPACK test routine (version 3.1) -- 00005 * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. 00006 * November 2006 00007 * 00008 * .. Scalar Arguments .. 00009 CHARACTER UPLO 00010 INTEGER ITYPE, KBAND, LDA, LDU, LDV, N 00011 * .. 00012 * .. Array Arguments .. 00013 DOUBLE PRECISION A( LDA, * ), D( * ), E( * ), RESULT( 2 ), 00014 $ TAU( * ), U( LDU, * ), V( LDV, * ), WORK( * ) 00015 * .. 00016 * 00017 * Purpose 00018 * ======= 00019 * 00020 * DSYT21 generally checks a decomposition of the form 00021 * 00022 * A = U S U' 00023 * 00024 * where ' means transpose, A is symmetric, U is orthogonal, and S is 00025 * diagonal (if KBAND=0) or symmetric tridiagonal (if KBAND=1). 00026 * 00027 * If ITYPE=1, then U is represented as a dense matrix; otherwise U is 00028 * expressed as a product of Householder transformations, whose vectors 00029 * are stored in the array "V" and whose scaling constants are in "TAU". 00030 * We shall use the letter "V" to refer to the product of Householder 00031 * transformations (which should be equal to U). 00032 * 00033 * Specifically, if ITYPE=1, then: 00034 * 00035 * RESULT(1) = | A - U S U' | / ( |A| n ulp ) *and* 00036 * RESULT(2) = | I - UU' | / ( n ulp ) 00037 * 00038 * If ITYPE=2, then: 00039 * 00040 * RESULT(1) = | A - V S V' | / ( |A| n ulp ) 00041 * 00042 * If ITYPE=3, then: 00043 * 00044 * RESULT(1) = | I - VU' | / ( n ulp ) 00045 * 00046 * For ITYPE > 1, the transformation U is expressed as a product 00047 * V = H(1)...H(n-2), where H(j) = I - tau(j) v(j) v(j)' and each 00048 * vector v(j) has its first j elements 0 and the remaining n-j elements 00049 * stored in V(j+1:n,j). 00050 * 00051 * Arguments 00052 * ========= 00053 * 00054 * ITYPE (input) INTEGER 00055 * Specifies the type of tests to be performed. 00056 * 1: U expressed as a dense orthogonal matrix: 00057 * RESULT(1) = | A - U S U' | / ( |A| n ulp ) *and* 00058 * RESULT(2) = | I - UU' | / ( n ulp ) 00059 * 00060 * 2: U expressed as a product V of Housholder transformations: 00061 * RESULT(1) = | A - V S V' | / ( |A| n ulp ) 00062 * 00063 * 3: U expressed both as a dense orthogonal matrix and 00064 * as a product of Housholder transformations: 00065 * RESULT(1) = | I - VU' | / ( n ulp ) 00066 * 00067 * UPLO (input) CHARACTER 00068 * If UPLO='U', the upper triangle of A and V will be used and 00069 * the (strictly) lower triangle will not be referenced. 00070 * If UPLO='L', the lower triangle of A and V will be used and 00071 * the (strictly) upper triangle will not be referenced. 00072 * 00073 * N (input) INTEGER 00074 * The size of the matrix. If it is zero, DSYT21 does nothing. 00075 * It must be at least zero. 00076 * 00077 * KBAND (input) INTEGER 00078 * The bandwidth of the matrix. It may only be zero or one. 00079 * If zero, then S is diagonal, and E is not referenced. If 00080 * one, then S is symmetric tri-diagonal. 00081 * 00082 * A (input) DOUBLE PRECISION array, dimension (LDA, N) 00083 * The original (unfactored) matrix. It is assumed to be 00084 * symmetric, and only the upper (UPLO='U') or only the lower 00085 * (UPLO='L') will be referenced. 00086 * 00087 * LDA (input) INTEGER 00088 * The leading dimension of A. It must be at least 1 00089 * and at least N. 00090 * 00091 * D (input) DOUBLE PRECISION array, dimension (N) 00092 * The diagonal of the (symmetric tri-) diagonal matrix. 00093 * 00094 * E (input) DOUBLE PRECISION array, dimension (N-1) 00095 * The off-diagonal of the (symmetric tri-) diagonal matrix. 00096 * E(1) is the (1,2) and (2,1) element, E(2) is the (2,3) and 00097 * (3,2) element, etc. 00098 * Not referenced if KBAND=0. 00099 * 00100 * U (input) DOUBLE PRECISION array, dimension (LDU, N) 00101 * If ITYPE=1 or 3, this contains the orthogonal matrix in 00102 * the decomposition, expressed as a dense matrix. If ITYPE=2, 00103 * then it is not referenced. 00104 * 00105 * LDU (input) INTEGER 00106 * The leading dimension of U. LDU must be at least N and 00107 * at least 1. 00108 * 00109 * V (input) DOUBLE PRECISION array, dimension (LDV, N) 00110 * If ITYPE=2 or 3, the columns of this array contain the 00111 * Householder vectors used to describe the orthogonal matrix 00112 * in the decomposition. If UPLO='L', then the vectors are in 00113 * the lower triangle, if UPLO='U', then in the upper 00114 * triangle. 00115 * *NOTE* If ITYPE=2 or 3, V is modified and restored. The 00116 * subdiagonal (if UPLO='L') or the superdiagonal (if UPLO='U') 00117 * is set to one, and later reset to its original value, during 00118 * the course of the calculation. 00119 * If ITYPE=1, then it is neither referenced nor modified. 00120 * 00121 * LDV (input) INTEGER 00122 * The leading dimension of V. LDV must be at least N and 00123 * at least 1. 00124 * 00125 * TAU (input) DOUBLE PRECISION array, dimension (N) 00126 * If ITYPE >= 2, then TAU(j) is the scalar factor of 00127 * v(j) v(j)' in the Householder transformation H(j) of 00128 * the product U = H(1)...H(n-2) 00129 * If ITYPE < 2, then TAU is not referenced. 00130 * 00131 * WORK (workspace) DOUBLE PRECISION array, dimension (2*N**2) 00132 * 00133 * RESULT (output) DOUBLE PRECISION array, dimension (2) 00134 * The values computed by the two tests described above. The 00135 * values are currently limited to 1/ulp, to avoid overflow. 00136 * RESULT(1) is always modified. RESULT(2) is modified only 00137 * if ITYPE=1. 00138 * 00139 * ===================================================================== 00140 * 00141 * .. Parameters .. 00142 DOUBLE PRECISION ZERO, ONE, TEN 00143 PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0, TEN = 10.0D0 ) 00144 * .. 00145 * .. Local Scalars .. 00146 LOGICAL LOWER 00147 CHARACTER CUPLO 00148 INTEGER IINFO, J, JCOL, JR, JROW 00149 DOUBLE PRECISION ANORM, ULP, UNFL, VSAVE, WNORM 00150 * .. 00151 * .. External Functions .. 00152 LOGICAL LSAME 00153 DOUBLE PRECISION DLAMCH, DLANGE, DLANSY 00154 EXTERNAL LSAME, DLAMCH, DLANGE, DLANSY 00155 * .. 00156 * .. External Subroutines .. 00157 EXTERNAL DGEMM, DLACPY, DLARFY, DLASET, DORM2L, DORM2R, 00158 $ DSYR, DSYR2 00159 * .. 00160 * .. Intrinsic Functions .. 00161 INTRINSIC DBLE, MAX, MIN 00162 * .. 00163 * .. Executable Statements .. 00164 * 00165 RESULT( 1 ) = ZERO 00166 IF( ITYPE.EQ.1 ) 00167 $ RESULT( 2 ) = ZERO 00168 IF( N.LE.0 ) 00169 $ RETURN 00170 * 00171 IF( LSAME( UPLO, 'U' ) ) THEN 00172 LOWER = .FALSE. 00173 CUPLO = 'U' 00174 ELSE 00175 LOWER = .TRUE. 00176 CUPLO = 'L' 00177 END IF 00178 * 00179 UNFL = DLAMCH( 'Safe minimum' ) 00180 ULP = DLAMCH( 'Epsilon' )*DLAMCH( 'Base' ) 00181 * 00182 * Some Error Checks 00183 * 00184 IF( ITYPE.LT.1 .OR. ITYPE.GT.3 ) THEN 00185 RESULT( 1 ) = TEN / ULP 00186 RETURN 00187 END IF 00188 * 00189 * Do Test 1 00190 * 00191 * Norm of A: 00192 * 00193 IF( ITYPE.EQ.3 ) THEN 00194 ANORM = ONE 00195 ELSE 00196 ANORM = MAX( DLANSY( '1', CUPLO, N, A, LDA, WORK ), UNFL ) 00197 END IF 00198 * 00199 * Compute error matrix: 00200 * 00201 IF( ITYPE.EQ.1 ) THEN 00202 * 00203 * ITYPE=1: error = A - U S U' 00204 * 00205 CALL DLASET( 'Full', N, N, ZERO, ZERO, WORK, N ) 00206 CALL DLACPY( CUPLO, N, N, A, LDA, WORK, N ) 00207 * 00208 DO 10 J = 1, N 00209 CALL DSYR( CUPLO, N, -D( J ), U( 1, J ), 1, WORK, N ) 00210 10 CONTINUE 00211 * 00212 IF( N.GT.1 .AND. KBAND.EQ.1 ) THEN 00213 DO 20 J = 1, N - 1 00214 CALL DSYR2( CUPLO, N, -E( J ), U( 1, J ), 1, U( 1, J+1 ), 00215 $ 1, WORK, N ) 00216 20 CONTINUE 00217 END IF 00218 WNORM = DLANSY( '1', CUPLO, N, WORK, N, WORK( N**2+1 ) ) 00219 * 00220 ELSE IF( ITYPE.EQ.2 ) THEN 00221 * 00222 * ITYPE=2: error = V S V' - A 00223 * 00224 CALL DLASET( 'Full', N, N, ZERO, ZERO, WORK, N ) 00225 * 00226 IF( LOWER ) THEN 00227 WORK( N**2 ) = D( N ) 00228 DO 40 J = N - 1, 1, -1 00229 IF( KBAND.EQ.1 ) THEN 00230 WORK( ( N+1 )*( J-1 )+2 ) = ( ONE-TAU( J ) )*E( J ) 00231 DO 30 JR = J + 2, N 00232 WORK( ( J-1 )*N+JR ) = -TAU( J )*E( J )*V( JR, J ) 00233 30 CONTINUE 00234 END IF 00235 * 00236 VSAVE = V( J+1, J ) 00237 V( J+1, J ) = ONE 00238 CALL DLARFY( 'L', N-J, V( J+1, J ), 1, TAU( J ), 00239 $ WORK( ( N+1 )*J+1 ), N, WORK( N**2+1 ) ) 00240 V( J+1, J ) = VSAVE 00241 WORK( ( N+1 )*( J-1 )+1 ) = D( J ) 00242 40 CONTINUE 00243 ELSE 00244 WORK( 1 ) = D( 1 ) 00245 DO 60 J = 1, N - 1 00246 IF( KBAND.EQ.1 ) THEN 00247 WORK( ( N+1 )*J ) = ( ONE-TAU( J ) )*E( J ) 00248 DO 50 JR = 1, J - 1 00249 WORK( J*N+JR ) = -TAU( J )*E( J )*V( JR, J+1 ) 00250 50 CONTINUE 00251 END IF 00252 * 00253 VSAVE = V( J, J+1 ) 00254 V( J, J+1 ) = ONE 00255 CALL DLARFY( 'U', J, V( 1, J+1 ), 1, TAU( J ), WORK, N, 00256 $ WORK( N**2+1 ) ) 00257 V( J, J+1 ) = VSAVE 00258 WORK( ( N+1 )*J+1 ) = D( J+1 ) 00259 60 CONTINUE 00260 END IF 00261 * 00262 DO 90 JCOL = 1, N 00263 IF( LOWER ) THEN 00264 DO 70 JROW = JCOL, N 00265 WORK( JROW+N*( JCOL-1 ) ) = WORK( JROW+N*( JCOL-1 ) ) 00266 $ - A( JROW, JCOL ) 00267 70 CONTINUE 00268 ELSE 00269 DO 80 JROW = 1, JCOL 00270 WORK( JROW+N*( JCOL-1 ) ) = WORK( JROW+N*( JCOL-1 ) ) 00271 $ - A( JROW, JCOL ) 00272 80 CONTINUE 00273 END IF 00274 90 CONTINUE 00275 WNORM = DLANSY( '1', CUPLO, N, WORK, N, WORK( N**2+1 ) ) 00276 * 00277 ELSE IF( ITYPE.EQ.3 ) THEN 00278 * 00279 * ITYPE=3: error = U V' - I 00280 * 00281 IF( N.LT.2 ) 00282 $ RETURN 00283 CALL DLACPY( ' ', N, N, U, LDU, WORK, N ) 00284 IF( LOWER ) THEN 00285 CALL DORM2R( 'R', 'T', N, N-1, N-1, V( 2, 1 ), LDV, TAU, 00286 $ WORK( N+1 ), N, WORK( N**2+1 ), IINFO ) 00287 ELSE 00288 CALL DORM2L( 'R', 'T', N, N-1, N-1, V( 1, 2 ), LDV, TAU, 00289 $ WORK, N, WORK( N**2+1 ), IINFO ) 00290 END IF 00291 IF( IINFO.NE.0 ) THEN 00292 RESULT( 1 ) = TEN / ULP 00293 RETURN 00294 END IF 00295 * 00296 DO 100 J = 1, N 00297 WORK( ( N+1 )*( J-1 )+1 ) = WORK( ( N+1 )*( J-1 )+1 ) - ONE 00298 100 CONTINUE 00299 * 00300 WNORM = DLANGE( '1', N, N, WORK, N, WORK( N**2+1 ) ) 00301 END IF 00302 * 00303 IF( ANORM.GT.WNORM ) THEN 00304 RESULT( 1 ) = ( WNORM / ANORM ) / ( N*ULP ) 00305 ELSE 00306 IF( ANORM.LT.ONE ) THEN 00307 RESULT( 1 ) = ( MIN( WNORM, N*ANORM ) / ANORM ) / ( N*ULP ) 00308 ELSE 00309 RESULT( 1 ) = MIN( WNORM / ANORM, DBLE( N ) ) / ( N*ULP ) 00310 END IF 00311 END IF 00312 * 00313 * Do Test 2 00314 * 00315 * Compute UU' - I 00316 * 00317 IF( ITYPE.EQ.1 ) THEN 00318 CALL DGEMM( 'N', 'C', N, N, N, ONE, U, LDU, U, LDU, ZERO, WORK, 00319 $ N ) 00320 * 00321 DO 110 J = 1, N 00322 WORK( ( N+1 )*( J-1 )+1 ) = WORK( ( N+1 )*( J-1 )+1 ) - ONE 00323 110 CONTINUE 00324 * 00325 RESULT( 2 ) = MIN( DLANGE( '1', N, N, WORK, N, 00326 $ WORK( N**2+1 ) ), DBLE( N ) ) / ( N*ULP ) 00327 END IF 00328 * 00329 RETURN 00330 * 00331 * End of DSYT21 00332 * 00333 END
1.7.3 
