TESTING/EIG/sget53.f

   1 *> \brief \b SGET53
   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 SGET53( A, LDA, B, LDB, SCALE, WR, WI, RESULT, INFO )
  12 *
  13 *       .. Scalar Arguments ..
  14 *       INTEGER            INFO, LDA, LDB
  15 *       REAL               RESULT, SCALE, WI, WR
  16 *       ..
  17 *       .. Array Arguments ..
  18 *       REAL               A( LDA, * ), B( LDB, * )
  19 *       ..
  20 *
  21 *
  22 *> \par Purpose:
  23 *  =============
  24 *>
  25 *> \verbatim
  26 *>
  27 *> SGET53  checks the generalized eigenvalues computed by SLAG2.
  28 *>
  29 *> The basic test for an eigenvalue is:
  30 *>
  31 *>                              | det( s A - w B ) |
  32 *>     RESULT =  ---------------------------------------------------
  33 *>               ulp max( s norm(A), |w| norm(B) )*norm( s A - w B )
  34 *>
  35 *> Two "safety checks" are performed:
  36 *>
  37 *> (1)  ulp*max( s*norm(A), |w|*norm(B) )  must be at least
  38 *>      safe_minimum.  This insures that the test performed is
  39 *>      not essentially  det(0*A + 0*B)=0.
  40 *>
  41 *> (2)  s*norm(A) + |w|*norm(B) must be less than 1/safe_minimum.
  42 *>      This insures that  s*A - w*B  will not overflow.
  43 *>
  44 *> If these tests are not passed, then  s  and  w  are scaled and
  45 *> tested anyway, if this is possible.
  46 *> \endverbatim
  47 *
  48 *  Arguments:
  49 *  ==========
  50 *
  51 *> \param[in] A
  52 *> \verbatim
  53 *>          A is REAL array, dimension (LDA, 2)
  54 *>          The 2x2 matrix A.
  55 *> \endverbatim
  56 *>
  57 *> \param[in] LDA
  58 *> \verbatim
  59 *>          LDA is INTEGER
  60 *>          The leading dimension of A.  It must be at least 2.
  61 *> \endverbatim
  62 *>
  63 *> \param[in] B
  64 *> \verbatim
  65 *>          B is REAL array, dimension (LDB, N)
  66 *>          The 2x2 upper-triangular matrix B.
  67 *> \endverbatim
  68 *>
  69 *> \param[in] LDB
  70 *> \verbatim
  71 *>          LDB is INTEGER
  72 *>          The leading dimension of B.  It must be at least 2.
  73 *> \endverbatim
  74 *>
  75 *> \param[in] SCALE
  76 *> \verbatim
  77 *>          SCALE is REAL
  78 *>          The "scale factor" s in the formula  s A - w B .  It is
  79 *>          assumed to be non-negative.
  80 *> \endverbatim
  81 *>
  82 *> \param[in] WR
  83 *> \verbatim
  84 *>          WR is REAL
  85 *>          The real part of the eigenvalue  w  in the formula
  86 *>          s A - w B .
  87 *> \endverbatim
  88 *>
  89 *> \param[in] WI
  90 *> \verbatim
  91 *>          WI is REAL
  92 *>          The imaginary part of the eigenvalue  w  in the formula
  93 *>          s A - w B .
  94 *> \endverbatim
  95 *>
  96 *> \param[out] RESULT
  97 *> \verbatim
  98 *>          RESULT is REAL
  99 *>          If INFO is 2 or less, the value computed by the test
 100 *>             described above.
 101 *>          If INFO=3, this will just be 1/ulp.
 102 *> \endverbatim
 103 *>
 104 *> \param[out] INFO
 105 *> \verbatim
 106 *>          INFO is INTEGER
 107 *>          =0:  The input data pass the "safety checks".
 108 *>          =1:  s*norm(A) + |w|*norm(B) > 1/safe_minimum.
 109 *>          =2:  ulp*max( s*norm(A), |w|*norm(B) ) < safe_minimum
 110 *>          =3:  same as INFO=2, but  s  and  w  could not be scaled so
 111 *>               as to compute the test.
 112 *> \endverbatim
 113 *
 114 *  Authors:
 115 *  ========
 116 *
 117 *> \author Univ. of Tennessee
 118 *> \author Univ. of California Berkeley
 119 *> \author Univ. of Colorado Denver
 120 *> \author NAG Ltd.
 121 *
 122 *> \date November 2011
 123 *
 124 *> \ingroup single_eig
 125 *
 126 *  =====================================================================
 127       SUBROUTINE SGET53( A, LDA, B, LDB, SCALE, WR, WI, RESULT, INFO )
 128 *
 129 *  -- LAPACK test routine (version 3.4.0) --
 130 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
 131 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
 132 *     November 2011
 133 *
 134 *     .. Scalar Arguments ..
 135       INTEGER            INFO, LDA, LDB
 136       REAL               RESULT, SCALE, WI, WR
 137 *     ..
 138 *     .. Array Arguments ..
 139       REAL               A( LDA, * ), B( LDB, * )
 140 *     ..
 141 *
 142 *  =====================================================================
 143 *
 144 *     .. Parameters ..
 145       REAL               ZERO, ONE
 146       PARAMETER          ( ZERO = 0.0, ONE = 1.0 )
 147 *     ..
 148 *     .. Local Scalars ..
 149       REAL               ABSW, ANORM, BNORM, CI11, CI12, CI22, CNORM,
 150      $                   CR11, CR12, CR21, CR22, CSCALE, DETI, DETR, S1,
 151      $                   SAFMIN, SCALES, SIGMIN, TEMP, ULP, WIS, WRS
 152 *     ..
 153 *     .. External Functions ..
 154       REAL               SLAMCH
 155       EXTERNAL           SLAMCH
 156 *     ..
 157 *     .. Intrinsic Functions ..
 158       INTRINSIC          ABS, MAX, SQRT
 159 *     ..
 160 *     .. Executable Statements ..
 161 *
 162 *     Initialize
 163 *
 164       INFO = 0
 165       RESULT = ZERO
 166       SCALES = SCALE
 167       WRS = WR
 168       WIS = WI
 169 *
 170 *     Machine constants and norms
 171 *
 172       SAFMIN = SLAMCH( 'Safe minimum' )
 173       ULP = SLAMCH( 'Epsilon' )*SLAMCH( 'Base' )
 174       ABSW = ABS( WRS ) + ABS( WIS )
 175       ANORM = MAX( ABS( A( 1, 1 ) )+ABS( A( 2, 1 ) ),
 176      $        ABS( A( 1, 2 ) )+ABS( A( 2, 2 ) ), SAFMIN )
 177       BNORM = MAX( ABS( B( 1, 1 ) ), ABS( B( 1, 2 ) )+ABS( B( 2, 2 ) ),
 178      $        SAFMIN )
 179 *
 180 *     Check for possible overflow.
 181 *
 182       TEMP = ( SAFMIN*BNORM )*ABSW + ( SAFMIN*ANORM )*SCALES
 183       IF( TEMP.GE.ONE ) THEN
 184 *
 185 *        Scale down to avoid overflow
 186 *
 187          INFO = 1
 188          TEMP = ONE / TEMP
 189          SCALES = SCALES*TEMP
 190          WRS = WRS*TEMP
 191          WIS = WIS*TEMP
 192          ABSW = ABS( WRS ) + ABS( WIS )
 193       END IF
 194       S1 = MAX( ULP*MAX( SCALES*ANORM, ABSW*BNORM ),
 195      $     SAFMIN*MAX( SCALES, ABSW ) )
 196 *
 197 *     Check for W and SCALE essentially zero.
 198 *
 199       IF( S1.LT.SAFMIN ) THEN
 200          INFO = 2
 201          IF( SCALES.LT.SAFMIN .AND. ABSW.LT.SAFMIN ) THEN
 202             INFO = 3
 203             RESULT = ONE / ULP
 204             RETURN
 205          END IF
 206 *
 207 *        Scale up to avoid underflow
 208 *
 209          TEMP = ONE / MAX( SCALES*ANORM+ABSW*BNORM, SAFMIN )
 210          SCALES = SCALES*TEMP
 211          WRS = WRS*TEMP
 212          WIS = WIS*TEMP
 213          ABSW = ABS( WRS ) + ABS( WIS )
 214          S1 = MAX( ULP*MAX( SCALES*ANORM, ABSW*BNORM ),
 215      $        SAFMIN*MAX( SCALES, ABSW ) )
 216          IF( S1.LT.SAFMIN ) THEN
 217             INFO = 3
 218             RESULT = ONE / ULP
 219             RETURN
 220          END IF
 221       END IF
 222 *
 223 *     Compute C = s A - w B
 224 *
 225       CR11 = SCALES*A( 1, 1 ) - WRS*B( 1, 1 )
 226       CI11 = -WIS*B( 1, 1 )
 227       CR21 = SCALES*A( 2, 1 )
 228       CR12 = SCALES*A( 1, 2 ) - WRS*B( 1, 2 )
 229       CI12 = -WIS*B( 1, 2 )
 230       CR22 = SCALES*A( 2, 2 ) - WRS*B( 2, 2 )
 231       CI22 = -WIS*B( 2, 2 )
 232 *
 233 *     Compute the smallest singular value of s A - w B:
 234 *
 235 *                 |det( s A - w B )|
 236 *     sigma_min = ------------------
 237 *                 norm( s A - w B )
 238 *
 239       CNORM = MAX( ABS( CR11 )+ABS( CI11 )+ABS( CR21 ),
 240      $        ABS( CR12 )+ABS( CI12 )+ABS( CR22 )+ABS( CI22 ), SAFMIN )
 241       CSCALE = ONE / SQRT( CNORM )
 242       DETR = ( CSCALE*CR11 )*( CSCALE*CR22 ) -
 243      $       ( CSCALE*CI11 )*( CSCALE*CI22 ) -
 244      $       ( CSCALE*CR12 )*( CSCALE*CR21 )
 245       DETI = ( CSCALE*CR11 )*( CSCALE*CI22 ) +
 246      $       ( CSCALE*CI11 )*( CSCALE*CR22 ) -
 247      $       ( CSCALE*CI12 )*( CSCALE*CR21 )
 248       SIGMIN = ABS( DETR ) + ABS( DETI )
 249       RESULT = SIGMIN / S1
 250       RETURN
 251 *
 252 *     End of SGET53
 253 *
 254       END