LAPACK  3.9.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.
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Date
December 2016

Definition at line 216 of file chet21.f.

216 *
217 * -- LAPACK test routine (version 3.7.0) --
218 * -- LAPACK is a software package provided by Univ. of Tennessee, --
219 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
220 * December 2016
221 *
222 * .. Scalar Arguments ..
223  CHARACTER UPLO
224  INTEGER ITYPE, KBAND, LDA, LDU, LDV, N
225 * ..
226 * .. Array Arguments ..
227  REAL D( * ), E( * ), RESULT( 2 ), RWORK( * )
228  COMPLEX A( LDA, * ), TAU( * ), U( LDU, * ),
229  $ V( LDV, * ), WORK( * )
230 * ..
231 *
232 * =====================================================================
233 *
234 * .. Parameters ..
235  REAL ZERO, ONE, TEN
236  parameter( zero = 0.0e+0, one = 1.0e+0, ten = 10.0e+0 )
237  COMPLEX CZERO, CONE
238  parameter( czero = ( 0.0e+0, 0.0e+0 ),
239  $ cone = ( 1.0e+0, 0.0e+0 ) )
240 * ..
241 * .. Local Scalars ..
242  LOGICAL LOWER
243  CHARACTER CUPLO
244  INTEGER IINFO, J, JCOL, JR, JROW
245  REAL ANORM, ULP, UNFL, WNORM
246  COMPLEX VSAVE
247 * ..
248 * .. External Functions ..
249  LOGICAL LSAME
250  REAL CLANGE, CLANHE, SLAMCH
251  EXTERNAL lsame, clange, clanhe, slamch
252 * ..
253 * .. External Subroutines ..
254  EXTERNAL cgemm, cher, cher2, clacpy, clarfy, claset,
255  $ cunm2l, cunm2r
256 * ..
257 * .. Intrinsic Functions ..
258  INTRINSIC cmplx, max, min, real
259 * ..
260 * .. Executable Statements ..
261 *
262  result( 1 ) = zero
263  IF( itype.EQ.1 )
264  $ result( 2 ) = zero
265  IF( n.LE.0 )
266  $ RETURN
267 *
268  IF( lsame( uplo, 'U' ) ) THEN
269  lower = .false.
270  cuplo = 'U'
271  ELSE
272  lower = .true.
273  cuplo = 'L'
274  END IF
275 *
276  unfl = slamch( 'Safe minimum' )
277  ulp = slamch( 'Epsilon' )*slamch( 'Base' )
278 *
279 * Some Error Checks
280 *
281  IF( itype.LT.1 .OR. itype.GT.3 ) THEN
282  result( 1 ) = ten / ulp
283  RETURN
284  END IF
285 *
286 * Do Test 1
287 *
288 * Norm of A:
289 *
290  IF( itype.EQ.3 ) THEN
291  anorm = one
292  ELSE
293  anorm = max( clanhe( '1', cuplo, n, a, lda, rwork ), unfl )
294  END IF
295 *
296 * Compute error matrix:
297 *
298  IF( itype.EQ.1 ) THEN
299 *
300 * ITYPE=1: error = A - U S U**H
301 *
302  CALL claset( 'Full', n, n, czero, czero, work, n )
303  CALL clacpy( cuplo, n, n, a, lda, work, n )
304 *
305  DO 10 j = 1, n
306  CALL cher( cuplo, n, -d( j ), u( 1, j ), 1, work, n )
307  10 CONTINUE
308 *
309  IF( n.GT.1 .AND. kband.EQ.1 ) THEN
310  DO 20 j = 1, n - 1
311  CALL cher2( cuplo, n, -cmplx( e( j ) ), u( 1, j ), 1,
312  $ u( 1, j-1 ), 1, work, n )
313  20 CONTINUE
314  END IF
315  wnorm = clanhe( '1', cuplo, n, work, n, rwork )
316 *
317  ELSE IF( itype.EQ.2 ) THEN
318 *
319 * ITYPE=2: error = V S V**H - A
320 *
321  CALL claset( 'Full', n, n, czero, czero, work, n )
322 *
323  IF( lower ) THEN
324  work( n**2 ) = d( n )
325  DO 40 j = n - 1, 1, -1
326  IF( kband.EQ.1 ) THEN
327  work( ( n+1 )*( j-1 )+2 ) = ( cone-tau( j ) )*e( j )
328  DO 30 jr = j + 2, n
329  work( ( j-1 )*n+jr ) = -tau( j )*e( j )*v( jr, j )
330  30 CONTINUE
331  END IF
332 *
333  vsave = v( j+1, j )
334  v( j+1, j ) = one
335  CALL clarfy( 'L', n-j, v( j+1, j ), 1, tau( j ),
336  $ work( ( n+1 )*j+1 ), n, work( n**2+1 ) )
337  v( j+1, j ) = vsave
338  work( ( n+1 )*( j-1 )+1 ) = d( j )
339  40 CONTINUE
340  ELSE
341  work( 1 ) = d( 1 )
342  DO 60 j = 1, n - 1
343  IF( kband.EQ.1 ) THEN
344  work( ( n+1 )*j ) = ( cone-tau( j ) )*e( j )
345  DO 50 jr = 1, j - 1
346  work( j*n+jr ) = -tau( j )*e( j )*v( jr, j+1 )
347  50 CONTINUE
348  END IF
349 *
350  vsave = v( j, j+1 )
351  v( j, j+1 ) = one
352  CALL clarfy( 'U', j, v( 1, j+1 ), 1, tau( j ), work, n,
353  $ work( n**2+1 ) )
354  v( j, j+1 ) = vsave
355  work( ( n+1 )*j+1 ) = d( j+1 )
356  60 CONTINUE
357  END IF
358 *
359  DO 90 jcol = 1, n
360  IF( lower ) THEN
361  DO 70 jrow = jcol, n
362  work( jrow+n*( jcol-1 ) ) = work( jrow+n*( jcol-1 ) )
363  $ - a( jrow, jcol )
364  70 CONTINUE
365  ELSE
366  DO 80 jrow = 1, jcol
367  work( jrow+n*( jcol-1 ) ) = work( jrow+n*( jcol-1 ) )
368  $ - a( jrow, jcol )
369  80 CONTINUE
370  END IF
371  90 CONTINUE
372  wnorm = clanhe( '1', cuplo, n, work, n, rwork )
373 *
374  ELSE IF( itype.EQ.3 ) THEN
375 *
376 * ITYPE=3: error = U V**H - I
377 *
378  IF( n.LT.2 )
379  $ RETURN
380  CALL clacpy( ' ', n, n, u, ldu, work, n )
381  IF( lower ) THEN
382  CALL cunm2r( 'R', 'C', n, n-1, n-1, v( 2, 1 ), ldv, tau,
383  $ work( n+1 ), n, work( n**2+1 ), iinfo )
384  ELSE
385  CALL cunm2l( 'R', 'C', n, n-1, n-1, v( 1, 2 ), ldv, tau,
386  $ work, n, work( n**2+1 ), iinfo )
387  END IF
388  IF( iinfo.NE.0 ) THEN
389  result( 1 ) = ten / ulp
390  RETURN
391  END IF
392 *
393  DO 100 j = 1, n
394  work( ( n+1 )*( j-1 )+1 ) = work( ( n+1 )*( j-1 )+1 ) - cone
395  100 CONTINUE
396 *
397  wnorm = clange( '1', n, n, work, n, rwork )
398  END IF
399 *
400  IF( anorm.GT.wnorm ) THEN
401  result( 1 ) = ( wnorm / anorm ) / ( n*ulp )
402  ELSE
403  IF( anorm.LT.one ) THEN
404  result( 1 ) = ( min( wnorm, n*anorm ) / anorm ) / ( n*ulp )
405  ELSE
406  result( 1 ) = min( wnorm / anorm, real( n ) ) / ( n*ulp )
407  END IF
408  END IF
409 *
410 * Do Test 2
411 *
412 * Compute U U**H - I
413 *
414  IF( itype.EQ.1 ) THEN
415  CALL cgemm( 'N', 'C', n, n, n, cone, u, ldu, u, ldu, czero,
416  $ work, n )
417 *
418  DO 110 j = 1, n
419  work( ( n+1 )*( j-1 )+1 ) = work( ( n+1 )*( j-1 )+1 ) - cone
420  110 CONTINUE
421 *
422  result( 2 ) = min( clange( '1', n, n, work, n, rwork ),
423  $ real( n ) ) / ( n*ulp )
424  END IF
425 *
426  RETURN
427 *
428 * End of CHET21
429 *
Here is the call graph for this function:
Here is the caller graph for this function:
clange
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:117
cgemm
subroutine cgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
CGEMM
Definition: cgemm.f:189
cher
subroutine cher(UPLO, N, ALPHA, X, INCX, A, LDA)
CHER
Definition: cher.f:137
clarfy
subroutine clarfy(UPLO, N, V, INCV, TAU, C, LDC, WORK)
CLARFY
Definition: clarfy.f:110
clacpy
subroutine clacpy(UPLO, M, N, A, LDA, B, LDB)
CLACPY copies all or part of one two-dimensional array to another.
Definition: clacpy.f:105
cunm2l
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:161
clanhe
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:126
cher2
subroutine cher2(UPLO, N, ALPHA, X, INCX, Y, INCY, A, LDA)
CHER2
Definition: cher2.f:152
lsame
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:55
claset
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:108
slamch
real function slamch(CMACH)
SLAMCH
Definition: slamch.f:70
cunm2r
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:161