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