LAPACK 3.11.0
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches

◆ 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.

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: