home *** CD-ROM | disk | FTP | other *** search
/ Geek Gadgets 1 / ADE-1.bin / ade-dist / octave-1.1.1p1-src.tgz / tar.out / fsf / octave / libcruft / lapack / dlanv2.f < prev    next >
Text File  |  1996-09-28  |  5KB  |  176 lines
  1.       SUBROUTINE DLANV2( A, B, C, D, RT1R, RT1I, RT2R, RT2I, CS, SN )
  2. *
  3. *  -- LAPACK auxiliary routine (version 2.0) --
  4. *     Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
  5. *     Courant Institute, Argonne National Lab, and Rice University
  6. *     September 30, 1994 
  7. *
  8. *     .. Scalar Arguments ..
  9.       DOUBLE PRECISION   A, B, C, CS, D, RT1I, RT1R, RT2I, RT2R, SN
  10. *     ..
  11. *
  12. *  Purpose
  13. *  =======
  14. *
  15. *  DLANV2 computes the Schur factorization of a real 2-by-2 nonsymmetric
  16. *  matrix in standard form:
  17. *
  18. *       [ A  B ] = [ CS -SN ] [ AA  BB ] [ CS  SN ]
  19. *       [ C  D ]   [ SN  CS ] [ CC  DD ] [-SN  CS ]
  20. *
  21. *  where either
  22. *  1) CC = 0 so that AA and DD are real eigenvalues of the matrix, or
  23. *  2) AA = DD and BB*CC < 0, so that AA + or - sqrt(BB*CC) are complex
  24. *  conjugate eigenvalues.
  25. *
  26. *  Arguments
  27. *  =========
  28. *
  29. *  A       (input/output) DOUBLE PRECISION
  30. *  B       (input/output) DOUBLE PRECISION
  31. *  C       (input/output) DOUBLE PRECISION
  32. *  D       (input/output) DOUBLE PRECISION
  33. *          On entry, the elements of the input matrix.
  34. *          On exit, they are overwritten by the elements of the
  35. *          standardised Schur form.
  36. *
  37. *  RT1R    (output) DOUBLE PRECISION
  38. *  RT1I    (output) DOUBLE PRECISION
  39. *  RT2R    (output) DOUBLE PRECISION
  40. *  RT2I    (output) DOUBLE PRECISION
  41. *          The real and imaginary parts of the eigenvalues. If the
  42. *          eigenvalues are both real, abs(RT1R) >= abs(RT2R); if the
  43. *          eigenvalues are a complex conjugate pair, RT1I > 0.
  44. *
  45. *  CS      (output) DOUBLE PRECISION
  46. *  SN      (output) DOUBLE PRECISION
  47. *          Parameters of the rotation matrix.
  48. *
  49. *  =====================================================================
  50. *
  51. *     .. Parameters ..
  52.       DOUBLE PRECISION   ZERO, HALF, ONE
  53.       PARAMETER          ( ZERO = 0.0D+0, HALF = 0.5D+0, ONE = 1.0D+0 )
  54. *     ..
  55. *     .. Local Scalars ..
  56.       DOUBLE PRECISION   AA, BB, CC, CS1, DD, P, SAB, SAC, SIGMA, SN1,
  57.      $                   TAU, TEMP
  58. *     ..
  59. *     .. External Functions ..
  60.       DOUBLE PRECISION   DLAPY2
  61.       EXTERNAL           DLAPY2
  62. *     ..
  63. *     .. Intrinsic Functions ..
  64.       INTRINSIC          ABS, SIGN, SQRT
  65. *     ..
  66. *     .. Executable Statements ..
  67. *
  68. *     Initialize CS and SN
  69. *
  70.       CS = ONE
  71.       SN = ZERO
  72. *
  73.       IF( C.EQ.ZERO ) THEN
  74.          GO TO 10
  75. *
  76.       ELSE IF( B.EQ.ZERO ) THEN
  77. *
  78. *        Swap rows and columns
  79. *
  80.          CS = ZERO
  81.          SN = ONE
  82.          TEMP = D
  83.          D = A
  84.          A = TEMP
  85.          B = -C
  86.          C = ZERO
  87.          GO TO 10
  88.       ELSE IF( (A-D).EQ.ZERO .AND. SIGN( ONE, B ).NE.
  89.      $   SIGN( ONE, C ) ) THEN
  90.          GO TO 10
  91.       ELSE
  92. *
  93. *        Make diagonal elements equal
  94. *
  95.          TEMP = A - D
  96.          P = HALF*TEMP
  97.          SIGMA = B + C
  98.          TAU = DLAPY2( SIGMA, TEMP )
  99.          CS1 = SQRT( HALF*( ONE+ABS( SIGMA ) / TAU ) )
  100.          SN1 = -( P / ( TAU*CS1 ) )*SIGN( ONE, SIGMA )
  101. *
  102. *        Compute [ AA  BB ] = [ A  B ] [ CS1 -SN1 ]
  103. *                [ CC  DD ]   [ C  D ] [ SN1  CS1 ]
  104. *
  105.          AA = A*CS1 + B*SN1
  106.          BB = -A*SN1 + B*CS1
  107.          CC = C*CS1 + D*SN1
  108.          DD = -C*SN1 + D*CS1
  109. *
  110. *        Compute [ A  B ] = [ CS1  SN1 ] [ AA  BB ]
  111. *                [ C  D ]   [-SN1  CS1 ] [ CC  DD ]
  112. *
  113.          A = AA*CS1 + CC*SN1
  114.          B = BB*CS1 + DD*SN1
  115.          C = -AA*SN1 + CC*CS1
  116.          D = -BB*SN1 + DD*CS1
  117. *
  118. *        Accumulate transformation
  119. *
  120.          TEMP = CS*CS1 - SN*SN1
  121.          SN = CS*SN1 + SN*CS1
  122.          CS = TEMP
  123. *
  124.          TEMP = HALF*( A+D )
  125.          A = TEMP
  126.          D = TEMP
  127. *
  128.          IF( C.NE.ZERO ) THEN
  129.             IF( B.NE.ZERO ) THEN
  130.                IF( SIGN( ONE, B ).EQ.SIGN( ONE, C ) ) THEN
  131. *
  132. *                 Real eigenvalues: reduce to upper triangular form
  133. *
  134.                   SAB = SQRT( ABS( B ) )
  135.                   SAC = SQRT( ABS( C ) )
  136.                   P = SIGN( SAB*SAC, C )
  137.                   TAU = ONE / SQRT( ABS( B+C ) )
  138.                   A = TEMP + P
  139.                   D = TEMP - P
  140.                   B = B - C
  141.                   C = ZERO
  142.                   CS1 = SAB*TAU
  143.                   SN1 = SAC*TAU
  144.                   TEMP = CS*CS1 - SN*SN1
  145.                   SN = CS*SN1 + SN*CS1
  146.                   CS = TEMP
  147.                END IF
  148.             ELSE
  149.                B = -C
  150.                C = ZERO
  151.                TEMP = CS
  152.                CS = -SN
  153.                SN = TEMP
  154.             END IF
  155.          END IF
  156.       END IF
  157. *
  158.    10 CONTINUE
  159. *
  160. *     Store eigenvalues in (RT1R,RT1I) and (RT2R,RT2I).
  161. *
  162.       RT1R = A
  163.       RT2R = D
  164.       IF( C.EQ.ZERO ) THEN
  165.          RT1I = ZERO
  166.          RT2I = ZERO
  167.       ELSE
  168.          RT1I = SQRT( ABS( B ) )*SQRT( ABS( C ) )
  169.          RT2I = -RT1I
  170.       END IF
  171.       RETURN
  172. *
  173. *     End of DLANV2
  174. *
  175.       END
  176.