Mercurial > octave

comparison libcruft/lapack/slae2.f @ 7789:82be108cc558

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
First attempt at single precision tyeps * * * corrections to qrupdate single precision routines * * * prefer demotion to single over promotion to double * * * Add single precision support to log2 function * * * Trivial PROJECT file update * * * Cache optimized hermitian/transpose methods * * * Add tests for tranpose/hermitian and ChangeLog entry for new transpose code
author David Bateman <dbateman@free.fr>
date Sun, 27 Apr 2008 22:34:17 +0200
parents
children
comparison
equal deleted inserted replaced
7788:45f5faba05a2 7789:82be108cc558
1 SUBROUTINE SLAE2( A, B, C, RT1, RT2 )
2 *
3 * -- LAPACK auxiliary routine (version 3.1) --
4 * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
5 * November 2006
6 *
7 * .. Scalar Arguments ..
8 REAL A, B, C, RT1, RT2
9 * ..
10 *
11 * Purpose
12 * =======
13 *
14 * SLAE2 computes the eigenvalues of a 2-by-2 symmetric matrix
15 * [ A B ]
16 * [ B C ].
17 * On return, RT1 is the eigenvalue of larger absolute value, and RT2
18 * is the eigenvalue of smaller absolute value.
19 *
20 * Arguments
21 * =========
22 *
23 * A (input) REAL
24 * The (1,1) element of the 2-by-2 matrix.
25 *
26 * B (input) REAL
27 * The (1,2) and (2,1) elements of the 2-by-2 matrix.
28 *
29 * C (input) REAL
30 * The (2,2) element of the 2-by-2 matrix.
31 *
32 * RT1 (output) REAL
33 * The eigenvalue of larger absolute value.
34 *
35 * RT2 (output) REAL
36 * The eigenvalue of smaller absolute value.
37 *
38 * Further Details
39 * ===============
40 *
41 * RT1 is accurate to a few ulps barring over/underflow.
42 *
43 * RT2 may be inaccurate if there is massive cancellation in the
44 * determinant A*C-B*B; higher precision or correctly rounded or
45 * correctly truncated arithmetic would be needed to compute RT2
46 * accurately in all cases.
47 *
48 * Overflow is possible only if RT1 is within a factor of 5 of overflow.
49 * Underflow is harmless if the input data is 0 or exceeds
50 * underflow_threshold / macheps.
51 *
52 * =====================================================================
53 *
54 * .. Parameters ..
55 REAL ONE
56 PARAMETER ( ONE = 1.0E0 )
57 REAL TWO
58 PARAMETER ( TWO = 2.0E0 )
59 REAL ZERO
60 PARAMETER ( ZERO = 0.0E0 )
61 REAL HALF
62 PARAMETER ( HALF = 0.5E0 )
63 * ..
64 * .. Local Scalars ..
65 REAL AB, ACMN, ACMX, ADF, DF, RT, SM, TB
66 * ..
67 * .. Intrinsic Functions ..
68 INTRINSIC ABS, SQRT
69 * ..
70 * .. Executable Statements ..
71 *
72 * Compute the eigenvalues
73 *
74 SM = A + C
75 DF = A - C
76 ADF = ABS( DF )
77 TB = B + B
78 AB = ABS( TB )
79 IF( ABS( A ).GT.ABS( C ) ) THEN
80 ACMX = A
81 ACMN = C
82 ELSE
83 ACMX = C
84 ACMN = A
85 END IF
86 IF( ADF.GT.AB ) THEN
87 RT = ADF*SQRT( ONE+( AB / ADF )**2 )
88 ELSE IF( ADF.LT.AB ) THEN
89 RT = AB*SQRT( ONE+( ADF / AB )**2 )
90 ELSE
91 *
92 * Includes case AB=ADF=0
93 *
94 RT = AB*SQRT( TWO )
95 END IF
96 IF( SM.LT.ZERO ) THEN
97 RT1 = HALF*( SM-RT )
98 *
99 * Order of execution important.
100 * To get fully accurate smaller eigenvalue,
101 * next line needs to be executed in higher precision.
102 *
103 RT2 = ( ACMX / RT1 )*ACMN - ( B / RT1 )*B
104 ELSE IF( SM.GT.ZERO ) THEN
105 RT1 = HALF*( SM+RT )
106 *
107 * Order of execution important.
108 * To get fully accurate smaller eigenvalue,
109 * next line needs to be executed in higher precision.
110 *
111 RT2 = ( ACMX / RT1 )*ACMN - ( B / RT1 )*B
112 ELSE
113 *
114 * Includes case RT1 = RT2 = 0
115 *
116 RT1 = HALF*RT
117 RT2 = -HALF*RT
118 END IF
119 RETURN
120 *
121 * End of SLAE2
122 *
123 END