LAPACK  3.5.0
LAPACK: Linear Algebra PACKage
 All Classes Files Functions Variables Typedefs Macros
chet21.f
Go to the documentation of this file.
1 *> \brief \b CHET21
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 CHET21( 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 * REAL D( * ), E( * ), RESULT( 2 ), RWORK( * )
20 * COMPLEX A( LDA, * ), TAU( * ), U( LDU, * ),
21 * $ V( LDV, * ), WORK( * )
22 * ..
23 *
24 *
25 *> \par Purpose:
26 * =============
27 *>
28 *> \verbatim
29 *>
30 *> CHET21 generally checks a decomposition of the form
31 *>
32 *> A = U S UC>
33 *> where * means conjugate transpose, A is hermitian, U is unitary, and
34 *> S is diagonal (if KBAND=0) or (real) symmetric tridiagonal (if
35 *> KBAND=1).
36 *>
37 *> If ITYPE=1, then U is represented as a dense matrix; otherwise U is
38 *> expressed as a product of Householder transformations, whose vectors
39 *> are stored in the array "V" and whose scaling constants are in "TAU".
40 *> We shall use the letter "V" to refer to the product of Householder
41 *> transformations (which should be equal to U).
42 *>
43 *> Specifically, if ITYPE=1, then:
44 *>
45 *> RESULT(1) = | A - U S U* | / ( |A| n ulp ) *andC> RESULT(2) = | I - UU* | / ( n ulp )
46 *>
47 *> If ITYPE=2, then:
48 *>
49 *> RESULT(1) = | A - V S V* | / ( |A| n ulp )
50 *>
51 *> If ITYPE=3, then:
52 *>
53 *> RESULT(1) = | I - UV* | / ( n ulp )
54 *>
55 *> For ITYPE > 1, the transformation U is expressed as a product
56 *> V = H(1)...H(n-2), where H(j) = I - tau(j) v(j) v(j)C> and each
57 *> vector v(j) has its first j elements 0 and the remaining n-j elements
58 *> stored in V(j+1:n,j).
59 *> \endverbatim
60 *
61 * Arguments:
62 * ==========
63 *
64 *> \param[in] ITYPE
65 *> \verbatim
66 *> ITYPE is INTEGER
67 *> Specifies the type of tests to be performed.
68 *> 1: U expressed as a dense unitary matrix:
69 *> RESULT(1) = | A - U S U* | / ( |A| n ulp ) *andC> RESULT(2) = | I - UU* | / ( n ulp )
70 *>
71 *> 2: U expressed as a product V of Housholder transformations:
72 *> RESULT(1) = | A - V S V* | / ( |A| n ulp )
73 *>
74 *> 3: U expressed both as a dense unitary matrix and
75 *> as a product of Housholder transformations:
76 *> RESULT(1) = | I - UV* | / ( n ulp )
77 *> \endverbatim
78 *>
79 *> \param[in] UPLO
80 *> \verbatim
81 *> UPLO is CHARACTER
82 *> If UPLO='U', the upper triangle of A and V will be used and
83 *> the (strictly) lower triangle will not be referenced.
84 *> If UPLO='L', the lower triangle of A and V will be used and
85 *> the (strictly) upper triangle will not be referenced.
86 *> \endverbatim
87 *>
88 *> \param[in] N
89 *> \verbatim
90 *> N is INTEGER
91 *> The size of the matrix. If it is zero, CHET21 does nothing.
92 *> It must be at least zero.
93 *> \endverbatim
94 *>
95 *> \param[in] KBAND
96 *> \verbatim
97 *> KBAND is INTEGER
98 *> The bandwidth of the matrix. It may only be zero or one.
99 *> If zero, then S is diagonal, and E is not referenced. If
100 *> one, then S is symmetric tri-diagonal.
101 *> \endverbatim
102 *>
103 *> \param[in] A
104 *> \verbatim
105 *> A is COMPLEX array, dimension (LDA, N)
106 *> The original (unfactored) matrix. It is assumed to be
107 *> hermitian, and only the upper (UPLO='U') or only the lower
108 *> (UPLO='L') will be referenced.
109 *> \endverbatim
110 *>
111 *> \param[in] LDA
112 *> \verbatim
113 *> LDA is INTEGER
114 *> The leading dimension of A. It must be at least 1
115 *> and at least N.
116 *> \endverbatim
117 *>
118 *> \param[in] D
119 *> \verbatim
120 *> D is REAL array, dimension (N)
121 *> The diagonal of the (symmetric tri-) diagonal matrix.
122 *> \endverbatim
123 *>
124 *> \param[in] E
125 *> \verbatim
126 *> E is REAL array, dimension (N-1)
127 *> The off-diagonal of the (symmetric tri-) diagonal matrix.
128 *> E(1) is the (1,2) and (2,1) element, E(2) is the (2,3) and
129 *> (3,2) element, etc.
130 *> Not referenced if KBAND=0.
131 *> \endverbatim
132 *>
133 *> \param[in] U
134 *> \verbatim
135 *> U is COMPLEX array, dimension (LDU, N)
136 *> If ITYPE=1 or 3, this contains the unitary matrix in
137 *> the decomposition, expressed as a dense matrix. If ITYPE=2,
138 *> then it is not referenced.
139 *> \endverbatim
140 *>
141 *> \param[in] LDU
142 *> \verbatim
143 *> LDU is INTEGER
144 *> The leading dimension of U. LDU must be at least N and
145 *> at least 1.
146 *> \endverbatim
147 *>
148 *> \param[in] V
149 *> \verbatim
150 *> V is COMPLEX array, dimension (LDV, N)
151 *> If ITYPE=2 or 3, the columns of this array contain the
152 *> Householder vectors used to describe the unitary matrix
153 *> in the decomposition. If UPLO='L', then the vectors are in
154 *> the lower triangle, if UPLO='U', then in the upper
155 *> triangle.
156 *> *NOTE* If ITYPE=2 or 3, V is modified and restored. The
157 *> subdiagonal (if UPLO='L') or the superdiagonal (if UPLO='U')
158 *> is set to one, and later reset to its original value, during
159 *> the course of the calculation.
160 *> If ITYPE=1, then it is neither referenced nor modified.
161 *> \endverbatim
162 *>
163 *> \param[in] LDV
164 *> \verbatim
165 *> LDV is INTEGER
166 *> The leading dimension of V. LDV must be at least N and
167 *> at least 1.
168 *> \endverbatim
169 *>
170 *> \param[in] TAU
171 *> \verbatim
172 *> TAU is COMPLEX array, dimension (N)
173 *> If ITYPE >= 2, then TAU(j) is the scalar factor of
174 *> v(j) v(j)* in the Householder transformation H(j) of
175 *> the product U = H(1)...H(n-2)
176 *> If ITYPE < 2, then TAU is not referenced.
177 *> \endverbatim
178 *>
179 *> \param[out] WORK
180 *> \verbatim
181 *> WORK is COMPLEX array, dimension (2*N**2)
182 *> \endverbatim
183 *>
184 *> \param[out] RWORK
185 *> \verbatim
186 *> RWORK is REAL array, dimension (N)
187 *> \endverbatim
188 *>
189 *> \param[out] RESULT
190 *> \verbatim
191 *> RESULT is REAL array, dimension (2)
192 *> The values computed by the two tests described above. The
193 *> values are currently limited to 1/ulp, to avoid overflow.
194 *> RESULT(1) is always modified. RESULT(2) is modified only
195 *> if ITYPE=1.
196 *> \endverbatim
197 *
198 * Authors:
199 * ========
200 *
201 *> \author Univ. of Tennessee
202 *> \author Univ. of California Berkeley
203 *> \author Univ. of Colorado Denver
204 *> \author NAG Ltd.
205 *
206 *> \date November 2011
207 *
208 *> \ingroup complex_eig
209 *
210 * =====================================================================
211  SUBROUTINE chet21( ITYPE, UPLO, N, KBAND, A, LDA, D, E, U, LDU, V,
212  $ ldv, tau, work, rwork, result )
213 *
214 * -- LAPACK test routine (version 3.4.0) --
215 * -- LAPACK is a software package provided by Univ. of Tennessee, --
216 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
217 * November 2011
218 *
219 * .. Scalar Arguments ..
220  CHARACTER uplo
221  INTEGER itype, kband, lda, ldu, ldv, n
222 * ..
223 * .. Array Arguments ..
224  REAL d( * ), e( * ), result( 2 ), rwork( * )
225  COMPLEX a( lda, * ), tau( * ), u( ldu, * ),
226  $ v( ldv, * ), work( * )
227 * ..
228 *
229 * =====================================================================
230 *
231 * .. Parameters ..
232  REAL zero, one, ten
233  parameter( zero = 0.0e+0, one = 1.0e+0, ten = 10.0e+0 )
234  COMPLEX czero, cone
235  parameter( czero = ( 0.0e+0, 0.0e+0 ),
236  $ cone = ( 1.0e+0, 0.0e+0 ) )
237 * ..
238 * .. Local Scalars ..
239  LOGICAL lower
240  CHARACTER cuplo
241  INTEGER iinfo, j, jcol, jr, jrow
242  REAL anorm, ulp, unfl, wnorm
243  COMPLEX vsave
244 * ..
245 * .. External Functions ..
246  LOGICAL lsame
247  REAL clange, clanhe, slamch
248  EXTERNAL lsame, clange, clanhe, slamch
249 * ..
250 * .. External Subroutines ..
251  EXTERNAL cgemm, cher, cher2, clacpy, clarfy, claset,
252  $ cunm2l, cunm2r
253 * ..
254 * .. Intrinsic Functions ..
255  INTRINSIC cmplx, max, min, real
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 = slamch( 'Safe minimum' )
274  ulp = slamch( 'Epsilon' )*slamch( '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( clanhe( '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*
298 *
299  CALL claset( 'Full', n, n, czero, czero, work, n )
300  CALL clacpy( cuplo, n, n, a, lda, work, n )
301 *
302  DO 10 j = 1, n
303  CALL cher( 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 = 1, n - 1
308  CALL cher2( cuplo, n, -cmplx( e( j ) ), u( 1, j ), 1,
309  $ u( 1, j-1 ), 1, work, n )
310  20 CONTINUE
311  END IF
312  wnorm = clanhe( '1', cuplo, n, work, n, rwork )
313 *
314  ELSE IF( itype.EQ.2 ) THEN
315 *
316 * ITYPE=2: error = V S V* - A
317 *
318  CALL claset( '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 clarfy( '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 clarfy( '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 = clanhe( '1', cuplo, n, work, n, rwork )
370 *
371  ELSE IF( itype.EQ.3 ) THEN
372 *
373 * ITYPE=3: error = U V* - I
374 *
375  IF( n.LT.2 )
376  $ RETURN
377  CALL clacpy( ' ', n, n, u, ldu, work, n )
378  IF( lower ) THEN
379  CALL cunm2r( '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 cunm2l( '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 = clange( '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, REAL( N ) ) / ( n*ulp )
404  END IF
405  END IF
406 *
407 * Do Test 2
408 *
409 * Compute UU* - I
410 *
411  IF( itype.EQ.1 ) THEN
412  CALL cgemm( '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( clange( '1', n, n, work, n, rwork ),
420  $ REAL( N ) ) / ( n*ulp )
421  END IF
422 *
423  RETURN
424 *
425 * End of CHET21
426 *
427  END