LAPACK: SRC/dlag2.f Source File

Go to the documentation of this file.
00001       SUBROUTINE DLAG2( A, LDA, B, LDB, SAFMIN, SCALE1, SCALE2, WR1,
00002      $                  WR2, WI )
00003 *
00004 *  -- LAPACK auxiliary routine (version 3.2) --
00005 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
00006 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
00007 *     November 2006
00008 *
00009 *     .. Scalar Arguments ..
00010       INTEGER            LDA, LDB
00011       DOUBLE PRECISION   SAFMIN, SCALE1, SCALE2, WI, WR1, WR2
00012 *     ..
00013 *     .. Array Arguments ..
00014       DOUBLE PRECISION   A( LDA, * ), B( LDB, * )
00015 *     ..
00016 *
00017 *  Purpose
00018 *  =======
00019 *
00020 *  DLAG2 computes the eigenvalues of a 2 x 2 generalized eigenvalue
00021 *  problem  A - w B, with scaling as necessary to avoid over-/underflow.
00022 *
00023 *  The scaling factor "s" results in a modified eigenvalue equation
00024 *
00025 *      s A - w B
00026 *
00027 *  where  s  is a non-negative scaling factor chosen so that  w,  w B,
00028 *  and  s A  do not overflow and, if possible, do not underflow, either.
00029 *
00030 *  Arguments
00031 *  =========
00032 *
00033 *  A       (input) DOUBLE PRECISION array, dimension (LDA, 2)
00034 *          On entry, the 2 x 2 matrix A.  It is assumed that its 1-norm
00035 *          is less than 1/SAFMIN.  Entries less than
00036 *          sqrt(SAFMIN)*norm(A) are subject to being treated as zero.
00037 *
00038 *  LDA     (input) INTEGER
00039 *          The leading dimension of the array A.  LDA >= 2.
00040 *
00041 *  B       (input) DOUBLE PRECISION array, dimension (LDB, 2)
00042 *          On entry, the 2 x 2 upper triangular matrix B.  It is
00043 *          assumed that the one-norm of B is less than 1/SAFMIN.  The
00044 *          diagonals should be at least sqrt(SAFMIN) times the largest
00045 *          element of B (in absolute value); if a diagonal is smaller
00046 *          than that, then  +/- sqrt(SAFMIN) will be used instead of
00047 *          that diagonal.
00048 *
00049 *  LDB     (input) INTEGER
00050 *          The leading dimension of the array B.  LDB >= 2.
00051 *
00052 *  SAFMIN  (input) DOUBLE PRECISION
00053 *          The smallest positive number s.t. 1/SAFMIN does not
00054 *          overflow.  (This should always be DLAMCH('S') -- it is an
00055 *          argument in order to avoid having to call DLAMCH frequently.)
00056 *
00057 *  SCALE1  (output) DOUBLE PRECISION
00058 *          A scaling factor used to avoid over-/underflow in the
00059 *          eigenvalue equation which defines the first eigenvalue.  If
00060 *          the eigenvalues are complex, then the eigenvalues are
00061 *          ( WR1  +/-  WI i ) / SCALE1  (which may lie outside the
00062 *          exponent range of the machine), SCALE1=SCALE2, and SCALE1
00063 *          will always be positive.  If the eigenvalues are real, then
00064 *          the first (real) eigenvalue is  WR1 / SCALE1 , but this may
00065 *          overflow or underflow, and in fact, SCALE1 may be zero or
00066 *          less than the underflow threshhold if the exact eigenvalue
00067 *          is sufficiently large.
00068 *
00069 *  SCALE2  (output) DOUBLE PRECISION
00070 *          A scaling factor used to avoid over-/underflow in the
00071 *          eigenvalue equation which defines the second eigenvalue.  If
00072 *          the eigenvalues are complex, then SCALE2=SCALE1.  If the
00073 *          eigenvalues are real, then the second (real) eigenvalue is
00074 *          WR2 / SCALE2 , but this may overflow or underflow, and in
00075 *          fact, SCALE2 may be zero or less than the underflow
00076 *          threshhold if the exact eigenvalue is sufficiently large.
00077 *
00078 *  WR1     (output) DOUBLE PRECISION
00079 *          If the eigenvalue is real, then WR1 is SCALE1 times the
00080 *          eigenvalue closest to the (2,2) element of A B**(-1).  If the
00081 *          eigenvalue is complex, then WR1=WR2 is SCALE1 times the real
00082 *          part of the eigenvalues.
00083 *
00084 *  WR2     (output) DOUBLE PRECISION
00085 *          If the eigenvalue is real, then WR2 is SCALE2 times the
00086 *          other eigenvalue.  If the eigenvalue is complex, then
00087 *          WR1=WR2 is SCALE1 times the real part of the eigenvalues.
00088 *
00089 *  WI      (output) DOUBLE PRECISION
00090 *          If the eigenvalue is real, then WI is zero.  If the
00091 *          eigenvalue is complex, then WI is SCALE1 times the imaginary
00092 *          part of the eigenvalues.  WI will always be non-negative.
00093 *
00094 *  =====================================================================
00095 *
00096 *     .. Parameters ..
00097       DOUBLE PRECISION   ZERO, ONE, TWO
00098       PARAMETER          ( ZERO = 0.0D+0, ONE = 1.0D+0, TWO = 2.0D+0 )
00099       DOUBLE PRECISION   HALF
00100       PARAMETER          ( HALF = ONE / TWO )
00101       DOUBLE PRECISION   FUZZY1
00102       PARAMETER          ( FUZZY1 = ONE+1.0D-5 )
00103 *     ..
00104 *     .. Local Scalars ..
00105       DOUBLE PRECISION   A11, A12, A21, A22, ABI22, ANORM, AS11, AS12,
00106      $                   AS22, ASCALE, B11, B12, B22, BINV11, BINV22,
00107      $                   BMIN, BNORM, BSCALE, BSIZE, C1, C2, C3, C4, C5,
00108      $                   DIFF, DISCR, PP, QQ, R, RTMAX, RTMIN, S1, S2,
00109      $                   SAFMAX, SHIFT, SS, SUM, WABS, WBIG, WDET,
00110      $                   WSCALE, WSIZE, WSMALL
00111 *     ..
00112 *     .. Intrinsic Functions ..
00113       INTRINSIC          ABS, MAX, MIN, SIGN, SQRT
00114 *     ..
00115 *     .. Executable Statements ..
00116 *
00117       RTMIN = SQRT( SAFMIN )
00118       RTMAX = ONE / RTMIN
00119       SAFMAX = ONE / SAFMIN
00120 *
00121 *     Scale A
00122 *
00123       ANORM = MAX( ABS( A( 1, 1 ) )+ABS( A( 2, 1 ) ),
00124      $        ABS( A( 1, 2 ) )+ABS( A( 2, 2 ) ), SAFMIN )
00125       ASCALE = ONE / ANORM
00126       A11 = ASCALE*A( 1, 1 )
00127       A21 = ASCALE*A( 2, 1 )
00128       A12 = ASCALE*A( 1, 2 )
00129       A22 = ASCALE*A( 2, 2 )
00130 *
00131 *     Perturb B if necessary to insure non-singularity
00132 *
00133       B11 = B( 1, 1 )
00134       B12 = B( 1, 2 )
00135       B22 = B( 2, 2 )
00136       BMIN = RTMIN*MAX( ABS( B11 ), ABS( B12 ), ABS( B22 ), RTMIN )
00137       IF( ABS( B11 ).LT.BMIN )
00138      $   B11 = SIGN( BMIN, B11 )
00139       IF( ABS( B22 ).LT.BMIN )
00140      $   B22 = SIGN( BMIN, B22 )
00141 *
00142 *     Scale B
00143 *
00144       BNORM = MAX( ABS( B11 ), ABS( B12 )+ABS( B22 ), SAFMIN )
00145       BSIZE = MAX( ABS( B11 ), ABS( B22 ) )
00146       BSCALE = ONE / BSIZE
00147       B11 = B11*BSCALE
00148       B12 = B12*BSCALE
00149       B22 = B22*BSCALE
00150 *
00151 *     Compute larger eigenvalue by method described by C. van Loan
00152 *
00153 *     ( AS is A shifted by -SHIFT*B )
00154 *
00155       BINV11 = ONE / B11
00156       BINV22 = ONE / B22
00157       S1 = A11*BINV11
00158       S2 = A22*BINV22
00159       IF( ABS( S1 ).LE.ABS( S2 ) ) THEN
00160          AS12 = A12 - S1*B12
00161          AS22 = A22 - S1*B22
00162          SS = A21*( BINV11*BINV22 )
00163          ABI22 = AS22*BINV22 - SS*B12
00164          PP = HALF*ABI22
00165          SHIFT = S1
00166       ELSE
00167          AS12 = A12 - S2*B12
00168          AS11 = A11 - S2*B11
00169          SS = A21*( BINV11*BINV22 )
00170          ABI22 = -SS*B12
00171          PP = HALF*( AS11*BINV11+ABI22 )
00172          SHIFT = S2
00173       END IF
00174       QQ = SS*AS12
00175       IF( ABS( PP*RTMIN ).GE.ONE ) THEN
00176          DISCR = ( RTMIN*PP )**2 + QQ*SAFMIN
00177          R = SQRT( ABS( DISCR ) )*RTMAX
00178       ELSE
00179          IF( PP**2+ABS( QQ ).LE.SAFMIN ) THEN
00180             DISCR = ( RTMAX*PP )**2 + QQ*SAFMAX
00181             R = SQRT( ABS( DISCR ) )*RTMIN
00182          ELSE
00183             DISCR = PP**2 + QQ
00184             R = SQRT( ABS( DISCR ) )
00185          END IF
00186       END IF
00187 *
00188 *     Note: the test of R in the following IF is to cover the case when
00189 *           DISCR is small and negative and is flushed to zero during
00190 *           the calculation of R.  On machines which have a consistent
00191 *           flush-to-zero threshhold and handle numbers above that
00192 *           threshhold correctly, it would not be necessary.
00193 *
00194       IF( DISCR.GE.ZERO .OR. R.EQ.ZERO ) THEN
00195          SUM = PP + SIGN( R, PP )
00196          DIFF = PP - SIGN( R, PP )
00197          WBIG = SHIFT + SUM
00198 *
00199 *        Compute smaller eigenvalue
00200 *
00201          WSMALL = SHIFT + DIFF
00202          IF( HALF*ABS( WBIG ).GT.MAX( ABS( WSMALL ), SAFMIN ) ) THEN
00203             WDET = ( A11*A22-A12*A21 )*( BINV11*BINV22 )
00204             WSMALL = WDET / WBIG
00205          END IF
00206 *
00207 *        Choose (real) eigenvalue closest to 2,2 element of A*B**(-1)
00208 *        for WR1.
00209 *
00210          IF( PP.GT.ABI22 ) THEN
00211             WR1 = MIN( WBIG, WSMALL )
00212             WR2 = MAX( WBIG, WSMALL )
00213          ELSE
00214             WR1 = MAX( WBIG, WSMALL )
00215             WR2 = MIN( WBIG, WSMALL )
00216          END IF
00217          WI = ZERO
00218       ELSE
00219 *
00220 *        Complex eigenvalues
00221 *
00222          WR1 = SHIFT + PP
00223          WR2 = WR1
00224          WI = R
00225       END IF
00226 *
00227 *     Further scaling to avoid underflow and overflow in computing
00228 *     SCALE1 and overflow in computing w*B.
00229 *
00230 *     This scale factor (WSCALE) is bounded from above using C1 and C2,
00231 *     and from below using C3 and C4.
00232 *        C1 implements the condition  s A  must never overflow.
00233 *        C2 implements the condition  w B  must never overflow.
00234 *        C3, with C2,
00235 *           implement the condition that s A - w B must never overflow.
00236 *        C4 implements the condition  s    should not underflow.
00237 *        C5 implements the condition  max(s,|w|) should be at least 2.
00238 *
00239       C1 = BSIZE*( SAFMIN*MAX( ONE, ASCALE ) )
00240       C2 = SAFMIN*MAX( ONE, BNORM )
00241       C3 = BSIZE*SAFMIN
00242       IF( ASCALE.LE.ONE .AND. BSIZE.LE.ONE ) THEN
00243          C4 = MIN( ONE, ( ASCALE / SAFMIN )*BSIZE )
00244       ELSE
00245          C4 = ONE
00246       END IF
00247       IF( ASCALE.LE.ONE .OR. BSIZE.LE.ONE ) THEN
00248          C5 = MIN( ONE, ASCALE*BSIZE )
00249       ELSE
00250          C5 = ONE
00251       END IF
00252 *
00253 *     Scale first eigenvalue
00254 *
00255       WABS = ABS( WR1 ) + ABS( WI )
00256       WSIZE = MAX( SAFMIN, C1, FUZZY1*( WABS*C2+C3 ),
00257      $        MIN( C4, HALF*MAX( WABS, C5 ) ) )
00258       IF( WSIZE.NE.ONE ) THEN
00259          WSCALE = ONE / WSIZE
00260          IF( WSIZE.GT.ONE ) THEN
00261             SCALE1 = ( MAX( ASCALE, BSIZE )*WSCALE )*
00262      $               MIN( ASCALE, BSIZE )
00263          ELSE
00264             SCALE1 = ( MIN( ASCALE, BSIZE )*WSCALE )*
00265      $               MAX( ASCALE, BSIZE )
00266          END IF
00267          WR1 = WR1*WSCALE
00268          IF( WI.NE.ZERO ) THEN
00269             WI = WI*WSCALE
00270             WR2 = WR1
00271             SCALE2 = SCALE1
00272          END IF
00273       ELSE
00274          SCALE1 = ASCALE*BSIZE
00275          SCALE2 = SCALE1
00276       END IF
00277 *
00278 *     Scale second eigenvalue (if real)
00279 *
00280       IF( WI.EQ.ZERO ) THEN
00281          WSIZE = MAX( SAFMIN, C1, FUZZY1*( ABS( WR2 )*C2+C3 ),
00282      $           MIN( C4, HALF*MAX( ABS( WR2 ), C5 ) ) )
00283          IF( WSIZE.NE.ONE ) THEN
00284             WSCALE = ONE / WSIZE
00285             IF( WSIZE.GT.ONE ) THEN
00286                SCALE2 = ( MAX( ASCALE, BSIZE )*WSCALE )*
00287      $                  MIN( ASCALE, BSIZE )
00288             ELSE
00289                SCALE2 = ( MIN( ASCALE, BSIZE )*WSCALE )*
00290      $                  MAX( ASCALE, BSIZE )
00291             END IF
00292             WR2 = WR2*WSCALE
00293          ELSE
00294             SCALE2 = ASCALE*BSIZE
00295          END IF
00296       END IF
00297 *
00298 *     End of DLAG2
00299 *
00300       RETURN
00301       END