git.maemo.org Git - opencv/blob - 3rdparty/lapack/dlae2.c

   1 #include "clapack.h"
   2
   3 /* Subroutine */ int dlae2_(doublereal *a, doublereal *b, doublereal *c__,
   4         doublereal *rt1, doublereal *rt2)
   5 {
   6     /* System generated locals */
   7     doublereal d__1;
   8
   9     /* Builtin functions */
  10     double sqrt(doublereal);
  11
  12     /* Local variables */
  13     doublereal ab, df, tb, sm, rt, adf, acmn, acmx;
  14
  15
  16 /*  -- LAPACK auxiliary routine (version 3.1) -- */
  17 /*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */
  18 /*     November 2006 */
  19
  20 /*     .. Scalar Arguments .. */
  21 /*     .. */
  22
  23 /*  Purpose */
  24 /*  ======= */
  25
  26 /*  DLAE2  computes the eigenvalues of a 2-by-2 symmetric matrix */
  27 /*     [  A   B  ] */
  28 /*     [  B   C  ]. */
  29 /*  On return, RT1 is the eigenvalue of larger absolute value, and RT2 */
  30 /*  is the eigenvalue of smaller absolute value. */
  31
  32 /*  Arguments */
  33 /*  ========= */
  34
  35 /*  A       (input) DOUBLE PRECISION */
  36 /*          The (1,1) element of the 2-by-2 matrix. */
  37
  38 /*  B       (input) DOUBLE PRECISION */
  39 /*          The (1,2) and (2,1) elements of the 2-by-2 matrix. */
  40
  41 /*  C       (input) DOUBLE PRECISION */
  42 /*          The (2,2) element of the 2-by-2 matrix. */
  43
  44 /*  RT1     (output) DOUBLE PRECISION */
  45 /*          The eigenvalue of larger absolute value. */
  46
  47 /*  RT2     (output) DOUBLE PRECISION */
  48 /*          The eigenvalue of smaller absolute value. */
  49
  50 /*  Further Details */
  51 /*  =============== */
  52
  53 /*  RT1 is accurate to a few ulps barring over/underflow. */
  54
  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. */
  59
  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. */
  63
  64 /* ===================================================================== */
  65
  66 /*     .. Parameters .. */
  67 /*     .. */
  68 /*     .. Local Scalars .. */
  69 /*     .. */
  70 /*     .. Intrinsic Functions .. */
  71 /*     .. */
  72 /*     .. Executable Statements .. */
  73
  74 /*     Compute the eigenvalues */
  75
  76     sm = *a + *c__;
  77     df = *a - *c__;
  78     adf = abs(df);
  79     tb = *b + *b;
  80     ab = abs(tb);
  81     if (abs(*a) > abs(*c__)) {
  82         acmx = *a;
  83         acmn = *c__;
  84     } else {
  85         acmx = *c__;
  86         acmn = *a;
  87     }
  88     if (adf > ab) {
  89 /* Computing 2nd power */
  90         d__1 = ab / adf;
  91         rt = adf * sqrt(d__1 * d__1 + 1.);
  92     } else if (adf < ab) {
  93 /* Computing 2nd power */
  94         d__1 = adf / ab;
  95         rt = ab * sqrt(d__1 * d__1 + 1.);
  96     } else {
  97
  98 /*        Includes case AB=ADF=0 */
  99
 100         rt = ab * sqrt(2.);
 101     }
 102     if (sm < 0.) {
 103         *rt1 = (sm - rt) * .5;
 104
 105 /*        Order of execution important. */
 106 /*        To get fully accurate smaller eigenvalue, */
 107 /*        next line needs to be executed in higher precision. */
 108
 109         *rt2 = acmx / *rt1 * acmn - *b / *rt1 * *b;
 110     } else if (sm > 0.) {
 111         *rt1 = (sm + rt) * .5;
 112
 113 /*        Order of execution important. */
 114 /*        To get fully accurate smaller eigenvalue, */
 115 /*        next line needs to be executed in higher precision. */
 116
 117         *rt2 = acmx / *rt1 * acmn - *b / *rt1 * *b;
 118     } else {
 119
 120 /*        Includes case RT1 = RT2 = 0 */
 121
 122         *rt1 = rt * .5;
 123         *rt2 = rt * -.5;
 124     }
 125     return 0;
 126
 127 /*     End of DLAE2 */
 128
 129 } /* dlae2_ */