numeric-linalg

Educational material on the SciPy implementation of numerical linear algebra algorithms

Name	Size	Mode
..
lapack/TESTING/EIG/dget53.f	7089B	-rw-r--r--

*> \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