LAPACK  3.4.1
LAPACK: Linear Algebra PACKage
cstt22.f
Go to the documentation of this file.
00001 *> \brief \b CSTT22
00002 *
00003 *  =========== DOCUMENTATION ===========
00004 *
00005 * Online html documentation available at 
00006 *            http://www.netlib.org/lapack/explore-html/ 
00007 *
00008 *  Definition:
00009 *  ===========
00010 *
00011 *       SUBROUTINE CSTT22( N, M, KBAND, AD, AE, SD, SE, U, LDU, WORK,
00012 *                          LDWORK, RWORK, RESULT )
00013 * 
00014 *       .. Scalar Arguments ..
00015 *       INTEGER            KBAND, LDU, LDWORK, M, N
00016 *       ..
00017 *       .. Array Arguments ..
00018 *       REAL               AD( * ), AE( * ), RESULT( 2 ), RWORK( * ),
00019 *      $                   SD( * ), SE( * )
00020 *       COMPLEX            U( LDU, * ), WORK( LDWORK, * )
00021 *       ..
00022 *  
00023 *
00024 *> \par Purpose:
00025 *  =============
00026 *>
00027 *> \verbatim
00028 *>
00029 *> CSTT22  checks a set of M eigenvalues and eigenvectors,
00030 *>
00031 *>     A U = U S
00032 *>
00033 *> where A is Hermitian tridiagonal, the columns of U are unitary,
00034 *> and S is diagonal (if KBAND=0) or Hermitian tridiagonal (if KBAND=1).
00035 *> Two tests are performed:
00036 *>
00037 *>    RESULT(1) = | U* A U - S | / ( |A| m ulp )
00038 *>
00039 *>    RESULT(2) = | I - U*U | / ( m ulp )
00040 *> \endverbatim
00041 *
00042 *  Arguments:
00043 *  ==========
00044 *
00045 *> \param[in] N
00046 *> \verbatim
00047 *>          N is INTEGER
00048 *>          The size of the matrix.  If it is zero, CSTT22 does nothing.
00049 *>          It must be at least zero.
00050 *> \endverbatim
00051 *>
00052 *> \param[in] M
00053 *> \verbatim
00054 *>          M is INTEGER
00055 *>          The number of eigenpairs to check.  If it is zero, CSTT22
00056 *>          does nothing.  It must be at least zero.
00057 *> \endverbatim
00058 *>
00059 *> \param[in] KBAND
00060 *> \verbatim
00061 *>          KBAND is INTEGER
00062 *>          The bandwidth of the matrix S.  It may only be zero or one.
00063 *>          If zero, then S is diagonal, and SE is not referenced.  If
00064 *>          one, then S is Hermitian tri-diagonal.
00065 *> \endverbatim
00066 *>
00067 *> \param[in] AD
00068 *> \verbatim
00069 *>          AD is REAL array, dimension (N)
00070 *>          The diagonal of the original (unfactored) matrix A.  A is
00071 *>          assumed to be Hermitian tridiagonal.
00072 *> \endverbatim
00073 *>
00074 *> \param[in] AE
00075 *> \verbatim
00076 *>          AE is REAL array, dimension (N)
00077 *>          The off-diagonal of the original (unfactored) matrix A.  A
00078 *>          is assumed to be Hermitian tridiagonal.  AE(1) is ignored,
00079 *>          AE(2) is the (1,2) and (2,1) element, etc.
00080 *> \endverbatim
00081 *>
00082 *> \param[in] SD
00083 *> \verbatim
00084 *>          SD is REAL array, dimension (N)
00085 *>          The diagonal of the (Hermitian tri-) diagonal matrix S.
00086 *> \endverbatim
00087 *>
00088 *> \param[in] SE
00089 *> \verbatim
00090 *>          SE is REAL array, dimension (N)
00091 *>          The off-diagonal of the (Hermitian tri-) diagonal matrix S.
00092 *>          Not referenced if KBSND=0.  If KBAND=1, then AE(1) is
00093 *>          ignored, SE(2) is the (1,2) and (2,1) element, etc.
00094 *> \endverbatim
00095 *>
00096 *> \param[in] U
00097 *> \verbatim
00098 *>          U is REAL array, dimension (LDU, N)
00099 *>          The unitary matrix in the decomposition.
00100 *> \endverbatim
00101 *>
00102 *> \param[in] LDU
00103 *> \verbatim
00104 *>          LDU is INTEGER
00105 *>          The leading dimension of U.  LDU must be at least N.
00106 *> \endverbatim
00107 *>
00108 *> \param[out] WORK
00109 *> \verbatim
00110 *>          WORK is COMPLEX array, dimension (LDWORK, M+1)
00111 *> \endverbatim
00112 *>
00113 *> \param[in] LDWORK
00114 *> \verbatim
00115 *>          LDWORK is INTEGER
00116 *>          The leading dimension of WORK.  LDWORK must be at least
00117 *>          max(1,M).
00118 *> \endverbatim
00119 *>
00120 *> \param[out] RWORK
00121 *> \verbatim
00122 *>          RWORK is REAL array, dimension (N)
00123 *> \endverbatim
00124 *>
00125 *> \param[out] RESULT
00126 *> \verbatim
00127 *>          RESULT is REAL array, dimension (2)
00128 *>          The values computed by the two tests described above.  The
00129 *>          values are currently limited to 1/ulp, to avoid overflow.
00130 *> \endverbatim
00131 *
00132 *  Authors:
00133 *  ========
00134 *
00135 *> \author Univ. of Tennessee 
00136 *> \author Univ. of California Berkeley 
00137 *> \author Univ. of Colorado Denver 
00138 *> \author NAG Ltd. 
00139 *
00140 *> \date November 2011
00141 *
00142 *> \ingroup complex_eig
00143 *
00144 *  =====================================================================
00145       SUBROUTINE CSTT22( N, M, KBAND, AD, AE, SD, SE, U, LDU, WORK,
00146      $                   LDWORK, RWORK, RESULT )
00147 *
00148 *  -- LAPACK test routine (version 3.4.0) --
00149 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
00150 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
00151 *     November 2011
00152 *
00153 *     .. Scalar Arguments ..
00154       INTEGER            KBAND, LDU, LDWORK, M, N
00155 *     ..
00156 *     .. Array Arguments ..
00157       REAL               AD( * ), AE( * ), RESULT( 2 ), RWORK( * ),
00158      $                   SD( * ), SE( * )
00159       COMPLEX            U( LDU, * ), WORK( LDWORK, * )
00160 *     ..
00161 *
00162 *  =====================================================================
00163 *
00164 *     .. Parameters ..
00165       REAL               ZERO, ONE
00166       PARAMETER          ( ZERO = 0.0E0, ONE = 1.0E0 )
00167       COMPLEX            CZERO, CONE
00168       PARAMETER          ( CZERO = ( 0.0E+0, 0.0E+0 ),
00169      $                   CONE = ( 1.0E+0, 0.0E+0 ) )
00170 *     ..
00171 *     .. Local Scalars ..
00172       INTEGER            I, J, K
00173       REAL               ANORM, ULP, UNFL, WNORM
00174       COMPLEX            AUKJ
00175 *     ..
00176 *     .. External Functions ..
00177       REAL               CLANGE, CLANSY, SLAMCH
00178       EXTERNAL           CLANGE, CLANSY, SLAMCH
00179 *     ..
00180 *     .. External Subroutines ..
00181       EXTERNAL           CGEMM
00182 *     ..
00183 *     .. Intrinsic Functions ..
00184       INTRINSIC          ABS, MAX, MIN, REAL
00185 *     ..
00186 *     .. Executable Statements ..
00187 *
00188       RESULT( 1 ) = ZERO
00189       RESULT( 2 ) = ZERO
00190       IF( N.LE.0 .OR. M.LE.0 )
00191      $   RETURN
00192 *
00193       UNFL = SLAMCH( 'Safe minimum' )
00194       ULP = SLAMCH( 'Epsilon' )
00195 *
00196 *     Do Test 1
00197 *
00198 *     Compute the 1-norm of A.
00199 *
00200       IF( N.GT.1 ) THEN
00201          ANORM = ABS( AD( 1 ) ) + ABS( AE( 1 ) )
00202          DO 10 J = 2, N - 1
00203             ANORM = MAX( ANORM, ABS( AD( J ) )+ABS( AE( J ) )+
00204      $              ABS( AE( J-1 ) ) )
00205    10    CONTINUE
00206          ANORM = MAX( ANORM, ABS( AD( N ) )+ABS( AE( N-1 ) ) )
00207       ELSE
00208          ANORM = ABS( AD( 1 ) )
00209       END IF
00210       ANORM = MAX( ANORM, UNFL )
00211 *
00212 *     Norm of U*AU - S
00213 *
00214       DO 40 I = 1, M
00215          DO 30 J = 1, M
00216             WORK( I, J ) = CZERO
00217             DO 20 K = 1, N
00218                AUKJ = AD( K )*U( K, J )
00219                IF( K.NE.N )
00220      $            AUKJ = AUKJ + AE( K )*U( K+1, J )
00221                IF( K.NE.1 )
00222      $            AUKJ = AUKJ + AE( K-1 )*U( K-1, J )
00223                WORK( I, J ) = WORK( I, J ) + U( K, I )*AUKJ
00224    20       CONTINUE
00225    30    CONTINUE
00226          WORK( I, I ) = WORK( I, I ) - SD( I )
00227          IF( KBAND.EQ.1 ) THEN
00228             IF( I.NE.1 )
00229      $         WORK( I, I-1 ) = WORK( I, I-1 ) - SE( I-1 )
00230             IF( I.NE.N )
00231      $         WORK( I, I+1 ) = WORK( I, I+1 ) - SE( I )
00232          END IF
00233    40 CONTINUE
00234 *
00235       WNORM = CLANSY( '1', 'L', M, WORK, M, RWORK )
00236 *
00237       IF( ANORM.GT.WNORM ) THEN
00238          RESULT( 1 ) = ( WNORM / ANORM ) / ( M*ULP )
00239       ELSE
00240          IF( ANORM.LT.ONE ) THEN
00241             RESULT( 1 ) = ( MIN( WNORM, M*ANORM ) / ANORM ) / ( M*ULP )
00242          ELSE
00243             RESULT( 1 ) = MIN( WNORM / ANORM, REAL( M ) ) / ( M*ULP )
00244          END IF
00245       END IF
00246 *
00247 *     Do Test 2
00248 *
00249 *     Compute  U*U - I
00250 *
00251       CALL CGEMM( 'T', 'N', M, M, N, CONE, U, LDU, U, LDU, CZERO, WORK,
00252      $            M )
00253 *
00254       DO 50 J = 1, M
00255          WORK( J, J ) = WORK( J, J ) - ONE
00256    50 CONTINUE
00257 *
00258       RESULT( 2 ) = MIN( REAL( M ), CLANGE( '1', M, M, WORK, M,
00259      $              RWORK ) ) / ( M*ULP )
00260 *
00261       RETURN
00262 *
00263 *     End of CSTT22
00264 *
00265       END
 All Files Functions