 LAPACK  3.10.0 LAPACK: Linear Algebra PACKage

◆ chet21()

 subroutine chet21 ( integer ITYPE, character UPLO, integer N, integer KBAND, complex, dimension( lda, * ) A, integer LDA, real, dimension( * ) D, real, dimension( * ) E, complex, dimension( ldu, * ) U, integer LDU, complex, dimension( ldv, * ) V, integer LDV, complex, dimension( * ) TAU, complex, dimension( * ) WORK, real, dimension( * ) RWORK, real, dimension( 2 ) RESULT )

CHET21

Purpose:
CHET21 generally checks a decomposition of the form

A = U S U**H

where **H means conjugate transpose, A is hermitian, U is unitary, and
S is diagonal (if KBAND=0) or (real) symmetric tridiagonal (if
KBAND=1).

If ITYPE=1, then U is represented as a dense matrix; otherwise U is
expressed as a product of Householder transformations, whose vectors
are stored in the array "V" and whose scaling constants are in "TAU".
We shall use the letter "V" to refer to the product of Householder
transformations (which should be equal to U).

Specifically, if ITYPE=1, then:

RESULT(1) = | A - U S U**H | / ( |A| n ulp ) and
RESULT(2) = | I - U U**H | / ( n ulp )

If ITYPE=2, then:

RESULT(1) = | A - V S V**H | / ( |A| n ulp )

If ITYPE=3, then:

RESULT(1) = | I - U V**H | / ( n ulp )

For ITYPE > 1, the transformation U is expressed as a product
V = H(1)...H(n-2),  where H(j) = I  -  tau(j) v(j) v(j)**H and each
vector v(j) has its first j elements 0 and the remaining n-j elements
stored in V(j+1:n,j).
Parameters
 [in] ITYPE ITYPE is INTEGER Specifies the type of tests to be performed. 1: U expressed as a dense unitary matrix: RESULT(1) = | A - U S U**H | / ( |A| n ulp ) and RESULT(2) = | I - U U**H | / ( n ulp ) 2: U expressed as a product V of Housholder transformations: RESULT(1) = | A - V S V**H | / ( |A| n ulp ) 3: U expressed both as a dense unitary matrix and as a product of Housholder transformations: RESULT(1) = | I - U V**H | / ( n ulp ) [in] UPLO UPLO is CHARACTER If UPLO='U', the upper triangle of A and V will be used and the (strictly) lower triangle will not be referenced. If UPLO='L', the lower triangle of A and V will be used and the (strictly) upper triangle will not be referenced. [in] N N is INTEGER The size of the matrix. If it is zero, CHET21 does nothing. It must be at least zero. [in] KBAND KBAND is INTEGER The bandwidth of the matrix. It may only be zero or one. If zero, then S is diagonal, and E is not referenced. If one, then S is symmetric tri-diagonal. [in] A A is COMPLEX array, dimension (LDA, N) The original (unfactored) matrix. It is assumed to be hermitian, and only the upper (UPLO='U') or only the lower (UPLO='L') will be referenced. [in] LDA LDA is INTEGER The leading dimension of A. It must be at least 1 and at least N. [in] D D is REAL array, dimension (N) The diagonal of the (symmetric tri-) diagonal matrix. [in] E E is REAL array, dimension (N-1) The off-diagonal of the (symmetric tri-) diagonal matrix. E(1) is the (1,2) and (2,1) element, E(2) is the (2,3) and (3,2) element, etc. Not referenced if KBAND=0. [in] U U is COMPLEX array, dimension (LDU, N) If ITYPE=1 or 3, this contains the unitary matrix in the decomposition, expressed as a dense matrix. If ITYPE=2, then it is not referenced. [in] LDU LDU is INTEGER The leading dimension of U. LDU must be at least N and at least 1. [in] V V is COMPLEX array, dimension (LDV, N) If ITYPE=2 or 3, the columns of this array contain the Householder vectors used to describe the unitary matrix in the decomposition. If UPLO='L', then the vectors are in the lower triangle, if UPLO='U', then in the upper triangle. *NOTE* If ITYPE=2 or 3, V is modified and restored. The subdiagonal (if UPLO='L') or the superdiagonal (if UPLO='U') is set to one, and later reset to its original value, during the course of the calculation. If ITYPE=1, then it is neither referenced nor modified. [in] LDV LDV is INTEGER The leading dimension of V. LDV must be at least N and at least 1. [in] TAU TAU is COMPLEX array, dimension (N) If ITYPE >= 2, then TAU(j) is the scalar factor of v(j) v(j)**H in the Householder transformation H(j) of the product U = H(1)...H(n-2) If ITYPE < 2, then TAU is not referenced. [out] WORK WORK is COMPLEX array, dimension (2*N**2) [out] RWORK RWORK is REAL array, dimension (N) [out] RESULT RESULT is REAL array, dimension (2) The values computed by the two tests described above. The values are currently limited to 1/ulp, to avoid overflow. RESULT(1) is always modified. RESULT(2) is modified only if ITYPE=1.

Definition at line 212 of file chet21.f.

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 = 2, 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 *
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:53
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 cgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
CGEMM
Definition: cgemm.f:187
real function clange(NORM, M, N, A, LDA, WORK)
CLANGE returns the value of the 1-norm, Frobenius norm, infinity-norm, or the largest absolute value ...
Definition: clange.f:115
real function clanhe(NORM, UPLO, N, A, LDA, WORK)
CLANHE returns the value of the 1-norm, or the Frobenius norm, or the infinity norm,...
Definition: clanhe.f:124
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 clarfy(UPLO, N, V, INCV, TAU, C, LDC, WORK)
CLARFY
Definition: clarfy.f:108
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 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
real function slamch(CMACH)
SLAMCH
Definition: slamch.f:68
Here is the call graph for this function:
Here is the caller graph for this function: