LAPACK 3.12.1
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  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, 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 = 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**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*
cgemm
subroutine cgemm(transa, transb, m, n, k, alpha, a, lda, b, ldb, beta, c, ldc)
CGEMM
Definition cgemm.f:188
cher2
subroutine cher2(uplo, n, alpha, x, incx, y, incy, a, lda)
CHER2
Definition cher2.f:150
cher
subroutine cher(uplo, n, alpha, x, incx, a, lda)
CHER
Definition cher.f:135
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:101
slamch
real function slamch(cmach)
SLAMCH
Definition slamch.f:68
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:113
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:122
clarfy
subroutine clarfy(uplo, n, v, incv, tau, c, ldc, work)
CLARFY
Definition clarfy.f:108
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:104
lsame
logical function lsame(ca, cb)
LSAME
Definition lsame.f:48
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:157
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:157
Here is the call graph for this function:
Here is the caller graph for this function: