TESTING/EIG/sstt22.f

   1 *> \brief \b SSTT22
   2 *
   3 *  =========== DOCUMENTATION ===========
   4 *
   5 * Online html documentation available at
   6 *            http://www.netlib.org/lapack/explore-html/
   7 *
   8 *  Definition:
   9 *  ===========
  10 *
  11 *       SUBROUTINE SSTT22( N, M, KBAND, AD, AE, SD, SE, U, LDU, WORK,
  12 *                          LDWORK, RESULT )
  13 *
  14 *       .. Scalar Arguments ..
  15 *       INTEGER            KBAND, LDU, LDWORK, M, N
  16 *       ..
  17 *       .. Array Arguments ..
  18 *       REAL               AD( * ), AE( * ), RESULT( 2 ), SD( * ),
  19 *      $                   SE( * ), U( LDU, * ), WORK( LDWORK, * )
  20 *       ..
  21 *
  22 *
  23 *> \par Purpose:
  24 *  =============
  25 *>
  26 *> \verbatim
  27 *>
  28 *> SSTT22  checks a set of M eigenvalues and eigenvectors,
  29 *>
  30 *>     A U = U S
  31 *>
  32 *> where A is symmetric tridiagonal, the columns of U are orthogonal,
  33 *> and S is diagonal (if KBAND=0) or symmetric tridiagonal (if KBAND=1).
  34 *> Two tests are performed:
  35 *>
  36 *>    RESULT(1) = | U' A U - S | / ( |A| m ulp )
  37 *>
  38 *>    RESULT(2) = | I - U'U | / ( m ulp )
  39 *> \endverbatim
  40 *
  41 *  Arguments:
  42 *  ==========
  43 *
  44 *> \param[in] N
  45 *> \verbatim
  46 *>          N is INTEGER
  47 *>          The size of the matrix.  If it is zero, SSTT22 does nothing.
  48 *>          It must be at least zero.
  49 *> \endverbatim
  50 *>
  51 *> \param[in] M
  52 *> \verbatim
  53 *>          M is INTEGER
  54 *>          The number of eigenpairs to check.  If it is zero, SSTT22
  55 *>          does nothing.  It must be at least zero.
  56 *> \endverbatim
  57 *>
  58 *> \param[in] KBAND
  59 *> \verbatim
  60 *>          KBAND is INTEGER
  61 *>          The bandwidth of the matrix S.  It may only be zero or one.
  62 *>          If zero, then S is diagonal, and SE is not referenced.  If
  63 *>          one, then S is symmetric tri-diagonal.
  64 *> \endverbatim
  65 *>
  66 *> \param[in] AD
  67 *> \verbatim
  68 *>          AD is REAL array, dimension (N)
  69 *>          The diagonal of the original (unfactored) matrix A.  A is
  70 *>          assumed to be symmetric tridiagonal.
  71 *> \endverbatim
  72 *>
  73 *> \param[in] AE
  74 *> \verbatim
  75 *>          AE is REAL array, dimension (N)
  76 *>          The off-diagonal of the original (unfactored) matrix A.  A
  77 *>          is assumed to be symmetric tridiagonal.  AE(1) is ignored,
  78 *>          AE(2) is the (1,2) and (2,1) element, etc.
  79 *> \endverbatim
  80 *>
  81 *> \param[in] SD
  82 *> \verbatim
  83 *>          SD is REAL array, dimension (N)
  84 *>          The diagonal of the (symmetric tri-) diagonal matrix S.
  85 *> \endverbatim
  86 *>
  87 *> \param[in] SE
  88 *> \verbatim
  89 *>          SE is REAL array, dimension (N)
  90 *>          The off-diagonal of the (symmetric tri-) diagonal matrix S.
  91 *>          Not referenced if KBSND=0.  If KBAND=1, then AE(1) is
  92 *>          ignored, SE(2) is the (1,2) and (2,1) element, etc.
  93 *> \endverbatim
  94 *>
  95 *> \param[in] U
  96 *> \verbatim
  97 *>          U is REAL array, dimension (LDU, N)
  98 *>          The orthogonal matrix in the decomposition.
  99 *> \endverbatim
 100 *>
 101 *> \param[in] LDU
 102 *> \verbatim
 103 *>          LDU is INTEGER
 104 *>          The leading dimension of U.  LDU must be at least N.
 105 *> \endverbatim
 106 *>
 107 *> \param[out] WORK
 108 *> \verbatim
 109 *>          WORK is REAL array, dimension (LDWORK, M+1)
 110 *> \endverbatim
 111 *>
 112 *> \param[in] LDWORK
 113 *> \verbatim
 114 *>          LDWORK is INTEGER
 115 *>          The leading dimension of WORK.  LDWORK must be at least
 116 *>          max(1,M).
 117 *> \endverbatim
 118 *>
 119 *> \param[out] RESULT
 120 *> \verbatim
 121 *>          RESULT is REAL array, dimension (2)
 122 *>          The values computed by the two tests described above.  The
 123 *>          values are currently limited to 1/ulp, to avoid overflow.
 124 *> \endverbatim
 125 *
 126 *  Authors:
 127 *  ========
 128 *
 129 *> \author Univ. of Tennessee
 130 *> \author Univ. of California Berkeley
 131 *> \author Univ. of Colorado Denver
 132 *> \author NAG Ltd.
 133 *
 134 *> \date November 2011
 135 *
 136 *> \ingroup single_eig
 137 *
 138 *  =====================================================================
 139       SUBROUTINE SSTT22( N, M, KBAND, AD, AE, SD, SE, U, LDU, WORK,
 140      $                   LDWORK, RESULT )
 141 *
 142 *  -- LAPACK test routine (version 3.4.0) --
 143 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
 144 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
 145 *     November 2011
 146 *
 147 *     .. Scalar Arguments ..
 148       INTEGER            KBAND, LDU, LDWORK, M, N
 149 *     ..
 150 *     .. Array Arguments ..
 151       REAL               AD( * ), AE( * ), RESULT( 2 ), SD( * ),
 152      $                   SE( * ), U( LDU, * ), WORK( LDWORK, * )
 153 *     ..
 154 *
 155 *  =====================================================================
 156 *
 157 *     .. Parameters ..
 158       REAL               ZERO, ONE
 159       PARAMETER          ( ZERO = 0.0E0, ONE = 1.0E0 )
 160 *     ..
 161 *     .. Local Scalars ..
 162       INTEGER            I, J, K
 163       REAL               ANORM, AUKJ, ULP, UNFL, WNORM
 164 *     ..
 165 *     .. External Functions ..
 166       REAL               SLAMCH, SLANGE, SLANSY
 167       EXTERNAL           SLAMCH, SLANGE, SLANSY
 168 *     ..
 169 *     .. External Subroutines ..
 170       EXTERNAL           SGEMM
 171 *     ..
 172 *     .. Intrinsic Functions ..
 173       INTRINSIC          ABS, MAX, MIN, REAL
 174 *     ..
 175 *     .. Executable Statements ..
 176 *
 177       RESULT( 1 ) = ZERO
 178       RESULT( 2 ) = ZERO
 179       IF( N.LE.0 .OR. M.LE.0 )
 180      $   RETURN
 181 *
 182       UNFL = SLAMCH( 'Safe minimum' )
 183       ULP = SLAMCH( 'Epsilon' )
 184 *
 185 *     Do Test 1
 186 *
 187 *     Compute the 1-norm of A.
 188 *
 189       IF( N.GT.1 ) THEN
 190          ANORM = ABS( AD( 1 ) ) + ABS( AE( 1 ) )
 191          DO 10 J = 2, N - 1
 192             ANORM = MAX( ANORM, ABS( AD( J ) )+ABS( AE( J ) )+
 193      $              ABS( AE( J-1 ) ) )
 194    10    CONTINUE
 195          ANORM = MAX( ANORM, ABS( AD( N ) )+ABS( AE( N-1 ) ) )
 196       ELSE
 197          ANORM = ABS( AD( 1 ) )
 198       END IF
 199       ANORM = MAX( ANORM, UNFL )
 200 *
 201 *     Norm of U'AU - S
 202 *
 203       DO 40 I = 1, M
 204          DO 30 J = 1, M
 205             WORK( I, J ) = ZERO
 206             DO 20 K = 1, N
 207                AUKJ = AD( K )*U( K, J )
 208                IF( K.NE.N )
 209      $            AUKJ = AUKJ + AE( K )*U( K+1, J )
 210                IF( K.NE.1 )
 211      $            AUKJ = AUKJ + AE( K-1 )*U( K-1, J )
 212                WORK( I, J ) = WORK( I, J ) + U( K, I )*AUKJ
 213    20       CONTINUE
 214    30    CONTINUE
 215          WORK( I, I ) = WORK( I, I ) - SD( I )
 216          IF( KBAND.EQ.1 ) THEN
 217             IF( I.NE.1 )
 218      $         WORK( I, I-1 ) = WORK( I, I-1 ) - SE( I-1 )
 219             IF( I.NE.N )
 220      $         WORK( I, I+1 ) = WORK( I, I+1 ) - SE( I )
 221          END IF
 222    40 CONTINUE
 223 *
 224       WNORM = SLANSY( '1', 'L', M, WORK, M, WORK( 1, M+1 ) )
 225 *
 226       IF( ANORM.GT.WNORM ) THEN
 227          RESULT( 1 ) = ( WNORM / ANORM ) / ( M*ULP )
 228       ELSE
 229          IF( ANORM.LT.ONE ) THEN
 230             RESULT( 1 ) = ( MIN( WNORM, M*ANORM ) / ANORM ) / ( M*ULP )
 231          ELSE
 232             RESULT( 1 ) = MIN( WNORM / ANORM, REAL( M ) ) / ( M*ULP )
 233          END IF
 234       END IF
 235 *
 236 *     Do Test 2
 237 *
 238 *     Compute  U'U - I
 239 *
 240       CALL SGEMM( 'T', 'N', M, M, N, ONE, U, LDU, U, LDU, ZERO, WORK,
 241      $            M )
 242 *
 243       DO 50 J = 1, M
 244          WORK( J, J ) = WORK( J, J ) - ONE
 245    50 CONTINUE
 246 *
 247       RESULT( 2 ) = MIN( REAL( M ), SLANGE( '1', M, M, WORK, M, WORK( 1,
 248      $              M+1 ) ) ) / ( M*ULP )
 249 *
 250       RETURN
 251 *
 252 *     End of SSTT22
 253 *
 254       END