LAPACK 3.12.0
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches
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 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, CHET21 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 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 REAL array, dimension (N)
124*> The diagonal of the (symmetric tri-) diagonal matrix.
125*> \endverbatim
126*>
127*> \param[in] E
128*> \verbatim
129*> E is REAL 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 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 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 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 array, dimension (2*N**2)
185*> \endverbatim
186*>
187*> \param[out] RWORK
188*> \verbatim
189*> RWORK is REAL array, dimension (N)
190*> \endverbatim
191*>
192*> \param[out] RESULT
193*> \verbatim
194*> RESULT is REAL 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 complex_eig
210*
211* =====================================================================
212 SUBROUTINE chet21( 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 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**H
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**H - 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**H - 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 U U**H - 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
subroutine chet21(itype, uplo, n, kband, a, lda, d, e, u, ldu, v, ldv, tau, work, rwork, result)
CHET21
Definition chet21.f:214
subroutine cgemm(transa, transb, m, n, k, alpha, a, lda, b, ldb, beta, c, ldc)
CGEMM
Definition cgemm.f:188
subroutine cher2(uplo, n, alpha, x, incx, y, incy, a, lda)
CHER2
Definition cher2.f:150
subroutine cher(uplo, n, alpha, x, incx, a, lda)
CHER
Definition cher.f:135
subroutine clacpy(uplo, m, n, a, lda, b, ldb)
CLACPY copies all or part of one two-dimensional array to another.
Definition clacpy.f:103
subroutine clarfy(uplo, n, v, incv, tau, c, ldc, work)
CLARFY
Definition clarfy.f:108
subroutine claset(uplo, m, n, alpha, beta, a, lda)
CLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values.
Definition claset.f:106
subroutine cunm2l(side, trans, m, n, k, a, lda, tau, c, ldc, work, info)
CUNM2L multiplies a general matrix by the unitary matrix from a QL factorization determined by cgeqlf...
Definition cunm2l.f:159
subroutine cunm2r(side, trans, m, n, k, a, lda, tau, c, ldc, work, info)
CUNM2R multiplies a general matrix by the unitary matrix from a QR factorization determined by cgeqrf...
Definition cunm2r.f:159