SRC/slanv2.f

   1 *> \brief \b SLANV2 computes the Schur factorization of a real 2-by-2 nonsymmetric matrix in standard form.
   2 *
   3 *  =========== DOCUMENTATION ===========
   4 *
   5 * Online html documentation available at
   6 *            http://www.netlib.org/lapack/explore-html/
   7 *
   8 *> \htmlonly
   9 *> Download SLANV2 + dependencies
  10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/slanv2.f">
  11 *> [TGZ]</a>
  12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/slanv2.f">
  13 *> [ZIP]</a>
  14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/slanv2.f">
  15 *> [TXT]</a>
  16 *> \endhtmlonly
  17 *
  18 *  Definition:
  19 *  ===========
  20 *
  21 *       SUBROUTINE SLANV2( A, B, C, D, RT1R, RT1I, RT2R, RT2I, CS, SN )
  22 *
  23 *       .. Scalar Arguments ..
  24 *       REAL               A, B, C, CS, D, RT1I, RT1R, RT2I, RT2R, SN
  25 *       ..
  26 *
  27 *
  28 *> \par Purpose:
  29 *  =============
  30 *>
  31 *> \verbatim
  32 *>
  33 *> SLANV2 computes the Schur factorization of a real 2-by-2 nonsymmetric
  34 *> matrix in standard form:
  35 *>
  36 *>      [ A  B ] = [ CS -SN ] [ AA  BB ] [ CS  SN ]
  37 *>      [ C  D ]   [ SN  CS ] [ CC  DD ] [-SN  CS ]
  38 *>
  39 *> where either
  40 *> 1) CC = 0 so that AA and DD are real eigenvalues of the matrix, or
  41 *> 2) AA = DD and BB*CC < 0, so that AA + or - sqrt(BB*CC) are complex
  42 *> conjugate eigenvalues.
  43 *> \endverbatim
  44 *
  45 *  Arguments:
  46 *  ==========
  47 *
  48 *> \param[in,out] A
  49 *> \verbatim
  50 *>          A is REAL
  51 *> \endverbatim
  52 *>
  53 *> \param[in,out] B
  54 *> \verbatim
  55 *>          B is REAL
  56 *> \endverbatim
  57 *>
  58 *> \param[in,out] C
  59 *> \verbatim
  60 *>          C is REAL
  61 *> \endverbatim
  62 *>
  63 *> \param[in,out] D
  64 *> \verbatim
  65 *>          D is REAL
  66 *>          On entry, the elements of the input matrix.
  67 *>          On exit, they are overwritten by the elements of the
  68 *>          standardised Schur form.
  69 *> \endverbatim
  70 *>
  71 *> \param[out] RT1R
  72 *> \verbatim
  73 *>          RT1R is REAL
  74 *> \endverbatim
  75 *>
  76 *> \param[out] RT1I
  77 *> \verbatim
  78 *>          RT1I is REAL
  79 *> \endverbatim
  80 *>
  81 *> \param[out] RT2R
  82 *> \verbatim
  83 *>          RT2R is REAL
  84 *> \endverbatim
  85 *>
  86 *> \param[out] RT2I
  87 *> \verbatim
  88 *>          RT2I is REAL
  89 *>          The real and imaginary parts of the eigenvalues. If the
  90 *>          eigenvalues are a complex conjugate pair, RT1I > 0.
  91 *> \endverbatim
  92 *>
  93 *> \param[out] CS
  94 *> \verbatim
  95 *>          CS is REAL
  96 *> \endverbatim
  97 *>
  98 *> \param[out] SN
  99 *> \verbatim
 100 *>          SN is REAL
 101 *>          Parameters of the rotation matrix.
 102 *> \endverbatim
 103 *
 104 *  Authors:
 105 *  ========
 106 *
 107 *> \author Univ. of Tennessee
 108 *> \author Univ. of California Berkeley
 109 *> \author Univ. of Colorado Denver
 110 *> \author NAG Ltd.
 111 *
 112 *> \date September 2012
 113 *
 114 *> \ingroup realOTHERauxiliary
 115 *
 116 *> \par Further Details:
 117 *  =====================
 118 *>
 119 *> \verbatim
 120 *>
 121 *>  Modified by V. Sima, Research Institute for Informatics, Bucharest,
 122 *>  Romania, to reduce the risk of cancellation errors,
 123 *>  when computing real eigenvalues, and to ensure, if possible, that
 124 *>  abs(RT1R) >= abs(RT2R).
 125 *> \endverbatim
 126 *>
 127 *  =====================================================================
 128       SUBROUTINE SLANV2( A, B, C, D, RT1R, RT1I, RT2R, RT2I, CS, SN )
 129 *
 130 *  -- LAPACK auxiliary routine (version 3.4.2) --
 131 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
 132 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
 133 *     September 2012
 134 *
 135 *     .. Scalar Arguments ..
 136       REAL               A, B, C, CS, D, RT1I, RT1R, RT2I, RT2R, SN
 137 *     ..
 138 *
 139 *  =====================================================================
 140 *
 141 *     .. Parameters ..
 142       REAL               ZERO, HALF, ONE
 143       PARAMETER          ( ZERO = 0.0E+0, HALF = 0.5E+0, ONE = 1.0E+0 )
 144       REAL               MULTPL
 145       PARAMETER          ( MULTPL = 4.0E+0 )
 146 *     ..
 147 *     .. Local Scalars ..
 148       REAL               AA, BB, BCMAX, BCMIS, CC, CS1, DD, EPS, P, SAB,
 149      $                   SAC, SCALE, SIGMA, SN1, TAU, TEMP, Z
 150 *     ..
 151 *     .. External Functions ..
 152       REAL               SLAMCH, SLAPY2
 153       EXTERNAL           SLAMCH, SLAPY2
 154 *     ..
 155 *     .. Intrinsic Functions ..
 156       INTRINSIC          ABS, MAX, MIN, SIGN, SQRT
 157 *     ..
 158 *     .. Executable Statements ..
 159 *
 160       EPS = SLAMCH( 'P' )
 161       IF( C.EQ.ZERO ) THEN
 162          CS = ONE
 163          SN = ZERO
 164          GO TO 10
 165 *
 166       ELSE IF( B.EQ.ZERO ) THEN
 167 *
 168 *        Swap rows and columns
 169 *
 170          CS = ZERO
 171          SN = ONE
 172          TEMP = D
 173          D = A
 174          A = TEMP
 175          B = -C
 176          C = ZERO
 177          GO TO 10
 178       ELSE IF( (A-D).EQ.ZERO .AND. SIGN( ONE, B ).NE.
 179      $   SIGN( ONE, C ) ) THEN
 180          CS = ONE
 181          SN = ZERO
 182          GO TO 10
 183       ELSE
 184 *
 185          TEMP = A - D
 186          P = HALF*TEMP
 187          BCMAX = MAX( ABS( B ), ABS( C ) )
 188          BCMIS = MIN( ABS( B ), ABS( C ) )*SIGN( ONE, B )*SIGN( ONE, C )
 189          SCALE = MAX( ABS( P ), BCMAX )
 190          Z = ( P / SCALE )*P + ( BCMAX / SCALE )*BCMIS
 191 *
 192 *        If Z is of the order of the machine accuracy, postpone the
 193 *        decision on the nature of eigenvalues
 194 *
 195          IF( Z.GE.MULTPL*EPS ) THEN
 196 *
 197 *           Real eigenvalues. Compute A and D.
 198 *
 199             Z = P + SIGN( SQRT( SCALE )*SQRT( Z ), P )
 200             A = D + Z
 201             D = D - ( BCMAX / Z )*BCMIS
 202 *
 203 *           Compute B and the rotation matrix
 204 *
 205             TAU = SLAPY2( C, Z )
 206             CS = Z / TAU
 207             SN = C / TAU
 208             B = B - C
 209             C = ZERO
 210          ELSE
 211 *
 212 *           Complex eigenvalues, or real (almost) equal eigenvalues.
 213 *           Make diagonal elements equal.
 214 *
 215             SIGMA = B + C
 216             TAU = SLAPY2( SIGMA, TEMP )
 217             CS = SQRT( HALF*( ONE+ABS( SIGMA ) / TAU ) )
 218             SN = -( P / ( TAU*CS ) )*SIGN( ONE, SIGMA )
 219 *
 220 *           Compute [ AA  BB ] = [ A  B ] [ CS -SN ]
 221 *                   [ CC  DD ]   [ C  D ] [ SN  CS ]
 222 *
 223             AA = A*CS + B*SN
 224             BB = -A*SN + B*CS
 225             CC = C*CS + D*SN
 226             DD = -C*SN + D*CS
 227 *
 228 *           Compute [ A  B ] = [ CS  SN ] [ AA  BB ]
 229 *                   [ C  D ]   [-SN  CS ] [ CC  DD ]
 230 *
 231             A = AA*CS + CC*SN
 232             B = BB*CS + DD*SN
 233             C = -AA*SN + CC*CS
 234             D = -BB*SN + DD*CS
 235 *
 236             TEMP = HALF*( A+D )
 237             A = TEMP
 238             D = TEMP
 239 *
 240             IF( C.NE.ZERO ) THEN
 241                IF( B.NE.ZERO ) THEN
 242                   IF( SIGN( ONE, B ).EQ.SIGN( ONE, C ) ) THEN
 243 *
 244 *                    Real eigenvalues: reduce to upper triangular form
 245 *
 246                      SAB = SQRT( ABS( B ) )
 247                      SAC = SQRT( ABS( C ) )
 248                      P = SIGN( SAB*SAC, C )
 249                      TAU = ONE / SQRT( ABS( B+C ) )
 250                      A = TEMP + P
 251                      D = TEMP - P
 252                      B = B - C
 253                      C = ZERO
 254                      CS1 = SAB*TAU
 255                      SN1 = SAC*TAU
 256                      TEMP = CS*CS1 - SN*SN1
 257                      SN = CS*SN1 + SN*CS1
 258                      CS = TEMP
 259                   END IF
 260                ELSE
 261                   B = -C
 262                   C = ZERO
 263                   TEMP = CS
 264                   CS = -SN
 265                   SN = TEMP
 266                END IF
 267             END IF
 268          END IF
 269 *
 270       END IF
 271 *
 272    10 CONTINUE
 273 *
 274 *     Store eigenvalues in (RT1R,RT1I) and (RT2R,RT2I).
 275 *
 276       RT1R = A
 277       RT2R = D
 278       IF( C.EQ.ZERO ) THEN
 279          RT1I = ZERO
 280          RT2I = ZERO
 281       ELSE
 282          RT1I = SQRT( ABS( B ) )*SQRT( ABS( C ) )
 283          RT2I = -RT1I
 284       END IF
 285       RETURN
 286 *
 287 *     End of SLANV2
 288 *
 289       END