LAPACK: TESTING/EIG/zhbt21.f Source File

Go to the documentation of this file.
00001       SUBROUTINE ZHBT21( UPLO, N, KA, KS, A, LDA, D, E, U, LDU, WORK,
00002      $                   RWORK, RESULT )
00003 *
00004 *  -- LAPACK test routine (version 3.1) --
00005 *     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
00006 *     November 2006
00007 *
00008 *     .. Scalar Arguments ..
00009       CHARACTER          UPLO
00010       INTEGER            KA, KS, LDA, LDU, N
00011 *     ..
00012 *     .. Array Arguments ..
00013       DOUBLE PRECISION   D( * ), E( * ), RESULT( 2 ), RWORK( * )
00014       COMPLEX*16         A( LDA, * ), U( LDU, * ), WORK( * )
00015 *     ..
00016 *
00017 *  Purpose
00018 *  =======
00019 *
00020 *  ZHBT21  generally checks a decomposition of the form
00021 *
00022 *          A = U S U*
00023 *
00024 *  where * means conjugate transpose, A is hermitian banded, U is
00025 *  unitary, and S is diagonal (if KS=0) or symmetric
00026 *  tridiagonal (if KS=1).
00027 *
00028 *  Specifically:
00029 *
00030 *          RESULT(1) = | A - U S U* | / ( |A| n ulp ) *and*
00031 *          RESULT(2) = | I - UU* | / ( n ulp )
00032 *
00033 *  Arguments
00034 *  =========
00035 *
00036 *  UPLO    (input) CHARACTER
00037 *          If UPLO='U', the upper triangle of A and V will be used and
00038 *          the (strictly) lower triangle will not be referenced.
00039 *          If UPLO='L', the lower triangle of A and V will be used and
00040 *          the (strictly) upper triangle will not be referenced.
00041 *
00042 *  N       (input) INTEGER
00043 *          The size of the matrix.  If it is zero, ZHBT21 does nothing.
00044 *          It must be at least zero.
00045 *
00046 *  KA      (input) INTEGER
00047 *          The bandwidth of the matrix A.  It must be at least zero.  If
00048 *          it is larger than N-1, then max( 0, N-1 ) will be used.
00049 *
00050 *  KS      (input) INTEGER
00051 *          The bandwidth of the matrix S.  It may only be zero or one.
00052 *          If zero, then S is diagonal, and E is not referenced.  If
00053 *          one, then S is symmetric tri-diagonal.
00054 *
00055 *  A       (input) COMPLEX*16 array, dimension (LDA, N)
00056 *          The original (unfactored) matrix.  It is assumed to be
00057 *          hermitian, and only the upper (UPLO='U') or only the lower
00058 *          (UPLO='L') will be referenced.
00059 *
00060 *  LDA     (input) INTEGER
00061 *          The leading dimension of A.  It must be at least 1
00062 *          and at least min( KA, N-1 ).
00063 *
00064 *  D       (input) DOUBLE PRECISION array, dimension (N)
00065 *          The diagonal of the (symmetric tri-) diagonal matrix S.
00066 *
00067 *  E       (input) DOUBLE PRECISION array, dimension (N-1)
00068 *          The off-diagonal of the (symmetric tri-) diagonal matrix S.
00069 *          E(1) is the (1,2) and (2,1) element, E(2) is the (2,3) and
00070 *          (3,2) element, etc.
00071 *          Not referenced if KS=0.
00072 *
00073 *  U       (input) COMPLEX*16 array, dimension (LDU, N)
00074 *          The unitary matrix in the decomposition, expressed as a
00075 *          dense matrix (i.e., not as a product of Householder
00076 *          transformations, Givens transformations, etc.)
00077 *
00078 *  LDU     (input) INTEGER
00079 *          The leading dimension of U.  LDU must be at least N and
00080 *          at least 1.
00081 *
00082 *  WORK    (workspace) COMPLEX*16 array, dimension (N**2)
00083 *
00084 *  RWORK   (workspace) DOUBLE PRECISION array, dimension (N)
00085 *
00086 *  RESULT  (output) DOUBLE PRECISION array, dimension (2)
00087 *          The values computed by the two tests described above.  The
00088 *          values are currently limited to 1/ulp, to avoid overflow.
00089 *
00090 *  =====================================================================
00091 *
00092 *     .. Parameters ..
00093       COMPLEX*16         CZERO, CONE
00094       PARAMETER          ( CZERO = ( 0.0D+0, 0.0D+0 ),
00095      $                   CONE = ( 1.0D+0, 0.0D+0 ) )
00096       DOUBLE PRECISION   ZERO, ONE
00097       PARAMETER          ( ZERO = 0.0D+0, ONE = 1.0D+0 )
00098 *     ..
00099 *     .. Local Scalars ..
00100       LOGICAL            LOWER
00101       CHARACTER          CUPLO
00102       INTEGER            IKA, J, JC, JR
00103       DOUBLE PRECISION   ANORM, ULP, UNFL, WNORM
00104 *     ..
00105 *     .. External Functions ..
00106       LOGICAL            LSAME
00107       DOUBLE PRECISION   DLAMCH, ZLANGE, ZLANHB, ZLANHP
00108       EXTERNAL           LSAME, DLAMCH, ZLANGE, ZLANHB, ZLANHP
00109 *     ..
00110 *     .. External Subroutines ..
00111       EXTERNAL           ZGEMM, ZHPR, ZHPR2
00112 *     ..
00113 *     .. Intrinsic Functions ..
00114       INTRINSIC          DBLE, DCMPLX, MAX, MIN
00115 *     ..
00116 *     .. Executable Statements ..
00117 *
00118 *     Constants
00119 *
00120       RESULT( 1 ) = ZERO
00121       RESULT( 2 ) = ZERO
00122       IF( N.LE.0 )
00123      $   RETURN
00124 *
00125       IKA = MAX( 0, MIN( N-1, KA ) )
00126 *
00127       IF( LSAME( UPLO, 'U' ) ) THEN
00128          LOWER = .FALSE.
00129          CUPLO = 'U'
00130       ELSE
00131          LOWER = .TRUE.
00132          CUPLO = 'L'
00133       END IF
00134 *
00135       UNFL = DLAMCH( 'Safe minimum' )
00136       ULP = DLAMCH( 'Epsilon' )*DLAMCH( 'Base' )
00137 *
00138 *     Some Error Checks
00139 *
00140 *     Do Test 1
00141 *
00142 *     Norm of A:
00143 *
00144       ANORM = MAX( ZLANHB( '1', CUPLO, N, IKA, A, LDA, RWORK ), UNFL )
00145 *
00146 *     Compute error matrix:    Error = A - U S U*
00147 *
00148 *     Copy A from SB to SP storage format.
00149 *
00150       J = 0
00151       DO 50 JC = 1, N
00152          IF( LOWER ) THEN
00153             DO 10 JR = 1, MIN( IKA+1, N+1-JC )
00154                J = J + 1
00155                WORK( J ) = A( JR, JC )
00156    10       CONTINUE
00157             DO 20 JR = IKA + 2, N + 1 - JC
00158                J = J + 1
00159                WORK( J ) = ZERO
00160    20       CONTINUE
00161          ELSE
00162             DO 30 JR = IKA + 2, JC
00163                J = J + 1
00164                WORK( J ) = ZERO
00165    30       CONTINUE
00166             DO 40 JR = MIN( IKA, JC-1 ), 0, -1
00167                J = J + 1
00168                WORK( J ) = A( IKA+1-JR, JC )
00169    40       CONTINUE
00170          END IF
00171    50 CONTINUE
00172 *
00173       DO 60 J = 1, N
00174          CALL ZHPR( CUPLO, N, -D( J ), U( 1, J ), 1, WORK )
00175    60 CONTINUE
00176 *
00177       IF( N.GT.1 .AND. KS.EQ.1 ) THEN
00178          DO 70 J = 1, N - 1
00179             CALL ZHPR2( CUPLO, N, -DCMPLX( E( J ) ), U( 1, J ), 1,
00180      $                  U( 1, J+1 ), 1, WORK )
00181    70    CONTINUE
00182       END IF
00183       WNORM = ZLANHP( '1', CUPLO, N, WORK, RWORK )
00184 *
00185       IF( ANORM.GT.WNORM ) THEN
00186          RESULT( 1 ) = ( WNORM / ANORM ) / ( N*ULP )
00187       ELSE
00188          IF( ANORM.LT.ONE ) THEN
00189             RESULT( 1 ) = ( MIN( WNORM, N*ANORM ) / ANORM ) / ( N*ULP )
00190          ELSE
00191             RESULT( 1 ) = MIN( WNORM / ANORM, DBLE( N ) ) / ( N*ULP )
00192          END IF
00193       END IF
00194 *
00195 *     Do Test 2
00196 *
00197 *     Compute  UU* - I
00198 *
00199       CALL ZGEMM( 'N', 'C', N, N, N, CONE, U, LDU, U, LDU, CZERO, WORK,
00200      $            N )
00201 *
00202       DO 80 J = 1, N
00203          WORK( ( N+1 )*( J-1 )+1 ) = WORK( ( N+1 )*( J-1 )+1 ) - CONE
00204    80 CONTINUE
00205 *
00206       RESULT( 2 ) = MIN( ZLANGE( '1', N, N, WORK, N, RWORK ),
00207      $              DBLE( N ) ) / ( N*ULP )
00208 *
00209       RETURN
00210 *
00211 *     End of ZHBT21
00212 *
00213       END