3 /* Subroutine */ int dlae2_(doublereal *a, doublereal *b, doublereal *c__,
4 doublereal *rt1, doublereal *rt2)
6 /* System generated locals */
9 /* Builtin functions */
10 double sqrt(doublereal);
13 doublereal ab, df, tb, sm, rt, adf, acmn, acmx;
16 /* -- LAPACK auxiliary routine (version 3.1) -- */
17 /* Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */
20 /* .. Scalar Arguments .. */
26 /* DLAE2 computes the eigenvalues of a 2-by-2 symmetric matrix */
29 /* On return, RT1 is the eigenvalue of larger absolute value, and RT2 */
30 /* is the eigenvalue of smaller absolute value. */
35 /* A (input) DOUBLE PRECISION */
36 /* The (1,1) element of the 2-by-2 matrix. */
38 /* B (input) DOUBLE PRECISION */
39 /* The (1,2) and (2,1) elements of the 2-by-2 matrix. */
41 /* C (input) DOUBLE PRECISION */
42 /* The (2,2) element of the 2-by-2 matrix. */
44 /* RT1 (output) DOUBLE PRECISION */
45 /* The eigenvalue of larger absolute value. */
47 /* RT2 (output) DOUBLE PRECISION */
48 /* The eigenvalue of smaller absolute value. */
53 /* RT1 is accurate to a few ulps barring over/underflow. */
55 /* RT2 may be inaccurate if there is massive cancellation in the */
56 /* determinant A*C-B*B; higher precision or correctly rounded or */
57 /* correctly truncated arithmetic would be needed to compute RT2 */
58 /* accurately in all cases. */
60 /* Overflow is possible only if RT1 is within a factor of 5 of overflow. */
61 /* Underflow is harmless if the input data is 0 or exceeds */
62 /* underflow_threshold / macheps. */
64 /* ===================================================================== */
66 /* .. Parameters .. */
68 /* .. Local Scalars .. */
70 /* .. Intrinsic Functions .. */
72 /* .. Executable Statements .. */
74 /* Compute the eigenvalues */
81 if (abs(*a) > abs(*c__)) {
89 /* Computing 2nd power */
91 rt = adf * sqrt(d__1 * d__1 + 1.);
92 } else if (adf < ab) {
93 /* Computing 2nd power */
95 rt = ab * sqrt(d__1 * d__1 + 1.);
98 /* Includes case AB=ADF=0 */
103 *rt1 = (sm - rt) * .5;
105 /* Order of execution important. */
106 /* To get fully accurate smaller eigenvalue, */
107 /* next line needs to be executed in higher precision. */
109 *rt2 = acmx / *rt1 * acmn - *b / *rt1 * *b;
110 } else if (sm > 0.) {
111 *rt1 = (sm + rt) * .5;
113 /* Order of execution important. */
114 /* To get fully accurate smaller eigenvalue, */
115 /* next line needs to be executed in higher precision. */
117 *rt2 = acmx / *rt1 * acmn - *b / *rt1 * *b;
120 /* Includes case RT1 = RT2 = 0 */