dget53.f

Go to the documentation of this file.

*> \brief \b DGET53
*
*  =========== DOCUMENTATION ===========
*
* Online html documentation available at
*            http://www.netlib.org/lapack/explore-html/
*
*  Definition:
*  ===========
*
*       SUBROUTINE DGET53( A, LDA, B, LDB, SCALE, WR, WI, RESULT, INFO )
*
*       .. Scalar Arguments ..
*       INTEGER            INFO, LDA, LDB
*       DOUBLE PRECISION   RESULT, SCALE, WI, WR
*       ..
*       .. Array Arguments ..
*       DOUBLE PRECISION   A( LDA, * ), B( LDB, * )
*       ..
*
*
*> \par Purpose:
*  =============
*>
*> \verbatim
*>
*> DGET53  checks the generalized eigenvalues computed by DLAG2.
*>
*> The basic test for an eigenvalue is:
*>
*>                              | det( s A - w B ) |
*>     RESULT =  ---------------------------------------------------
*>               ulp max( s norm(A), |w| norm(B) )*norm( s A - w B )
*>
*> Two "safety checks" are performed:
*>
*> (1)  ulp*max( s*norm(A), |w|*norm(B) )  must be at least
*>      safe_minimum.  This insures that the test performed is
*>      not essentially  det(0*A + 0*B)=0.
*>
*> (2)  s*norm(A) + |w|*norm(B) must be less than 1/safe_minimum.
*>      This insures that  s*A - w*B  will not overflow.
*>
*> If these tests are not passed, then  s  and  w  are scaled and
*> tested anyway, if this is possible.
*> \endverbatim
*
*  Arguments:
*  ==========
*
*> \param[in] A
*> \verbatim
*>          A is DOUBLE PRECISION array, dimension (LDA, 2)
*>          The 2x2 matrix A.
*> \endverbatim
*>
*> \param[in] LDA
*> \verbatim
*>          LDA is INTEGER
*>          The leading dimension of A.  It must be at least 2.
*> \endverbatim
*>
*> \param[in] B
*> \verbatim
*>          B is DOUBLE PRECISION array, dimension (LDB, N)
*>          The 2x2 upper-triangular matrix B.
*> \endverbatim
*>
*> \param[in] LDB
*> \verbatim
*>          LDB is INTEGER
*>          The leading dimension of B.  It must be at least 2.
*> \endverbatim
*>
*> \param[in] SCALE
*> \verbatim
*>          SCALE is DOUBLE PRECISION
*>          The "scale factor" s in the formula  s A - w B .  It is
*>          assumed to be non-negative.
*> \endverbatim
*>
*> \param[in] WR
*> \verbatim
*>          WR is DOUBLE PRECISION
*>          The real part of the eigenvalue  w  in the formula
*>          s A - w B .
*> \endverbatim
*>
*> \param[in] WI
*> \verbatim
*>          WI is DOUBLE PRECISION
*>          The imaginary part of the eigenvalue  w  in the formula
*>          s A - w B .
*> \endverbatim
*>
*> \param[out] RESULT
*> \verbatim
*>          RESULT is DOUBLE PRECISION
*>          If INFO is 2 or less, the value computed by the test
*>             described above.
*>          If INFO=3, this will just be 1/ulp.
*> \endverbatim
*>
*> \param[out] INFO
*> \verbatim
*>          INFO is INTEGER
*>          =0:  The input data pass the "safety checks".
*>          =1:  s*norm(A) + |w|*norm(B) > 1/safe_minimum.
*>          =2:  ulp*max( s*norm(A), |w|*norm(B) ) < safe_minimum
*>          =3:  same as INFO=2, but  s  and  w  could not be scaled so
*>               as to compute the test.
*> \endverbatim
*
*  Authors:
*  ========
*
*> \author Univ. of Tennessee
*> \author Univ. of California Berkeley
*> \author Univ. of Colorado Denver
*> \author NAG Ltd.
*
*> \ingroup double_eig
*
*  =====================================================================
      SUBROUTINE dget53( A, LDA, B, LDB, SCALE, WR, WI, RESULT, INFO )
*
*  -- LAPACK test routine --
*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
*
*     .. Scalar Arguments ..
      INTEGER            INFO, LDA, LDB
      DOUBLE PRECISION   RESULT, SCALE, WI, WR
*     ..
*     .. Array Arguments ..
      DOUBLE PRECISION   A( LDA, * ), B( LDB, * )
*     ..
*
*  =====================================================================
*
*     .. Parameters ..
      DOUBLE PRECISION   ZERO, ONE
      parameter( zero = 0.0d0, one = 1.0d0 )
*     ..
*     .. Local Scalars ..
      DOUBLE PRECISION   ABSW, ANORM, BNORM, CI11, CI12, CI22, CNORM,
     $                   CR11, CR12, CR21, CR22, CSCALE, DETI, DETR, S1,
     $                   SAFMIN, SCALES, SIGMIN, TEMP, ULP, WIS, WRS
*     ..
*     .. External Functions ..
      DOUBLE PRECISION   DLAMCH
      EXTERNAL           dlamch
*     ..
*     .. Intrinsic Functions ..
      INTRINSIC          abs, max, sqrt
*     ..
*     .. Executable Statements ..
*
*     Initialize
*
      info = 0
      result = zero
      scales = scale
      wrs = wr
      wis = wi
*
*     Machine constants and norms
*
      safmin = dlamch( 'Safe minimum' )
      ulp = dlamch( 'Epsilon' )*dlamch( 'Base' )
      absw = abs( wrs ) + abs( wis )
      anorm = max( abs( a( 1, 1 ) )+abs( a( 2, 1 ) ),
     $        abs( a( 1, 2 ) )+abs( a( 2, 2 ) ), safmin )
      bnorm = max( abs( b( 1, 1 ) ), abs( b( 1, 2 ) )+abs( b( 2, 2 ) ),
     $        safmin )
*
*     Check for possible overflow.
*
      temp = ( safmin*bnorm )*absw + ( safmin*anorm )*scales
      IF( temp.GE.one ) THEN
*
*        Scale down to avoid overflow
*
         info = 1
         temp = one / temp
         scales = scales*temp
         wrs = wrs*temp
         wis = wis*temp
         absw = abs( wrs ) + abs( wis )
      END IF
      s1 = max( ulp*max( scales*anorm, absw*bnorm ),
     $     safmin*max( scales, absw ) )
*
*     Check for W and SCALE essentially zero.
*
      IF( s1.LT.safmin ) THEN
         info = 2
         IF( scales.LT.safmin .AND. absw.LT.safmin ) THEN
            info = 3
            result = one / ulp
            RETURN
         END IF
*
*        Scale up to avoid underflow
*
         temp = one / max( scales*anorm+absw*bnorm, safmin )
         scales = scales*temp
         wrs = wrs*temp
         wis = wis*temp
         absw = abs( wrs ) + abs( wis )
         s1 = max( ulp*max( scales*anorm, absw*bnorm ),
     $        safmin*max( scales, absw ) )
         IF( s1.LT.safmin ) THEN
            info = 3
            result = one / ulp
            RETURN
         END IF
      END IF
*
*     Compute C = s A - w B
*
      cr11 = scales*a( 1, 1 ) - wrs*b( 1, 1 )
      ci11 = -wis*b( 1, 1 )
      cr21 = scales*a( 2, 1 )
      cr12 = scales*a( 1, 2 ) - wrs*b( 1, 2 )
      ci12 = -wis*b( 1, 2 )
      cr22 = scales*a( 2, 2 ) - wrs*b( 2, 2 )
      ci22 = -wis*b( 2, 2 )
*
*     Compute the smallest singular value of s A - w B:
*
*                 |det( s A - w B )|
*     sigma_min = ------------------
*                 norm( s A - w B )
*
      cnorm = max( abs( cr11 )+abs( ci11 )+abs( cr21 ),
     $        abs( cr12 )+abs( ci12 )+abs( cr22 )+abs( ci22 ), safmin )
      cscale = one / sqrt( cnorm )
      detr = ( cscale*cr11 )*( cscale*cr22 ) -
     $       ( cscale*ci11 )*( cscale*ci22 ) -
     $       ( cscale*cr12 )*( cscale*cr21 )
      deti = ( cscale*cr11 )*( cscale*ci22 ) +
     $       ( cscale*ci11 )*( cscale*cr22 ) -
     $       ( cscale*ci12 )*( cscale*cr21 )
      sigmin = abs( detr ) + abs( deti )
      result = sigmin / s1
      RETURN
*
*     End of DGET53
*
      END