LAPACK: TESTING/EIG/zstt22.f Source File

Go to the documentation of this file.
00001       SUBROUTINE ZSTT22( N, M, KBAND, AD, AE, SD, SE, U, LDU, WORK,
00002      $                   LDWORK, 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       INTEGER            KBAND, LDU, LDWORK, M, N
00010 *     ..
00011 *     .. Array Arguments ..
00012       DOUBLE PRECISION   AD( * ), AE( * ), RESULT( 2 ), RWORK( * ),
00013      $                   SD( * ), SE( * )
00014       COMPLEX*16         U( LDU, * ), WORK( LDWORK, * )
00015 *     ..
00016 *
00017 *  Purpose
00018 *  =======
00019 *
00020 *  ZSTT22  checks a set of M eigenvalues and eigenvectors,
00021 *
00022 *      A U = U S
00023 *
00024 *  where A is Hermitian tridiagonal, the columns of U are unitary,
00025 *  and S is diagonal (if KBAND=0) or Hermitian tridiagonal (if KBAND=1).
00026 *  Two tests are performed:
00027 *
00028 *     RESULT(1) = | U* A U - S | / ( |A| m ulp )
00029 *
00030 *     RESULT(2) = | I - U*U | / ( m ulp )
00031 *
00032 *  Arguments
00033 *  =========
00034 *
00035 *  N       (input) INTEGER
00036 *          The size of the matrix.  If it is zero, ZSTT22 does nothing.
00037 *          It must be at least zero.
00038 *
00039 *  M       (input) INTEGER
00040 *          The number of eigenpairs to check.  If it is zero, ZSTT22
00041 *          does nothing.  It must be at least zero.
00042 *
00043 *  KBAND   (input) INTEGER
00044 *          The bandwidth of the matrix S.  It may only be zero or one.
00045 *          If zero, then S is diagonal, and SE is not referenced.  If
00046 *          one, then S is Hermitian tri-diagonal.
00047 *
00048 *  AD      (input) DOUBLE PRECISION array, dimension (N)
00049 *          The diagonal of the original (unfactored) matrix A.  A is
00050 *          assumed to be Hermitian tridiagonal.
00051 *
00052 *  AE      (input) DOUBLE PRECISION array, dimension (N)
00053 *          The off-diagonal of the original (unfactored) matrix A.  A
00054 *          is assumed to be Hermitian tridiagonal.  AE(1) is ignored,
00055 *          AE(2) is the (1,2) and (2,1) element, etc.
00056 *
00057 *  SD      (input) DOUBLE PRECISION array, dimension (N)
00058 *          The diagonal of the (Hermitian tri-) diagonal matrix S.
00059 *
00060 *  SE      (input) DOUBLE PRECISION array, dimension (N)
00061 *          The off-diagonal of the (Hermitian tri-) diagonal matrix S.
00062 *          Not referenced if KBSND=0.  If KBAND=1, then AE(1) is
00063 *          ignored, SE(2) is the (1,2) and (2,1) element, etc.
00064 *
00065 *  U       (input) DOUBLE PRECISION array, dimension (LDU, N)
00066 *          The unitary matrix in the decomposition.
00067 *
00068 *  LDU     (input) INTEGER
00069 *          The leading dimension of U.  LDU must be at least N.
00070 *
00071 *  WORK    (workspace) COMPLEX*16 array, dimension (LDWORK, M+1)
00072 *
00073 *  LDWORK  (input) INTEGER
00074 *          The leading dimension of WORK.  LDWORK must be at least
00075 *          max(1,M).
00076 *
00077 *  RWORK   (workspace) DOUBLE PRECISION array, dimension (N)
00078 *
00079 *  RESULT  (output) DOUBLE PRECISION array, dimension (2)
00080 *          The values computed by the two tests described above.  The
00081 *          values are currently limited to 1/ulp, to avoid overflow.
00082 *
00083 *  =====================================================================
00084 *
00085 *     .. Parameters ..
00086       DOUBLE PRECISION   ZERO, ONE
00087       PARAMETER          ( ZERO = 0.0D0, ONE = 1.0D0 )
00088       COMPLEX*16         CZERO, CONE
00089       PARAMETER          ( CZERO = ( 0.0D+0, 0.0D+0 ),
00090      $                   CONE = ( 1.0D+0, 0.0D+0 ) )
00091 *     ..
00092 *     .. Local Scalars ..
00093       INTEGER            I, J, K
00094       DOUBLE PRECISION   ANORM, ULP, UNFL, WNORM
00095       COMPLEX*16         AUKJ
00096 *     ..
00097 *     .. External Functions ..
00098       DOUBLE PRECISION   DLAMCH, ZLANGE, ZLANSY
00099       EXTERNAL           DLAMCH, ZLANGE, ZLANSY
00100 *     ..
00101 *     .. External Subroutines ..
00102       EXTERNAL           ZGEMM
00103 *     ..
00104 *     .. Intrinsic Functions ..
00105       INTRINSIC          ABS, DBLE, MAX, MIN
00106 *     ..
00107 *     .. Executable Statements ..
00108 *
00109       RESULT( 1 ) = ZERO
00110       RESULT( 2 ) = ZERO
00111       IF( N.LE.0 .OR. M.LE.0 )
00112      $   RETURN
00113 *
00114       UNFL = DLAMCH( 'Safe minimum' )
00115       ULP = DLAMCH( 'Epsilon' )
00116 *
00117 *     Do Test 1
00118 *
00119 *     Compute the 1-norm of A.
00120 *
00121       IF( N.GT.1 ) THEN
00122          ANORM = ABS( AD( 1 ) ) + ABS( AE( 1 ) )
00123          DO 10 J = 2, N - 1
00124             ANORM = MAX( ANORM, ABS( AD( J ) )+ABS( AE( J ) )+
00125      $              ABS( AE( J-1 ) ) )
00126    10    CONTINUE
00127          ANORM = MAX( ANORM, ABS( AD( N ) )+ABS( AE( N-1 ) ) )
00128       ELSE
00129          ANORM = ABS( AD( 1 ) )
00130       END IF
00131       ANORM = MAX( ANORM, UNFL )
00132 *
00133 *     Norm of U*AU - S
00134 *
00135       DO 40 I = 1, M
00136          DO 30 J = 1, M
00137             WORK( I, J ) = CZERO
00138             DO 20 K = 1, N
00139                AUKJ = AD( K )*U( K, J )
00140                IF( K.NE.N )
00141      $            AUKJ = AUKJ + AE( K )*U( K+1, J )
00142                IF( K.NE.1 )
00143      $            AUKJ = AUKJ + AE( K-1 )*U( K-1, J )
00144                WORK( I, J ) = WORK( I, J ) + U( K, I )*AUKJ
00145    20       CONTINUE
00146    30    CONTINUE
00147          WORK( I, I ) = WORK( I, I ) - SD( I )
00148          IF( KBAND.EQ.1 ) THEN
00149             IF( I.NE.1 )
00150      $         WORK( I, I-1 ) = WORK( I, I-1 ) - SE( I-1 )
00151             IF( I.NE.N )
00152      $         WORK( I, I+1 ) = WORK( I, I+1 ) - SE( I )
00153          END IF
00154    40 CONTINUE
00155 *
00156       WNORM = ZLANSY( '1', 'L', M, WORK, M, RWORK )
00157 *
00158       IF( ANORM.GT.WNORM ) THEN
00159          RESULT( 1 ) = ( WNORM / ANORM ) / ( M*ULP )
00160       ELSE
00161          IF( ANORM.LT.ONE ) THEN
00162             RESULT( 1 ) = ( MIN( WNORM, M*ANORM ) / ANORM ) / ( M*ULP )
00163          ELSE
00164             RESULT( 1 ) = MIN( WNORM / ANORM, DBLE( M ) ) / ( M*ULP )
00165          END IF
00166       END IF
00167 *
00168 *     Do Test 2
00169 *
00170 *     Compute  U*U - I
00171 *
00172       CALL ZGEMM( 'T', 'N', M, M, N, CONE, U, LDU, U, LDU, CZERO, WORK,
00173      $            M )
00174 *
00175       DO 50 J = 1, M
00176          WORK( J, J ) = WORK( J, J ) - ONE
00177    50 CONTINUE
00178 *
00179       RESULT( 2 ) = MIN( DBLE( M ), ZLANGE( '1', M, M, WORK, M,
00180      $              RWORK ) ) / ( M*ULP )
00181 *
00182       RETURN
00183 *
00184 *     End of ZSTT22
00185 *
00186       END