LAPACK: SRC/slanv2.f Source File

Go to the documentation of this file.
00001       SUBROUTINE SLANV2( A, B, C, D, RT1R, RT1I, RT2R, RT2I, CS, SN )
00002 *
00003 *  -- LAPACK auxiliary routine (version 3.2.2) --
00004 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
00005 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
00006 *     June 2010
00007 *
00008 *     .. Scalar Arguments ..
00009       REAL               A, B, C, CS, D, RT1I, RT1R, RT2I, RT2R, SN
00010 *     ..
00011 *
00012 *  Purpose
00013 *  =======
00014 *
00015 *  SLANV2 computes the Schur factorization of a real 2-by-2 nonsymmetric
00016 *  matrix in standard form:
00017 *
00018 *       [ A  B ] = [ CS -SN ] [ AA  BB ] [ CS  SN ]
00019 *       [ C  D ]   [ SN  CS ] [ CC  DD ] [-SN  CS ]
00020 *
00021 *  where either
00022 *  1) CC = 0 so that AA and DD are real eigenvalues of the matrix, or
00023 *  2) AA = DD and BB*CC < 0, so that AA + or - sqrt(BB*CC) are complex
00024 *  conjugate eigenvalues.
00025 *
00026 *  Arguments
00027 *  =========
00028 *
00029 *  A       (input/output) REAL            
00030 *  B       (input/output) REAL            
00031 *  C       (input/output) REAL            
00032 *  D       (input/output) REAL            
00033 *          On entry, the elements of the input matrix.
00034 *          On exit, they are overwritten by the elements of the
00035 *          standardised Schur form.
00036 *
00037 *  RT1R    (output) REAL 
00038 *  RT1I    (output) REAL            
00039 *  RT2R    (output) REAL            
00040 *  RT2I    (output) REAL            
00041 *          The real and imaginary parts of the eigenvalues. If the
00042 *          eigenvalues are a complex conjugate pair, RT1I > 0.
00043 *
00044 *  CS      (output) REAL            
00045 *  SN      (output) REAL            
00046 *          Parameters of the rotation matrix.
00047 *
00048 *  Further Details
00049 *  ===============
00050 *
00051 *  Modified by V. Sima, Research Institute for Informatics, Bucharest,
00052 *  Romania, to reduce the risk of cancellation errors,
00053 *  when computing real eigenvalues, and to ensure, if possible, that
00054 *  abs(RT1R) >= abs(RT2R).
00055 *
00056 *  =====================================================================
00057 *
00058 *     .. Parameters ..
00059       REAL               ZERO, HALF, ONE
00060       PARAMETER          ( ZERO = 0.0E+0, HALF = 0.5E+0, ONE = 1.0E+0 )
00061       REAL               MULTPL
00062       PARAMETER          ( MULTPL = 4.0E+0 )
00063 *     ..
00064 *     .. Local Scalars ..
00065       REAL               AA, BB, BCMAX, BCMIS, CC, CS1, DD, EPS, P, SAB,
00066      $                   SAC, SCALE, SIGMA, SN1, TAU, TEMP, Z
00067 *     ..
00068 *     .. External Functions ..
00069       REAL               SLAMCH, SLAPY2
00070       EXTERNAL           SLAMCH, SLAPY2
00071 *     ..
00072 *     .. Intrinsic Functions ..
00073       INTRINSIC          ABS, MAX, MIN, SIGN, SQRT
00074 *     ..
00075 *     .. Executable Statements ..
00076 *
00077       EPS = SLAMCH( 'P' )
00078       IF( C.EQ.ZERO ) THEN
00079          CS = ONE
00080          SN = ZERO
00081          GO TO 10
00082 *
00083       ELSE IF( B.EQ.ZERO ) THEN
00084 *
00085 *        Swap rows and columns
00086 *
00087          CS = ZERO
00088          SN = ONE
00089          TEMP = D
00090          D = A
00091          A = TEMP
00092          B = -C
00093          C = ZERO
00094          GO TO 10
00095       ELSE IF( (A-D).EQ.ZERO .AND. SIGN( ONE, B ).NE.
00096      $   SIGN( ONE, C ) ) THEN
00097          CS = ONE
00098          SN = ZERO
00099          GO TO 10
00100       ELSE
00101 *
00102          TEMP = A - D
00103          P = HALF*TEMP
00104          BCMAX = MAX( ABS( B ), ABS( C ) )
00105          BCMIS = MIN( ABS( B ), ABS( C ) )*SIGN( ONE, B )*SIGN( ONE, C )
00106          SCALE = MAX( ABS( P ), BCMAX )
00107          Z = ( P / SCALE )*P + ( BCMAX / SCALE )*BCMIS
00108 *
00109 *        If Z is of the order of the machine accuracy, postpone the
00110 *        decision on the nature of eigenvalues
00111 *
00112          IF( Z.GE.MULTPL*EPS ) THEN
00113 *
00114 *           Real eigenvalues. Compute A and D.
00115 *
00116             Z = P + SIGN( SQRT( SCALE )*SQRT( Z ), P )
00117             A = D + Z
00118             D = D - ( BCMAX / Z )*BCMIS
00119 *
00120 *           Compute B and the rotation matrix
00121 *
00122             TAU = SLAPY2( C, Z )
00123             CS = Z / TAU
00124             SN = C / TAU
00125             B = B - C
00126             C = ZERO
00127          ELSE
00128 *
00129 *           Complex eigenvalues, or real (almost) equal eigenvalues.
00130 *           Make diagonal elements equal.
00131 *
00132             SIGMA = B + C
00133             TAU = SLAPY2( SIGMA, TEMP )
00134             CS = SQRT( HALF*( ONE+ABS( SIGMA ) / TAU ) )
00135             SN = -( P / ( TAU*CS ) )*SIGN( ONE, SIGMA )
00136 *
00137 *           Compute [ AA  BB ] = [ A  B ] [ CS -SN ]
00138 *                   [ CC  DD ]   [ C  D ] [ SN  CS ]
00139 *
00140             AA = A*CS + B*SN
00141             BB = -A*SN + B*CS
00142             CC = C*CS + D*SN
00143             DD = -C*SN + D*CS
00144 *
00145 *           Compute [ A  B ] = [ CS  SN ] [ AA  BB ]
00146 *                   [ C  D ]   [-SN  CS ] [ CC  DD ]
00147 *
00148             A = AA*CS + CC*SN
00149             B = BB*CS + DD*SN
00150             C = -AA*SN + CC*CS
00151             D = -BB*SN + DD*CS
00152 *
00153             TEMP = HALF*( A+D )
00154             A = TEMP
00155             D = TEMP
00156 *
00157             IF( C.NE.ZERO ) THEN
00158                IF( B.NE.ZERO ) THEN
00159                   IF( SIGN( ONE, B ).EQ.SIGN( ONE, C ) ) THEN
00160 *
00161 *                    Real eigenvalues: reduce to upper triangular form
00162 *
00163                      SAB = SQRT( ABS( B ) )
00164                      SAC = SQRT( ABS( C ) )
00165                      P = SIGN( SAB*SAC, C )
00166                      TAU = ONE / SQRT( ABS( B+C ) )
00167                      A = TEMP + P
00168                      D = TEMP - P
00169                      B = B - C
00170                      C = ZERO
00171                      CS1 = SAB*TAU
00172                      SN1 = SAC*TAU
00173                      TEMP = CS*CS1 - SN*SN1
00174                      SN = CS*SN1 + SN*CS1
00175                      CS = TEMP
00176                   END IF
00177                ELSE
00178                   B = -C
00179                   C = ZERO
00180                   TEMP = CS
00181                   CS = -SN
00182                   SN = TEMP
00183                END IF
00184             END IF
00185          END IF
00186 *
00187       END IF
00188 *
00189    10 CONTINUE
00190 *
00191 *     Store eigenvalues in (RT1R,RT1I) and (RT2R,RT2I).
00192 *
00193       RT1R = A
00194       RT2R = D
00195       IF( C.EQ.ZERO ) THEN
00196          RT1I = ZERO
00197          RT2I = ZERO
00198       ELSE
00199          RT1I = SQRT( ABS( B ) )*SQRT( ABS( C ) )
00200          RT2I = -RT1I
00201       END IF
00202       RETURN
00203 *
00204 *     End of SLANV2
00205 *
00206       END