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

◆ zhet21()

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 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  and whose scaling constants are in .
!> We shall use the letter  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, 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)**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*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.

Definition at line 212 of file zhet21.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 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**H
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**H - 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**H - 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 U U**H - 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*
zgemm
subroutine zgemm(transa, transb, m, n, k, alpha, a, lda, b, ldb, beta, c, ldc)
ZGEMM
Definition zgemm.f:188
zher2
subroutine zher2(uplo, n, alpha, x, incx, y, incy, a, lda)
ZHER2
Definition zher2.f:150
zher
subroutine zher(uplo, n, alpha, x, incx, a, lda)
ZHER
Definition zher.f:135
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:101
dlamch
double precision function dlamch(cmach)
DLAMCH
Definition dlamch.f:69
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:113
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,...
Definition zlanhe.f:122
zlarfy
subroutine zlarfy(uplo, n, v, incv, tau, c, ldc, work)
ZLARFY
Definition zlarfy.f:108
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:104
lsame
logical function lsame(ca, cb)
LSAME
Definition lsame.f:48
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:157
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:157
Here is the call graph for this function:
Here is the caller graph for this function: