226 SUBROUTINE zhpt21( ITYPE, UPLO, N, KBAND, AP, D, E, U, LDU, VP,
227 $ TAU, WORK, RWORK, RESULT )
235 INTEGER ITYPE, KBAND, LDU, N
238 DOUBLE PRECISION D( * ), E( * ), RESULT( 2 ), RWORK( * )
239 COMPLEX*16 AP( * ), TAU( * ), U( LDU, * ), VP( * ),
246 DOUBLE PRECISION ZERO, ONE, TEN
247 parameter( zero = 0.0d+0, one = 1.0d+0, ten = 10.0d+0 )
248 DOUBLE PRECISION HALF
249 parameter( half = 1.0d+0 / 2.0d+0 )
250 COMPLEX*16 CZERO, CONE
251 parameter( czero = ( 0.0d+0, 0.0d+0 ),
252 $ cone = ( 1.0d+0, 0.0d+0 ) )
257 INTEGER IINFO, J, JP, JP1, JR, LAP
258 DOUBLE PRECISION ANORM, ULP, UNFL, WNORM
259 COMPLEX*16 TEMP, VSAVE
263 DOUBLE PRECISION DLAMCH, ZLANGE, ZLANHP
265 EXTERNAL lsame, dlamch, zlange, zlanhp, zdotc
272 INTRINSIC dble, dcmplx, max, min
284 lap = ( n*( n+1 ) ) / 2
286 IF( lsame( uplo,
'U' ) )
THEN
294 unfl = dlamch(
'Safe minimum' )
295 ulp = dlamch(
'Epsilon' )*dlamch(
'Base' )
299 IF( itype.LT.1 .OR. itype.GT.3 )
THEN
300 result( 1 ) = ten / ulp
308 IF( itype.EQ.3 )
THEN
311 anorm = max( zlanhp(
'1', cuplo, n, ap, rwork ), unfl )
316 IF( itype.EQ.1 )
THEN
320 CALL zlaset(
'Full', n, n, czero, czero, work, n )
321 CALL zcopy( lap, ap, 1, work, 1 )
324 CALL zhpr( cuplo, n, -d( j ), u( 1, j ), 1, work )
327 IF( n.GT.1 .AND. kband.EQ.1 )
THEN
329 CALL zhpr2( cuplo, n, -dcmplx( e( j ) ), u( 1, j ), 1,
330 $ u( 1, j-1 ), 1, work )
333 wnorm = zlanhp(
'1', cuplo, n, work, rwork )
335 ELSE IF( itype.EQ.2 )
THEN
339 CALL zlaset(
'Full', n, n, czero, czero, work, n )
343 DO 40 j = n - 1, 1, -1
344 jp = ( ( 2*n-j )*( j-1 ) ) / 2
346 IF( kband.EQ.1 )
THEN
347 work( jp+j+1 ) = ( cone-tau( j ) )*e( j )
349 work( jp+jr ) = -tau( j )*e( j )*vp( jp+jr )
353 IF( tau( j ).NE.czero )
THEN
356 CALL zhpmv(
'L', n-j, cone, work( jp1+j+1 ),
357 $ vp( jp+j+1 ), 1, czero, work( lap+1 ), 1 )
358 temp = -half*tau( j )*zdotc( n-j, work( lap+1 ), 1,
360 CALL zaxpy( n-j, temp, vp( jp+j+1 ), 1, work( lap+1 ),
362 CALL zhpr2(
'L', n-j, -tau( j ), vp( jp+j+1 ), 1,
363 $ work( lap+1 ), 1, work( jp1+j+1 ) )
367 work( jp+j ) = d( j )
372 jp = ( j*( j-1 ) ) / 2
374 IF( kband.EQ.1 )
THEN
375 work( jp1+j ) = ( cone-tau( j ) )*e( j )
377 work( jp1+jr ) = -tau( j )*e( j )*vp( jp1+jr )
381 IF( tau( j ).NE.czero )
THEN
384 CALL zhpmv(
'U', j, cone, work, vp( jp1+1 ), 1, czero,
386 temp = -half*tau( j )*zdotc( j, work( lap+1 ), 1,
388 CALL zaxpy( j, temp, vp( jp1+1 ), 1, work( lap+1 ),
390 CALL zhpr2(
'U', j, -tau( j ), vp( jp1+1 ), 1,
391 $ work( lap+1 ), 1, work )
394 work( jp1+j+1 ) = d( j+1 )
399 work( j ) = work( j ) - ap( j )
401 wnorm = zlanhp(
'1', cuplo, n, work, rwork )
403 ELSE IF( itype.EQ.3 )
THEN
409 CALL zlacpy(
' ', n, n, u, ldu, work, n )
410 CALL zupmtr(
'R', cuplo,
'C', n, n, vp, tau, work, n,
411 $ work( n**2+1 ), iinfo )
412 IF( iinfo.NE.0 )
THEN
413 result( 1 ) = ten / ulp
418 work( ( n+1 )*( j-1 )+1 ) = work( ( n+1 )*( j-1 )+1 ) - cone
421 wnorm = zlange(
'1', n, n, work, n, rwork )
424 IF( anorm.GT.wnorm )
THEN
425 result( 1 ) = ( wnorm / anorm ) / ( n*ulp )
427 IF( anorm.LT.one )
THEN
428 result( 1 ) = ( min( wnorm, n*anorm ) / anorm ) / ( n*ulp )
430 result( 1 ) = min( wnorm / anorm, dble( n ) ) / ( n*ulp )
438 IF( itype.EQ.1 )
THEN
439 CALL zgemm(
'N',
'C', n, n, n, cone, u, ldu, u, ldu, czero,
443 work( ( n+1 )*( j-1 )+1 ) = work( ( n+1 )*( j-1 )+1 ) - cone
446 result( 2 ) = min( zlange(
'1', n, n, work, n, rwork ),
447 $ dble( n ) ) / ( n*ulp )
subroutine zaxpy(n, za, zx, incx, zy, incy)
ZAXPY
subroutine zcopy(n, zx, incx, zy, incy)
ZCOPY
subroutine zgemm(transa, transb, m, n, k, alpha, a, lda, b, ldb, beta, c, ldc)
ZGEMM
subroutine zhpmv(uplo, n, alpha, ap, x, incx, beta, y, incy)
ZHPMV
subroutine zhpr2(uplo, n, alpha, x, incx, y, incy, ap)
ZHPR2
subroutine zhpr(uplo, n, alpha, x, incx, ap)
ZHPR
subroutine zlacpy(uplo, m, n, a, lda, b, ldb)
ZLACPY copies all or part of one two-dimensional array to another.
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.
subroutine zupmtr(side, uplo, trans, m, n, ap, tau, c, ldc, work, info)
ZUPMTR
subroutine zhpt21(itype, uplo, n, kband, ap, d, e, u, ldu, vp, tau, work, rwork, result)
ZHPT21