LAPACK  3.8.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 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, 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)* 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 213 of file chet21.f.

213 *
214 * -- LAPACK test routine (version 3.7.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 * December 2016
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 *
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
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
subroutine cher(UPLO, N, ALPHA, X, INCX, A, LDA)
CHER
Definition: cher.f:137
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
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:55
real function slamch(CMACH)
SLAMCH
Definition: slamch.f:69
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
subroutine cher2(UPLO, N, ALPHA, X, INCX, Y, INCY, A, LDA)
CHER2
Definition: cher2.f:152
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, or the element of largest absolute value of a complex Hermitian matrix.
Definition: clanhe.f:126
subroutine cgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
CGEMM
Definition: cgemm.f:189
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
subroutine clarfy(UPLO, N, V, INCV, TAU, C, LDC, WORK)
CLARFY
Definition: clarfy.f:110
Here is the call graph for this function:
Here is the caller graph for this function: