LAPACK 3.12.0
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches
zhet22.f
Go to the documentation of this file.
1*> \brief \b ZHET22
2*
3* =========== DOCUMENTATION ===========
4*
5* Online html documentation available at
6* http://www.netlib.org/lapack/explore-html/
7*
8* Definition:
9* ===========
10*
11* SUBROUTINE ZHET22( ITYPE, UPLO, N, M, KBAND, A, LDA, D, E, U, LDU,
12* V, LDV, TAU, WORK, RWORK, RESULT )
13*
14* .. Scalar Arguments ..
15* CHARACTER UPLO
16* INTEGER ITYPE, KBAND, LDA, LDU, LDV, M, N
17* ..
18* .. Array Arguments ..
19* DOUBLE PRECISION D( * ), E( * ), RESULT( 2 ), RWORK( * )
20* COMPLEX*16 A( LDA, * ), TAU( * ), U( LDU, * ),
21* $ V( LDV, * ), WORK( * )
22* ..
23*
24*
25*> \par Purpose:
26* =============
27*>
28*> \verbatim
29*>
30*> ZHET22 generally checks a decomposition of the form
31*>
32*> A U = U S
33*>
34*> where A is complex Hermitian, the columns of U are orthonormal,
35*> and S is diagonal (if KBAND=0) or symmetric tridiagonal (if
36*> KBAND=1). If ITYPE=1, then U is represented as a dense matrix,
37*> otherwise the U is expressed as a product of Householder
38*> transformations, whose vectors are stored in the array "V" and
39*> whose scaling constants are in "TAU"; we shall use the letter
40*> "V" to refer to the product of Householder transformations
41*> (which should be equal to U).
42*>
43*> Specifically, if ITYPE=1, then:
44*>
45*> RESULT(1) = | U**H A U - S | / ( |A| m ulp ) and
46*> RESULT(2) = | I - U**H U | / ( m ulp )
47*> \endverbatim
48*
49* Arguments:
50* ==========
51*
52*> \verbatim
53*> ITYPE INTEGER
54*> Specifies the type of tests to be performed.
55*> 1: U expressed as a dense orthogonal matrix:
56*> RESULT(1) = | A - U S U**H | / ( |A| n ulp ) *and
57*> RESULT(2) = | I - U U**H | / ( n ulp )
58*>
59*> UPLO CHARACTER
60*> If UPLO='U', the upper triangle of A will be used and the
61*> (strictly) lower triangle will not be referenced. If
62*> UPLO='L', the lower triangle of A will be used and the
63*> (strictly) upper triangle will not be referenced.
64*> Not modified.
65*>
66*> N INTEGER
67*> The size of the matrix. If it is zero, ZHET22 does nothing.
68*> It must be at least zero.
69*> Not modified.
70*>
71*> M INTEGER
72*> The number of columns of U. If it is zero, ZHET22 does
73*> nothing. It must be at least zero.
74*> Not modified.
75*>
76*> KBAND INTEGER
77*> The bandwidth of the matrix. It may only be zero or one.
78*> If zero, then S is diagonal, and E is not referenced. If
79*> one, then S is symmetric tri-diagonal.
80*> Not modified.
81*>
82*> A COMPLEX*16 array, dimension (LDA , N)
83*> The original (unfactored) matrix. It is assumed to be
84*> symmetric, and only the upper (UPLO='U') or only the lower
85*> (UPLO='L') will be referenced.
86*> Not modified.
87*>
88*> LDA INTEGER
89*> The leading dimension of A. It must be at least 1
90*> and at least N.
91*> Not modified.
92*>
93*> D DOUBLE PRECISION array, dimension (N)
94*> The diagonal of the (symmetric tri-) diagonal matrix.
95*> Not modified.
96*>
97*> E DOUBLE PRECISION array, dimension (N)
98*> The off-diagonal of the (symmetric tri-) diagonal matrix.
99*> E(1) is ignored, E(2) is the (1,2) and (2,1) element, etc.
100*> Not referenced if KBAND=0.
101*> Not modified.
102*>
103*> U COMPLEX*16 array, dimension (LDU, N)
104*> If ITYPE=1, this contains the orthogonal matrix in
105*> the decomposition, expressed as a dense matrix.
106*> Not modified.
107*>
108*> LDU INTEGER
109*> The leading dimension of U. LDU must be at least N and
110*> at least 1.
111*> Not modified.
112*>
113*> V COMPLEX*16 array, dimension (LDV, N)
114*> If ITYPE=2 or 3, the lower triangle of this array contains
115*> the Householder vectors used to describe the orthogonal
116*> matrix in the decomposition. If ITYPE=1, then it is not
117*> referenced.
118*> Not modified.
119*>
120*> LDV INTEGER
121*> The leading dimension of V. LDV must be at least N and
122*> at least 1.
123*> Not modified.
124*>
125*> TAU COMPLEX*16 array, dimension (N)
126*> If ITYPE >= 2, then TAU(j) is the scalar factor of
127*> v(j) v(j)**H in the Householder transformation H(j) of
128*> the product U = H(1)...H(n-2)
129*> If ITYPE < 2, then TAU is not referenced.
130*> Not modified.
131*>
132*> WORK COMPLEX*16 array, dimension (2*N**2)
133*> Workspace.
134*> Modified.
135*>
136*> RWORK DOUBLE PRECISION array, dimension (N)
137*> Workspace.
138*> Modified.
139*>
140*> RESULT DOUBLE PRECISION array, dimension (2)
141*> The values computed by the two tests described above. The
142*> values are currently limited to 1/ulp, to avoid overflow.
143*> RESULT(1) is always modified. RESULT(2) is modified only
144*> if LDU is at least N.
145*> Modified.
146*> \endverbatim
147*
148* Authors:
149* ========
150*
151*> \author Univ. of Tennessee
152*> \author Univ. of California Berkeley
153*> \author Univ. of Colorado Denver
154*> \author NAG Ltd.
155*
156*> \ingroup complex16_eig
157*
158* =====================================================================
159 SUBROUTINE zhet22( ITYPE, UPLO, N, M, KBAND, A, LDA, D, E, U, LDU,
160 $ V, LDV, TAU, WORK, RWORK, RESULT )
161*
162* -- LAPACK test routine --
163* -- LAPACK is a software package provided by Univ. of Tennessee, --
164* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
165*
166* .. Scalar Arguments ..
167 CHARACTER UPLO
168 INTEGER ITYPE, KBAND, LDA, LDU, LDV, M, N
169* ..
170* .. Array Arguments ..
171 DOUBLE PRECISION D( * ), E( * ), RESULT( 2 ), RWORK( * )
172 COMPLEX*16 A( LDA, * ), TAU( * ), U( LDU, * ),
173 $ v( ldv, * ), work( * )
174* ..
175*
176* =====================================================================
177*
178* .. Parameters ..
179 DOUBLE PRECISION ZERO, ONE
180 parameter( zero = 0.0d0, one = 1.0d0 )
181 COMPLEX*16 CZERO, CONE
182 parameter( czero = ( 0.0d0, 0.0d0 ),
183 $ cone = ( 1.0d0, 0.0d0 ) )
184* ..
185* .. Local Scalars ..
186 INTEGER J, JJ, JJ1, JJ2, NN, NNP1
187 DOUBLE PRECISION ANORM, ULP, UNFL, WNORM
188* ..
189* .. External Functions ..
190 DOUBLE PRECISION DLAMCH, ZLANHE
191 EXTERNAL dlamch, zlanhe
192* ..
193* .. External Subroutines ..
194 EXTERNAL zgemm, zhemm, zunt01
195* ..
196* .. Intrinsic Functions ..
197 INTRINSIC dble, max, min
198* ..
199* .. Executable Statements ..
200*
201 result( 1 ) = zero
202 result( 2 ) = zero
203 IF( n.LE.0 .OR. m.LE.0 )
204 $ RETURN
205*
206 unfl = dlamch( 'Safe minimum' )
207 ulp = dlamch( 'Precision' )
208*
209* Do Test 1
210*
211* Norm of A:
212*
213 anorm = max( zlanhe( '1', uplo, n, a, lda, rwork ), unfl )
214*
215* Compute error matrix:
216*
217* ITYPE=1: error = U**H A U - S
218*
219 CALL zhemm( 'L', uplo, n, m, cone, a, lda, u, ldu, czero, work,
220 $ n )
221 nn = n*n
222 nnp1 = nn + 1
223 CALL zgemm( 'C', 'N', m, m, n, cone, u, ldu, work, n, czero,
224 $ work( nnp1 ), n )
225 DO 10 j = 1, m
226 jj = nn + ( j-1 )*n + j
227 work( jj ) = work( jj ) - d( j )
228 10 CONTINUE
229 IF( kband.EQ.1 .AND. n.GT.1 ) THEN
230 DO 20 j = 2, m
231 jj1 = nn + ( j-1 )*n + j - 1
232 jj2 = nn + ( j-2 )*n + j
233 work( jj1 ) = work( jj1 ) - e( j-1 )
234 work( jj2 ) = work( jj2 ) - e( j-1 )
235 20 CONTINUE
236 END IF
237 wnorm = zlanhe( '1', uplo, m, work( nnp1 ), n, rwork )
238*
239 IF( anorm.GT.wnorm ) THEN
240 result( 1 ) = ( wnorm / anorm ) / ( m*ulp )
241 ELSE
242 IF( anorm.LT.one ) THEN
243 result( 1 ) = ( min( wnorm, m*anorm ) / anorm ) / ( m*ulp )
244 ELSE
245 result( 1 ) = min( wnorm / anorm, dble( m ) ) / ( m*ulp )
246 END IF
247 END IF
248*
249* Do Test 2
250*
251* Compute U**H U - I
252*
253 IF( itype.EQ.1 )
254 $ CALL zunt01( 'Columns', n, m, u, ldu, work, 2*n*n, rwork,
255 $ result( 2 ) )
256*
257 RETURN
258*
259* End of ZHET22
260*
261 END
subroutine zgemm(transa, transb, m, n, k, alpha, a, lda, b, ldb, beta, c, ldc)
ZGEMM
Definition zgemm.f:188
subroutine zhemm(side, uplo, m, n, alpha, a, lda, b, ldb, beta, c, ldc)
ZHEMM
Definition zhemm.f:191
subroutine zhet22(itype, uplo, n, m, kband, a, lda, d, e, u, ldu, v, ldv, tau, work, rwork, result)
ZHET22
Definition zhet22.f:161
subroutine zunt01(rowcol, m, n, u, ldu, work, lwork, rwork, resid)
ZUNT01
Definition zunt01.f:126