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

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