LAPACK: TESTING/EIG/dget53.f Source File

Go to the documentation of this file.
00001       SUBROUTINE DGET53( A, LDA, B, LDB, SCALE, WR, WI, RESULT, INFO )
00002 *
00003 *  -- LAPACK test routine (version 3.1) --
00004 *     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
00005 *     November 2006
00006 *
00007 *     .. Scalar Arguments ..
00008       INTEGER            INFO, LDA, LDB
00009       DOUBLE PRECISION   RESULT, SCALE, WI, WR
00010 *     ..
00011 *     .. Array Arguments ..
00012       DOUBLE PRECISION   A( LDA, * ), B( LDB, * )
00013 *     ..
00014 *
00015 *  Purpose
00016 *  =======
00017 *
00018 *  DGET53  checks the generalized eigenvalues computed by DLAG2.
00019 *
00020 *  The basic test for an eigenvalue is:
00021 *
00022 *                               | det( s A - w B ) |
00023 *      RESULT =  ---------------------------------------------------
00024 *                ulp max( s norm(A), |w| norm(B) )*norm( s A - w B )
00025 *
00026 *  Two "safety checks" are performed:
00027 *
00028 *  (1)  ulp*max( s*norm(A), |w|*norm(B) )  must be at least
00029 *       safe_minimum.  This insures that the test performed is
00030 *       not essentially  det(0*A + 0*B)=0.
00031 *
00032 *  (2)  s*norm(A) + |w|*norm(B) must be less than 1/safe_minimum.
00033 *       This insures that  s*A - w*B  will not overflow.
00034 *
00035 *  If these tests are not passed, then  s  and  w  are scaled and
00036 *  tested anyway, if this is possible.
00037 *
00038 *  Arguments
00039 *  =========
00040 *
00041 *  A       (input) DOUBLE PRECISION array, dimension (LDA, 2)
00042 *          The 2x2 matrix A.
00043 *
00044 *  LDA     (input) INTEGER
00045 *          The leading dimension of A.  It must be at least 2.
00046 *
00047 *  B       (input) DOUBLE PRECISION array, dimension (LDB, N)
00048 *          The 2x2 upper-triangular matrix B.
00049 *
00050 *  LDB     (input) INTEGER
00051 *          The leading dimension of B.  It must be at least 2.
00052 *
00053 *  SCALE   (input) DOUBLE PRECISION
00054 *          The "scale factor" s in the formula  s A - w B .  It is
00055 *          assumed to be non-negative.
00056 *
00057 *  WR      (input) DOUBLE PRECISION
00058 *          The real part of the eigenvalue  w  in the formula
00059 *          s A - w B .
00060 *
00061 *  WI      (input) DOUBLE PRECISION
00062 *          The imaginary part of the eigenvalue  w  in the formula
00063 *          s A - w B .
00064 *
00065 *  RESULT  (output) DOUBLE PRECISION
00066 *          If INFO is 2 or less, the value computed by the test
00067 *             described above.
00068 *          If INFO=3, this will just be 1/ulp.
00069 *
00070 *  INFO    (output) INTEGER
00071 *          =0:  The input data pass the "safety checks".
00072 *          =1:  s*norm(A) + |w|*norm(B) > 1/safe_minimum.
00073 *          =2:  ulp*max( s*norm(A), |w|*norm(B) ) < safe_minimum
00074 *          =3:  same as INFO=2, but  s  and  w  could not be scaled so
00075 *               as to compute the test.
00076 *
00077 *  =====================================================================
00078 *
00079 *     .. Parameters ..
00080       DOUBLE PRECISION   ZERO, ONE
00081       PARAMETER          ( ZERO = 0.0D0, ONE = 1.0D0 )
00082 *     ..
00083 *     .. Local Scalars ..
00084       DOUBLE PRECISION   ABSW, ANORM, BNORM, CI11, CI12, CI22, CNORM,
00085      $                   CR11, CR12, CR21, CR22, CSCALE, DETI, DETR, S1,
00086      $                   SAFMIN, SCALES, SIGMIN, TEMP, ULP, WIS, WRS
00087 *     ..
00088 *     .. External Functions ..
00089       DOUBLE PRECISION   DLAMCH
00090       EXTERNAL           DLAMCH
00091 *     ..
00092 *     .. Intrinsic Functions ..
00093       INTRINSIC          ABS, MAX, SQRT
00094 *     ..
00095 *     .. Executable Statements ..
00096 *
00097 *     Initialize
00098 *
00099       INFO = 0
00100       RESULT = ZERO
00101       SCALES = SCALE
00102       WRS = WR
00103       WIS = WI
00104 *
00105 *     Machine constants and norms
00106 *
00107       SAFMIN = DLAMCH( 'Safe minimum' )
00108       ULP = DLAMCH( 'Epsilon' )*DLAMCH( 'Base' )
00109       ABSW = ABS( WRS ) + ABS( WIS )
00110       ANORM = MAX( ABS( A( 1, 1 ) )+ABS( A( 2, 1 ) ),
00111      $        ABS( A( 1, 2 ) )+ABS( A( 2, 2 ) ), SAFMIN )
00112       BNORM = MAX( ABS( B( 1, 1 ) ), ABS( B( 1, 2 ) )+ABS( B( 2, 2 ) ),
00113      $        SAFMIN )
00114 *
00115 *     Check for possible overflow.
00116 *
00117       TEMP = ( SAFMIN*BNORM )*ABSW + ( SAFMIN*ANORM )*SCALES
00118       IF( TEMP.GE.ONE ) THEN
00119 *
00120 *        Scale down to avoid overflow
00121 *
00122          INFO = 1
00123          TEMP = ONE / TEMP
00124          SCALES = SCALES*TEMP
00125          WRS = WRS*TEMP
00126          WIS = WIS*TEMP
00127          ABSW = ABS( WRS ) + ABS( WIS )
00128       END IF
00129       S1 = MAX( ULP*MAX( SCALES*ANORM, ABSW*BNORM ),
00130      $     SAFMIN*MAX( SCALES, ABSW ) )
00131 *
00132 *     Check for W and SCALE essentially zero.
00133 *
00134       IF( S1.LT.SAFMIN ) THEN
00135          INFO = 2
00136          IF( SCALES.LT.SAFMIN .AND. ABSW.LT.SAFMIN ) THEN
00137             INFO = 3
00138             RESULT = ONE / ULP
00139             RETURN
00140          END IF
00141 *
00142 *        Scale up to avoid underflow
00143 *
00144          TEMP = ONE / MAX( SCALES*ANORM+ABSW*BNORM, SAFMIN )
00145          SCALES = SCALES*TEMP
00146          WRS = WRS*TEMP
00147          WIS = WIS*TEMP
00148          ABSW = ABS( WRS ) + ABS( WIS )
00149          S1 = MAX( ULP*MAX( SCALES*ANORM, ABSW*BNORM ),
00150      $        SAFMIN*MAX( SCALES, ABSW ) )
00151          IF( S1.LT.SAFMIN ) THEN
00152             INFO = 3
00153             RESULT = ONE / ULP
00154             RETURN
00155          END IF
00156       END IF
00157 *
00158 *     Compute C = s A - w B
00159 *
00160       CR11 = SCALES*A( 1, 1 ) - WRS*B( 1, 1 )
00161       CI11 = -WIS*B( 1, 1 )
00162       CR21 = SCALES*A( 2, 1 )
00163       CR12 = SCALES*A( 1, 2 ) - WRS*B( 1, 2 )
00164       CI12 = -WIS*B( 1, 2 )
00165       CR22 = SCALES*A( 2, 2 ) - WRS*B( 2, 2 )
00166       CI22 = -WIS*B( 2, 2 )
00167 *
00168 *     Compute the smallest singular value of s A - w B:
00169 *
00170 *                 |det( s A - w B )|
00171 *     sigma_min = ------------------
00172 *                 norm( s A - w B )
00173 *
00174       CNORM = MAX( ABS( CR11 )+ABS( CI11 )+ABS( CR21 ),
00175      $        ABS( CR12 )+ABS( CI12 )+ABS( CR22 )+ABS( CI22 ), SAFMIN )
00176       CSCALE = ONE / SQRT( CNORM )
00177       DETR = ( CSCALE*CR11 )*( CSCALE*CR22 ) -
00178      $       ( CSCALE*CI11 )*( CSCALE*CI22 ) -
00179      $       ( CSCALE*CR12 )*( CSCALE*CR21 )
00180       DETI = ( CSCALE*CR11 )*( CSCALE*CI22 ) +
00181      $       ( CSCALE*CI11 )*( CSCALE*CR22 ) -
00182      $       ( CSCALE*CI12 )*( CSCALE*CR21 )
00183       SIGMIN = ABS( DETR ) + ABS( DETI )
00184       RESULT = SIGMIN / S1
00185       RETURN
00186 *
00187 *     End of DGET53
00188 *
00189       END