LAPACK  3.6.1
LAPACK: Linear Algebra PACKage
subroutine zhet21 ( integer  ITYPE,
character  UPLO,
integer  N,
integer  KBAND,
complex*16, dimension( lda, * )  A,
integer  LDA,
double precision, dimension( * )  D,
double precision, dimension( * )  E,
complex*16, dimension( ldu, * )  U,
integer  LDU,
complex*16, dimension( ldv, * )  V,
integer  LDV,
complex*16, dimension( * )  TAU,
complex*16, dimension( * )  WORK,
double precision, dimension( * )  RWORK,
double precision, dimension( 2 )  RESULT 
)

ZHET21

Purpose: 
 ZHET21 generally checks a decomposition of the form

    A = U S UC>
 where * 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* | / ( |A| n ulp ) *andC>    RESULT(2) = | I - UU* | / ( n ulp )

 If ITYPE=2, then:

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

 If ITYPE=3, then:

    RESULT(1) = | I - UV* | / ( 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)C> 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* | / ( |A| n ulp )   *andC>             RESULT(2) = | I - UU* | / ( n ulp )

          2: U expressed as a product V of Housholder transformations:
             RESULT(1) = | A - V S V* | / ( |A| n ulp )

          3: U expressed both as a dense unitary matrix and
             as a product of Housholder transformations:
             RESULT(1) = | I - UV* | / ( 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, ZHET21 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*16 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 DOUBLE PRECISION array, dimension (N)
          The diagonal of the (symmetric tri-) diagonal matrix.
[in]E
          E is DOUBLE PRECISION 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*16 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*16 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*16 array, dimension (N)
          If ITYPE >= 2, then TAU(j) is the scalar factor of
          v(j) v(j)* 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*16 array, dimension (2*N**2)
[out]RWORK
          RWORK is DOUBLE PRECISION array, dimension (N)
[out]RESULT
          RESULT is DOUBLE PRECISION 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
November 2011

Definition at line 213 of file zhet21.f.

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  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*
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 = 1, 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* - 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* - 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 UU* - 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 *
zlacpy
subroutine zlacpy(UPLO, M, N, A, LDA, B, LDB)
ZLACPY copies all or part of one two-dimensional array to another.
Definition: zlacpy.f:105
dlamch
double precision function dlamch(CMACH)
DLAMCH
Definition: dlamch.f:65
zlanhe
double precision function zlanhe(NORM, UPLO, N, A, LDA, WORK)
ZLANHE returns the value of the 1-norm, or the Frobenius norm, or the infinity norm, or the element of largest absolute value of a complex Hermitian matrix.
Definition: zlanhe.f:126
zunm2l
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:161
zgemm
subroutine zgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
ZGEMM
Definition: zgemm.f:189
zlaset
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:108
zher2
subroutine zher2(UPLO, N, ALPHA, X, INCX, Y, INCY, A, LDA)
ZHER2
Definition: zher2.f:152
zlange
double precision function zlange(NORM, M, N, A, LDA, WORK)
ZLANGE returns the value of the 1-norm, Frobenius norm, infinity-norm, or the largest absolute value ...
Definition: zlange.f:117
zunm2r
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:161
zher
subroutine zher(UPLO, N, ALPHA, X, INCX, A, LDA)
ZHER
Definition: zher.f:137
zlarfy
subroutine zlarfy(UPLO, N, V, INCV, TAU, C, LDC, WORK)
ZLARFY
Definition: zlarfy.f:110
lsame
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:55

Here is the call graph for this function:

Here is the caller graph for this function: