LAPACK: TESTING/EIG/dstt22.f Source File

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