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