SRC/slag2.f

   1 *> \brief \b SLAG2 computes the eigenvalues of a 2-by-2 generalized eigenvalue problem, with scaling as necessary to avoid over-/underflow.
   2 *
   3 *  =========== DOCUMENTATION ===========
   4 *
   5 * Online html documentation available at
   6 *            http://www.netlib.org/lapack/explore-html/
   7 *
   8 *> \htmlonly
   9 *> Download SLAG2 + dependencies
  10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/slag2.f">
  11 *> [TGZ]</a>
  12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/slag2.f">
  13 *> [ZIP]</a>
  14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/slag2.f">
  15 *> [TXT]</a>
  16 *> \endhtmlonly
  17 *
  18 *  Definition:
  19 *  ===========
  20 *
  21 *       SUBROUTINE SLAG2( A, LDA, B, LDB, SAFMIN, SCALE1, SCALE2, WR1,
  22 *                         WR2, WI )
  23 *
  24 *       .. Scalar Arguments ..
  25 *       INTEGER            LDA, LDB
  26 *       REAL               SAFMIN, SCALE1, SCALE2, WI, WR1, WR2
  27 *       ..
  28 *       .. Array Arguments ..
  29 *       REAL               A( LDA, * ), B( LDB, * )
  30 *       ..
  31 *
  32 *
  33 *> \par Purpose:
  34 *  =============
  35 *>
  36 *> \verbatim
  37 *>
  38 *> SLAG2 computes the eigenvalues of a 2 x 2 generalized eigenvalue
  39 *> problem  A - w B, with scaling as necessary to avoid over-/underflow.
  40 *>
  41 *> The scaling factor "s" results in a modified eigenvalue equation
  42 *>
  43 *>     s A - w B
  44 *>
  45 *> where  s  is a non-negative scaling factor chosen so that  w,  w B,
  46 *> and  s A  do not overflow and, if possible, do not underflow, either.
  47 *> \endverbatim
  48 *
  49 *  Arguments:
  50 *  ==========
  51 *
  52 *> \param[in] A
  53 *> \verbatim
  54 *>          A is REAL array, dimension (LDA, 2)
  55 *>          On entry, the 2 x 2 matrix A.  It is assumed that its 1-norm
  56 *>          is less than 1/SAFMIN.  Entries less than
  57 *>          sqrt(SAFMIN)*norm(A) are subject to being treated as zero.
  58 *> \endverbatim
  59 *>
  60 *> \param[in] LDA
  61 *> \verbatim
  62 *>          LDA is INTEGER
  63 *>          The leading dimension of the array A.  LDA >= 2.
  64 *> \endverbatim
  65 *>
  66 *> \param[in] B
  67 *> \verbatim
  68 *>          B is REAL array, dimension (LDB, 2)
  69 *>          On entry, the 2 x 2 upper triangular matrix B.  It is
  70 *>          assumed that the one-norm of B is less than 1/SAFMIN.  The
  71 *>          diagonals should be at least sqrt(SAFMIN) times the largest
  72 *>          element of B (in absolute value); if a diagonal is smaller
  73 *>          than that, then  +/- sqrt(SAFMIN) will be used instead of
  74 *>          that diagonal.
  75 *> \endverbatim
  76 *>
  77 *> \param[in] LDB
  78 *> \verbatim
  79 *>          LDB is INTEGER
  80 *>          The leading dimension of the array B.  LDB >= 2.
  81 *> \endverbatim
  82 *>
  83 *> \param[in] SAFMIN
  84 *> \verbatim
  85 *>          SAFMIN is REAL
  86 *>          The smallest positive number s.t. 1/SAFMIN does not
  87 *>          overflow.  (This should always be SLAMCH('S') -- it is an
  88 *>          argument in order to avoid having to call SLAMCH frequently.)
  89 *> \endverbatim
  90 *>
  91 *> \param[out] SCALE1
  92 *> \verbatim
  93 *>          SCALE1 is REAL
  94 *>          A scaling factor used to avoid over-/underflow in the
  95 *>          eigenvalue equation which defines the first eigenvalue.  If
  96 *>          the eigenvalues are complex, then the eigenvalues are
  97 *>          ( WR1  +/-  WI i ) / SCALE1  (which may lie outside the
  98 *>          exponent range of the machine), SCALE1=SCALE2, and SCALE1
  99 *>          will always be positive.  If the eigenvalues are real, then
 100 *>          the first (real) eigenvalue is  WR1 / SCALE1 , but this may
 101 *>          overflow or underflow, and in fact, SCALE1 may be zero or
 102 *>          less than the underflow threshold if the exact eigenvalue
 103 *>          is sufficiently large.
 104 *> \endverbatim
 105 *>
 106 *> \param[out] SCALE2
 107 *> \verbatim
 108 *>          SCALE2 is REAL
 109 *>          A scaling factor used to avoid over-/underflow in the
 110 *>          eigenvalue equation which defines the second eigenvalue.  If
 111 *>          the eigenvalues are complex, then SCALE2=SCALE1.  If the
 112 *>          eigenvalues are real, then the second (real) eigenvalue is
 113 *>          WR2 / SCALE2 , but this may overflow or underflow, and in
 114 *>          fact, SCALE2 may be zero or less than the underflow
 115 *>          threshold if the exact eigenvalue is sufficiently large.
 116 *> \endverbatim
 117 *>
 118 *> \param[out] WR1
 119 *> \verbatim
 120 *>          WR1 is REAL
 121 *>          If the eigenvalue is real, then WR1 is SCALE1 times the
 122 *>          eigenvalue closest to the (2,2) element of A B**(-1).  If the
 123 *>          eigenvalue is complex, then WR1=WR2 is SCALE1 times the real
 124 *>          part of the eigenvalues.
 125 *> \endverbatim
 126 *>
 127 *> \param[out] WR2
 128 *> \verbatim
 129 *>          WR2 is REAL
 130 *>          If the eigenvalue is real, then WR2 is SCALE2 times the
 131 *>          other eigenvalue.  If the eigenvalue is complex, then
 132 *>          WR1=WR2 is SCALE1 times the real part of the eigenvalues.
 133 *> \endverbatim
 134 *>
 135 *> \param[out] WI
 136 *> \verbatim
 137 *>          WI is REAL
 138 *>          If the eigenvalue is real, then WI is zero.  If the
 139 *>          eigenvalue is complex, then WI is SCALE1 times the imaginary
 140 *>          part of the eigenvalues.  WI will always be non-negative.
 141 *> \endverbatim
 142 *
 143 *  Authors:
 144 *  ========
 145 *
 146 *> \author Univ. of Tennessee
 147 *> \author Univ. of California Berkeley
 148 *> \author Univ. of Colorado Denver
 149 *> \author NAG Ltd.
 150 *
 151 *> \date June 2016
 152 *
 153 *> \ingroup realOTHERauxiliary
 154 *
 155 *  =====================================================================
 156       SUBROUTINE SLAG2( A, LDA, B, LDB, SAFMIN, SCALE1, SCALE2, WR1,
 157      $                  WR2, WI )
 158 *
 159 *  -- LAPACK auxiliary routine (version 3.6.1) --
 160 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
 161 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
 162 *     June 2016
 163 *
 164 *     .. Scalar Arguments ..
 165       INTEGER            LDA, LDB
 166       REAL               SAFMIN, SCALE1, SCALE2, WI, WR1, WR2
 167 *     ..
 168 *     .. Array Arguments ..
 169       REAL               A( LDA, * ), B( LDB, * )
 170 *     ..
 171 *
 172 *  =====================================================================
 173 *
 174 *     .. Parameters ..
 175       REAL               ZERO, ONE, TWO
 176       PARAMETER          ( ZERO = 0.0E+0, ONE = 1.0E+0, TWO = 2.0E+0 )
 177       REAL               HALF
 178       PARAMETER          ( HALF = ONE / TWO )
 179       REAL               FUZZY1
 180       PARAMETER          ( FUZZY1 = ONE+1.0E-5 )
 181 *     ..
 182 *     .. Local Scalars ..
 183       REAL               A11, A12, A21, A22, ABI22, ANORM, AS11, AS12,
 184      $                   AS22, ASCALE, B11, B12, B22, BINV11, BINV22,
 185      $                   BMIN, BNORM, BSCALE, BSIZE, C1, C2, C3, C4, C5,
 186      $                   DIFF, DISCR, PP, QQ, R, RTMAX, RTMIN, S1, S2,
 187      $                   SAFMAX, SHIFT, SS, SUM, WABS, WBIG, WDET,
 188      $                   WSCALE, WSIZE, WSMALL
 189 *     ..
 190 *     .. Intrinsic Functions ..
 191       INTRINSIC          ABS, MAX, MIN, SIGN, SQRT
 192 *     ..
 193 *     .. Executable Statements ..
 194 *
 195       RTMIN = SQRT( SAFMIN )
 196       RTMAX = ONE / RTMIN
 197       SAFMAX = ONE / SAFMIN
 198 *
 199 *     Scale A
 200 *
 201       ANORM = MAX( ABS( A( 1, 1 ) )+ABS( A( 2, 1 ) ),
 202      $        ABS( A( 1, 2 ) )+ABS( A( 2, 2 ) ), SAFMIN )
 203       ASCALE = ONE / ANORM
 204       A11 = ASCALE*A( 1, 1 )
 205       A21 = ASCALE*A( 2, 1 )
 206       A12 = ASCALE*A( 1, 2 )
 207       A22 = ASCALE*A( 2, 2 )
 208 *
 209 *     Perturb B if necessary to insure non-singularity
 210 *
 211       B11 = B( 1, 1 )
 212       B12 = B( 1, 2 )
 213       B22 = B( 2, 2 )
 214       BMIN = RTMIN*MAX( ABS( B11 ), ABS( B12 ), ABS( B22 ), RTMIN )
 215       IF( ABS( B11 ).LT.BMIN )
 216      $   B11 = SIGN( BMIN, B11 )
 217       IF( ABS( B22 ).LT.BMIN )
 218      $   B22 = SIGN( BMIN, B22 )
 219 *
 220 *     Scale B
 221 *
 222       BNORM = MAX( ABS( B11 ), ABS( B12 )+ABS( B22 ), SAFMIN )
 223       BSIZE = MAX( ABS( B11 ), ABS( B22 ) )
 224       BSCALE = ONE / BSIZE
 225       B11 = B11*BSCALE
 226       B12 = B12*BSCALE
 227       B22 = B22*BSCALE
 228 *
 229 *     Compute larger eigenvalue by method described by C. van Loan
 230 *
 231 *     ( AS is A shifted by -SHIFT*B )
 232 *
 233       BINV11 = ONE / B11
 234       BINV22 = ONE / B22
 235       S1 = A11*BINV11
 236       S2 = A22*BINV22
 237       IF( ABS( S1 ).LE.ABS( S2 ) ) THEN
 238          AS12 = A12 - S1*B12
 239          AS22 = A22 - S1*B22
 240          SS = A21*( BINV11*BINV22 )
 241          ABI22 = AS22*BINV22 - SS*B12
 242          PP = HALF*ABI22
 243          SHIFT = S1
 244       ELSE
 245          AS12 = A12 - S2*B12
 246          AS11 = A11 - S2*B11
 247          SS = A21*( BINV11*BINV22 )
 248          ABI22 = -SS*B12
 249          PP = HALF*( AS11*BINV11+ABI22 )
 250          SHIFT = S2
 251       END IF
 252       QQ = SS*AS12
 253       IF( ABS( PP*RTMIN ).GE.ONE ) THEN
 254          DISCR = ( RTMIN*PP )**2 + QQ*SAFMIN
 255          R = SQRT( ABS( DISCR ) )*RTMAX
 256       ELSE
 257          IF( PP**2+ABS( QQ ).LE.SAFMIN ) THEN
 258             DISCR = ( RTMAX*PP )**2 + QQ*SAFMAX
 259             R = SQRT( ABS( DISCR ) )*RTMIN
 260          ELSE
 261             DISCR = PP**2 + QQ
 262             R = SQRT( ABS( DISCR ) )
 263          END IF
 264       END IF
 265 *
 266 *     Note: the test of R in the following IF is to cover the case when
 267 *           DISCR is small and negative and is flushed to zero during
 268 *           the calculation of R.  On machines which have a consistent
 269 *           flush-to-zero threshold and handle numbers above that
 270 *           threshold correctly, it would not be necessary.
 271 *
 272       IF( DISCR.GE.ZERO .OR. R.EQ.ZERO ) THEN
 273          SUM = PP + SIGN( R, PP )
 274          DIFF = PP - SIGN( R, PP )
 275          WBIG = SHIFT + SUM
 276 *
 277 *        Compute smaller eigenvalue
 278 *
 279          WSMALL = SHIFT + DIFF
 280          IF( HALF*ABS( WBIG ).GT.MAX( ABS( WSMALL ), SAFMIN ) ) THEN
 281             WDET = ( A11*A22-A12*A21 )*( BINV11*BINV22 )
 282             WSMALL = WDET / WBIG
 283          END IF
 284 *
 285 *        Choose (real) eigenvalue closest to 2,2 element of A*B**(-1)
 286 *        for WR1.
 287 *
 288          IF( PP.GT.ABI22 ) THEN
 289             WR1 = MIN( WBIG, WSMALL )
 290             WR2 = MAX( WBIG, WSMALL )
 291          ELSE
 292             WR1 = MAX( WBIG, WSMALL )
 293             WR2 = MIN( WBIG, WSMALL )
 294          END IF
 295          WI = ZERO
 296       ELSE
 297 *
 298 *        Complex eigenvalues
 299 *
 300          WR1 = SHIFT + PP
 301          WR2 = WR1
 302          WI = R
 303       END IF
 304 *
 305 *     Further scaling to avoid underflow and overflow in computing
 306 *     SCALE1 and overflow in computing w*B.
 307 *
 308 *     This scale factor (WSCALE) is bounded from above using C1 and C2,
 309 *     and from below using C3 and C4.
 310 *        C1 implements the condition  s A  must never overflow.
 311 *        C2 implements the condition  w B  must never overflow.
 312 *        C3, with C2,
 313 *           implement the condition that s A - w B must never overflow.
 314 *        C4 implements the condition  s    should not underflow.
 315 *        C5 implements the condition  max(s,|w|) should be at least 2.
 316 *
 317       C1 = BSIZE*( SAFMIN*MAX( ONE, ASCALE ) )
 318       C2 = SAFMIN*MAX( ONE, BNORM )
 319       C3 = BSIZE*SAFMIN
 320       IF( ASCALE.LE.ONE .AND. BSIZE.LE.ONE ) THEN
 321          C4 = MIN( ONE, ( ASCALE / SAFMIN )*BSIZE )
 322       ELSE
 323          C4 = ONE
 324       END IF
 325       IF( ASCALE.LE.ONE .OR. BSIZE.LE.ONE ) THEN
 326          C5 = MIN( ONE, ASCALE*BSIZE )
 327       ELSE
 328          C5 = ONE
 329       END IF
 330 *
 331 *     Scale first eigenvalue
 332 *
 333       WABS = ABS( WR1 ) + ABS( WI )
 334       WSIZE = MAX( SAFMIN, C1, FUZZY1*( WABS*C2+C3 ),
 335      $        MIN( C4, HALF*MAX( WABS, C5 ) ) )
 336       IF( WSIZE.NE.ONE ) THEN
 337          WSCALE = ONE / WSIZE
 338          IF( WSIZE.GT.ONE ) THEN
 339             SCALE1 = ( MAX( ASCALE, BSIZE )*WSCALE )*
 340      $               MIN( ASCALE, BSIZE )
 341          ELSE
 342             SCALE1 = ( MIN( ASCALE, BSIZE )*WSCALE )*
 343      $               MAX( ASCALE, BSIZE )
 344          END IF
 345          WR1 = WR1*WSCALE
 346          IF( WI.NE.ZERO ) THEN
 347             WI = WI*WSCALE
 348             WR2 = WR1
 349             SCALE2 = SCALE1
 350          END IF
 351       ELSE
 352          SCALE1 = ASCALE*BSIZE
 353          SCALE2 = SCALE1
 354       END IF
 355 *
 356 *     Scale second eigenvalue (if real)
 357 *
 358       IF( WI.EQ.ZERO ) THEN
 359          WSIZE = MAX( SAFMIN, C1, FUZZY1*( ABS( WR2 )*C2+C3 ),
 360      $           MIN( C4, HALF*MAX( ABS( WR2 ), C5 ) ) )
 361          IF( WSIZE.NE.ONE ) THEN
 362             WSCALE = ONE / WSIZE
 363             IF( WSIZE.GT.ONE ) THEN
 364                SCALE2 = ( MAX( ASCALE, BSIZE )*WSCALE )*
 365      $                  MIN( ASCALE, BSIZE )
 366             ELSE
 367                SCALE2 = ( MIN( ASCALE, BSIZE )*WSCALE )*
 368      $                  MAX( ASCALE, BSIZE )
 369             END IF
 370             WR2 = WR2*WSCALE
 371          ELSE
 372             SCALE2 = ASCALE*BSIZE
 373          END IF
 374       END IF
 375 *
 376 *     End of SLAG2
 377 *
 378       RETURN
 379       END