LAPACK  3.10.0
LAPACK: Linear Algebra PACKage
zhet21.f
Go to the documentation of this file.
1 *> \brief \b ZHET21
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 ZHET21( ITYPE, UPLO, N, KBAND, A, LDA, D, E, U, LDU, V,
12 * LDV, TAU, WORK, RWORK, RESULT )
13 *
14 * .. Scalar Arguments ..
15 * CHARACTER UPLO
16 * INTEGER ITYPE, KBAND, LDA, LDU, LDV, 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 *> ZHET21 generally checks a decomposition of the form
31 *>
32 *> A = U S U**H
33 *>
34 *> where **H means conjugate transpose, A is hermitian, U is unitary, and
35 *> S is diagonal (if KBAND=0) or (real) symmetric tridiagonal (if
36 *> KBAND=1).
37 *>
38 *> If ITYPE=1, then U is represented as a dense matrix; otherwise U is
39 *> expressed as a product of Householder transformations, whose vectors
40 *> are stored in the array "V" and whose scaling constants are in "TAU".
41 *> We shall use the letter "V" to refer to the product of Householder
42 *> transformations (which should be equal to U).
43 *>
44 *> Specifically, if ITYPE=1, then:
45 *>
46 *> RESULT(1) = | A - U S U**H | / ( |A| n ulp ) and
47 *> RESULT(2) = | I - U U**H | / ( n ulp )
48 *>
49 *> If ITYPE=2, then:
50 *>
51 *> RESULT(1) = | A - V S V**H | / ( |A| n ulp )
52 *>
53 *> If ITYPE=3, then:
54 *>
55 *> RESULT(1) = | I - U V**H | / ( n ulp )
56 *>
57 *> For ITYPE > 1, the transformation U is expressed as a product
58 *> V = H(1)...H(n-2), where H(j) = I - tau(j) v(j) v(j)**H and each
59 *> vector v(j) has its first j elements 0 and the remaining n-j elements
60 *> stored in V(j+1:n,j).
61 *> \endverbatim
62 *
63 * Arguments:
64 * ==========
65 *
66 *> \param[in] ITYPE
67 *> \verbatim
68 *> ITYPE is INTEGER
69 *> Specifies the type of tests to be performed.
70 *> 1: U expressed as a dense unitary matrix:
71 *> RESULT(1) = | A - U S U**H | / ( |A| n ulp ) and
72 *> RESULT(2) = | I - U U**H | / ( n ulp )
73 *>
74 *> 2: U expressed as a product V of Housholder transformations:
75 *> RESULT(1) = | A - V S V**H | / ( |A| n ulp )
76 *>
77 *> 3: U expressed both as a dense unitary matrix and
78 *> as a product of Housholder transformations:
79 *> RESULT(1) = | I - U V**H | / ( n ulp )
80 *> \endverbatim
81 *>
82 *> \param[in] UPLO
83 *> \verbatim
84 *> UPLO is CHARACTER
85 *> If UPLO='U', the upper triangle of A and V will be used and
86 *> the (strictly) lower triangle will not be referenced.
87 *> If UPLO='L', the lower triangle of A and V will be used and
88 *> the (strictly) upper triangle will not be referenced.
89 *> \endverbatim
90 *>
91 *> \param[in] N
92 *> \verbatim
93 *> N is INTEGER
94 *> The size of the matrix. If it is zero, ZHET21 does nothing.
95 *> It must be at least zero.
96 *> \endverbatim
97 *>
98 *> \param[in] KBAND
99 *> \verbatim
100 *> KBAND is INTEGER
101 *> The bandwidth of the matrix. It may only be zero or one.
102 *> If zero, then S is diagonal, and E is not referenced. If
103 *> one, then S is symmetric tri-diagonal.
104 *> \endverbatim
105 *>
106 *> \param[in] A
107 *> \verbatim
108 *> A is COMPLEX*16 array, dimension (LDA, N)
109 *> The original (unfactored) matrix. It is assumed to be
110 *> hermitian, and only the upper (UPLO='U') or only the lower
111 *> (UPLO='L') will be referenced.
112 *> \endverbatim
113 *>
114 *> \param[in] LDA
115 *> \verbatim
116 *> LDA is INTEGER
117 *> The leading dimension of A. It must be at least 1
118 *> and at least N.
119 *> \endverbatim
120 *>
121 *> \param[in] D
122 *> \verbatim
123 *> D is DOUBLE PRECISION array, dimension (N)
124 *> The diagonal of the (symmetric tri-) diagonal matrix.
125 *> \endverbatim
126 *>
127 *> \param[in] E
128 *> \verbatim
129 *> E is DOUBLE PRECISION array, dimension (N-1)
130 *> The off-diagonal of the (symmetric tri-) diagonal matrix.
131 *> E(1) is the (1,2) and (2,1) element, E(2) is the (2,3) and
132 *> (3,2) element, etc.
133 *> Not referenced if KBAND=0.
134 *> \endverbatim
135 *>
136 *> \param[in] U
137 *> \verbatim
138 *> U is COMPLEX*16 array, dimension (LDU, N)
139 *> If ITYPE=1 or 3, this contains the unitary matrix in
140 *> the decomposition, expressed as a dense matrix. If ITYPE=2,
141 *> then it is not referenced.
142 *> \endverbatim
143 *>
144 *> \param[in] LDU
145 *> \verbatim
146 *> LDU is INTEGER
147 *> The leading dimension of U. LDU must be at least N and
148 *> at least 1.
149 *> \endverbatim
150 *>
151 *> \param[in] V
152 *> \verbatim
153 *> V is COMPLEX*16 array, dimension (LDV, N)
154 *> If ITYPE=2 or 3, the columns of this array contain the
155 *> Householder vectors used to describe the unitary matrix
156 *> in the decomposition. If UPLO='L', then the vectors are in
157 *> the lower triangle, if UPLO='U', then in the upper
158 *> triangle.
159 *> *NOTE* If ITYPE=2 or 3, V is modified and restored. The
160 *> subdiagonal (if UPLO='L') or the superdiagonal (if UPLO='U')
161 *> is set to one, and later reset to its original value, during
162 *> the course of the calculation.
163 *> If ITYPE=1, then it is neither referenced nor modified.
164 *> \endverbatim
165 *>
166 *> \param[in] LDV
167 *> \verbatim
168 *> LDV is INTEGER
169 *> The leading dimension of V. LDV must be at least N and
170 *> at least 1.
171 *> \endverbatim
172 *>
173 *> \param[in] TAU
174 *> \verbatim
175 *> TAU is COMPLEX*16 array, dimension (N)
176 *> If ITYPE >= 2, then TAU(j) is the scalar factor of
177 *> v(j) v(j)**H in the Householder transformation H(j) of
178 *> the product U = H(1)...H(n-2)
179 *> If ITYPE < 2, then TAU is not referenced.
180 *> \endverbatim
181 *>
182 *> \param[out] WORK
183 *> \verbatim
184 *> WORK is COMPLEX*16 array, dimension (2*N**2)
185 *> \endverbatim
186 *>
187 *> \param[out] RWORK
188 *> \verbatim
189 *> RWORK is DOUBLE PRECISION array, dimension (N)
190 *> \endverbatim
191 *>
192 *> \param[out] RESULT
193 *> \verbatim
194 *> RESULT is DOUBLE PRECISION array, dimension (2)
195 *> The values computed by the two tests described above. The
196 *> values are currently limited to 1/ulp, to avoid overflow.
197 *> RESULT(1) is always modified. RESULT(2) is modified only
198 *> if ITYPE=1.
199 *> \endverbatim
200 *
201 * Authors:
202 * ========
203 *
204 *> \author Univ. of Tennessee
205 *> \author Univ. of California Berkeley
206 *> \author Univ. of Colorado Denver
207 *> \author NAG Ltd.
208 *
209 *> \ingroup complex16_eig
210 *
211 * =====================================================================
212  SUBROUTINE zhet21( ITYPE, UPLO, N, KBAND, A, LDA, D, E, U, LDU, V,
213  $ LDV, TAU, WORK, RWORK, RESULT )
214 *
215 * -- LAPACK test routine --
216 * -- LAPACK is a software package provided by Univ. of Tennessee, --
217 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
218 *
219 * .. Scalar Arguments ..
220  CHARACTER UPLO
221  INTEGER ITYPE, KBAND, LDA, LDU, LDV, N
222 * ..
223 * .. Array Arguments ..
224  DOUBLE PRECISION D( * ), E( * ), RESULT( 2 ), RWORK( * )
225  COMPLEX*16 A( LDA, * ), TAU( * ), U( LDU, * ),
226  $ v( ldv, * ), work( * )
227 * ..
228 *
229 * =====================================================================
230 *
231 * .. Parameters ..
232  DOUBLE PRECISION ZERO, ONE, TEN
233  parameter( zero = 0.0d+0, one = 1.0d+0, ten = 10.0d+0 )
234  COMPLEX*16 CZERO, CONE
235  parameter( czero = ( 0.0d+0, 0.0d+0 ),
236  $ cone = ( 1.0d+0, 0.0d+0 ) )
237 * ..
238 * .. Local Scalars ..
239  LOGICAL LOWER
240  CHARACTER CUPLO
241  INTEGER IINFO, J, JCOL, JR, JROW
242  DOUBLE PRECISION ANORM, ULP, UNFL, WNORM
243  COMPLEX*16 VSAVE
244 * ..
245 * .. External Functions ..
246  LOGICAL LSAME
247  DOUBLE PRECISION DLAMCH, ZLANGE, ZLANHE
248  EXTERNAL lsame, dlamch, zlange, zlanhe
249 * ..
250 * .. External Subroutines ..
251  EXTERNAL zgemm, zher, zher2, zlacpy, zlarfy, zlaset,
252  $ zunm2l, zunm2r
253 * ..
254 * .. Intrinsic Functions ..
255  INTRINSIC dble, dcmplx, max, min
256 * ..
257 * .. Executable Statements ..
258 *
259  result( 1 ) = zero
260  IF( itype.EQ.1 )
261  $ result( 2 ) = zero
262  IF( n.LE.0 )
263  $ RETURN
264 *
265  IF( lsame( uplo, 'U' ) ) THEN
266  lower = .false.
267  cuplo = 'U'
268  ELSE
269  lower = .true.
270  cuplo = 'L'
271  END IF
272 *
273  unfl = dlamch( 'Safe minimum' )
274  ulp = dlamch( 'Epsilon' )*dlamch( 'Base' )
275 *
276 * Some Error Checks
277 *
278  IF( itype.LT.1 .OR. itype.GT.3 ) THEN
279  result( 1 ) = ten / ulp
280  RETURN
281  END IF
282 *
283 * Do Test 1
284 *
285 * Norm of A:
286 *
287  IF( itype.EQ.3 ) THEN
288  anorm = one
289  ELSE
290  anorm = max( zlanhe( '1', cuplo, n, a, lda, rwork ), unfl )
291  END IF
292 *
293 * Compute error matrix:
294 *
295  IF( itype.EQ.1 ) THEN
296 *
297 * ITYPE=1: error = A - U S U**H
298 *
299  CALL zlaset( 'Full', n, n, czero, czero, work, n )
300  CALL zlacpy( cuplo, n, n, a, lda, work, n )
301 *
302  DO 10 j = 1, n
303  CALL zher( cuplo, n, -d( j ), u( 1, j ), 1, work, n )
304  10 CONTINUE
305 *
306  IF( n.GT.1 .AND. kband.EQ.1 ) THEN
307  DO 20 j = 2, n - 1
308  CALL zher2( cuplo, n, -dcmplx( e( j ) ), u( 1, j ), 1,
309  $ u( 1, j-1 ), 1, work, n )
310  20 CONTINUE
311  END IF
312  wnorm = zlanhe( '1', cuplo, n, work, n, rwork )
313 *
314  ELSE IF( itype.EQ.2 ) THEN
315 *
316 * ITYPE=2: error = V S V**H - A
317 *
318  CALL zlaset( 'Full', n, n, czero, czero, work, n )
319 *
320  IF( lower ) THEN
321  work( n**2 ) = d( n )
322  DO 40 j = n - 1, 1, -1
323  IF( kband.EQ.1 ) THEN
324  work( ( n+1 )*( j-1 )+2 ) = ( cone-tau( j ) )*e( j )
325  DO 30 jr = j + 2, n
326  work( ( j-1 )*n+jr ) = -tau( j )*e( j )*v( jr, j )
327  30 CONTINUE
328  END IF
329 *
330  vsave = v( j+1, j )
331  v( j+1, j ) = one
332  CALL zlarfy( 'L', n-j, v( j+1, j ), 1, tau( j ),
333  $ work( ( n+1 )*j+1 ), n, work( n**2+1 ) )
334  v( j+1, j ) = vsave
335  work( ( n+1 )*( j-1 )+1 ) = d( j )
336  40 CONTINUE
337  ELSE
338  work( 1 ) = d( 1 )
339  DO 60 j = 1, n - 1
340  IF( kband.EQ.1 ) THEN
341  work( ( n+1 )*j ) = ( cone-tau( j ) )*e( j )
342  DO 50 jr = 1, j - 1
343  work( j*n+jr ) = -tau( j )*e( j )*v( jr, j+1 )
344  50 CONTINUE
345  END IF
346 *
347  vsave = v( j, j+1 )
348  v( j, j+1 ) = one
349  CALL zlarfy( 'U', j, v( 1, j+1 ), 1, tau( j ), work, n,
350  $ work( n**2+1 ) )
351  v( j, j+1 ) = vsave
352  work( ( n+1 )*j+1 ) = d( j+1 )
353  60 CONTINUE
354  END IF
355 *
356  DO 90 jcol = 1, n
357  IF( lower ) THEN
358  DO 70 jrow = jcol, n
359  work( jrow+n*( jcol-1 ) ) = work( jrow+n*( jcol-1 ) )
360  $ - a( jrow, jcol )
361  70 CONTINUE
362  ELSE
363  DO 80 jrow = 1, jcol
364  work( jrow+n*( jcol-1 ) ) = work( jrow+n*( jcol-1 ) )
365  $ - a( jrow, jcol )
366  80 CONTINUE
367  END IF
368  90 CONTINUE
369  wnorm = zlanhe( '1', cuplo, n, work, n, rwork )
370 *
371  ELSE IF( itype.EQ.3 ) THEN
372 *
373 * ITYPE=3: error = U V**H - I
374 *
375  IF( n.LT.2 )
376  $ RETURN
377  CALL zlacpy( ' ', n, n, u, ldu, work, n )
378  IF( lower ) THEN
379  CALL zunm2r( 'R', 'C', n, n-1, n-1, v( 2, 1 ), ldv, tau,
380  $ work( n+1 ), n, work( n**2+1 ), iinfo )
381  ELSE
382  CALL zunm2l( 'R', 'C', n, n-1, n-1, v( 1, 2 ), ldv, tau,
383  $ work, n, work( n**2+1 ), iinfo )
384  END IF
385  IF( iinfo.NE.0 ) THEN
386  result( 1 ) = ten / ulp
387  RETURN
388  END IF
389 *
390  DO 100 j = 1, n
391  work( ( n+1 )*( j-1 )+1 ) = work( ( n+1 )*( j-1 )+1 ) - cone
392  100 CONTINUE
393 *
394  wnorm = zlange( '1', n, n, work, n, rwork )
395  END IF
396 *
397  IF( anorm.GT.wnorm ) THEN
398  result( 1 ) = ( wnorm / anorm ) / ( n*ulp )
399  ELSE
400  IF( anorm.LT.one ) THEN
401  result( 1 ) = ( min( wnorm, n*anorm ) / anorm ) / ( n*ulp )
402  ELSE
403  result( 1 ) = min( wnorm / anorm, dble( n ) ) / ( n*ulp )
404  END IF
405  END IF
406 *
407 * Do Test 2
408 *
409 * Compute U U**H - I
410 *
411  IF( itype.EQ.1 ) THEN
412  CALL zgemm( 'N', 'C', n, n, n, cone, u, ldu, u, ldu, czero,
413  $ work, n )
414 *
415  DO 110 j = 1, n
416  work( ( n+1 )*( j-1 )+1 ) = work( ( n+1 )*( j-1 )+1 ) - cone
417  110 CONTINUE
418 *
419  result( 2 ) = min( zlange( '1', n, n, work, n, rwork ),
420  $ dble( n ) ) / ( n*ulp )
421  END IF
422 *
423  RETURN
424 *
425 * End of ZHET21
426 *
427  END
subroutine zher(UPLO, N, ALPHA, X, INCX, A, LDA)
ZHER
Definition: zher.f:135
subroutine zher2(UPLO, N, ALPHA, X, INCX, Y, INCY, A, LDA)
ZHER2
Definition: zher2.f:150
subroutine zgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
ZGEMM
Definition: zgemm.f:187
subroutine zhet21(ITYPE, UPLO, N, KBAND, A, LDA, D, E, U, LDU, V, LDV, TAU, WORK, RWORK, RESULT)
ZHET21
Definition: zhet21.f:214
subroutine zlacpy(UPLO, M, N, A, LDA, B, LDB)
ZLACPY copies all or part of one two-dimensional array to another.
Definition: zlacpy.f:103
subroutine zlarfy(UPLO, N, V, INCV, TAU, C, LDC, WORK)
ZLARFY
Definition: zlarfy.f:108
subroutine zlaset(UPLO, M, N, ALPHA, BETA, A, LDA)
ZLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values.
Definition: zlaset.f:106
subroutine zunm2l(SIDE, TRANS, M, N, K, A, LDA, TAU, C, LDC, WORK, INFO)
ZUNM2L multiplies a general matrix by the unitary matrix from a QL factorization determined by cgeqlf...
Definition: zunm2l.f:159
subroutine zunm2r(SIDE, TRANS, M, N, K, A, LDA, TAU, C, LDC, WORK, INFO)
ZUNM2R multiplies a general matrix by the unitary matrix from a QR factorization determined by cgeqrf...
Definition: zunm2r.f:159