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