3 /* Subroutine */ int slae2_(real *a, real *b, real *c__, real *rt1, real *rt2)
5 /* System generated locals */
8 /* Builtin functions */
9 double sqrt(doublereal);
12 real ab, df, tb, sm, rt, adf, acmn, acmx;
15 /* -- LAPACK auxiliary routine (version 3.1) -- */
16 /* Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */
19 /* .. Scalar Arguments .. */
25 /* SLAE2 computes the eigenvalues of a 2-by-2 symmetric matrix */
28 /* On return, RT1 is the eigenvalue of larger absolute value, and RT2 */
29 /* is the eigenvalue of smaller absolute value. */
35 /* The (1,1) element of the 2-by-2 matrix. */
38 /* The (1,2) and (2,1) elements of the 2-by-2 matrix. */
41 /* The (2,2) element of the 2-by-2 matrix. */
43 /* RT1 (output) REAL */
44 /* The eigenvalue of larger absolute value. */
46 /* RT2 (output) REAL */
47 /* The eigenvalue of smaller absolute value. */
52 /* RT1 is accurate to a few ulps barring over/underflow. */
54 /* RT2 may be inaccurate if there is massive cancellation in the */
55 /* determinant A*C-B*B; higher precision or correctly rounded or */
56 /* correctly truncated arithmetic would be needed to compute RT2 */
57 /* accurately in all cases. */
59 /* Overflow is possible only if RT1 is within a factor of 5 of overflow. */
60 /* Underflow is harmless if the input data is 0 or exceeds */
61 /* underflow_threshold / macheps. */
63 /* ===================================================================== */
65 /* .. Parameters .. */
67 /* .. Local Scalars .. */
69 /* .. Intrinsic Functions .. */
71 /* .. Executable Statements .. */
73 /* Compute the eigenvalues */
80 if (dabs(*a) > dabs(*c__)) {
88 /* Computing 2nd power */
90 rt = adf * sqrt(r__1 * r__1 + 1.f);
91 } else if (adf < ab) {
92 /* Computing 2nd power */
94 rt = ab * sqrt(r__1 * r__1 + 1.f);
97 /* Includes case AB=ADF=0 */
102 *rt1 = (sm - rt) * .5f;
104 /* Order of execution important. */
105 /* To get fully accurate smaller eigenvalue, */
106 /* next line needs to be executed in higher precision. */
108 *rt2 = acmx / *rt1 * acmn - *b / *rt1 * *b;
109 } else if (sm > 0.f) {
110 *rt1 = (sm + rt) * .5f;
112 /* Order of execution important. */
113 /* To get fully accurate smaller eigenvalue, */
114 /* next line needs to be executed in higher precision. */
116 *rt2 = acmx / *rt1 * acmn - *b / *rt1 * *b;
119 /* Includes case RT1 = RT2 = 0 */